Simulations of incompressible fluid flows by a least squares finite element method

Abstract

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-omega and u-p-tau 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-omega and u-p-tau 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.

Navier-Stokes equations; large eddy simulation; least-squares finite element; fluid flows

TECHNICAL PAPERS

Simulations of incompressible fluid flows by a least squares finite element method

V. D. Pereira; 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; [e-mail]; [e-mail]

ABSTRACT

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-w and u-p-t 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-w and u-p-t 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

The 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-w and u-p-t 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-t, formulation of Navier-Stokes equations by Ding & Tsang (2001). We didn't find in the literature that we had access the u-p-w 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.

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

= dimensionless pressure

Re = Reynolds number

t = Time

u = component of dimensional velocity in the x_i - axis direction

u_i= component of velocity in the x_i - axis direction

u₀= reference velocity

U = u/u₀ - dimensionless component of velocity in the X-axis direction

U_i = dimensionless component of velocity in the X_i -axis direction

n = component of dimensional velocity in the y -axis direction

V = n/u₀ - dimensionless component of velocity in the Y -axis direction

X = x/L dimensionless X coordinate

x_i = i^th- axis in Cartesian coordinates

X_i = x_i / L = i^th dimensionless coordinate

Y = y / L dimensionless Y coordinate

Greek Symbols

a = index that indicates local node number inside an element

d_ij = Krönecker delta

q = time discretization parameter

µ = dynamic viscosity

µ_t = dynamic eddy viscosity

r = density

f = any scalar or vector variable

F = nodal variable in elements

y = stream function

w_j = vorticity around the j-axis

Superscripts

n = variable evaluated at time t

n+1 = variable evaluated at time t + Dt

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:

where r is the fluid density, u is the velocity vector with components u_i, is the pressure, µ is the dynamic viscosity and f is the body forces vector with components f_i.

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: , the Navier-Stokes can be rewritten as, now in tensorial notation

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

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

Now, after modeling of the sub-grid scale effects the dimensionless form of the Eqs. (6)-(8) is as follow

The dimensionless variables in Eqs. (10)-(12) are defined in function of the characteristic parameters of length L and velocity as u₀

The eddy viscosity is calculated according to the Smagorinsky model in the form

where C_s is the constant of Smagorinsky and D is the filter width defined as: D = (DxDyDz)^1/3 for 3D or D = (DxDy)^1/2 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

where f = [U, V, P, W]^T is the vector of unknown variables, f = [S_u, S_v,0,0]^T is the vector of independent terms and now L is a matrix differential operator defined as

The source terms S_u and S_v are:

and 0 < q < 1 is a parameter of time discretization.

The variable j in finite element methods, for equal order interpolation of all degrees of freedom, can be interpolated inside an element in the form:

where N_j is the interpolation function associated to the node j of the element and N_e 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

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 , the functional, in the whole domain divided in Nelem elements, can be calculated as follow

The minimization of the functional requires that dJ(Fⁿ⁺¹) = 0, which results in the following matrix system:

Now, in Eq. (21), F is the global vector of nodal values. The global matrix K is assembled from the element matrices,

where A_e 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

in which

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 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.

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.

Formulation u-p-t

Now it is considered the u-p-t formulation. Applying a filtering process to the Navier-Stokes equation (2.1) result the following equations:

where the viscous and the sub-grid-scale stresses are defined as

The sub-grid-scale stress tensor is modeled through the model of Smagorinsky resulting

The modeled Navier-Stokes equations, in u-p-t formulation, now have the form

In Cartesian coordinates results the following system of equations for 2D problems:

In the Equations (33)-(38) are the components of velocity in the axis x,y, respectively; is pressure field is the stress tensor and µ_e = µ + µ_t 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-t 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.

Backward-Facing Step Flow

In this section some results for a laminar backward-facing step flow are presented also using the u-p-t 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: , where h is the step height, r_w is the half spacing of the small channel and in this case r_w = h/2. 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.

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-w formulation. The four-node elements with three-point Gaussian quadrature not lead to good 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.

Paper accepted June, 2005.

Technical Editor: Aristeu da Silveira Neto.

• 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. 297319.
• 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. 843859.
• 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 24^th 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. 4369.

Publication Dates

• Publication in this collection
  05 Sept 2005
• Date of issue
  Sept 2005