Spatial direct numerical simulation of transitional boundary layer over compliant surfaces

Computers & Fluids 34 (2005) 1062–1095 www.elsevier.com/locate/compfluid

Spatial direct numerical simulation of transitional boundary layer over compliant surfaces Z. Wang a, K.S. Yeo

b,*

, B.C. Khoo

b,c

a

c

Temasek Laboratories, National University of Singapore, 1 Engineering Drive 2, Singapore 117576, Singapore b Department of Mechanical Engineering, National University of Singapore, 1 Engineering Drive 2, Singapore 117576, Singapore Singapore MIT Alliance, National University of Singapore, 1 Engineering Drive 2, Singapore 117576, Singapore Received 13 August 2003; received in revised form 27 May 2004; accepted 30 August 2004 Available online 24 December 2004

Abstract An iterative implicit fractional step method is developed and employed for the simulation of transitional boundary layer over compliant surfaces. The three-dimensional perturbation Navier–Stokes equations are discretized by curvilinear finite volumes on a collocated grid system. A multigrid procedure is used for computations associated with the pressure-Poisson equation, while simulation is carried out by MPI-based parallel computation with domain decomposition. Results in the literature for oblique linear waves and the non-linear breakdown of a wave triad over a rigid wall are repeated to check the accuracy of the codes that had been developed. Oblique linear TS (Tollmien–Schlichting) waves over finite-length compliant membranes are generally found to coexist with CIFI (compliance induced flow instability) or FISI (flow induced surface instabilities) waves. The latter waves usually possess longer wavelengths and thus propagate at a larger oblique wave angles than the TS waves. Simulation reveals that compliant surfaces may slow down the development of secondary instabilities during the early stages of laminar–turbulent transition. However, during the later stages of fundamental wave breakdown, interactions with CIFI and edge-generated waves may increase the amplitudes of the original 2D and 3D TS waves, leading to an earlier breakdown on compliant surfaces. Linear interaction between the flow and compliant membranes has been assumed.
2004 Elsevier Ltd. All rights reserved.

*

Corresponding author. E-mail address: [email protected] (K.S. Yeo).

0045-7930/$ - see front matter
2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2004.08.005

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

1063

1. Introduction Many theoretical studies have found that compliant surfaces may stabilize two-dimensional Tollmien–Schlichting (TS) waves in a Blasius boundary layer, thus extending the linear stage of disturbance wave growth and delaying the onset of non-linear boundary-layer transition. In addition to the Tollmien–Schlichting instability (TSI) waves, new types of instability waves were discovered in these studies. All the waves arising from the interaction of the flow and the compliant surface can be classified into three categories according to the energy needed to activate or excite the disturbance [1,2]. The Class A waves, for which the TS wave is an example, are characterized by negative activation energy. The class B instability waves, for which the travelingwave flutter instability is an example, require positive activation energy for their excitation and growth. For the Class C instability (KH instability for instance) the activation energy is nearly zero. Later, Carpenter and Garrad [3] coined the term Ôflow-induced surface instabilityÕ (abbreviated FISI) to describe all waves that have their direct origin in the compliant quality of the wall to distinguish them from the TSI, which may exist even when the wall is rigid. The FISI was also termed as Ôcompliance-induced flow instabilityÕ (CIFI) by Yeo [4]. This latter terminology emphasizes the essentially passive quality of the compliant walls that he investigated. Stability analysis of Yeo [4] found that three-dimensional (3D) oblique wave modes might grow faster than their 2D counterparts for some compliant boundaries. This is particularly so for 3D TS modes on soft compliant surfaces. This conclusion was further confirmed by Joslin et al. [5] and Yeo [6]. For a compliant surface to be effective in delaying transition, it is necessary for both 2D and 3D modes to be brought under proper regulation. Three-dimensional modes also play a leading role in the transition of the boundary layer. Upon gaining sufficient amplitude, oblique wave pairs may interact non-linearly with selected 2D modes to form rapidly amplifying wave triads. The selections are governed by resonance principle. These wave triads are the immediate precursors of the chaotic fluctuations that characterize boundary-layer breakdown to turbulence in boundary layer on a rigid surface. All of these are fairly well understood at least for Blasius boundary layer on a rigid surface. However, the roles of 3D waves in boundary-layer stability, non-linear interaction and wave breakdown over compliant surfaces have yet to be more fully examined. Joslin and Morris [7] applied HerbertÕs secondary instability theory [8] to boundary-layer flow over isotropic and anisotropic plate surfaces. They found that the optimized isotropic wall does not have an adverse effect on the non-linear stability of the flow relative to the rigid wall case. The direct numerical simulation (DNS) of 3D non-linear wave interactions in boundary layers over compliant surfaces was first carried out by Metcalfe et al. [9]. Linear wall boundary conditions were employed in these studies. Their results show that non-linear secondary instabilities could arise and cause the flow to become unstable when it is predicted to be stable by the linear theory. Both of these studies employed the temporal model of instability, which assumes the waves are uniform spatially and grow only with time. The temporal postulation offers simplicity in both mathematical and computational modeling. In temporal DNS, the computational domain is usually limited to just a few streamwise waves under the condition spatial periodicity. This leads to simplified boundary conditions as well as considerable savings in computational cost. However, the temporal approach is subject to a number of limitations when applied to the modeling of wave

1064

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

systems over compliant surfaces. The divergence of the boundary layer implies that the uniformwave assumption is at best an approximation. The temporal wave model cannot handle the strong streamwise inhomogeneities presented by edge boundary conditions of finite length compliant surfaces—practical compliant walls are necessarily finite in length. The typical TS waves on compliant surfaces frequently coexist with CIFI waves. The temporal wave model cannot accommodate multiple waves if these have incommensurate wavelengths. For the above reasons, a more general spatial treatment is needed to study waves on compliant surfaces. Two dimensional spatial–temporal numerical simulations of unstable linear waves on finite compliant surfaces without the above limitations have been carried out by Davies and Carpenter [10] and Wiplier and Ehrenstein [11]. A spatial–temporal DNS model for three-dimensional waves in boundary layers over finite-length compliant surfaces is developed in the present study. The model is needed for investigating 3D secondary instabilities and other nonlinear transition processes in boundary layers over such surfaces. The present 3D model employs the flow variables in their primitive form, whereas Fasel et al. [12] and Rist and Fasel [13] had applied a vorticity–velocity formulation for the governing flow equations. A variant of the vorticity–velocity method was developed by Davies and Carpenter [14] for simulating disturbance evolution in 3D boundary layer over compliant surfaces, and applied to boundary-layer flows over rigid and compliant rotating disks. Unlike the customary 3D vorticity–velocity formulation in which six governing equations are usually required, only three governing equations are solved in their highly efficient method. The present computation uses an iterative fractional step procedure with curvilinear finite-volume discretization on a non-staggered grid system. The buffer domain technique of Liu and Liu [15] is employed to reduce wave reflection from the outflow boundary. The study is restricted here to linear membrane surfaces, which are modeled numerically by a high-order difference scheme. Implicit in the adoption of the wall model is the assumption that the surface waves are small in their amplitude. Consequently, the interaction of the flow and the membrane surface is assumed to be essentially linear although the perturbation dynamics of the flow is fully non-linear. Lucey et al. [16] used a non-linear membrane model in their recent study of largeamplitude static divergence waves. Results on the propagation of 3D linear oblique waves and the non-linear interactions of wave triads leading to subharmonic and fundamental breakdown of the boundary layer are presented. Despite the assumption of linear flow–wall interaction, compliant surfaces are found to have a fairly strong influence on the development of non-linear waves and their breakdown.

2. Fully implicit fractional step method The fractional step method for solving the Navier–Stokes equations was originally developed by Chorin [17]. It was popularized by Kim and Moin [18] and has now become the most widely used approach for DNS. In this paper, the fractional step method is applied to the three-dimensional (3D) incompressible perturbation Navier–Stokes equations. A pressure correction method is used in the time splitting strategy to enforce the divergence-free flow condition. To handle the fluid–structure interaction problem, a fully implicit iterative variant of the fractional step method is employed in our study. The computational domain is presented in Fig. 1.

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095 y

1065

z x Free stream Outflowbuffer domain

Inflow

membrane

ribbon

Fig. 1. Three-dimensional (3D) computational domain of transitional boundary layer over a 2D membrane.

2.1. Governing equations The Cartesian tensor notation is adopted for the 3D incompressible perturbation Navier– Stokes equations: The continuity equation oui ¼ 0; oxi the momentum equations


 ouk o o 1 ouk op þ l ðul uk þ uk ul þ ul uk Þ ¼ l
k l ox ox Re ox ox ot

ð1Þ

ðk ¼ 1; 2; 3Þ:

ð2Þ

