12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

EFFECT OF STREAMWISE SURFACE UNDULATIONS ON THENONLINEAR STABILITY OF CROSSFLOW INSTABILITIESS. Westerbeek
Faculty of Aerospace Engineering
TU Delft
Leeghwaterstraat 42, 2628 CA Delft
S.H.J.Westerbeek@tudelft.nl

M. Kotsonis
Faculty of Aerospace Engineering
TU Delft
Leeghwaterstraat 42, 2628 CA Delft
M.Kotsonis@tudelft.nl

ABSTRACT

importance of surface wavelength and phase. The existence of
surface undulations is significant to the stability of the flow andcan result in both positive (stabilization) and negative (destabilization) effects on the dominant mode. The parameter space islarge and requires an efficient (nonlinear) stability calculation
framework to be analyzed thoroughly.
The Nonlinear Parabolized Stability Equations (NPSE)
are commonly used to calculate the stability of growing boundary layers. Nevertheless, their application is limited due tothe assumption of slowly varying perturbation shape. This assumption appears to be valid for smooth surfaces with limited
localized surface gradients as shown by Park & Park (2013)
and Park & Oh (2020). The PSE have therein been shown tobe effective in predicting the stability around localized curvature for 2-dimensional boundary layers. Franco Sumariva et al.
(2020) show that sharp geometric features break this assumption and require a fully elliptic stability method instead. Additionally, the step size limitation, inherent in the PSE methodology, is mentioned as a crucial limitation for stability oversmooth humps. A relaxation of the step size restriction thatwould allow for finer streamwise discretization is not considered explicitly in that work. Thomas et al. (2016) analyzed
the linear stability of a swept wing featuring sunsoidal surface
deformations with LPSE and Linear Harmonic Navier-Stokes(LHNS), thus establishing some validity of the PSE assumption in undulated surface boundary layer cases, even for caseswith small separation regions in the basic state.
Based on the afformentioned studies, the PSE are used
for this problem due to the unrivaled speed by which nonlinear simulations can be performed. This allows for a parametric study to be performed that considers also the shape of theundulations. This study aims to analyze the effect of undulation shape on the stability of three-dimensional boundary layers and show the effectiveness of the PSE framework for thistype of problem.

The nonlinear stability of three-dimensional boundary
layers over various undulated surfaces was calculated using the generalized Nonlinear Parabolized Stability Equations(NPSE). The results are compared with a flat plate configuration to assess the effect of the undulation shape on the stabilityof the boundary layer. It was found that the effect of surface
undulations is significant and should not be ignored when performing stability analysis. All undulation shapes considered inthis work showed a destabilization of the primary mode and theassociated harmonics. The stability of the boundary layer was
directly affected by the amplitude of the undulations, while
their respective shape did not meaningfully affect the evolution of the crossflow instabilities within the parameter range
considered in this work.

Introduction
The laminar-to-turbulent transition of boundary layers directly affects the skin friction drag experienced by an aircraft
wing resulting in increased fuel costs and associated emissions. At the same time, the manufacturing of wings is subject
to limitations resulting in necessary junctions (gaps and steps)
and bolts on the wing surface. Furthermore, plate rolling techniques and operational loads affect the surface smoothness of
the wing plate itself, by inherently creating long-wavelength
streamwise waviness. Despite the presence of such geometrical deviations from an ideal wing and their long known importance to the transition scenario (see e.g. Holmes et al. (1985);
Nayfeh et al. (1988) ), stability calculations are often still performed assuming a perfectly smooth surface. Accounting for
finer physical details in the stability calculations of boundary
layers is crucial for the design of realistic laminar wings that
aim to delay the onset of turbulence.
In one of the earliest works concerning the stabilityof three-dimensional boundary layers over undulated plates,
Masad (1996) found a destabilizing effect that was positively
correlated with the amplitude of the undulations. Thomas et al.
(2016) instead found a stabilization of the primary mode forsome cases with sinusoidal surface undulations and cites theNondimensionalization
In this work, all quantities are presented nondimensionally. Distances are normalized by the boundary layer thick-1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

H

U

W

y
VLE

Figure 1.

z

η
x

z

ξ

L

Schematic of a sinusoidally undulated swept plate and the computational domain.

ness, taken as δ99 at the inflow, i.e. the reference length denoted lre f = 7.71 · 10−4 m. The inflow velocity U0 = 15.1 m/s
is used to normalize velocities. A global Reynolds number
is defined using these constants and the kinematic viscosityU ∗lν = 1.51 · 10−5 m/s2 as: Re = 0 ν re f .

