Towards the Understanding of Humpback Whale Tubercles: Linear Stability Analysis of a Wavy Flat Plate

: The results from a temporal linear stability analysis of a subsonic boundary layer over a ﬂat plate with a straight and wavy leading edge are presented in this paper for a swept and un-swept plate. For the wavy leading-edge case, an extensive study on the e ﬀ ects of the amplitude and wavelength of the waviness was performed. Our results show that the wavy leading edge increases the critical Reynolds number for both swept and un-swept plates. For the un-swept plate, increasing the leading-edge amplitude increased the critical Reynolds number, while changing the leading-edge wavelength had no e ﬀ ect on the mean ﬂow and hence the ﬂow stability. For the swept plate, a local analysis at the leading-edge peak showed that increasing the leading-edge amplitude increased the critical Reynolds number asymptotically, while the leading-edge wavelength required optimization. A global analysis was subsequently performed across the span of the swept plate, where smaller leading-edge wavelengths produced relatively constant critical Reynolds number proﬁles that were larger than those of the straight leading edge, while larger leading-edge wavelengths produced oscillating critical Reynolds number proﬁles. It was also found that the most ampliﬁed wavenumber was not a ﬀ ected by the wavy leading-edge geometry and hence independent of the waviness.


Introduction
The aerodynamic, hydrodynamic, and acoustic performance of an airfoil with a wavy leading edge (LE) has been a subject of research for the past 30 years. This airfoil geometry was inspired by the pectoral fin of a humpback whale (Megaptera novaeangliae), which features large LE protuberances that are presumed to aid in maneuverability. One feeding behavior of a humpback whale involves swimming in circles of relatively tight radii to create a bubble net, corralling prey for efficient consumption. This banking maneuver is efficiently performed if the pectoral fin has a large stall angle combined with increased lift and decreased drag to reduce energy expenditure. Each of these performance characteristics has been investigated in the literature. Early studies of a humpback whale pectoral fin with LE protuberances were performed by Bushnell and Moore [1], Fish and Battle [2], and Miklosovic et al. [3], which suggest that the protuberances improve hydrodynamic maneuverability by increasing the stall angle and flattening out the lift curve, at the expense of reducing the maximum lift. Drag reduction was also observed due to LE tubercules. Parametric experiments by Johari et al. [4] and aerodynamic modeling by van Nierop et al. [5] were performed for several sinusoidal LE amplitudes from 2.5% to 12% chord and two wavelengths of 25% and 50% chord. These studies showed that increasing the LE amplitude increases the stall angle, while changing the LE

Geometry
The geometry for the general case of a swept flat plate is shown in Figure 1. For the swept plate with crossflow velocities, two coordinate systems are used, one a body fixed system in the chord and span directions xc and yc, and a second a rotated coordinate system xs and ys, shown rotated in line with the local streamline. For the simplified case of zero sweep, the velocities are only in the chord direction and the rotated coordinates are not used. The geometry for a swept flat plate with a wavy leading edge is shown in Figure 2. The plate chord is set to c = 0.5 m to match the DLR experimental conditions, with the amplitude, A, and wavelength, λ, of the sinusoidal profile defined as a percentage of the chord. The peak of the LE wavelength is located at λ, and the trough is at 1.5 λ. The stability results for the swept flat plate are presented in terms of the 99.9% boundary layer thickness Reynolds number. The results for the zero-sweep case are presented in the nondimensional variables typically used in the Blasius equation. To ensure an equivalent wetted surface area of the flat plate with straight and wavy LEs, the chord for the wavy LE was extended by a length equal to one half of the LE amplitude A. All the chord fractions for the wavy LE plate were normalized using a chord equal to c + A/2.

Geometry
The geometry for the general case of a swept flat plate is shown in Figure 1. For the swept plate with crossflow velocities, two coordinate systems are used, one a body fixed system in the chord and span directions xc and yc, and a second a rotated coordinate system xs and ys, shown rotated in line with the local streamline. For the simplified case of zero sweep, the velocities are only in the chord direction and the rotated coordinates are not used. The geometry for a swept flat plate with a wavy leading edge is shown in Figure 2. The plate chord is set to c = 0.5 m to match the DLR experimental conditions, with the amplitude, A, and wavelength, λ, of the sinusoidal profile defined as a percentage of the chord. The peak of the LE wavelength is located at λ, and the trough is at 1.5 λ. The stability results for the swept flat plate are presented in terms of the 99.9% boundary layer thickness Reynolds number. The results for the zero-sweep case are presented in the nondimensional variables typically used in the Blasius equation. To ensure an equivalent wetted surface area of the flat plate with straight and wavy LEs, the chord for the wavy LE was extended by a length equal to one half of the LE amplitude A. All the chord fractions for the wavy LE plate were normalized using a chord equal to c + A/2.

Linear Stability Equations
The solution method for the three-dimensional linear stability equations follows Mack [29]. The nonlinear perturbation equations were derived from the nondimensional incompressible Navier-Stokes equations, ∂u ∂x + ∂v ∂y + ∂w by introducing the perturbed velocities and pressures, u = U + u , v = V + v , w = W + w and p = P + p , and subtracting the original non-perturbed equation. Here, u, v, w are the fluid velocities in the x, y, z coordinates, respectively, with the z coordinate across the shear layer. The Reynolds number is R; p is the pressure; U, V, W, P are the local base flow velocities and pressure; and u , v , w , p are the perturbation velocities and pressure. The boundary layer thickness and the local freestream velocity were used as the reference length and velocity scales, respectively. By assuming small disturbances and a locally parallel base flow, the equations were simplified to a linear and separable system of partial differential equations. This allows the use of a separable function for the disturbances, where e, φ, h, π are the disturbance amplitudes corresponding to u , v , w , p , respectively; α and β are disturbance wavenumbers in the x and y directions, respectively; and ω is the frequency. For the temporal analysis, the frequency is a complex number, where the imaginary component is the amplification factor. The amplification factors can be positive, zero, or negative, indicating an unstable, neutrally stable, or stable disturbance, respectively. Substituting (2) into the linearized, separable perturbation equations gives the nondimensional continuity and momentum perturbation equations, now a system of ordinary differential equations: In Equation (3), the prime, , and double prime, ", indicate first and second derivatives with respect to z, respectively. The boundary conditions for Equation (3) are no slip at the wall, and that the disturbances go to zero in the freestream, The shooting method was used to satisfy the flat plate boundary condition and solve iteratively for the eigenvalues. The compound matrix method of Ng and Reid [26,43] was used to relieve the stiffness of the equations.

