V. D. Pereira and J. B. Campos Silva / Vol. XXVII, No. 3, July-September 2005 ABCM 274 V. D. Pereira and J. B. Campos Silva UNESP Ilha Solteira Faculdade de Engenharia Departamento de Engenharia Mecânica Av. Brasil Sul, 56 15385-000 Ilha Solteira, SP. Brazil vdvanco@yahoo.com.br jbcampos@dem.feis.unesp.br Simulations of Incompressible Fluid Flows by a Least Squares Finite Element Method In this work simulations of incompressible fluid flows have been done by a Least Squares Finite Element Method (LSFEM) using velocity-pressure-vorticity and velocity-pressure- stress formulations, named u-p-ω and u-p-τ formulations respectively. These formulations are preferred because the resulting equations are partial differential equations of first order, which is convenient for implementation by LSFEM. The main purposes of this work are the numerical computation of laminar, transitional and turbulent fluid flows through the application of large eddy simulation (LES) methodology using the LSFEM. The Navier- Stokes equations in u-p-ω and u-p-τ formulations are filtered and the eddy viscosity model of Smagorinsky is used for modeling the sub-grid-scale stresses. Some benchmark problems are solved for validate the numerical code and the preliminary results are presented and compared with available results from the literature. Keywords: Navier-Stokes equations, large eddy simulation, least-squares finite element, fluid flows Introduction 1The finite element method (FEM) is one of the most used techniques for numerical solution of partial differential equations in engineering and applied sciences. The mathematical foundation of the finite element method can be based on the weight residual method (WRM), Finlayson, (1972), which originate different formulations according to the weight functions used. The main versions of the FEM are the Bubnov-Galerkin, Petrov-Galerkin, Collocation, Sub-domain and Least-Squares formulations. Another classification underlining the variational principle considers three major groups: the Rayleigh-Ritz method, The Galerkin Method and the Least-squares method. For convection dominated problems the Galerkin-based methods present often spurious oscillation of the solutions (Jiang, 1998). In recent works, Romão et al. (2003) and Romão (2004) applied different versions of the finite element method for convection-diffusion problems. Several authors have investigated the LSFEM for solution of incompressible and compressible fluid flows. Jiang (1998) presented a list of such works. Winterscheidt & Surana (1994) also have applied p-versions of least-squares finite element method for fluid flows. The least squares have also been used for stabilization of the Galerkin finite element method. Jansen (1999) presented a work for computing turbulent flows in unstructured-grids. In the present work, the least- squares finite element method (LSFEM), Jiang (1998), has been implemented for simulation of laminar, transitional and turbulent fluid flows using velocity-pressure-vorticity, u-p-ω and u-p-τ formulations which result in a first-order partial differential system of equations. The filtered Navier-Stokes equations with the Smagorinsky model of eddy viscosity are solved for simulations of the lid-driven cavity flow and the backward-facing step flow as benchmark problems. Which is of known of the authors the large-eddy simulation (LES) was implemented with the velocity- pressure-stress, u-p-τ, formulation of Navier-Stokes equations by Ding & Tsang (2001). We didn’t find in the literature that we had access the u-p-ω formulation with LES for fluid flow simulation in the way presented here. Although, turbulence is a 3D phenomenon, in this preliminary work, only 2D simulations have been considered to understand the behavior of the LSFEM in the proposed formulations and due to the computational capacity available. Paper accepted June, 2005. Technical Editor: Aristeu da Silveira Neto. Jiang (1998) enumerated several features of the LSFEM, among them: universality, efficiency, robustness, optimality, concurrent simulation of multiple physics and general-purpose coding. Jiang also claimed that no upwind technique is necessary for numerical calculation of convection dominated problems, because the resulting matrix systems of equations from the LSFEM application are always symmetrical and positive-definite. For Galerkin type solutions of fluid dynamics, upwind generally has to be applied for stability. In this work, the lid-driven cavity and the backward-facing step flows have been solved with linear and quadratic quadrilateral elements for investigation of different interpolation functions, and some values of the Smagorinsky constant have been also tried. The refine of the meshes has been investigated, too. The results are compared with results from other authors. Beyond of this introduction, the paper covers some aspects of formulation of the proposed model, presents some results, discussions, conclusions and references. Nomenclature k = turbulent kinetic energy L = length of reference p = pressure 2 0u ppP o ρ − = = dimensionless pressure Re = Reynolds number t = Time u = component of dimensional velocity in the ix - axis direction iu = component of velocity in the ix - axis direction 0u = reference velocity ouuU = - dimensionless component of velocity in the X-axis direction iU = dimensionless component of velocity in the iX -axis direction v = component of dimensional velocity in the y -axis direction 0/ uvV = - dimensionless component of velocity in the Y -axis direction LxX = = dimensionless X coordinate ix = ith- axis in Cartesian coordinates Simulations of Incompressible Fluid Flows by a … J. of the Braz. Soc. of Mech. Sci. & Eng. Copyright © 2005 by ABCM July-September 2005, Vol. XXVII, No. 3 / 275 LxX ii = = ith dimensionless coordinate LyY = = dimensionless Y coordinate Greek Symbols α = index that indicates local node number inside an element ijδ = Krönecker delta θ = time discretization parameter µ = dynamic viscosity tµ = dynamic eddy viscosity ρ = density φ = any scalar or vector variable Φ = nodal variable in elements ψ = stream function jω = vorticity around the j-axis Superscripts n = variable evaluated at time t n+1 = variable evaluated at time tt ∆+ Subscripts i = direction of the axis in the system of coordinates or component j = direction of the axis in the system of coordinates or component k = direction of the axis in the system of coordinates or component Formulation Governing Equations The Navier-Stokes equations for incompressible transient fluid flows in vector notation can be written as follow: fu uuu =∇++⎟ ⎠ ⎞ ⎜ ⎝ ⎛ •+ ∂ ∂ 2µρ p t ∇∇ (1) 0=•u∇ (2) where ρ is the fluid density, u is the velocity vector with components iu , p is the pressure, µ is the dynamic viscosity and f is the body forces vector with components if . The Equation (1) is a second order partial differential equation and this is not the most appropriated form for solution by LSFEM. The LSFEM generally is applied for first order differential equations. However, second order partial differential can be transformed in first order system by using auxiliary variables. This is another advantage of the least-squares method: the direct calculation of secondary variables that have physical interpretation such as heat and mass fluxes, stresses and vorticity. According to Brodkey (1967), using vectorial identities: uuuuuu ××−•=• ∇∇∇ 2/)( and )()(2 uuu ××−•=∇ ∇∇∇∇ , the Navier-Stokes can be rewritten as, now in tensorial notation i j k ijk i i kjijk i f xx up u t u = ∂ ∂ + ∂ +∂ +⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − ∂ ∂ ω µεωερ )2/( 2 ; (3) 0= ∂ ∂ i i x u ; (4) j k ijki x u ∂ ∂ ε=ω . (5) For application of the large-eddy simulation methodology, the equations must be filtered for separation of the large scales from the sub-grid scales. So, the large scales are simulated by solution of the equations for the filtered variables after modeling the sub-grid scales terms that come from the filter process. Chidambaram (1998) presented different filter functions for LES. The filtered Eq. (3)–(5) are of the form ( ) ikjkjijk j k ijk i ii i i kjijk i fuu x x ]/)uu(p[ x uu t u =ω−ωρε− ∂ ω∂ µε+ + ∂ −ρ+∂ +⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ +ωε− ∂ ∂ ρ 2 2 1 222 (6) 0= ∂ ∂ i i x u ; (7) j k ijki x u ∂ ∂ = εω . (8) The differences between Eqs. (3) and (6) are the additional term to the pressure and the fourth term of left hand side of Eq. (6) that originated from the convection term of the Navier-Stokes equations. These terms correspond to the turbulent kinetic energy and the vorticity of the sub-grid scales respectively. The purpose of this work is the modeling of the fourth term, by analogy with the modeling of the sub-grid scale stresses in the conventional formulation of the Navier-Stokes equations. So, it is defined the sub-grid scale effects and the turbulent pressure as j k tkjkj x uu ∂ ∂ −=− ω µωωρ )( ; 2/)( 22 iit uupp −+= ρ . (9) Now, after modeling of the sub-grid scale effects the dimensionless form of the Eqs. (6)-(8) is as follow i j k tijk i t j i j i S XReX P X UU t U = ∂ Ω∂ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ν+ε+ ∂ ∂ + ∂ ∂ + ∂ ∂ 1 ; (10) 0= ∂ ∂ i i X U ; (11) j k ijki X U ∂ ∂ ε=Ω . (12) The dimensionless variables in Eqs. (10)-(12) are defined in function of the characteristic parameters of length L and velocity 0u as L x X i i = ; 0u u U i i = ; 2 0 0 u pp P t t ρ − = ; 0 * / uL tt = ; ( ) 2/1 2 0 * 2 klkls t t SS L C Lu ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∆ == ν ν ; 0u Li i ω =Ω ; µ ρ Lu0Re = . The eddy viscosity is calculated according to the Smagorinsky model in the form V. D. Pereira and J. B. Campos Silva / Vol. XXVII, No. 3, July-September 2005 ABCM 276 ( ) ( ) 212 2 / klkls * t ssC ∆=ν ; ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ = i j j i ij x u x us 2 1 (13) where sC is the constant of Smagorinsky and ∆ is the filter width defined as: ( ) 3/1zyx ∆∆∆=∆ for 3D or ( ) 2/1yx∆∆=∆ for 2D geometry. First Order System and Interpolation of Variables The first order system (10)-(12) for 2D problems, after discretizing the transient term can be written in a compact form as nn fL =φ +1 ; (14) where [ ]T,P,V,U Ω=φ is the vector of unknown variables, [ ]Tvu ,,S,Sf 00= is the vector of independent terms and now L is a matrix differential operator defined as ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ + ∆ ∂ ∂ ∂ ∂ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ + ∆ = 1 0 - 0 0 Y - Y 1 0 0 1 e e XY X Xy V x U YXY V X U L νθ τ νθ τ (15) The source terms uS and vS are: ( ) n e n Xu YX P Y UV X UU t UfS ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ Ω∂ µ+ ∂ ∂ + ∂ ∂ + ∂ ∂ θ−− ∆ += 1 ; (16) ( ) n e n Yv XY P Y VV X VU t VfS ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ Ω∂ − ∂ ∂ + ∂ ∂ + ∂ ∂ −− ∆ += µθ1 . (17) and 10 ≤≤ θ is a parameter of time discretization. The variable φ in finite element methods, for equal order interpolation of all degrees of freedom, can be interpolated inside an element in the form: ∑ = ⎪ ⎪ ⎭ ⎪ ⎪ ⎬ ⎫ ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ Ω =φ Ne j j j e )t( P V U )Y,X(N)t,Y,X( 1 ; (18) where jN is the interpolation function associated to the node j of the element and Ne is the number of nodes. It has been pointed out that the LSFEM is not subjected to LBB (Ladyzhenskaya-Babuska- Brezzi) condition like the Galerkin-based method, Jiang (1998), Winterscheidt & Surana (1994). Least-Squares Finite Element Method Substituting Eq. (18) in Eq. (14) a residual vector can be defined inside an element as nn fLNR −Φ= +1 . (19) The application of LSFEM consists in the minimization of a functional defined as the integral of the squared residuals. If inside an element ones define a functional as ( ) ∫ •=Φ + eA Tn e RdARJ 1 , the functional, in the whole domain divided in Nelem elements, can be calculated as follow ( ) ∑ ∫ = + •=Φ Nelem e eA Tn RdARJ 1 1 . (20) The minimization of the functional requires that ( ) 01 =Φ +nJδ , which results in the following matrix system: FK =Φ (21) Now, in Eq. (21), Φ is the global vector of nodal values. The global matrix K is assembled from the element matrices, dALNLNLNLNLNLNK Ne T A Nee e ,,,,,, 2121∫= (22) where eA is the area of the finite element, T denotes the transpose and the global vector F is assembled with the contribution of the element vectors fdALNLNLNF T A Nee e ∫= ,,, 21 ; (23) in which ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ + ∆ ∂ ∂ ∂ ∂ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ + ∆ = j jj jj jjjjj jjjjj j N X N Y N N X N X NN y N V x N U N Y N X N Y N V X N U N LN 0 - 0 0 Y - Y 0 0 e e µθ τ µθ τ (24) Results and Discussions for Velocity-Pressure-Vorticity Formulation Lid-Driven Cavity flow In order to validate the proposed model, the classical benchmark problem of the lid-driven cavity flow was numerically simulated. The convective terms are linearized by successive substitutions in each time step and the linear system resulting is solved by the frontal method (Taylor & Hughes, 1981). However, a solver like the preconditioned conjugate gradient method should be preferred, once the matrix from the LSFEM is always symmetric and positive- definite. The element matrices are obtained using three-point for linear quadrilateral elements and nine-point Gaussian quadrature for the nine-node finite element. The Reynolds numbers of 100, 400 1,000 and 10,000 were considered. It has been tested some values for the constant of Smagorinsky and the mesh refinement. The Figure 1 shows the geometry and boundary conditions. With the aid of secondary variables all boundary conditions are of Dirichlet type Simulations of Incompressible Fluid Flows by a … J. of the Braz. Soc. of Mech. Sci. & Eng. Copyright © 2005 by ABCM July-September 2005, Vol. XXVII, No. 3 / 277 and in most cases no boundary condition is necessary for the vorticity, Jiang (1992). Figures 2 to 5 show the results for U velocity component at the mid of the cavity for different low Reynolds numbers and a mesh of 100 x 100 bilinear elements and mesh of 30 x 30 quadratic elements, without turbulence model. The four-node elements give poor results and the nine-node element give better results, but no so good yet. Most probably, the refinement of the meshes could lead to results in better agreement with results from literature. (1,0) (0,1) (1,1) Figure 1. Cavity geometry and boundary conditions. -0,4 -0,2 0,0 0,2 0,4 0,6 0,8 1,0 1,2 0,0 0,2 0,4 0,6 0,8 1,0 Re = 400 (present work) Re = 400 (Ghia et al., 1982) Y U Figure 2. Velocity U at X = 0.5, linear element. -0,4 -0,2 0,0 0,2 0,4 0,6 0,8 1,0 1,2 0,0 0,2 0,4 0,6 0,8 1,0 Re = 400 ( present work) Re = 400 (Ghia et al.,1982) Y U Figure 3. Velocity U at X = 0.5, quadratic element. -0,4 -0,2 0,0 0,2 0,4 0,6 0,8 1,0 1,2 0,0 0,2 0,4 0,6 0,8 1,0 Re = 1000 (present work) Re = 1000 (Ghia et al., 1982) Y U Figure 4. Velocity U at X = 0.5, linear element. -0,4 -0,2 0,0 0,2 0,4 0,6 0,8 1,0 1,2 0,0 0,2 0,4 0,6 0,8 1,0 Re = 1000 (present work) Re = 1000 (Ghia et al., 1982) Y U Figure 5. Velocity U at X = 0.5, quadratic element. V. D. Pereira and J. B. Campos Silva / Vol. XXVII, No. 3, July-September 2005 ABCM 278 Figure 6 shows the component U of velocity for different Reynolds numbers for a mesh of 30 x 30 quadratic element and Cs = 0,1612. When the Reynolds increases the results are not in so good agreement with the results of Ghia et al. (1982). Figure 7 shows the component V for the same cases of Figure 6 and a similar behavior of the velocity field is observed. The Figures 8 and 9 show that when the mesh is refined, the agreement between the results is quite good. The mesh for the results of Figures 6 and 7 is of 60 by 60 elements or 121 by 121 grid points. For a high Reynolds, Figure 10 shows yet some discrepancy between the results. This difference may due to the constant of Smagorinsky adopted and a non- sufficient refinement of the mesh. We are investigating these possibilities. Results for the stream function are shown in Figures 11 to 14 for Reynolds numbers of 100, 400 1000 and 10,000 and turbulence model in the simulations. The lengths of the secondary vortexes are in agreement to the results from the literature. -0,4 -0,2 0,0 0,2 0,4 0,6 0,8 1,0 1,2 0,0 0,2 0,4 0,6 0,8 1,0 Re = 100 (present work) Re = 400 (present work) Re = 1000 (present work) Re = 100 (Ghia et al., 1982) Re = 400 (Ghia et al., 1982) Re = 1000 (Ghia et al., 1982) Y U Figure 6. U velocity at X = 0.5 for different Reynolds with turbulence model. 0,0 0,2 0,4 0,6 0,8 1,0 -0,7 -0,6 -0,5 -0,4 -0,3 -0,2 -0,1 0,0 0,1 0,2 0,3 0,4 0,5 0,6 X Re = 100 (present work) Re = 400 (present work) Re = 1000 (present work) Re = 100 (Ghia et al., 1982) Re = 400 (Ghia et al., 1982) Re = 1000 (Ghia et al., 1982) V Figure 7. V velocity at Y = 0.5 for different Reynolds with turbulence model. -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.0 0.2 0.4 0.6 0.8 1.0 Re = 400 (present work) Re = 400 (Ghia et al., 1982) Y U Figure 8. U velocity at X = 0.5 for Reynolds 400 with turbulence model and a more refined mesh. -0,4 -0,2 0,0 0,2 0,4 0,6 0,8 1,0 1,2 0,0 0,2 0,4 0,6 0,8 1,0 Re = 1000 (present work) Re = 1000 (Ghia et al., 1982) Y U Figure 9. U velocity at X = 0.5 for Reynolds 1000 with turbulence model and a more refined mesh. -0,6 -0,4 -0,2 0,0 0,2 0,4 0,6 0,8 1,0 1,2 0,0 0,2 0,4 0,6 0,8 1,0 Re = 10000 (present work) Re = 10000 (Ghia et al., 1982) Y U Figure 10. U velocity at X = 0.5 for Reynolds 10,000 with turbulence model. Simulations of Incompressible Fluid Flows by a … J. of the Braz. Soc. of Mech. Sci. & Eng. Copyright © 2005 by ABCM July-September 2005, Vol. XXVII, No. 3 / 279 1 2 3 3 4 4 4 5 5 5 6 6 6 6 7 7 7 7 8 8 8 8 8 9 9 9 9 9 9 10 10 10 10 10 10 11 11 11 11 11 11 11 12 12 12 12 12 12 12 13 13 13 13 13 13 1313 14 14 14 14 14 14 1414 X Y 0 0.25 0.5 0.75 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Level SF 14 -1.33832E-06 13 -1.29465E-05 12 -0.000144182 11 -0.00163398 10 -0.0065731 9 -0.0141146 8 -0.0253852 7 -0.0380777 6 -0.0507703 5 -0.0634629 4 -0.0761555 3 -0.0888481 2 -0.0951944 1 -0.0998343 Figure 11. Stream function for Re = 100 with turbulence model. 2 3 3 4 4 5 5 6 6 6 7 7 7 7 8 8 8 8 9 9 9 9 9 10 10 10 10 10 10 11 11 11 11 11 11 11 12 12 12 12 12 12 12 12 13 13 13 13 13 13 13 13 14 14 14 14 14 14 14 14 15 15 15 15 15 15 15 15 15 16 16 17 17 18 19 X Y 0 0.25 0.5 0.75 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Level SF 19 0.000391569 18 0.000217252 17 3.88613E-05 16 2.39935E-06 15 -6.92722E-08 14 -1.03687E-05 13 -8.01651E-05 12 -0.00128279 11 -0.00603213 10 -0.0189839 9 -0.0319357 8 -0.0448875 7 -0.0578393 6 -0.0707911 5 -0.0837429 4 -0.0902188 3 -0.0966947 2 -0.099686 1 -0.102323 Figure 12. Stream function for Re = 400 with turbulence model. 2 3 4 4 5 5 56 6 6 7 7 7 7 8 8 8 8 9 9 9 9 9 10 10 10 10 10 10 11 11 11 11 11 11 12 12 12 12 12 12 1213 13 13 13 13 13 13 13 14 14 14 14 14 14 14 14 15 15 15 15 16 16 16 17 17 18 19 20 21 X Y 0 0.25 0.5 0.75 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Level SF 21 0.0012879 20 0.00125246 19 0.000934794 18 0.00064222 17 0.000285514 16 5.58351E-05 15 1.48062E-06 14 -4.82521E-05 13 -0.000620296 12 -0.00497407 11 -0.0175007 10 -0.0300273 9 -0.0425539 8 -0.0550805 7 -0.0676071 6 -0.0801337 5 -0.086397 4 -0.0926603 3 -0.0966589 2 -0.0980666 1 -0.0987134 Figure 13. Stream function for Re = 1,000 with turbulence model. 3 3 4 4 5 5 5 6 6 6 6 7 7 7 7 8 8 8 8 8 9 9 9 9 9 10 10 10 10 10 10 1011 11 11 11 11 11 11 12 12 12 12 12 12 12 12 13 13 13 14 14 14 15 X Y 0 0.25 0.5 0.75 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Level SF 16 0.00369097 15 0.00233594 14 0.000713377 13 2.19354E-05 12 -0.000217243 11 -0.00609176 10 -0.016061 9 -0.0359996 8 -0.0559382 7 -0.0758768 6 -0.0958154 5 -0.115754 4 -0.135693 3 -0.145662 2 -0.153399 1 -0.1555 Figure 14. Stream Function for Re = 10,000 with turbulence model. Formulation u-p-τ Now it is considered the u-p-τ formulation. Applying a filtering process to the Navier-Stokes equation (2.1) result the following equations: i j sgs ij j L ij ij jii f xxx p x uu t u + ∂ ∂ − ∂ ∂ + ∂ ∂ −= ∂ ∂ + ∂ ∂ ττ ρ 1 ; (25) 0= ∂ ∂ i i x u . (26) where the viscous and the sub-grid-scale stresses are defined as ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ = i j j iL ij x u x uντ ; (27) ijij sgs ij uuuu −=τ . (28) The sub-grid-scale stress tensor is modeled through the model of Smagorinsky resulting ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ −= i j j i t sgs ij x u x uk ντ 3 2 . (29) The modeled Navier-Stokes equations, in u-p-τ formulation, now have the form ( ) i j ij ij i j i f xx kp x u u t u + ∂ ∂ + ∂ +∂ −= ∂ ∂ + ∂ ∂ τρ 3/2/ ; (30) ( ) ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ += i j j i tij x u x u νντ ; (31) 0= ∂ ∂ i i x u . (32) V. D. Pereira and J. B. Campos Silva / Vol. XXVII, No. 3, July-September 2005 ABCM 280 In Cartesian coordinates results the following system of equations for 2D problems: x xyxx f yxx p y uv x uu t u =⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ − ∂ ∂ +⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ + ∂ ∂ ττρρ ; (33) y yyxy f yxy p y vv x vu t v =⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ − ∂ ∂ +⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ + ∂ ∂ ττ ρρ ; (34) 02 = ∂ ∂ − x u exx µτ ; (35) 0=⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ − x v y u exy µτ ; (36) 02 = ∂ ∂ − y v eyy µτ ; (37) 0= ∂ ∂ + ∂ ∂ y v x u . (38) In the Equations (33)-(38) vu , are the components of velocity in the axis yx, , respectively; p is pressure field ijτ is the stress tensor and te µµµ += is dynamic effective viscosity. Results and Discussions for Velocity-Pressure-Stress Formulation Lid-Driven Cavity flow In this section are presented some results from tests of the u-p-τ formulation. Initially, were tried quadrilateral elements to discretized the domain. The results don’t converged to the expected solution, mainly with the increase of the Reynolds number. The results presented in the following, Fig. 17-18, were obtained using nine-node elements without turbulence model, the aim here was test this kind of formulation that can useful for simulation of Non- Newtonian fluid flows. For the Reynolds numbers presented the agreement between the results of the present work and those ones from Ghia et al. (1982) is quite good. The mesh here was of 121 by 121 grid points while Ghia et al. used 129 by 129 grid points. -0,4 -0,2 0,0 0,2 0,4 0,6 0,8 1,0 0,0 0,2 0,4 0,6 0,8 1,0 Re = 100 (present work) Re = 100 (Ghia et al., 1982) Y U Figure 15. Velocity U at X = 0.5; linear element. 0,0 0,2 0,4 0,6 0,8 1,0 -0,3 -0,2 -0,1 0,0 0,1 0,2 Re = 100 (present work) Re = 100 (Ghia et al., 1982) V X Figure 16. Velocity V at Y = 0.5; linear element. -0,4 0,0 0,4 0,8 1,2 0,0 0,2 0,4 0,6 0,8 1,0 Re = 400 (present work) Re = 400 (Ghia et al., 1982) Y U Figure 17. Velocity U at X = 0.5; linear element. 0,0 0,2 0,4 0,6 0,8 1,0 -0,5 -0,4 -0,3 -0,2 -0,1 0,0 0,1 0,2 0,3 0,4 Re = 400 (present work) Re = 400 (Ghia et al., 1982) V X Figure 18. Velocity V at Y = 0.5; linear element. Simulations of Incompressible Fluid Flows by a … J. of the Braz. Soc. of Mech. Sci. & Eng. Copyright © 2005 by ABCM July-September 2005, Vol. XXVII, No. 3 / 281 Backward-Facing Step Flow In this section some results for a laminar backward-facing step flow are presented also using the u-p-τ formulation for Reynolds of 100, 200 and 400. The Reynolds numbers has been based on the step height. The ratio of expansion in this case is of 1:2. At the entrance of the channel was imposed a parabolic velocity profile of the form: [ ]22 0 /)(/)(22/3 ww rhyrhyuu −−−= , where h is the step height, wr is the half spacing of the small channel and in this case 2/hrw = . Figures 19-21 show profiles of velocity at some stations along the channel and Figures 22-24 shows the streamlines for Re of 100, 200 and 400 respectively. The behavior of the flow has been simulated satisfactorily. The non-dimensional reattachment lengths are approximately, 4.2, 6.2 and 8 times step height for Re = 100, 200, 400 respectively. The length reattachment results of Barber & Fonty (2003) for a similar flow and the same Reynolds numbers are of about 3.5, 5.5 and 8.4 times step height. For the Reynolds number of 400 a second vortex next the upper wall appears. Some authors say that for Reynolds greater than 400 three-dimensional effects appear in the flow. -0,2 0,0 0,2 0,4 0,6 0,8 1,0 1,2 1,4 1,6 0,0 0,5 1,0 1,5 2,0 Re 100 Re = 0,74 Re = 1,88 Re = 2,65 Re = 10,69 Re = 28,74 Y U Figure 19. Profiles of velocity at some stations along the channel. -0,2 0,0 0,2 0,4 0,6 0,8 1,0 1,2 1,4 1,6 0,0 0,5 1,0 1,5 2,0 Re = 200 X = 0,74 X = 1,88 X = 2,65 X = 10,69 X = 28,74 Y U Figure 20. Profiles of velocity at some stations along the channel. -0,2 0,0 0,2 0,4 0,6 0,8 1,0 1,2 1,4 1,6 0,0 0,5 1,0 1,5 2,0 Re = 400 X = 0,74 X = 1,88 X = 2,65 X = 10,69 X = 28,74 Y U Figure 21. Profiles of velocity at some stations along the channel. X Y 0 10 20 30 40 0 0.5 1 1.5 2 Figure 22. Streamlines for Re = 100. X Y 0 10 20 30 40 0 0.25 0.5 0.75 1 1.25 1.5 1.75 2 Figure 23. Streamlines for Re = 200. X Y 0 10 20 30 40 0 0.5 1 1.5 2 Figure 24. Streamlines for Re = 400. Conclusions A least-squares finite element method with eddy viscosity model of Smagorinsky has been implemented in this work for simulation of Navier-Stokes equations, in u-p-ω formulation. The four-node elements with three-point Gaussian quadrature not lead to good V. D. Pereira and J. B. Campos Silva / Vol. XXVII, No. 3, July-September 2005 ABCM 282 results. The nine-node elements with nine-point Gaussian quadrature presented better results with the refinement of the mesh. For low Reynolds number even a coarse mesh produces good agreement. Since, the interest is to simulate high Reynolds flows, more investigation is still necessary for improvement of the model. Some points to be investigated are the constant of Smagorinsky and filter width. Some aspects of the method were elucidated, but several enhancements can still be done. In the LSFEM the interpolation functions and the rules of integration must be choose with care, this is another point of investigation. According to Jiang (1992) troubles with interpolation functions can be overcome with the use of the p-version least-squares method. Since turbulence is a three-dimensional phenomenon, cases of 3D geometry shall be treated in future works. These aspects pointed as possible causes of discrepancies between the results are nowadays being investigated. Acknowledgements The authors would like to acknowledge the Fundação de Amparo a Pesquisa do Estado de São Paulo = FAPESP for the computational resources (Proc. 03/11737-2) and the CAPES for the financial support. References Barber, R. W. & Fonty, A., 2003. Comparison of vortex-element and finite-volume simulations of low Reynolds number flow over a confined backward-facing step. CFD2003 – Vortex Methods, May, 28-30, Vancouver, Canada. Brodkey, R. S., 1967. The Phenomena of Fluid Motions. Dover Publications, Inc. Chidambaram, N., 1998. Colocated-grid Finite Volume Formulation for the Large Eddy Simulation of Incompressible and Compressible Turbulent Flows. M.Sc. Thesis. Graduate College, Department of Mechanical Engineering, Iowa State University, Ames, Iowa, USA. Ding, X. & Tsang, T. H., 2001. Large eddy simulation of turbulent flows by a least-squares finite element method. International Journal for Numerical Methods in Fluids, vol. 37, pp. 297–319. Finlayson, B. A.,1972. The Method of Weighted Residuals and Variational Principles. Academic Press. Ghia, U., Ghia, K.N. & Shin, C.T., 1982. High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method. Journal of Computational Physics, v. 5, n. 48, p.387-411. Jansen, K. E., 1999. A stabilized finite element method for computing turbulence. Comput. Methods Appl. Mech. Engrg, v. 174, p.299-317. Jiang, B-N., 1992. A least squares finite element method for incompressible Navier-Stokes problems. International Journal for Numerical Methods in Fluids, vol. 14, pp. 843–859. Jiang, B-N., 1998. The Least-Squares Finite Element Method: Theory and Applications in Computational Fluid Dynamics and Electromagnetics. Springer-Verlag. Romão, E.C., Campos Silva, J. B. & Aparecido, J.B., 2003. Variantes do Método de Elementos Finitos para Solução de Problemas Convectivo- Difusivos Unidimensionais. In Proceedings of 24th Iberain Latin-American Congress on Computational Methods in Engineering – CILAMCE 2003 (paper CIL 198-31), October 29-31, 2003, Ouro Preto, Minas Gerais, Brasil. Romão, E. C., 2004. Variantes do Método dos Elementos Finitos para Solução de Fenômenos Convectivo-Difusivos Bidimensionais. Dissertação de Mestrado. UNESP, Faculdade de Engenharia de Ilha Solteira, Depto. de Engenharia Mecânica. Taylor, C. & Hughes, T.G., 1981. Finite Element Programming of the Navier-Stokes Equations. Pineridge Press Limited, Swansea, U.K. 244 p. Winterscheidt, D. & Surana, K. S., 1994. p-version least squares finite element formulations for two-dimensional, incompressible fluid flow. International Journal for Numerical Methods in Fluids, vol. 18, pp. 43–69.