hyw (x) = r(x) · h · Norm tan−1syw (x) = r(x) · h ·Problem Statement
The stability of the boundary layers over swept plateswith respect to the freestream flow direction is considered. A
schematic of the problem is shown in figure 1. The virtual
leading edge is shown upstream of the computational stability
domain meaning that a pre-existing boundary layer and crossflow instability is present at the inflow. The computational domain L = 470 · lre f (0.3m). The undulation starts at a distanceof 27·lre f (2 cm) downstream of the inflow, ramping up to the
maximum amplitude at a distance of 162 lre f (12.5 cm) downstream of the inflow. This is explained in more detail later inthis section. The spanwise component of the external flow isconstant, W∞ = −1.24 ·U∞ , the streamwise velocity is a function of the leading edge orthogonal dimension, x, following
the external velocity distribution of the experiments of RiusVidales & Kotsonis (2021) as described in Casacuberta et al.
(2021).

isin(s)
k + 1 − cos(s)1 + b2sin(s)
1 + b2 · sin2 (s)

(1)

(2)

Here, Norm indicates a normalization of the shape and(x−S)2πs(x) = λ, where λw is the streamwise wavelength of thew
undulations, S is a constant offset before which the plate is flat,h is the undulation height,
skew parameter.
 and k is a constantThe ramp, r(x) = tanh 20 · lre f (x − S) , ensures a smooth introduction of the waviness. Note that the cases described by
k = 1000 and b = 0 are the same case as these shapes approximate a sine function.
The amplitude can be adjusted through the parameter halthough it is maintained at h = 0.3059 · lre f in this work. Still,
it will be shown that this small undulation amplitude has a noticeable effect on the stability of the base flow. The base flow
does not feature any separation for the cases presented here.
The geometric shape of the undulation can be adjusted to be
more square by varying the parameter b in equation 2. The
resulting surfaces are shown below in figure 2.

The transition scenario on swept-wings is generally governed by the evolution of stationary crossflow instabilities dueto the low freestream turbulence found in free-flight conditions
as noted by Bippes & Müller (1990) and Deyhle & Bippes(1996). In this work, therefore, no travelling crossflow modes
are considered and for all modes ω = 0. For the flat plate case,it was found that the mode specified by the spanwise wavelength of λz = 7.5mm is the most unstable mode and thus thismode will be introduced at the inflow for all cases. It is possible that the undulations affect the stability of the boundary
layer in such a way that other modes govern the transition scenario instead. A more robust transition analysis would therefore require a nonlinear N-envelope analysis that is consideredout of the scope of this work however.

Methodology
The stability of swept undulated plates is considered nonlinearly using the NPSE methodology. The NPSE are derived for generalized coordinates, parabolizing the equationsin the body-fitted grid along the streamwise direction. The
basic state is found using the steady incompressible laminar
Navier Stokes formulation in the commercial finite element
solver COMSOL. The details for both the base flow solver and
PSE methods are presented in this section.
In COMSOL, the inflow is prescribed by an analytical
Falkner-Scan-Cooke solution. At the wall, the no-slip condition is imposed and a pressure description of a favourablepressure gradient is prescribed at the top boundary as foundin Casacuberta et al. (2021). This top boundary is placed at aheight, H, which through preliminary investigations is foundto have no influence on the results (H = 51.9 ∗ lre f ). The outlet prescribes a constant static pressure condition equal to thestatic pressure found at the top boundary at that streamwise
location. Two finite elements are present in the spanwise direction, z, with periodic boundary conditions. Velocities and

The nonlinear stability of several variations of streamwise surface undulations is considered. After a flat plate region downstream of the inflow, a growing sinusoidal function is altered to either contain a skew or a flat top to assess the effect of the undulations’ shape. For the undulations,the shape is described by algebraic functions subject to either the skew parameter k(= [3, 5, 7, 1000]), or shape parameter b(= [0, 1, 1.5, 2]). These parameters appear in the surface
shape functions as:2

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

(a)
0.2

0

-0.2

(b)
0.2

0

-0.2

100

150

200

250

300

Figure 2. Surface undulations (not to scale) for the (a) skewed sinusoidals with k =[1000 7 5 3] and (b) squared sinusoidals withb =[2.5 2 1.5 0] in orange, blue, purple and green respectively.

pressures are calculated using second-order elements of which
there exist 78 in the wall-normal direction clustered to the walland 877 equispaced in the streamwise direction, using a direct
solver procedure.







L=