Base Flow Models
Solutions of the Blasius equation, the FSC equations, and the three-dimensional incompressible BLE were used as the input to the stability code. Similarity solutions of the Blasius and FSC equations for a straight LE plate were used for the validation of the stability code for swept and un-swept plates.
The FSC equations can be used to approximate an infinite span swept flat plate with a favorable pressure gradient and straight LE, with the pressure gradient controlled by the dimensionless Hartree parameter β H and the three-dimensionality controlled by the local sweep angle θ. The chordwise and spanwise components of the FSC model flow are and the functions f and g are the solutions to the following FSC equations: with boundary conditions In the above equations, U 0c and V 0c are the chord and span velocities, respectively. It is noted that setting β H = θ = 0 results in the Blasius equation. The FSC equations were solved using the shooting method, where the freestream boundary conditions were met by iterating on the unspecified initial conditions f (0) and g (0) to drive f (∞) and g(∞) to unity.
The boundary layer approximation to the incompressible Navier-Stokes equations for steady flow is given in Equation (9). These are parabolic equations and can be marched in the chord and span directions with the known free stream velocities U and V and specified initial conditions. For an infinite span flat plate, the flow is three-dimensional, but its velocities are independent of the span-wise coordinate y; hence, the bracketed terms in Equation (9) are neglected.
Equation (9) is subject to the following boundary conditions: For the wavy LE plate, free boundary conditions in the plane of the plate surface were applied when the flow was in between LE peaks and not contacting the plate surface. The BLE were solved using an implicit finite difference discretization and uniform grid, with the tri-diagonal matrix algorithm for matrix inversion [44].

Results
This section begins with the validation of the base flow and stability calculations by comparison with published numerical solutions and experimental data. After the base flow, the stability results for Fluids 2020, 5, 212 7 of 23 a zero-sweep and swept flat plate are presented. Most of our results are presented in nondimensional variables unless otherwise specified. For the wavy leading-edge plate, a uniform notation is adopted similar to that found in the literature [8][9][10] to indicate the wavelength and amplitude of the waviness. For example, "aAwW" means a% chord Amplitude and w% chord Wavelength with a and w varying between 0 and 40 in this paper. As an example, a 10A30W wavy leading-edge plate indicates a 10% chord amplitude and a 30% chord wavelength. A straight leading edge is represented by "00AwW", indicating a 0% amplitude and any wavelength. In addition to these wavy leading-edge notations, the notation "aAwW lL" is also used to indicate the spanwise integration extent. Namely, "10A30W 2L" refers to 10% chord amplitude and 30% chord wavelength integrated to two wavelengths.

Base Flow Validation
The FSC equations, infinite span BLE, and full 3-D BLE integrations were first solved for the case of a zero pressure gradient and a 45 • swept plate, with all the results matching published values for the Blasius profile. Next, the FSC and infinite span BLE were compared to experimental data from the DLR [33] for a swept plate with a strong favorable pressure gradient, at an 80% chord location. An effective sweep angle of 42.5 • was used in the calculations based on the DLR data. For the BLE code, the spanwise flow component of the freestream velocity was set to a constant 13 m/s, and the freestream velocity component in the chord direction increased linearly from 4.9 to 19.1 m/s to match the measured velocity gradient. Fischer and Dallmann [45] provide values for the FSC parameters from 10% to 90% chord, which were used to set up the FSC code. The infinite span BLE code requires an initial velocity profile to begin the marching procedure. The initial conditions at 10% chord were provided from the FSC integrator and marched forward in the BLE code to 80% chord. Figures 3  and 4 show the comparison of the experimental and calculated velocity profiles aligned with the local freestream vector. Good agreement was found between the FSC calculations of Fischer and Dallmann [45] and our results, and between the BLE integrations and experimental data. between 0 and 40 in this paper. As an example, a 10A30W wavy leading-edge plate indicates a 10% chord amplitude and a 30% chord wavelength. A straight leading edge is represented by "00AwW", indicating a 0% amplitude and any wavelength. In addition to these wavy leading-edge notations, the notation "aAwW lL" is also used to indicate the spanwise integration extent. Namely, "10A30W 2L" refers to 10% chord amplitude and 30% chord wavelength integrated to two wavelengths.

