Numerical simulation and dynamic mode decomposition analysis of the flow past wings with spanwise waviness

This paper presents an investigation on the compressible flow past spanwise wavy wings with an incoming flow of and . The scale adaptive simulation is employed here to solve the Favre-averaged equations. Three wings, with the ratio of wavelengths (λ) to chord length (c) respectively equaling to 1/2 (W4), 1/3 (W6), 1/4 (W8), are aerodynamically investigated. The aerodynamic force and separation were compared. The results obtained show that the lift-to-drag ratio of W-8 is 12.7% and 14.4% higher than W-4 and W-6 because of the larger low-pressure regions, while the lift and drag coefficient curves of W-8 oscillate more strongly. Leading edge separated bubbles and dilation-type recirculation regions are observed in the flow field structures. Spanwise momentum transportation and local airflow downwash effect are generated by the flow around the separation regions. Spanwise momentum transportation and downwash effect play important roles in the flow control mechanism of spanwise wavy wings. Three components of three-dimensional velocity vectors are modally decomposed by dynamic mode decomposition method to reveal the flow control mechanisms and vortex structures. The main modes characterizing the motions of arched vortices and shear layer vortices have been obtained. Strouhal number (St) is adopted to define the periodic oscillations of vortex shedding. The St of the modes are related to the scale of the vortex structures in the flow field. Small-scale eddy motions have been decomposed in high-frequency modes, and thus the modes of the flow field are more meticulous for higher St.


