An alternative finite element formulation for determination of streamlines in two-dimensional problems

An Alternative Finite Element Formulation for Determination of Streamlines in Two-Dimensional Problems

Sérgio Frey

Thermal Sciences and Energy Systems Group (GESTE), Mechanical Engineering Department

Universidade Federal do Rio Grande do Sul, Rua Sarmento Leite, 425 -- 90050-170 Porto Alegre/RS, Brazil

[e-mail] - http://www.mecanica.ufrgs.br

Maria Laura Martins-Costa

Laboratory of Theoretical and Applied Mechanics (LMTA), Mechanical Engineering Department

Universidade Federal Fluminense, Rua Passo da Pátria, 156 -- 24210-240 Niterói/RJ, Brazil

[e-mail]

José Henrique Carneiro de Araujo

Laboratory of Theoretical and Applied Mechanics (LMTA) , Computational Sciences Department

Universidade Federal Fluminense, Rua Passo da Pátria, 156 -- 24210-240 Niterói/RJ, Brazil

[e-mail]

It is well known that the numerical solutions of incompressible viscous flows are of great importance in Fluid Dynamics. The graphics output capabilities of their computational codes have revolutionized the communication of ideas to the non-specialist public. In general those codes include, in their hydrodynamic features, the visualization of flow streamlines - essentially a form of contour plot showing the line patterns of the flow - and the magnitudes and orientations of their velocity vectors. However, the standard finite element formulation to compute streamlines suffers from the disadvantage of requiring the determination of boundary integrals, leading to cumbersome implementations at the construction of the finite element code. In this article, we introduce an efficient way - via an alternative variational formulation - to determine the streamlines for fluid flows, which does not need the computation of contour integrals. In order to illustrate the good performance of the alternative formulation proposed, we capture the streamlines of three viscous models: Stokes, Navier-Stokes and Viscoelastic flows.

Keywords: Streamlines, Variational Formulation, Finite Element Method, Navier-Stokes Flows, Viscoelastic Liquids.

Introduction

It is well known that the numerical solutions of incompressible viscous flows, even in two-dimensional cases, are of great importance in Fluid Dynamics. As in the pre-processing phase of these solutions, a huge amount of development work has recently taken place in their graphical post-processing, in order to explore all the potential of the numerical approximations of the fluid flows. The graphics output capabilities of the codes of Computational Fluid Dynamics - which in general include visualization of the flow streamlines - have revolutionized the communication of ideas to the non-specialist public. Streamlines - essentially a form of contour plot showing lines of the flow - are better than velocity vectors (the other standard way of depicting fluid flows) at representing the flow direction, although they tell little about how fast the flow is moving.

However the standard finite element formulation to compute these streamlines suffers from the disadvantage of requiring the determination of boundary integrals. In this article, we introduce an efficient way to determine the streamlines for such flows without the need of computing line integrals. That computation leads to either the implementation of two distinct mappings in the finite element code or to the use of additional pointer arrays to identify the nodes on the domain contour and to calculate unit tangent vectors on those nodes. It is worth to note that the former alternative - namely, the implementation of one mapping for the two-dimensional elements in the interior domain and other for the one-dimensional elements in the boundary domain - compared with the latter one, allows saving some computational time. Besides, adequate numerical techniques must be implemented to compute these integrals (Hughes, 1987).

Mechanical Modeling

A mechanical body B is a set whose elements x are denoted particles or material particles. A one-to-one continuous mapping of this set onto a region of the Euclidean space e is called a configuration of the body. The point c=X(x ) of e is the place occupied by the particle x , and x = X ^-1 (c), the particle whose place in e is x. The configuration of a body may be equivalently described in terms of the position vector x of the point x Î e with respect to the origin O: x = c (x ). A motion of a body is a one-parameter family of configurations,

(1)

where t represents the time. The streamlines for a given time t form a family of curves to which the velocity field is everywhere tangent at time t, so that their parametric equations are solutions of following system of equations,