Base Flow Validation
The FSC equations, infinite span BLE, and full 3-D BLE integrations were first solved for the case of a zero pressure gradient and a 45° swept plate, with all the results matching published values for the Blasius profile. Next, the FSC and infinite span BLE were compared to experimental data from the DLR [33] for a swept plate with a strong favorable pressure gradient, at an 80% chord location. An effective sweep angle of 42.5° was used in the calculations based on the DLR data. For the BLE code, the spanwise flow component of the freestream velocity was set to a constant 13 m/s, and the freestream velocity component in the chord direction increased linearly from 4.9 to 19.1 m/s to match the measured velocity gradient. Fischer and Dallmann [45] provide values for the FSC parameters from 10% to 90% chord, which were used to set up the FSC code. The infinite span BLE code requires an initial velocity profile to begin the marching procedure. The initial conditions at 10% chord were provided from the FSC integrator and marched forward in the BLE code to 80% chord. Figures 3 and  4 show the comparison of the experimental and calculated velocity profiles aligned with the local freestream vector. Good agreement was found between the FSC calculations of Fischer and Dallmann [45] and our results, and between the BLE integrations and experimental data.  Next, the full 3-D BLE were compared for an infinite span flat plate with wavy and straight LEs under the DLR freestream conditions. To properly compare the infinite span straight LE and wavy LE boundary layers, the BLE for the wavy LE geometry were integrated spanwise across successive LE peaks until the flow profile at each LE peak did not change with further spanwise integration. The infinite span boundary layer profile was used as the initial condition at the wavy LE peak to start the integration. Figure 5 shows the crossflow boundary layer profiles for a wavy LE of 10% chord amplitude and 7.5% chord wavelength (10A075W) at different span locations. It is seen that the flow develops from the infinite span straight LE profile until it reaches a limiting profile at four wavelengths (10A075W 4L). Subsequent integrations to 5L and 6L show no change in the base flow profile. Furthermore, the same grid was used to integrate the infinite span solution spanwise across six wavelengths for a straight LE (labeled 00A075W 6L), and it is seen that the integration preserves the initial condition velocity profile, indicating no significant discretization errors. A similar result is found in Figure 6 for the wavy LE with a 10% chord amplitude and 30% chord wavelength (10A30W), except that the boundary layer reaches the limiting profile at one LE wavelength.  Next, the full 3-D BLE were compared for an infinite span flat plate with wavy and straight LEs under the DLR freestream conditions. To properly compare the infinite span straight LE and wavy LE boundary layers, the BLE for the wavy LE geometry were integrated spanwise across successive LE peaks until the flow profile at each LE peak did not change with further spanwise integration. The infinite span boundary layer profile was used as the initial condition at the wavy LE peak to start the integration. Figure 5 shows the crossflow boundary layer profiles for a wavy LE of 10% chord amplitude and 7.5% chord wavelength (10A075W) at different span locations. It is seen that the flow develops from the infinite span straight LE profile until it reaches a limiting profile at four wavelengths (10A075W 4L). Subsequent integrations to 5L and 6L show no change in the base flow profile. Furthermore, the same grid was used to integrate the infinite span solution spanwise across six wavelengths for a straight LE (labeled 00A075W 6L), and it is seen that the integration preserves the initial condition velocity profile, indicating no significant discretization errors. A similar result is found in Figure 6 for the wavy LE with a 10% chord amplitude and 30% chord wavelength (10A30W), except that the boundary layer reaches the limiting profile at one LE wavelength.   Next, the full 3-D BLE were compared for an infinite span flat plate with wavy and straight LEs under the DLR freestream conditions. To properly compare the infinite span straight LE and wavy LE boundary layers, the BLE for the wavy LE geometry were integrated spanwise across successive LE peaks until the flow profile at each LE peak did not change with further spanwise integration. The infinite span boundary layer profile was used as the initial condition at the wavy LE peak to start the integration. Figure 5 shows the crossflow boundary layer profiles for a wavy LE of 10% chord amplitude and 7.5% chord wavelength (10A075W) at different span locations. It is seen that the flow develops from the infinite span straight LE profile until it reaches a limiting profile at four wavelengths (10A075W 4L). Subsequent integrations to 5L and 6L show no change in the base flow profile. Furthermore, the same grid was used to integrate the infinite span solution spanwise across six wavelengths for a straight LE (labeled 00A075W 6L), and it is seen that the integration preserves the initial condition velocity profile, indicating no significant discretization errors. A similar result is found in Figure 6 for the wavy LE with a 10% chord amplitude and 30% chord wavelength (10A30W), except that the boundary layer reaches the limiting profile at one LE wavelength.

Stability Validation
Prior to the wavy LE calculations, the stability code was validated by comparing our results with the benchmark solution of Grosch and Orszag [24] for Blasius flow eigenvalues at a single Reynolds number and wavenumber. Table 1 shows that a four-digit accuracy is achieved at a low grid resolution, and a six-digit accuracy is achieved at a higher grid resolution. Eigenvalues were Fluids 2020, 5, 212 9 of 23 also calculated at different Reynolds numbers and compared to the results from Mack [29] in Table 2, with good agreement.

Stability Validation
Prior to the wavy LE calculations, the stability code was validated by comparing our results with the benchmark solution of Grosch and Orszag [24] for Blasius flow eigenvalues at a single Reynolds number and wavenumber. Table 1 shows that a four-digit accuracy is achieved at a low grid resolution, and a six-digit accuracy is achieved at a higher grid resolution. Eigenvalues were also calculated at different Reynolds numbers and compared to the results from Mack [29] in Table 2, with good agreement. Next, the stability equations were solved for an infinite span swept flat plate with a straight LE using an FSC base flow and compared to the experimental measurements of Nitschke-Kowsky and Bippes [33]. The most unstable stationary mode parallel to the leading edge was calculated and compared to measurements at two different wind tunnel turbulence intensities, along with published numerical solutions from Nitschke-Kowsky and Bippes [33] and Collier et al. [46], who both used FSC similarity solutions for the base flow. Good agreement among the codes and experimental data is shown in Figure 7.  Next, the stability equations were solved for an infinite span swept flat plate with a straight LE using an FSC base flow and compared to the experimental measurements of Nitschke-Kowsky and Bippes [33]. The most unstable stationary mode parallel to the leading edge was calculated and compared to measurements at two different wind tunnel turbulence intensities, along with published numerical solutions from Nitschke-Kowsky and Bippes [33] and Collier et al. [46], who both used FSC similarity solutions for the base flow. Good agreement among the codes and experimental data is shown in Figure 7.

Zero-Sweep Flat Plate-Base Flow and Stability
The results for the base flow and stability of the boundary layer over a flat plate with a zero pressure gradient and zero sweep relative to the freestream were calculated for a straight and wavy LE. The base flow was in the chord direction with zero crossflow, with T-S waves present in the boundary layer. The analysis began by calculating base flows for different LE geometries using the three-dimensional BLE. Insights from the base flow analysis were used to simplify the stability analysis.