In the above equations, Re is the Reynolds number based on the freestream velocity of mean flow (U1), the displacement thickness (dx0) of the mean boundary layer at the location where the ribbon exciter is placed, and the kinematic viscosity (m) of the fluid. The x1, x2 and x3 denote the streamwise (x), wall normal (y) and spanwise (z) coordinates of the Cartesian frame respectively. The velocity components u1, u2 and u3 in the x-, y- and z-directions are used interchangeably with u, v and w, respectively. The u and  u indicate the perturbation and base flow velocity components respectively. 2.2. Time splitting strategy The employment of second-order backward-Euler scheme requires that the following momentum and mass conservative equations, Eqs. (3) and (4), are satisfied at each time step,
 3umþ1
4umk þ um
1 o o 1 oumþ1 opmþ1 k k mþ1 mþ1 mþ1 mþ1 k ; ð3Þ
þ l ðul uk þ uk ul þ uk ul Þ ¼ l ox ox Re oxl 2Dt oxk oumþ1 k ¼ 0; oxk where m + 1 indicates the current time step.

ð4Þ

1066

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

Different time splitting strategies can be developed based on the different ways Eq. (3) can be decomposed. In the present iterative fractional step method, a pressure correction scheme is adopted. Two levels of iterations are carried out at each time step, termed the outer iteration and the inner iteration. Outer iterations between the momentum and the pressure correction equations are performed to fully couple those equations at each time step. Each outer iteration begins with the solving of the momentum equations to obtain the intermediate velocity field uk :
 3uk
4umk þ um
1 o o 1 ouk opðrÞ     k : ð5Þ
þ l ðul uk þ uk ul þ uk ul Þ ¼ l 2Dt oxk ox ox Re oxl Here (r) indicates number of outer iterations. p(r) is the perturbation pressure calculated from the previous outer iteration. At each time step, the initial perturbation pressure p(0) is derived from the perturbation pressure at the previous time level pm. For each outer iteration (r) ! (r + 1), Eq. (5) is iterated a number of times to yield u*, with u(r) as the initial guess value. With u* thus obtained, several iterations are then performed on the pressure correction equation: o op 3 ouk ¼ : oxk oxk 2Dt oxk

ð6Þ

to give correction pressure field p. The intermediate velocity field uk and the pressure correction p are then used to correct the velocity and pressure fields: ðrþ1Þ

uk

¼ uk


2Dt op ; 3 oxk

pðrþ1Þ ¼ pðrÞ þ p:

ð7Þ ð8Þ

At each time step from m ! m + 1, the outer iteration involving Eqs. (5)–(8) is repeated until the convergence criteria for the perturbation velocity and pressure fields are satisfied simultaðrþ1Þ ¼ uk and neously. The field values at the m + 1 time level are then updated as follows: umþ1 k m+1 (r+1) =p . p

3. Numerical algorithm 3.1. Coordinate transformations Tensor notation [19] is adopted to describe the coordinate transformation. In a curvilinear nonorthogonal coordinate system, two different frames of basis vectors are usually defined. The covariant basis vectors are tangent to the coordinate lines: eðiÞ ¼

oxk ik ¼ J ki ik ; oni

while the contravariant basis vectors are normal to the coordinate surfaces.

ð9Þ

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

1067

oni ik ¼ J ik ik : ð10Þ oxk These two vectors are parallel if and only if the curvilinear coordinates are orthogonal. The chain rule connects the derivatives in the curvilinear and the Cartesian coordinate systems: eðiÞ ¼

o/ oxj o/ o/ ¼ J ji j ; i ¼ i j ox on on ox

o/ onj o/ j o/ ¼ i j ¼ Ji j : i ox ox on on

ð11a,bÞ

The covariant, contravariant and mixed metric tensor components are given by, gij ¼ eðiÞ  eðjÞ ¼ J ki J kj ;

i

j

gij ¼ eðiÞ  eðjÞ ¼ J k J k ;

gij ¼ eðiÞ  eðjÞ ¼ dij :

ð12a,b,cÞ

The face area vector of a control volume generated by two of the three vectors is given by AðkÞ ¼ eðlÞ  eðmÞ ;

ð13Þ

where k, l, and m are cyclic. The Jacobian determinant of the coordinate transformation is defined by  ox1 ox2 ox3    1  on on1 on1     ox1 ox2 ox3  i ð14Þ J ¼ detðJ j Þ ¼  on2 on2 on2  ¼ eð1Þ  ðeð2Þ  eð3Þ Þ ¼ eð1Þ  Að1Þ :     1  ox3 ox23 ox33  on

on

on

In general AðiÞ ¼ J eðiÞ

i

or Aik ¼ J J k :

ð15Þ

The above equations connect quantities in the physical and computational domains. If a control volume has unit dimensions in the computational space, the corresponding volume in the physical space has the area A(i) and volume J. A velocity vector U can be expressed in terms of its Cartesian, covariant, or contravariant components as follows: U ¼ ui ii ¼ U 0i eðiÞ ¼ U 0i eðiÞ ;

ð16Þ

while its normal flux through a surface is given by b i ¼ U  AðiÞ ¼ uj Ai : U j

ð17Þ

The relation between the Cartesian velocity components and the normal flux components arises from uj ¼

J ij i b : U J

ð18Þ

3.2. Finite volume discretization Finite volume in a collocated grid system is adopted in this study. Rhie and ChowÕs [20] method is employed to suppress spurious check-board fluctuations in the pressure solutions.

1068

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

3.2.1. Discretization of the momentum equations Compared with the normal Navier–Stokes equations, two more linear convective terms ul uk and ul uk , which represent the coupling between the base flow and perturbation velocities, are to be found in the governing equation (2). The first linear term represents the transport effects of base flow on the perturbation velocity components. The sum of this term and the non-linear term uluk forms the total convective term that is treated implicitly in our calculations. On the other hand, the second linear term ul uk , which represents the transport effects of the perturbation flow on base flow velocity components, is arranged into the source term on the RHS. Therefore, the integral of the momentum equation (5) over a volume dV in the physical plane can be written as,  Z   3uk
4umk þ um
1 o 1 o2 uk    k þ l ðul uk þ uk ul Þ
 dV ox Re oxl oxl 2Dt dV  Z  opðrÞ o 
k
l ðuk ul Þ  dV : ð19Þ ¼ ox ox dV In our simulation, the finite volumes have unit dimension in the computational coordinate frame. Therefore, using the Gauss divergence theorem and Eq. (17), Eq. (19) may be written as:

  5 5 ðrÞ X X 3uk
4umk þ um
1 n n b i j op i k uk Þn ; þ ð
1Þ ðQl Al Þn ¼
Ak j
ð
1Þ ð U ð20Þ J 2Dt on P n¼0 P n¼0 where n = 0, 1, 2, 3, 4, 5 mark the six faces of the control volume e (east), w (west), n (north), s (south), t (top), b (bottom). The superscript index i = int(n/2) + 1, while the subscript P denotes the center of the finite volume. In above equation,    i 1 ouk i i  bi þ Q bi ; b Þu
1 Gij ouk ¼ Q bi þ U ul Þuk
¼ ð U ð21Þ A Ql Al ¼ ðul þ  o no l k j l Re ox Re on where Aik Ajk ; G ¼ J ij

 i b i ¼ ðU b Þu
1 Gii ouk ; bi þ U Q o k Re oni
 i 1 ij ouk b : Q no ¼
G Re onj i6¼j

ð22Þ ð23Þ ð24Þ

At this juncture, a general discretized equation may be obtained by applying a suitable spatial scheme,
 5 ðrÞ X X j op   b i uk Þ ; aP ukP ¼ anb uknb þ b
Ak j þ ð
1Þn ð U ð25Þ n on P n¼0 P P
4um þum
1 ou 1 3J Gij onkj Þi6n¼j
J ð k2Dt k Þ with i = int(n/2) + 1; aP ¼ anb þ 2Dt and anb dewhere b ¼ 5n¼0 ð
1Þn ðRe pend on the selected spatial scheme.

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

1069

The second-order central-difference scheme is used for spatial discretization here. Thus the orthogonal flux at the east and west faces of the physical volume are given by, 1

b Þ ¼ AðP e ÞDe ðu
u Þ þ F e u ; ðQ kP kE kP o e

ð26Þ

1

b Þ ¼ AðP w ÞDw ðu
u Þ þ F w u ; ðQ kW kP kP o w

ð27Þ

where AðP Þ ¼ 1
0:5jP j þ maxf0;
P g; Fn Pn ¼ ; Dn

i

ð28Þ


i

b Þ; b þU F n ¼ ðU n

Dn ¼

Gii Re  Dni

 ;

ðno summation for iÞ:

ð29Þ

n

Therefore, the coefficients aE and aW are given by aE ¼ AðP e ÞDe ;

aW ¼ Að
P w ÞDw :

ð30Þ

3.2.2. Discretization of pressure correction equation The pressure correction equation (6), in the following integral form, Z Z o op 3 ouk  dV ¼  dV k k 2Dt dV oxk dV ox ox

