- Email: [email protected]
Large eddy simulation of turbulent flows in external flow field using three-step FEM–BEM model
Engineering Analysis with Boundary Elements 30 (2006) 564–576 www.elsevier.com/locate/enganabound
Large eddy simulation of turbulent flows in external flow field using three-step FEM–BEM model D.L. Young a,*, T.I. Eldho b, J.T. Chang a a
Department of Civil Engineering and Hydrotech Research Institute, National Taiwan University, 1, Sec. 4, Roosvelt Road, Taipei 10617, Taiwan b Department of Civil Engineering, Indian Institute of Technology, Mumbai 400 076, India Received 28 March 2005; accepted 2 February 2006 Available online 19 April 2006
Abstract An innovative computational model is presented for the large eddy simulation (LES) modeling of multi-dimensional unsteady turbulent flow problems in external flow field. Based on the LES principles, the model uses a pressure projection method to solve the Navier–Stokes equations in transient condition. The turbulent motion is simulated by Smagorinsky sub-grid scale (SGS) eddy viscosity model. The momentum equation of the flow motion is solved using a three-step finite element method (FEM). The external flow field is simulated using a boundary element method (BEM) by solving a pressure Poisson equation that assumes the pressure as zero at the infinity. Through extracting the boundary effects on a specified finite computational domain, the model is able to solve the infinite boundary value problems. The present model is used to simulate the flows past a two-dimensional square rib and a three-dimensional cube at high Reynolds number. The simulation results are found to be reasonable and comparable with other models available in the literature even for coarse meshes. q 2006 Elsevier Ltd. All rights reserved. Keywords: Large eddy simulation; Turbulent flows; Multi-dimensional external flows; Navier–Stokes equations; Three-step finite element method; Boundary element method
1. Introduction Three-dimensional turbulent, external flows and their interaction with structures play a vital role in the design of high-rise buildings, suspended bridges and aircrafts. Many researchers have investigated inherently three-dimensional separation zones around structures in a surface boundary layer. Experimental investigations carried out to simulate these external flow phenomena, provided limited information for the entire flow field. Hence numerical models are essential for predicting the turbulent flow fields around buildings and other complex flows. One difficulty in simulating these flows lies in effectively setting up the boundary conditions and also the computational domain to solve external flows with infinite domain. In these types of problems, the vortex shedding around the structure plays an important role in characterizing various types of fluid-induced structural vibrations [1]. To develop a computational model for this class of turbulent flow problems, it becomes necessary to deal with the unsteady Navier–Stokes * Corresponding author. Tel./fax: C886 2 2362 6114. E-mail address: [email protected] (D.L. Young).
0955-7997/$ - see front matter q 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2006.02.004
equations with a turbulence closure form such as the LES model in a multi-dimensional infinite domain. To study the turbulent flows in external flow field, numerical simulation can be carried out using two different concepts to provide a solution [2]. In the first concept, the Reynolds averaged Navier–Stokes (RANS) equations are solved for the mean flow field together with a chosen turbulence closure model. In the second concept, known as the large eddy simulation (LES), the basic Navier–Stokes equations (in primitive variables) are solved and filtered over the minimum grid cells (2D) or volumes (3D), which can be resolved on a given computer. A sub-grid scale (SGS) turbulence closure model can be used to account for the non-resolvable turbulent motion on the grid-scale flow field. Smagorinsky–Lilly model [3] is the most commonly used SGS turbulence closure model. As shown later when the mesh is fine enough, the LES closure model will approach the direct numerical simulation (DNS) model, which still becomes very expensive to modeling even for modern advanced computers. Deardorff [4] first used Smagorinsky’s LES closure model to investigate the numerical simulation of three-dimensional turbulent flows. The model has been commonly used in channel flow and bluff body flow studies (e.g. [5–7]). In the LES modeling, the fluctuating motion of turbulence is treated
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
directly except only neglecting the fine-scale fluctuations contributing to the energy dissipation. The LES concept allows us to evaluate the Reynolds shear stresses on the mesh-scale, and the SGS contributions are negligibly away from rigid walls. Hence the LES model makes it possible to provide information for modelers who used the first concept of the RANS model. The LES modeling also provides a threedimensional and time dependent database, which can be used in further studies of complex interactions of large-scale structures in turbulent flows around obstacles. Furthermore, LES closure model has the capability to freeze the flow and yield velocities, pressures and various turbulence parameters immediately for the entire flow field. Kobayashi et al. [8] obtained solutions in two-dimensional flows, from a LES model for a periodic arrangement of square ribs. Murakami et al. [7] first time applied the LES concept in three-dimensional cases to calculate the flow over a periodic arrangement of cubes, in a simulated atmospheric boundary layer through applying boundary conditions in the main flow direction. The LES closure model will be used in the present paper on two and three dimensions for the simulation of turbulent flows over surface mounted obstacles. The sub-grid scale motions are simulated by using the SGS eddy viscosity model. In general, numerical methods such as finite volume method (FVM), finite difference method (FDM), finite element method (FEM) and boundary element method (BEM) are used to solve multi-dimensional turbulent flows. Song and He [9] used FVM to study three-dimensional wind flow around a tall building by solving weakly compressible flow equations with the Smagorinsky’s sub-grid-scale turbulence model. The numerical simulations of the unsteady incompressible Navier–Stokes equations using FDM have been reported in many studies (for example, Murakami et al. [7], Rail and Moin [10]). Numerous researchers have explored the different applications of FEM [11] and BEM [12] for the solutions of numerous forms of the Navier–Stokes equations in multi-dimensions. Various forms of FEM formulations are available in the literature, for the solutions of incompressible viscous flows. Commonly used Galerkin schemes have limitations to effectively deal with the convective terms, and hence other forms of FEM like Petrov–Galerkin formulations [13] and Taylor–Galerkin schemes [14] were developed. Jiang and Kawahara [15] developed a three-step finite element formulation for the solution of an unsteady incompressible viscous flow based on the Taylor–Galerkin scheme. Different forms of BEM such as direct BEM, indirect BEM and dual reciprocity BEM [12] are also used for the solutions of fluid dynamics problems. In the present investigation, a new computational procedure is developed for the LES closure model of two- and threedimensional turbulent flow problems in external flow field. The model is based on the pressure projection method [16–19] of the unsteady Navier–Stokes equations in the LES terms of pressure and velocity. The momentum equation of the flow domain is solved using an explicit three-step FEM.
565
The pressure Poisson equation for the external flow field is simulated by the BEM. Different external flow problems such as flow over an immersed body and other fluid-structure interaction problems can be solved using the present model through coupling of BEM and three-step FEM within the LES framework. Utilization of BEM for the solution of the pressure Poisson equation helps us to handle the infinite domain of the external flow problem from a finite discrete domain, as the fundamental solutions used in the BEM formulations automatically satisfy the conditions at infinity [16–19]. The use of three-step FEM for the solution of the Navier–Stokes equations is also effective to deal with the convection dominant flows. The developed model has been applied to simulate two multi-dimensional turbulent flows over surface mounted obstacles. A two-dimensional turbulent flow over a square rib and a three-dimensional turbulent flow around a cube were simulated in the Reynolds number range of 4.2!105. The results are in good agreement with other model results even for very coarse meshes used in the simulations. This salient feature of simulation of high-Reynolds number flows using very coarse mesh was also achieved by the previous laminar flow works [18]. After presenting the governing equations and the details of the LES concept, the numerical formulations using coupled three-step FEM and BEM model are briefly described. Then the solution procedures and numerical results for twodimensional and three-dimensional turbulent flow problems are described, and followed by a few concluding remarks. 2. Governing equations The governing equations for transient viscous flows can be described by the continuity and momentum equations in conservative form and are written in Cartesian tensor notation as follows vui Z0 vxi
(1)
vui vui uj 1 vp v2 ui C ZK Cy r vxi vt vxj vxj vxj
(2)
where ui, iZ1,2,3 is the three components of the velocity vector, p is the pressure, y is the kinematic viscosity, t is the time and r is the mass density. Here it is noted that the pressure terms in the momentum equations include the vertical body force, since free surface effects are not taken into consideration. 2.1. Description of large eddy simulation In the present investigation, generally used LES formulation is applied in the model (e.g. [4,7]). In the LES formulation, any physical quantity f is decomposed into two parts f Z f C f 0
(3)
566
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
where f is the resolvable-scale component and f 0 is the sub-grid scale component. f is defined with a filter function G(x,y,z,x 0 ,y 0 ,z 0 ) as follows: ððð f ðx;y;z;tÞ Z 1 Gðx;y;z;x 0 ;y 0 ;z 0 Þf ðx 0 ;y 0 ;z 0 ;tÞdV 0 (4) Dv
eij Z
vu j vu i C vxi vxj
The characteristics filter width can be defined as, D Z ðDx Dy Dz Þ1=3
(12)
Here a sharp spectral filter function [4,7] is used as G(x,y,z,x 0 ,y 0 ,z 0 ) which can be represented as
If D is encompassed by the inertia sub-range scale, from dimensional analysis we can get
Gðx;y;z;x 0 ;y 0 ;z 0 Þ
ysgs Z ðCs DÞ4=3 31=3
Z 1jxKx 0 j!
Dx Dy Dz ; jyKy 0 j! ; jzKz 0 j! 2 2 2
(5)
Gðx;y;z;x 0 ;y 0 ;z 0 Þ Dx Dy Dz ; jyKy 0 jR ; jzKz 0 jR Z 0jxKx jR 2 2 2 0
(6)
where Dx, Dy and Dz are the computational grid size in x-, yand z-directions, respectively. In Eqs. (1) and (2), using Eq. (3) and converting ui as ui Z u i C ui0 and imposing Eq. (4), the filtered equations can be written as follows, vu i Z0 vxi
(7)
vu i vu i uj 1 vp v vu C ZK C Lij C Cij C Rij C y i r vxi vxj vt vxj vxj
(8)
where Lij Z u i uj Kui u j ; Cij ZKui uj0 Kui0 u j ; Rij ZKui0 uj0 In the above terms, the primes are deviations from local grid-volume means, and the term ui0 uj0 thus represents the subgrid-scale (SGS) Reynolds stress tensor. In the present investigation, Reynolds averaging assumption is used, namely, u i uj0 C ui0 u j C u i uj Kui u j Z 0
(9)
Then Eq. (8) can be written as vu i vu i uj 1 vp^ C ZK r vxi vt vxj v 1 0 0 vu i 0 0 C Kui uj C ui uj dij C y vxj 3 vxj
where
where Cs is a numerical parameter called the Smagorinsky constant, 3 is the rate of dissipation of turbulence kinetic energy within the sub-grid-scale. Assuming that production equals dissipation, 3 can be approximated by, 3 zKui0 uj0
vu i vxj
(14)
Substituting Eq. (14) into Eq. (13) and using Eq. (11), ysgs can be expressed as, 1=2 eij (15) ysgs Z ðCs DÞ2 S; S Z 2 The value of Cs depends on the property of the flow field and on the size of the computational mesh. In the present model, the variation of Cs is selected as given by Yoshizawa [20] Cs DS Z 1KCA SK2 Dt Cs0
(16)
where ðD=DtÞZ ðv=vtÞC u$V and CAZ1.8, Cs0Z0.16. The value of Cs generally varies [20] from 0.1 (e.g. for channel flows) to 0.23 (e.g. for decaying turbulence). Using Eq. (11) in Eq. (10), we can finally write the governing equations in the LES form as, vui Z0 vxi
(17)
vui vui u j 1 vp^ v2 ui v vu C ZK Cy C ysgs i r vxi vt vxj vxj vxj vxj vxj
(18)
The following dimensionless forms of the variables are used in the present investigation xi Z xi =L; (10)
where p^ Z pC ðr=3Þui0 uj0 dij and dij is the Kronecker delta. The viscous sub-layer at the wall is not explicitly treated. In the Eq. (10), u i and p^ are the unknown variables whose values are determined at each node of the computational domain. In Eq. (10), the remaining unknowns are the SGS Reynolds stresses ui0 uj0 Kð1=3Þui0 uj0 dij , which are modeled using the SGS eddy viscosity ysgs as 1 ui0 uj0 K ui0 uj0 dij ZKysgs eij 3
(13)
(11)
u i Z u i =u0 ;
t Z tu0 =L;
p Z p^ =ru20
(19)
where L is the characteristic length and u0 is maximum flow velocity. Now Eqs. (17) and (18) can be written as (after dropping the asterisk from the dimensionless variables for brevity from now on) vui Z0 vxi vui vui u j vp 1 v2 ui v C ZK C C vxi Re vxj vxj vxj vt vxj
(20)
1 vu i Resgs vxj
(21)
where ReZLu0/y and ResgsZLu0/ysgs. It is worthwhile to notice that if the inertia sub-range scale is small enough, then RewResgs, and the LES model will approach the DNS model.
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
567
Fig. 1. Computational domain for a sample external turbulent flow problem. (1) Flow separation; (2) Von Karman vortex; (3) Horseshoe vortex; (4) Tip vortex.
2.2. Boundary conditions For the solutions of transient multi-dimensional turbulent external flow problems, appropriate initial and boundary conditions should be prescribed. The initial conditions for the LES Navier–Stokes problems consist of specifying the values of velocity at the initial time as u i ðx;y;z;0Þ Z u i0 ðx;y;zÞ
(22)
v Z 0;
w Z 0
(23)
As shown in Fig. 1, at the inlet plane of the finite flow domain, a uniform flow profile was fixed as, v Z 0; u Z U;
w Z 0
(24)
As far as external flows are concerned, the outer boundaries are located at the infinity and the pressure at infinity is assumed to be zero. Due to the limitations of computational facilities, it is assumed in numerical computations, that the computational domain is limited to a finite region. It is assumed outside the finite region, that the flow is uniform, continuity equation is valid, convective acceleration can be neglected and the pressure gradient in the flow direction is zero [21]. Hence it is assumed that 2
V p Z 0
Cd Z
Fd ; 1=2ru20 A
Fd Z
# p cos f dA C # t sin f dA
Cl Z
(25)
is valid out of the finite computational domain where the fluid is assumed to be undisturbed. As a consequence, it is assumed that the outflow boundary condition of flow field is satisfied by the radiation boundary condition. A sample computational domain for the type of problems considered in the present investigation is shown in Fig. 1. The boundary condition of the fixed body in the flow is set as no-slip boundary and the boundary condition of the pressure will be extracted from the Navier–Stokes equations in the LES form of Eq. (21).
Fl Z
Fl ; 1=2ru20 A
s
St Z
s
s
For the external flow problem considered in the present study, no-slip boundary conditions are imposed on the solid wall surface, that is, u Z 0;
The coefficient of drag and the coefficient of lift on the solid body and the Strouhal number St are found from the following equations, fL u0
(26)
(27)
s
# p sin f dA C # t cos f dA s
s
s
(28)
s
where u0 is the characteristic fluid velocity, Fd is the drag force, r is the mass density, L is the characteristic dimension, Fl is the lift force, f is the frequency of the oscillation, ts is the shear force acting on the body, p s is the pressure acting on the body, A is projected area, f is the angle between the direction of flow and the normal to the surface element. The wall static pressure is found from the following equation Cp Z
ðpKp 0Þ 1=2ru20
(29)
where p0 is the reference pressure. 3. Numerical formulations In the present model as mentioned earlier, a coupled three-step FEM–BEM approach is used in the solution of the governing differential equations in the LES formulation. The numerical formulation is briefly described in this section. 3.1. Three-step FEM formulation In the present model, the mass-momentum Eq. (21) is approximated using an explicit three-step FEM based on
568
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
a Taylor series expansion in time and standard Galerkin FEM in space, or the so called Taylor–Galerkin method as proposed by Jiang and Kawahara [15]. Applying the three-step FEM scheme [15] to Eq. (21), the three-step scheme can be obtained in two phases. The viscous and convective terms are obtained in the first phase, and the pressure terms are derived in the second phase. In the first phase, from Eq. (21), the following equations are obtained as far as time discretization is concerned, Step 1 vuni u nj vpn 1 v2 u ni u inC1=3 Ku ni ZK K C Re vx2j Dt=3 vxj vxi v 1 vuni C vxj Resgs vxj
Mij
Ku nj u nC1=2 1 j ZKAnC1=3 KB ij p nj K Sij ujnC1=3 u nC1=3 ij j Re Dt=2 1 1 1 nC1=3 R C R inC1=3 K S u nC1=3 C Re Resgs ij j Resgs i (35)
For step 3, Mij
u j Ku nj 1 1 ZKAijnC1=2 ujnC1=2 K Sij ujnC1=2 C R inC1=2 Re Re Dt K
(30)
Eqs. (33) can be discretized as
ð
vu inCð1=3Þ u jnCð1=3Þ vp n 1 v2 u inCð1=3Þ ZK K C Re vxj vxi vx2j v 1 vu nC1=3 i C vxj Resgs vxj
Mij Z Ni Nj dU;
(37)
U
R ni Z nCð1=2Þ nCð1=2Þ u j
v C vxj
vxj
1 v2 u inCð1=2Þ Re vx2j nC1=2
(32)
(33)
u jnC1=3 Ku nj Dt=3
1 n R Resgs i
Ni
U
n
vu dS; vn
Z R nC1=2 i
vNi vNj dU; vxk vxk
R inC1=3 Z
ð Ni
vu nC1=3 dS; vn
vU
ð Ni
nC1=2
vu vn
dS
vU
Spatial discretization of Eqs. (30)–(33) are performed by the standard Galerkin criterion using the four-point bilinear elements in two-dimensions and eight-point tri-linear elements in three-dimensions [22]. Here it should be noted that the discretized equation for the u component only is described in the following. The other two components v and w can be derived in a similar way. The resulting finite element equations can be expressed as: For step 1,
ZKAnij unj KB ij p nj K
Sij Z
vU
1 vu i Resgs vxj
vpnC1 vxi
ð
vNj dU; vxk
ð
C
where u i are the apparent velocities from which the velocities in the present time step can be derived as, Z ui KDt u nC1 i
n vNj Z Ni Nk uk dU; vxk ð
U
B ij Z Ni
Step 3 vui u i Ku ni ZK Dt
Anij
U
ð
(31)
C
(36)
where
u inCð1=2Þ Ku ni
Mij
1 1 nC1=2 S unC1=2 C R Resgs ij j Resgs i
Mij ujnC1 Z Mij u j KDtB ij p nC1 j
Step 2
Dt=2
For step 2,
1 1 1 S un C R ni K S un Re ij j Re Resgs ij j (34)
in which Ni, Nj and Nk are the shape functions. After assembling the system and applying the initial and boundary conditions described in Eqs. (22)–(24), the system of equations are solved using the Jacobi iteration. In the Jacobi iteration, we initially assume some values for unknown variables and solve the linear system of equations iteratively using the values from the previous iteration. If a linear system of equations of the form MXZR is to be solved, the system is transformed as kC1 Z RKðMKMÞX k MX
(38)
where M is the consistent coefficient matrix, X is the unknown is vector, R is the known vector, k is the iteration number and M a lumped matrix. Generally three iterations can achieve satisfactory accuracy [15]. Before calculating the velocity in the present time step using Eq. (33), the pressure and its derivatives are to be solved, in the second phase. Combining the continuity Eq. (20) and taking the gradient of Eq. (33) and using the pressure projection method [16–19,23], the pressure Poisson equation is derived to correct the velocity equation as, 1 vu i 2 nC1 V p Z (39) Dt vxi
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
In the present model, the pressure Poisson equation is solved using the BEM formulation, given in the following section using the boundary conditions derived from the solution of Eq. (25) on the infinite domain from the finite discrete domain. The second phase (derivation of pressure terms) is being necessarily implicit obtained from the first phase, and hence an iterative procedure is necessary between the solutions of unC1 ; vnC1 ; w nC1 and p nC1 . 3.2. BEM formulation for pressure Poisson equation Consider the Poisson type pressure equation in p and ui Eq. (39), 1 vui V2 p Z Zb (40) Dt vxj with pressure boundary conditions as, on G1 ;
p Z p0
q Z
vp0 vn
on G2
(41)
where n is the unit outward normal vector and G indicates the boundary of the domain. Here an iterative scheme is used such that the velocity ui is known in the current iteration and time step from the previous step by solving the Navier–Stokes equations in LES form. A weighting function p* can now be introduced such that it has continuous first derivatives within the domain. The following weighted residual statement can now be written ð ð ð qÞp ^ dGK ðpK pÞq ^ dG ðV2 pKbÞp dU Z ðqK (42) G2
U
G1
vp=vn where qZ and q Z vp =vn. Let p* be the fundamental solution of the Laplace equation in three dimensions, represented as p*Z1/(4pr), and in two dimensions as, p*Z(Kln(r)/2p), where r is the distance from the collocation point (k) to other field points (i) given as ffi (in pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi three-dimensions) rZ ðxk Kxi Þ2 C ðyk Kyi Þ2 C ðzk Kzi Þ2 . Now applying Green’s second identity theorem to Eq. (42) and using the standard boundary element procedure [24], we can get the boundary integral equation as ð ð ð dG C bp dU Z qp dG Ci p i C pq (43) G
U
G
vp=vn; where Ci, is the Green’s constant and qZ q Z vp =vn. In Eq. (43), we have boundary integrals and domain integrals. In the present model, domain integration is done by subdividing it into a series of internal cells on each of which a numerical integration is performed. In threedimensions, bilinear elements are used for the boundary discretization and three-dimensional isoparametric quadrilateral volumes are used for the internal discretization. In twodimensions, linear line elements are used for the boundary discretization and two-dimensional isoparametric quadrilateral cells are used for the internal discretization. The details of the element properties, shape functions, coordinate
569
transformation and numerical integration used here are described in Brebbia et al. [24]. If the domain is discretized into M internal cells, then domain integral can be written as " # ð M NI X X wk ðbp Þk Ue (44) Di Z bp dU Z eZ1
U
kZ1
where the integral has been approximated by a summation over different cells or volumes (e varies from 1 to M), wk are the Gauss integration weights. The function (bp*)k needs to be evaluated at integration points k on each cell (k varies from 1 to NI, where NI is the total number of integration points on each cell) and Ue is the volume of cell e. The term Di is the result of the numerical integration and is different for each position i of the boundary nodes. Assuming that the boundary of the domain is discretized into NE linear elements with N nodes, Eq. (43) can be discretized and written in matrix form as, Ci p i C
N X
H ij pj C Di Z
jZ1
N X
Gij q j
(45)
jZ1
Combining the effect of the constant term C with the H matrix, we can write the system of matrix as, H p C D Z Gq
(46)
In Eq. (46), the boundary conditions are introduced and the known values are taken to the right hand side to form a system of linear equations of the form AX Z F
(47)
where X is a vector of unknown boundary values of p and q, and F is a known vector. Eq. (47) is solved by using the Gauss elimination scheme and all the boundary values will be then known. Once this is done, it is possible to calculate the internal values of p or its derivatives. The values of p are calculated at any internal point using the Eq. (43) that can be written in condensed form as: ð ð ð dGK pq dGK bp dU p i Z qp (48) G
G
U
The same discretization is used for the boundary integrals, that is, p i Z
N X jZ1
Gij qj K
N X
H ij p j KDi
(49)
jZ1
where i is an internal point. It should be noted that the same BEM formulation without any domain integral given above (Eq. (43)) is valid for the solution of the Laplace equation given in Eq. (25). Eq. (25) is solved for the infinite domain using the same discretization from the discrete finite domain. The main advantage of using BEM in the solution of the pressure Poisson equation is the effectiveness of BEM to deal with the infinite domain of the external flow problem from
570
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
a finite discrete domain. Here only the boundary condition of pressure on the finite computational domain is known which is used to solve the infinite domain problem. In BEM, the fundamental solutions used in the formulations automatically satisfy the conditions at infinity, and hence it is easier to solve the infinite domain problem from the finite discrete domain [16–19,24]. 4. Solution procedures As mentioned in the previous sections, here an iterative scheme is used in the solution of the Navier–Stokes equations in the LES form. The computational procedures adopted in the present model include the following iterative steps: For the time step nZ1, 0 and solve the (1) Assume at infinite domain, pressure pZ pressure Laplace equation (Eq. (25)) outside the computational domain and pressure Poisson equation (Eq. (39)) inside the computational domain and then get the pressure boundary conditions on the boundaries of the computational domain. These equations are solved simultaneously by utilizing the compatibility and equilibrium conditions of continuity of pressure and flux at the common boundary [16–19,24]. In the problem considered, if U1 represents the computational domain and U2 represents outside the computational domain joined by interface GI, then the compatibility and equilibrium conditions along the interface G1 can be written as p1I Z p2I ; q1I C q2I Z 0
(50)
where p1I ; q1I refer to the pressures and fluxes on the interface GI for the region l (lZ1,2). (2) Solution of the LES form of Navier–Stokes equations using three-step FEM and pressure projection method. † Solve for the unknown apparent velocity values from Eq. (36). † Calculate the pressure distribution for the current time step from the pressure Poisson equation (Eq. (39)) using BEM. † Determine the new velocity values by solving Eq. (37). (3) Check for convergence of the velocity and pressure components in the present iteration, for example, ju nkC1 Ku nk j % 0:00001 junk j
uniform CFL numerical stability condition of the scheme. Various numerical experiments on two-dimensional incompressible viscous flow problem showed that the model is stable and hence a for 0%Cr%1 (CrZCourant numberZuDt=Dx) large time step can be used. A detailed stability analysis for the three-step FEM scheme, used in the present study can be found in Jiang and Kawahara [15]. Here it is noted that, the solution of the pressure Poisson Eq. (39) using BEM includes domain integrals. We have to solve a large system of equations, while finding out the internal pressure distribution using BEM with a fine mesh. As the system matrices are fully populated, the computational costs would increase. Numerical investigations showed that the computing efficiency could be considerably increased by solving the boundary pressure values on the computational domain using the BEM, while the FEM is used to get the internal pressure values. As FEM discretization for the internal domain is already done and matrix properties are available, it is more efficient to use FEM to determine the internal pressure distribution from the boundary pressure already found from using BEM. Both the system matrices of BEM and FEM are coupled using the compatibility and equilibrium conditions as described earlier. 5. Model applications To verify the feasibility of the new model, the proposed LES three-step FEM–BEM model has been applied on two test problems of turbulent flows over surface mounted obstacles. In the Reynolds number range of 4.2!105, a two-dimensional turbulent flow over a square rib and a three-dimensional turbulent flow around a cube were simulated. For the convenience of graphical illustration purposes the bar over physical properties are dropped from now on. 5.1. Turbulent flow over a square rib The three-step FEM–BEM model has been applied on a surface mounted square rib initially and compared the results with the experimental measurement of Crabb et al. [25] and the numerical model with k–3 turbulence closure model results of Benodekar et al. [26]. Fig. 2 shows the form of calculation domain with mesh discretization. The inlet boundary is located
If convergence criterion is satisfied, then we proceed to the next step, otherwise go to step 1. (4) In the successive time step, we use the velocity and pressure components from the previous time step as initial conditions and the new initial guesses, and use the iterative procedure, steps 2–3. The procedure is repeated until the prescribed time step is reached. Detailed numerical investigations showed that the present model with the three-step FEM for the momentum equations is more stable than other FEM schemes like two-step scheme or Lax–Wendroff schemes due to the third-order accuracy and
Fig. 2. Computational domain for flow over a square rib with discretization.
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
Fig. 3. Mean velocity field for flow over a square rib, ReZ420,000.
sufficiently upstream of the obstacle so that its presence is not felt and the outlet is located sufficiently downstream of the point of reattachment as shown in Fig. 2 (by trial and error). The computational domain is chosen with x/DZ26.5 and y/DZ8 (where D is the height of the square rib). The values of the various dimensions and boundary conditions are as given in Crabb et al. [25] and Benodekar et al. [26]. The flow field consists of three zones of flow separation, which are located in the front, behind, and the top of the obstacle. No-slip boundary conditions are imposed along the walls, and uZu0, vZ0 are prescribed along the inlet boundary. At the inlet plane, a power-law boundary layer profile for a smooth wall is imposed on axial velocity as in the measurements of Crabb et al. [25], given by y a u Z (51) u0 d where a is a constant taken as 0.25, d is the boundary layer thickness taken as half the thickness of the rib, u is the mean velocity and u0 is the maximum flow velocity. In the present analysis, the turbulent flow is considered at a Reynolds number of 420,000. For the analysis, the domain considered is discretized using 3204 bilinear quadrilateral elements and 3344 nodes, as shown in Fig. 2. Here a coarse mesh is used to show the effectiveness of the present model to yield reasonable results even for coarser mesh. A time step of 0.025 (nondimensional) is used in the present analysis.
571
Figs. 3–6 show the mean velocity values of the flow field, pressure distributions, streamlines and the vorticity field, respectively, for the turbulent flow at a total time of tZ300. As discerned from Fig. 6, the predicted reattachment length is 10.5 D, which is less than the measured value of 12.3 D. As shown in Eqs. (13)–(16) and also in [20], the reattachment length is dependent on Cs and the mesh selected for the simulation. The under prediction of the reattachment length can be attributed to the coarse mesh used in the simulation. The predicted height of re-circulation of 1.8 D agrees very well with the value inferred from Crabb et al.’s data [25]. As observed in the measurements, the present predictions indicate the presence of a small counterrotating eddy behind the obstacle. Figs. 7 and 8, show the coefficients of drag and lift on the obstacle with reference to time. These two figures display conspicuously the transient nature of the LES turbulence modeling. These two graphs also reveal that there are large scale oscillations captured in the mean flow which could be used for comparison by an unsteady RANS approach. The mean u-velocity profiles at a total time of tZ300 at different axial locations at upstream and downstream are compared, with the results of the LDA measurements of Crabb et al. [25] and the numerical prediction (FVM) of Benodekar et al. [26] in Figs. 9–11. The agreement is quite close at x/DZK0.9, 9.5 and 12.5, especially for the LDA measured values. There is some discrepancy in the upstream and downstream regions close to the rib. The calculated mean wall static pressure distribution at a total time of tZ300 is compared with the LDA measurements of Crabb et al. [25] and the numerical prediction (FVM) of Benodekar et al. [26] as shown in Fig. 12. The correspondence with the LDA measurement is good at upstream of the rib, beginning and end of the downstream. Despite an underprediction occurs at middle of downstream for unclear reasons, the overall comparison with the other predictions is reasonable. Further analysis of computed data of the instantaneous and mean profiles of the velocity, pressure, streamline, and vorticity distribution reveals that fluctuating motions of turbulent mean flows are achieved by the present LES turbulence model. Even in the two-dimensional case studies,
Fig. 4. Mean pressure field for flow over a square rib, ReZ420,000.
572
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
Fig. 5. Mean streamline distribution for flow over a square rib, ReZ420,000.
Fig. 9. Comparison of velocity predictions and measurements on upstream of the square rib.
Fig. 6. Mean vorticity distribution for flow over a square rib, ReZ420,000.
Fig. 7. Variation of drag coefficient with time for flow over a square rib.
Fig. 10. Comparison of velocity predictions and measurements on downstream of the square rib.
this unsteady phenomenon also reflects on the coefficients of drag and lift as a function of time. 5.2. Turbulent flow around a cube
Fig. 8. Variation of lift coefficient with time for flow over a square rib.
The large eddy simulation by three-step FEM–BEM model has been applied for the simulation to the case of flow around a surface mounted cube. The unsteady flow field around a cube is very complex and much more difficult to predict than the above described two-dimensional case. In addition to the vortex structure in the two-dimensional case, there are more complex vortex structures, consisting of four major vortex systems [9]. They are: (1) the secondary currents generated by the flow
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
573
Fig. 14. Velocity field for flow around a cube on x–z plane at yZ0, tZ160.
Fig. 11. Comparison of velocity predictions and measurements on top of the square rib. Fig. 15. Pressure field for flow around a cube on x–z plane at yZ0, tZ160.
Fig. 12. Mean wall static pressure distribution for flow over a square rib. Fig. 16. Vorticity distribution for flow around a cube on x–z plane at yZ0, tZ160.
Fig. 13. Computational domain for flow around a cube with discretization.
Fig. 17. Pressure field for flow around a cube on y–z plane at xZ0, tZ160.
574
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
Fig. 18. Vorticity distribution for flow around a cube on y–z plane at xZ0, tZ160.
separation on the top of the cube and the large circulation cavity behind the cube; (2) von Karman vortex street with axes in the vertical direction; (3) Horseshoe vortex around the structure near the ground; and (4) the twin axial vortices issuing from the side edges of the top of cube. Moreover, the three-dimensional flow always appears to be at transient state. The form of calculation domain with mesh discretization is shown in Fig. 13. As in Fig. 13 (by trial and error) the inlet boundary is located sufficiently upstream of the obstacle so that
its presence is not felt and the outlet is located sufficiently downstream of the point of reattachment. The computational domain in the present study corresponds to x/dZ22.5, y/dZ8 and z/dZ7 (d is the size of the cube). No-slip boundary conditions are imposed along the walls, and uZu0, vZ0, wZ0 are prescribed along the inlet boundary. At the inlet plane, a power-law boundary layer profile for a smooth wall is imposed on axial velocity as in the case of Eq. (51), in which d is the boundary layer thickness taken as half the thickness of the cube. In the present numerical analysis, the turbulent flow is considered at a Reynolds number of 420,000. For the analysis, the domain is discretized into a mesh of 61!37!29, using 58,752 tri-linear cubic elements and 64,001 nodes, as in Fig. 13. Here a fine mesh is used to capture the various details of turbulence. A time step of 0.05 is used in the analysis. As mentioned earlier, the turbulent flow over a cube contains four types of vortices. The time averaged flow patterns from the present numerical investigations are shown in Fig. 1, which are almost same as the earlier numerical as well as experimental results [7,9,27]. As in Fig. 1, the shear layer develops into a large recirculation cavity region 1. Along the obstacle height, the vertical vortices 2 induced by von Karman vortex shedding occur. The horseshoe-shaped vortices 3 occur around the base of the cube and the tip vortex 4 occurs away downstream. Figs. 14–16 represent the velocity vectors, pressure field and vorticity distribution, respectively, at the instant of time tZ160
Fig. 19. Velocity field for flow around a cube on x–y plane at zZ0.5, tZ160.
Fig. 20. Vorticity distribution for flow around a cube on x–y plane at zZ0.5, tZ160.
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
(non-dimensional) at yZ0.0 in the x–z plane. The pressure field and vorticity distributions at the instant time tZ160 (nondimensional) at xZ0.0 in the y–z plane, are represented by Figs. 17 and 18, respectively. Figs. 19 and 20 represent, respectively, the velocity field and vorticity distribution at the instant of time tZ160 at zZ0.5 in the x–y plane. When compared the velocity field to the two-dimensional case, the reattachment length behind the obstacle in Fig. 14 is much shorter than that in Fig. 3 and the separation is much more significant. For the three-dimensional cube case, the front separation is inclined as the flow passes around both sides of the cube, forming the horseshoe vortex structure. It can also be seen that the flow, separating after the front edge, almost no longer reattaches on the top of the cube. Although we do not present a quantitative comparison of various parameters with other LES models, a qualitative comparison of flow and pressure fields with the FVM or FDM results of He and Song [28], Yu and Kareem [29], and Murakami et al. [7] all gives encouraging overall agreement considering the coarse nature of the mesh used. Further observation of the computational results demonstrates the complex characteristics of unsteady and threedimensional nature of the LES turbulence model. The flow structures vary temporally as well as locally. The turbulent features are vividly simulated in the present study, especially for the dynamics of vortex shedding, vortex stretching, and interaction between tip, Karman and horseshoe vortices. These properties all relate to the mechanism of flow separation. They are very difficult to numerical modeling or experimental measurement, due to the unsteady and three-dimensional nature. 6. Concluding remarks A novel computational model has been developed to solve two- and three-dimensional turbulent flow problems in external flow field with large eddy simulation (LES) concept by coupling three-step FEM and BEM. The pressure projection model based on LES form of the Navier–Stokes equations is able to solve the infinite boundary value problems by extracting the boundary effects on a specified finite computational domain. In the proposed model, the Navier–Stokes equations in the LES form are solved using a three-step FEM, and the Poisson type pressure equations are solved using BEM. By coupling the three-step FEM and BEM, the model is able to handle infinite domain problems. The model has been applied to solve a twodimensional turbulent flow over a square rib, and a threedimensional turbulent flow around a surface mounted cube. Both applications give reasonable and comparable results and further show the feasibility of the proposed model as comparing with other models in the literature even for smaller number of grids. Acknowledgements The work reported in this paper was supported by the National Science Council of Taiwan under the Grant no. NSC 93-E-2611-002-001. It is greatly appreciated.
575
References [1] Blevins RD. Flow-induced vibration. 2nd ed. New York: Van Nostrand Reinhold; 1990. [2] Baetke F, Werner H. Numerical simulation of turbulence flow over surface-mounted obstacles with sharp edges and corners. J Wind Eng Ind Aerod 1990;35:129–47. [3] Smagorinsky JS. General circulation experiments with the primitive equations. Part 1. Basic experiments. Mon Weather Rev 1963;91:99–164. [4] Deardorff JW. A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers. J Fluid Mech 1970;41(2):453–80. [5] Schumann U. Sub-grid scale model for finite difference simulation of turbulent flows in plane channels and annuli. J Comput Phys 1975;18: 376–404. [6] Moin P, Kim J. Numerical investigation of turbulent channel flow. J Fluid Mech 1982;118:341–77. [7] Murakami S, Mochida A, Hibi K. Three-dimensional simulation of air flow around a cubic model by means of large eddy simulation. J Wind Eng Ind Aerod 1987;25:291–305. [8] Kobayashi T, Ishiara T, Kano M. Prediction of turbulent flow in twodimensional channel with turbulence promoters, III— improvement of large eddy simulation and formation of streaklines. JSME Bull 1985; 28(246):2948–53. [9] Song CCS, He J. Computation of wind flow around a tall building and the large-scale vortex structure. J Wind Eng Ind Aerod 1993;46–47:219–28. [10] Rail MM, Moin P. Direct simulations of turbulent flow using finite difference schemes. J Comput Phys 1991;96:15–53. [11] Reddy JN, Gartling DK. The finite element method in heat transfer and fluid dynamics. London: CRC Press; 1992. [12] Power H, Wrobel LC. Boundary integral methods in fluid mechanics. Southampton: Computational Mechanics Publications; 1995. [13] Brookes AN, Hughes TJR. Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations. Comput Meth Appl Mech Eng 1982;32:199–259. [14] Lohner R, Morgan K, Zienkiewics OC. The solution of nonlinear hyperbolic equations systems by the finite element method. Int J Numer Methods Fluids 1984;4:1043–63. [15] Jiang CB, Kawahara M. The analysis of unsteady incompressible flows by a three-step finite element method. Int J Numer Methods Fluids 1993;16: 793–811. [16] Young DL, Chang JT, Eldho TI. A coupled BEM and arbitrary Lagrangian–Eulerian FEM model for the solution of two-dimensional laminar flows in external flow fields. Int J Numer Meth Eng 2001;(9): 1053–77. [17] Young DL, Huang JL, Eldho TI. Simulation of laminar vortex shedding flow past cylinders using a coupled BEM and FEM model. Comput Meth Appl Mech Eng 2001;190:5975–98. [18] Young DL, Huang JL, Eldho TI. Numerical simulation of high-Reynolds number flow around circular cylinders by a three-step FEM–BEM model. Int J Numer Methods Fluids 2001;37:657–89. [19] Young DL, Chang JT, Eldho TI. Solution of three-dimensional unsteady external flow using a coupled arbitrary Lagrangian FEM–BEM model. Eng Anal Bound Elem 2004;28:711–23. [20] Yoshizawa A. Eddy-viscosity-type subgrid-scale model with a variable Smagorinsky coefficient and its relationship with the oneequation model in large eddy simulation. J Phys Fluids 1991;A3: 2007–9. [21] Debler WA. Fluid mechanics fundamentals. Englewood Cliffs, NJ: Prentice Hall; 1990. [22] Reddy JN. Finite element method. New York: McGraw-Hill; 1993. [23] Young DL, Liao CB, Sheen HJ. Computations of recirculation zones of a confined annular swirling flow. Int J Numer Methods Fluids 1999;29: 791–810. [24] Brebbia CA, Telles JCF, Wrobel LC. Boundary element techniques— theory and applications in engineering. Berlin: Springer; 1984.
576
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
[25] Crabb D, Durao DFG, Whitelaw JH. Velocity characteristics in the vicinity of a two dimensional rib Proceedings of fourth Brazilian congress on mechanical engineering, Florianopolis, Brazil, December; 1977. p. 415–29. [26] Benodekar RW, Goddard AJH, Gosman AD, Issac RI. Numerical prediction of turbulent flow over surface-mounted ribs. AIAA J 1985; 23(3):359–66.
[27] Castro IP, Robins AG. The flow around a surface mounted cube in uniform and turbulent streams. J Fluid Mech 1977;79(2):307–35. [28] He J, Song CCS. Computation of turbulent shear flow over surface-mounted obstacle. J Eng Mech ASCE 1992;118(11): 2282–97. [29] Yu D, Kareem A. Parametric study of flow around rectangular prisms using LES. J Wind Eng Ind Aerod 1998;77/78:653–62.
Large eddy simulation of turbulent flows in external flow field using three-step FEM–BEM model D.L. Young a,*, T.I. Eldho b, J.T. Chang a a
Department of Civil Engineering and Hydrotech Research Institute, National Taiwan University, 1, Sec. 4, Roosvelt Road, Taipei 10617, Taiwan b Department of Civil Engineering, Indian Institute of Technology, Mumbai 400 076, India Received 28 March 2005; accepted 2 February 2006 Available online 19 April 2006
Abstract An innovative computational model is presented for the large eddy simulation (LES) modeling of multi-dimensional unsteady turbulent flow problems in external flow field. Based on the LES principles, the model uses a pressure projection method to solve the Navier–Stokes equations in transient condition. The turbulent motion is simulated by Smagorinsky sub-grid scale (SGS) eddy viscosity model. The momentum equation of the flow motion is solved using a three-step finite element method (FEM). The external flow field is simulated using a boundary element method (BEM) by solving a pressure Poisson equation that assumes the pressure as zero at the infinity. Through extracting the boundary effects on a specified finite computational domain, the model is able to solve the infinite boundary value problems. The present model is used to simulate the flows past a two-dimensional square rib and a three-dimensional cube at high Reynolds number. The simulation results are found to be reasonable and comparable with other models available in the literature even for coarse meshes. q 2006 Elsevier Ltd. All rights reserved. Keywords: Large eddy simulation; Turbulent flows; Multi-dimensional external flows; Navier–Stokes equations; Three-step finite element method; Boundary element method
1. Introduction Three-dimensional turbulent, external flows and their interaction with structures play a vital role in the design of high-rise buildings, suspended bridges and aircrafts. Many researchers have investigated inherently three-dimensional separation zones around structures in a surface boundary layer. Experimental investigations carried out to simulate these external flow phenomena, provided limited information for the entire flow field. Hence numerical models are essential for predicting the turbulent flow fields around buildings and other complex flows. One difficulty in simulating these flows lies in effectively setting up the boundary conditions and also the computational domain to solve external flows with infinite domain. In these types of problems, the vortex shedding around the structure plays an important role in characterizing various types of fluid-induced structural vibrations [1]. To develop a computational model for this class of turbulent flow problems, it becomes necessary to deal with the unsteady Navier–Stokes * Corresponding author. Tel./fax: C886 2 2362 6114. E-mail address: [email protected] (D.L. Young).
0955-7997/$ - see front matter q 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2006.02.004
equations with a turbulence closure form such as the LES model in a multi-dimensional infinite domain. To study the turbulent flows in external flow field, numerical simulation can be carried out using two different concepts to provide a solution [2]. In the first concept, the Reynolds averaged Navier–Stokes (RANS) equations are solved for the mean flow field together with a chosen turbulence closure model. In the second concept, known as the large eddy simulation (LES), the basic Navier–Stokes equations (in primitive variables) are solved and filtered over the minimum grid cells (2D) or volumes (3D), which can be resolved on a given computer. A sub-grid scale (SGS) turbulence closure model can be used to account for the non-resolvable turbulent motion on the grid-scale flow field. Smagorinsky–Lilly model [3] is the most commonly used SGS turbulence closure model. As shown later when the mesh is fine enough, the LES closure model will approach the direct numerical simulation (DNS) model, which still becomes very expensive to modeling even for modern advanced computers. Deardorff [4] first used Smagorinsky’s LES closure model to investigate the numerical simulation of three-dimensional turbulent flows. The model has been commonly used in channel flow and bluff body flow studies (e.g. [5–7]). In the LES modeling, the fluctuating motion of turbulence is treated
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
directly except only neglecting the fine-scale fluctuations contributing to the energy dissipation. The LES concept allows us to evaluate the Reynolds shear stresses on the mesh-scale, and the SGS contributions are negligibly away from rigid walls. Hence the LES model makes it possible to provide information for modelers who used the first concept of the RANS model. The LES modeling also provides a threedimensional and time dependent database, which can be used in further studies of complex interactions of large-scale structures in turbulent flows around obstacles. Furthermore, LES closure model has the capability to freeze the flow and yield velocities, pressures and various turbulence parameters immediately for the entire flow field. Kobayashi et al. [8] obtained solutions in two-dimensional flows, from a LES model for a periodic arrangement of square ribs. Murakami et al. [7] first time applied the LES concept in three-dimensional cases to calculate the flow over a periodic arrangement of cubes, in a simulated atmospheric boundary layer through applying boundary conditions in the main flow direction. The LES closure model will be used in the present paper on two and three dimensions for the simulation of turbulent flows over surface mounted obstacles. The sub-grid scale motions are simulated by using the SGS eddy viscosity model. In general, numerical methods such as finite volume method (FVM), finite difference method (FDM), finite element method (FEM) and boundary element method (BEM) are used to solve multi-dimensional turbulent flows. Song and He [9] used FVM to study three-dimensional wind flow around a tall building by solving weakly compressible flow equations with the Smagorinsky’s sub-grid-scale turbulence model. The numerical simulations of the unsteady incompressible Navier–Stokes equations using FDM have been reported in many studies (for example, Murakami et al. [7], Rail and Moin [10]). Numerous researchers have explored the different applications of FEM [11] and BEM [12] for the solutions of numerous forms of the Navier–Stokes equations in multi-dimensions. Various forms of FEM formulations are available in the literature, for the solutions of incompressible viscous flows. Commonly used Galerkin schemes have limitations to effectively deal with the convective terms, and hence other forms of FEM like Petrov–Galerkin formulations [13] and Taylor–Galerkin schemes [14] were developed. Jiang and Kawahara [15] developed a three-step finite element formulation for the solution of an unsteady incompressible viscous flow based on the Taylor–Galerkin scheme. Different forms of BEM such as direct BEM, indirect BEM and dual reciprocity BEM [12] are also used for the solutions of fluid dynamics problems. In the present investigation, a new computational procedure is developed for the LES closure model of two- and threedimensional turbulent flow problems in external flow field. The model is based on the pressure projection method [16–19] of the unsteady Navier–Stokes equations in the LES terms of pressure and velocity. The momentum equation of the flow domain is solved using an explicit three-step FEM.
565
The pressure Poisson equation for the external flow field is simulated by the BEM. Different external flow problems such as flow over an immersed body and other fluid-structure interaction problems can be solved using the present model through coupling of BEM and three-step FEM within the LES framework. Utilization of BEM for the solution of the pressure Poisson equation helps us to handle the infinite domain of the external flow problem from a finite discrete domain, as the fundamental solutions used in the BEM formulations automatically satisfy the conditions at infinity [16–19]. The use of three-step FEM for the solution of the Navier–Stokes equations is also effective to deal with the convection dominant flows. The developed model has been applied to simulate two multi-dimensional turbulent flows over surface mounted obstacles. A two-dimensional turbulent flow over a square rib and a three-dimensional turbulent flow around a cube were simulated in the Reynolds number range of 4.2!105. The results are in good agreement with other model results even for very coarse meshes used in the simulations. This salient feature of simulation of high-Reynolds number flows using very coarse mesh was also achieved by the previous laminar flow works [18]. After presenting the governing equations and the details of the LES concept, the numerical formulations using coupled three-step FEM and BEM model are briefly described. Then the solution procedures and numerical results for twodimensional and three-dimensional turbulent flow problems are described, and followed by a few concluding remarks. 2. Governing equations The governing equations for transient viscous flows can be described by the continuity and momentum equations in conservative form and are written in Cartesian tensor notation as follows vui Z0 vxi
(1)
vui vui uj 1 vp v2 ui C ZK Cy r vxi vt vxj vxj vxj
(2)
where ui, iZ1,2,3 is the three components of the velocity vector, p is the pressure, y is the kinematic viscosity, t is the time and r is the mass density. Here it is noted that the pressure terms in the momentum equations include the vertical body force, since free surface effects are not taken into consideration. 2.1. Description of large eddy simulation In the present investigation, generally used LES formulation is applied in the model (e.g. [4,7]). In the LES formulation, any physical quantity f is decomposed into two parts f Z f C f 0
(3)
566
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
where f is the resolvable-scale component and f 0 is the sub-grid scale component. f is defined with a filter function G(x,y,z,x 0 ,y 0 ,z 0 ) as follows: ððð f ðx;y;z;tÞ Z 1 Gðx;y;z;x 0 ;y 0 ;z 0 Þf ðx 0 ;y 0 ;z 0 ;tÞdV 0 (4) Dv
eij Z
vu j vu i C vxi vxj
The characteristics filter width can be defined as, D Z ðDx Dy Dz Þ1=3
(12)
Here a sharp spectral filter function [4,7] is used as G(x,y,z,x 0 ,y 0 ,z 0 ) which can be represented as
If D is encompassed by the inertia sub-range scale, from dimensional analysis we can get
Gðx;y;z;x 0 ;y 0 ;z 0 Þ
ysgs Z ðCs DÞ4=3 31=3
Z 1jxKx 0 j!
Dx Dy Dz ; jyKy 0 j! ; jzKz 0 j! 2 2 2
(5)
Gðx;y;z;x 0 ;y 0 ;z 0 Þ Dx Dy Dz ; jyKy 0 jR ; jzKz 0 jR Z 0jxKx jR 2 2 2 0
(6)
where Dx, Dy and Dz are the computational grid size in x-, yand z-directions, respectively. In Eqs. (1) and (2), using Eq. (3) and converting ui as ui Z u i C ui0 and imposing Eq. (4), the filtered equations can be written as follows, vu i Z0 vxi
(7)
vu i vu i uj 1 vp v vu C ZK C Lij C Cij C Rij C y i r vxi vxj vt vxj vxj
(8)
where Lij Z u i uj Kui u j ; Cij ZKui uj0 Kui0 u j ; Rij ZKui0 uj0 In the above terms, the primes are deviations from local grid-volume means, and the term ui0 uj0 thus represents the subgrid-scale (SGS) Reynolds stress tensor. In the present investigation, Reynolds averaging assumption is used, namely, u i uj0 C ui0 u j C u i uj Kui u j Z 0
(9)
Then Eq. (8) can be written as vu i vu i uj 1 vp^ C ZK r vxi vt vxj v 1 0 0 vu i 0 0 C Kui uj C ui uj dij C y vxj 3 vxj
where
where Cs is a numerical parameter called the Smagorinsky constant, 3 is the rate of dissipation of turbulence kinetic energy within the sub-grid-scale. Assuming that production equals dissipation, 3 can be approximated by, 3 zKui0 uj0
vu i vxj
(14)
Substituting Eq. (14) into Eq. (13) and using Eq. (11), ysgs can be expressed as, 1=2 eij (15) ysgs Z ðCs DÞ2 S; S Z 2 The value of Cs depends on the property of the flow field and on the size of the computational mesh. In the present model, the variation of Cs is selected as given by Yoshizawa [20] Cs DS Z 1KCA SK2 Dt Cs0
(16)
where ðD=DtÞZ ðv=vtÞC u$V and CAZ1.8, Cs0Z0.16. The value of Cs generally varies [20] from 0.1 (e.g. for channel flows) to 0.23 (e.g. for decaying turbulence). Using Eq. (11) in Eq. (10), we can finally write the governing equations in the LES form as, vui Z0 vxi
(17)
vui vui u j 1 vp^ v2 ui v vu C ZK Cy C ysgs i r vxi vt vxj vxj vxj vxj vxj
(18)
The following dimensionless forms of the variables are used in the present investigation xi Z xi =L; (10)
where p^ Z pC ðr=3Þui0 uj0 dij and dij is the Kronecker delta. The viscous sub-layer at the wall is not explicitly treated. In the Eq. (10), u i and p^ are the unknown variables whose values are determined at each node of the computational domain. In Eq. (10), the remaining unknowns are the SGS Reynolds stresses ui0 uj0 Kð1=3Þui0 uj0 dij , which are modeled using the SGS eddy viscosity ysgs as 1 ui0 uj0 K ui0 uj0 dij ZKysgs eij 3
(13)
(11)
u i Z u i =u0 ;
t Z tu0 =L;
p Z p^ =ru20
(19)
where L is the characteristic length and u0 is maximum flow velocity. Now Eqs. (17) and (18) can be written as (after dropping the asterisk from the dimensionless variables for brevity from now on) vui Z0 vxi vui vui u j vp 1 v2 ui v C ZK C C vxi Re vxj vxj vxj vt vxj
(20)
1 vu i Resgs vxj
(21)
where ReZLu0/y and ResgsZLu0/ysgs. It is worthwhile to notice that if the inertia sub-range scale is small enough, then RewResgs, and the LES model will approach the DNS model.
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
567
Fig. 1. Computational domain for a sample external turbulent flow problem. (1) Flow separation; (2) Von Karman vortex; (3) Horseshoe vortex; (4) Tip vortex.
2.2. Boundary conditions For the solutions of transient multi-dimensional turbulent external flow problems, appropriate initial and boundary conditions should be prescribed. The initial conditions for the LES Navier–Stokes problems consist of specifying the values of velocity at the initial time as u i ðx;y;z;0Þ Z u i0 ðx;y;zÞ
(22)
v Z 0;
w Z 0
(23)
As shown in Fig. 1, at the inlet plane of the finite flow domain, a uniform flow profile was fixed as, v Z 0; u Z U;
w Z 0
(24)
As far as external flows are concerned, the outer boundaries are located at the infinity and the pressure at infinity is assumed to be zero. Due to the limitations of computational facilities, it is assumed in numerical computations, that the computational domain is limited to a finite region. It is assumed outside the finite region, that the flow is uniform, continuity equation is valid, convective acceleration can be neglected and the pressure gradient in the flow direction is zero [21]. Hence it is assumed that 2
V p Z 0
Cd Z
Fd ; 1=2ru20 A
Fd Z
# p cos f dA C # t sin f dA
Cl Z
(25)
is valid out of the finite computational domain where the fluid is assumed to be undisturbed. As a consequence, it is assumed that the outflow boundary condition of flow field is satisfied by the radiation boundary condition. A sample computational domain for the type of problems considered in the present investigation is shown in Fig. 1. The boundary condition of the fixed body in the flow is set as no-slip boundary and the boundary condition of the pressure will be extracted from the Navier–Stokes equations in the LES form of Eq. (21).
Fl Z
Fl ; 1=2ru20 A
s
St Z
s
s
For the external flow problem considered in the present study, no-slip boundary conditions are imposed on the solid wall surface, that is, u Z 0;
The coefficient of drag and the coefficient of lift on the solid body and the Strouhal number St are found from the following equations, fL u0
(26)
(27)
s
# p sin f dA C # t cos f dA s
s
s
(28)
s
where u0 is the characteristic fluid velocity, Fd is the drag force, r is the mass density, L is the characteristic dimension, Fl is the lift force, f is the frequency of the oscillation, ts is the shear force acting on the body, p s is the pressure acting on the body, A is projected area, f is the angle between the direction of flow and the normal to the surface element. The wall static pressure is found from the following equation Cp Z
ðpKp 0Þ 1=2ru20
(29)
where p0 is the reference pressure. 3. Numerical formulations In the present model as mentioned earlier, a coupled three-step FEM–BEM approach is used in the solution of the governing differential equations in the LES formulation. The numerical formulation is briefly described in this section. 3.1. Three-step FEM formulation In the present model, the mass-momentum Eq. (21) is approximated using an explicit three-step FEM based on
568
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
a Taylor series expansion in time and standard Galerkin FEM in space, or the so called Taylor–Galerkin method as proposed by Jiang and Kawahara [15]. Applying the three-step FEM scheme [15] to Eq. (21), the three-step scheme can be obtained in two phases. The viscous and convective terms are obtained in the first phase, and the pressure terms are derived in the second phase. In the first phase, from Eq. (21), the following equations are obtained as far as time discretization is concerned, Step 1 vuni u nj vpn 1 v2 u ni u inC1=3 Ku ni ZK K C Re vx2j Dt=3 vxj vxi v 1 vuni C vxj Resgs vxj
Mij
Ku nj u nC1=2 1 j ZKAnC1=3 KB ij p nj K Sij ujnC1=3 u nC1=3 ij j Re Dt=2 1 1 1 nC1=3 R C R inC1=3 K S u nC1=3 C Re Resgs ij j Resgs i (35)
For step 3, Mij
u j Ku nj 1 1 ZKAijnC1=2 ujnC1=2 K Sij ujnC1=2 C R inC1=2 Re Re Dt K
(30)
Eqs. (33) can be discretized as
ð
vu inCð1=3Þ u jnCð1=3Þ vp n 1 v2 u inCð1=3Þ ZK K C Re vxj vxi vx2j v 1 vu nC1=3 i C vxj Resgs vxj
Mij Z Ni Nj dU;
(37)
U
R ni Z nCð1=2Þ nCð1=2Þ u j
v C vxj
vxj
1 v2 u inCð1=2Þ Re vx2j nC1=2
(32)
(33)
u jnC1=3 Ku nj Dt=3
1 n R Resgs i
Ni
U
n
vu dS; vn
Z R nC1=2 i
vNi vNj dU; vxk vxk
R inC1=3 Z
ð Ni
vu nC1=3 dS; vn
vU
ð Ni
nC1=2
vu vn
dS
vU
Spatial discretization of Eqs. (30)–(33) are performed by the standard Galerkin criterion using the four-point bilinear elements in two-dimensions and eight-point tri-linear elements in three-dimensions [22]. Here it should be noted that the discretized equation for the u component only is described in the following. The other two components v and w can be derived in a similar way. The resulting finite element equations can be expressed as: For step 1,
ZKAnij unj KB ij p nj K
Sij Z
vU
1 vu i Resgs vxj
vpnC1 vxi
ð
vNj dU; vxk
ð
C
where u i are the apparent velocities from which the velocities in the present time step can be derived as, Z ui KDt u nC1 i
n vNj Z Ni Nk uk dU; vxk ð
U
B ij Z Ni
Step 3 vui u i Ku ni ZK Dt
Anij
U
ð
(31)
C
(36)
where
u inCð1=2Þ Ku ni
Mij
1 1 nC1=2 S unC1=2 C R Resgs ij j Resgs i
Mij ujnC1 Z Mij u j KDtB ij p nC1 j
Step 2
Dt=2
For step 2,
1 1 1 S un C R ni K S un Re ij j Re Resgs ij j (34)
in which Ni, Nj and Nk are the shape functions. After assembling the system and applying the initial and boundary conditions described in Eqs. (22)–(24), the system of equations are solved using the Jacobi iteration. In the Jacobi iteration, we initially assume some values for unknown variables and solve the linear system of equations iteratively using the values from the previous iteration. If a linear system of equations of the form MXZR is to be solved, the system is transformed as kC1 Z RKðMKMÞX k MX
(38)
where M is the consistent coefficient matrix, X is the unknown is vector, R is the known vector, k is the iteration number and M a lumped matrix. Generally three iterations can achieve satisfactory accuracy [15]. Before calculating the velocity in the present time step using Eq. (33), the pressure and its derivatives are to be solved, in the second phase. Combining the continuity Eq. (20) and taking the gradient of Eq. (33) and using the pressure projection method [16–19,23], the pressure Poisson equation is derived to correct the velocity equation as, 1 vu i 2 nC1 V p Z (39) Dt vxi
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
In the present model, the pressure Poisson equation is solved using the BEM formulation, given in the following section using the boundary conditions derived from the solution of Eq. (25) on the infinite domain from the finite discrete domain. The second phase (derivation of pressure terms) is being necessarily implicit obtained from the first phase, and hence an iterative procedure is necessary between the solutions of unC1 ; vnC1 ; w nC1 and p nC1 . 3.2. BEM formulation for pressure Poisson equation Consider the Poisson type pressure equation in p and ui Eq. (39), 1 vui V2 p Z Zb (40) Dt vxj with pressure boundary conditions as, on G1 ;
p Z p0
q Z
vp0 vn
on G2
(41)
where n is the unit outward normal vector and G indicates the boundary of the domain. Here an iterative scheme is used such that the velocity ui is known in the current iteration and time step from the previous step by solving the Navier–Stokes equations in LES form. A weighting function p* can now be introduced such that it has continuous first derivatives within the domain. The following weighted residual statement can now be written ð ð ð qÞp ^ dGK ðpK pÞq ^ dG ðV2 pKbÞp dU Z ðqK (42) G2
U
G1
vp=vn where qZ and q Z vp =vn. Let p* be the fundamental solution of the Laplace equation in three dimensions, represented as p*Z1/(4pr), and in two dimensions as, p*Z(Kln(r)/2p), where r is the distance from the collocation point (k) to other field points (i) given as ffi (in pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi three-dimensions) rZ ðxk Kxi Þ2 C ðyk Kyi Þ2 C ðzk Kzi Þ2 . Now applying Green’s second identity theorem to Eq. (42) and using the standard boundary element procedure [24], we can get the boundary integral equation as ð ð ð dG C bp dU Z qp dG Ci p i C pq (43) G
U
G
vp=vn; where Ci, is the Green’s constant and qZ q Z vp =vn. In Eq. (43), we have boundary integrals and domain integrals. In the present model, domain integration is done by subdividing it into a series of internal cells on each of which a numerical integration is performed. In threedimensions, bilinear elements are used for the boundary discretization and three-dimensional isoparametric quadrilateral volumes are used for the internal discretization. In twodimensions, linear line elements are used for the boundary discretization and two-dimensional isoparametric quadrilateral cells are used for the internal discretization. The details of the element properties, shape functions, coordinate
569
transformation and numerical integration used here are described in Brebbia et al. [24]. If the domain is discretized into M internal cells, then domain integral can be written as " # ð M NI X X wk ðbp Þk Ue (44) Di Z bp dU Z eZ1
U
kZ1
where the integral has been approximated by a summation over different cells or volumes (e varies from 1 to M), wk are the Gauss integration weights. The function (bp*)k needs to be evaluated at integration points k on each cell (k varies from 1 to NI, where NI is the total number of integration points on each cell) and Ue is the volume of cell e. The term Di is the result of the numerical integration and is different for each position i of the boundary nodes. Assuming that the boundary of the domain is discretized into NE linear elements with N nodes, Eq. (43) can be discretized and written in matrix form as, Ci p i C
N X
H ij pj C Di Z
jZ1
N X
Gij q j
(45)
jZ1
Combining the effect of the constant term C with the H matrix, we can write the system of matrix as, H p C D Z Gq
(46)
In Eq. (46), the boundary conditions are introduced and the known values are taken to the right hand side to form a system of linear equations of the form AX Z F
(47)
where X is a vector of unknown boundary values of p and q, and F is a known vector. Eq. (47) is solved by using the Gauss elimination scheme and all the boundary values will be then known. Once this is done, it is possible to calculate the internal values of p or its derivatives. The values of p are calculated at any internal point using the Eq. (43) that can be written in condensed form as: ð ð ð dGK pq dGK bp dU p i Z qp (48) G
G
U
The same discretization is used for the boundary integrals, that is, p i Z
N X jZ1
Gij qj K
N X
H ij p j KDi
(49)
jZ1
where i is an internal point. It should be noted that the same BEM formulation without any domain integral given above (Eq. (43)) is valid for the solution of the Laplace equation given in Eq. (25). Eq. (25) is solved for the infinite domain using the same discretization from the discrete finite domain. The main advantage of using BEM in the solution of the pressure Poisson equation is the effectiveness of BEM to deal with the infinite domain of the external flow problem from
570
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
a finite discrete domain. Here only the boundary condition of pressure on the finite computational domain is known which is used to solve the infinite domain problem. In BEM, the fundamental solutions used in the formulations automatically satisfy the conditions at infinity, and hence it is easier to solve the infinite domain problem from the finite discrete domain [16–19,24]. 4. Solution procedures As mentioned in the previous sections, here an iterative scheme is used in the solution of the Navier–Stokes equations in the LES form. The computational procedures adopted in the present model include the following iterative steps: For the time step nZ1, 0 and solve the (1) Assume at infinite domain, pressure pZ pressure Laplace equation (Eq. (25)) outside the computational domain and pressure Poisson equation (Eq. (39)) inside the computational domain and then get the pressure boundary conditions on the boundaries of the computational domain. These equations are solved simultaneously by utilizing the compatibility and equilibrium conditions of continuity of pressure and flux at the common boundary [16–19,24]. In the problem considered, if U1 represents the computational domain and U2 represents outside the computational domain joined by interface GI, then the compatibility and equilibrium conditions along the interface G1 can be written as p1I Z p2I ; q1I C q2I Z 0
(50)
where p1I ; q1I refer to the pressures and fluxes on the interface GI for the region l (lZ1,2). (2) Solution of the LES form of Navier–Stokes equations using three-step FEM and pressure projection method. † Solve for the unknown apparent velocity values from Eq. (36). † Calculate the pressure distribution for the current time step from the pressure Poisson equation (Eq. (39)) using BEM. † Determine the new velocity values by solving Eq. (37). (3) Check for convergence of the velocity and pressure components in the present iteration, for example, ju nkC1 Ku nk j % 0:00001 junk j
uniform CFL numerical stability condition of the scheme. Various numerical experiments on two-dimensional incompressible viscous flow problem showed that the model is stable and hence a for 0%Cr%1 (CrZCourant numberZuDt=Dx) large time step can be used. A detailed stability analysis for the three-step FEM scheme, used in the present study can be found in Jiang and Kawahara [15]. Here it is noted that, the solution of the pressure Poisson Eq. (39) using BEM includes domain integrals. We have to solve a large system of equations, while finding out the internal pressure distribution using BEM with a fine mesh. As the system matrices are fully populated, the computational costs would increase. Numerical investigations showed that the computing efficiency could be considerably increased by solving the boundary pressure values on the computational domain using the BEM, while the FEM is used to get the internal pressure values. As FEM discretization for the internal domain is already done and matrix properties are available, it is more efficient to use FEM to determine the internal pressure distribution from the boundary pressure already found from using BEM. Both the system matrices of BEM and FEM are coupled using the compatibility and equilibrium conditions as described earlier. 5. Model applications To verify the feasibility of the new model, the proposed LES three-step FEM–BEM model has been applied on two test problems of turbulent flows over surface mounted obstacles. In the Reynolds number range of 4.2!105, a two-dimensional turbulent flow over a square rib and a three-dimensional turbulent flow around a cube were simulated. For the convenience of graphical illustration purposes the bar over physical properties are dropped from now on. 5.1. Turbulent flow over a square rib The three-step FEM–BEM model has been applied on a surface mounted square rib initially and compared the results with the experimental measurement of Crabb et al. [25] and the numerical model with k–3 turbulence closure model results of Benodekar et al. [26]. Fig. 2 shows the form of calculation domain with mesh discretization. The inlet boundary is located
If convergence criterion is satisfied, then we proceed to the next step, otherwise go to step 1. (4) In the successive time step, we use the velocity and pressure components from the previous time step as initial conditions and the new initial guesses, and use the iterative procedure, steps 2–3. The procedure is repeated until the prescribed time step is reached. Detailed numerical investigations showed that the present model with the three-step FEM for the momentum equations is more stable than other FEM schemes like two-step scheme or Lax–Wendroff schemes due to the third-order accuracy and
Fig. 2. Computational domain for flow over a square rib with discretization.
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
Fig. 3. Mean velocity field for flow over a square rib, ReZ420,000.
sufficiently upstream of the obstacle so that its presence is not felt and the outlet is located sufficiently downstream of the point of reattachment as shown in Fig. 2 (by trial and error). The computational domain is chosen with x/DZ26.5 and y/DZ8 (where D is the height of the square rib). The values of the various dimensions and boundary conditions are as given in Crabb et al. [25] and Benodekar et al. [26]. The flow field consists of three zones of flow separation, which are located in the front, behind, and the top of the obstacle. No-slip boundary conditions are imposed along the walls, and uZu0, vZ0 are prescribed along the inlet boundary. At the inlet plane, a power-law boundary layer profile for a smooth wall is imposed on axial velocity as in the measurements of Crabb et al. [25], given by y a u Z (51) u0 d where a is a constant taken as 0.25, d is the boundary layer thickness taken as half the thickness of the rib, u is the mean velocity and u0 is the maximum flow velocity. In the present analysis, the turbulent flow is considered at a Reynolds number of 420,000. For the analysis, the domain considered is discretized using 3204 bilinear quadrilateral elements and 3344 nodes, as shown in Fig. 2. Here a coarse mesh is used to show the effectiveness of the present model to yield reasonable results even for coarser mesh. A time step of 0.025 (nondimensional) is used in the present analysis.
571
Figs. 3–6 show the mean velocity values of the flow field, pressure distributions, streamlines and the vorticity field, respectively, for the turbulent flow at a total time of tZ300. As discerned from Fig. 6, the predicted reattachment length is 10.5 D, which is less than the measured value of 12.3 D. As shown in Eqs. (13)–(16) and also in [20], the reattachment length is dependent on Cs and the mesh selected for the simulation. The under prediction of the reattachment length can be attributed to the coarse mesh used in the simulation. The predicted height of re-circulation of 1.8 D agrees very well with the value inferred from Crabb et al.’s data [25]. As observed in the measurements, the present predictions indicate the presence of a small counterrotating eddy behind the obstacle. Figs. 7 and 8, show the coefficients of drag and lift on the obstacle with reference to time. These two figures display conspicuously the transient nature of the LES turbulence modeling. These two graphs also reveal that there are large scale oscillations captured in the mean flow which could be used for comparison by an unsteady RANS approach. The mean u-velocity profiles at a total time of tZ300 at different axial locations at upstream and downstream are compared, with the results of the LDA measurements of Crabb et al. [25] and the numerical prediction (FVM) of Benodekar et al. [26] in Figs. 9–11. The agreement is quite close at x/DZK0.9, 9.5 and 12.5, especially for the LDA measured values. There is some discrepancy in the upstream and downstream regions close to the rib. The calculated mean wall static pressure distribution at a total time of tZ300 is compared with the LDA measurements of Crabb et al. [25] and the numerical prediction (FVM) of Benodekar et al. [26] as shown in Fig. 12. The correspondence with the LDA measurement is good at upstream of the rib, beginning and end of the downstream. Despite an underprediction occurs at middle of downstream for unclear reasons, the overall comparison with the other predictions is reasonable. Further analysis of computed data of the instantaneous and mean profiles of the velocity, pressure, streamline, and vorticity distribution reveals that fluctuating motions of turbulent mean flows are achieved by the present LES turbulence model. Even in the two-dimensional case studies,
Fig. 4. Mean pressure field for flow over a square rib, ReZ420,000.
572
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
Fig. 5. Mean streamline distribution for flow over a square rib, ReZ420,000.
Fig. 9. Comparison of velocity predictions and measurements on upstream of the square rib.
Fig. 6. Mean vorticity distribution for flow over a square rib, ReZ420,000.
Fig. 7. Variation of drag coefficient with time for flow over a square rib.
Fig. 10. Comparison of velocity predictions and measurements on downstream of the square rib.
this unsteady phenomenon also reflects on the coefficients of drag and lift as a function of time. 5.2. Turbulent flow around a cube
Fig. 8. Variation of lift coefficient with time for flow over a square rib.
The large eddy simulation by three-step FEM–BEM model has been applied for the simulation to the case of flow around a surface mounted cube. The unsteady flow field around a cube is very complex and much more difficult to predict than the above described two-dimensional case. In addition to the vortex structure in the two-dimensional case, there are more complex vortex structures, consisting of four major vortex systems [9]. They are: (1) the secondary currents generated by the flow
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
573
Fig. 14. Velocity field for flow around a cube on x–z plane at yZ0, tZ160.
Fig. 11. Comparison of velocity predictions and measurements on top of the square rib. Fig. 15. Pressure field for flow around a cube on x–z plane at yZ0, tZ160.
Fig. 12. Mean wall static pressure distribution for flow over a square rib. Fig. 16. Vorticity distribution for flow around a cube on x–z plane at yZ0, tZ160.
Fig. 13. Computational domain for flow around a cube with discretization.
Fig. 17. Pressure field for flow around a cube on y–z plane at xZ0, tZ160.
574
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
Fig. 18. Vorticity distribution for flow around a cube on y–z plane at xZ0, tZ160.
separation on the top of the cube and the large circulation cavity behind the cube; (2) von Karman vortex street with axes in the vertical direction; (3) Horseshoe vortex around the structure near the ground; and (4) the twin axial vortices issuing from the side edges of the top of cube. Moreover, the three-dimensional flow always appears to be at transient state. The form of calculation domain with mesh discretization is shown in Fig. 13. As in Fig. 13 (by trial and error) the inlet boundary is located sufficiently upstream of the obstacle so that
its presence is not felt and the outlet is located sufficiently downstream of the point of reattachment. The computational domain in the present study corresponds to x/dZ22.5, y/dZ8 and z/dZ7 (d is the size of the cube). No-slip boundary conditions are imposed along the walls, and uZu0, vZ0, wZ0 are prescribed along the inlet boundary. At the inlet plane, a power-law boundary layer profile for a smooth wall is imposed on axial velocity as in the case of Eq. (51), in which d is the boundary layer thickness taken as half the thickness of the cube. In the present numerical analysis, the turbulent flow is considered at a Reynolds number of 420,000. For the analysis, the domain is discretized into a mesh of 61!37!29, using 58,752 tri-linear cubic elements and 64,001 nodes, as in Fig. 13. Here a fine mesh is used to capture the various details of turbulence. A time step of 0.05 is used in the analysis. As mentioned earlier, the turbulent flow over a cube contains four types of vortices. The time averaged flow patterns from the present numerical investigations are shown in Fig. 1, which are almost same as the earlier numerical as well as experimental results [7,9,27]. As in Fig. 1, the shear layer develops into a large recirculation cavity region 1. Along the obstacle height, the vertical vortices 2 induced by von Karman vortex shedding occur. The horseshoe-shaped vortices 3 occur around the base of the cube and the tip vortex 4 occurs away downstream. Figs. 14–16 represent the velocity vectors, pressure field and vorticity distribution, respectively, at the instant of time tZ160
Fig. 19. Velocity field for flow around a cube on x–y plane at zZ0.5, tZ160.
Fig. 20. Vorticity distribution for flow around a cube on x–y plane at zZ0.5, tZ160.
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
(non-dimensional) at yZ0.0 in the x–z plane. The pressure field and vorticity distributions at the instant time tZ160 (nondimensional) at xZ0.0 in the y–z plane, are represented by Figs. 17 and 18, respectively. Figs. 19 and 20 represent, respectively, the velocity field and vorticity distribution at the instant of time tZ160 at zZ0.5 in the x–y plane. When compared the velocity field to the two-dimensional case, the reattachment length behind the obstacle in Fig. 14 is much shorter than that in Fig. 3 and the separation is much more significant. For the three-dimensional cube case, the front separation is inclined as the flow passes around both sides of the cube, forming the horseshoe vortex structure. It can also be seen that the flow, separating after the front edge, almost no longer reattaches on the top of the cube. Although we do not present a quantitative comparison of various parameters with other LES models, a qualitative comparison of flow and pressure fields with the FVM or FDM results of He and Song [28], Yu and Kareem [29], and Murakami et al. [7] all gives encouraging overall agreement considering the coarse nature of the mesh used. Further observation of the computational results demonstrates the complex characteristics of unsteady and threedimensional nature of the LES turbulence model. The flow structures vary temporally as well as locally. The turbulent features are vividly simulated in the present study, especially for the dynamics of vortex shedding, vortex stretching, and interaction between tip, Karman and horseshoe vortices. These properties all relate to the mechanism of flow separation. They are very difficult to numerical modeling or experimental measurement, due to the unsteady and three-dimensional nature. 6. Concluding remarks A novel computational model has been developed to solve two- and three-dimensional turbulent flow problems in external flow field with large eddy simulation (LES) concept by coupling three-step FEM and BEM. The pressure projection model based on LES form of the Navier–Stokes equations is able to solve the infinite boundary value problems by extracting the boundary effects on a specified finite computational domain. In the proposed model, the Navier–Stokes equations in the LES form are solved using a three-step FEM, and the Poisson type pressure equations are solved using BEM. By coupling the three-step FEM and BEM, the model is able to handle infinite domain problems. The model has been applied to solve a twodimensional turbulent flow over a square rib, and a threedimensional turbulent flow around a surface mounted cube. Both applications give reasonable and comparable results and further show the feasibility of the proposed model as comparing with other models in the literature even for smaller number of grids. Acknowledgements The work reported in this paper was supported by the National Science Council of Taiwan under the Grant no. NSC 93-E-2611-002-001. It is greatly appreciated.
575
References [1] Blevins RD. Flow-induced vibration. 2nd ed. New York: Van Nostrand Reinhold; 1990. [2] Baetke F, Werner H. Numerical simulation of turbulence flow over surface-mounted obstacles with sharp edges and corners. J Wind Eng Ind Aerod 1990;35:129–47. [3] Smagorinsky JS. General circulation experiments with the primitive equations. Part 1. Basic experiments. Mon Weather Rev 1963;91:99–164. [4] Deardorff JW. A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers. J Fluid Mech 1970;41(2):453–80. [5] Schumann U. Sub-grid scale model for finite difference simulation of turbulent flows in plane channels and annuli. J Comput Phys 1975;18: 376–404. [6] Moin P, Kim J. Numerical investigation of turbulent channel flow. J Fluid Mech 1982;118:341–77. [7] Murakami S, Mochida A, Hibi K. Three-dimensional simulation of air flow around a cubic model by means of large eddy simulation. J Wind Eng Ind Aerod 1987;25:291–305. [8] Kobayashi T, Ishiara T, Kano M. Prediction of turbulent flow in twodimensional channel with turbulence promoters, III— improvement of large eddy simulation and formation of streaklines. JSME Bull 1985; 28(246):2948–53. [9] Song CCS, He J. Computation of wind flow around a tall building and the large-scale vortex structure. J Wind Eng Ind Aerod 1993;46–47:219–28. [10] Rail MM, Moin P. Direct simulations of turbulent flow using finite difference schemes. J Comput Phys 1991;96:15–53. [11] Reddy JN, Gartling DK. The finite element method in heat transfer and fluid dynamics. London: CRC Press; 1992. [12] Power H, Wrobel LC. Boundary integral methods in fluid mechanics. Southampton: Computational Mechanics Publications; 1995. [13] Brookes AN, Hughes TJR. Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations. Comput Meth Appl Mech Eng 1982;32:199–259. [14] Lohner R, Morgan K, Zienkiewics OC. The solution of nonlinear hyperbolic equations systems by the finite element method. Int J Numer Methods Fluids 1984;4:1043–63. [15] Jiang CB, Kawahara M. The analysis of unsteady incompressible flows by a three-step finite element method. Int J Numer Methods Fluids 1993;16: 793–811. [16] Young DL, Chang JT, Eldho TI. A coupled BEM and arbitrary Lagrangian–Eulerian FEM model for the solution of two-dimensional laminar flows in external flow fields. Int J Numer Meth Eng 2001;(9): 1053–77. [17] Young DL, Huang JL, Eldho TI. Simulation of laminar vortex shedding flow past cylinders using a coupled BEM and FEM model. Comput Meth Appl Mech Eng 2001;190:5975–98. [18] Young DL, Huang JL, Eldho TI. Numerical simulation of high-Reynolds number flow around circular cylinders by a three-step FEM–BEM model. Int J Numer Methods Fluids 2001;37:657–89. [19] Young DL, Chang JT, Eldho TI. Solution of three-dimensional unsteady external flow using a coupled arbitrary Lagrangian FEM–BEM model. Eng Anal Bound Elem 2004;28:711–23. [20] Yoshizawa A. Eddy-viscosity-type subgrid-scale model with a variable Smagorinsky coefficient and its relationship with the oneequation model in large eddy simulation. J Phys Fluids 1991;A3: 2007–9. [21] Debler WA. Fluid mechanics fundamentals. Englewood Cliffs, NJ: Prentice Hall; 1990. [22] Reddy JN. Finite element method. New York: McGraw-Hill; 1993. [23] Young DL, Liao CB, Sheen HJ. Computations of recirculation zones of a confined annular swirling flow. Int J Numer Methods Fluids 1999;29: 791–810. [24] Brebbia CA, Telles JCF, Wrobel LC. Boundary element techniques— theory and applications in engineering. Berlin: Springer; 1984.
576
D.L. Young et al. / Engineering Analysis with Boundary Elements 30 (2006) 564–576
[25] Crabb D, Durao DFG, Whitelaw JH. Velocity characteristics in the vicinity of a two dimensional rib Proceedings of fourth Brazilian congress on mechanical engineering, Florianopolis, Brazil, December; 1977. p. 415–29. [26] Benodekar RW, Goddard AJH, Gosman AD, Issac RI. Numerical prediction of turbulent flow over surface-mounted ribs. AIAA J 1985; 23(3):359–66.
[27] Castro IP, Robins AG. The flow around a surface mounted cube in uniform and turbulent streams. J Fluid Mech 1977;79(2):307–35. [28] He J, Song CCS. Computation of turbulent shear flow over surface-mounted obstacle. J Eng Mech ASCE 1992;118(11): 2282–97. [29] Yu D, Kareem A. Parametric study of flow around rectangular prisms using LES. J Wind Eng Ind Aerod 1998;77/78:653–62.