Zero-Sweep Flat Plate-Base Flow and Stability
The results for the base flow and stability of the boundary layer over a flat plate with a zero pressure gradient and zero sweep relative to the freestream were calculated for a straight and wavy LE. The base flow was in the chord direction with zero crossflow, with T-S waves present in the boundary layer. The analysis began by calculating base flows for different LE geometries using the three-dimensional BLE. Insights from the base flow analysis were used to simplify the stability analysis.

Base Flow
For the zero-sweep flat plate, the freestream conditions used the same resultant velocity magnitude as the DLR wind tunnel experiment for a swept flat plate, with the modification of the zero pressure gradient and zero sweep, reducing the base flow to the Blasius boundary layer. For the wavy LE geometry, four parameters were varied: the LE wavelength, LE amplitude, chord location, and span location. Four different base flow studies were performed by (1) changing the span location for a fixed chord location, LE amplitude, and LE wavelength; (2) changing the span and chord location for a fixed LE amplitude and LE wavelength; (3) changing the span location, chord location, and LE wavelength for a fixed LE amplitude; and (4) changing the chord location and LE amplitude for a fixed span location and LE wavelength. Figure 8 shows the boundary layer profiles at different span locations between the LE peaks for the straight and wavy LEs at a fixed 20% chord location, fixed wavy LE amplitude, and wavelengths of 10% and 30% chord, respectively. For the wavy LE, the boundary layer at the LE peak, labeled as 10A30W 1.0L, is identical to the Blasius profile. Deviations from the Blasius profile are seen at span locations between the LE peaks. The boundary layer thins across the span to a minimum at the trough and is symmetric about the trough.

Base Flow
For the zero-sweep flat plate, the freestream conditions used the same resultant velocity magnitude as the DLR wind tunnel experiment for a swept flat plate, with the modification of the zero pressure gradient and zero sweep, reducing the base flow to the Blasius boundary layer. For the wavy LE geometry, four parameters were varied: the LE wavelength, LE amplitude, chord location, and span location. Four different base flow studies were performed by (1) changing the span location for a fixed chord location, LE amplitude, and LE wavelength; (2) changing the span and chord location for a fixed LE amplitude and LE wavelength; (3) changing the span location, chord location, and LE wavelength for a fixed LE amplitude; and (4) changing the chord location and LE amplitude for a fixed span location and LE wavelength. Figure 8 shows the boundary layer profiles at different span locations between the LE peaks for the straight and wavy LEs at a fixed 20% chord location, fixed wavy LE amplitude, and wavelengths of 10% and 30% chord, respectively. For the wavy LE, the boundary layer at the LE peak, labeled as 10A30W 1.0L, is identical to the Blasius profile. Deviations from the Blasius profile are seen at span locations between the LE peaks. The boundary layer thins across the span to a minimum at the trough and is symmetric about the trough. The boundary layer thicknesses down the chord and at different span locations for the 10A30W wavy LE are shown in Figure 9. At the wavy LE peak, the boundary layer thickness is equivalent everywhere to the Blasius boundary layer thickness. At span locations between the LE peaks, the The boundary layer thicknesses down the chord and at different span locations for the 10A30W wavy LE are shown in Figure 9. At the wavy LE peak, the boundary layer thickness is equivalent everywhere to the Blasius boundary layer thickness. At span locations between the LE peaks, the boundary layer is thinner and asymptotically trends toward the Blasius boundary layer thickness as the chord location is increased. The boundary layer thicknesses down the chord and at different span locations for the 10A30W wavy LE are shown in Figure 9. At the wavy LE peak, the boundary layer thickness is equivalent everywhere to the Blasius boundary layer thickness. At span locations between the LE peaks, the boundary layer is thinner and asymptotically trends toward the Blasius boundary layer thickness as the chord location is increased. Next, the boundary layer profiles for the wavy LE were calculated at different chord and span locations for wavelengths of 30%, 15%, and 7.5% at a fixed 10% amplitude (10A30 W, 10A15 W, and 10A075 W). It was found that the boundary layer profile plotted in both dimensional and nondimensional units at a given span location was independent of wavelength. Figure 10 shows the boundary layer thickness across the chord for three different LE wavelengths at 1.2 L and 1.5 L span locations as an example of the independence of the boundary layer profile with the LE wavelength. Next, the boundary layer profiles for the wavy LE were calculated at different chord and span locations for wavelengths of 30%, 15%, and 7.5% at a fixed 10% amplitude (10A30 W, 10A15 W, and 10A075 W). It was found that the boundary layer profile plotted in both dimensional and nondimensional units at a given span location was independent of wavelength. Figure 10 shows the boundary layer thickness across the chord for three different LE wavelengths at 1.2 L and 1.5 L span locations as an example of the independence of the boundary layer profile with the LE wavelength. Finally, the effect of the LE amplitude was investigated by calculating boundary layer profiles at different chord locations for amplitudes of 5%, 10%, and 20% chord, with a fixed 30% wavelength and fixed span location at the wavy LE trough. As shown in Figure 11, increasing the amplitude at the trough has a similar effect on the boundary layer thickness as the results in Figure 9 show; Finally, the effect of the LE amplitude was investigated by calculating boundary layer profiles at different chord locations for amplitudes of 5%, 10%, and 20% chord, with a fixed 30% wavelength and fixed span location at the wavy LE trough. As shown in Figure 11, increasing the amplitude at the trough has a similar effect on the boundary layer thickness as the results in Figure 9 show; increasing the amplitude decreases the boundary layer thickness towards the LE, and the thicknesses tend towards the Blasius solution as the chord location increases. Finally, the effect of the LE amplitude was investigated by calculating boundary layer profiles at different chord locations for amplitudes of 5%, 10%, and 20% chord, with a fixed 30% wavelength and fixed span location at the wavy LE trough. As shown in Figure 11, increasing the amplitude at the trough has a similar effect on the boundary layer thickness as the results in Figure 9 show; increasing the amplitude decreases the boundary layer thickness towards the LE, and the thicknesses tend towards the Blasius solution as the chord location increases. In summary, the boundary layer at the LE peak is the same as the Blasius profile and thins across the span to a minimum at the trough. The boundary layer profiles are symmetric about the trough, and the boundary layer profiles across the wavy LE span are independent of the wavy LE wavelength, while the boundary layer thins as the LE amplitude increases. Finally, the deviation of the base flow from the Blasius profile is diminished with an increasing chord location; therefore, one can conclude that the LE geometry only modifies the base flow close to the LE. In summary, the boundary layer at the LE peak is the same as the Blasius profile and thins across the span to a minimum at the trough. The boundary layer profiles are symmetric about the trough, and the boundary layer profiles across the wavy LE span are independent of the wavy LE wavelength, while the boundary layer thins as the LE amplitude increases. Finally, the deviation of the base flow from the Blasius profile is diminished with an increasing chord location; therefore, one can conclude that the LE geometry only modifies the base flow close to the LE.