To analyze the stability of the basic state, a body-fitted coordinate system is introduced with the ξ -axis at η = 0 being locally parallel and normal to the wall surface. The η-coordinate
is constructed such that the grid is locally orthogonal at thewall as is required for parabolization along the ξ -axis. After aperturbation analysis is performed on the generalized NavierStokes equations, the nonlinear solution ansatz is introduced:

L + ∂U
∂x

∂U
∂y

∂V
∂x

L + ∂V
∂y

∂W
∂x

∂W
∂y

i( ξs α(ξ )dξ ∗ +βn z−ωm t)

q0 (ξ , η, z,t) = ∑ ∑ q̂m,n (ξ , η)e

+ c.c

(3)
Here q0 is the perturbation state vector [u0 , v0 , w0 , p0 ]T describedby a truncated sum of modes denoted by m and n for angularfrequency ω, and spanwise wavelength, β , respectively. The
streamwise wavenumber, α, is complex with the real part describing the wavenumber and the imaginary part the growthrate. The ∗ denotes an integration variable. Lastly, t is thetime, i is the imaginary unit and c.c denotes the complex conjugate.



 i 2

− Re (ξx + ξy2 )
0
0
0




i
2
2

0
− Re (ξx + ξy )
0
0



N=


i
2 + ξ 2 ) 0