Introduction
Humpback whales peocoral flippers are approximately 550 cm in length and are shaped like bird wings. The front edge of the flippers has irregular tumor-like tubercles and the rear edge has wavy cuts, making humpback whales one of the world's largest animals to possess amazing agility. A study by Fish and Battle (1995) showed that the raised nodules on the front edge of the humpback whales peocoral flippers can delay stall and maintain a high lift coefficient at high angles of attack to provide the necessary centripetal force for turning maneuvers. In recent years, based on the flippers of humpback whales, the bionic wavy modified wings have attracted great attention as new passive methods for flow separation control, include the wings with leading edge waviness, trailing edge waviness and spanwise waviness, of which wavy trailing edge is often used in noise control related research.
The water tunnel force measurements of NACA0020 full-span and half-span wings with or without wavy leading edges were conducted by Miklosovic et al. (2004, CONTACT Jian-Hong Sun jhsun@nuaa.edu.cn 2007). Experiments show that for a half-span wing with a Reynolds number of about 5.0 × 10 5 , the tubercles of the leading edge delay the stall angle of attack, and the maximum lift coefficient is increased by 6%, which significantly improves the stall performance. For a fullspan model with a Reynolds number of about 2.7 × 10 5 , the leading edge tubercles cause the wing to stall prematurely, which is only effective in the post-stall phase. Carreira and Kobayashi (2008) used detached eddy simulation to verify the effect of the wavy tubercles on the stall characteristics of the whale peocoral flipper model. The effect of wavy leading edge on hydrodynamic characteristics was investigated by Yoon et al. (2011) through solving Reynolds-averaged Navier-Stokes equations (RANS). The study shows that under the condition of small angles of attack, the wavy front has negligible influence on the flow structure and the lift drag coefficient; for the flow at high angles of attack, the wavy leading edge causes the flow to produce an unfavorable pressure gradient in the trough area and separate the flow. Stanway (2008) conducted an experiment similar to Miklosovic's half-span model. By analyzing the experimental results, it is believed that the Reynolds number is an important factor in determining whether the wavy tubercles improve the aerodynamic performance of the wing. The large eddy simulations (LES) of the cavitating flow around wavy leading edge hydrofoil conducted by Pendar et al. (2020) exhibit a significantly reduced level of unsteady fluctuations in aerodynamic forces. Experimental study of wavy leading edge wing was carried out by Mehraban and Djavareshkian (2021) under ground effect. The wavy leading edge wing shows an excellent performance in the poststall region. An aerodynamic optimization study is conducted by Lu et al. (2021) on the wings with wavy leading edge. The increments of the stall angles for optimal wings run up to as much as 33.3%. Based on a large number of experiments and numerical simulations, the structure of the flow field for the wavy modified wings has been extensively studied. Cai et al. (2017) and Natarajan et al. (2014) carried out the particle image velocimetry (PIV) experiment, the combination experiment of oil flow, and the wire experiment of the wavy leading edge wing flow field, respectively. The results show that as the angles of attack increase, the fluid separation occurs in advance at the trough section of the wing, and the flow separation is delayed at the peak section. Gross and Fasel (2013), Lin et al. (2013), Prez-Torr and Kim (2017) and Skillen and Favier (2015) observed similar separation flow characteristics through LES and direct numerical simulation (DNS) methods respectively. Stanway (2008) observed a series of counterrotating streamwise vortices induced by the wavy tubercles in water tunnel experiment research. Hansen et al. (2011Hansen et al. ( , 2016 further studied the shape and development of the similar streamwise vortices in the wind tunnel experiment through the hydrogen bubble visualization method and PIV. The development and dissipation process of streamwise vortices in the downstream were studied using PIV by Wei et al. (2015). Arai et al. (2010), Carreira and Kobayashi (2008) and Rostamzadeh et al. (2014) captured the streamwise vortices in their numerical investigations. The flow field around the spanwise wavy wings and the wavy leading edge wings will show periodic and non-periodic flow fields under certain conditions. Custodio (2007) first described the periodic phenomenon of flow around a wavy leading edge wing in the dye experiment. Attached flow and separation vortices appear alternately in adjacent troughs on the upper surface of the wing. This phenomenon has been verified in the simulation results of Cai et al. (2015), Melo and Camara (2013) and Zhao et al. (2017). Custodio (2007) found in the wire experiment that when the wing has a wingspan of eight times the leading edge wavelength, the flow presents a more complex non-periodic phenomenon. Local stalled vortices can be observed at some trough positions, while at other trough positions the flow remains unstalled and is accompanied by a complicated spanwise flow. Wei et al. (2015) used PIV to study the development and evolution of the streamwise vortices of the wing with 16 waves and found that the flow vortex pairs are relatively stable at a small angle of attack, and the adjacent vortex pairs gradually merge and dissipate with the increase of the angle of attack, showing non-periodic characteristics in space. In terms of numerical calculations, the non-periodic flow phenomena was found by Dropkin et al. (2012) and Rostamzadeh et al. (2017) in their investigations based on unsteady RANS.
Due to the complex flow field structure and flow phenomena, there is still no consensus on the flow control mechanism of the bionic modified wing. Miklosovic et al. (2007) believed that the tubercles played a role similar to the vortex generators in delaying the occurrence of stalls by enhancing the momentum exchange between the boundary layer and the main flow. The flow around wings with spanwise waviness were investigated by Serson and Meneghini (2015); Serson et al. (2017) for Reynolds numbers of 1.0 × 10 4 and 5.0 × 10 4 . It was believed that the spanwise flow induced by the tubercles of the leading edge weakened the suction peak of the wave crests, so as to maintain fluid attached. Van et al. (2007) believed that the downwash velocity of the streamwise vortex pairs reduced the effective angle of attack of the wave crest section and delayed the flow separation at high angles of attack. Cai et al. (2017) studied the influence of single and multiple tubercles on the control of stall characteristics through wind tunnel tests and numerical methods and indicated that the partition effect of high-momentum attached flow and the downwash from the local stall sections are the control mechanisms of the tubercles. Prez-Torr and Kim (2017) believed that the aerodynamic characteristics of the wavy leading edge were attributed to the reduced spanwise coherence of the flow, which weakened the strength of the Karman vortex shedding. Hansen et al. (2011) believed that the tubercles of the leading edge changed the pressure distribution along the span of the wing, so that the crests had a higher pressure and a lower pressure gradient, thus delaying the flow separation.
In recent years, the reduced-order model (ROM) of the unsteady flow field has been widely used in the study of the characteristics of the unsteady flow field. Among them, the proper orthogonal decomposition (POD) and the dynamic mode decomposition (DMD) are two typical modal decomposition methods. Chen et al. (2010) and Poplingher and Raveh (2019) used POD and DMD methods to study the flow field modes and shock wave motion characteristics in transonic buffeting. Kou and Zhang (2017) used the DMD method to analyze the characteristic modes of the flow around a cylinder. The flow characteristics of wavy leading edge wing in static and dynamic stalls was investigated by Zhao et al. (2020Zhao et al. ( , 2021. The results showed that compared with the POD method, the DMD method could decompose the flow field according to the frequency characteristics, which is beneficial to the dynamic analysis and research of the complex flow field. When a plane is in a stall condition, the region of flow separation covers the entire upper surface, and the aerodynamic force changes over time violently. Wings are prone to flutter, which can cause structural fatigue or even damage. Some progress has been made in active flow control (Kim & Kim, 2021), but it is relatively difficult to implement in practice due to its complex structure. Passive flow control based on geometry modifications may have brighter application foreground. Although several numerical and experimental studies on the bionic wavy modified wings were carried out, the corresponding flow control effects and mechanisms of different wavelengths are rarely get attention, especially for high Reynolds number conditions (Re > 10 6 ). The investigation of wavelength for spanwise wavy wings in compressible flow with high Reynolds number is necessary to be conducted. In this paper, the numerical simulations on spanwise wavy wings with different wavelengths at Re = 10 7 were carried out, and the aerodynamic characteristics, flow field structure characteristics, and turbulence characteristics were studied to reveal the flow control effects and mechanisms. The three components of the three-dimensional flow field velocity vector have been modally decomposed by the DMD method, and the spatial distribution, frequency characteristics, and vortex structures of the modes are analyzed.

Scale adaptive simulation (SAS)
Although, large eddy simulation (LES) is able to simulate the separated flow accurately (Zahiri & Roohi, 2021). It is rarely used to solve high Reynolds number problems for the higher computational cost. The SAS concept is based on the introduction of the von Karman length scale into the turbulence scale equation. The information provided by the von Karman length scale allows SAS models to dynamically adjust to resolved structures in an unsteady Reynolds-Averaged Navier-Stokes (URANS) simulation resulting in an LES-like behavior in unsteady regions. Simultaneously, the model provides standard RANS capabilities in stable flow regions (Younsi et al., 2008). The SAS approach could acquire the accurate prediction for unsteady separated flow with the lower cost. Therefore Spalart-Allmaras (SA) one-equation turbulence model is utilized to establish the SAS approach.
In the Cartesian coordinate system, the three-dimensional compressible Favre-averaged Navier-Stokes equations include mass, momentum, and energy equations: Here:'−' and '∼ ' represent the direct and densityweighted filtering operations, respectively.ρ,ũ i , t,p,Ẽ denote the density, velocity component, time, pressure and total energy, respectively. The viscous stress tensor τ ij and heat fluxq i can be written as Here, μ is the molecular viscosity coefficient, calculated by Sutherlands law Sutherland (1893); γ is the specific heat ratio of the air;T is the temperature; Pr is the molecular Prandtl number; the strain rate tensorS ij = 0.5(∂ũ i /∂x j + ∂ũ j /∂x i ). The total energyẼ can be given by the relationship: In Equations (2) and (3), the unclosed terms τ ij,t and q i,t can be processed by the turbulent eddy viscosity hypothesis: Here, μ T is the turbulent viscosity coefficient; Pr T is the turbulence Prandtl number, which is equal to 0.92 in SAS.
In the practical operation, the wall distance d can be modified by Xu et al. (2019): Here, κ is the von Karman constant (κ = 0.41). The von Karman length scales can be written as The numerical code based on the finite volume method is used in the investigation. Implicit approximate factorization method and sub-iteration technique are used for temporal advancement to achieve the second-order accuracy. The traditional second-order central scheme is used to discretize the viscous terms. The second-order central/upwind hybrid scheme is proposed to discretize the convection terms. The code has been verified in our previous investigations (Xu & Ni, 2018;Xu et al., 2019). Numerical results of axisymmetric base, circular cylinder, and supersonic boundary layer are agreed with the experimental data.
In this paper, the experimental results of flow past a cylinder (Rodriguez, 1984) are used to verify the credibility of the grids for compressible separation flows. The spanwise length is three times the diameter of the cylinder D, the incoming Mach number is 0.55, and the Reynolds number based on the diameter of cylinder is Re = 2.56 × 10 5 . The number of grid points in the azimuthal, radial, and spanwise direction are 481, 241, and 121. The first off-wall grid spacing in viscous wall units is about 1. At least 30 grid points are in the boundary layer region, where the grid sizes in viscous wall units are X + max < 80, Y + max < 50, Z + max < 200. The both sides of the model are the periodic boundaries and the remaining boundaries are the pressure farfields. The time step is 0.005D/a ∞ ( a ∞ is the local speed of sound), the number of sub-iteration steps is 20. The CFL number is 1. The flow field information sampled after the flow field is fully evolved with 30 sampling intervals, and the sampling time length is more than 150c/a ∞ . The computational meshes are shown in Figure 1.
The Strouhal number ) based on the simulated time series of the lift coefficient is 0.20, which is consistent with experimental results. Figure 2 shows that the time-averaged dimensionless pressure on the surface of the numerical results is agreed with experimental results, indicates that the calculation method has sufficient credibility.

Dynamic mode decomposition (DMD)
The DMD method is a Reduced-Order Model (ROM) proposed by Schmid and Sesterhenn (2010) based on characteristic frequency. It could be used to process highdimensional sequential measurements, extracts coherent  structures, isolates dynamic behavior, and reduces complex evolution processes to their dominant features and essential components.
The specific method is to obtain snapshots of the flow field from time 1 to N through experiments and numerical methods, and divide them into two matrices of 1 ∼ N-1 and 2 ∼ N, as shown in Equations (11) and (12). Suppose there is a linear mapping A such that Y = AX: Carry out singular value decomposition on matrix X: A is the optimal low-dimensional approximate matrix of A, The eigenvalue ofÃ is part of A, that is, the Ritz eigenvalue. Calculate the linear map A and its eigenvectors ω and eigenvalues λ [ω, λ] = eig(A) Define the dynamic modes as The growth rate and frequency of mode j are defined as Define the energy of the mode as the sum of the squares of the signal amplitudes at all moments and points. The signal at a certain moment is equal to the initial signal multiplied by the characteristic value, which can be calculated as

Data reduction
The Reynolds number Re c based on the chord length is expressed as follows: where ρ ∞ is the density of the incoming flow, U ∞ is the velocity at the farfield, μ ∞ is the flow dynamic viscosity. The Mach number of incoming flow Ma ∞ is written as where a is the velocity of sound.
The lift coefficient (C L ) and drag coefficient (C D ) can be defined as Here, F L and F D are the lift and drag, respectively. The reference area is defined as S = 1. The pressure coefficient (C P ) is written as Here, P is the local pressure and P ∞ is the pressure of farfield.

Parameterization of the wavy wing geometry
The wings with spanwise waviness can be regarded as the combination of the other two types of wings. NACA0012 with no trailing edge thickness is used for simulation at M ∞ = 0.5 and Re = 10 7 based on chord length (c = 1 m). It is a typical subsonic airfoil commonly used in flow control investigation. The geometrical models are obtained by stretching NACA0012 along the spanwise direction of a cosine curve, and the cosine curve is as in Equation (25). The aspect ratio is 2, and the flight angle of attack is set to 12 deg. The wavelength of the curves are selected as λ = c/2, c/3, c/4, respectively. The ratio of wave height h to chord length is 0.1, and cosine waves of 4, 6, and 8 periods are distributed along the span of the wings, as shown in Figure 3. The parameters of mesh (seen in Figure 4) and time steps used for the calculation are the same as for the validation case.

Force behaviors
Compared with W-4 and W-6, the time-averaged lift-todrag ratio of W-8 increases by 12.7% and 14.4% (seen in Table 1). The standard deviations of C L and C D for W-8 are significantly higher than that for W-4 and W-6. Figure 5 shows the change of C L and C D with time for different wavelength wings. The C L and C D curves of the three wings show oscillations of different magnitudes with time, indicating that the flow fields exhibit unsteady properties. The aerodynamic of W-8 oscillates more strongly than W-4 and W-6. The lift power spectrum density for different wavelength wings is obtained through Fourier transform, as  Considering that the aerodynamic performance of wing is closely related to the recirculation regions in the flow field, Figure 7 shows the distribution of time-averaged recirculation regions for different wavelength wings. The recirculation regions are shown in a     (2017) and Rostamzadeh et al. (2017). The flow field of W-6 exhibits periodic characteristics of three times the wavelength, and a pair of separated bubbles skewed to the center wave crest is formed between the two dilation-type recirculation regions. The size of recirculation regions on the upper surface of W-8 is much lower. The distribution of these regions exhibits irregularities due to the inter-transition of separation bubbles and expansion recirculation regions. Figure 8 shows the contour of pressure coefficient and limiting streamlines. The separation occurred at the troughs, and a relatively high negative pressure is observed in the leading separated bubbles of W-4 and W-6 wings. It is mainly due to the flow being deflected by the tubercles at the leading edge and moving around it.  . The low pressure regions of multiple troughs at the W-8 leading edge communicate with each other and are significantly larger than that of W-4 and W-6. Only one dilation-type recirculation region is formed in the corresponding span, and thus the W-8 wing can maintain high lift at AOA = 12 deg. The vortex core locations in the recirculation regions move upward as the wave number goes up with shorter wavelength.

Flow structures
To better understand the control mechanism of wings with spanwise waviness, the three-dimensional vortex structures of the time-averaged flow fields are displayed by the Q criterion, as shown in Figure 9 (Q = 0.5( 2 − S 2 ), Q = 0.5). It can be seen that the deflected bubbles of the W-6 wing fall off downstream and evolve to form a streamwise counter-rotating vortex pairs (CVPs), which are considered to be an important reason for the separation control of the bionic wavy modified wings by Hansen et al. (2016). For the W-4 flow field, the separated bubbles did not evolve to form the CVPs, but achieved a flow control effect similar to that of W-6, indicating that CVPs may only appear in the flow field of a wavy modified wing at certain specific wavelengths, and does not constitute the main mechanism of the flow control effect. The flow field structure  of W-8 is more complex. The flow vortex pairs merge and dissipate downstream of the separation bubble, and the scale of the spanwise vortex structure is between the separation bubble and the dilation-type recirculation region. Figure 10 shows the spanwise pressure gradient distributions and streamlines (the transparent surface is the isosurface V x = −0.001V ∞ ). The flow field structure of the leading separated bubbles and the dilationtype recirculation regions changes the pressure gradient of the leeward side surface, drives the spanwise flow of the airflow. The momentum transport generated by the spanwise flow is one of the reasons to keep the airflow reattached downstream of the leading separated bubbles. As W-8 has a greater waveheight-wavelength ratio, the spanwise pressure gradient of the upper surface is more densely alternated between positive and negative, which will increase the spanwise disturbance of the airflow, resulting in a more complicated flow field structure of the W-8 wing. Analyzing the vortex dynamics characteristics of the flow field can further understand the flow control mechanism of wings with spanwise waviness, and the surface is one of the important sources of vorticity. Boundary vorticity flux (BVF) represents the rate at which the vorticity is generated on the non-slip wall and diffuses into the external flow, and it can also be used to represent the pressure bulging effects and shear effects of the fluid on the wall. For high Reynolds number conditions, BVF can be approximated as Equation 26 (Wu et al., 2007): Figure 11 shows that the streamwise component of BVF is higher in the leading edge of both sides of the crest and corresponding downstream regions, indicating that the tubercles have a strong pressure bulging effect and shear effect on the fluid. Comparing the three types of wings, the streamwise BVF of W-4 is much smaller, and the streamwise BVF of the dashed areas on W-6 and W-8 has several higher peaks, indicating that the streamwise vorticity is generated from the wall and diffuses into the flow field rapidly. The dashed areas in Figure 11 are consistent with the position where the CVPs generate on the W-6 and W-8 wings, as we can see from Figure 9. The Lamb vector divergence distributions reflects the transport of momentum and vorticity in the flow field (Hamman et al., 2008). The generalized Lamb vector (L = ω × V − T∇s) serves as the most commonly encountered mechanism for the coupling between the compressing and shearing processes in the interior of the flow field. The appearance of T∇s implies that both shearing and compressing processes are coupled with thermodynamic process (Wu et al., 2007). For moderate M ∞ , like M ∞ = 0.5 considered here, entropy gradients are small as the solution is smooth. The Lamb vector divergence is, therefore, approximately given by ∇ · L = V · ∇ × ω − ω · ω. ∇ · L > 0 represents that the straining motion caused by the flexion product is dominant; ∇ · L < 0 represents that the vorticity-bearing motion caused by the entropy is dominant; the positive and negative values alternately represent the change of the momentum transport mechanism in the flow field (Xu et al., 2020). Figure 12 shows the distributions of Lamb vector divergence on the characteristic section. SLACE A, D, and H are the chordal cross-sections of the dilationtype recirculation regions of the three types of wings. They are a two-layer structure with positive and negative arrangements, corresponding to the existence of different momentum transport mechanisms on each sides of the shear layer: positive values correspond to the strong strain rate side is dominated by straining motion, and a negative value corresponds to the side of the vorticity zone, dominated by vorticity bearing motion. SLICE C and F are a three-layer structure with a positive and negative arrangement in the cross section of the separated bubbles. Different from the shear layer in the dilationtype recirculation regions, the separated bubbles have strong straining motion for fluid near the wall. SLACE B and E are the attached flow areas between the separated bubbles and the dilation-type recirculation regions, and ∇ · L are positive along the flow direction, indicating that the momentum transport in this areas is dominated by straining motion, which is consistent with the phenomenon shown in Figure 9. SLACE G and I are the attached flow cross-sections between two separated bubbles, showing a weak positive and negative two-layer structure.
Take the normal section of the incoming flow at the half chord (x = 0.5c) and observe the z-direction velocity distribution of the airflow, as shown in Figure 13. Van et al. (2007) believed that the recirculation regions reduced the effective angle of attack at the local wing. The change in the effective angle of attack caused the upwash and downwash of the airflow. The segment with a smaller effective angle of attack has a lower circulation, thereby inducing a strong downwash airflow in the section with a higher effective angle of attack and circulation, as shown in Figure 13. Therefore, the upwash airflow is formed in the dilation-type recirculation regions, and the downwash airflow is formed downstream of the separation bubbles, so that the airflow reattaches downstream of the separation bubble. W-4 and W-6 wings upwash flow field and downwash flow field have similar span area distribution. The W-8's downwash area is significantly larger than the first two types of wings, so the W-8 has smaller recirculation regions and a higher lift-to-drag ratio.

DMD-based characterization of the Favre-averaged model resolved unsteady flow component
To further study the frequency characteristics and modal characteristics of the wings flow field, 600 fully evolved continuous flow snapshots were selected of the W-6 wing, and the three components of the three-dimensional flow field velocity vector have been modally decomposed by the DMD method. Figure 14 shows the distribution of the real and imaginary parts of the modal eigenvalues. The points in the unit circle indicate the stable modes, on the unit circle indicate the periodic modes, and outside the unit circle indicate the divergent modes. It can be seen from the figure that most of the modes decomposed by DMD are on the unit circle, which are periodic modes, indicating that the flow field is in a critically stable state. Figure 15 shows the energy spectra obtained by DMD. Red marks the top five modes sorted by energy, as well as the four main high-frequency modes. M1 is the firstorder mode (St M1 = 0) of the W-6 wing flow field, which has the most energy and represents the time-averaged flow field. Figure 16 shows the vortex structure of M1 mode through the isosurface Q = 1 × 10 −6 . It can be seen that the leading edge separation vortex (LSV) and the dilation-type recirculation regions are alternately distributed at the positions of the troughs. The airflow  downstream of the LSV is induced to form the counterrotating vortex pairs (CVPs). Due to the speed difference on both sides, the shear layer instability occurs at the front edge of the dilation-type recirculation regions, which induces the generation of shear-layer vortex (SLV), and the SLV moves downstream and merges to form the Arch vortex (AV). The symmetrical secondary vortex pairs (SVPs) are formed under the induced action of the AV, the shearing effect of the airflow on the lower surface and the induction of the trailing edge protrusion.
The turbulent characteristics of the flow field can be revealed by the pressure pulsation information of the detection point. Six points were selected to set the pressure signal probes Point 1 to Point 6, marked in Figure 17, and power spectral density of the pressure signals was analyzed. Figure 18 shows the power spectral density profiles of pressure signals at six points. Different slope relationships reveal the propagation characteristics of the solid wall pressure disturbance in the flow field. Point 1 reflects the characteristics of the pressure signal at the shear layer of the dilation-type recirculation region, indicates the characteristic frequency St ∼ 1.4, which is associated with the SLV motions. Point 2 and Point 3 reflect the pressure signal characteristics of the SLV and CVPs with the characteristic frequency St ≈ 1.16 ∼ 1.36 and its harmonics. Point 4 is located at the SV of the trailing edge with the characteristic frequency St ≈ 0.33; Point 5 is located at the side of the AV with St ≈ 0.33, and a high frequency signal of St ≈ 1.33 appears due to the disturbance of the CVPs. Point 6 is located at the top of the AV with St ≈ 0.4, therefore the characteristic frequency of the AV downward shedding movement is St ≈ 0.33 ∼ 0.4.
At Points 1-6, the lower-frequency regime has a slope relationship of St 0.4 , which is related to the turbulence structure in the outer region of the boundary layer (Panton & Linebarger, 1974); an St −1 scaling can be observed at the intermediate frequencies of Points 1, 4-6, related   to the vortex structure in the logarithm law region of the boundary layer; the St −5/3 scaling at Point 1 reveals that the calculation results seem to reach the inertial subrange scales (Na & Moin, 1998). At points 2 and 3, the spectrums with St −0.7 and St −2 scaling are presented. Tsuji and Fransson (2007) assumed that the contribution to the 0.7 power law is seen to originate from positions up to the end of the log region. The slope of −2 indicates that the wall-pressure field has not yet recovered from the separation (Na & Moin, 1998). The power spectral density profiles of pressure pulsation at Points 2-6 all show a scaling of St −3 . Simpson et al. (1987) pointed out that the pressure in the adverse pressure gradient boundary layer and the separated turbulent boundary layer will have the −3 power slope attenuation relationship.
It can be seen from Figure 19 that St M2 = 0, indicating that M2 is a steady mode; St M3 = 1.72 is similar to the characteristic St number of the shear layer, which may be related to the SLV motions. In the M4 and M5 diagrams, there are several arched vortices arranged in an orderly manner downstream, indicating that the M4 and M5 modes are related to the shedding evolution of the AV in the flow field, and St M4 = 0.32 is consistent with the characteristic St number of Point 4. Figure 20 shows the high-frequency modes occupying the main energy. From the frequency characteristics and the vortex structures of the modes, it can be seen that the high-frequency modes characterize the small-scale eddy motions in the flow field, and their characteristic St numbers are related to the scale of the vortex structures. The modes of the flow field are more meticulous for the higher St number. It also shows that the DMD method can obtain different scale vortex motion modes of the flow field based on the characteristic frequency.

Conclusion
In this paper, the SAS method is used to numerically investigate the compressible flow around wings with spanwise waviness for an incoming flow with M ∞ = 0.5 and Re = 10 7 . Various fundamental mechanisms dictating the complex flow behaviors are summarized as follows.
The leading edge separated bubbles and the dilationtype recirculation regions are alternately distributed at the upper surface of the wings, which changes the pressure gradient of the flow field. The pressure gradient drives the spanwise flow of the airflow and generates spanwise momentum transport. At the same time, the upwash airflow is formed in the dilation-type recirculation regions, and the downwash airflow is formed downstream of the separation bubbles. The spanwise momentum transport of the airflow and the downwash effect of the partial wing section keep the airflow reattached downstream of the separated bubbles. Spanwise momentum transportation and downwash effect play important roles in the flow control mechanism of the wings with spanwise waviness.
The low-pressure regions of multiple troughs at the W-8 leading edge communicate with each other and are significantly larger than that of the W-4 and W-6 wings; at the same time, the airflow downwash area increases to make the flow separation on upper surface reduced. Therefore, compared with W-4 and W-6, the lift-to-drag ratio of W-8 increases by 12.7% and 14.4%, and its lift and drag coefficient curves oscillate more strongly.
The three components of the three-dimensional flow field velocity vector have been modally decomposed by the DMD method. The vortex structures of the flow field based on different frequencies are presented and analyzed. In the power spectral density profiles of W-6 pressure probes, the separation bubble and its downstream show a relationship of St 0.4 , St −0.7 , St −2 , and St −3 scalings, and the dilation-type recirculation region shows St 0.4 , St −1 , and St −3 scalings. The main modes characterizing the motions of AV and SLV have been obtained through the DMD method. The St numbers of the modes are related to the scale of the vortex structures in the flow field. Small-scale eddy motions have been decomposed in high-frequency modes, and thus the modes of the flow field are more meticulous for the higher St number. It also shows that the DMD method could obtain different scale vortex motion modes of the flow field based on the characteristic frequency.
Although wavy modification makes part of the flow attached to the wing in stall condition. It is still difficult to completely eliminate flow separation. The combination of active flow control or optimization with wavy modification is worthy of further investigation.

Disclosure statement
No potential conflict of interest was reported by the author(s).