Stability
A parametric study of the boundary layer stability was performed for the straight and wavy LE flat plate. From the base flow results in the previous section, it was found that the LE wavelength has no effect on the base flow profile, and it was therefore not a parameter in the stability analysis. The stability results are presented in the form of neutral curves and critical Reynolds numbers, and require base flows at different chord locations. Thus, the two parameters in the stability analysis were the span location and LE amplitude.
The effect of the span location for a fixed 10% LE amplitude is shown in Figure 12. The neutral curves for the wavy LE peak and Blasius flow of the straight LE are equivalent, with a critical Reynolds number of approximately 302, which is in good agreement with Jordinson [22]. Both the critical Reynolds number and the range of unstable wavenumbers increase from the peak to the trough. The neutral curves are symmetric about the LE trough as expected from the symmetry of the base flows.
The effects of varying the LE amplitude at a fixed span location were calculated. Figure 13 shows that the effect of doubling the amplitude from 10% to 20% had a similar trend to varying the span location for a fixed amplitude, as expected from the base flow study. Both the critical Reynolds number and range of unstable wavenumbers were greater for the 20A case than the 10A case. The effects of the span location and LE amplitude on the critical Reynolds number are more clearly shown in Figure 14. The critical Reynolds number profile follows an oscillating pattern with the span location.
The maximum critical Reynolds number appears at the trough and 0.5 and 1.5 wavelengths, and is increased for a larger LE amplitude. base flows.
The effects of varying the LE amplitude at a fixed span location were calculated. Figure 13 shows that the effect of doubling the amplitude from 10% to 20% had a similar trend to varying the span location for a fixed amplitude, as expected from the base flow study. Both the critical Reynolds number and range of unstable wavenumbers were greater for the 20A case than the 10A case. The effects of the span location and LE amplitude on the critical Reynolds number are more clearly shown in Figure 14. The critical Reynolds number profile follows an oscillating pattern with the span location. The maximum critical Reynolds number appears at the trough and 0.5 and 1.5 wavelengths, and is increased for a larger LE amplitude.   In summary, the stability analysis for the zero-sweep plate with a zero pressure gradient and wavy LE produced similar conclusions to those of the base flow study. The neutral curve at the LE wave peak was equal to that of the Blasius profile. Both the critical Reynolds number and range of unstable wavenumbers were increased across the span to a maximum at the trough, and the neutral curves were symmetric about the trough. The wavy LE neutral curves all tended towards the Blasius neutral curve as the Reynolds number was increased. The most general conclusion of the stability analysis is that the leading-edge wave geometry leads to an increase in the critical Reynolds number above that of the straight LE, with the trough of the LE wave having the highest critical Reynolds number. In summary, the stability analysis for the zero-sweep plate with a zero pressure gradient and wavy LE produced similar conclusions to those of the base flow study. The neutral curve at the LE wave peak was equal to that of the Blasius profile. Both the critical Reynolds number and range of unstable wavenumbers were increased across the span to a maximum at the trough, and the neutral curves were symmetric about the trough. The wavy LE neutral curves all tended towards the Blasius neutral curve as the Reynolds number was increased. The most general conclusion of the stability analysis is that the leading-edge wave geometry leads to an increase in the critical Reynolds number above that of the straight LE, with the trough of the LE wave having the highest critical Reynolds number. unstable wavenumbers were increased across the span to a maximum at the trough, and the neutral curves were symmetric about the trough. The wavy LE neutral curves all tended towards the Blasius neutral curve as the Reynolds number was increased. The most general conclusion of the stability analysis is that the leading-edge wave geometry leads to an increase in the critical Reynolds number above that of the straight LE, with the trough of the LE wave having the highest critical Reynolds number.

Swept Flat Plate-Base Flow and Stability
A linear stability analysis for stationary crossflow disturbances was performed on an infinite span swept flat plate with straight and wavy LEs. The analysis was built around a DLR experiment [37] on crossflow stability using an effective sweep angle of 42.5°. The experiment considered here used a displacement body placed above a swept flat plate with a straight LE in a wind tunnel to create a strong favorable pressure gradient in the chord direction, stabilizing the T-S disturbances. The

Swept Flat Plate-Base Flow and Stability
A linear stability analysis for stationary crossflow disturbances was performed on an infinite span swept flat plate with straight and wavy LEs. The analysis was built around a DLR experiment [37] on crossflow stability using an effective sweep angle of 42.5 • . The experiment considered here used a displacement body placed above a swept flat plate with a straight LE in a wind tunnel to create a strong favorable pressure gradient in the chord direction, stabilizing the T-S disturbances. The combination of sweep and a favorable pressure gradient created a crossflow within the boundary layer. In the DLR experiment, end plates were placed on the swept flat plate to approximate infinite span flow conditions along the center of the plate. A parametric study of the flat plate with a wavy LE is presented. A systematic analysis of the base flow across the wavy LE was first performed to aid in the parametric study of the boundary layer stability. The results are first presented as a local analysis at a single span location, the wavy LE peak, and then a global analysis at span locations between successive peaks.