(2)

or

(3)

with a representing an arbitrary parameter measured along the curves, e_ijk the alternator symbol (Germain and Muller, 1986) and the time t a constant. In two-dimensional flows the streamlines represent curves along which the stream-function is a constant.

Initially, we must introduce the stream-function j which aims to describe the streamlines and velocity scales at representative points in the flow. It is obtained by relating the streamlines to the principle of mass conservation. By a two-dimensional motion, we mean here one such that in some coordinate system the velocity field has only two nonzero components. Let us further restrict our discussion to incompressible fluids, so that the equation of continuity reduces to (Landau and Lifchitz, 1971),

(4)

where u is a velocity field and W a bounded open of R² with a regular boundary G . Consider a motion such that in rectangular cartesian coordinates,

(5)

so we may write eq.(4) as

(6)

Assuming the required regularity for the function j with j Î C² (W ), where C² (W ) stands for the space of continuous functions with continuous first and second derivatives

(7)

comparing eqs.(6) and (7) we see that we may define a stream-function j such that

(8)

Equation (8) gives rise to

(9)

The advantage of employing a stream-function j is that in this way the equation of continuity is identically satisfied for the flow described by eq.(5). On the other hand, on the boundary G we have

(10)

with n and t standing respectively for the unit outward normal and tangent vectors to the boundary G while e_i and e_j represent unit vectors of the cartesian system.

Now, combining eqs.(8), (9) and (10), we may form a boundary value problem to determine the streamlines of a flow: For a given velocity field u = ( u_x (x) , u_y (x)), find a scalar j = j (x) such that

(11)

The problem defined by eq.(11) has a solution only to within an additive constant if the following relation is verified (Ciarlet, 1978),

(12)

The main subject of the present work is to present a Galerkin formulation for the streamline problem stated in eq.(11) in such a way that undesirable integrations on the boundary G are not required. As examples, streamlines for Stokes, Navier-Stokes and Viscoelastic flows are obtained by employing a finite element method.

Variational Formulations

Assuming the required regularity (Ciarlet, 1978) for both trial function j and test function y with j Î H¹ and y Î H¹, we may construct a classical variational formulation for eq.(11) as follows,

(13)

in which Gauss divergence theorem (Germain and Muller, 1986) has been employed.

Now, according to eq.(13), we may state a symmetric variational formulation of the problem defined by eq.(11) as: Given u = (u_x (x), u_y (x)) : W ® R² , find j = j (x) Î W such that

(14)

in which the functional set W is defined by,

(15)

with H¹(W ) denoting the Sobolev space of functions with square-integrable value and first-derivatives in W , (Ciarlet, 1978)

(16)

and L²(W ) standing for the space of square-integrable functions in W ,

(17)

The variational formulation defined by eq.(14) suffers from the disadvantage of demanding the computation of boundary integrations. This shortcoming may require - in order to save computational time - the implementation of two distinct mappings in the finite element code, namely, one for interior domain W employing two-dimensional elements and another with one-dimensional elements in order to approximate the boundary domain G . Besides the implementation of two mappings - which is not mandatory - additional pointer arrays must be employed in order to identify the nodes on the domain contour and unit tangent vectors must be calculated. Furthermore, it is necessary to implement adequate numerical techniques to evaluate these integrals.

In order to overcome this computational shortcoming - leading to cumbersome implementations at the construction of the finite element code - it is convenient to present an alternative formulation to eq.(14). Noting that

(18)

with k being the z-Cartesian director, and introducing the vectorial identity (Germain and Muller, 1986) we may rewrite eq.(14) as

(19)

From Germain and Muller (1986), we may extract the further identity,

(20)

Then, combining eqs.(19) and (20), we obtain

(21)

Thus the following alternative variational formulation for the problem described in eq.(11) runs as follows: Given u = (u_x (x), u_y (x)) : W ® R² , find j = j (x) Î W such that