0
0
−
(ξ
y


Re x



The resulting equations, in generalized coordinates, arepresented in system form as:

0

∂α
∂ q̂
+N
q̂ = r
∂ξ
∂ξ

0





 0 M 0 ξy 



M=


0 0 M 0ξx ξy 0 0

m n

Lq̂ + MM 0 0 ξx
0 ηy D1 + iαξy 





L
iβ



iαξx + ηx D1 iαξy + ηy D1 iβ



R ξe

0 iαξx + ηx D1



(4)

0

0

−u ∂∂ ξu ξx − u ∂∂ηu ηx − v ∂∂ ξu ξy − v ∂∂ηu ηy − iβ uw0

−u ∂ v ξ − u ∂ v η − v ∂ v ξ − v ∂ v η − iβ vw ∂ξ x

∂η x
∂η y
∂ξ y


r = −u ∂ w ξx − u ∂ w ηx − v ∂ w ξy − v ∂ w ηy − iβ w2  ,
 ∂ξ

∂η
∂η
∂ξ



Where L, M and N are:03

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

10 0

10 -1

10 -2

10 -3

10 -4

100

150

200

250

300

350

100

150

200

250

300

350

Figure 3. Amplitude evolution of the 1st (thick line), 2nd (dashed), 3rd (dash-dotted), 4th (thin line) and 5th (thin dash-dotted)
harmonic and the mean flow distortion (dotted) over swept plates with (a) skewed surface undulations k =[1000 7 5 3] and (b) squaresinusoidals with b =[0 1 1.5 2] in orange, blue, red and green respectively. The reference CFI evolution over a flat plate is shown inblack.
With D1 and D2 derivative matrices andvalid, a linear validation was also performed neglecting the∂p, as will be shown in the following section.
∂ξGuaranteeing the numerical stability of the NPSE is difficult for the discussed cases, though it can be improved inseveral ways if a fine enough base flow simulation is used. In
the present work, some additional measures are taken to thisgoal, by adding the mean flow distortion to the base flow after
a transformation of the derivatives to the physical coordinate
frame. This effectively reduces the weight of the right-handside, as the effect of the mean flow distortion is accounted forin the left-hand side rather than through modal interactions as
was suggested by Zhao et al. (2016). Additionally, the nonlinear forcing term is subjected to an under-relaxation factorof 0.2 in the nonlinear convergence iteration. Lastly, the convergence procedure is split in two steps as shown by Park &Park (2011), specifically: First, for a given forcing term, theα is converged for all modes. Then, the forcing term is updated and the process is repeated until, again, the streamwise
wavenumber is converged.
The NPSE are discretized in 400 stages in the streamwise
direction (∆ξ = 0.88) using a first-order backward implicit
Euler scheme. In the wall-normal direction spectral colocation is employed incorporating 80 Chebyshev polynomials onan equal number of collocation points that are clustered near
the wall. The initial condition for the primary mode is created by solving the local eigenvalue problem for the perturbation velocities, pressure and the complex-valued streamwise
wavenumber. The desired nondimensional initial amplitude,
A0 = 1.75 · 10−3 based on the maximum of the streamwise
perturbation velocity, is then imposed on the primary mode

L = −iω + iUαξx +Uηx D1 + iαV ξy +V ηy D1
1
1+ iβW + (α 2 (ξx2 + ξy2 ) + β 2 ) − (ηx2 + ηy2 )D2
Re
Re
1
1
− 2iα(ηy ξy + ηx ξx )D1 − iα(ξyy + ξxx )ReRe
1− (ηxx + ηyy )D1
(5a)
Re
and
M = Uξx +V ξy −
−

1
2iα 2
(ξx + ξy2 ) − (ξxx + ξyy )ReRe2(ηy ξy + ηx ξx )
D1
Re

(6a)

It must be noted that the base flow derivatives are stillpresent in these equations as derivatives with respect to thephysical xy-coordinate system as they are constants to thestability problem. Furthermore, the last column of M canbe neglected, effectively multiplying the streamwise pressure
derivative ∂∂ ξp by 0 for all modes as suggested by Li & Malik(1996), in order to relax the step size restriction of the NPSE.
Although this step can be performed reliably in most flat platestability simulations, more complex geometries require a careful consideration as the streamwise pressure gradient of theperturbations might not be negligible. To ensure this step is4

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

4

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

2
0
4

1.6

2
1.4

0
0

10

20

0

10

20

0

10

20

0

10

20
1.2

4

(i)

1

2
0.8

0
80
4

100

120

140

160

180

200

220
0.6

(j)

2

0.4

0
240
4

260

280

300

320

340

360

380

0.2
0

(k)

2
0
400

420

440

460

480

500

520

540

Figure 4. Fully reconstructed velocity fields based on the sum of base flow and instability modes. Streamwise velocity in the yz-plane
for the flat plate case (a-d) and the undulated plate (e-h) with k = 1000 at x = 182 (a,e), x = 305 (b,f), x = 367(c,g) and x =427 (d,h)
with the y = 0 representing the local wall location. Contours of streamwise velocity in the xy-plane over the entire streamwise domain(i,j,k).

and its complex conjugate. The nonlinear terms are incorporated explicitly in the equations through a forcing term on theright-hand side. For all stability simulations, five harmonics
and the mean flow distortion are considered while higher harmonics are suppressed even if they are predicted to cross thenonlinear introduction threshold of A = 1 · 10−08 .
The base flows capture all features and the base flow profiles and derivatives were ensured to be smooth throughout thedomain, however, the base flow was not subjected to a systematic grid convergence study and it is hypothesized that thenumerical stability and convergence of the NPSE can be significantly improved through the use of finer base flows.

ward the validity of the slowly varying flow assumption forboth the linear and nonlinear system.
Figure 5 shows a comparison between the linear amplitude evolution as calculated per PSE and HNS for the case withthe strongest wall gradients (b = 1.5), proving the validity ofneglecting the higher order derivatives in the PSE derivation.

Results
Base flow simulations were performed for all undulated
plates of section and the reference swept flat plate case. The
basic states were then used to perform nonlinear stability calculations using the NPSE. The base flows show no flow separation due to the mild amplitude of the waviness.
The amplitude evolution is shown separately for theskewed sinusoidal and squared sinusoidal cases in figure 3(a)and (b) respectively. In all reported cases, the growth of thecrossflow instability is enhanced by the surface undulations.
The local effect of the surface undulations is clearly visible forall cases, although the differences between the cases are minor indicating that the shape is not critically important to thestability response of the boundary layer within the investigatedparameter space.

Validation of the PSE Assumption in mildly undulating flows
The PSE assumption of slowly varying flow is validated
for the primary mode by comparing the linear PSE result withresults from a harmonic Navier-Stokes (HNS) simulation that
does not neglect higher-order derivatives. This validation is assumed to be independent of the harmonic considered as it applies to the shape function, rather than the total mode. Therefore, a linear match between the two methods would point to-5

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
for wall-undulations is too large to cover in this work. The
effects of wavelength, phase, shape and amplitude need
to be considered. This work could only tackle the effectof undulation shape limited by the PSE assumptions. All
cases caused an accumulated destabilization of the crossflow
instabilities which could lead to a significant advancement of
the transition front. Consequently, it is concluded that surface
undulations have to be accounted for in stability analysis if
present at amplitudes of a similar scale as presented here.

100

10-1

10-2

100

150

200

250

300

350

REFERENCES
Bippes, H & Müller, B 1990 Disturbance growth in an unstable three-dimensional boundary layer. In Numerical andPhysical Aspects of Aerodynamic Flows IV, pp. 345–358.
Springer.
Casacuberta, Jordi, Hickel, Stefan & Kotsonis, Marios 2021
Mechanisms of interaction between stationary crossflow instabilities and forward-facing steps. In AIAA Scitech 2021
Forum, p. 0854.
Deyhle, H & Bippes, H 1996 Disturbance growth in an unstable three-dimensional boundary layer and its dependence on
environmental conditions. Journal of Fluid Mechanics 316,
73–113.
Franco Sumariva, Juan Alberto, Hein, Stefan & Valero, Eusebio 2020 On the influence of two-dimensional humproughness on laminar–turbulent transition. Physics of Fluids 32 (3), 034102.
Holmes, Bruce J, Obara, Clifford J, Martin, Glenn L & Domack, Christopher S 1985 Manufacturing tolerances fornatural laminar flow airframe surfaces. SAE Transactions
pp. 522–531.
Li, Fei & Malik, Mujeeb R 1996 On the nature of pse approximation. Theoretical and Computational Fluid Dynamics 8 (4), 253–273.
Masad, Jamal A 1996 Effect of surface waviness on transitionin three-dimensional boundary-layer flow. Tech. Rep..
Nayfeh, Ali H, Ragab, Saad A & Al-Maaitah, Ayman A 1988
Effect of bulges on the stability of boundary layers. The
Physics of fluids 31 (4), 796–806.
Park, Donghun & Oh, Sejong 2020 Effect of shape of twodimensional smooth hump on boundary layer instability.
International Journal of Aeronautical and Space Sciences
21 (4), 906–923.
Park, Donghun & Park, Seung O 2013 Linear and non-linear
stability analysis of incompressible boundary layer over a
two-dimensional hump. Computers & Fluids 73, 80–96.
Park, Dong-Hun & Park, Seung-O 2011 Nonlinear stability
analysis of boundary layers by using nonlinear parabolizedstabiltiy equations. Journal of the Korean Society for Aeronautical & Space Sciences 39 (9), 805–815.
Rius-Vidales, Alberto F & Kotsonis, M 2021 Impact of aforward-facing step on the development of crossflow instability. Journal of Fluid Mechanics 924.
Thomas, Christian, Mughal, Shahid M, Gipon, Matthew, Ashworth, Richard & Martinez-Cava, Alejandro 2016 Stabilityof an infinite swept wing boundary layer with surface waviness. AIAA Journal 54 (10), 3024–3038.
Zhao, Lei, Zhang, Cun-bo, Liu, Jian-xin & Luo, Ji-sheng 2016
Improved algorithm for solving nonlinear parabolized stability equations. Chinese Physics B 25 (8), 084701.

Figure 5. Linear amplitude evolution of the primary crossflow instability mode (λz = 9.7 · łre f , ω = 0) as calculated per
PSE (line) and HLNS (dashed) for the case b = 1.5.

It should be noted that the validity of NPSE, judged hereby the ability of the NPSE to converge in the presented domain
and providing a linear match with HNS results, is limited anddepends strongly on the geometry (shape and amplitude) of
the wall. For the current case, a greater undulation amplitude(h > 0.3059 · lre f ), stronger skew (k < 3) or more square shape(b > 3) could affect the accuracy of the results.
In figure 4, streamwise velocity contours taken from thefully reconstructed flowfield are depicted. In figures 4(a-d) the
evolution over a flat plate is shown, which can be directly compared to the results of 4(e-h) below, showing the contours at thesame streamwise location for the sinusoidal plate. The comparison highlights the deformation of the developed crossflow
instability in the strongly nonlinear regime for the undulated
plate. In figures 4(i-j) xy-planes of the streamwise velocity
contours are shown demonstrating the relative scale of the surface undulations compared to the incoming perturbation. The
streamwise wavelength of the undulation is approximately six
times that of the perturbation and the height is a just under
half the vortex core height. The effect of the mild undulations
is visible in the contours of the crossflow vortex which is inaccordance with the discovered effect on the amplitude evolution.

Conclusion
The effect of surface undulations on the nonlinear evolution of crossflow instabilities was considered for several wall
surface undulation shapes.
The local effect of the wall shape variations was foundto be strong and can be directly explained by considering theeffect on the pressure gradient and curvature. Any undulationconsidered in this work strongly affects the local stability ofthe primary mode directly. Across all undulation shapes a rising wall expectedly promotes the growth of the incoming CFI
as a result of the pressure gradient becoming more favourable.
After the apex, the pressure gradient becomes more adverse
and the growth of crossflow instabilities is instead reduced.
The behaviour of higher harmonics is found to be more
affected by the nonlinear forcing, governed by the primarymode, than by the linear effect of the wall shape present in
the base flow and curvature terms.
As mentioned in the introduction, the parameter space

6