Local Analysis
The fully three-dimensional BLE were solved for a swept flat plate with straight and wavy LEs for one set of freestream conditions. Calculations were performed at a single span location, i.e., the peak, for the local analysis.
The freestream velocity profiles U(x) and V(x) were chosen to match the DLR experiment, with initial U(x = 0) = 4.9 m/s. This allowed for non-zero velocities to be present everywhere in the flow, except for the plate surface, and facilitate solution by the boundary layer code without the requirement of a self-similar FSC or Hiemenz flow initial condition.
For a wavy LE, the flow profile will approach that of a straight LE plate when the LE wavelength approaches both zero and infinity. It is therefore expected that an optimum LE wavelength between zero and infinity may exist that maximizes the wavy LE modification of the base flow. To estimate the optimized LE wavelength, base flow calculations were performed over a range of LE wavelengths at a fixed LE wave amplitude of 10% chord and a fixed 10% chord location, and compared to the infinite span solution for a straight LE. Figure 15 plots the crossflow peak velocity ratio, defined as the peak crossflow velocity normal to the freestream resultant vector for a wavy LE divided by the same velocity for the straight LE. It is seen in Figure 15 that the velocity modification factor approaches unity for a zero wavelength and for large wavelengths, and that there exists a wavelength near 25% chord that minimizes the crossflow velocity. at a fixed LE wave amplitude of 10% chord and a fixed 10% chord location, and compared to the infinite span solution for a straight LE. Figure 15 plots the crossflow peak velocity ratio, defined as the peak crossflow velocity normal to the freestream resultant vector for a wavy LE divided by the same velocity for the straight LE. It is seen in Figure 15 that the velocity modification factor approaches unity for a zero wavelength and for large wavelengths, and that there exists a wavelength near 25% chord that minimizes the crossflow velocity. The stability of the wavy LE was next investigated by calculating neutral curves for different LE geometries. First, the effect of different LE amplitudes at a fixed wavelength was investigated. Based on the results from Figure 15, a wavelength of 25% chord was chosen, and amplitudes of 5%, 10%, 20%, and 40% of the plate chord were used to calculate the neutrally stable eigenvalues. Figure 16 shows the neutral curves for each case, all of which retained the same general shape. The critical Reynolds number nearly reached a limiting value at a 20% amplitude. It is noted that each neutral  The stability of the wavy LE was next investigated by calculating neutral curves for different LE geometries. First, the effect of different LE amplitudes at a fixed wavelength was investigated. Based on the results from Figure 15, a wavelength of 25% chord was chosen, and amplitudes of 5%, 10%, 20%, and 40% of the plate chord were used to calculate the neutrally stable eigenvalues. Figure 16 shows the neutral curves for each case, all of which retained the same general shape. The critical Reynolds number nearly reached a limiting value at a 20% amplitude. It is noted that each neutral curve approached the straight LE neutral curve at a sufficiently high Reynolds number, and that all the neutral curves exhibited instability over a large range of wavenumbers.  Next, the effect of different LE wavelengths from 15% to 40% chord at a fixed 10% amplitude is shown in Figure 17. The 30% chord wavelength had the highest critical Reynolds number. The 25% chord wavelength, not shown in the figure, had a similar shape and nearly the same critical Reynolds number as the 30% wavelength. The 40% chord wavelength approached the limiting case of the straight LE. Figure 18 shows the critical Reynolds number for a range of wavelengths and amplitudes. The asymptotic behavior of different amplitudes at a fixed 25% chord wavelength is evident, along with the presence of both a local and global maximum critical Reynolds number at a fixed 10% amplitude. For the local analysis at the LE peak, the effect of the LE wavelength is an optimization problem, while the effect of increasing the LE amplitude is asymptotic.  Next, the effect of different LE wavelengths from 15% to 40% chord at a fixed 10% amplitude is shown in Figure 17. The 30% chord wavelength had the highest critical Reynolds number. The 25% chord wavelength, not shown in the figure, had a similar shape and nearly the same critical Reynolds number as the 30% wavelength. The 40% chord wavelength approached the limiting case of the straight LE. Figure 18 shows the critical Reynolds number for a range of wavelengths and amplitudes. The asymptotic behavior of different amplitudes at a fixed 25% chord wavelength is evident, along with the presence of both a local and global maximum critical Reynolds number at a fixed 10% amplitude. For the local analysis at the LE peak, the effect of the LE wavelength is an optimization problem, while the effect of increasing the LE amplitude is asymptotic.
Next, the effect of different LE wavelengths from 15% to 40% chord at a fixed 10% amplitude is shown in Figure 17. The 30% chord wavelength had the highest critical Reynolds number. The 25% chord wavelength, not shown in the figure, had a similar shape and nearly the same critical Reynolds number as the 30% wavelength. The 40% chord wavelength approached the limiting case of the straight LE. Figure 18 shows the critical Reynolds number for a range of wavelengths and amplitudes. The asymptotic behavior of different amplitudes at a fixed 25% chord wavelength is evident, along with the presence of both a local and global maximum critical Reynolds number at a fixed 10% amplitude. For the local analysis at the LE peak, the effect of the LE wavelength is an optimization problem, while the effect of increasing the LE amplitude is asymptotic.  Plots of the temporal amplification factor for a Reynolds number of 2300 at a fixed 10% amplitude and for a range of wavelengths are shown in Figure 19. The most amplified wavenumber is nearly constant for all the LE geometries, showing that the LE waviness has little effect on the most amplified wavenumber. Plots of the temporal amplification factor for a Reynolds number of 2300 at a fixed 10% amplitude and for a range of wavelengths are shown in Figure 19. The most amplified wavenumber is nearly constant for all the LE geometries, showing that the LE waviness has little effect on the most amplified wavenumber.