(22)

Finite Element Modelling

Consider a bounded domain W Ì R² with a regular boundary G and let C_h be a partition of W into elements consisting of convex quadrilaterals, which is performed in the usual way (Ciarlet, 1978),

(23)

Now it is convenient to define the following finite element spaces,

(24)

where P_k stands for polynomial space of degree k ³ 0. Based upon definition (24) and supposing that the finite element approximation j _hÎ W_h admits the representation (Hughes, 1987)

(25)

where f _hÎ W_h and j ⁰_{h = j} ⁰ on x_{0Î W} , a finite element approximation of the Galerkin variational problem defined in eq.(22) may be stated as: Given u = (u_x (x), u_y (x)) : W ® R² , find j _h = j _h (x) Î W_h, j _h (x) = f _h (x) + j ⁰_h (x) with f _h (x) Î W_h such that

(26)

Let us now express the finite element approximation as a combination of known shape functions N_A and unknown degrees-of-freedom y _A (Ciarlet, 1978),

(27)

with h -{x₀} denoting the complement of the Dirichlet prescribed node {x₀} in the set of all global nodes h Thus, h -{x₀} represents the set of nodal points at which y _h is to be determined, namely, the number of nodes h -{x₀} is equal to the number of equations. Substituting eq.(27) in the Galerkin problem expressed in eq.(26), we obtain

(28)

where

(29)

As the eqs.(28)-(29) must hold for all test functions y _h, we get that the coefficients y _A in eq.(27) may assume any arbitrary value and therefore we shall impose G_A =0, " AÎ h -{x₀}, resulting in the following matrix problem associated to the Galerkin formulation (eq.(26)): Given the stiffness matrix K and load vector F, Find the vector d such that,

(30)

where

(31)

with the array D introduced to assign to a global node A or B the corresponding global equation numbers P or Q and n_eq standing for the number of algebraic equations (the dimension of the complement h -{x₀}). In order to save computational memory, it is more convenient to handle the global matrices K and F of system (30)-(31) in terms of element matrices k^K and f^K as follows,

(32)

in which A represents the assembly operator (Hughes, 1987) to form K and F from the element arrays k^K and f^K, 1 £ a, b £ n_en with n_en being the number of element nodes, n_el the number of mesh elements and j ^K_b = j ⁰ at the prescribed global node x₀ and, otherwise, equals to zero.

Some Applications

In this section we present some two-dimensional simulations of the streamlines for viscous incompressible flows employing the alternative finite element formulation introduced in eq.(26). All computations have employed the finite element code FEM and graphics post-processor VIEW, both codes under development at Laboratory of Theoretical and Applied Mechanics (LMTA) of Fluminense Federal University.

Newtonian Incompressible Flows

Initially, two Newtonian models in fluid dynamics for the leaky cavity and the stepward channel benchmarks were considered: first the linear Stokes model and further the incompressible Navier-Stokes one, which may be described by the following set of equations,

(33)

where the pair (u, p) represents velocity and pressure fields, n the fluid kinematical viscosity, e (u) the symmetrical part of Ñ u tensor, f the body force acting on the flow and s is Cauchy stress tensor following the Newtonian hypothesis (Germain and Muller, 1986)

(34)

with I denoting the identity tensor, G _g the region of the boundary G on which essential (Dirichlet) conditions are imposed and G _h the region of G subjected to natural (Neumann) ones,

(35)

The Newtonian leaky cavity. The geometry and boundary conditions of the problem are sketched in Fig.1: the domain W is such that -0.5£ x, y £ 0.5 and represents a biunity cavity subjected to a flow on its superior boundary - expressed by Dirichlet boundary conditions,

(36)

(37)

