|
|
A FINITE ELEMENT FORMULATION FOR THE RESOLUTION
OF THE UNSTEADY INCOMPRESSIBLE VISCOUS FLOW FOR LOW REYNOLDS NUMBERS.
*P. VELLANDO, **J. PUERTAS, **J.
BONILLO, *J. FE.
ETS. de Ingenieros de Caminos,
Canales y Puertos.*Dpto. de Métodos Matemáticos y de Representación. **Dpto. de
Tecnología de la Construcción.
Universidad de La Coruña. Campus
de Elviña. 15192 La Coruña, Spain.
Tel: (34) 981 167000, fax:
(34) 981 167170. E-mail: vellando@iccp.udc.es
ABSTRACT
The following paper shows a Finite Element formulation
for the resolution of the -local and convective acceleration terms including-
Navier-Stokes equations, which gives analytical response to the problem of
viscous, incompressible, unsteady flows. The integration of the resulting
non-linear system of first order ordinary differential equations, is made upon
a successive approximation algorithm together with an implicit backward time
integrating scheme. The interpolation of the spatial domain is made in terms of
a Q1/P0 pair (bilinear velocity-constant pressure). The usage of a Bubnov
Galerkin formulation in the process of obtaining a weak form implies that flows
of a certain velocity need the employment of a very refined spatial mesh so as
to avoid numerical instability. For high Reynolds numbers the convection term
becomes predominant compared to the diffussion term and a different algorithm
(SPGU, GLS), should be introduced. Finally the developed program is checked
over some of the most commonly used flow tests and its results on velocity and
pressure are shown.
GOVERNING EQUATIONS
The unsteady incompressible Navier-Stokes equations
are given by:
in
in
(1.1)
for 
, (where
is a specified time), together with the initial and boundary
conditions:
in
on
(1.2)
where u is
the velocity, p is the pressure,
t is
the time, 
is the cinematic viscosity, f is the body force per unit
mass and 
are the gradient and laplacian tensor operators. The problem
consists in finding 
and 
in a time-space domain 
x(0,T).
FINITE ELEMENT FORMULATION
Let V be the subspace of
D(
) satisfying the incompressibility constraint: 
, let H be the closure of V in 
. If 
is sufficiently well
behaved then 
.
If we apply the weighted-residual argument by taking
the scalar product of the momentum equation times an arbitrary test function v satisfying 
, that is for 
,
(2.1)
Applying the divergence theorem and taking into
account that 
in 
, we arrive to:
(2.2)
for all admissible functions 
.
Even with 
in 
, we will maintain the pressure term so as to be able to use
a mixed formulation rather than the penalized one. This allows us to keep as
variables the pressure unknowns.Although it produces a certain increase in
computational cost, the penalty parameter formulation is known to be the cause
of loss of accuracy for small values and for holding up the convergence of the
solution for too large ones.
In the same way, making the scalar product of the
continuity equation with an arbitrary test function 
, the weak expression turns into:
(2.3)
The velocity and pressure fields may be approximated
now as piecewise polynomials on the discretization by (Q1/P0) Lagrange-Type
elements, so they are expressed in terms of the trial functions, 
,
.
(2.4)
Replacing (2.4) into (2.2) and (2.3):
(2.5)
(2.6)
For all admissible test functions 
and 
.
Introducing the expansions for 
, and integrating, the following non-linear ordinary
differential equations system is obtained.
(2.7)
In expanded matrix form:
(2.8)
where:
(2.9)
In order to turn equation (2.8) into
a linear system of equations we are going to approximate the non-linear term 
, by an iterative scheme known as successive-approximation:
(2.10)
Thus (2.7) is written:
(2.11)
Where
is:
(2.12)
Finally, the time integration based on an implicit
backward scheme, results in the following non-differential linear system:
(2.13)
Which characterises the unsteady viscous
incompressible flow.
DISCRETIZATION
The Q1/P0 pair means a
quadratic first order approximation for the velocity field and a constant
pressure approximation for each basic element. The shape functions for the
velocity field Ni are expressed in terms of the local coordinates 
.
(3.1)
The derivatives of the form functions with respect to
the global axis coordinates must be calculated so as to be able to constitute
the basic matrices, the change of coordinates leads to:
(3.2)
Combining the former expressions, the matrices M, A, Bx,
By, C may be now
expressed in terms of the derivatives with respect to the local axis
coordinates. So, the matrix Bx (for instance) being:
(3.3)
The integration of expression (3.3) and that of the
other matrices will be carried out, following a numerical 4-point
two-dimensional Gauss integration, where the Gauss points are in 
, and the weighting function H=1.
(3.4)
NUMERICAL EXAMPLES
A FORTRAN code has been developed making use of the
previously explained formulation. The program has been checked with the well
known backward step and driven cavity flow tests and the results for velocity
and pressure unknowns are shown bellow
The following figure shows the evolution in time of
the behaviour of the streamlines for the backward step problem in which the
eddy takes form for flow decreasing boundary conditions.
|
|
|
|
|
|
Fig 1. Backward step
streamlines varying in time.
The results of the velocity field, both for the
backward step and cavity flow are plotted bellow
|
|
|
Fig 2. Velocity fields.
Pressure
field numerical results in iso-lines and 3D graph
|
|
|
Fig 3 Pressure fields for the
cavity flow problem.
Fig 4 Pressure field for the
backward (16x60) step flow
CONCLUSIONS AND FURTHER
DEVELOPMENTS.
The results may be considered accurate enough for
small Reynolds numbers, with its magnitude being a function of the mesh grade
of refinement. Once the Reynolds
number becomes of considerable magnitude (and so the convective term gets
bigger compared with the diffusive term), the upwinding Petrov Galerking
formulation together with the inclusion of an artificial streamline diffusion
term (SUPG), should provide a computationaly speaking inexpensive and
convenient approach. Next, a turbulence-considering scheme and a pattern being
able to model transitions in flow, would be implemented.
REFERENCES
[1] I Christie, O. F. Griffith, A. R. Mitchel, O. C. Zienkiewicz. Finite
element methods for second-order differential equations with significant first
order derivatives. Int. J. Numer. Methods Engrg. 10 (1976) 1389-1396.
[2] A. N. Brooks, T. J. R. Hughes,
Streamline upwind/Petrov-Galerkin formulation for convection dominated flows
with particular emphasis on the incompressible Navier-Stokes equation, Comput.
Methods Appl. Mech. Engrg. 32 (1982) 199-259.
[3] G. F. Carey, J. T. Oden, Finite
Elements, Prentice Hall, (1983).
[4] , T. J. R. Hughes, M. Mallet, A.
Mizukami, A new finite element formulation for computational fluid dynamics: II
Beyod SUPG, Comput. Methods Appl. Mech. Engrg. 54 (1986) 341-355.
[5] M. D. Gunzburger, "Finite Element Methods for
Viscous Incompressible Flows, AcademicPress,(1989).
[6] 0. Pironneau, "FEM for Fluids". Willey, 1989.
[7] , T. J. R. Hughes, L. P. Franca,
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 (1989) 173-189.
[8] O. C. Zienkiewicz, R. L. Taylor,
The Finite Element Method, Mc Graw Hill, (1991).
[9] L. P. Franca, S. L. Frey
Stabilized finite element methods: II The incompressible Navier-Stokes
equation. Comput. Methods Appl. Mech. Engrg. 99 (1992) 209-233.
[10] S. K. Hannani, M. Stanislas, P.
Dupont, Incompressible Navier-Stokes computations with SUPG and GLS
formulations-A comparison study. Comput. Methods Appl. Mech. Engrg. 124 (1995)
153-170.
[11] H. G. Choi, H. Choi, J. Y. Yoo.
A fractional four-step finite element formulation of the unsteady
incompressible Navier Stokes equations using SUPG and linear equal-order
element methods. Comput. Methods Appl. Mech. Engrg. 143 (1997) 333-348.