Global Analysis
A global stability analysis was performed for different LE geometries by calculating the critical Reynolds numbers and maximum amplification factors at span locations between successive LE peaks. The effect of the LE wavelength is presented first. Figure 20 shows the critical Reynolds numbers for wavy LE geometries at a fixed 10% amplitude and 30%, 15%, and 7.5% wavelengths (10A30W, 10A15W, and 10A075W, respectively). The critical Reynolds numbers are plotted against a nondimensionalized LE wavelength, with integer values of the nondimensionalized wavelength corresponding to the LE wave peak, and values of 0.5 and 1.5 corresponding to the wave trough. The 10A30W geometry exhibits strong oscillations across the span, with a peak critical Reynolds number near 1.0 (the dotted vertical line shows where the local analysis was carried out, Figure 17) and 2.0 and a low critical Reynolds number at or slightly below the straight leading-edge case near 0.25 and 1.25. The effect of sweep is characterized by the steep decline in the critical Reynolds number past the peak near the wavelength of 1.0 and a slow gradual rise past the wavelength of 1.25. By contrast, for the 10A15W and 10A075W geometries, the critical Reynolds number had less pronounced oscillations across the span, with both cases having similar critical Reynolds number when averaged over the LE wavelength. Additionally, it is worth noting that the maximum critical Reynolds number for both cases-10A15W and 10A075W-occurred near the trough: wavelengths 0.5 and 1.5. On average, the critical Reynolds numbers for all three cases were higher than the straight leading-edge case. Plots of the temporal amplification factor for a Reynolds number of 2300 at a fixed 10% amplitude and for a range of wavelengths are shown in Figure 19. The most amplified wavenumber is nearly constant for all the LE geometries, showing that the LE waviness has little effect on the most amplified wavenumber.