The boundary-valued problem defined by eqs.(33)-(34), (36)-(37) was computed in Franca and Frey (1992) via a stabilized finite element method employing an equal-order serendipity elements (Q2/Q2) in order to approximate the pair velocity/pressure. The employed method dealt with the transient and non-linear features of the model via a predictor/multi-corrector algorithm based upon the general trapezoidal method (Hughes, 1987). The dynamical viscosity is n = 400^-1 for the Navier-Stokes application and n ® ¥ for the Stokes one.

Substituting the velocity field established by Franca and Frey (1992) in eq.(26), the streamlines for both Stokes and steady-state Navier-Stokes problems - sketched in Fig.2 - were obtained employing a uniform finite element mesh with 16x16 Lagragean biquadratic (Q2) elements. The Stokes problem shows symmetric streamline patterns - in relation to y-axis, while the non-symmetric streamline features for the nonlinear flow were also adequately captured in the Navier-Stokes model. A good agreement with the results presented by Gresho and Chan (1990) and Tezduyar et al. (1990) has been verified.

The Newtonian backward step flow. The problem statement and its boundary conditions are shown in Fig.3: the domain represents a flat stepwise divergent channel subjected to Dirichlet boundary conditions. At the channel entrance and exit, we have imposed parabolic flow profiles preserving the mass conservation. On its walls, the classical no-slip condition has been employed. We take Re=50 with respect to the inflow section. In order to obtain the velocity field governed by the system (33)-(34) and subjected to the above-mentioned boundary conditions, we again recur to the stabilized method introduced in Franca and Frey (1992) with a Q2/Q2 serendipity equal-order pair of elements.

So, employing the velocity field obtained by Franca and Frey (1992) for the above mentioned case in eq.(26), the streamlines for the Navier-Stokes steady-state model are obtained. They are sketched in Fig.4, for a uniform finite element mesh with 224 Lagragean biquadratic (Q2) elements. In the computations shown in this figure, we may observe a well-defined re-circulation at the step region as found in Macedo (1995), which was absent in the linear Stokes model.

A Non-Newtonian Sudden Contraction

Secondly, an incompressible steady-state non-Newtonian flow through a four to one abrupt contraction channel will be considered. This flow may be modeled by the following boundary-value problem,

(38)

where u, p, e (u) and f are defined as in eq.(33), G_g and G_h are defined as in eq.(35), s₁ is the part of the extra stress tensor s defined by (Marchal and Crochet, 1987)

(39)

h ₁ and h ₂ dynamical viscosities, l a relaxation time, n unit outward normal and the boundary G - _g is expressed by

(40)

The velocity field, solution of the system (38), was obtained from Frey et.al. (1999), employing a three-field finite element method. In a partition C_h of W defined by quadrilateral elements, the velocity field is approximated by continuous piecewise biquadratic elements and the pressure field by discontinuous piecewise linear ones. The extra-stress tensor is approximated by continuous piecewise bilinear elements in direct sum with twelve bubble tensors in each element as defined in quoted article. The geometry and boundary conditions of the problem are sketched in Fig.5: the domain is a plane channel with a sudden 4:1 contraction. The boundary conditions for velocity and s₁ tensor, at the channel entrance, are fully developed profiles while, at its exit, a fully developed profile was considered for the velocity field only - taking into account the mass conservation. In the numerical tests, the authors selected a viscosity ratio h ₂/(h ₁+h ₂)=0.11 - usual in Oldroyd-B simulations (Ruas and Carneiro de Araujo, 1992), measured the Deborah number (Marchal and Crochet, 1987) at downstream fully developed region and made use of an incremental Newton's method with a continuation strategy on the stress relaxation time (Marchal and Crochet, 1987) to solve the set of non-linear equations.

Fig.6 shows details of the flow streamlines near the contraction for Deborah number (De) equal to 1.2 and 15. As it can be seen in the figure, the alternative finite element formulation defined by eq.(26) was able to capture a well-defined back-flow pattern for the streamlines nearby the channel contraction. These results are qualitatively in accordance with experimental visualizations introduced in Evans and Walters (1986).

Conclusions

