Revision as of 16:54, 26 April 2016 (view source) Test1 (talk | contribs) (Tag: Visual edit) ← Older edit		Latest revision as of 12:40, 7 October 2016 (view source) Scipedia (talk | contribs)
(3 intermediate revisions by one other user not shown)
Line 1:		Line 1:
−	==FIC/FEM FORMULATION WITH MATRIX STABILIZING TERMS FOR INCOMPRESSIBLE FLOWS AT LOW AND HIGH REYNOLDS NUMBERS==
−
−	'''E. Oñate, A. Valls '''and'''  J. García'''
−
−	{|
−	|-
−	International Center for Numerical Methods in Engineering
−	|-
−	| Universidad Politécnica de Cataluña
−	|-
−	| Gran Capitán s/n, 08034 Barcelona, Spain
−	|-
−	| E–mail: [mailto:onate@cimne.upc.edu onate@cimne.upc.edu]
−	|-
−	| Web page: [http://www.cimne.upc.es www.cimne.upc.es]
−	|}
−
−	==Abstract==
−
−	We present a general formulation for incompressible fluid flow analysis using the finite element method (FEM). The necessary stabilization for dealing with convective effects and the incompressibility condition are introduced via the Finite Calculus (FIC) method using a matrix form of the stabilization parameters. This allows to model a wide range of fluid flow problems for low and high Reynolds numbers flows without introducing a turbulence model.  Examples of application to the analysis of incompressible flows with moderate and large Reynolds numbers are presented.
−
−	'''keywords''' Stabilized formulation, incompressible fluid, finite calculus, finite element method, high Reynold's numbers, turbulence model.
−
	==1 INTRODUCTION==		==1 INTRODUCTION==

---

Latest revision as of 12:40, 7 October 2016

1 INTRODUCTION

Much effort has been spent in developing the so called stabilized numerical methods overcoming the two main sources of instability in incompressible flow analysis, namely those originated by the high values of the convective terms and those induced by the difficulty in satisfying the incompressibility condition.

The solution of above problems in the context of the finite element method (FEM) has been attempted in a number of ways. The first attempts to correct the underdiffusive character of the Galerkin FEM for high convection flows were based in adding some kind of artificial viscosity terms to the standard Galerkin equations 1.

A popular way to overcome the problems with the incompressibility constraint is by introducing a pseudo-compressibility in the flow and using implicit and explicit algorithms developed for this kind of problems such as artificial compressibility schemes 3 and preconditioning techniques 6. State of the art FEM schemes for fluid flow analysis with good stabilization properties for the convective and incompressibility terms are based in Petrov-Galerkin (PG) techniques. The background of PG methods are the non-centred (upwind) schemes for computing the first derivatives of the convective operator in finite difference and finite volume methods 7. A general class of stabilized PG FEM has been recently developed where the standard Galerkin variational form of the momentum and mass balance equations is extended with adequate residual-based terms in order to achieve a stabilized numerical scheme. References 9 list some of the more popular stabilized FEM of this kind. A review of many of these methods can be found in 1.

In this paper a stabilized FEM for incompressible flows is derived via a finite calculus (FIC) approach 32. The FIC method is based in invoking the balance of fluxes in a fluid domain of finite size. This introduces naturally additional terms in the classical differential equations of momentum and mass balance of infinitesimal fluid mechanics which are a function of characteristic length dimensions related to the element size in the discretized problem. The FIC terms in the modified governing equations provide the necessary stabilization to the discrete equations obtained via the standard Galerkin FEM. The FIC/FEM formulation allows to use low order finite elements (such as linear triangles and tetrahedra) with equal order approximations for the velocity and pressure variables.

The FIC/FEM formulation has proven to be very effective for the solution of a wide class of problems, such as convection-diffusion 32 and convection-diffusion-reaction 40 involving arbitrary high gradients, incompressible flow problems accounting for free surface effects and fluid-structure interaction situations 32 and quasi and fully incompressible problems in solid mechanics 52.

The FIC equations for incompressible flow derived in previous works of the authors assumed that the dimensions of the domain where the momentum conservations law was enforced remain the same independently of the direction along which balance of momentum is imposed. As a consequence, each of the resulting FIC momentum equations contain the same characteristic dimensions which can be grouped in a characteristic distance vector. In this paper, a refined FIC momentum equations are derived by accepting that the dimensions of the momentum balance domain are different for each of the momentum equations. This introduces a matrix form of the characteristic distances and of the corresponding FIC terms which have better intrinsic stabilization properties.

The idea of a matrix form of the stabilization parameters is close to the element-matrix-based and element-vector-based stabilization parameters proposed in 56 where different intrinsic time parameters wsere defined separately for each degree of freedom of the equation system.

Stabilized FEM have been successfully used in the past to solve a wide range of fluid mechanics problems. The intrinsice dissipative properties of the stabilization terms (which can interpreted as an additional viscosity) typically suffices to yield good results for low and moderate values of the Reynolds number (${\textstyle Re}$). For high values of ${\textstyle Re}$ most stabilized FEM fail to provide physically meaningful results and the numerical solution is often unstable or inaccurate. The introduction of a turbulence model is mandatory in order to obtain meaningful results in these cases.

The relationship between the additional dissipation introduced by the turbulence model and the intrinsic dissipative properties of stabilized FEM is an open topic which is attracting increasing attention in the CFD community. It is clear that both remedies (the turbulence model and the stabilization terms) play a similar role in the numerical solution, i.e. that of ensuring a solution which is ``physically sound´´ and as accurate as possible.

It is our belief that the matrix stabilization terms introduced by the FIC/FEM formulation here presented allow to model accurately high ${\textstyle Re}$ number flows without the need of introducing any turbulence model. The background of this belief originates in the positive experiences in the application of a very similar formulation for solving advection-diffusion and advection-diffusion-reaction problems with arbitrary sharp gradients without introducing any transverse dissipation terms 39. The extension of these ideas to the Navier-Stokes equations described here provides a straightforward procedure for solving a wide class of flow problems from low to high Reynolds numbers, as demonstrated by the good results presented in the paper.

The layout of the paper is the following. In the next section the FIC equations for incompressible flows with matrix stabilization terms are presented. The finite element discretization is introduced and the resulting matrix equations are detailed. A fractional step scheme for the transient solution is detailed. Examples of applications to the 2D analysis of flows passing a backward facing step and a cylinder at different Reynolds numbers are presented.

2 GENERAL FIC EQUATIONS FOR VISCOUS INCOMPRESSIBLE FLOW

The FIC governing equations for a viscous incompressible fluid can be written in an Eulerian frame of reference as

Momentum

${\displaystyle r_{m_{i}}-{\underline {{1 \over 6}h_{ij}{\partial r_{m_{i}} \over \partial x_{j}}}}{1 \over 2}h_{ij}{\partial r_{t_{z}} \over \partial x_{j}}=34\qquad {\hbox{in }}\Pi }$

(1)

Mass balance

${\displaystyle r_{d}-{\underline {{1 \over 2}h_{j}{\partial r_{d} \over \partial x_{j}}}}{1 \over 2}h_{j}{\partial r_{d} \over \partial x_{j}}=0\qquad {\hbox{in }}\Omega }$

(2)

where

${\displaystyle r_{m_{i}}=\rho \left({\partial u_{i} \over \partial t}+u_{j}{\partial u_{i} \over \partial x_{j}}\right)+{\partial p \over \partial x_{i}}-{\partial s_{ij} \over \partial x_{j}}-b_{i}}$

(3a)

${\displaystyle r_{d}={\partial u_{i} \over \partial x_{i}}\qquad i,j=1,n_{d}}$

(3b)

Above ${\textstyle \Omega }$ is the analysis domain, ${\textstyle n_{d}}$ is the number of space dimensions (${\textstyle n_{d}=2}$ for 2D problems), ${\textstyle u_{i}}$ is the velocity along the ith global axis, ${\textstyle \rho }$ is the (constant) density of the fluid, ${\textstyle p}$ is the absolute pressure (defined positive in compression), ${\textstyle b_{i}}$ are the body forces and ${\textstyle s_{ij}}$ are the viscous deviatoric stresses related to the viscosity ${\textstyle \mu }$ by the standard expression

${\displaystyle s_{ij}=2\mu \left({\dot {\varepsilon }}_{ij}-\delta _{ij}{1 \over 3}{\partial u_{k} \over \partial x_{k}}\right)}$

where ${\textstyle \delta _{ij}}$ is the Kronecker delta and the strain rates ${\textstyle {\dot {\varepsilon }}_{ij}}$ are

${\displaystyle {\dot {\varepsilon }}_{ij}={1 \over 2}\left({\partial u_{i} \over \partial x_{j}}+{\partial u_{j} \over \partial x_{i}}\right)}$

(4)

The FIC boundary conditions are

${\displaystyle n_{j}\sigma _{ij}-t_{i}+{\underline {{1 \over 2}h_{ij}n_{j}r_{m_{i}}}}{1 \over 2}h_{ij}n_{j}r_{m_{i}}=0\quad {\hbox{on }}\Gamma _{t}}$

(6a)

${\displaystyle u_{j}-u_{j}^{p}=0\quad {\hbox{on }}\Gamma _{u}}$

(6b)

and the initial condition is ${\textstyle u_{j}=u_{j}^{0}}$ for ${\textstyle t=t_{0}}$.

Summation convention for repeated indices in products and derivatives is used unless otherwise specified.

In Eqs.(13) and (14) ${\textstyle t_{i}}$ and ${\textstyle u_{j}^{p}}$ are surface tractions and prescribed displacements on the boundaries ${\textstyle \Gamma _{t}}$ and ${\textstyle \Gamma _{u}}$, respectively, ${\textstyle n_{j}}$ are the components of the unit normal vector to the boundary and ${\textstyle \sigma _{ij}}$ are the total stresses given by ${\textstyle \sigma _{ij}=s_{ij}-\delta _{ij}p}$.

The ${\textstyle h_{ij}}$ and ${\textstyle h_{j}}$ are characteristic distances of the domain where balance of momentum and mass is enforced. In Eq.(7) these lengths define the domain where equilibrium of boundary tractions is established 32. In the discretized problem the characteristic distances become of the order of the typical element dimensions. Note that by making these distances equal to zero the standard infinitessimal form of the fluid mechanics equations is recovered 1.

Eqs.(1)–(6) are the starting point for deriving stabilized FEM for solving the incompressible Navier-Stokes equations. The underlined FIC terms in Eq.(1) are essential to overcome the numerical instabilities due to the convective terms in the momentum equations, whereas the underlined terms in Eq.(2) take care of the instabilities due to the incompressibility constraint. An interesting feature of the FIC formulation is that it allows to use equal order interpolation for the velocity and pressure variables 43.

Remark 1. In previous work of the authors the characteristic distances in the momentum equations had a vector form, i.e. the FIC momentum equations were written as

${\displaystyle r_{i}-{1 \over 2}h_{j}{\partial r_{i} \over \partial x_{j}}=0}$

(7a)

or

${\displaystyle r_{i}-{1 \over 2}{h}_{j}\nabla r_{i}=0}$

(7b)

where (for 2D problems) ${\textstyle {h}=[h_{1},h_{2}]^{T}}$ is the characteristic length vector 32.

The difference of Eqs.(7) with Eq.(1) is that the characteristic distances have now a matrix form, i.e. the expanded form of the momentum equations (1) is (for 2D problems)

${\displaystyle {\begin{array}{l}r_{1}-\displaystyle {1 \over 2}\left(h_{11}\displaystyle {\partial r_{1} \over \partial x_{1}}+h_{12}\displaystyle {\partial r_{1} \over \partial x_{2}}\right)=0\\r_{2}-\displaystyle {1 \over 2}\left(h_{21}\displaystyle {\partial r_{1} \over \partial x_{1}}+h_{22}\displaystyle {\partial r_{1} \over \partial x_{2}}\right)=0\end{array}}}$

(8)

The rationale of Eqs.(8) is briefly explained in the Appendix.

The matrix of stabilization parameters H is defined as (for 2D problems)

${\displaystyle {H}=\left[{\begin{matrix}h_{11}&h_{12}\\h_{21}&h_{22}\\\end{matrix}}\right]}$

(9)

Remark 2. Note that the characteristic distances in the FIC mass conservation equation (2) have a vector form. As mentioned above distances ${\textstyle h_{1}}$ and ${\textstyle h_{2}}$ in Eq.(2) (for 2D problems) denote the dimensions of the domain where balance of mass is globally enforced (see Appendix). This is a basic difference with the momentum equations where the momentum balance law is applied along each global coordinate direction.

2.1 Stabilized integral forms

From the momentum equations it can be obtained 43

${\displaystyle {\partial r_{d} \over \partial x_{i}}\simeq {h_{ii} \over 2a_{i}}{\partial r_{m_{i}} \over \partial x_{j}}\quad ,\quad {\hbox{no sum in }}i}$

(10a)

where

${\displaystyle a_{i}={2\mu \over 3}+{\rho u_{i}h_{i} \over 2}\quad ,\quad {\hbox{no sum in }}i}$

(10b)

Substituting Eq.(15) into Eq.(8) and retaining the terms involving the derivatives of ${\textstyle r_{m_{i}}}$ with respect to ${\textstyle x_{i}}$ only, leads to the following alternative expression for the stabilized mass balance equation

${\displaystyle r_{d}-\sum \limits _{i=1}^{n_{d}}\tau _{i}{\partial r_{m_{i}} \over \partial x_{i}}=0}$

with

${\displaystyle \tau _{i}=\left({8\mu \over 3h_{ii}{\bar {h}}_{i}}+{2\rho u_{i} \over h_{ii}}\right)^{-1}}$

(11)

The ${\textstyle \tau _{i}}$'s in Eq.(11) when multiplied by the density are equivalent to the intrinsic time parameters, seen extensively in the stabilization literature. The interest of Eq.(11) is that it introduces the first space derivatives of the momentum equations into the mass balance equation. These terms have intrinsic good stability properties as explained next.

The weighted residual form of the momentum and mass balance equations (Eqs.(1) and (11)) is written as

${\displaystyle \int _{\Omega }\delta u_{i}\left[r_{m_{i}}-{h_{ij} \over 2}{\partial r_{m_{i}} \over \partial x_{j}}\right]d\Omega +\int _{\Gamma _{t}}\delta u_{i}(\sigma _{ij}n_{j}-t_{i}+{h_{j} \over 2}n_{j}r_{m_{i}})d\Gamma =0}$

(12)

${\displaystyle \int _{\Omega }q\left[r_{d}-\sum \limits _{i=1}^{n_{d}}\tau _{i}{\partial r_{m_{i}} \over \partial x_{i}}\right]d\Omega =0}$

(13)

where ${\textstyle \delta u_{i}}$ and ${\textstyle q}$ are arbitrary weighting functions representing virtual velocities and virtual pressure fields. Integrating by parts the ${\textstyle r_{m_{i}}}$ terms in Eqs.(13) and (14) leads to

${\displaystyle \int _{\Omega }\delta u_{i}r_{m_{i}}d\Omega +\int _{\Gamma _{t}}\delta u_{i}(\sigma _{ij}n_{j}-t_{i})d\Gamma +\int _{\Omega }{h_{ij} \over 2}{\partial \delta u_{i} \over \partial x_{j}}r_{m_{i}}d\Omega =0}$

(15a)

${\displaystyle \int _{\Omega }qr_{d}d\Omega +\int _{\Omega }\left[\sum \limits _{i=1}^{n_{d}}\tau _{i}{\partial q \over \partial x_{i}}r_{m_{i}}\right]d\Omega -\int _{\Gamma }\left[\sum \limits _{i=1}^{n_{d}}q\tau _{i}n_{i}r_{m_{i}}\right]d\Gamma =0}$

(15b)

We will neglect hereonwards the third integral in Eq.(15b) by assuming that ${\textstyle r_{m_{i}}}$ is negligible on the boundaries. The deviatoric stresses and the pressure terms in the first integral of Eq.(15a) are integrated by parts in the usual manner. The resulting momentum and mass balance equations are

${\displaystyle {\begin{array}{r}\displaystyle \int _{\Omega }\left[\delta u_{i}\rho \left({\partial u_{i} \over \partial t}+u_{j}{\partial {u_{i}} \over \partial x_{j}}\right)+{\partial \delta u_{i} \over \partial x_{j}}\left(\mu {\partial u_{i} \over \partial x_{j}}-\delta _{ij}p\right)\right]d\Omega -\int _{\Omega }\delta u_{i}b_{i}d\Omega -\qquad \\\displaystyle -\int _{\Gamma _{t}}\delta u_{i}t_{i}d\Gamma +\int _{\Omega }{h_{ij} \over 2}{\partial \delta u_{i} \over \partial x_{j}}r_{m_{i}}d\Omega =0\qquad \end{array}}}$

(16a)

${\displaystyle \int _{\Omega }q{\partial u_{i} \over \partial x_{i}}d\Omega +\int _{\Omega }\left[\sum \limits _{i=1}^{n_{d}}\tau _{i}{\partial q \over \partial x_{i}}r_{m_{i}}\right]d\Omega =0}$

(16b)

In the derivation of the viscous term in Eq.(16a) we have used the following identity holding for incompressible fluids (prior to the integration by parts)

${\displaystyle {\partial s_{ij} \over \partial x_{j}}=2\mu {\partial \varepsilon _{ij} \over \partial x_{j}}=\mu {\partial ^{2}u_{i} \over \partial x_{j}\partial x_{j}}}$

2.2 Convective and pressure gradient projections

The computation of the residual terms are simplified if we introduce the convective and pressure gradient projections ${\textstyle c_{i}}$ and ${\textstyle \pi _{i}}$, respectively defined as

${\displaystyle {\begin{array}{l}\displaystyle c_{i}=r_{m_{i}}-\rho u_{j}{\partial u_{i} \over \partial x_{j}}\\\displaystyle \pi _{i}=r_{m_{i}}-{\partial p \over \partial x_{i}}\end{array}}}$

(17)

We can express ${\textstyle r_{m_{i}}}$ in Eqs.(16a) and (16b) in terms of ${\textstyle c_{i}}$ and ${\textstyle \pi _{i}}$, respectively which then become additional variables. The system of integral equations is now augmented in the necessary number of equations by imposing that the residual ${\textstyle r_{m_{i}}}$ vanishes (in average sense) for both forms given by Eqs.(18). This gives the final system of governing equation as:

${\displaystyle \int _{\Omega }\left[\delta u_{i}\rho \left({\partial u_{i} \over \partial t}+u_{j}{\partial {u_{i}} \over \partial x_{j}}\right)+{\partial \delta u_{i} \over \partial x_{j}}\left(\mu {\partial u_{i} \over \partial x_{j}}-\delta _{ij}p\right)\right]d\Omega -\int _{\Omega }\delta u_{i}b_{i}d\Omega -}$ ${\displaystyle -\int _{\Gamma _{t}}\delta u_{i}t_{i}d\Gamma +\int _{\Omega }{h_{ik} \over 2}{\partial (\delta u_{i}) \over \partial x_{k}}\left(\rho u_{j}{\partial {u_{i}} \over \partial x_{j}}+c_{i}\right)d\Omega =0}$

(18)

${\displaystyle \int _{\Omega }q{\partial u_{i} \over \partial x_{i}}d\Omega +\int _{\Omega }\sum \limits _{i=1}^{n_{d}}\tau _{i}{\partial q \over \partial x_{i}}\left({\partial p \over \partial x_{i}}+\pi _{i}\right)d\Omega =0}$

(19)

${\displaystyle \int _{\Omega }\delta c_{i}\rho \left(\rho u_{j}{\partial {u_{i}} \over \partial x_{j}}+c_{i}\right)d\Omega =0\qquad {\hbox{no sum in }}i}$

(20)

${\displaystyle \int _{\Omega }\delta \pi _{i}\tau _{i}\left({\partial p \over \partial x_{i}}+\pi _{i}\right)d\Omega =0\qquad {\hbox{no sum in }}i}$

(21)

with ${\textstyle i,j,k=1,n_{d}}$. In Eqs.(21) and (22) ${\textstyle \delta c_{i}}$ and ${\textstyle \delta \pi _{i}}$ are appropriate weighting functions and the ${\textstyle \rho }$ and ${\textstyle \tau _{i}}$ weights are introduced for convenience.

We note that accounting for the convective and pressure gradient projections enforces the consistency of the formulation as it ensures that the stabilization terms in Eqs.(19) and (20) have a residual form which vanishes for the ``exact´´ solution. Neglecting these terms can reduce the accuracy of the numerical solution and it makes the formulation more sensitive to the value of the stabilization parameters as shown in references 51.

3 FINITE ELEMENT DISCRETIZATION

We choose ${\textstyle C^{\circ }}$ continuous linear interpolations of the velocities, the pressure, the convection projections ${\textstyle c_{i}}$ and the pressure gradient projections ${\textstyle \pi _{i}}$ over 3-noded triangles (2D) and 4-noded tetrahedra (3D). The linear interpolations are written as

${\displaystyle {\begin{array}{l}\displaystyle u_{i}=N^{k}{\bar {u}}_{i}^{k}\quad ,\quad p=N^{k}{\bar {p}}^{k}\\\displaystyle c_{i}=N^{k}{\bar {c}}_{i}^{k}\quad ,\quad \pi _{i}=N^{k}{\bar {\pi }}_{i}^{k}\end{array}}}$

(22)

where the sum goes over the number of nodes of each element ${\textstyle n}$ (${\textstyle n=3/4}$ for triangles/tetrahedra), ${\textstyle {\bar {(\cdot )}}^{k}}$ denotes the nodal variables and ${\textstyle N^{k}}$ are the linear shape functions 1.

Substituting the approximations (23) into Eqs.(19)–(22) and choosing the Galerking form with ${\textstyle \delta u_{i}=q=\delta c_{i}=\delta \pi _{i}=N^{i}}$ leads to following system of discretized equations

${\displaystyle \displaystyle {M}{\dot {\bar {u}}}+{H}{\bar {u}}-{G}{\bar {p}}+{C}{\bar {c}}={f}}$

(24a)

${\displaystyle \displaystyle {G}^{T}{\bar {u}}+{\hat {L}}{\bar {p}}+{Q}{\bar {\pi {\hbox{ }}}}={0}}$

(24b)

${\displaystyle \displaystyle {\hat {C}}{\bar {u}}+{M}{\bar {c}}={0}}$

(24c)

${\displaystyle \displaystyle {Q}^{T}{\bar {p}}+{\hat {M}}{\bar {\pi {\hbox{ }}}}={0}}$

(24d)

where

${\displaystyle {H}={A}+{K}+{\hat {K}}}$

If we denote the node indexes with superscripts ${\textstyle a,b}$, the space indices with subscripts ${\textstyle i,j}$, the element contributions to the components of the arrays involved in these equations are (${\textstyle i,j=1,3}$ for 3D problems)

${\displaystyle {\begin{array}{l}\displaystyle M_{ij}^{ab}=\left(\int _{\Omega ^{e}}\rho N^{a}N^{b}d\Omega \right)\delta _{ij}\quad ,\quad A_{ij}^{ab}=\left(\int _{\Omega ^{e}}\rho N^{a}({u}^{T}{\nabla }N^{b})d\Omega \right)\delta _{ij}\\\\\displaystyle {K}_{ij}^{ab}=\left(\int _{\Omega ^{e}}\mu {\nabla {\hbox{ }}}^{T}N^{a}{\nabla {\hbox{ }}}N^{b}d\Omega \right)\delta _{ij}\quad ,\quad {\nabla {\hbox{ }}}=\left[{\partial \over \partial x_{1}},{\partial \over \partial x_{2}},{\partial \over \partial x_{3}}\right]^{T}\\\\\displaystyle {\hat {K}}_{ij}^{ab}=\left({1 \over 2}\int _{\Omega ^{e}}h_{ij}{\partial N^{a} \over \partial x_{j}}(\rho {u}^{T}{\nabla {\hbox{ }}}N^{b})d\Omega \right)\delta _{ij}\quad ,\quad {G}_{i}^{ab}=\int _{\Omega ^{e}}{\partial N^{a} \over \partial x_{i}}N^{b}d\Omega \\\\\displaystyle {C}=\left[{\begin{matrix}{C}_{1}\\{C}_{2}\\{C}_{3}\\\end{matrix}}\right]\quad ,\quad {C}_{i}^{ab}={1 \over 2}\int _{\Omega ^{e}}h_{ij}{\partial N^{a} \over \partial x_{j}}N^{b}d\Omega \\\end{array}}}$

(25)

${\displaystyle {\begin{array}{l}\displaystyle {\hat {L}}^{ab}=\int _{\Omega ^{e}}({\nabla {\hbox{ }}}^{T}N^{a})[\tau ]{\nabla {\hbox{ }}}N^{b}d\Omega \quad ,\quad [\tau ]=\left[{\begin{matrix}\tau _{1}&0&0\\0&\tau _{2}&0\\0&0&\tau _{3}\\\end{matrix}}\right]\\\\\displaystyle {Q}=[{Q}_{1},{Q}_{2},{Q}_{3}]\quad ,\quad \displaystyle Q_{i}^{ab}=\int _{\Omega ^{e}}\tau _{i}{\partial N^{a} \over \partial x_{i}}N^{b}d\Omega \quad \quad {\hbox{no sum in }}i\end{array}}}$

(26)

${\displaystyle {\begin{array}{l}\displaystyle {\hat {C}}=[{\hat {C}}_{1},{\hat {C}}_{2},{\hat {C}}_{3}]\quad ,\quad \displaystyle {\hat {C}}_{1}^{ab}={\hat {C}}_{2}^{ab}={\hat {C}}_{3}^{ab}=\int _{\Omega ^{e}}\rho ^{2}N^{a}({u}^{T}{\nabla {\hbox{ }}}N^{b})d\Omega \\\\\displaystyle {\hat {M}}_{ij}^{ab}=\left(\int _{\Omega ^{e}}\tau _{i}N^{a}N^{b}d\Omega \right)\delta _{ij}\quad ,\quad \displaystyle {f}_{i}^{a}=\int _{\Omega ^{e}}N^{a}f_{i}d\Omega +\int _{\Gamma ^{e}}N^{a}t_{i}d\Gamma \end{array}}}$

It is understood that all the arrays are matrices (except ${\textstyle {f}}$ which is a vector) whose components are obtained by grouping together the left indices in the previous expressions (${\textstyle a}$ and possibly ${\textstyle i}$) and the right indices (${\textstyle b}$ and possibly ${\textstyle j}$).

Note that the stabilization matrix ${\textstyle {\hat {K}}}$ in Eq.(25) adds additional orthotropic diffusivity terms of value ${\textstyle \rho \displaystyle {h_{ij}u_{l} \over 2}}$.

The overall stabilization terms introduced by the FIC formulation above presented have the intrinsic capacity to ensure physically sound numerical solutions for a wide spectrum of Reynolds numbers without the need of introducing additional turbulence modelling terms. This interesting property is validated in the solution of the examples presented in a next section.

3.1 Transient solution scheme

The solution in time of the system of Eqs.(24) can be written in general form as

${\displaystyle {M}\displaystyle {1 \over \Delta t}({\bar {u}}^{n+1}-{\bar {u}}^{n})+{H}^{n+\theta }{\bar {u}}^{n+\theta }-{G}{\bar {p}}^{n+\theta }+{C}^{n+\theta }{\bar {c}}^{n+\theta }={f}^{n+\theta }}$

(27a)

${\displaystyle {G}^{T}{\bar {u}}^{n+\theta }+{\hat {L}}^{n+\theta }{\bar {p}}^{n+\theta }+{Q}{\bar {\pi {\hbox{ }}}}^{n+\theta }={0}}$

(27b)

${\displaystyle {\hat {C}}^{n+\theta }{\bar {u}}^{n+\theta }+{M}{\bar {c}}^{n+\theta }={0}}$

(27c)

${\displaystyle {G}^{T}{\bar {p}}^{n+\theta }+{\hat {M}}^{n+\theta }{\bar {\pi {\hbox{ }}}}^{n+\theta }={0}}$

(27d)

where ${\textstyle {H}^{n+\theta }={H}({u}^{n+\theta })}$, etc and the parameter ${\textstyle \theta \in [0,1]}$. The direct monolitic solution of Eqs.(27) is possible using an adequate iterative scheme 52. However, in our work we have used the fractional step method described next.

4 FRACTIONAL STEP METHOD

A fractional step scheme is derived by noting that the discretized momentum equation (27a) can be split into the two following equations

${\displaystyle {M}\displaystyle {1 \over \Delta t}({\tilde {u}}^{n+1}-{\bar {u}}^{n})+{H}^{n+\theta }{\bar {u}}^{n+\theta }-\alpha {G}{\bar {p}}^{n}+{C}^{n+\theta }{\bar {c}}^{n+\theta }={f}^{n+\theta }}$

(28a)

${\displaystyle {M}\displaystyle {1 \over \Delta t}({\bar {u}}^{n+1}-{\tilde {u}}^{n+1})-{G}({\bar {p}}^{n+1}-\alpha {\bar {p}}^{n})={0}}$

(28b)

In Eqs.(28) ${\textstyle {\tilde {u}}^{n+1}}$ is a predicted value of the velocity at time ${\textstyle n+1}$ and ${\textstyle \alpha }$ is a variable whose values of interest are zero and one. For ${\textstyle \alpha =0}$ (first order scheme) the splitting error is of order ${\textstyle 0(\Delta t)}$, whereas for ${\textstyle \alpha =1}$ (second order scheme) the error is of order ${\textstyle 0(\Delta t^{2})}$ 52. We have chosen ${\textstyle \alpha =1}$ for the solution of the examples presented in the paper.

Eqs.(28) are completed with the following three equations emanating from Eqs.(27b-d)

${\displaystyle {G}^{T}{\bar {u}}^{n+1}+{\hat {L}}^{n}{\bar {p}}^{n+1}+{Q}{\bar {\pi {\hbox{ }}}}^{n}={0}}$

(29a)

${\displaystyle {\hat {C}}^{n+1}{\bar {u}}^{n+1}+{M}{\bar {c}}^{n+1}={0}}$

(29b)

${\displaystyle {Q}^{T}{\bar {p}}^{n+1}+{\hat {M}}^{n+1}{\bar {\pi {\hbox{ }}}}^{n+1}={0}}$

(29c)

The value of ${\textstyle {\bar {u}}^{n+1}}$ obtained from Eq.(28b) is substituted into Eq.(29a) to give

${\displaystyle {G}^{T}{\tilde {u}}^{n+1}+\Delta t{G}^{T}{M}^{-1}{G}({\bar {p}}^{n+1}-\alpha {\bar {p}}^{n})+{\hat {L}}^{n}{p}^{n+1}+{Q}{\bar {\pi {\hbox{ }}}}^{n}={0}}$

The product ${\textstyle {G}^{T}{M}^{-1}{G}}$ can be approximated by a laplacian matrix, i.e.

${\displaystyle {G}^{T}{M}^{-1}{G}={1 \over \rho }{L}\quad {\hbox{with }}L^{ab}=\int _{\Omega ^{e}}{\nabla {\hbox{ }}}^{T}N^{a}{\nabla {\hbox{ }}}N^{b}d\Omega }$

(30)

where ${\textstyle L^{ab}}$ are the element contributions to ${\textstyle {L}}$.

The steps of the fractional step scheme are:

Step 1

Eq.(28a) is linearized as

${\displaystyle {M}{{\tilde {u}}^{n+1}-{\bar {u}}^{n} \over \Delta t}+{\tilde {H}}^{n+\theta }{\tilde {u}}^{n+\theta }-\alpha {G}{\bar {p}}^{n}+{\tilde {C}}^{n+\theta }{\bar {c}}^{n}={\bar {f}}^{n+\theta }}$

(31)

where ${\textstyle {\tilde {u}}^{n+\theta }=\theta {\tilde {u}}^{n+1}+(1-\theta ){\bar {u}}^{n}}$, ${\textstyle {\tilde {H}}^{n+\theta }={H}({\tilde {u}}^{n+\theta })}$, and ${\textstyle {\tilde {C}}^{n+\theta }={C}({\tilde {u}}^{n+\theta })}$. We have chosen in our computation ${\textstyle \theta =0}$. For this value, the fractional nodal velocities ${\textstyle {\tilde {u}}^{n+1}}$ can be explicitely computed from Eq.(32) by

${\displaystyle {\tilde {u}}^{n+1}={\bar {u}}^{n}-\Delta t{M}_{d}^{-1}[{\tilde {H}}^{n}{\bar {u}}^{n}-\alpha {G}{\bar {p}}^{n}+{\tilde {C}}^{n}{\bar {c}}^{n}-{\bar {f}}^{n}]}$

(32)

where ${\textstyle {M}_{d}}$ is the lumped diagonal form of M.

Step 2 Compute ${\textstyle {\bar {p}}^{n+1}}$ from Eq.(30) as

${\displaystyle {\bar {p}}^{n+1}=-[{\hat {L}}^{n}+{\Delta t \over \rho }{L}]^{-1}[{G}^{T}{\tilde {u}}^{n+1}-\alpha {\Delta t \over \rho }{L}{\bar {p}}^{n}+{Q}{\bar {\pi {\hbox{ }}}}^{n}]}$

(33)

Step 3 Compute ${\textstyle {\bar {u}}^{n+1}}$ explicitly from Eq.(28a) as

${\displaystyle {\bar {u}}^{n+1}={\tilde {u}}^{n+1}+\Delta t{M}_{d}^{-1}{G}({\bar {p}}^{n+1}-\alpha {\bar {p}}^{n})}$

(34)

Step 4 Compute ${\textstyle {\bar {c}}^{n+1}}$ explicitly from Eq.(29b) as

${\displaystyle {\bar {c}}^{n+1}=-{M}_{d}^{-1}{\hat {C}}^{n+1}{\bar {u}}^{n+1}}$

(35)

Step 5 Compute ${\textstyle {\bar {\pi {\hbox{ }}}}^{n+1}}$ explicitly from Eq.(29c) as

${\displaystyle {\bar {\pi {\hbox{ }}}}^{n+1}=-{\hat {M}}_{d}^{-1}{Q}^{T}{\bar {p}}^{n+1}}$

(36)

Above algorithm has improved stabilization properties versus the standard segregation methods due to the introduction of the laplacian matrix ${\textstyle {\hat {L}}}$ in Eq.(34) which emanates from the FIC stabilization terms.

The boundary conditions are applied as follows. No condition is applied in the computation of the fractional velocities ${\textstyle {\tilde {u}}^{n+1}}$ in Eq.(33). The prescribed velocities at the boundary are applied when solving for ${\textstyle {\bar {u}}^{n+1}}$ in the step 3. The prescribed pressures at the boundary are imposed by making ${\textstyle {\bar {p}}^{n}}$ equal to the prescribed pressure values.

5 STOKES FLOW

The formulation for a Stokes flow can be readily obtained simply by neglecting the convective terms in the general Navier-Stokes formulation. Consequently, the convective stabilization terms and the convective projection variables are not larger necessary. Also the intrinsic time parameters ${\textstyle \tau _{i}}$ take now the simpler form (see Eq.(12)):

${\displaystyle \tau _{i}={3h_{ii}h_{i} \over 8\mu }}$

(37)

The resulting discretized system of equations can be written as (see Eqs.(29))

${\displaystyle {\begin{array}{l}\displaystyle {M}{\dot {\bar {u}}}+{K}{\bar {u}}-{G}{\bar {p}}={f}\\\displaystyle {G}^{T}{\bar {u}}+{\hat {L}}{\bar {p}}+{Q}{\bar {\pi {\hbox{ }}}}={0}\\\displaystyle {Q}^{T}{\bar {p}}+{\hat {M}}{\bar {\pi {\hbox{ }}}}={0}\end{array}}}$

(38)

The fractional step algorithm of the previous section can now be implemented. We note that convergence of the predictor-corrector scheme is now faster due to the absence of the non linear convective terms in the momentum equation.

The steady-state form of Eqs.(39) can be expressed in matrix form as

${\displaystyle \left[{\begin{matrix}{K}&-{G}&{0}\\-{G}^{T}&-{\hat {L}}&-{Q}\\{0}&-{Q}^{T}&-{\hat {M}}\\\end{matrix}}\right]\left\{{\begin{matrix}{\bar {u}}\\{\bar {p}}\\{\bar {\pi {\hbox{ }}}}\\\end{matrix}}\right\}=\left\{{\begin{matrix}{f}\\{0}\\{0}\\\end{matrix}}\right\}}$

(39)

The system is symmetric and always positive definite and therefore leads to a non singular solution. This property holds for any interpolation function chosen for ${\textstyle {\bar {u}},{\bar {p}}}$ and ${\textstyle {\bar {\pi {\hbox{ }}}}}$, therefore overcoming the Babuŝka-Brezzi (BB) restrictions 1.

A reduced velocity-pressure formulation can be obtained by eliminating the pressure gradient projection variables ${\textstyle {\bar {\pi {\hbox{ }}}}}$ from the last equation to give

${\displaystyle \left[{\begin{matrix}{K}&-{G}\\-{G}^{T}&-({\hat {L}}-{Q}{\hat {M}}^{-1}{Q}^{T})\\\end{matrix}}\right]\left\{{\begin{matrix}{\bar {u}}\\{\bar {p}}\\\end{matrix}}\right\}=\left\{{\begin{matrix}{f}\\{0}\\\end{matrix}}\right\}}$

(40)

The reduction process is simplified by using a diagonal form of matrix ${\textstyle {\hat {M}}}$. Applications of this scheme to incompressible solid mechanics problems are reported in 52.

6 COMPUTATION OF THE CHARACTERISTIC DISTANCES

The computation of the stabilization parameters is a crucial issue as they affect both the stability and accuracy of the numerical solution. The different procedures to compute the stabilization parameters are typically based on the study of simplified forms of the stabilized equations. Contributions to this topic are reported in 11. Despite the relevance of the problem there still lacks a general method to compute the stabilization parameters for all the range of flow situations.

Recent work of the authors in the application of the FIC/FEM formulation to convection-diffusion problems with sharp arbitrary gradients 39 has shown that the stabilizing FIC terms take the form of a simple orthotropic diffusion if the balance equation is written in the principal curvature directions of the solution. Excellent results were reported in 39 by computing first the characteristic length distances along the principal curvature directions, followed by a standard transformation of the distances to global axes. The resulting stabilized finite element equations capture the high gradient zones in the vicinity of the domain edges (boundary layers) as well as the sharp gradients appearing randomly in the interior of the domain 39. The FIC/FEM thus reproduces the best features of both the so called transverse (cross-wind) dissipation or shock capturing methods 58.

The numerical computations are simplified without apparent loss of accuracy if the main principal curvature direction of the solution at each element point is approximated by the direction of the gradient vector at the element center. The second principal direction (for 2D problems) is taken in the orthogonal direction to the gradient. For linear triangles and quadrilaterals these directions are assumed to be constant within the element 39.

Above simple scheme has been extended in this work for the computation of the characteristic distances ${\textstyle h_{ij}}$ for the momentum equations. As for the length parameters ${\textstyle h_{i}}$ in the mass conservation equation, the simplest assumption ${\textstyle h_{i}=h_{ii}}$ has been taken. Details of the algorithm for computing ${\textstyle h_{ij}}$ are given next (the method is explained for 2D problems although it is readily extendible to 3D problems).

() Definition of the principal curvature direction ${\displaystyle {\vec {\xi }}_{1}^{i}}$ along the gradient of ${\displaystyle u_{i}}$.

Definition of the element characteristic distances ${\displaystyle l_{i1}}$ and ${\displaystyle l_{i2}}$ corresponding to the ${\displaystyle i}$th momentum equation.

For the ${\textstyle i}$-th momentum balance equation and every step of the fractional step method described in Section 4:

1. A coordinate system ${\textstyle {\vec {\xi }}_{1}^{i}}$, ${\textstyle {\vec {\xi }}_{2}^{i}}$ is defined at each element point such that ${\textstyle {\vec {\xi }}_{i}^{1}}$ is aligned with the gradient of ${\textstyle u_{i}}$ (${\textstyle {\vec {\xi }}_{1}^{i}={\vec {\nabla u_{i}}}}$) and ${\textstyle {\vec {\xi }}_{2}^{i}}$ is orthogonal to ${\textstyle {\vec {\xi }}_{1}^{i}}$ in anticlockwise sense (Figure 1). The angle that ${\textstyle {\vec {\xi }}_{1}^{i}}$ forms with the global ${\textstyle x_{1}}$ axis is defined as ${\textstyle \alpha _{i}}$. Recall that upper and lower index ${\textstyle i}$ denotes the ${\textstyle i}$th momentum equation.
2. The element characteristic distances ${\textstyle l_{i1}}$ and ${\textstyle l_{i2}}$ are defined as the maximum projections of the element sides along the ${\textstyle {\vec {\xi }}_{1}^{i}}$ and ${\textstyle {\vec {\xi }}_{2}^{i}}$ axes, respectively (Figure 2).
3. The characteristic distances ${\textstyle h_{i1}}$ and ${\textstyle h_{i2}}$ are computed as

   ${\displaystyle \left\{{\begin{matrix}h_{i1}\\h_{i2}\\\end{matrix}}\right\}=\left[{\begin{matrix}c_{i}&-s_{i}\\s_{i}&c_{i}\\\end{matrix}}\right]\left\{{\begin{matrix}{\bar {h}}_{i1}\\{\bar {h}}_{i2}\\\end{matrix}}\right\}\quad ,\quad i=1,2}$

   (41)

with ${\textstyle c_{i}=\cos \alpha _{i},s_{i}=\sin \alpha _{i}}$ and the local distances ${\textstyle {\bar {h}}_{i1}}$ and ${\textstyle {\bar {h}}_{i2}}$ are

${\displaystyle {\bar {h}}_{i1}=\left(\coth {\bar {\gamma }}_{ij}-{1 \over {\bar {\gamma }}_{ij}}\right)l_{ij}\quad ,\quad \gamma _{ij}={{\bar {u}}_{j}l_{ij} \over 2\mu }\quad j=1,2}$

(42)

where ${\textstyle {\bar {u}}_{1}}$ and ${\textstyle {\bar {u}}_{2}}$ are the components of the velocity vector along the local axes ${\textstyle {\vec {\xi }}_{1}^{i}}$ and ${\textstyle {\vec {\xi }}_{2}^{i}}$, respectively (Figure 1).

7 EXAMPLES

The examples were solved with the Tdyn code where the formulation here presented has been implemented. The Tdyn code can be downloaded from the webpage given in [58].

7.1 Backwards facing step at high Reynolds number

Figure 3 shows the geometry of the standard backwards facing step problem. The boundary conditions were the following: ${\textstyle u_{1}=1}$ and ${\textstyle u_{2}=0}$ were taken at the entry while ${\textstyle p=0}$ was assumed at the exit. Slipping conditions were assumed at the rest of the vertical and horizontal walls. A value of the kinematic viscosity ${\textstyle \nu =\displaystyle {\mu \over \rho }=2.1\times 10^{-5}}$ was taken giving a Reynolds number of ${\textstyle Re=\displaystyle {u_{free}H \over \nu }=47619}$ for ${\textstyle H=1}$ and ${\textstyle u_{free}=1}$.

Figure 3 also shows the relatively coarse mesh chosen of 30850 three-noded triangular elements and 15426 nodes. The contours of the horizontal and vertical velocities and details of the velocity vectors are shown in Figures 4 and 5, respectively. Figure 6 shows the distribution of the horizontal velocity along the bottom line starting from the vertical wall of the step. The point where the horizontal velocity changes sign indicates the end of the recirculation area.

The length of the circulation area computed from Figure 6 is 6.79. This value compares very well (3.2%) with the experimental value reported by Kim et al. [59] (see Table 1). The FIC/FEM results are remarkably accurate in comparison with other results reported in the literature obtained using ${\textstyle K-\varepsilon }$ and ${\textstyle K}$-tau turbulence models [60–63]. We note again that the FIC/FEM formulation does not include any additional turbulence terms.

Error creating thumbnail: File missing

Backwards facing step. Geometry and finite element mesh of 30850 three-noded triangles. Mesh detail at the vicinity of the step.

Backwards facing step. Geometry and finite element mesh of 30850 three-noded triangles. Mesh detail at the vicinity of the step.

File:Testeando 2016d-Eje6 1Fig3.png

Backwards facing step. Contours of horizontal (above) and vertical velocities.

Backwards facing step. Contours of horizontal (above) and vertical velocities.

Error creating thumbnail: File missing

Backwards facing step. Velocity vectors and recirculation distance ${\displaystyle D}$.

Distribution of the horizontal velocity along the bottom line starting from the vertical wall of the step. The circle shows the end of the recirculation region.

Table. 1 Backwards facing step. Length of the recirculation distance ${\displaystyle D}$ for ${\displaystyle Re=47619}$. Comparison of the FIC/FEM result with experimental data and with numerical results obtained using different turbulence models.

Model	Length ${\displaystyle D}$	Error range	Average error
Exp. [59]	6.0–7.0
K-${\textstyle \varepsilon }$ [60]	5.2	13.3–26%	19.6%
K-${\textstyle \varepsilon }$ [61]	5.88	2–16%	9%
K-${\textstyle \varepsilon }$ [62]	6.0	0–14%	7%
K-${\textstyle \varepsilon }$ [63]	6.2	13.7–11.4%	12.6%
K-Tau [63]	6.82	13.7-2.5%	8.1%
FIC/FEM	6.71	11.8-4.1%	7.9%

7.2 Flow past a cylinder. Computation of the Strouhal unstability

Figure 7 shows the geometry for the analysis of the flow past a cylinder of unit diameter (${\textstyle D}$). A unit horizontal velocity is prescribed at the inlet boundary and at the two horizontal walls. Zero pressure is prescribed at the outlet boundary. The dimensions of the analysis domain are ${\textstyle 27\times 27}$ units. The origin of the coordinate system has been sampled at the center of the cylinder located at a distance of 13.1 units from the entry wall. Zero velocity is prescribed at the cylinder wall. The kinematic viscosity is ${\textstyle \nu =0.01}$. Figure 8 shows the mesh of 91316 three-noded elements used for the computation. A detail of the mesh in the vicinity of the cylinder is also shown in Figure 8.

Error creating thumbnail: File missing

Flow past a cylinder of unit diameter. Analysis domain and boundary conditions.

Flow past a cylinder. Mesh of 91316 three-noded triangles used for the computations.

Error creating thumbnail: File missing

Flow past a cylinder, ${\displaystyle Re=100}$. Contour of the velocity vector modulus for ${\displaystyle t=100}$ secs.

The problem has been analyzed first for a value of the horizontal velocity at the entry of ${\textstyle u_{1}=1}$ giving a Reynolds number of ${\textstyle Re=100}$. Figures 9 and 10 respectively show the velocity modulus contours and the velocity vectors for ${\textstyle t=100}$secs.

Figure 11 shows images of the trajectory of a substance over a band of 2.45 units transported at the entry across the flow for ${\textstyle t=100}$secs. The picture shows clearly the oscillatory nature of the flow.

Figure 12 shows the oscillations of the horizontal velocity at the point ${\textstyle A}$ with coordinates (6.7, -1.02) for ${\textstyle t=100}$secs. The Strouhal number computed from the shedding frequency ${\textstyle n}$ as ${\textstyle S={nD \over \vert {u}\vert }}$ is ${\textstyle S=0.1702}$. This number compares very well with the experimental result available in the literature (see Figure 13).

The same problem was analyzed for a value of the kinematic viscosity ${\textstyle \nu =0.001}$ giving ${\textstyle Re=1000}$. The same mesh of 91316 linear triangles of Figure 7 was used. Figures 14–16 show respectively the velocity modulus contours, the velocity vectors in the vicinity of the cylinder for ${\textstyle t=100}$ secs. and the trajectories of a substance transported across the flow. Figure 17 finally shows the oscillations of the horizontal velocity at point ${\textstyle A}$. The computed value of the Strouhal number in this case was ${\textstyle S=0.2103}$. This value again coincides well with the reported experimental data (see Figure 13).

It is a well known fact that for ${\textstyle Re>300}$ the flow past a cylinder exhibits 3D features. In [64] results from 2D and 3D computation were compared for ${\textstyle Re=300}$ and 800. While 3D features were observed even at ${\textstyle Re=300}$ and more so at ${\textstyle Re=800}$, there were no large discrepances between the global flow parameters (such as drag, lift and Strouhal number) obtained from 2D and 3D computations. These conclusions justify the results of the 2D computations presented in the paper.

Error creating thumbnail: File missing

Flow past a cylinder, ${\displaystyle Re=100}$. Velocity vectors for ${\displaystyle t=100}$ secs.

Error creating thumbnail: File missing

Error creating thumbnail: File missing

Flow past a cylinder, ${\displaystyle Re=100}$. Trajectories of a substance over a band of 2.45 units at the entry transported across the flow for ${\displaystyle t=100}$ secs.

Error creating thumbnail: File missing

Flow past a cylinder, ${\displaystyle Re=100}$. Oscillations with time of the horizontal velocity at the point with coordinates ${\displaystyle A}$ (6.7–1.02).

Flow past a cylinder. Experimental (thick line) and computed (${\displaystyle {\blacklozenge }}$) values of the Strouhal number ${\displaystyle S}$ in terms of the Reynolds number. Experimental values taken from [2].

Error creating thumbnail: File missing

Flow past a cylinder, ${\displaystyle Re=1000}$. Contour of the velocity vector module for ${\displaystyle t=100}$ secs.

Error creating thumbnail: File missing

Flow past a cylinder, ${\displaystyle Re=1000}$. Velocity vectors for ${\displaystyle t=100}$ secs.

Error creating thumbnail: File missing

Error creating thumbnail: File missing

Flow past a cylinder, ${\displaystyle Re=1000}$. Trajectories of a substance over a band of 2.45 units at the entry transported across the flow for ${\displaystyle t=100}$ secs.

Error creating thumbnail: File missing

Flow past a cylinder, ${\displaystyle Re=1000}$. Oscillations with time of the horizontal velocity at the point with coordinates ${\displaystyle A}$ (6.7–1.02).

8 CONCLUSIONS

The finite calculus (FIC) form of the fluid mechanics equations is a good starting point for deriving stabilized FEM for solving a variety of incompressible fluid flow problems. The matrix stabilization terms introduced by the FIC formulation here presented allow to obtain physically sound solutions in the presence of sharp gradients occuring for high Reynolds numbers without the need of introducing a turbulence model. Good numerical solutions have been obtained in the 2D examples solved with relatively coarse meshes for moderate and high values of the Reynolds number. These preliminary results reinforce our idea that the stabilization terms introduced by the FIC formulation suffice to provide good results for problems for which turbulence models are required using alternative numerical methods. These results also confirm the close link between the stabilized methods and turbulence models, which surely will be the object of much research in the near future.

ACKNOWLEDGEMENTS

The authors thank Prof. S.R. Idelsohn and C.A. Felippa for many useful discussions.

APPENDIX

The FIC momentum equations in two dimensions (2D) are obtained by expressing the balance of momentum along the horizontal and vertical directions in the finite domains shown in Figures A.1 and A.2, respectively.

Error creating thumbnail: File missing

Finite domain where balance of momentum is imposed along the horizontal direction

Error creating thumbnail: File missing

Finite domain where balance of momentum is imposed along the vertical direction

The balance equation is written for each finite domain as

${\displaystyle \sum f_{i}d\Omega ={\partial \over \partial t}\int _{\Omega }\rho u_{i}d\Omega +\int _{\Gamma }(\rho u_{i}){u}^{T}{n}d\Gamma \quad i=1,2}$

(43)

where ${\textstyle f_{i}}$ includes the forces due to the stresses acting on the boundary of the balance domain and the body forces per unit area (Figures A.1 and A.2).

Expressing the values of the momentum and force terms at the end point of the balance domain in terms the values at an arbitrary point (such as the corner point A) using higher order Taylor expansions and retaining second order tems gives after some algebra 32 the FIC momentum equations along the ${\textstyle i}$th coordinate direction as

${\displaystyle r_{i}-{1 \over 2}h_{ij}{\partial r_{i} \over \partial x_{j}}=0\quad \quad i,j=1,2}$

(44)

with

${\displaystyle r_{i}:=\rho \left({\partial u_{i} \over \partial t}+u_{j}{\partial u_{i} \over \partial x_{j}}\right)+{\partial p \over \partial x_{i}}-{\partial \sigma _{ij} \over \partial x_{i}}-b_{i}}$

(45)

with ${\textstyle \sigma _{ij}=\tau _{ij}-p\delta _{ij}}$ where ${\textstyle \tau _{ij}}$ and ${\textstyle p}$ are the deviatoric stresses and the pressure, respectively.

Note that distance ${\textstyle h_{12}}$ is arbitrary when writting the balance of momentum along the ${\textstyle x_{1}}$ direction. The same applies for the distance ${\textstyle h_{21}}$ when deriving the balance equation along the ${\textstyle x_{2}}$ direction. Thus, in general, ${\textstyle h_{12}\not =h_{21}}$ and this explains the matrix form of the FIC momentum equations.

The FIC mass balance equation is obtained by invoking the balance of mass in the finite domain of Figure A.3

${\displaystyle \int _{\Gamma }\rho {u}^{T}{n}d\Gamma =0}$

(46)

Finite domain where balance of mass is enforced

Expanding the values of ${\textstyle \rho u_{i}}$ at the corner points in terms of the value at an arbitrary point gives of the mass balance domain the FIC mass balance equation as 32

${\displaystyle {\partial u_{i} \over \partial x_{i}}-{1 \over 2}h_{j}{\partial \over \partial x_{j}}\left({\partial u_{i} \over \partial x_{i}}\right)=0}$

Note that a matrix form of the characteristic distances is not obtained in this case as the mass balance equation expresses the conservation of the mass in the whole domain ABCD of Figure A.3 with dimensions ${\textstyle h_{1}}$ and ${\textstyle h_{2}}$. Distances ${\textstyle h_{1}}$ and ${\textstyle h_{2}}$ should be taken in general different from distances ${\textstyle h_{ij}}$ defining the domain where balance of momentum is enforced. In our computations we have however assumed that ${\textstyle h_{1}=h_{11}}$ and ${\textstyle h_{2}=h_{22}}$ for simplicity.

BIBLIOGRAPHY

[1] O.C. Zienkiewicz, R.C. Taylor and P. Nietharasu. The finite element method, 6th Edition, Volume 3, Elsevier, 2005.

[2] D.W. Kelly, S. Nakazawa and O.C. Zienkiewicz, ``A note on anisotropic balancing dissipation in the finite element method approximation to convective diffusion problems´´, Int. J. Num. Meth. Engng., 15, 1705-11, 1980.

[3] W.R. Briley, S.S. Neerarambam and D.L. Whitfield, ``Multigrid algorithm for three-dimensional incompressible high-Reynolds number turbulent flows´´, AIAA Journal, 33 (1), 2073–2079, 1995.

[4] J. Peraire, K. Morgan and J. Peiro, ``The simulation of 3D incompressible flows using unstructured grids´´, In Frontiers of Computational Fluid Dynamics, Caughey DA and Hafez MM. (eds), Chapter 16, J. Wiley, 1994.

[5] C. Sheng, L.K. Taylor and D.L. Whitfield, ``Implicit lower-upper/approximate-factorization schemes for incompressible flows´´ Journal of Computational Physics, 128 (1), 32–42, 1996.

[6] M. Storti, N. Nigro and S.R. Idelsohn, ``Steady state incompressible flows using explicit schemes with an optimal local preconditioning´´, Comput. Meth. Appl. Mech. Engrg., 124, 231–252, 1995.

[7] C. Hirsch, Numerical computation of internal and external flow, J. Wiley, Vol. 1 1988, Vol. 2, 1990.

[8] J.C. Heinrich, P.S. Hayakorn and O.C. Zienkiewicz, ``An upwind finite element scheme for two dimensional convective transport equations´´, Int. J. Num. Meth. Engng., 11, 131–143, 1977.

[9] S.R. Idelsohn, M. Storti and N. Nigro, ``Stability analysis of mixed finite element formulation with special mention to stabilized equal-order interpolations´´ Int. J. for Num. Meth. in Fluids, 20, 1003-1022, 1995.

[10] S.R. Idelsohn, N. Nigro, M. Storti and G. Buscaglia, ``A Petrov-Galerkin formulation for advection-reaction-diffusion´´, Comput. Meth. Appl. Mech. Engrg., 136, 27–46, 1996.

[11] A. Brooks and T.J.R. Hughes, ``Streamline upwind/Petrov-Galerkin formulation for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations´´, Comput. Methods Appl. Mech. Engrg., 32, 199–259, 1982.

[12] T.J.R. Hughes and T.E. Tezduyar, ``Finite element methods for first-order hyperbolic systems with particular emphasis on the compressible Euler equations´´, Comput. Meth. Appl. Mech. Engrg., 45, 217–284, 1984.

[13] T.E. Tezduyar and Y.J. Park, ``Discontinuity capturing finite element formulations for nonlinear convection-diffusion-reaction equations", Comput. Meth. Appl. Mech. Engrg., 59, 307–325, 1986.

[14] T.J.R. Hughes and M. Mallet, ``A new finite element formulations for computational fluid dynamics: III. The generalized streamline operator for multidimensional advective-diffusive systems´´, Comput Methods Appl. Mech. Engrg., 58, pp. 305–328, 1986.

[15] T.J.R. Hughes, L.P. Franca and M. Balestra, ``A new finite element formulation for computational fluid dynamics. V Circumventing the Babuska-Brezzi condition: A stable Petrov-Galerkin formulation of the Stokes problem accomodating equal order interpolations´´, Comput. Methods Appl. Mech. Engrg., 59, 85–99, 1986.

[16] L.P. Franca and S.L. Frey, ``Stabilized finite element methods: II. The incompressible Navier-Stokes equations´´, Comput. Method Appl. Mech. Engrg., Vol. 99, pp. 209–233, 1992.

[17] T.J.R. Hughes, G. Hauke and K. Jansen, ``Stabilized finite element methods in fluids: Inspirations, origins, status and recent developments´´, in: Recent Developments in Finite Element Analysis. A Book Dedicated to Robert L. Taylor, T.J.R. Hughes, E. Oñate and O.C. Zienkiewicz (Eds.), (International Center for Numerical Methods in Engineering, Barcelona, Spain, pp. 272–292, 1994.

[18] M.A. Cruchaga and E. Oñate, ``A finite element formulation for incompressible flow problems using a generalized streamline operator´´, Comput. Methods in Appl. Mech. Engrg., 143, 49–67, 1997.

[19] M.A. Cruchaga and E. Oñate, ``A generalized streamline finite element approach for the analysis of incompressible flow problems including moving surfaces´´, Comput. Methods in Appl. Mech. Engrg., 173, 241–255, 1999.

[20] T.J.R. Hughes, L.P. Franca and G.M. Hulbert, ``A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective-diffusive equations´´, Comput. Methods Appl. Mech. Engrg., 73, pp. 173–189, 1989.

[21] T.E. Tezduyar, S. Mittal, S.E. Ray and R. Shih, ``Incompressible flow computations with stabilized bilinear and linear equal order interpolation velocity–pressure elements´´, Comput. Methods Appl. Mech. Engrg., 95, 221–242, 1992.

[22] J. Donea, ``A Taylor-Galerkin method for convective transport problems´´, Int. J. Num. Meth. Engng., 20, 101–119, 1984.

[23] O. Pironneau, ``On the transport-diffusion algorithm and its applications to the Navier-Stokes equations´´, Numer. Math., 38, 309, 1982.

[24] R. Löhner, K. Morgan, O.C. Zienkiewicz, ``The solution of non-linear hyperbolic equation systems by the finite element method´´, Int. J. Num. Meth. in Fluids, 4, 1043, 1984.

[25] R. Codina, M. Vazquez and O.C. Zienkiewicz, ``A general algorithm for compressible and incompressible flow - Part III. The semi-implicit form´´ Int. J. Num. Meth. in Fluids, 27, 13–32, 1998.

[26] R. Codina and O.C. Zienkiewicz, ``CBS versus GLS stabilization of the incompressible Navier-Stokes equations and the role of the time step as stabilization parameter´´, Communications Num. Meth. Engng., 18 (2), 99–112, 2002.

[27] R. Codina and J. Blasco, ``Stabilized finite element method for the transient Navier-Stokes equations based on a pressure gradient operator´´, Comput. Methods in Appl. Mech. Engrg., 182, 277–301, 2000.

[28] T.J.R. Hughes, ``Multiscale phenomena: Green functions, subgrid scale models, bubbles and the origins of stabilized methods´´, Comput. Methods Appl. Mech. Engrg, Vol. 127, pp. 387–401, 1995.

[29] F. Brezzi, L.P. Franca, T.J.R. Hughes and A. Russo, ``${\displaystyle b=\int g}$´´, Comput. Methods Appl. Mech. Engrg., 145, 329–339, 1997.

[30] R. Codina, ``Stabilized finite element approximation of transient incompressible flows using orthogonal subscales´´, Comput. Methods Appl. Mech. Engrg., 191, 4295–4321, 2002.

[31] J. Donea and A. Huerta, ``Finite element method for flow problems´´, J. Wiley, 2003.

[32] E. Oñate, ``Derivation of stabilized equations for advective-diffusive transport and fluid flow problems´´, Comput. Meth. Appl. Mech. Engng., Vol. 151, pp. 233–267, (1998).

[33] E. Oñate, ``Possibilities of finite calculus in computational mechanics´´, Int. J. Num. Meth. Engng., 60(1), 255–281, 2004.

[34] E. Oñate, J. García and S. Idelsohn, ``Computation of the stabilization parameter for the finite element solution of advective-diffusive problems´´, Int. J. Num. Meth. Fluids, Vol. 25, pp. 1385–1407, (1997).

[35] E. Oñate, J. García and S. Idelsohn, ``An alpha-adaptive approach for stabilized finite element solution of advective-diffusive problems with sharp gradients´´, New Adv. in Adaptive Comp. Met. in Mech., P. Ladeveze and J.T. Oden (Eds.), Elsevier, (1998).

[36] E. Oñate and M. Manzan, ``A general procedure for deriving stabilized space-time finite element methods for advective-diffusive problems´´, Int. J. Num. Meth. Fluids, 31, 203–221, 1999.

[37] E. Oñate and M. Manzan, ``Stabilization techniques for finite element analysis of convection diffusion problems´´, in Computational Analysis of Heat Transfer, G. Comini and B. Sunden (Eds.), WIT Press, Southampton, 2000.

[38] E. Oñate and S. Idelsohn, A mesh free finite point method for advective-diffusive transport and fluid flow problems, Computational Mechanics, 21, 283–292, 1998.

[39] E. Oñate, F. Zárate and S.R. Idelsohn, ``Finite element formulation for convective-diffusive problems with sharp gradients using finite calculus´´, Comput. Meth. Appl. Mech. Engng., 2005.

[40] E. Oñate, J. Miquel and G. Hauke, ``Stabilized formulation for the advection-diffusion-absorption equation using finite calculus and linear finite elements´´, Comput. Meth. Appl. Mech. Engng., 2005.

[41] E. Oñate, J. Miquel and F. Zárate, ``Stabilized solution of the multidimensional advection-diffusion-absorption equation using linear finite elements´´, Comput. Meth. Appl. Mech. Engng., 2005.

[42] C.A. Felippa and E. Oñate, ``Nodally exact Ritz discretizations of 1D diffusion-absorption and Helmholtz equations by variational FIC and modified equation methods´´, Computational Mechanics, accepted for publication, 2005.

[43] E. Oñate, ``A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation´´, Comp. Meth. Appl. Mech. Engng., 182, 1–2, 355–370, 2000.

[44] E. Oñate, C. Sacco and S.R. Idelsohn, ``A finite point method for incompressible flow problems´´, Computing and Visualization in Science, 2, 67–75, 2000.

[45] E. Oñate and J. García, ``A finite element method for fluid-structure interaction with surface waves using a finite calculus formulation´´, Comput. MethodsAppl. Mech. Engrg., 191, 635–660, (2001).

[46] J. García and E. Oñate, ``An unstructured finite element solver for ship hydrodynamic problems´´, Journal of Applied Mechanics, 70, 18–26, January 2003.

[47] S.R. Idelsohn, E. Oñate and F. Del Pin, ``A lagrangian meshless finite element method applied to fluid-structure interaction problems´´, Computer and Structures, 81, 655–671, 2003.

[48] E. Oñate, S.R. Idelsohn, F. Del Pin and R. Aubry, ``The particle finite element method. An overview´´, Int. J. Comput. Meth., 1 (2), 267–307, 2004.

[49] S.R. Idelsohn, E. Oñate and F. Del Pin, ``The Particle Finite Element Method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves´´, Int. J. Num. Meth. Engng., 61, 964–989, 2004.

[50] E. Oñate, J. García and S.R. Idelsohn, ``Ship hydrodynamics´´, in Encyclopedia of Computational Mechanics, E. Stein, R. de Borst and T.J.R. Hughes (Eds), J. Wiley Sons Ltd, Vol. 3, Chapter 18, pp. 579–607, 2004.

[51] E. Oñate, J. García, S.R. Idelsohn and F. Del Pin, ``FIC formulations for finite element analysis of incompressible flows. Eulerian, ALE and Lagrangian approaches´´, Comput. Meth. Appl. Mech. Engng., 2006.

[52] E. Oñate, J. Rojek, R.L. Taylor and O.C. Zienkiewicz, ``Finite calculus formulation for incompressible solids using linear triangles and tetrahedra´´, Int. J. Num. Meth. Engng., 59, 1473–1500, 2004.

[53] E. Oñate, R.L. Taylor, O.C. Zienkiewicz and J. Rojek, ``A residual correction method based on finite calculus´´, Engineering Computations, 20, (5/6), 629–658, 2003.

[54] T.E. Tezduyar and Y. Osawa, ``Finite element stabilization parameters computed from element matrices and vectors´´, Computer Methods in Applied Mechanics and Engineering, 190, 411–430, 2000.

[55] T.E. Tezduyar, "Computation of moving boundaries and interfaces and stabilization parameters", International Journal for Numerical Methods in Fluids, 43 (2003) 555-575.

[56] T.J.R. Hughes and M. Mallet, ``A new finite element formulations for computational fluid dynamics: IV. A discontinuity capturing operator for multidimensional advective-diffusive system, Comput. Methods Appl. Mech. Engrg., 58, 329–336, 1986.

[57] R. Codina, ``A discontinuity-capturing crosswind dissipation for the finite element solution of the convection-diffusion equation´´, Comput. Methods Appl. Mech. Engrg., 110, 325–342, 1993.

[58] Tdyn. A finite element code for fluid-dynamic analysis, COMPASS Ingeniería y Sistemas SA, www.compassis.com, 2005.

[59] J.J. Kim, ``Investigation of separation and reattachment of turbulent shear layer: flow over a backward facing step´´. PhD Thesis, Standford Univeristy, 1978.

[60] N.N. Mansour and T. Morel, In AFOSR-HTTM Stanford Conference on Complex Turbulent Flows III, Stanford University, CA, pp. 1418–1423, 1981.

[61] A. Pollard, In AFOSR-HTTM Stanford Conference on Complex Turbulent Flows III, Stanford University, CA, pp. 1486–1490, 1981.

[62] A.M. Abdelmengui, S.Y. Goh, J. Ilegbusi and D.B. Spalding, ``Predictions of complex turbulent flows using the PHOENICS computer code´´, In AFOSR-HTTM Stanford Conference on Complex Turbulent Flows III, Stanford University, CA, pp. 1521–1526, 1981.

[63] F. Illinca, J.-F. Hetu and D. Pelletier, ``A Unified Finite Element algorithm for two-equations models of turbulence´´. Computers and Fluids, Vol. 27, No. 3, 291–310, 1998.

[64] V. Kalro and T. Tezduyar, ``Parallel 3D computation of unsteady flows around circular cylinders´´, Parallel Computing, 23, 1235–1248, 1997.