ð31Þ

is also discretized by the finite volume scheme. Applying the divergence theorem, we have
 5 5 X 3 X n ij op b i Þ with i ¼ intðn=2Þ þ 1: ð
1Þ G ¼ ð
1Þn ð U ð32Þ n j on n 2Dt n¼0 n¼0 In the primitive-variable formulation, checkerboard oscillations typically occur in the solution pressure field when the flow variables are defined on a collocated grid (velocities and pressure defined at the same grid location) [21]. Close coupling between the pressure and velocity fields is necessary to suppress this undesirable phenomenon. Many momentum interpolation methods have been developed to achieve this. Among them, the Rhie and ChowÕs [20] interpolation method is the most widely used,
 

 
 i i 2Dt ii op op b b G
: ð33Þ U e ¼ ð U e Þave
i 3 on e ave oni e ave Eq. (33) is used to calculate the convective mass flux at the control volume interface. The subscript ave denotes weighted linear interpolation from neighboring grid points. The second-order central-difference scheme is used to discretize Eq. (32) to yield X anb pnb þ b; ð34Þ aP pP ¼ where aE ¼

G11 e ; Dn1e

aW ¼

G11 w ; Dn1w

aP ¼

X nb

anb ;

ð35Þ

1070

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

and b¼

i6¼j 5 5
X 3 X op ð
1Þn ðU i Þn
Gij j 2Dt n¼0 on n n¼0

with i ¼ intðn=2Þ þ 1:

ð36Þ

3.2.3. Corrections of pressure and velocity fields Finally, the pressure and velocity fields are updated using, pðrþ1Þ ¼ pðrÞ þ p; ðrþ1Þ

uk

¼ uk


2Dt Aik op : 3 J oni

ð37Þ ð38Þ

3.3. Boundary treatment A rectangular coordinate system with grid stretching along the y-axis is utilized in this study, where the geometry is relatively simple. Grid stretching in the y-direction allows us to simultaneously set the freestream boundary far from the wall and to achieve good vertical spatial resolution at the wall. The following transformation functions are used: xðn1 Þ ¼ n1 ; yðn2 Þ ¼

y max cn2 n2max c þ y max ðn2max
n2 Þ

ð39Þ ð40Þ

following Liu and Liu [15]. In above equations, superscript 2 in n2 indicates the wall-normal direction. The perturbation velocities ui are assumed to decay to zero far from the wall. The values of ui at the free stream boundary are hence ui ¼ 0

ði ¼ 1; 2; 3Þ:

ð41Þ

The perturbation velocity is set to zero at the inflow boundary. The flow is excited by a spanwise pulsating line source, of specified frequency content, located on the wall a short distance from the inflow boundary. This is intended to approximate the action of vibrating ribbons used in experiments. The upstream influence of the exciting line source is assumed to be negligible at the inflow boundary. Periodic boundary conditions are applied in the spanwise (z) direction. The perturbation velocities ui = 0(i = 1, 2, 3) at the wall for a rigid surface. The wall boundary conditions for non-rigid membrane surfaces are described in the next section. A buffer domain is introduced at the outflow boundary to reduce upstream wave reflections, following Liu and Liu [15]. Inside the buffer region, two buffer functions are employed, which act on the diffusion terms of governing equations. The first buffer function b(n1) is applied to the streamwise diffusion terms to gradually parabolize the governing equations in the downstream direction. The second buffer function bRe(n1) acts on all the diffusive terms to reduce the Reynolds

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

1071

number in the buffer domain gradually to a value below the critical Reynolds number so as to damp out all the perturbation modes. These treatments yield:  1   1  bðn ÞbRe ðn1 Þ  G11 bðn ÞbRe ðn1 Þ  G11 ; D ¼ ; e Re  Dn1 Re  Dn1 w e     bReðn1 ÞG22 bRe ðn1 Þ  G22 Ds ¼ ; Dn ¼ ; Re  Dn2 s Re  Dn2 n     bReðn1 ÞG33 bRe ðn1 Þ  G33 Db ¼ ; Dt ¼ ; Re  Dn3 b Re  Dn3 t

Dw ¼

ð42Þ