In this article, an efficient way to determine streamlines for fluid flows without the computation of contour integrals has been introduced. The alternative Galerkin formulation proposed has the computational advantage of not requiring the implementation of distinct mappings in finite element codes nor the use of additional pointer arrays to identify the nodes on the domain contour nor the cumbersome computation of these integrals. In order to illustrate the good performance of the alternative formulation proposed in eq.(26), we have captured streamlines for Stokes, Navier-Stokes and Viscoelastic flows. The method proposed was sufficiently accurate to describe even secondary flows, as back-flow re-circulations in the Navier-Stokes model - for both leaky cavity and backward step flows - and sudden contraction flows for Viscoelastic liquids.

Acknowledgments

The authors S. Frey, M.L. Martins-Costa and J.H. Carneiro de Araujo gratefully acknowledge the financial support provided, respectively, by the agency CNPq through grants n^o 350747/93-8 and 300404/91-3 and by the agency FAPERJ through grant n^o E-26/171316/97.

Manuscript received: May 1998, Technical Editor: Angela Ourívio Nieckele.

• Ciarlet, P.G., 1978, "The Finite Element Method for Elliptic Problems", North-Holland, Amsterdam.
• Evans, R.E., and Walters, K., 1986, "Flow Characteristics Associated with Abrupt Changes in Geometry in the Case of Highly Elastic Liquids", J. Non-Newtonian Fluid Mech., Vol.20, pp.11-29, 1986.
• Franca, L.P., and Frey, S., 1992, "Stabilized Finite Element Methods: II. The Incompressible Navier-Stokes Equations", Comput. Methods Appl. Mech. Engrg., Vol. 99, pp.209-233.
• Frey, S., Carneiro de Araujo, J.H. and Silva Ramos, M.A., 1999, "Finite Element Approximation for the Heat Transfer in Incompressible Oldroyd-B Flows", J. Braz. Soc. Mech. Sci., Vol. 21/1, pp.29-40.
• Germain, P. and Muller, P., 1986, "Introduction `a la Mécanique des Milieux Continus:", Masson Editions, Paris.
• Gresho, P.M. and Chan, S.T., 1990, "On Theory of Semi-Implicit Projection for Viscous Incompressible Flow and its Implementation via a Finite Element Method that Also Introduces a Nearly Consistent Mass Matrix. Part 2: Implementation", Int. J. Numer. Methods Fluids, Vol.11, pp.621-659.
• Hughes, T.J.R., 1987, "The Finite Element Method: Linear Static and Dynamic Finite Element Analysis", Prentice-Hall, Singapore.
• Landau, L. and Lifchitz, E., 1971, "Mécanique des Fluides", Mir Editions, Moscou.
• Macedo, A.P., 1995, "Application of Totally Stabilized Finite Element Methods - GLS - to the Numerical Simulation of Laminar and Turbulent Flows", (in Portuguese) Master Dissertation - Mech. Engrg. Dpt., University of Brasília, DF, Brazil.
• Marchal, J.M. and Crochet, M., 1987, "A New Mixed Finite Element for Calculating Viscoelastic Flow", J. Non-Newtonian Fluid Mech., Vol.26, pp.77-117.
• Ruas, V. and Carneiro de Araujo, J.H., 1992, "Un Método Mejorado de Segundo Orden para la Simulatión de Flujo Viscoelástico con Elementos Finitos Quadrilaterales", Rev. Intern. Métodos Numéricos para Cálculo y Diseño en Ingeniería, Vol.8/1, pp.77-85.
• Tezduyar, T.E., Shih, R., Mittal, S. and Ray, S.E., 1990, "Incompressible Flow Using Stabilized Bilinear and Linear Equal-Order-Interpolation Velocity-Pressure Elements", Preprint UMSI, No.90/165, University of Minnesota, Minneapolis, USA.

Publication Dates

• Publication in this collection
  15 Dec 2000
• Date of issue
  2000

History

• Received
  May 1998