Global Analysis
A global stability analysis was performed for different LE geometries by calculating the critical Reynolds numbers and maximum amplification factors at span locations between successive LE peaks. The effect of the LE wavelength is presented first. Figure 20 shows the critical Reynolds numbers for wavy LE geometries at a fixed 10% amplitude and 30%, 15%, and 7.5% wavelengths (10A30W, 10A15W, and 10A075W, respectively). The critical Reynolds numbers are plotted against a nondimensionalized LE wavelength, with integer values of the nondimensionalized wavelength corresponding to the LE wave peak, and values of 0.5 and 1.5 corresponding to the wave trough. The 10A30W geometry exhibits strong oscillations across the span, with a peak critical Reynolds number near 1.0 (the dotted vertical line shows where the local analysis was carried out, Figure 17  The maximum amplification factors at a fixed 0.124 chord fraction are shown in Figure 21. This chord fraction was chosen since it was far enough down the chord that unstable disturbances were present, but not so far down the chord that the effects of the wavy LE on the base flow were diminished. The maximum amplification factor for the 10A30W geometry had an oscillating profile, with the amplification being smaller everywhere than the straight LE case. The 10A15W and 10A075W maximum amplification factor profiles were more constant and less amplified than the straight LE geometry. The 10A075W maximum amplification factor was larger in magnitude everywhere than the 10A15W case, trending towards the straight LE maximum amplification value.  diminished. The maximum amplification factor for the 10A30W geometry had an oscillating profile, with the amplification being smaller everywhere than the straight LE case. The 10A15W and 10A075W maximum amplification factor profiles were more constant and less amplified than the straight LE geometry. The 10A075W maximum amplification factor was larger in magnitude everywhere than the 10A15W case, trending towards the straight LE maximum amplification value. The maximum amplification factors at a fixed 0.124 chord fraction are shown in Figure 21. This chord fraction was chosen since it was far enough down the chord that unstable disturbances were present, but not so far down the chord that the effects of the wavy LE on the base flow were diminished. The maximum amplification factor for the 10A30W geometry had an oscillating profile, with the amplification being smaller everywhere than the straight LE case. The 10A15W and 10A075W maximum amplification factor profiles were more constant and less amplified than the straight LE geometry. The 10A075W maximum amplification factor was larger in magnitude everywhere than the 10A15W case, trending towards the straight LE maximum amplification value. The full amplification curves at a 0.124 chord fraction are shown in Figure 22 for different span locations of the 10A30W geometry and the straight LE. It is seen that the wavenumber of the most amplified mode was nearly constant at all span locations. Similar to the local analysis of Figure 19, the LE waviness had little effect on the most amplified wavenumber across the span. The full amplification curves at a 0.124 chord fraction are shown in Figure 22 for different span locations of the 10A30W geometry and the straight LE. It is seen that the wavenumber of the most amplified mode was nearly constant at all span locations. Similar to the local analysis of Figure 19, the LE waviness had little effect on the most amplified wavenumber across the span. The global effect of the LE amplitude was calculated for 30% and 7.5% LE wavelengths. Figure  23 shows the critical Reynolds numbers for the 05A30W and 10A30W geometries. For the 5% and 10% amplitudes, the profiles oscillated across the span. Increasing the amplitude had a beneficial effect on the stability for all cases, and the maximum global critical Reynolds number was increased as the amplitude increased. Portions of the 5% and 10% amplitude cases were less stable than the straight LE geometry. Figure 24 shows the global effect of LE amplitude for the 05A075W and 10A075W geometries. Each of the 7.5% wavelength profiles were more constant along the span than the 30% wavelength profiles, and the boundary layer for the 7.5% wavelength was more stable than the straight LE everywhere for 5% and 10% amplitudes, in contrast to the 30% wavelength geometry. The global effect of the LE amplitude was calculated for 30% and 7.5% LE wavelengths. Figure 23 shows the critical Reynolds numbers for the 05A30W and 10A30W geometries. For the 5% and 10% amplitudes, the profiles oscillated across the span. Increasing the amplitude had a beneficial effect on the stability for all cases, and the maximum global critical Reynolds number was increased as the amplitude increased. Portions of the 5% and 10% amplitude cases were less stable than the straight LE geometry. Figure 24 shows the global effect of LE amplitude for the 05A075W and 10A075W geometries. Each of the 7.5% wavelength profiles were more constant along the span than the 30% wavelength profiles, and the boundary layer for the 7.5% wavelength was more stable than the straight LE everywhere for 5% and 10% amplitudes, in contrast to the 30% wavelength geometry. As the amplitude was decreased for the 7.5% wavelength case, the nearly constant critical Reynolds number profiles approached the straight LE value with little oscillation. The global effect of the LE amplitude was calculated for 30% and 7.5% LE wavelengths. Figure  23 shows the critical Reynolds numbers for the 05A30W and 10A30W geometries. For the 5% and 10% amplitudes, the profiles oscillated across the span. Increasing the amplitude had a beneficial effect on the stability for all cases, and the maximum global critical Reynolds number was increased as the amplitude increased. Portions of the 5% and 10% amplitude cases were less stable than the straight LE geometry. Figure 24 shows the global effect of LE amplitude for the 05A075W and 10A075W geometries. Each of the 7.5% wavelength profiles were more constant along the span than the 30% wavelength profiles, and the boundary layer for the 7.5% wavelength was more stable than the straight LE everywhere for 5% and 10% amplitudes, in contrast to the 30% wavelength geometry. As the amplitude was decreased for the 7.5% wavelength case, the nearly constant critical Reynolds number profiles approached the straight LE value with little oscillation.   In summary, for the swept flat plate, the wavy LE had stability benefits verses a straight LE for a range of LE wavelengths and amplitudes. In general, the wavy LE increased the critical Reynolds number across the span and reduced the temporal amplification factor. The most amplified wavenumber was relatively constant across the span and was nearly equal to that of the straight LE. Smaller wavelengths seem to have a stability advantage, in that the critical Reynolds number was increased above the straight LE value for all the span locations and all the LE amplitudes tested.

Discussion
This work performed a temporal linear stability analysis for an incompressible fluid over a flat In summary, for the swept flat plate, the wavy LE had stability benefits verses a straight LE for a range of LE wavelengths and amplitudes. In general, the wavy LE increased the critical Reynolds number across the span and reduced the temporal amplification factor. The most amplified wavenumber was relatively constant across the span and was nearly equal to that of the straight LE. Smaller wavelengths seem to have a stability advantage, in that the critical Reynolds number was increased above the straight LE value for all the span locations and all the LE amplitudes tested.

Discussion
This work performed a temporal linear stability analysis for an incompressible fluid over a flat plate with a straight and wavy LE, for both zero and 42.5 • sweep. The plate geometry and freestream conditions were chosen to match the experiment performed at the DLR for the swept plate analysis. The freestream conditions for the zero-sweep plate were modified from the DLR conditions by using a zero pressure gradient in the chord direction. The three-dimensional BLE were used to calculate the base flows for straight and wavy LE geometries. The disturbances for the zero-sweep plate consisted of traveling waves in the chord direction, while stationary disturbances in the crossflow direction were calculated for the swept plate. For both swept and un-swept plates, the modification of the base flow by the LE wave was most significant near the leading edge. The wavy LE base flow tended toward the straight LE base flow as chord location increased. For both the swept and zero-sweep geometries, the wavy LE resulted in higher critical Reynolds numbers and decreased the temporal amplification factor of the most amplified wavenumber at chord locations near the LE.
For the zero-sweep flat plate with a zero pressure gradient, the LE wave modified the base flow profile with respect to the Blasius boundary layer by thinning the boundary layer from the LE peak to the trough, with a minimum boundary layer thickness at the trough and symmetric boundary layer profiles with respect to the trough. The Blasius boundary layer profile occurred at the LE peak. The LE wavelength had no effect on the boundary layer profiles, while increasing the LE amplitude increased the critical Reynolds number at a fixed span location. Compared to the straight LE, the wavy LE geometry increased the critical Reynolds number and decreased the temporal amplification at all the span locations between the LE peaks, with the trough being the most stable span location.
For the swept flat plate, a local linear stability analysis was first performed at the wavy LE peak. Unlike that for the zero-sweep flat plate, the base flow for the swept flat plate was dependent on both the LE wavelength and amplitude. The effect of the LE wavelength was an optimization problem, while the effect of increasing the LE amplitude was asymptotic, with the 10A30W geometry being chosen as the most stable LE geometry at the peak since it had the largest critical Reynolds number. The most amplified wavenumber was also nearly constant and approximately equal to that of the straight LE.
A global stability analysis for the swept plate was performed at span locations between the LE peaks. In general, the wavy LE increases the critical Reynolds number across the span and reduces the temporal amplification factor. The most amplified wavenumber for the swept plate is nearly constant and equal to that of the straight LE plate across the span and chord. For larger LE wavelengths, the critical Reynolds number profiles oscillated across the span, while the smaller LE wavelengths had more constant span profiles. For a fixed LE amplitude, the largest LE wavelength of 30% chord produced the largest maximum critical Reynolds number, but due to the oscillating profile, the stability at portions of the span was the same or slightly worse than that for the straight LE. The smallest wavelengths of 15% and 7.5% chord produced more constant critical Reynolds number profiles across the span, being more stable than the straight LE everywhere. Varying the amplitude at a fixed 7.5% wavelength produced relatively constant span profiles at all amplitudes, with the average critical Reynolds number approaching that of the straight LE as the amplitude was decreased. For the 30% wavelength case, all the LE amplitudes produced an oscillating profile across the span. Amplitudes of 5% and 10% chord and a 30% wavelength had portions of the span with a smaller critical Reynolds number than the straight LE.