Computers & Fluids 31 (2002) 683–693 www.elsevier.com/locate/compfluid
Numerical simulation of separating boundary-layer flow Matthieu Marquillie, Uwe Ehrenstein * Laboratoire J.-A. Dieudonn e, Universit e de Nice-Sophia Antipolis, Parc Valrose, F-06108 Nice Cedex 2, France Received 1 August 2001
Abstract A two-dimensional boundary-layer flow over a bump has been studied computationally by solving the incompressible Navier–Stokes equations. The numerical solution procedure uses a mapping transforming the physical coordinate system into a Cartesian one. The resulting system is solved using a mixed finite differences––Chebyshev collocation discretization. We compare the influence matrix technique with a fractional time-step procedure used to enforce the divergence-free condition. Numerical experiments are performed for bumps with different aspect ratios. The separation structure is investigated for increasing Reynolds number and self-induced vortex shedding is reproduced numerically. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Boundary-layer flow; Separating flow; Fractional-step method
1. Introduction The behavior of a boundary layer growing along a surface is highly influenced by the pressure distribution. For adverse pressure gradients the flow may separate. The flow over a backward facing step as a prototype of a flow geometry leading to separation has found much attention during the last decades from both the experimental and numerical point of view. Data for laminar separating flow may be found for instance in Ref. [1], among others. Numerical investigations of the stability of separating flow over a backward-facing step have been reported by Kaiktsis et al. [2,3], whereas for instance Le et al. [4] reported turbulent flow simulations for this prototype flow geometry. An alternative way to produce an adverse pressure gradient required for the formation of a separation bubble is the use of suction through the upper boundary. The influence of the resulting
*
Corresponding author. Tel.: +33-4-9207-6257; fax: +33-4-9351-7974. E-mail address: [email protected] (U. Ehrenstein).
0045-7930/02/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 3 0 ( 0 1 ) 0 0 0 7 4 - 3
684
M. Marquillie, U. Ehrenstein / Computers & Fluids 31 (2002) 683–693
Fig. 1. Geometry of the bump: (- - -) h ¼ 1, chord-to-height ratio Lc =h ¼ 15; (––) h ¼ 2, chord-to-height ratio Lc =h ¼ 10.
external adverse pressure gradient on separation has been numerically investigated by Pauley et al. [5,6]. In a very recent experimental investigation [7], separation of a laminar boundary layer on the plate is triggered imposing an adverse pressure gradient by an opposite contoured wall with suction. In the present numerical investigation a bump has been considered. The bump is formed by a convex surface with at the front and at the rear two additional concave regions in order to achieve continuity of curvature between the bump and the adjacent flat plate. The geometry (reminiscent of that of a wing) of the bump is sketched in Fig. 1 (for two different aspect ratios). This bump (with more or less similar aspect ratios) has been designed for measurements of fully developed turbulent flow [8] and it is suitable for investigation of both curvature and pressure-gradient effects. In the laminar regime the flow will first be accelerated and separation will then occur at the rear of the bump. In the present work a numerical method is presented capable of computing the separating boundary layer over the bump. Section 2 is devoted to the description of the numerical solution procedure. Section 3 contains some numerical experiments for two-dimensional flow, focusing first on comparison between the fractional-step (FS) method and the influence matrix (IM) technique. Then, steady-state solutions are investigated with respect to the recirculation length as function of the Reynolds number. Increasing the Reynolds number beyond a critical value, selfinduced vortex shedding is produced. Some discussion and a outlook is given in Section 4.
2. Problem definition and solution procedure We consider an incompressible fluid flow of viscosity m past a bump. The bump has been designed for the flow to transit continuously from the bump to the adjacent plane wall. The lower limit of the domain is hence given by the graph of the bump gðxÞ and the Navier–Stokes system has to be solved in the domain xa 6 x 6 xb , g ðx Þ 6 y < 1. The equations are made dimensionless using the displacement thickness sffiffiffiffiffiffiffiffi m xa ; c ¼ 17;208; ð1Þ da ¼ c U1 and hence at inflow xa as reference length, the reference velocity being the velocity at infinity U1 the Reynolds number is
M. Marquillie, U. Ehrenstein / Computers & Fluids 31 (2002) 683–693
Re ¼
U1 da : m
685
ð2Þ
2.1. Mapping and discretization Rather to write the Navier–Stokes system in curvilinear coordinates we merely transform the partial differential operators using the mapping (the barred quantities are the physical coordinates) t ¼ t;
x ¼ x;
y ¼ y gðxÞ:
ð3Þ
Hence the gradient, Laplacian and curl operators are transformed and ~¼r ~þG ~g ; r
¼ D þ Lg ; D
~ ðr ~ ~ ~ ðr ~ ~ r uÞ ¼ r uÞ þ ~ Rg ð~ uÞ
ð4Þ
with
2 2 og o o2 g o og o2 og o ; ; 0 ; Lg ¼ 2 2 þ ox oy ox oy ox oxoy ox oy 2 2 2 og o v og o u ~ Rg ð~ uÞ ¼ ; Lg v ; ox oy 2 ox oy 2 ~g ¼ G
ð5Þ
while the normal vector along the bump is 1 ~ n ¼ ð0; 1Þ and ~ ng ¼ n ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n þ~ ng Þ with ~ 2ffi ð~ 1 þ og ox
og ;0 : ox
ð6Þ
Hence, the physical domain is transformed into a Cartesian one in the computational variables ðx; yÞ. The two-dimensional Navier–Stokes system then writes o~ u 1 ~p G ~g p þ 1 D~ ~ Þ~ ~g Þ~ u ¼ r þ ð~ u r u þ ð~ u G u þ Lg~ u; ot Re Re
ð7Þ
~ ~ ~g ~ r u ¼ G u;
ð8Þ
which has to be solved in the domain xa 6 x 6 xb and 0 6 y < 1. The Poisson equation ou ov ou ov Dp ¼ Lg p þ J ð u; vÞ with J ðu; vÞ ¼ 2 ox oy oy ox
ð9Þ
for the pressure is obtained applying the divergence operator to the momentum equations together with Eq. (8). For space discretization fourth-order central finite differences are used for the second derivatives in the streamwise x-direction. All first derivatives of the flow quantities appear explicitly in the time-advancing scheme and the first derivatives in x are discretized using eighthorder finite differences. Chebyshev collocation is used in the wall-normal y-direction, the unbounded domain 0 6 y < 1 being mapped on the finite domain n 2 ½1; 1 using the algebraic transformation
686
M. Marquillie, U. Ehrenstein / Computers & Fluids 31 (2002) 683–693
y¼
ymax Lð1 þ nÞ : 2L þ ymax ð1 nÞ
ð10Þ
For time-integration implicit second-order backward Euler differencing is used; the Cartesian part of the Laplacian is taken implicitly whereas an explicit second-order Adams–Bashforth scheme is used for the operators depending on g as well as for the nonlinear convective terms. The resulting Poisson equations are solved efficiently using the matrix-diagonalization technique. (The second-derivative operator in y is diagonalized once for all and the system is hence transformed into a series of one-dimensional Helmholtz-like equations in x which can be solved efficiently using the pentadiagonal structure of the matrices.) A similar algorithm has previously been used for the computation of a boundary-layer flow interacting with a passive compliant coating [9]. In this latter work the IM-technique has been used in order to enforce the divergence-free condition using an appropriate Dirichlet-boundary condition for the pressure, solution of the Poisson Eq. (9). The FS-method, which has been meeting an increasing success during the last two decades and has given rise to more or less equivalent formulation for solving the incompressible Navier–Stokes equations (cf. Refs. [10–12]), has been adapted to the problem under investigation. 2.2. Fractional-step method The equations once discretized in time, the following system has to be solved at each time step ~ pnþ1 þ f n;n1 ðD 3sÞ~ unþ1 ¼ Rer Dpnþ1 ¼ ½Lg p þ J ðu; vÞ
n;n1
with s ¼
Re ; 2 Dt
;
ð11Þ ð12Þ
~ ~ r unþ1 ¼ 0;
ð13Þ
~ Þ~ ~g Þ~ f n;n1 ¼ 4s~ un þ s~ un1 ½Lg ~ u r u þ ð~ u G u n;n1 : u n;n1 þ Re½ð~
ð14Þ
with n;n1
Here the superscript n, n 1 means an explicit Adams–Bashforth differencing with ½
¼ 2½
n ½
n1 . At the lower boundary the non-slip condition is considered whereas at infinity uniform flow in streamwise direction is recovered setting u ¼ 1 and v ¼ 0. At inflow (upstream of the bump) the Blasius profile is imposed with u ¼ Ub and v ¼ 0. At outflow, a convective condition is used as advocated for instance by Pauley et al. [5] who used a similar condition for the computation of separating flow. The outflow condition is Z y o~ u o~ u 1 uðxb ; yÞ dy; ð15Þ þ Uc ¼ 0 with Uc ¼ ot ox y 0 with Uc a time-dependent mean convective velocity at outflow xb . For unsteady flow calculations characterized by vortex shedding, the propagation speed of the vortices may provide a suitable value for Uc [6]. Based on our numerical simulation results, the upper bound y of the integral has been chosen for the streamwise component of the velocity at outflow uðxb ; y Þ to be equal to 0:5, that is half the value of uniform flow at infinity. For unsteady flow the integration bound y
M. Marquillie, U. Ehrenstein / Computers & Fluids 31 (2002) 683–693
687
fluctuates in time (Uc is computed in the time-marching algorithm explicitly from values at the previous time step). It was found that the resulting convective velocity allowed the vortices to exit the domain without distortion. Although numerous papers have been devoted to the different formulations of the FS-method (or projection method) we briefly summarize the formulation used for the present computations. Indeed, in our formulation we had to take into account the presence of the metric terms depending on g, that is the graph of the bump. Knowing the flow quantities at previous time steps n, n 1, an intermediate pressure is computed via n;n1
Dp ¼ ½J ðu; vÞ Lg p
ð16Þ
;
with the Neumann-type boundary condition " 3~ unþ1 4~ un þ ~ un1 ~ ~ Þ~ ~g Þ~ ~g p rp ~ ð~ u r u þ ð~ u G n¼ uþ G 2 Dt n;n1 # 1 ~ ~ ~ ~ p n;n1 ~ uÞ þ ~ Rg ð~
ð~ n þ~ ng Þ ½r uÞÞ ng : þ ðr ðr Re
ð17Þ
Here we update at each time step the boundary condition for the pressure in order to minimize the effect of erroneous numerical boundary-layers induced by splitting methods, as suggested by Karniadakis et al. [12]. In a quite recent paper Hugues and Randriamampianina [13] showed that to update the pressure-gradient at the boundary provides indeed a time accuracy of the same order as that of the underlying temporal scheme. The intermediate pressure once obtained one recovers an intermediate velocity field solving ~ p þ f n;n1 ðD 3sÞ~ u ¼ Rer
with ~ uC ¼ ~ unþ1 C ;
ð18Þ
C being the boundary. Our numerical experiments showed that for seak of stability (in particular for steeper bumps) it is much preferable to solve during the projection step the pressure correction / ¼ pnþ1 p , with ~ / ¼ 3 ð~ ~ ~ r u Þ; r unþ1 ¼ 0; ð19Þ unþ1 ~ 2 Dt ~ containing both the Cartesian and the metric part. implicitly in the physical gradient operator r Applying the physical divergence operator to Eq. (19), using the incompressibility condition, and imposing homogeneous Neumann boundary conditions for the pressure correction, one recovers the following iteration sequence 3 ~ ~ ~ /kþ1 ~ ~g /kþ1 ~ ~ /k ~ r ~ u þ Gg ~ u Lg /k with r nþG n g ¼ r ng : ð20Þ D/kþ1 ¼ 2Dt In practice only one iteration proved to be necessary in order to recover a pressure correction up to second order in time which is the overall precision of the time-marching algorithm. Hence, the new pressure as well as a divergence-free velocity field is obtained writing ~ u unþ1 ¼ ~
2 Dt ~ ~g /Þ ðr/ þ G 3
and pnþ1 ¼ p þ /:
ð21Þ
688
M. Marquillie, U. Ehrenstein / Computers & Fluids 31 (2002) 683–693
3. Results and discussion First the bump with a dimensionless height of one (that is equal to the displacement thickness at inflow) has been considered. The corresponding graph is sketched as the broken line in Fig. 1, the chord-to-height ratio being approximately Lc =h ¼ 15. The inflow (where the Blasius profile is described) is at a distance of about x ¼ 6 upstream of the bump. 3.1. Fractional-step method versus influence matrix technique Two different approaches for computing the pressure ensuring a divergence-free velocity field have been compared: the IM-technique [14] and the FS-method described in Section 2. For time integration up to a steady-state solution a time step Dt ¼ 0:02 has been considered. A more or less long computational domain in the streamwise direction from inflow xa ¼ 0 to outflow xb has been considered. Almost all the computations have been performed using 500 discretization points in the streamwise direction whereas the wall-normal direction has been discretized using 76 collocation points. The stretching factor L in Eq. (10) has been set to L ¼ 6 when the IM-technique has been used, whereas a higher value L ¼ 14 proved to be appropriate for the FS-algorithm. An upper bound ymax ¼ 100 (respectively ymax ¼ 80) has been considered for the IM-approach (respectively the FS-approach). One comparison for an inflow Reynolds number Re ¼ 600 is depicted in Fig. 2, where the wallpressure distribution is shown. As expected for this steady-state solution, both the IM- and FSmethod exhibit the same pressure distribution along the bump. Increasing the computational domain up to xb ¼ 150 (for this latter case 750 points have been used in x) did not modify the pressure curve, as shown by the dotted line in Fig. 2. The pressure gradient changes from adverse (at the short concave part of the bump nearby inflow) to favorable at the convex surface and then
Fig. 2. Pressure distribution along the bump and at the wall, Re ¼ 600; IM: influence matrix technique; FS-method with outflow at xb ¼ 100: FS1, with xb ¼ 150: FS2.
M. Marquillie, U. Ehrenstein / Computers & Fluids 31 (2002) 683–693
689
becomes adverse, the minimum of the pressure being located at the summit of the bump. As expected, the pressure gradient tends to zero along the adjacent flat plate, for the larger computational domain. The flow is slightly accelerated at the outflow of the computational domain due to the convective boundary condition (15). However, this numerical condition hardly affects the pressure gradient upstream of the outflow boundary. In order to compare both the FS- and IM-approach for a time-dependent flow, the separated steady state at Re ¼ 600 has been perturbed using blowing suction at the summit of the bump. For this purpose the time-dependent boundary condition centered at x0 vðx0 þ qÞ ¼ AgðqÞ sinð2pftÞ
with gðqÞ ¼ aðq5 2q3 þ qÞ; 1 6 q 6 1;
ð22Þ
has been considered. The function g continuously fits the blowing-suction profile and the adjacent homogeneous boundary condition (with a a constant for the maximum amplitude to be A). For the computations, the frequency f ¼ 0:01 has been considered, with a blowing-suction amplitude A ¼ 0:1. Hence, a time-periodic fluid flow is triggered and vortex-shedding occurs, as shown in Fig. 3 where instantaneous streamlines for equally spaced time intervals over one period are depicted. The computational domain reaches from 0 6 x 6 60 and the time evolution for the pressure at the wall for three different x locations is shown in Fig. 4. The computations have been performed for both the FS- and IM-approach using an identical time-step Dt ¼ 0:01. The pressure difference pðx; 0Þ pð0; 0Þ is shown (for comparison, at each location the pressure curves start at zero for t ¼ 0 in Fig. 4), at x ¼ 15, that is slightly downstream the summit of the bump, at the center x ¼ 30 of the computational domain and at x ¼ 45. Using the IM-technique a Dirichlet-boundary condition for the pressure is computed, equivalent to the divergence-free condition. This method provides the pressure with the same time accuracy as the time-marching algorithm. As shown in Fig. 4 the FS-approach predicts an almost identical time-dependent wall-pressure distribution (despite some small differences in the overall amplitude, the frequency being identical). As discussed in Ref. [13], to update the Neumann
Fig. 3. Instantaneous streamlines over one period in time, for periodic vortex shedding triggered by blowing-suction at the summit of the bump, at Re ¼ 600.
690
M. Marquillie, U. Ehrenstein / Computers & Fluids 31 (2002) 683–693
Fig. 4. Time history for pressure distribution at the wall, flow triggered by blowing suction, at Re ¼ 600, at (from top to bottom) x ¼ 15, 30, 45: (––) FS-approach; (- - -) IM-technique.
condition for the pressure according to Eq. (17) at each time step improves the FS-approach with regard to time accuracy of the pressure. We emphasize that the comparison is made for a physically meaningful flow field, which in particular is affected by the numerical boundary conditions imposed at inflow and at outflow. The FS-method applied to the present flow geometry hence provides reliable results and in the following this approach is used in order to numerically describe the separated flow at the rear of the bump. 3.2. Steady separation As shown by Kaiktsis et al. [3], the asymptotic flow states for a two-dimensional flow over a backward-facing step are stable for Reynolds numbers (based on the height of the step) up to 2500, in absence of any external excitation. Here we considered first the bump with the dimensionless height of 1 and one may expect to recover time-independent states with separation for increasing Reynolds numbers, by analogy with the results of Kaiktsis et al. Asymptotic steady states could indeed be recovered (for these computations the FS-method has been used), for increasing Reynolds number from Re ¼ 300 to Re ¼ 1300. Considering a second bump depicted as the solid line in Fig. 1 twice as high and with a smaller chord-to-height ratio Lc =h ¼ 10, similar steady-state solutions with separation could be found and corresponding streamlines are shown in Fig. 5. The Reynolds numbers are Re ¼ 300, 450, formed with the displacement thickness at inflow. Using the bump height h as reference length the corresponding Reynolds numbers Reh are Reh ¼ 600, 900.
M. Marquillie, U. Ehrenstein / Computers & Fluids 31 (2002) 683–693
691
Fig. 5. Streamlines for steady separation, bump with h ¼ 2 and chord-to-height ratio Lc =h ¼ 10 (outflow at xb ¼ 100), for increasing Reynolds number Re ¼ 300, 450.
Fig. 6. Normalized recirculation length xr =h as function of the Reynolds number Reh based on bump height h: (þ) bump with h ¼ 1 (Lc =h ¼ 15) and ( ) bump with h ¼ 2 (Lc =h ¼ 10).
One hence may wonder whether the recirculation length is independent of the specific bump shape considered. For this purpose the recirculation length has been normalized using the bump height as reference length and the results are shown in Fig. 6 as function of Reh . The recirculation length for the bump with h ¼ 1 depends almost linearly on Reh . The two results for the higher bump are slightly above those for the smaller bump with the same Reh values, which means that the recirculation length also depends on the geometry. 3.3. Unsteady separation bubble For the different bump geometries considered one may expect a critical Reynolds number for the separation bubble to become unstable, the instability being characterized by vortex shedding. Focusing on the bump with height h ¼ 2 (and Lc =h ¼ 10), we considered a sufficiently long computational domain for the flow field inside the separation bubble to be independent of the outflow conditions. Considering a domain 0 6 x 6 190 (for these computations 950 points have been used in x and 92 collocation points in y) we increased the inflow Reynolds number and for
692
M. Marquillie, U. Ehrenstein / Computers & Fluids 31 (2002) 683–693
Fig. 7. Instantaneous streamlines in the domain 0 6 x 6 120, bump with height h ¼ 2 (Lc =h ¼ 10), at Re ¼ 600.
Fig. 8. Streamwise velocity at ðx; yÞ ¼ ð120; 0:4Þ as function of time, bump with height h ¼ 2 (Lc =h ¼ 10), at Re ¼ 600.
Re P 500 the steady state, for a long time integration, progressively became unstable. Finally, for Re ¼ 600 no transiently stable flow state could be recovered. Some instantaneous streamlines for this latter case are shown in Fig. 7. While we did not attempt to accurately locate the critical Reynolds number it should be in the vicinity of Re ¼ 550. Finally, the streamwise velocity at ðx; yÞ ¼ ð120; 0:4Þ, as function of time, is depicted in Fig. 8, illustrating the unsteadiness of the flow field.
4. Concluding remarks A numerical solution procedure suitable for the computation of a separating boundary-layer flow has been presented. Here the adverse pressure-gradient leading to the detachment of the flow is induced at the rear of a bump. Steady states with separation could be recovered for increasing Reynolds number. While there is evidence of vortex shedding above critical Reynolds numbers, the stability of the flow with recirculation has to be assessed by forthcoming numerical experiments. The FS-method used proved to be robust and capable of reproducing the pressure distribution along the bump, which is responsible for the flow structure with separation. The FS-approach necessitates less computer storage compared to the IM-technique, which is particularly important when extending the solution procedure to three-dimensional problems; it proved also to be less CPU consuming (with a gain of about 20%). The numerical technique used may be easily generalized by taking into account three-dimensionality. A three-dimensional version of the algorithm, the geometry being homogeneous in the third (spanwise) coordinate, is under development. The ultimate goal being to contribute to the understanding of the instabilities and transition in separating boundary-layer type flows.
M. Marquillie, U. Ehrenstein / Computers & Fluids 31 (2002) 683–693
693
Acknowledgements Parts of the computations have been performed on NEC-SX5 of the IDRIS, France.
References [1] Sinha SN, Gupta AK, Oberai MM. Laminar separating flow over backsteps and cavities. Part I: backsteps. AIAA J 1981;19:1527–30. [2] Kaiktsis L, Karniadakis GE, Orszag SA. Onset of three-dimensionality, equilibria, and early transition in flow over a backward-facing step. J Fluid Mech 1991;231:501–28. [3] Kaiktsis L, Karniadakis GE, Orszag SA. Unsteadiness and convective instabilities in two-dimensional flow over a backward-facing step. J Fluid Mech 1996;321:157–87. [4] Le H, Moin P, Kim J. Direct numerical simulation of turbulent flow over a backward-facing step. J Fluid Mech 1997;330:349–74. [5] Pauley LL, Moin P, Reynolds WC. The structure of two-dimensional separation. J Fluid Mech 1990;220:397–411. [6] Pauley LL. Structure of local pressure-driven three-dimensional transient boundary-layer separation. AIAA J 1994;32:997–1005. [7] H€aggmark CP, Bakchinov AA, Alfredsson PH. Experiments on a two-dimensional laminar separation bubble. Philos Trans R Soc London A 2000;358:3198–205. [8] Bernard A, Foucaut JM, Dupont P, Stanislas M. Control of a decelerating boundary-layer. Part 1: characterization of the flow. AIAA J, in press. [9] Wiplier O, Ehrenstein U. Numerical simulation of linear and nonlinear disturbance evolution in a boundary layer with compliant walls. J Fluids Struct 2000;14:157–82. [10] Kim J, Moin P. Application of a fractional-step method to incompressible Navier–Stokes equations. J Comput Phys 1985;59:308–23. [11] Streett CL, Macaraeg MG. Spectral multi-domain for large-scale fluid dynamic simulations. Appl Numer Math 1989/1990;6:123–39. [12] Karniadakis GE, Israeli M, Orszag SA. High-order splitting methods for the incompressible Navier–Stokes equations. J Comput Phys 1991;97:414–43. [13] Hugues S, Randriamampianina A. An improved projection scheme applied to pseudospectral methods for the incompressible Navier–Stokes equations. Int J Numer Methods Fluids 1998;28:501–21. [14] Ehrenstein U, Peyret R. A Chebyshev collocation method for the Navier–Stokes equations with applications to double-diffusive convection. Int J Numer Methods Fluids 1989;9:427–52.
Numerical simulation of separating boundary-layer flow Matthieu Marquillie, Uwe Ehrenstein * Laboratoire J.-A. Dieudonn e, Universit e de Nice-Sophia Antipolis, Parc Valrose, F-06108 Nice Cedex 2, France Received 1 August 2001
Abstract A two-dimensional boundary-layer flow over a bump has been studied computationally by solving the incompressible Navier–Stokes equations. The numerical solution procedure uses a mapping transforming the physical coordinate system into a Cartesian one. The resulting system is solved using a mixed finite differences––Chebyshev collocation discretization. We compare the influence matrix technique with a fractional time-step procedure used to enforce the divergence-free condition. Numerical experiments are performed for bumps with different aspect ratios. The separation structure is investigated for increasing Reynolds number and self-induced vortex shedding is reproduced numerically. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Boundary-layer flow; Separating flow; Fractional-step method
1. Introduction The behavior of a boundary layer growing along a surface is highly influenced by the pressure distribution. For adverse pressure gradients the flow may separate. The flow over a backward facing step as a prototype of a flow geometry leading to separation has found much attention during the last decades from both the experimental and numerical point of view. Data for laminar separating flow may be found for instance in Ref. [1], among others. Numerical investigations of the stability of separating flow over a backward-facing step have been reported by Kaiktsis et al. [2,3], whereas for instance Le et al. [4] reported turbulent flow simulations for this prototype flow geometry. An alternative way to produce an adverse pressure gradient required for the formation of a separation bubble is the use of suction through the upper boundary. The influence of the resulting
*
Corresponding author. Tel.: +33-4-9207-6257; fax: +33-4-9351-7974. E-mail address: [email protected] (U. Ehrenstein).
0045-7930/02/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 3 0 ( 0 1 ) 0 0 0 7 4 - 3
684
M. Marquillie, U. Ehrenstein / Computers & Fluids 31 (2002) 683–693
Fig. 1. Geometry of the bump: (- - -) h ¼ 1, chord-to-height ratio Lc =h ¼ 15; (––) h ¼ 2, chord-to-height ratio Lc =h ¼ 10.
external adverse pressure gradient on separation has been numerically investigated by Pauley et al. [5,6]. In a very recent experimental investigation [7], separation of a laminar boundary layer on the plate is triggered imposing an adverse pressure gradient by an opposite contoured wall with suction. In the present numerical investigation a bump has been considered. The bump is formed by a convex surface with at the front and at the rear two additional concave regions in order to achieve continuity of curvature between the bump and the adjacent flat plate. The geometry (reminiscent of that of a wing) of the bump is sketched in Fig. 1 (for two different aspect ratios). This bump (with more or less similar aspect ratios) has been designed for measurements of fully developed turbulent flow [8] and it is suitable for investigation of both curvature and pressure-gradient effects. In the laminar regime the flow will first be accelerated and separation will then occur at the rear of the bump. In the present work a numerical method is presented capable of computing the separating boundary layer over the bump. Section 2 is devoted to the description of the numerical solution procedure. Section 3 contains some numerical experiments for two-dimensional flow, focusing first on comparison between the fractional-step (FS) method and the influence matrix (IM) technique. Then, steady-state solutions are investigated with respect to the recirculation length as function of the Reynolds number. Increasing the Reynolds number beyond a critical value, selfinduced vortex shedding is produced. Some discussion and a outlook is given in Section 4.
2. Problem definition and solution procedure We consider an incompressible fluid flow of viscosity m past a bump. The bump has been designed for the flow to transit continuously from the bump to the adjacent plane wall. The lower limit of the domain is hence given by the graph of the bump gðxÞ and the Navier–Stokes system has to be solved in the domain xa 6 x 6 xb , g ðx Þ 6 y < 1. The equations are made dimensionless using the displacement thickness sffiffiffiffiffiffiffiffi m xa ; c ¼ 17;208; ð1Þ da ¼ c U1 and hence at inflow xa as reference length, the reference velocity being the velocity at infinity U1 the Reynolds number is
M. Marquillie, U. Ehrenstein / Computers & Fluids 31 (2002) 683–693
Re ¼
U1 da : m
685
ð2Þ
2.1. Mapping and discretization Rather to write the Navier–Stokes system in curvilinear coordinates we merely transform the partial differential operators using the mapping (the barred quantities are the physical coordinates) t ¼ t;
x ¼ x;
y ¼ y gðxÞ:
ð3Þ
Hence the gradient, Laplacian and curl operators are transformed and ~¼r ~þG ~g ; r
¼ D þ Lg ; D
~ ðr ~ ~ ~ ðr ~ ~ r uÞ ¼ r uÞ þ ~ Rg ð~ uÞ
ð4Þ
with
2 2 og o o2 g o og o2 og o ; ; 0 ; Lg ¼ 2 2 þ ox oy ox oy ox oxoy ox oy 2 2 2 og o v og o u ~ Rg ð~ uÞ ¼ ; Lg v ; ox oy 2 ox oy 2 ~g ¼ G
ð5Þ
while the normal vector along the bump is 1 ~ n ¼ ð0; 1Þ and ~ ng ¼ n ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n þ~ ng Þ with ~ 2ffi ð~ 1 þ og ox
og ;0 : ox
ð6Þ
Hence, the physical domain is transformed into a Cartesian one in the computational variables ðx; yÞ. The two-dimensional Navier–Stokes system then writes o~ u 1 ~p G ~g p þ 1 D~ ~ Þ~ ~g Þ~ u ¼ r þ ð~ u r u þ ð~ u G u þ Lg~ u; ot Re Re
ð7Þ
~ ~ ~g ~ r u ¼ G u;
ð8Þ
which has to be solved in the domain xa 6 x 6 xb and 0 6 y < 1. The Poisson equation ou ov ou ov Dp ¼ Lg p þ J ð u; vÞ with J ðu; vÞ ¼ 2 ox oy oy ox
ð9Þ
for the pressure is obtained applying the divergence operator to the momentum equations together with Eq. (8). For space discretization fourth-order central finite differences are used for the second derivatives in the streamwise x-direction. All first derivatives of the flow quantities appear explicitly in the time-advancing scheme and the first derivatives in x are discretized using eighthorder finite differences. Chebyshev collocation is used in the wall-normal y-direction, the unbounded domain 0 6 y < 1 being mapped on the finite domain n 2 ½1; 1 using the algebraic transformation
686
M. Marquillie, U. Ehrenstein / Computers & Fluids 31 (2002) 683–693
y¼
ymax Lð1 þ nÞ : 2L þ ymax ð1 nÞ
ð10Þ
For time-integration implicit second-order backward Euler differencing is used; the Cartesian part of the Laplacian is taken implicitly whereas an explicit second-order Adams–Bashforth scheme is used for the operators depending on g as well as for the nonlinear convective terms. The resulting Poisson equations are solved efficiently using the matrix-diagonalization technique. (The second-derivative operator in y is diagonalized once for all and the system is hence transformed into a series of one-dimensional Helmholtz-like equations in x which can be solved efficiently using the pentadiagonal structure of the matrices.) A similar algorithm has previously been used for the computation of a boundary-layer flow interacting with a passive compliant coating [9]. In this latter work the IM-technique has been used in order to enforce the divergence-free condition using an appropriate Dirichlet-boundary condition for the pressure, solution of the Poisson Eq. (9). The FS-method, which has been meeting an increasing success during the last two decades and has given rise to more or less equivalent formulation for solving the incompressible Navier–Stokes equations (cf. Refs. [10–12]), has been adapted to the problem under investigation. 2.2. Fractional-step method The equations once discretized in time, the following system has to be solved at each time step ~ pnþ1 þ f n;n1 ðD 3sÞ~ unþ1 ¼ Rer Dpnþ1 ¼ ½Lg p þ J ðu; vÞ
n;n1
with s ¼
Re ; 2 Dt
;
ð11Þ ð12Þ
~ ~ r unþ1 ¼ 0;
ð13Þ
~ Þ~ ~g Þ~ f n;n1 ¼ 4s~ un þ s~ un1 ½Lg ~ u r u þ ð~ u G u n;n1 : u n;n1 þ Re½ð~
ð14Þ
with n;n1
Here the superscript n, n 1 means an explicit Adams–Bashforth differencing with ½
¼ 2½
n ½
n1 . At the lower boundary the non-slip condition is considered whereas at infinity uniform flow in streamwise direction is recovered setting u ¼ 1 and v ¼ 0. At inflow (upstream of the bump) the Blasius profile is imposed with u ¼ Ub and v ¼ 0. At outflow, a convective condition is used as advocated for instance by Pauley et al. [5] who used a similar condition for the computation of separating flow. The outflow condition is Z y o~ u o~ u 1 uðxb ; yÞ dy; ð15Þ þ Uc ¼ 0 with Uc ¼ ot ox y 0 with Uc a time-dependent mean convective velocity at outflow xb . For unsteady flow calculations characterized by vortex shedding, the propagation speed of the vortices may provide a suitable value for Uc [6]. Based on our numerical simulation results, the upper bound y of the integral has been chosen for the streamwise component of the velocity at outflow uðxb ; y Þ to be equal to 0:5, that is half the value of uniform flow at infinity. For unsteady flow the integration bound y
M. Marquillie, U. Ehrenstein / Computers & Fluids 31 (2002) 683–693
687
fluctuates in time (Uc is computed in the time-marching algorithm explicitly from values at the previous time step). It was found that the resulting convective velocity allowed the vortices to exit the domain without distortion. Although numerous papers have been devoted to the different formulations of the FS-method (or projection method) we briefly summarize the formulation used for the present computations. Indeed, in our formulation we had to take into account the presence of the metric terms depending on g, that is the graph of the bump. Knowing the flow quantities at previous time steps n, n 1, an intermediate pressure is computed via n;n1
Dp ¼ ½J ðu; vÞ Lg p
ð16Þ
;
with the Neumann-type boundary condition " 3~ unþ1 4~ un þ ~ un1 ~ ~ Þ~ ~g Þ~ ~g p rp ~ ð~ u r u þ ð~ u G n¼ uþ G 2 Dt n;n1 # 1 ~ ~ ~ ~ p n;n1 ~ uÞ þ ~ Rg ð~
ð~ n þ~ ng Þ ½r uÞÞ ng : þ ðr ðr Re
ð17Þ
Here we update at each time step the boundary condition for the pressure in order to minimize the effect of erroneous numerical boundary-layers induced by splitting methods, as suggested by Karniadakis et al. [12]. In a quite recent paper Hugues and Randriamampianina [13] showed that to update the pressure-gradient at the boundary provides indeed a time accuracy of the same order as that of the underlying temporal scheme. The intermediate pressure once obtained one recovers an intermediate velocity field solving ~ p þ f n;n1 ðD 3sÞ~ u ¼ Rer
with ~ uC ¼ ~ unþ1 C ;
ð18Þ
C being the boundary. Our numerical experiments showed that for seak of stability (in particular for steeper bumps) it is much preferable to solve during the projection step the pressure correction / ¼ pnþ1 p , with ~ / ¼ 3 ð~ ~ ~ r u Þ; r unþ1 ¼ 0; ð19Þ unþ1 ~ 2 Dt ~ containing both the Cartesian and the metric part. implicitly in the physical gradient operator r Applying the physical divergence operator to Eq. (19), using the incompressibility condition, and imposing homogeneous Neumann boundary conditions for the pressure correction, one recovers the following iteration sequence 3 ~ ~ ~ /kþ1 ~ ~g /kþ1 ~ ~ /k ~ r ~ u þ Gg ~ u Lg /k with r nþG n g ¼ r ng : ð20Þ D/kþ1 ¼ 2Dt In practice only one iteration proved to be necessary in order to recover a pressure correction up to second order in time which is the overall precision of the time-marching algorithm. Hence, the new pressure as well as a divergence-free velocity field is obtained writing ~ u unþ1 ¼ ~
2 Dt ~ ~g /Þ ðr/ þ G 3
and pnþ1 ¼ p þ /:
ð21Þ
688
M. Marquillie, U. Ehrenstein / Computers & Fluids 31 (2002) 683–693
3. Results and discussion First the bump with a dimensionless height of one (that is equal to the displacement thickness at inflow) has been considered. The corresponding graph is sketched as the broken line in Fig. 1, the chord-to-height ratio being approximately Lc =h ¼ 15. The inflow (where the Blasius profile is described) is at a distance of about x ¼ 6 upstream of the bump. 3.1. Fractional-step method versus influence matrix technique Two different approaches for computing the pressure ensuring a divergence-free velocity field have been compared: the IM-technique [14] and the FS-method described in Section 2. For time integration up to a steady-state solution a time step Dt ¼ 0:02 has been considered. A more or less long computational domain in the streamwise direction from inflow xa ¼ 0 to outflow xb has been considered. Almost all the computations have been performed using 500 discretization points in the streamwise direction whereas the wall-normal direction has been discretized using 76 collocation points. The stretching factor L in Eq. (10) has been set to L ¼ 6 when the IM-technique has been used, whereas a higher value L ¼ 14 proved to be appropriate for the FS-algorithm. An upper bound ymax ¼ 100 (respectively ymax ¼ 80) has been considered for the IM-approach (respectively the FS-approach). One comparison for an inflow Reynolds number Re ¼ 600 is depicted in Fig. 2, where the wallpressure distribution is shown. As expected for this steady-state solution, both the IM- and FSmethod exhibit the same pressure distribution along the bump. Increasing the computational domain up to xb ¼ 150 (for this latter case 750 points have been used in x) did not modify the pressure curve, as shown by the dotted line in Fig. 2. The pressure gradient changes from adverse (at the short concave part of the bump nearby inflow) to favorable at the convex surface and then
Fig. 2. Pressure distribution along the bump and at the wall, Re ¼ 600; IM: influence matrix technique; FS-method with outflow at xb ¼ 100: FS1, with xb ¼ 150: FS2.
M. Marquillie, U. Ehrenstein / Computers & Fluids 31 (2002) 683–693
689
becomes adverse, the minimum of the pressure being located at the summit of the bump. As expected, the pressure gradient tends to zero along the adjacent flat plate, for the larger computational domain. The flow is slightly accelerated at the outflow of the computational domain due to the convective boundary condition (15). However, this numerical condition hardly affects the pressure gradient upstream of the outflow boundary. In order to compare both the FS- and IM-approach for a time-dependent flow, the separated steady state at Re ¼ 600 has been perturbed using blowing suction at the summit of the bump. For this purpose the time-dependent boundary condition centered at x0 vðx0 þ qÞ ¼ AgðqÞ sinð2pftÞ
with gðqÞ ¼ aðq5 2q3 þ qÞ; 1 6 q 6 1;
ð22Þ
has been considered. The function g continuously fits the blowing-suction profile and the adjacent homogeneous boundary condition (with a a constant for the maximum amplitude to be A). For the computations, the frequency f ¼ 0:01 has been considered, with a blowing-suction amplitude A ¼ 0:1. Hence, a time-periodic fluid flow is triggered and vortex-shedding occurs, as shown in Fig. 3 where instantaneous streamlines for equally spaced time intervals over one period are depicted. The computational domain reaches from 0 6 x 6 60 and the time evolution for the pressure at the wall for three different x locations is shown in Fig. 4. The computations have been performed for both the FS- and IM-approach using an identical time-step Dt ¼ 0:01. The pressure difference pðx; 0Þ pð0; 0Þ is shown (for comparison, at each location the pressure curves start at zero for t ¼ 0 in Fig. 4), at x ¼ 15, that is slightly downstream the summit of the bump, at the center x ¼ 30 of the computational domain and at x ¼ 45. Using the IM-technique a Dirichlet-boundary condition for the pressure is computed, equivalent to the divergence-free condition. This method provides the pressure with the same time accuracy as the time-marching algorithm. As shown in Fig. 4 the FS-approach predicts an almost identical time-dependent wall-pressure distribution (despite some small differences in the overall amplitude, the frequency being identical). As discussed in Ref. [13], to update the Neumann
Fig. 3. Instantaneous streamlines over one period in time, for periodic vortex shedding triggered by blowing-suction at the summit of the bump, at Re ¼ 600.
690
M. Marquillie, U. Ehrenstein / Computers & Fluids 31 (2002) 683–693
Fig. 4. Time history for pressure distribution at the wall, flow triggered by blowing suction, at Re ¼ 600, at (from top to bottom) x ¼ 15, 30, 45: (––) FS-approach; (- - -) IM-technique.
condition for the pressure according to Eq. (17) at each time step improves the FS-approach with regard to time accuracy of the pressure. We emphasize that the comparison is made for a physically meaningful flow field, which in particular is affected by the numerical boundary conditions imposed at inflow and at outflow. The FS-method applied to the present flow geometry hence provides reliable results and in the following this approach is used in order to numerically describe the separated flow at the rear of the bump. 3.2. Steady separation As shown by Kaiktsis et al. [3], the asymptotic flow states for a two-dimensional flow over a backward-facing step are stable for Reynolds numbers (based on the height of the step) up to 2500, in absence of any external excitation. Here we considered first the bump with the dimensionless height of 1 and one may expect to recover time-independent states with separation for increasing Reynolds numbers, by analogy with the results of Kaiktsis et al. Asymptotic steady states could indeed be recovered (for these computations the FS-method has been used), for increasing Reynolds number from Re ¼ 300 to Re ¼ 1300. Considering a second bump depicted as the solid line in Fig. 1 twice as high and with a smaller chord-to-height ratio Lc =h ¼ 10, similar steady-state solutions with separation could be found and corresponding streamlines are shown in Fig. 5. The Reynolds numbers are Re ¼ 300, 450, formed with the displacement thickness at inflow. Using the bump height h as reference length the corresponding Reynolds numbers Reh are Reh ¼ 600, 900.
M. Marquillie, U. Ehrenstein / Computers & Fluids 31 (2002) 683–693
691
Fig. 5. Streamlines for steady separation, bump with h ¼ 2 and chord-to-height ratio Lc =h ¼ 10 (outflow at xb ¼ 100), for increasing Reynolds number Re ¼ 300, 450.
Fig. 6. Normalized recirculation length xr =h as function of the Reynolds number Reh based on bump height h: (þ) bump with h ¼ 1 (Lc =h ¼ 15) and ( ) bump with h ¼ 2 (Lc =h ¼ 10).
One hence may wonder whether the recirculation length is independent of the specific bump shape considered. For this purpose the recirculation length has been normalized using the bump height as reference length and the results are shown in Fig. 6 as function of Reh . The recirculation length for the bump with h ¼ 1 depends almost linearly on Reh . The two results for the higher bump are slightly above those for the smaller bump with the same Reh values, which means that the recirculation length also depends on the geometry. 3.3. Unsteady separation bubble For the different bump geometries considered one may expect a critical Reynolds number for the separation bubble to become unstable, the instability being characterized by vortex shedding. Focusing on the bump with height h ¼ 2 (and Lc =h ¼ 10), we considered a sufficiently long computational domain for the flow field inside the separation bubble to be independent of the outflow conditions. Considering a domain 0 6 x 6 190 (for these computations 950 points have been used in x and 92 collocation points in y) we increased the inflow Reynolds number and for
692
M. Marquillie, U. Ehrenstein / Computers & Fluids 31 (2002) 683–693
Fig. 7. Instantaneous streamlines in the domain 0 6 x 6 120, bump with height h ¼ 2 (Lc =h ¼ 10), at Re ¼ 600.
Fig. 8. Streamwise velocity at ðx; yÞ ¼ ð120; 0:4Þ as function of time, bump with height h ¼ 2 (Lc =h ¼ 10), at Re ¼ 600.
Re P 500 the steady state, for a long time integration, progressively became unstable. Finally, for Re ¼ 600 no transiently stable flow state could be recovered. Some instantaneous streamlines for this latter case are shown in Fig. 7. While we did not attempt to accurately locate the critical Reynolds number it should be in the vicinity of Re ¼ 550. Finally, the streamwise velocity at ðx; yÞ ¼ ð120; 0:4Þ, as function of time, is depicted in Fig. 8, illustrating the unsteadiness of the flow field.
4. Concluding remarks A numerical solution procedure suitable for the computation of a separating boundary-layer flow has been presented. Here the adverse pressure-gradient leading to the detachment of the flow is induced at the rear of a bump. Steady states with separation could be recovered for increasing Reynolds number. While there is evidence of vortex shedding above critical Reynolds numbers, the stability of the flow with recirculation has to be assessed by forthcoming numerical experiments. The FS-method used proved to be robust and capable of reproducing the pressure distribution along the bump, which is responsible for the flow structure with separation. The FS-approach necessitates less computer storage compared to the IM-technique, which is particularly important when extending the solution procedure to three-dimensional problems; it proved also to be less CPU consuming (with a gain of about 20%). The numerical technique used may be easily generalized by taking into account three-dimensionality. A three-dimensional version of the algorithm, the geometry being homogeneous in the third (spanwise) coordinate, is under development. The ultimate goal being to contribute to the understanding of the instabilities and transition in separating boundary-layer type flows.
M. Marquillie, U. Ehrenstein / Computers & Fluids 31 (2002) 683–693
693
Acknowledgements Parts of the computations have been performed on NEC-SX5 of the IDRIS, France.
References [1] Sinha SN, Gupta AK, Oberai MM. Laminar separating flow over backsteps and cavities. Part I: backsteps. AIAA J 1981;19:1527–30. [2] Kaiktsis L, Karniadakis GE, Orszag SA. Onset of three-dimensionality, equilibria, and early transition in flow over a backward-facing step. J Fluid Mech 1991;231:501–28. [3] Kaiktsis L, Karniadakis GE, Orszag SA. Unsteadiness and convective instabilities in two-dimensional flow over a backward-facing step. J Fluid Mech 1996;321:157–87. [4] Le H, Moin P, Kim J. Direct numerical simulation of turbulent flow over a backward-facing step. J Fluid Mech 1997;330:349–74. [5] Pauley LL, Moin P, Reynolds WC. The structure of two-dimensional separation. J Fluid Mech 1990;220:397–411. [6] Pauley LL. Structure of local pressure-driven three-dimensional transient boundary-layer separation. AIAA J 1994;32:997–1005. [7] H€aggmark CP, Bakchinov AA, Alfredsson PH. Experiments on a two-dimensional laminar separation bubble. Philos Trans R Soc London A 2000;358:3198–205. [8] Bernard A, Foucaut JM, Dupont P, Stanislas M. Control of a decelerating boundary-layer. Part 1: characterization of the flow. AIAA J, in press. [9] Wiplier O, Ehrenstein U. Numerical simulation of linear and nonlinear disturbance evolution in a boundary layer with compliant walls. J Fluids Struct 2000;14:157–82. [10] Kim J, Moin P. Application of a fractional-step method to incompressible Navier–Stokes equations. J Comput Phys 1985;59:308–23. [11] Streett CL, Macaraeg MG. Spectral multi-domain for large-scale fluid dynamic simulations. Appl Numer Math 1989/1990;6:123–39. [12] Karniadakis GE, Israeli M, Orszag SA. High-order splitting methods for the incompressible Navier–Stokes equations. J Comput Phys 1991;97:414–43. [13] Hugues S, Randriamampianina A. An improved projection scheme applied to pseudospectral methods for the incompressible Navier–Stokes equations. Int J Numer Methods Fluids 1998;28:501–21. [14] Ehrenstein U, Peyret R. A Chebyshev collocation method for the Navier–Stokes equations with applications to double-diffusive convection. Int J Numer Methods Fluids 1989;9:427–52.