where ( tanhðL 1

bðn Þ ¼

bRe ðn1 Þ ¼

1 total
n Þ tanhðLbuffer Þ

þ 1  10
20 ;

1; 0 6 n1 6 Loriginal 8 < c ðn1
Loriginal Þ 2 þ 1; L :

tanhðLbuffer Þ2

1;

Loriginal 6 n1 6 Ltotal

original

;

6 n1 6 Ltotal

ð43Þ

ð44Þ

1

0 6 n 6 Loriginal :

At the outflow boundary of the buffer domain, o2 ui ¼ 0 ði ¼ 1; 2; 3Þ: ox2

ð45Þ

4. The membrane equation The spring-backed isotropic tensioned membrane with elastic foundation and damping is used in this study. The dimensional governing equation of this compliant surface is,
2   2   o2 g  og o g  o g þ k  g ¼
pw ;
T þ ð46Þ m 2 þ d ot ot ox2 oz2 where g* denotes the dimensional vertical displacement of the compliant surface. m* and T* are the surface mass density and the tension of the membrane respectively; while k* and d* are the elastic and damping constants of the foundation respectively. The superscript * denotes dimensional quantity here. pw is the flow disturbance pressure at wall. The dimensional membrane parameters are non-dimensionalized by the free stream velocity U 1 , fluid density q*, and a reference length scale L*: mðLÞ ¼

m ; q L

d ðLÞ ¼

d ; q U 1

T ðLÞ ¼

T ;  q U 2 1L

k ðLÞ ¼

k  L : q U 2 1

ð47a,b,c,dÞ

1072

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

The properties of the membrane and its foundation are specified in terms of these non-dimensional parameters to permit comparison between different compliant walls. For the purpose of computation, however, the reference length employed is the displacement thickness d0 of the mean boundary layer at the location of the exciter. The properties of the membrane and foundation in the two length scales are related by: ReðLÞ ReðLÞ Re ; ð48a,b,c,dÞ ; d ¼ d ðLÞ ; T ¼ T ðLÞ ; k ¼ k ðLÞ m ¼ mðLÞ ReðLÞ Re Re where ReðLÞ ¼ L  U 1 =m and Re ¼ d0 U 1 =m are the Reynolds numbers based on the membrane reference length scale L* and the computational (flow) reference length scale d0 respectively. Finally, the non-dimensional membrane equation in the computational length scale is:
2  o2 g og o g o2 g þ ð49Þ þ kg ¼
pw : m 2 þd
T ot ot ox2 oz2 The finite membrane spans (xcs, xce) in the streamwise direction. Hence at the two ends of the membrane, the non-dimensional displacements are gðxcs Þ ¼ gðxce Þ ¼ 0:

ð50Þ

Periodic boundary conditions are applied in the spanwise direction for consistency with the periodic character of the perturbation in that direction. The perturbation fluid velocity (uw, vw, ww) at the wall is coupled to the displacement of the surface by ( 0; x 6 xs or x P xe ; ð51Þ uw ¼ ou
g oy ; xs < x < xe ( vw ¼

0;

x 6 xs or x P xe

og ; ot

xs < x < xe

ww ¼ 0;

;

ð52Þ ð53Þ

which may be derived via linear approximation. Finally, the backward-Euler scheme is adopted for the temporal integration of the membrane equation, while fourth-order accurate central-difference scheme is employed for spatial discretization. The use of fourth-order method for structure is not strictly necessary, but we find it improves the convergence for the coupled problem.

5. Computational implementation In our DNS project, geometric multigrid (GMG) algorithm [22] is selected and applied for solving the discretized pressure correction equations. A 3D multigrid Poisson solver is written and applied to the iteration of Eq. (6). Its application greatly accelerates the overall speed of the unsteady computation. 3D ADI solver is used as the smoother. To save memory, flow variables, grid parameters and coefficients of the discretized equations are stored as 1D arrays in the main program and re-declared as 3D arrays in the subroutines, where most of the operations are

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

1073

performed. The size of these 1D arrays and the starting addresses of the 3D arrays for each grid level must be carefully calculated in order to ensure the correct mapping of variables between the main program and the subroutines. In our simulation, the computational domain is uniformly partitioned into 32 blocks in the streamwise direction. After domain decomposition, a multilevel grid is constructed inside each block. Interface communications between adjacent blocks is facilitated by the employment of ghost volumes. One level of overlap is used in the present study. MPI [23] subroutines are employed to handle boundary data transfers. New data types are constructed to reflect the data structures at the six different interfaces of each block. Using these new data types, communication time among different CPUs is tremendously reduced. Consequently, the scalability of current code is effectively increased. Finally, the computational procedure for fluid–structure interaction comprises seven sub-steps: (1) Guess values for all the variables; such as u, v, w, p etc. are obtained from the solutions at the previous time step. (2) The boundary conditions are updated and velocities at the wall (uw, vw, ww) are calculated from the membrane displacement Eqs. (51)–(53). (3) The momentum equations (Eq. (5)) are relaxed a fixed number of times to obtain uk from the pressure field p(r). (4) The pressure-correction equation (Eq. (6)) is iterated upon a fixed number of times to obtain p. (5) The pressure, mass flow rate and Cartesian velocity components are corrected. (6) The membrane equation is iterated upon a fixed number of times to obtain the membrane displacement. (7) These new values are used in a new outer iteration cycle. Steps (2)–(7) are repeated until convergence criteria for flow variables and wall displacements are satisfied at each time step. Around 14–16 outer iterations are usually needed at each time step to reduce the L2 residual errors by 10
2. The final relative L2 residuals of Eqs. (25), (34) and (49) are of the order 10
8 or smaller at each time step. Each run takes about 26.6–60 h of CPU time on a cluster of 32 Pentium III 1.0 GHz processors with 1 GB RAM per processor, depending on the grid resolutions and the total number of time steps. The curvilinear finite-volume formulation and domain decomposition feature of the present code combine to yield a geometrically flexible code with good conservation properties. As the computational geometry involved in the present study is relatively simple, involving only coordinate stretching in one direction, not all of these useful features are fully exploited here.

6. Computational results The solver described in the preceding sections is first validated by comparing with existing rigidwall results. The linear wave computation of Fasel et al. [12] and non-linear growth data from Saric et al. [24]Õs experiments, as presented by Liu and Liu [15], are used in this comparison. Several cases of 3D linear wave evolution over compliant membranes are studied thereafter. Finally,

1074

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

the secondary instabilities and the late stages of subharmonic and fundamental breakdown are investigated. All computations employ 32 computational streamwise blocks and time step equal to 1/320 of the TS wave period, unless otherwise indicated. 6.1. Code validation The 3-D waves in the present linear study are driven by a sinusoidal spanwise line source of the following form located at the wall: u2 ¼ A3D sinðb0 z
x0 tÞ;

at x ¼ x0 ;

y ¼ 0;

ð54Þ

where b0 and x0 are the spanwise wavenumber and the wave frequency respectively. The specific case of a 3D wave driven by a source with b0 = 0.10994 and x0 = 0.054971 (or F = 100) placed at x0 = 185.644 (Re = 549.71) is examined by Fasel et al. [12] in some detail. This case is simulated here by a computational domain that spans x = [124.430, 600.6828], y = [0, 52.76] and z = [
28.575, 28.575]. A total of 513 · 65 · 33 grid points are used, with the stretch parameter in the normal direction being set at c = 1.8. A small source amplitude of A3D = 10
4 is used to ensure that the amplitude of the excited 3-D wave is within the linear range, and simulation is carried out for 15 TS-wave periods. Figs. 2–4 compare our computational results with corresponding results scaled from Fasel et al. [12]. Fig. 2 shows the three components of velocity perturbations at the location Re = 770. The agreement is excellent throughout the range, including the zero point of u-perturbation. Fig. 3 shows a comparison of local streamwise wavenumber ar with Fasel et al. and linear stability values. Our results (at y = 2.25) track the linear stability curve closely over the full range of Reynolds number up from Re = 650. The oscillation in the results below Re = 650 is due to evanescent waves in the proximity of the perturbation source, which is also evident in Fasel 1 0.9 u

0.8

Amplitude

0.7 0.6 0.5 w

0.4 0.3 0.2

v

0.1 0

0

2

4

6

8

10

y (in current scale)

Fig. 2. Amplitude distributions for the three disturbances velocity components at Re = 770 on a rigid flat plate (present results ——,– – – –, –  –  –; Fasel et al. [12] s, h,
).

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095 0.015

Linear stability theory umax (present study) umax (Fasel et al. [12]) vmax (present study)

0.01

vmax (Fasel et al. [12]) wmax (present study) wmax (Fasel et al. [12])

0.005

αi

1075

0

-0.005

-0.01

-0.015 600

700

800

900

Re

Fig. 3. Comparison of growth rates ai with Fasel et al. [12]Õs results and linear stability theory.

0.18

0.17

Linear Stability Theory present study Fasel et al 1991

αr

0.16

0.15

0.14

0.13 550

600

650

700

750

800

850

900

950

Re

Fig. 4. Comparison of wavenumber ar with Fasel et al. [12]Õs results and linear stability theory.

et al.Õs result. Fasel et al.Õs result is based on 7 TS-wave cycles; their data above Re = 800 are not applicable to the comparison as they reflect the remnants of the leading wave packet. In our longer computation, the leading waves have passed out of the computational domain and our results are in general agreement with the linear stability curve up to Re = 900. Comparison of local spatial growth rates ai is given in Fig. 4. For reasons noted above, comparison with Fasel et al. is

1076

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

valid only in the range of Re from 650 to 800. There is generally good agreement between our results and those of Fasel et al. for both ar and ai. We also note that the numerically derived values of ar and aiare fairly sensitive to the spatial distribution of wave amplitude, whose accurate determination in simulation generally requires a significant passage of time for the initial transient to die away. Numerical differentiation to obtain ar and ai also tend to amplify any errors. Some of the observed differences between our results and those of Fasel et al. may be due to these sensitivities. Differences from linear stability prediction could of course be due to non-parallel effects not accounted for in the calculation of eigenvalues. The solver is also verified for non-linear wave simulation by comparison with the numerical results of Liu and Liu [15] and the experimental results of Saric et al. [24] for the fundamental breakdown of a resonant wave triad. The propagating wave triad is generated by a spanwise line source at the wall having the form of u2 ¼ A2D sinð
x2D tÞ þ A3D sinðb0 z
x3D tÞ þ A3D sinð
b0 z
x3D tÞ

at x ¼ x0 ; y ¼ 0;

ð55Þ

where A2D and A3D denote the excitation amplitudes of the primary 2D wave and 3D oblique wave pair respectively. The x2D and x3D are the corresponding excitation frequencies of the 2D and 3D waves. b0 indicates the spanwise wavenumber of the 3D mode. For the present case, the wave source is located at Re = 1221.77 (x0 = 412.609) with x2D = x3D = 0.0928 and spanwise wavenumber b0 = 0.02451. The computational domain extends over x = [327.559, 854.447], y = [0, 50] and z = [
25.635, 76.906] and is covered by a grid with 1025 · 65 · 33 node points and c = 1.8. The excitation wave amplitudes of A2D = 0.0181 and A3D = 0.00129 are selected after extensive pre-simulation trails to give perturbation wave amplitudes at excitation location similar to those of Liu and Liu and Saric et al., and the simulation is carried out for 20 TS-wave cycles. The u-perturbation amplitude along the streamwise direction at y = 0.4 is compared with the results of Saric et al. [24] and Liu and Liu [15] in Fig. 5. Good overall agreement in both the magnitude and amplification rate of the u-perturbation amplitude is obtained. The small deviation is due at least in part to the difficulty of exactly matching the initial 2D and 3D amplitude of the non-linear waves produced by the different methods of wave generation involved in the three works. The above simulations of 3D linear wave evolution and non-linear wave breakdown on a rigid wall affirms the veracity of the current code. 120 100 Liu & Liu, [15] Saric et al., [24] Present study

u'/u'0

80 60 40 20 0 600

700

800

900

1000

1/2

(Re*x)

Fig. 5. Comparisons of u amplitude along streamwise direction.

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

1077

6.2. 3D linear waves over finite membrane surfaces Six cases of linear 3D wave propagation over membrane surfaces are examined in the present study. Cases 1–4 are performed to show the effects of membrane parameters on linear oblique waves. Cases 5 and 6 are specially selected to show that under certain conditions a 3D wave can grow faster than its 2D counterpart of the same frequency. The parameters for the six cases are presented in Table 1(Panels A–C). In Table 1(Panel A), xin, xout, ymin, ymax, zmin and zmax indicate the boundaries of the computational domain. The xr denotes the location of the wave source, while xcs and xce mark the locations of leading and trailing edges of the finite membrane. The perturbation parameters are tabulated in Table 1(Panel B). These are scaled based on the displacement thickness of the boundary layer at the location (xr) of the perturbation source. Table 1(Panel C) summarizes the properties of the surfaces employed in the six cases. A grid of 513 · 65 · 33 is used for Cases 1–4. The evolving waves are simulated for 20 TS-wave cycles. The long simulation is essential for one to capture the complete wave system sustained by the compliant membrane, including the full effects of wave reflection/generation at its leading and trailing edges. In this regard, direct spatial numerical simulation will be seen to offer a more complete evaluation of the true stability performance of a real compliant surface than could be deduced from eigenvalue studies. Table 1 Parameters of the computational domain (Panel A), perturbation source (Panel B) and properties of the membrane surface (Panel C) Case

Computational domain xr

Panel A 1–4 135.08 5–6 1281.30

xin

xout

ymin

ymax

zmin

zmax

xcs

xce

56.56 1246.67

847.65 1731.49

0 0

>50 >50


62.83
41.40

62.83 41.40

203.4 1306.3

700.6 1671.7

Panel B Perturbation parameters 1–4 5 6

Re

b0

x2D

x3D

A2D

A3D

400 3794 3794

0.05 0 0.07588

0 0.01897 0

0.0624 0 0.01897

0 0.0004 0

0.00004 0 0.0004

Panel C Membrane parameters m(L)

T(L)

d(L)

k(L)

Re(L)

1 2 3 4

1 1 1 1

10 20 10 10

0.1 0.1 0.1 0.5

0 0 1 0

580 580 580 580

5 6

1 1

4 4

0.02 0.02

0.1 0.1

580 580

1078

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095 0.15

0.12

ω

0.09

0.06

0.03

0

0

500

1000

1500

2000

Re

Fig. 6. Neutral curves for computational cases. Rigid wall (  ), Case 1 (——), Case 2 (– – –), Case 3 (–    –) and Case 4 (–  –).

The neutral curves of Tollmien–Schlichting instability for Cases 1–4 and the rigid wall are presented in Fig. 6. The u-perturbation velocity of Case 1 over the compliant membrane at y = 0.3 is plotted in Fig. 7(a). Besides the TS waves, CIFI waves could also be generated inside the boundary layer according to linear stability theory. These CIFI waves usually possess a different streamwise wavenumber and therefore travel at a different oblique angle than the TS wave. The waves are responsible for the mild wavy modulation of the primary TS wave that may be detected in its crests in Fig. 7(a). The same modulation of the wave crests is even more clearly evident in the surface displacement plot of the membrane surface in Fig. 7(b). This is because the amplitude of CIFI waves, whose existence is associated with the presence of the compliant membrane, tends to be relatively larger at the compliant surface. The wave crests for the corresponding oblique TS wave over a rigid surface (Fig. 8), on the other hand, are straight and free from such modulation. In Fig. 9(a), a comparison of the u-perturbation velocity of Case 1 and corresponding rigid-wall result is made. The rigid-wall wave has an amplitude minimum at x = 343 (Re = 637.4) and a maximum at x = 476 (Re = 750.9), which correspond to its lower and upper neutral stability branches respectively as given in Fig. 6. While the same minimum and maximum point (associated with the neutral TS branches) also exist for the membrane surface, it is more difficult to pinpoint precisely their locations along the membrane surface because of amplitude modulation by co-existing CIFI waves. There is a clearly amplification of the TS wave between the two turning points. Thereafter, the TS wave decays monotonically over the compliant membrane, so that its amplitude is significantly smaller than the same wave over the rigid wall. This behavior of the wave over the compliant membrane is consistent qualitatively with the smaller unstable region indicated by its neutral curve (compared to the rigid wall) in Fig. 6, and it shows that 3D oblique TS wave can be stabilized by this finite membrane despite the presence of CIFI. Fig. 9(b) presents the wavenumber spectrum of membrane displacement for Case 1. Three peaks are observed, which are

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

Flow

dire

ctio

1079

n

(a) 2E-05 1E-05 100

0 u

200 300

-1E-05

400 500

X

600

-2E-05 -50

700

Z 0 50

800

Flow

dire

ctio

n

(b)

1E-05 100

η

0

200 300 400

-1E-05 -2E-05 -50

z 0 50

500 600

x

700 800

Fig. 7. (a) Perturbation velocity u at y = 0.3 (Case 1). (b) Wall displacement (Case 1).

Flow

n

ctio

dire

2E-05 1E-05 100

0 u

200

-1E-05 -2E-05 -50

Z 0 50

400 600

500

300

X

700 800

Fig. 8. Perturbation velocity u at y = 0.3 over a rigid wall for perturbation parameters corresponding to Case 1.

annotated as TSI, CIFI(1) and CIFI(2). The real wavenumbers for the TSI and CIFI(1) in Fig. 9(b) are 0.161 and 0.030. These are in reasonably good agreement with the real part of eigenvalues calculated from linear stability theory (aTSI = 0.1639 + 0.0138i and aCIFI(1) = 0.0349 + 0.0009i). The poorer match between the spectrum-analyzed wavenumber and the stability eigenvalue (real part only) for the CIFI(1) wave may be attributed to the fact that the compliant membrane contains only a few CIFI(1) wave periods and the wavenumber of the long CIFI(1) wave may vary

1080

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095 5E-06 4E-06 3E-06 2E-06 1E-06 0 -1E-06 -2E-06 -3E-06 -4E-06 -5E-06

u

(a)

200

1

(b)

400

x

600

800

TS I

FFT(η)

0.8 0.6

CI FI (1 )

0.4 CI FI (2 )

0.2 0 0

0.25

0.5

0.75

1

Fig. 9. (a) Comparisons of perturbation velocity u of Case 1 (——) and rigid wall (– – –) (y = 0.3, z = 0). (b) Wavenumber spectrum of wall displacement of Case 1 at z = 0.

(albeit slowly) along the membrane. The CIFI(2) peak in Fig. 9(b) may be due to evanescent (damped) waves generated at the edges. Its eigenvalue has not been identified due to the lack of sufficiently accurate guess value. In Cases 2 and 3, we consider the effects increased membrane tension T(L) and foundation stiffness k(L) have on the oblique TS wave respectively. In Fig. 10(a), we can see that the region of spatial amplification (unstable region) is now increased as a consequence of an increase in the membrane tension from T(L) = 10 to 20. Furthermore, the rate of wave decay after the maximum amplification point is also slower for Case 2. These indicate that increased membrane tension is destabilizing for TS waves; which is consistent with the neutral stability curves presented in Fig. 6 for the two cases. At the level of the wall (y = 0), Fig. 10(b) shows that the membrane displacements associated with the primary TS wave and the modulating CIFI mode are reduced by the increase in membrane tension. The same type of effect may also be seen when the foundation stiffness is increased from k(L) = 0 (Case 1) to 1.0 (Case 3) in Fig. 11; i.e. the enlargement of the region of positive spatial amplification (Fig. 11(a)), and the reduction in membrane displacement (Fig. 11(b)). The spatial growth of the TS mode, marking its destabilization, is also clearly evident in the membrane displacement for Case 3. Significantly, the modulating effect by the CIFI, which is so marked for Case 1, is entirely absent from the membrane displacement for Case 3. The absence of any distinct CIFI mode for Case 3 is confirmed by the wall displacement spectrum presented in Fig. 11(c). This shows that an increase in the stiffness of the foundation may be able to suppress the CIFI or wall-related modes. The above results also show that although increase in membrane tension and/or foundation stiffness (both resulting in a stiffer surface) may result in more rapid spatial growth (destabilizing) for the primary TS mode, they actually reduce the

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095 (a)

u

5E-06 4E-06 3E-06 2E-06 1E-06 0 -1E-06 -2E-06 -3E-06 -4E-06 -5E-06

1081

200 4E-06 3E-06

400

x

600

800

600

800

(b)

2E-06 η

1E-06 0

-1E-06 -2E-06 -3E-06 -4E-06

200

400

x

Fig. 10. (a) Comparisons of perturbation velocity u of Case 2 (——) and Case 1 (– – –) (y = 0.3, z = 0). (b) Wall displacement of Case 2 (——) and Case 1 (– – –) (z = 0).

amplitude of the wave sustained at the surface. Conversely, this means that a soft TS-stabilizing compliant membrane may suffer from having large amplitude surface waves; which may cause it to become more susceptible to non-linear flow–surface interaction effects, even in the absence of CIFI. Finally, the effect of wall damping is considered. Fig. 12(a) and (b) show the effects increased wall damping has on the primary TS wave at y = 0.3 in the flow and on the membrane surface. The increase in damping produces a small increase in the amplification rate of the TS wave and thus damping may be regarded as being destabilizing. This behavior is also in general accord with the neutral curves presented in Fig. 6. The CIFI-modulation is also reduced by the damping (Fig. 12(c)), indicating that damping helps to suppress the CIFI modes. According to SquireÕs theorem [25], the most unstable temporal mode for a rigid wall is twodimensional, and that to any unstable 3D mode, there is an unstable 2D mode at a lower Reynolds number. SquireÕs theorem is not valid, however, if the bounding surface is non-rigid; the reason being that a compliant surface will appear stiffer to an oblique wave (than to the associated 2D wave with b = 0) because of the reduced mean flow that the oblique wave rides upon [6]. A 2D TS wave on a compliant surface thus experiences a stronger stabilization effect than corresponding oblique TS waves. Indeed, for a soft surface, the critical Reynolds number may belong to a 3D TS mode. Consequently there will be situations in which 3D oblique TS waves may grow faster temporally or spatially than their 2D counterparts. An example is given in Fig. 13(a) for a 2D (b0 = 0) and 3D (b0 = 0.0759) TS wave modes driven at the same frequency and amplitude from the same perturbation source location. For these cases, the time step is set equal to 1/640 TS period and

1082

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095 5E-06 4E-06 3E-06 2E-06 1E-06 0 -1E-06 -2E-06 -3E-06 -4E-06 -5E-06 u

(a)

4E-06 3E-06 2E-06 1E-06 0 -1E-06 -2E-06 -3E-06 -4E-06

200

400

200

400

x

600

800

600

800

FFT (η)

η

(b)

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

(c)

0

x

TSI

0.2

0.4

α

0.6

0.8

1

Fig. 11. (a) Comparisons of perturbation velocity u of Case 3 (——) and Case 1 (– – –) (y = 0.3, z = 0). (b) Wall displacement of Case 3 (——) and Case 1 (– – –) (z = 0). (c) Wavenumber spectrum of wall displacement of Case 3 at z = 0.

spatial simulation is halted shortly after 5 wave cycles to prevent excessive interference from coexisting CIFI modes and waves reflected/scattered from the edges of the finite membrane. Fig. 13(a) shows that the 2D TS wave initially decays. It then propagates at almost constant amplitude before appearing to amplify near the end of the Re range. The 3D oblique mode, on the other hand, amplifies spatially over the whole Re range. Despite careful attempt to minimize interference from other waves through selection of parameters, there is still the hint of a co-existing long CIFI wave in the 2D wave profile, as marked by its descending mean value. The behavior of the 2D and 3D modes is consistent with the lower neutral curves for these modes in Fig. 13(b). At this juncture it may be appropriate to highlight the role played by the edges of the finite membrane as revealed through our simulation of the linear waves. To begin with, the leading edge of the membrane is a principal source of CIFI or wall-related waves that are observed in the simulations. The co-existing CIFI that we had referred to earlier are usually triggered as the primary TS wave, generated by the upstream line source, passes over the leading edge. In the simulations,

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095 5E-06

1083

(a)

4E-06 3E-06 2E-06 u

1E-06 0

-1E-06 -2E-06 -3E-06 -4E-06 -5E-06

4E-06

200

400

200

400

x

600

800

600

800

(b)

3E-06 2E-06

η

1E-06 0

-1E-06 -2E-06 -3E-06 -4E-06

1

(c)

0.9

x

TSI

0.8

FFT(η)

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0

0. 2

0. 4

α

0. 6

0. 8

1

Fig. 12. (a) Comparisons of perturbation velocity u of Case 4 (——) and Case 1 (– – –) (y = 0.3, z = 0). (b) Wall displacement of Case 4 (——) and Case 1 (– – –) (z = 0). (c) Wavenumber spectrum of wall displacement of Case 4 at z = 0.

the leading edge is a point of sharp property changes from a rigid to a compliant surface. Subject to the forcing action of the passing TS wave, the leading edge responds by generating its own waves. These waves are associated with the dynamic quality of the compliant surface (modified by presence of the fluid), and are distinct from the surface deflection that reflects the direct driving action of the TS wave over the surface of the membrane. They are the coupled flow–wall eigenmodes that we have termed CIFI, or FISI by Carpenter and co-workers. In this process, the leading edge facilitates the transfer of energy from the TS wave to the CIFI wave. Because of the Class A or negative-energy character of the TS wave, the TS wave may sometimes experience a local increase in its amplitude as it passes over the leading edge. Energy transfer is usually not possible between linear waves because of their mutual independence. The singular nature of the leading edge (boundary condition), however, allows local coupling between the TS and CIFI waves to take place. A CIFI wave thus produced at the leading edge may amplify or decay as it propagates forward, depending on the properties of the compliant membrane. In our simulations, the CIFI wave typically travels ahead of the TS wave because of its larger wave speed. Upon arriving at the

1084

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095 0.0002

(a) 0.00015 2D

0.0001

3D oblique

u

5E-05 0

-5E-05 -0.0001 -0.00015 -0.0002 1300 0.03

1400

x

1500

1600

1700

(b)

0.02

2D

ω

3D oblique

0.01

0 3000

Re

4000

5000

Fig. 13. (a) Perturbation velocity u at y = 0.3 and z = 0; (b) Lower branches of neutral curves; Case 5 (—–) and Case 6 (– – –).

trailing edge of the membrane, the ensuing interaction between incident and reflected/scattered CIFI waves may result in a quasi-standing or a more complex wave system forming on the membrane. The interaction between the incident TS wave and reflected waves from the trailing edge may also result in large amplitude deflections of the membrane surface near the trailing edge. These may spread/propagate upstream to cover the entire compliant surface when conditions for linear wave coalescence/resonance are satisfied, leading to absolute instability such as those seen in Wiplier and EhrensteinÕs [26] 2D wave simulation. It is difficult to fully anticipate and quantify the wave effects due to the edges without direct simulation. Over the central portion of the compliant membrane, the behaviour of the TS waves conforms broadly to linear stability prediction, provided it is not masked by the simultaneous presence of strong CIFI modes. These aspects of 3D linear waves, not unexpectedly, are similar to what one would see for 2D waves. 6.3. Secondary instabilities and wave breakdown Four cases (Cases 7–10) of wave triad (see Eq. (55)) interaction are simulated to investigate the occurrence of secondary instabilities and non-linear wave breakdown over membrane surfaces.

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

1085

The finite membrane chosen for our study is an optimized membrane (m(L) = 12.0, T(L) = 6.96, d(L) = 0.0, k(L) = 0.0 and Re(L) = 661) obtained from the temporal DNS study of Metcalfe et al. [9]. The details for the computational domain and perturbation source are summarized in Table 2. The computational grid is set at 513 · 65 · 65 with a vertical stretch factor of c = 1.8, and simulations are carried out for 10 cycles of the primary 2D TS waves. Cases 7 and 8 are designed to study the occurrence of secondary instabilities. The excitation amplitude A3D of the oblique wave pair is chosen to be much smaller than the amplitude A2D of the primary 2D wave. The focus of these cases is therefore on the development of small 3D oblique waves growing on the 2D primary waves over a compliant surface. Case 7 examines the case for subharmonic instability. The spatial development of the primary 2D wave mode (1, 0) and leading subharmonic 3D wave mode (1/2, 1) over the compliant membrane are plotted in Fig. 14 together with corresponding results for a rigid surface. Here, we have followed the mode designation scheme of Fasel et al. [12] where (n, k) denotes a mode whose frequency is nx2D and spanwise wavenumber kb0. Unlike the smooth and monotonic increase in amplitude of the 2D primary wave on a rigid surface, the amplitude of the 2D primary mode (1, 0) for the compliant membrane displays an oscillatory behavior. The oscillatory behavior of the (1, 0) amplitude curve may be explained by the interference between the primary 2D TS wave (a2DTSI = 0.2247 + 0.001392i) and a co-existing CIFI wave (a2DCIFI = 0.1291
0.0031869i) of small amplitude. The real wavenumber difference Dar  0.0956 yields about 2.7 wave periods, which correspond roughly to what is observed in Fig. 14. Such oscillation of the amplitude curve is also apparent in the study of 2D linear waves over viscoelastic compliant layers by Wang [27]. A similar but milder oscillation may also be perceived in the amplitude curve for the subharmonic oblique (1/2, 1) mode for the compliant membrane. Putting CIFI-related modulation aside, more importantly we note that the growth of the subharmonic (1/2, 1) mode over the compliant

Table 2 Parameters of the computational domain (Panel A), perturbation source (Panel B) for nonlinear wave simulations Case

Computational domain xr

xin

xout

ymin

ymax

zmin

zmax

xcs

xce

Panel A 7 8

247.20 303.94

223.09 279.20

472.36 566.20

0 0

>50 >50


12.99
31.42

38.98 31.42

260.0 318.9

441.6 530.0

9 10

247.20 303.94

223.09 279.20

472.36 566.20

0 0

>50 >50


12.99
31.42

38.98 31.42

260.0 318.9

441.6 530.0

Panel B Perturbation parameter Re

b0

x2D

x3D

A2D

A3D

7 8

732 900

0.2418 0.1

0.0909 0.0774

0.04545 0.0774

0.0372 0.0372

0.000124 0.000124

9 10

732 900

0.2418 0.1

0.0909 0.0774

0.04545 0.0774

0.0372 0.0372

0.0124 0.0124

1086

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095 10-1

mode (1,0)

|u|

10-2

10-3

10-4

10-5

mode (1/2,1)

300

350

400

x

Fig. 14. Streamwise development of the Fourier amplitudes of selected wave modes at y = 0.26 during subharmonic secondary instability (membrane ——, rigid wall –  –  –).

membrane is significantly smaller than the growth of the same mode over the rigid surface. This is despite the fact that the average amplitudes of the primary 2D waves are about the same over both surfaces. In Case 8, we turn our attention to the spatial growth of fundamental oblique modes. The results are summarized in Fig. 15. Here again we note the smooth and steady growth in the amplitude of the 2D primary (1, 0) mode on the rigid wall and the modulated amplitude behavior of the compliant wall curve. The fundamental oblique mode (1, 1) for the rigid wall also sees a steady if accelerating growth in amplitude, while the same mode for the compliant membrane experiences a much slower average rate of amplitude growth over the computed range. However, compared to the earlier subhamonic (1/2, 1) mode of Case 7 (Fig. 14), the fundamental secondary (1, 1) mode here possesses a distinctly stronger oscillatory component with about four wave periods. This can again be explained on the basis of wave interference between the (1, 1) mode (a3DTSI = 0.1899 + 0.001678i) and a coexisting 3D CIFI mode (a3DCIFI = 0.06470
0.0002173i) over the membrane surface. The wavenumber difference Dar  0.1252 yields a wavelength of 50.2, which represents about 4 waves over the domain shown. The above results show that the favorable or stabilizing effect of surface compliance on linear TS waves also extends to small-amplitude 3D secondary waves riding on 2D primary TS waves of significant amplitude. They reaffirm similar results of previous theoretical study by Joslin et al. [7] and computational study by Metcalfe et al. [9]. Unlike the present work, both of these studies employed the simpler temporal mode model. While such models are well suited for investigating the basic properties of secondary instabilities, they could not be used to provide an accurate assessment of the performance of the more practical walls, which are necessarily of finite extent and possess edges that could reflect or scatter on-coming waves.

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

1087

10-1 mode (1,0)

|u|

10-2

10-3

10-4

10-5

mode (1,1)

350

400

450

500

x

Fig. 15. Streamwise development of the Fourier amplitudes of selected wave modes at y = 0.26 during fundamental secondary instability (membrane ——, rigid wall –  –  –).

In the next two cases (Cases 9 and 10), we will investigate the behavior of the wave triads when the amplitude of the participating oblique wave modes are a significant fraction of the primary 2D wave. The large amplitude of the oblique waves, a hundred folds increase from that in Cases 7 and 8, is chosen to greatly enhance wave–triad interaction. The purpose is to approximately model events that one may expect to see during the late stages of subharmonic and fundamental secondary instabilities as well as the wave breakdown that follows. The same finite membrane as Cases 7 and 8 is used here and linear flow–membrane interaction is assumed for simplicity. Fig. 16 shows the spatial–temporal development of the subharmonic wave triad over a rigid wall and the compliant membrane in term of the u-perturbation velocity contours at y = 0.46. Distinct patterns of staggered K vortices that characterize subharmonic wave instability are formed as the perturbation propagates downstream. The streamwise interval between the K-vortices is seen to be larger on the compliant membrane. Furthermore, the K vortices become evident earlier in its development over the rigid wall, suggesting that the secondary waves associated with its appearance are more highly developed relative to its counterparts on the compliant membrane. The larger streamwise interval of the K vortices over the compliant membrane is found to be linked to the longer wavelength of its primary 2D wave. The further consequence of this is the more stretched or elongated appearance of the K vortices over the compliant membrane. Comparison of contour patterns for the two surfaces at time t = 6T and 8T, where T is the period of the primary 2D wave, shows quite clearly the greater intensity of the K vortices propagating over the rigid wall. Fig. 17 shows the amplitude growth of the primary (1, 0) and principal subharmonic (1/2, 1) wave components over the two surfaces. The amplitude curve for the primary 2D mode over the rigid wall is no longer the smooth ascending curve that it was in Fig. 14. This is clearly due to the stronger non-linear interactions that occur between the primary and subharmonic modes

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095 30

30

20

20

Z

Z

1088

10

10

0

0

-10

-10

X

400

(a)

t=2T

30

30

20

20

Z

Z

300

10

0

-10

-10 400

(b)

t=4T

30

30

20

20

Z

Z

X

10

0

-10

-10 400

(c)

t=6T

30

30

20

20

Z

Z

X

10

0

-10

-10

X

400

(d)

300

X

400

300

X

400

300

X

400

10

0 300

400

10

0 300

X

10

0 300

300

t=8T

Fig. 16. Perturbation velocity u contours at different times during subharmonic breakdown over a rigid wall (left) and the compliant membrane (right).

in the present case, where a much larger subharmonic excitation has been applied. For the same reason, the variation of the 2D primary (1, 0) wave amplitude over the compliant membrane also differs somewhat from that of the earlier subharmonic case (Case 7, Fig. 14). Our greatest interest, however, concerns the principal subharmonic component (1/2, 1). Here, we observe that the growth of the subharmonic mode over the rigid wall is considerably stronger than its growth over the compliant membrane, thus showing that the membrane has a stabilizing effect on the principal subharmonic modes even when the latter have acquired non-linearly significant amplitude. This amplitude analysis of the results substantiates what we had noted qualitatively of the development of the K-vortices (u-velocity contours) over the two surfaces in Fig. 16. Fig. 18 depicts the evolution of the perturbation on the compliant membrane itself. The leading 2D TS wave has arrived at the trailing edge of the membrane at about t  4T. At the trailing edge of the membrane, the wave is reflected. At about t  4T, subharmonic (3D) waves are just becoming visible as small-amplitude spanwise modulation of the second and third wave crests. The small raised humps along the wave crests seem to correlate with the locations of the K-vortices in Fig. 16. A graphic animation of the process shows that these (modulated structures) propagate downstream with the forward motion of the wave crests, until they eventually interact with reflected

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

1089

10-1

mode (1,0)

10

-2

|u|

mode (1/2,1)

10-3

10-4

10-5

300

350

400

x

Fig. 17. Streamwise development of the Fourier amplitudes of selected wave modes at y = 0.26 during subharmonic breakdown (membrane ——, rigid wall –  –  –).

waves towards the end of the membrane. On the whole, the waves on the membrane are predominantly two-dimensional. Three-dimensional or spanwise variation of surface displacement becomes more prominent towards the downstream half of the membrane, but the variation remains relatively weaker than the 2D wave. The maximum membrane displacement is less than 1.3% of the local boundary-layer thickness based on 0.995U1. We hence believe that the linearized wall boundary conditions that we had used do not introduce significant errors into these results. The fundamental breakdown of the wave triad over the membrane surface is now considered. Fig. 19 depicts the spatial–temporal development of the wave triad over the compliant membrane and a rigid wall in term of the u-perturbation velocity contours at y = 0.46. Instead of staggered K-vortices, streamwise aligned K-vortices are formed on both surfaces. From the distribution of the contours, we can see that the developing 2D primary wave over the membrane possesses larger amplitude than its counterpart over the rigid wall. At t = 5T, the first K-vortex becomes visible on both surfaces; with the vortex on the membrane occurring a little further downstream. The K-vortex on the compliant membrane is clearly of stronger intensity. The second K-vortex is formed by t = 6T (a third may also be seen on the rigid wall). At this time, preliminary signs of vortex breakdown are already evident in the second K-vortex on the compliant wall and the third K-vortex on the rigid wall. The breakdown appears to be slightly more advanced on the compliant membrane (which may have prevented the formation of a third K-vortex). General breakdown of the wave structure into small fluctuations occurs thereafter (see t = 7T) on both surfaces. Fig. 20 shows the development of wave amplitudes over the two surfaces. Similar to the subharmonic case (Case 9, Fig. 17), the amplitude of the 2D primary waves (1, 0) for both walls exhibits significant non-linear modification (compare with the small-amplitude oblique-wave case in

1090

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

275 300 325

η

0

350 375

-0.05

Z

X

400

0 20

425 40

450

(a)

t=2T

275 300 325

η

0

350 375

-0.05 400

0

Z

X

425

20 40

450

(b)

t=4T

275 300 325

η

0

350 375

-0.05 400

0

Z

X

425

20 40

450

(c)

t=6T

275 300 325

η

0

350 375

-0.05

400

0

Z

X

425

20 40

450

(d)

t=8T

Fig. 18. Membrane displacement g during subharmonic breakdown.

Fig. 15) as a result of their interaction with the large-amplitude oblique (1, 1) modes. The largeramplitude of the (1, 0) mode for the compliant membrane confirms an earlier observation about

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

400

X

500

(a)

30 20 10 0 -10 -20 -30

400

X

500

(b)

t=6T 30 20 10 0 -10 -20 -30

300

400

X

500

(c)

400

300

400

300

400

300

400

X

500

X

500

t=7T

X

500

Z

30 20 10

Z

30 20 10 0 -10 -20 -30

300

Z

300

Z

30 20 10 0 -10 -20 -30

t=5T

Z

300

Z

30 20 10 0 -10 -20 -30

Z

30 20 10 0 -10 -20 -30

Z

30 20 10 0 -10 -20 -30

1091

0 -10 -20 -30

300

400

X

500

(d)

t=8T

X

500

Fig. 19. Perturbation velocity u contours at different times during fundamental breakdown over a rigid wall (left) and the compliant membrane (right).

the distribution of the u-perturbation velocity contours in Fig. 19. The fundamental (1, 1) mode over the compliant membrane is also larger in amplitude than the same over the rigid wall. This may account for the stronger intensity of the first K-vortex on the membrane and the earlier breakdown of the K-vortex structure. The evolving displacement of the membrane surface is given in Fig. 21. The dimple-like depressions on the surface at t = 3T and 4T mark the locations just ahead of where a K-vortex is starting to form. A K-shaped impression may also be seen in the crest of the wave near x = 450 at t = 5T, which corresponds to the location of the first distinct K-vortex over the membrane (see u-contours in Fig. 19). The surface snapshots from t = 5T–7T depict the development of the K-vortex into a turbulent spot. Strong reflection of the primary 2D wave from the trailing edge of the membrane may be noted from t = 6T. Interaction between the oncoming and reflected waves leads to large displacement of surface (g < 0.07db.l.) that is outside the linear range. Such large surface displacement would most probably promote stronger non-linear interaction between the flow and the surface leading to even more rapid wave breakdown into turbulence. In summary, we note that finite compliant membranes with the appropriate properties are able to slow down the development of small-amplitude subharmonic and fundamental secondary instabilities; which is in agreement with the predictions of HerbertÕs secondary instability theory

1092

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095 10-1

10-2

mode (1,0)

|u|

mode (1,1)

10-3

10

-4

10-5

350

400

450

500

x

Fig. 20. Streamwise development of the Fourier amplitudes of selected wave modes at y = 0.26 during fundamental breakdown (membrane ——, rigid wall –  –  –).

[8] employed by Joslin et al. [7]. We also find that the finite compliant membrane investigated continues to be able to slow down the growth of oblique subharmonic (1/2, 1) waves of significantly large amplitude, but to be less successful with regards to large-amplitude fundamental (1, 1) mode.

7. Conclusions A numerical model based on a primitive-variables formulation of the perturbation Navier– Stokes equations is described. The numerical method and its treatment of flow-wall interaction may be regarded as an iterative variant of the implicit fractional step method. The governing equations are spatially discretized by a finite volume procedure on a transformed collocated (non-staggered) grid. A Rhie and Chow [20] type scheme is used to suppress spurious checkerboard fluctuations of the pressure field and solution. A multigrid procedure is used to solve the pressure-Poisson equation. The numerical model is implemented for parallel computation using domain decomposition and the Message Passing Interface (MPI) protocol. The generalized numerical model is developed with applications involving more complex geometries in mind, and its full potential is not exploited in the present study. The computational model is applied to study the spatial–temporal development of linear and non-linear waves in Blasius boundary layer over rigid and finite-length compliant membranes. Results thus obtained for oblique linear waves are in accord with results of linear stability theory. Our study verifies that finite-length compliant membrane is able to reduce the downstream growth of oblique linear TS waves. The TS wave is usually accompanied by co-existing CIFI mode(s)

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

1093

Fig. 21. Membrane displacement g during fundamental breakdown.

produced at the leading edge of the membrane. The CIFI mode propagates at different angle from the TS wave and its presence is most evident in the displacement of the membrane surface. The edges play a non-trivial role in determining the eventual wave system on a finite-length membrane, which may only be analyzed via a full spatial and temporal simulation, such as we have attempted here in 3D.

1094

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

The computational model is also applied to investigate the secondary instabilities of TS wave modes on a finite membrane. As linear flow-wall conditions are assumed in this model, a fairly stiff membrane is chosen for the non-linear study to keep the membrane displacement small. The finite membrane is found to reduce the spatial growth rates of the principal subharmonic (1/2, 1) and fundamental (1, 1) waves during the earlier stages of the instabilities, i.e. when the amplitudes of these waves are small. The membrane continues to exert a positive stabilizing effect on the subharmonic mode (1/2, 1) even when the amplitude of the oblique wave pair is fairly large. However, the membrane appears to be less successful with regards to a fundamental mode (1, 1) of the same amplitude. Here, slightly earlier breakdown of the K-vortex and its degeneration into a turbulent-like spot is observed compared with the same over a rigid surface. Large displacement of the membrane surface near the trailing edge of the membrane was also observed during the simulation of the fundamental breakdown, which suggests that the actual breakdown process could be more severe than represented by our use of linear flow-wall boundary conditions. Acknowledgments The financial support of a research scholarship from the National University of Singapore (NUS) and a top-up grant (C-014-000-109-531) from the Defense Science and Technology Agency (DSTA) is gratefully acknowledged. This work is also partially supported by Temasek Laboratories, National University of Singapore. References [1] Benjamin TB. The three-fold classification of unstable disturbances in flexible surfaces bounding inviscid flows. J Fluid Mech 1963;16:436. [2] Riley JJ, Gad-el-Hak M, Metcalfe RW. Compliant coatings. Ann Rev Fluid Mech 1988;20:393. [3] Carpenter PW, Garrad AD. The hydrodynamics stability of flow over Kramer-type compliant surfaces. Part 1. Tollmien–Schlichting instabilities. J Fluid Mech 1985;155:465. [4] Yeo KS. The stability of flow over flexible surfaces, PhD thesis, University of Cambridge, 1986. [5] Joslin RD, Morris PJ, Carpenter PW. The role of three-dimensional instabilities in compliant wall boundary-layer transition. AIAA J 1991;29:1603. [6] Yeo KS. The three-dimensional stability of boundary-layer flow over compliant walls. J Fluid Mech 1992;238:537. [7] Joslin RD, Morris PJ. The effect of compliant walls on secondary instabilities in boundary layer transition. AIAA J 1992;30(2):332. [8] Herbert T. Secondary instability of plane channel flow to subharmonic three-dimensional disturbances. Phys Fluids 1983;26(4):871. [9] Metcalfe RW, Battistoni F, Orzo S, Ekeroot J. Evolution of boundary layer flow over a compliant wall during transition to turbulence, in boundary layer transition and control, Cambridge UK. Proc Royal Aeronaut Soc 1991:36.1. [10] Davies C, Carpenter PW. Numerical simulation of the evolution of Tollmien–Schlichting waves over finite compliant panels. J Fluid Mech 1997;335:361. [11] Wiplier O, Ehrenstein U. Numerical simulation of linear and nonlinear disturbance evolution in a boundary layer with complaint walls. J Fluid Struct 2000;14:157. [12] Fasel HF, Rist U, Konzelmann U. Numerical investigation of the three-dimensional development in boundarylayer transition. AIAA J 1990;28(1):29.

Z. Wang et al. / Computers & Fluids 34 (2005) 1062–1095

1095

[13] Rist U, Fasel H. Direct numerical simulation of controlled transition in a flat-plate boundary layer. J Fluid Mech 1995;298:211. [14] Davies C, Carpenter PW. A Novel velocity–vorticity formulation of the Navier–Stokes equations with applications to boundary layer disturbance evolution. J Comput Phys 2001;172:119. [15] Liu C, Liu Z. Multigrid mapping box relaxation for simulation of the whole process of flow transition in 3D boundary layer. J Comput Phys 1995;119(2):325. [16] Lucey AD, Cafolla GJ, Carpenter PW, Yang M. The nonlinear hydroelastic behaviour of flexible walls. J Fluid Struct 1997;11:717. [17] Chorin AJ. On the convergence of discrete approximations to the Navier–Stokes equations. Math Comput 1969;23:509. [18] Kim J, Moin P. Application of a fractional-step method to incompressible Navier–Stokes equations. J Comput Phys 1985;59:308. [19] Melaaen MC. Calculation of fluid-flows with staggered and non-staggered curvilinear non-orthogonal grid—the theory. Numer Heat Transfer Part B—Fundamentals 1992;21(1):1. [20] Rhie CM, Chow WL. Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA J 1983;21(11):1525. [21] Patankar SV. Numerical heat transfer and fluid flow. Hemisphere Pub. Corp.; 1980. [22] Wesseling P, Oosterlee CW. Geometric multigrid with applications to computational fluid dynamics. J Comput Appl Math 2001;128:311. [23] Marc S, Steve O, Steven H, David W, Jack D. MPI: The complete reference. The MIT Press; 1996. [24] Saric WS, Kozlov VV, Levchenko CYa., 1984, AIAA Paper No. 84-0007 (unpublished). [25] Squire HB. On the stability of three-dimensional distribution of viscous fluid between parallel walls. Proc Roy Soc, Ser A 1933;142:621. [26] Wiplier O, Ehrenstein U. On the absolute instability in a boundary-layer flow with compliant coatings. Euro J Mech B—Fluids 2001;20:127. [27] Wang Z. Computational simulation of unsteady boundary layer over compliant surfaces, Ph.D. thesis, National University of Singapore, 2004.