- Email: [email protected]
Boundary domain integral method for high Reynolds viscous fluid flows in complex planar geometries
Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220 www.elsevier.com/locate/cma
Boundary domain integral method for high Reynolds viscous fluid flows in complex planar geometries M. Hribersˇek *, L. Sˇkerget Faculty of Mechanical Engineering, Laboratory of Computational Fluid Dynamics, Institute of Power, Process and Environmental Engineering, University of Maribor, Smetanova 17, SI-2000 Maribor, Slovenia Received 9 February 2004; received in revised form 5 July 2004; accepted 2 November 2004
Abstract The article presents new developments in boundary domain integral method (BDIM) for computation of viscous fluid flows, governed by the Navier–Stokes equations. The BDIM algorithm uses velocity–vorticity formulation and is based on Poisson velocity equation for flow kinematics. This results in accurate determination of boundary vorticity values, a crucial step in constructing an accurate numerical algorithm for computation of flows in complex geometries, i.e. geometries with sharp corners. The domain velocity computations are done by the segmentation technique using large segments. After solving the kinematics equation the vorticity transport equation is solved using macro-element approach. This enables the use of macro-element based diffusion–convection fundamental solution, a key factor in assuring accuracy of computations for high Reynolds value laminar flows. The versatility and accuracy of the proposed numerical algorithm is shown for several test problems, including the standard driven cavity together with the driven cavity flow in an L shaped cavity and flow in a Z shaped channel. The values of Reynolds number reach 10,000 for driven cavity and 7500 for L shaped driven cavity, whereas the Z shaped channel flow is computed up to Re = 400. The comparison of computational results shows that the developed algorithm is capable of accurate resolution of flow fields in complex geometries. 2004 Elsevier B.V. All rights reserved. Keywords: Boundary element algorithms; Incompressible viscous fluid; Navier–Stokes equations; Velocity–vorticity formulation; Segmentation technique; Driven cavity flow
*
Corresponding author. Tel.: +386 62 220 7772; fax: +386 62 220 7990. E-mail address: [email protected] (M. Hribersˇek).
0045-7825/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2004.11.002
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4197
1. Introduction The origins of the boundary element method (BEM) are strongly connected with computation of linear or weakly non-linear field problems, which can be computed by means of boundary-only discretization of computational domain. In case of flows of viscous fluids at higher values of Re number, the flow phenomena is strongly non-linear and boundary-only discretization had to be paired with domain discretization as well. In the context of BEM related methods for laminar viscous fluid flows several successful attempts have already been made, e.g. [2,5,9,10,12] among others. An excellent survey of these approaches can be found in [8]. These numerical approaches were based on different forms of Navier–Stokes equations, representing the frame for the solution of viscous flow problems. Different techniques of capturing non-linear domain effects were developed, including internal cells, macro-elements and dual and multiple reciprocity methods. Different approaches were successful also in cases of higher values of Re number flows, although the highest values, like for the case of driven cavity, computed with other established approximation methods (Re = 10,000), were not easily obtained [3] (Re = 10,000) or were restricted to lower or moderately high Re number values [10] (up to Re = 1000) and [9] (up to Re = 5000). In boundary domain integral method (BDIM), velocity–vorticity formulation was used for the solution of Navier–Stokes equations [2]. Its main advantage is an implicit computation of boundary vorticity values, which presents a solid basis for development of a successful computational technique. The use of different types of fundamental solutions, especially macro-element based diffusion–convection fundamental solution [5] led to development of flexible and accurate computational methods [3,7] able to accurately resolve also higher value Re number flows. In obtaining a divergence-free numerical solution for all types of geometry a special attention has to be paid to a proper numerical treatment of governing equations [20]. In case of BDIM, vector potential formulation [4] and vector–velocity formulation of flow kinematics [5] were already used. In the first case in flow kinetics the parabolic-diffusion fundamental solution was used [2] and in the second the diffusion–convection fundamental solution was used. In [11] it was shown that both approaches result in the equal final integral statement. In order to lower computational cost of the method subdomain technique was used with vector–potential formulation [4] although it was restricted to non-star arrangement of subdomains and to segmentation of flow kinematics. On the other hand, the velocity vector formulation of flow kinematics [5] allowed the use of macro-element based subdomain technique in both flow kinematics and flow kinetics, but it lacked the conservation properties when applied to complex geometries. In order to overcome these drawbacks the paper presents a further development of BDIM technique, combining the velocity Poisson equation for flow kinematics and segmentation of flow kinematics together with macro-element based computation of flow kinetics with the use of diffusion–convection fundamental solution. The result is a stable numerical technique, able of accurate computation of high value Re flows in complex geometries. The new computational algorithm can be viewed as a development of the BDIM method [4], where the flow kinetics is computed in a macro-element based domain discretization with the use of diffusion–convection fundamental solution. Due to this the BDIM computational algorithm has to be modified in order to combine segmentation technique of flow kinematics with macro-element approach used in flow kinetics. Additionally, the frozen integrals and matrix technique has to be developed in order to lower computational demands of flow kinetics. The paper is organized as follows: Section 2 deals with the set of governing equations in its primitive and velocity–vorticity form, followed by Section 3 describing transformation of partial differential forms of governing equations to its integral forms, representing the starting point for discretization by boundary elements and domain cells, described in Section 4. In Section 5 the use of subdomain and macro-element technique in conjunction with segmentation technique in flow kinematics is explained. Section 6 deals with computational aspects of the derived algorithm, especially with the frozen integrals and frozen matrix
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4198
approach, and Section 7 concentrates on computational results for selected test cases: driven cavity flow, flow in L shaped driven cavity and flow in a Z shaped channel. The paper closes with conclusions based on derived numerical algorithm and obtained numerical results.
2. Governing equations When developing a numerical method for solution of viscous incompressible fluid flows the Navier– Stokes set of equations serves as the starting point for derivation of the final algebraic equations, suitable for numerical solution. For the numerical solution of the Navier–Stokes equation set we choose velocity– vorticity formulation, which has several advantages when used in the context of BDIM. When using this formulation the dynamics of a viscous incompressible fluid is partitioned into its kinematic and kinetic aspect through the use of derived vector vorticity field function xi(rj, t). Vorticity is obtained as a curl of the velocity field vi(rj, t), namely xi ¼ eijk
ovk oxj
ð1Þ
and thus xi is solenoidal, oxi ¼ 0: oxi
ð2Þ
By applying the curl operator to the vorticity definition and by using the continuity constraint for the incompressible flows the following vector elliptic Poisson equation for velocity vector is obtained: o2 v i oxk þ eijk ¼ 0: oxj oxj oxj
ð3Þ
Eq. (3) represents the kinematics of an incompressible fluid motion, expressing the compatibility and restriction conditions between velocity and vorticity field functions. The kinetic aspect is governed by the parabolic diffusion–convection vorticity transport equation obtained by applying the curl operator to the momentum equation, oxi oxi oxj vi o2 x i þ vj ¼ þm : ot oxj oxj oxj oxj
ð4Þ
The present investigation deals with constant material properties, kinematic viscosity m. One of the advantages of using velocity–vorticity formulation lies in removing pressure gradient terms from the solution process, resulting in a higher numerical stability of the computational scheme. The pressure does not appear in the solution procedure and has no influence on the velocity field, a fact that is of course valid only for incompressible fluid approximation. The pressure can be computed after velocity field is being solved. For the solution of the Poisson pressure equation the Neumann conditions are specified over the boundary of the domain, but at least in one point the pressure value has to be specified. This could well be the point on the inlet or outlet, thus the computed pressure field will correspond to the specified boundary condition and computed velocity field. In the case of compressible fluid the pressure is a thermodynamic variable and has direct influence on the density of the fluid, which further influences velocity and vorticity fields. The velocity–vorticity BEM method is currently being extended to the solution of this class of problems. For the case of the 2D flow, which is considered in the paper, the vorticity vector has just one component, perpendicular to the plane of the flow, simplifying Eqs. (3) and (4) to
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4199
o2 v i ox þ eij ¼ 0; oxj oxj oxj
ð5Þ
ox ox o2 x þ vj ¼m : ot oxj oxj oxj
ð6Þ
3. Integral representations The unique advantage of the BDIM originates from the application of GreenÕs fundamental solutions as particular weighting functions. Since they only consider the linear transport phenomenon, an appropriate selection of a linear differential operator L½ is of main importance in establishing stable and accurate singular integral representations corresponding to the original differential conservation equations. All different conservation models can be written in the following general form: L½u þ b ¼ 0;
ð7Þ
where the operator L½ can be either elliptic or parabolic, u(rj, t) is an arbitrary field function, and the nonhomogeneous term b(rj, t) is applied for non-linear transport effects or pseudo body forces. In the following, integral representations will be derived using different fundamental solutions based on particular type of governing equation under consideration, i.e. elliptic Poisson partial differential equation (PDE) for flow kinematics and elliptic diffusion–convection PDE for vorticity transport. 3.1. Integral representations for flow kinematics When transforming the elliptic Poisson PDE (3) into algebraic form the first step is the choice of the linear differential operator L½, which in the case of (3) reads as L½ ¼
o2 ðÞ ; oxj oxj
ð8Þ
allowing us to rewrite Eq. (3) as L ½ v i þ bi ¼
o2 v i þ bi ¼ 0; oxj oxj
ð9Þ
where bi ¼ eijk
oxk : oxj
ð10Þ
There are two possibilities in deriving the final integral statement of flow kinematics, suitable for the discretization and consequent use in the computational BDIM algorithm. The first one was presented in [5] and allowed the use of macro-element based discretization. Despite the low computational demands of this approach, resulting in sparse flow kinematics matrices, the formulation was limited to simple geometries with prescribed no-slip boundary conditions. In [1] an excellent explanation for this behaviour can be found noting that the solution of flow kinematics equation in [5] does not necessarily satisfy both continuity constraint and conservation equation (2), although vice versa is always true. In order to exclude non-conservative solutions additional enforcement of continuity constraint in equation has to be undertaken, as was shown recently in [11]. This formulation presents the basis for the flow kinematics integral equations, used with the presented BDIM computational algorithm. The derivation of
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4200
the integral form starts with the choice of uw as the Laplace fundamental solution, application of the GreenÕs theorems for scalar functions and additional implementation of the continuity constraint and vorticity definition, the final integral form of flow kinematics is obtained [11] which for the planar case reads as Z Z Z ouH ouH ouH dC ¼ dC eij x vj dX ð11Þ cðnÞvi ðnÞ þ vi on os oxj C C X and for the general 3D case Z Z Z ~ H ~ ~ H ~ ~ H Þ dX: cðnÞ~ vðnÞ þ ðru nÞ~ v dC ¼ ðru nÞ ~ v dC þ ð~ x ru C
C
ð12Þ
X
Eqs. (11) and (12) represent the starting point for discretization by means of boundary elements and domain cells. It can be used for evaluating unknown velocities in the interior of computational domain, whereas for determination of values on the boundary of computational domain (mainly vorticities) the discretization leads to the singular system matrix [2]. Therefore for the boundary, when the unknowns are the boundary vorticity values or the tangential velocity component values, one has to use the tangential component of the vector, Eq. (12), namely Z Z Z H H ~ ~ ~ H Þ dX cðnÞ~ nðnÞ ~ vðnÞ þ ~ nðnÞ ðru ~ nÞ~ v dC ¼ ~ nðnÞ ðru ~ nÞ ~ v dC þ ~ nðnÞ ð~ x ru C
C
X
ð13Þ in order to obtain the appropriate non-singular implicit system of equations. When the normal velocity components to the boundary are unknown, the normal form has to be employed: Z Z Z ~ H ~ ~ H ~ ~ H Þ dX: cðnÞ~ nðnÞ ~ vðnÞ þ ~ nðnÞ ðru nÞ~ v dC ¼ ~ nðnÞ ðru nÞ ~ v dC þ ~ nðnÞ ð~ x ru ð14Þ C
C
X
3.2. Integral representations for flow kinetics We discretize
du dt
in Eq. (4) by the backward Euler scheme approximation, given by
ou uF uF 1 ; ot Dt
ð15Þ
where Dt is the time step and F denotes the time step index, the linear differential operator can be written in the following diffusion-convective non-homogeneous PDE form, namely L½ ¼ m
o2 ðÞ oðÞ ðÞ vj ; oxj oxj oxj Dt
ð16Þ
where the velocity field vj was decomposed into an average constant vector vj and a perturbed one ^vj , i.e. vj ¼ vj þ ^vj . The following can now be stated: L½xi þ bi ¼ m
o2 xi oxi xi vj þ bi oxj oxj oxj Dt
ð17Þ
with bi ¼ ^vj
oxi oxj vi xi;F 1 : þ þ oxj oxj Dt
The second term on the right-hand side is zero in case of 2D flows considered in the present paper.
ð18Þ
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4201
As fundamental solution u* the diffusion–convection fundamental solution is implemented, being the solution of equation m
o2 u ou 1 u þ dðn; sÞ ¼ 0; þ vj oxj oxj oxj Dt
ð19Þ
given for the planar case as 1 ðvj rj Þ K 0 ðlrÞ exp ; 2p 2m 2 v 1 l2 ¼ ; þ 2m m Dt u ¼
ð20Þ
ð21Þ
where K0 is modified Bessel functions of the second kind, v is the magnitude of the average part of the velocity vector, s is the field point and n the source point. By applying the GreenÕs theorem [5] the following integral statement is obtained: Z Z Z Z Z ou ox 1 1 ou 1 u dC dC ¼ cðnÞxðnÞ þ x xvn u dC þ x^vj dX þ xF 1 u dX ð22Þ m C m X mDt X on oxj C C on with vn = vjnj. The velocity decomposition is necessary since the fundamental solution of Eq. (16) is known only for the linear diffusion-convective PDE with constant coefficients, i.e. known average velocity, the result of the evaluation of flow kinematics equations in the current BDIM iteration step.
4. Algebraic representations of integral equations The set of integral forms of governing equations, Eqs. (11) and (22), is transformed to its algebraic form first by dividing boundary C into boundary elements Ce and domain X into domain cells Xc. Second, geometry, the field functions vi, x and their derivatives are approximated by the use of quadratic interpolation functions [5]. Due to the use of macro-element approach in discretization of flow kinetics equations different positions of collocation nodes on boundary elements for flow kinematics and flow kinetics have to be selected, leading to two types of implemented boundary elements, Fig. 1: • Continuous quadratic elements: collocation nodes coincide with geometry nodes, this discretization is used for boundary elements in flow kinematics. • Discontinuous quadratic elements: collocation nodes are shifted from corners into the interior of the boundary element, thus avoiding the problem of obtaining an overdetermined system of equations in case of macro-element discretization approach [6] where these elements are used.
Fig. 1. Discrete model: continuous and discontinuous quadratic boundary elements (left), continuous quadratic domain cell (right).
4202
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
The detailed description of the implemented macro-element approach is given in Section 5. The use of continuous boundary elements in flow kinematics discretization ensures continuity of boundary vorticity values at sharp corners, what is of great advantage over the previous BDIM algorithm [5] where discontinuous boundary elements were used for both, flow kinematics and flow kinetics computations. After collocating over all nodes one obtains the final system of algebraic equations. In case of flow kinematics the matrix form of the final system of equations reads as ½H e fvx g ¼ ½H t fvy g ½Dy fxg;
ð23Þ
½H e fvy g ¼ ½H t fvx g þ ½Dx fxg:
ð24Þ
Forms (23) and (24) are discretized forms of Eq. (11) and are used for computation of velocity values for domain nodes. For computation of unknowns on the boundaries of computational domain equations (13) and (14) have to be used. After assembling equations written for all source points, the final matrix statements for tangential and normal forms read as ð½H 1 þ ½H t1 Þfvt g ¼ ð½H 2 ½H t2 Þfvn g þ ½Dn fxg;
ð25Þ
ð½H 1 þ ½H t1 Þfvn g ¼ ð½H t2 ½H 2 Þfvt g þ ½Dt fxg:
ð26Þ
In the case of vorticity kinetics the final matrix form of equations is obtained after discretization of Eq. (22) and reads as ox 1 1 ½DfxgF 1 : ½H fxg ¼ ½G ð27Þ þ ð½Dj ½C j Þ½^vj fxg þ on m mDt The elements of matrix [Cj] are computed using continuous boundary elements, which is necessary since both [Dj] and [Cj] are the result of transformation of the first term on the right-hand side bi of flow kinetics, Eq. (18), into the boundary and domain integral parts. In Eqs. (23)–(27), all values of field functions correspond to the actual time step (denoted also as F) and (F 1) denotes the previous time step values.
5. Subdomain and macro-element technique One of the main drawbacks of BEMs were always large matrices, that one had to store in-core or out-ofcore and corresponding large computing times, partly due to computationally expensive integration approaches and partly due to large computing times needed for solution of resulting sets of linear equations. In both cases, considerable saving in computational time and memory requirements are possible by applying subdomain technique to original set of integral equations. The original idea of the subdomain technique is to decompose computational domain into a few relatively large non-overlapping subdomains, interconnected by suitable interface conditions [8]. In comparison with other solution techniques the subdomain technique could be denoted as a direct Domain Decomposition Method (DDM), similar to Schur complement type DDM method [13]. In the context of BDIM this was successfully implemented in [4]. From the numerical approach point of view the technique could be applied to any kind of integral equation however it was shown [4] that in the case of flow kinematics computation it does seriously affect convergence of the coupled velocity–vorticity iteration scheme. This is attributed to the loss of interconnectivity between computational nodes, as the influence of fundamental solution reaches far beyond one domain cell or boundary element. Boundary unknowns, mainly unknown vorticity values on solid boundaries, must therefore be computed from the original one-domain computational mesh.
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4203
Since domain values of velocities are computed explicitly from Eqs. (23) and (24) it is possible to perform computation of these values in a segmentation mesh. The segmentation mesh is composed of several segments, which are geometrically identical to subdomains. However, segmentation technique has to be distinguished from the subdomain technique since segments are independent of each other and no interface conditions have to be imposed. In order to get accurate results the boundary values on segment boundaries have to be computed explicitly from the one-domain mesh. This is the main idea of segmentation technique, presented in [4]. In [4], the unknowns in vorticity transport equation were computed using standard subdomain technique, therefore the computational mesh used for segmentation technique could be used also for the assembly of vorticity transport matrices of integrals. The fundamental solution, used with this approach did not account for any field function values and the numerical scheme was not very stable for higher values of Reynolds number. When the subdomain technique is applied in its limit version, i.e. one subdomain is equal to one internal cell with corresponding boundary elements, the macro-element approach is derived. A macro-element is in general a small subdomain and requires application of implicit interface (I) conditions between neighboring macro-elements (for example 1 and 2), i.e. 1
2
xjI ¼ xjI ; 1 2 ox ox m ¼ m ; on I on I
ð28Þ
in the assembly of the final system of equations. The macro-element approach was used with the flow kinematics velocity–vector form in the BDIM algorithm [5] and tested for various test examples with Dirichlet boundary conditions. The approach was unfortunately limited to regular geometries of computational domains. On the other hand, application of cell or macro-element based diffusion–convection fundamental solution, developed in [5], showed that the numerical stability and accuracy of computational algorithm can be substantially increased with its implementation. Additionally, it is very important that the second main advantage of the macro-element discretization are sparse matrices of integrals and sparse system matrix, generally containing 25 entries in a row. In the computational scheme of flow kinetics (vorticity and heat transport) in [11] all the nodes (domain and boundary) had to be included into the system matrices in order to be able to compute higher Re number value flows, resulting in large fully populated system matrices. When compared with standard fully populated BEM system matrices, the macro-element discretization of flow kinetics is therefore clearly the best possible solution in the construction of a fast computational algorithm. The present paper presents a novel numerical approach to the BDIM which combines segmentation technique for flow kinematics with macro-element approach for flow kinetics i.e. vorticity transport. In this way the algorithm combines conservation properties of vector–potential formulation of flow kinematics and accurate computation of diffusion–convection transport phenomena, performed by macro-element based computation of flow kinetics. The new numerical approach is composed of the following steps: 1. Computation of boundary vorticity and velocity values from flow kinematics equation, discretized by the continuous full domain discretization, Fig. 2 (left), presented here in a final matrix form of system of equations obtained after implementation of boundary conditions and rearrangement of Eqs. (26) and (25): ½AC fxgC ¼ bCþX :
ð29Þ
4204
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
Fig. 2. Sequence of computational meshes.
2. Computation of velocity values on segment boundaries: fvx gS i ¼ ½H t fvy gC ½Dy fxgCþX ;
ð30Þ
fvy gS i ¼ ½H t fvx gC þ ½Dx fxgCþX :
ð31Þ
3. Computation of velocity values inside the segments, Fig. 2 (middle): fvx gXi ¼ ½H t fvy gCi ½Dy fxgCi þXi ;
ð32Þ
fvy gXi ¼ ½H t fvx gCi þ ½Dx fxgCi þXi :
ð33Þ
4. Transformation of velocity and vorticity values from segmentation mesh (Xi) to macro-element mesh (XME): p1 : fvi gXi ! fvi gXME ;
ð34Þ
p2 : fxgC ! fxgCME :
ð35Þ
5. Computation of domain vorticity values and vorticity fluxes from macro-element mesh, presented here in a final matrix form of system of equations obtained after implementation of boundary and interface conditions and rearrangement of Eqs. (27): ½ACME þXME fxgCME þXME ¼ bCME þXME :
ð36Þ
6. Transformation of vorticity values from macro-element mesh (XME) to segmentation mesh (Xi): r : fxgXME ! fxgXi :
ð37Þ
6. Computational aspects of the method 6.1. Inner iteration loop The kinematics and kinetics of the fluid flow motion can be seen as two interlaced problems in a computational procedure. For the general time dependent flows two iteration loops are necessary in order to obtain a solution. The first loop, the time advancing loop or the outer loop, is using known information about the motion at an instant time level for computation of values at subsequent time level, here achieved by the approximation (15). The other iteration loop, the inner loop, serves as a means of computing a divergence free velocity and vorticity field values. In the BDIM context, results of computation of flow kinetics (36) are values of vorticities inside the computational domain, and the result of computation of flow kinematics (29) are vorticity
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4205
values at the boundary of computational domain. Vorticity solution in the whole domain must therefore be searched in an iterative manner, where, with the use of under-relaxation parameter (H), new domain vorticities are computed: iþ1
fxgXME ¼ H iþ1 fxgXME þ ð1 HÞ i fxgXME :
ð38Þ
The inner loop is stopped after the following convergence check for a predefined error (usually 106, Nc = number of computational nodes) is satisfied: PN c iþ1 2 xj i xj Þ j¼1 ð ð39Þ C¼ PN c iþ1 2 : xj Þ j¼1 ð These two iteration loops constitute a computational loop which advances the solution from initial time level to a final time level. 6.2. Recomputation of matrices of integrals of flow kinetics Since the value of velocity vector is changing during each inner iteration formally this would also mean that the entries in matrices of integrals of kinetics must change, as fundamental solution (20) includes the average part of the velocity vector, vj , in a macro-element. This would strongly increase the computational cost of the method therefore the average part of velocity vector is frozen until the relative change in its magnitude has reached a certain level. In practice this means the average velocity in fundamental solution (20) and hence the integrals of a macro-element are updated after the criterion jvj jvj ð40Þ jvj P int is satisfied in the macro-element. Since recomputation of integrals of a macro-element leads to the need of updating the flow kinetics system matrix, this would lead to the increase of overall computation time and could to some extent cancel out positive effects of freezing the velocity in the fundamental solution. To overcome this problem a simple but effective strategy is to implement the ratio of number of macro-elements Nnew, where the criterion (40) was reached, to the overall number of macro-element NME, N new P int N ME
ð41Þ
as the criterion for recomputation of integrals and hence recomputation of system matrix entries. The effect of implementing criteria (40) and (41) can be seen in Fig. 3. Although the full update of integrals of flow kinetics in each inner loop iteration does result in the fastest convergence in the early stage of iterating, the overall convergence for different int is quite similar and leads to the conclusion, that freezing the integrals for several inner iterations does not significantly influence the overall stability of the BDIM iteration procedure. The value of int = 0.05 was therefore chosen and used in all computations in the test examples section. 6.3. Recomputation of system matrix of flow kinetics In order to construct a computationally favourable BDIM algorithm we have to deal also with the process of reformation of system matrix of flow kinetics. The need for the reformation of this matrix comes from the change in integrals but also from including the convection contributions in the entries of the system matrix of flow kinetics, corresponding to unknown vorticity values, i.e.
4206
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
0 -0.5 -1
log(C)
-1.5 epsint=0% epsint=5% epsint=10% epsint=25%
-2 -2.5 -3 -3.5 -4
10
20
30
40
50 60 iteration
70
80
90
100
Fig. 3. Convergence of BDIM scheme for different thresholds (epsint = int) for recomputation of integrals of kinetics. Test case of driven cavity at Re = 100 with non-equidistant 10 · 10 mesh and Dt = 1016.
1 1 ½ACME þXME ¼ ½H þ ½C j ½vj ½Dj ½^vj m m with the corresponding right-hand side ox 1 ½DfxgF 1 : fbg ¼ ½G þ on mDt
ð42Þ
ð43Þ
In this way the system matrix contains the majority of information and the right-hand side of system of equations consists mainly of the known boundary conditions and vorticity contribution from the previous time step. Unfortunately, this means that the flow kinetics system matrix has to be updated in each inner iteration due to the change of vj and hence ^vj . By freezing velocity values of terms in the system matrix, say vj,f, for several inner loop iterations, the system matrix remains constant, 1 1 ½ACME þXME ¼ ½H þ ½C j ½vj;f ½Dj ½vj;f vj m m and only a change in the vector of the right-hand side has to be made: ox 1 1 ½DfxgF 1 : fbg ¼ ½G þ ð½Dj ½C j Þ½vj vj;f fxg þ on m mDt
ð44Þ
ð45Þ
In the case of no freezing, the vj = vj,f and the second term on the right-hand side is zero. In the case of freezing velocities in the system matrix the update of integration is synchronized with the freezing process. An example for this can be seen from comparison of percentage of recomputed integrals of flow kinetics for different numbers of frozen iterations (nfroz), Fig. 4. Selecting larger values of frozen iterations do increase the percentage of newly integrated macro-elements, but the overall number of recomputed integrals is still significantly lower than for the case of update of velocity field for system matrix in each inner iteration (nfroz = 1).
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4207
% of integrated macro elements
100
80
60 nfroz=1 nfroz=10 nfroz=25
40
20
20
40
60
80
100 120 iteration
140
160
180
200
Fig. 4. Comparison of percentage of macro-elements where new integrals were computed based on relative change of 5% in velocity field. Test case as in Fig. 3.
The comparison of convergence of the C error for different numbers of frozen iterations is shown in Fig. 5. Differences in convergence of BDIM in the early stage of inner iteration loop result from the fact that the velocity field, used for the computation of initial integrals of kinetics, represents the solution to the potential flow in the computational domain, as x = 0 is the initial guess and hence the first solution of flow kinematics solves the equation o2 v i ¼ 0: oxj oxj
ð46Þ
0
log(C)
-1
-2
nfroz=1 nfroz=10 nfroz=25
-3
-4
-5 20
40
60
80 100 iteration
120
140
160
180
Fig. 5. Comparison of inner iteration loop convergence for different freezing steps. Test case as in Fig. 3.
4208
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
This is the reason for slower decrease of error C when increasing the number of frozen iterations, as shown in Fig. 5. After the update of velocity field in the integrals of flow kinetics the convergence significantly increases, although in the later stage of inner iteration loop the convergence tends to be similar regardless of the number of frozen iterations. 6.4. Use of diffusion–convection fundamental solution for large cell Re number values The argument of the exponential function in the diffusion–convection fundamental solution (20), pffiffiffiffiffiffiffiffiffi vME AME xe ; ð47Þ 2m where AME is the area of a macro-element, can present numerical problems due to limitations of computer processor in handling very large numbers (overflow). For a given fluid the xe depends on macro-element size and velocity magnitude. In order to avoid computational problems either mesh size or velocity magnitude should be decreased. With a given computational mesh the second possibility is applicable, therefore macro-element based maximum velocity magnitude is defined as 2mxe;max vME;max ¼ pffiffiffiffiffiffiffiffiffi ; AME
ð48Þ
where the xe,max is the largest argument that does not result in an overflow during computation. If the actual magnitude of the average velocity vector in a macro-element exceeds the value (48), local scaling is applied, i.e. vME;max ¼ /MEvME ;
ð49Þ
where /ME is a scaling factor. The scaling is performed simultaneously with recomputation of integrals of flow kinetics, 6.2, and the difference ð1 /ME ÞvME of average part of velocity is treated similar to frozen iteration approach, described in Section 6.2. An interesting observation can be made here. The argument (47) can be equated to half of a macroelement Re number, pffiffiffiffiffiffiffiffiffi vME AME 1 ¼ ReME ð50Þ 2 2m pffiffiffiffiffiffiffiffiffi with characteristic length L ¼ AME and characteristic velocity v ¼ vME applied in definition of ReME. For the AMD based 32-bit personal computer, on which all the computations were made, the maximum macroelement Re number was ReME = 2 · 308 = 616. 6.5. Iterative solvers for solution of resulting systems of equations To lower computational cost of the method, iterative Krylov subspace type solvers are used for the solution of resulting systems of equations: • diagonally preconditioned CGS or GMRES solver for flow kinematics system of equations [4], • level based ILU preconditioned CGS or GMRES solver for flow kinetics system of equations [7]. When compared with block-oriented ILU preconditioning, used in [4], the macro-element discretization allows more efficient use of level based ILU preconditioned, developed in [7]. Freezing the velocities in system matrix of kinetics has a consequence of freezing the preconditioned, so the same preconditioned is used
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4209
for computation of all frozen iterations. Thus, freezing of velocity field in integrals and system matrix of kinetics greatly decreases computational cost of the proposed algorithm. The developed BDIM algorithm has also lower memory demands as the original method [4], which is due to the macro-element discretization of flow kinetics instead of using large subdomains (i.e. segments) as done in [4]. There are further parts of BDIM that need improvements regarding memory demands, especially flow kinematics with its fully occupied system matrix. Wavelet transform techniques present here an effective possibility, a work that is already under development for the BDIM.
7. Test problems The presence of non-smooth boundaries can significantly affect the fluid flow conditions as compared with the flows in confined rectangular cavities. Abrupt changes in velocity profiles, caused by geometric irregularities, are connected with large values of vorticity, therefore an accurate determination of boundary vorticity values is a key step towards accurate solution of the governing set of equations. Except for the first test case, the driven cavity problem, the selected test cases have abrupt changes in geometry. In order to accurately resolve high gradients of velocities and vorticities, highly stretched computational meshes have to be used [19] and the following examples show the application of such meshes in BDIM. 7.1. Flow in a square driven cavity Fluid flow in a square cavity is one of the most difficult benchmark problems in flows of incompressible fluids, which is currently used for testing different numerical schemes for solution of Navier–Stokes equations [14–17,20]. The imposed velocity field at the top moving wall of the cavity drives a large recirculation region inside the cavity. With increasing moving wall velocity and hence the Re number value, additional smaller recirculation zones appear in the corners of the cavity. The complexity of the flow field in connection with a simple geometry makes this problem ideal for testing CFD codes. On the other side, the driven cavity problem can be found in a variety of important engineering applications, like flows in mixing vessels equipped with disk turbines, flows in polymer extruders and environmental flows. The test case was computed on two mesh densities, Fig. 6. The coarse mesh consisted of 30 · 30 macroelements with ratio 4 between the shortest and the longest element side, and 4 large segments for computation of flow kinematics. The fine mesh consisted of 40 · 40 macro-elements with ratio 20, and also 4 segments. The comparison of results between the two meshes showed that the coarse mesh gave accurate results up to Re = 5000, but was unable to give a converged solution for Re = 10,000. For Re = 10,000, the fine mesh was used. A check of mesh convergence with a 60 · 60 ratio 20 mesh with 4 subdomains
Fig. 6. Computational meshes for driven cavity: 30 · 30 with ratio 4 (left) and 40 · 40 macro-elements with ratio 20 (right).
4210
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
showed no improvement in the quality of results, therefore the 40 · 40 mesh was used for all computations of Re = 5000 and higher. All computations were performed by time stepping procedure, where the typical time steps were Dt = 1 for Re 1000 and Re = 5000 and Dt = 0.1 for Re = 7500 and Re = 10,000 (Figs. 7–13). For computation of
Fig. 7. Streamlines and vorticity isolines for Re = 1000.
Fig. 8. Streamlines and vorticity isolines for Re = 5000.
Fig. 9. Streamlines and vorticity isolines for Re = 7500.
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4211
Fig. 10. Streamlines and vorticity isolines for Re = 10,000.
higher value Re number flows starting values of lower Re number value computations were used. Underrelaxation parameter was 0.01 for all the computed cases. In Table 1 a comparison of primary vortex strength values with results of other authors is presented. 7.2. Flow in L shaped driven cavity In the flow in L shaped driven cavity case, the motion of fluid in the cavity is caused by imposed velocity field at the top of the cavity, which, due to the action of viscous forces, drives the recirculating flows inside the lower right and upper left cavity. In contrast to the standard driven cavity flow, the presence of a step inside the cavity causes the flow field to form three large recirculation zones, two in the upper part of the cavity and one in the lower, smaller part of the cavity (Fig. 14). The coarse computational mesh consisted of 300 macro-elements and the dense one consisted of 1200 (three times 20 · 20) macro-elements. In both cases 3 segments were used for segmentation of flow kinematics (Fig. 15).
1
Ghia Re=1000 BDIM Re=1000 Ghia Re=5000 BDIM Re=5000 Ghia Re=7500 BDIM Re=7500 Ghia Re=10,000 BDIM Re=10,000
0.8
y
0.6
0.4
0.2
0 -0.6
-0.4
-0.2
0
0.2 Vx(0.5,y)
0.4
0.6
Fig. 11. Comparison of vx velocity profiles at x = 0.5.
0.8
1
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4212
0.6
0.4
Vy(0.5,y)
0.2
0 Ghia Re=1000 BDIM Re=1000 Ghia Re=5000 BDIM Re=5000 Ghia Re=7500 BDIM Re=7500 Ghia Re=10,000 BDIM Re=10,000
-0.2
-0.4
-0.6 0
0.2
0.4
0.6
0.8
1
x
Fig. 12. Comparison of vy velocity profiles at y = 0.5.
200
0
-200
-400 Ghia Re=1000 BDIM Re=1000 Ghia Re=5000 BDIM Re=5000 Ghia Re=7500 BDIM Re=7500 Ghia Re=10,000 BDIM Re=10,000
y
-600
-800
-1000
-1200
-1400
-1600 0
0.2
0.4
0.6
0.8
1
w(x,1.0)
Fig. 13. Comparison of vorticity profile at the moving boundary.
The computations were performed for Re number values in the range 1000–7500 (Figs. 16–19). Here, the Reynolds number value was defined as (h stands for the total height of the cavity wall), Re ¼
vtop h : m
ð51Þ
As seen from the computation results the irregular geometry of the cavity causes formation of three large vortices (Table 1). The main vortex is at the right top part of the cavity and generally exhibits the characteristics of the main cavity formed in the driven cavity test case except for the lack of smaller corner recirculation regions, which is due to extension of the cavity to the left and to the bottom from the sharp corner.
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4213
Table 1 Comparison of primary vortex strength values Method
Re number value 1000
5000
7500
10,000
Wmin, present, BDIM Wmin [14] FDM Wmin [17] FEM Wmin [9] Poly-BEM Wmin [15] FDM Wmin [3] BEM–BDIM
0.1187 0.1179 0.1189 0.1193 0.1157 0.113
0.1199 0.1189 – 0.1221 0.1123 –
0.1245 0.1199 0.1224 – 0.1052 –
0.1307 0.1197 0.1224 – – 0.109
FDM—finite difference method, FEM—finite element method, BEM—boundary element method, BDIM—boundary domain integral method, Poly-BEM—poly-region boundary element method.
Fig. 14. Presentation of L shaped driven cavity problem.
Fig. 15. Computational mesh for L shaped driven cavity: 300 macro-elements with ratio 4 (left) and 1200 macro-elements with ratio 10 (right).
As a consequence of this large vortex two additional vortices appear in the computational domain, one left and the other at the bottom of the main vortex. They are rotating in the direction opposite to the main vortex.
4214
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
Fig. 16. Streamlines and vorticity isolines for Re = 1000.
Fig. 17. Streamlines and vorticity isolines for Re = 3200.
Fig. 18. Streamlines and vorticity isolines for Re = 5000.
Interesting flow behaviour can be observed on the top of the left vortex, where a thin stratified layer between the moving top wall and the left vortex is formed. It is becoming thinner as the value of Re number is increasing and more and more it exhibits characteristics of a semi closed recirculation region, identified also in the vorticity value increase in the middle of the moving top wall, Fig. 22.
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4215
Fig. 19. Streamlines and vorticity isolines for Re = 7500.
1
0.8
y
0.6
Re=1000 Re=3200 Re=5000 Re=7500
0.4
0.2
0 -0.6
-0.4
-0.2
0
0.2 Vx(0.75,y)
0.4
0.6
0.8
1
Fig. 20. Comparison of vx velocity profiles at x = 0.75.
Each of the additional large vortices has its own recirculation regions at the corners, the left vortex one and the bottom vortex two recirculation regions. They are becoming larger as the Re number value is increased. In the case of the bottom vortex results for Re = 7500 show that a new recirculation region is beginning to form on the right wall between the main and the bottom vortex. The comparison of velocity profiles along the lines x = 0.75 and y = 0.75, presented in Figs. 20 and 21 shows that the main characteristics of the driven cavity flow problem are replicated also in all of the three vortices: almost linear velocity profiles and constant vorticity values inside the vortex regions together with large velocity gradients and hence large vorticity values at vortex boundaries. 7.3. Flow in a Z shaped channel The last test example deals with the flow in a channel that has two bends and the total geometry has a Z-like form. The test case is interesting because of the sudden expansion over the first bend and sudden
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4216
0.6
0.4
Vy(x,0.75)
0.2
0
-0.2 Re=1000 Re=3200 Re=5000 Re=7500
-0.4
-0.6
-0.8 0
0.2
0.4
0.6
0.8
1
x
Fig. 21. Comparison of vy velocity profiles at y = 0.75.
0
-100
vorticity
-200
-300 Re=1000 Re=3200 Re=5000 Re=7500
-400
-500
-600 0
0.2
0.4
0.6
0.8
1
x
Fig. 22. Comparison of vorticity profiles at the moving boundary.
contraction over the second bend which influence the main flow path. Both abrupt changes in the geometry have ratio 2. The test case is interesting also because of mixed type boundary conditions, as the flow conditions at the outlet boundary are not known in advance and are a result of the numerical simulation. The geometry of the channel is presented in Fig. 23. The dimensions correspond to those of [18]. The Re number is defined as Re ¼
vh ; m
ð52Þ
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4217
Fig. 23. Z shaped channel geometry, boundary conditions.
where v is the average velocity at the inlet of the channel and h is the height of the channel at its inlet. In our case v ¼ 1 and h = 1 was chosen. Boundary conditions were prescribed as follows: at the inlet a fully developed laminar flow profile was prescribed, at the outlet the developed flow was assumed by vt = 0 and ox ¼ 0, and the rest of boundaries on was prescribed as no-slip solid wall, i.e. vx = 0, vy = 0. The result of applying the described outlet boundary conditions leads to mixed type boundary conditions in flow kinematics as well as in flow kinetics. In flow kinematics, vorticity is unknown all over the boundary C except for the outlet, where velocities normal to the outlet boundary are unknowns. In vorticity transport equation, the unknowns at the boundary are vorticity fluxes except for the outlet, where vorticities are unknowns. This combination of boundary conditions is another test for the numerical stability of the proposed computational method. The computational domain consisted of 896 macro-elements and seven segments (I–VII, Fig. 23). Computations were performed for Re = 100, 200 and 400. The time step and the under-relaxation values in computations were 0.25 and 0.01, respectively. Computation of higher value Re number flows was done in the same way as in all other test cases, i.e. by using results of the lower Re number value flow as initial conditions for the next higher Re number value. Computational results, presented in Figs. 24 and 25 show that three main recirculation regions are formed regardless of the value of Reynolds number values under consideration. The first one is in the upper right corner, the second is in the lower left corner and the third is formed behind the second bend on the upper wall of the contracted region. For the case of the highest computed Re number value, Re = 400, additional recirculation regions occur. The first additional region is formed along the left wall of the expanded region and is a consequence of the stronger recirculation in the lower left corner. The length of this region is clearly visible from Fig. 26, where vorticity values on the solid wall at x = 2.0 are plotted for various Re number values. The second rather small additional recirculation region appears behind the large recirculation region in the lower left corner. From Fig. 26 one can see that this region is starting to form already at Re = 200, although it becomes pronounced at Re = 400. The length of this region is 0.3 units for Re = 400 and 0.19 units for Re = 200. Also, high vorticity gradients at the sharp corner of the expansion can be observed, although because of the quadratic interpolation functions used these values show oscillations in the vicinity of the corner. The third additional recirculation region is a consequence of redirection of the mean flow after it hits the lower solid wall at the beginning of the contracted region. This region forms on the lower solid wall of the contracted region and is comparable in length with the recirculation region behind the contraction corner. The length of this region for Re = 400 is 1.73 units.
4218
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
Fig. 24. Z shaped channel flow at Re = 200: streamlines (top) and vorticity isolines (bottom).
Fig. 25. Z shaped channel flow at Re = 400: streamlines (top) and vorticity isolines (bottom).
Results, obtained for Re = 200 are in good agreement with results in [18], where QUICK finite differencing numerical scheme was used.
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4219
4.5
4
3.5
3
Re=100 Re=200 Re=400
y
2.5
2
1.5
1
0.5
0 0
2
4
6
8
10
vorticity on solid wall at (2.0,y)
Fig. 26. Comparison of vorticity x values along the solid wall at x = 2.0.
8. Conclusions A new computational algorithm for the numerical solution of incompressible fluid flow on basis of BEMs has been presented. The algorithm is based on the BDIM and combines two different approximation techniques, approximation of flow kinematics by means of boundary elements and internal cells and approximation of flow kinetics by means of macro-elements. The first ensures divergence free numerical solution and the latter, in connection with the diffusion–convection fundamental solution, delivers accurate solutions to diffusion–convection transport phenomena at higher values of Re number. In order to decrease computational cost of the method, segmentation technique is applied in the flow kinematics computation. Computed test examples of flow in complex geometries show that the developed algorithm is capable of accurate solutions of high Re value flows. Since the BDIM presents an extension of the BEM it is straightforward to couple the method with the established BEM methods and to make use of the advantages BEM has to offer.
References [1] J.C. Wu, J.F. Thompson, Numerical solution of time dependent incompressible Navier–Stokes equations using an integrodifferential formulation, Comput. Fluids 1 (1973) 197–215. [2] L. Skerget, A. Alujevic, C.A. Brebbia, G. Kuhn, Natural and forced convection simulation using the velocity–vorticity approach, in: C.A. Brebbia (Ed.), Topics in Boundary Element Research, vol. 5, Springer-Verlag, Berlin, 1989 (Chapter 4). [3] Z. Rek, L. Skerget, Boundary element method for steady 2D high-Reynolds-number flow, Int. J. Numer. Methods Fluids 19 (4) (1994) 343–362. [4] M. Hribersˇek, L. Skerget, Iterative methods in solving Navier–Stokes equations by the boundary element method, Int. J. Numer. Methods Engrg. 39 (1996) 115–139. [5] L. Sˇkerget, M. Hribersˇek, G. Kuhn, Computational fluid dynamics by boundary domain integral method, Int. J. Numer. Methods Engrg. 46 (1999) 1291–1311. [6] M. Ramsak, L. Sˇkerget, Mixed boundary elements for laminar flows, Int. J. Numer. Methods Fluids 31 (1999) 861–877. [7] M. Hribersˇek, L. Sˇkerget, Fast boundary-domain integral algorithm for computation of incompressible fluid flow problems, Int. J. Numer. Methods Fluids 31 (1999) 891–907. [8] L.C. Wrobel, Applications in thermo-fluids and acoustics. The Boundary Element Method, vol. 1, Wiley, New York, 2002.
4220
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
[9] M.M. Grigoriev, G.F. Dargush, A poly-region boundary element method for incompressible viscous fluid flows, Int. J. Numer. Methods Engrg. 46 (1999) 1127–1158. [10] W.F. Florez, H. Power, Comparison between continuous and discontinuous boundary elements in the multidomain dual reciprocity method for the solution of the two-dimensional Navier–Stokes equations, Engrg. Anal. Bound. Elem. 25 (2001) 57–69. ˇ unicˇ, Natural convection flows in complex cavities by BEM, Int. J. Numer. Methods Heat Fluid [11] L. Sˇkerget, M. Hribersˇek, Z. Z Flow 13 (5/6) (2003) 720–736. [12] M.S. Ingber, A vorticity method for the solution of natural convection flows in enclosures, Int. J. Numer. Methods Heat Fluid Flow 13 (5/6) (2003) 655–671. [13] W. Hackbusch, Iterative Lo¨sung Grosser Schwachbesetzter Gleichungssysteme, B.G. Teubner, Stuttgart, 1991. [14] U. Ghia, K.N. Ghia, C.T. Shin, High-Re solutions for incompressible flow using the Navier–Stokes equations and a multigrid method, J. Comput. Phys. 48 (1982) 387–411. [15] O. Goyon, High-Reynolds number solutions of Navier–Stokes equations using incremental unknowns, Comput. Methods Appl. Mech. Engrg. 130 (1996) 319–335. [16] C.H. Liu, D.Y.C. Leung, Development of a finite element solution for the unsteady Navier–Stokes equations using projection method and fractional H-scheme, Comput. Methods Appl. Mech. Engrg. 190 (2001) 4301–4317. [17] E. Barragy, G.F. Carey, Stream function–vorticity driven cavity solution using p finite elements, Comput. Fluids 26 (5) (1997) 453–468. [18] C.Y. Perng, R.L. Street, A coupled multigrid-domain-splitting technique for simulating incompressible flows in geometrically complex domains, Int. J. Numer. Methods Fluids 13 (1991) 269–286. [19] P. Sockol, Multigrid solution of the Navier–Stokes equations on highly stretched grids, Int. J. Numer. Methods Engrg. 17 (1993) 543–566. [20] S. Turek, Tools for simulating non-stationary incompressible flow via discretely divergence free finite element models, Int. J. Numer. Methods Fluids 18 (1994) 71–105.
Boundary domain integral method for high Reynolds viscous fluid flows in complex planar geometries M. Hribersˇek *, L. Sˇkerget Faculty of Mechanical Engineering, Laboratory of Computational Fluid Dynamics, Institute of Power, Process and Environmental Engineering, University of Maribor, Smetanova 17, SI-2000 Maribor, Slovenia Received 9 February 2004; received in revised form 5 July 2004; accepted 2 November 2004
Abstract The article presents new developments in boundary domain integral method (BDIM) for computation of viscous fluid flows, governed by the Navier–Stokes equations. The BDIM algorithm uses velocity–vorticity formulation and is based on Poisson velocity equation for flow kinematics. This results in accurate determination of boundary vorticity values, a crucial step in constructing an accurate numerical algorithm for computation of flows in complex geometries, i.e. geometries with sharp corners. The domain velocity computations are done by the segmentation technique using large segments. After solving the kinematics equation the vorticity transport equation is solved using macro-element approach. This enables the use of macro-element based diffusion–convection fundamental solution, a key factor in assuring accuracy of computations for high Reynolds value laminar flows. The versatility and accuracy of the proposed numerical algorithm is shown for several test problems, including the standard driven cavity together with the driven cavity flow in an L shaped cavity and flow in a Z shaped channel. The values of Reynolds number reach 10,000 for driven cavity and 7500 for L shaped driven cavity, whereas the Z shaped channel flow is computed up to Re = 400. The comparison of computational results shows that the developed algorithm is capable of accurate resolution of flow fields in complex geometries. 2004 Elsevier B.V. All rights reserved. Keywords: Boundary element algorithms; Incompressible viscous fluid; Navier–Stokes equations; Velocity–vorticity formulation; Segmentation technique; Driven cavity flow
*
Corresponding author. Tel.: +386 62 220 7772; fax: +386 62 220 7990. E-mail address: [email protected] (M. Hribersˇek).
0045-7825/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2004.11.002
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4197
1. Introduction The origins of the boundary element method (BEM) are strongly connected with computation of linear or weakly non-linear field problems, which can be computed by means of boundary-only discretization of computational domain. In case of flows of viscous fluids at higher values of Re number, the flow phenomena is strongly non-linear and boundary-only discretization had to be paired with domain discretization as well. In the context of BEM related methods for laminar viscous fluid flows several successful attempts have already been made, e.g. [2,5,9,10,12] among others. An excellent survey of these approaches can be found in [8]. These numerical approaches were based on different forms of Navier–Stokes equations, representing the frame for the solution of viscous flow problems. Different techniques of capturing non-linear domain effects were developed, including internal cells, macro-elements and dual and multiple reciprocity methods. Different approaches were successful also in cases of higher values of Re number flows, although the highest values, like for the case of driven cavity, computed with other established approximation methods (Re = 10,000), were not easily obtained [3] (Re = 10,000) or were restricted to lower or moderately high Re number values [10] (up to Re = 1000) and [9] (up to Re = 5000). In boundary domain integral method (BDIM), velocity–vorticity formulation was used for the solution of Navier–Stokes equations [2]. Its main advantage is an implicit computation of boundary vorticity values, which presents a solid basis for development of a successful computational technique. The use of different types of fundamental solutions, especially macro-element based diffusion–convection fundamental solution [5] led to development of flexible and accurate computational methods [3,7] able to accurately resolve also higher value Re number flows. In obtaining a divergence-free numerical solution for all types of geometry a special attention has to be paid to a proper numerical treatment of governing equations [20]. In case of BDIM, vector potential formulation [4] and vector–velocity formulation of flow kinematics [5] were already used. In the first case in flow kinetics the parabolic-diffusion fundamental solution was used [2] and in the second the diffusion–convection fundamental solution was used. In [11] it was shown that both approaches result in the equal final integral statement. In order to lower computational cost of the method subdomain technique was used with vector–potential formulation [4] although it was restricted to non-star arrangement of subdomains and to segmentation of flow kinematics. On the other hand, the velocity vector formulation of flow kinematics [5] allowed the use of macro-element based subdomain technique in both flow kinematics and flow kinetics, but it lacked the conservation properties when applied to complex geometries. In order to overcome these drawbacks the paper presents a further development of BDIM technique, combining the velocity Poisson equation for flow kinematics and segmentation of flow kinematics together with macro-element based computation of flow kinetics with the use of diffusion–convection fundamental solution. The result is a stable numerical technique, able of accurate computation of high value Re flows in complex geometries. The new computational algorithm can be viewed as a development of the BDIM method [4], where the flow kinetics is computed in a macro-element based domain discretization with the use of diffusion–convection fundamental solution. Due to this the BDIM computational algorithm has to be modified in order to combine segmentation technique of flow kinematics with macro-element approach used in flow kinetics. Additionally, the frozen integrals and matrix technique has to be developed in order to lower computational demands of flow kinetics. The paper is organized as follows: Section 2 deals with the set of governing equations in its primitive and velocity–vorticity form, followed by Section 3 describing transformation of partial differential forms of governing equations to its integral forms, representing the starting point for discretization by boundary elements and domain cells, described in Section 4. In Section 5 the use of subdomain and macro-element technique in conjunction with segmentation technique in flow kinematics is explained. Section 6 deals with computational aspects of the derived algorithm, especially with the frozen integrals and frozen matrix
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4198
approach, and Section 7 concentrates on computational results for selected test cases: driven cavity flow, flow in L shaped driven cavity and flow in a Z shaped channel. The paper closes with conclusions based on derived numerical algorithm and obtained numerical results.
2. Governing equations When developing a numerical method for solution of viscous incompressible fluid flows the Navier– Stokes set of equations serves as the starting point for derivation of the final algebraic equations, suitable for numerical solution. For the numerical solution of the Navier–Stokes equation set we choose velocity– vorticity formulation, which has several advantages when used in the context of BDIM. When using this formulation the dynamics of a viscous incompressible fluid is partitioned into its kinematic and kinetic aspect through the use of derived vector vorticity field function xi(rj, t). Vorticity is obtained as a curl of the velocity field vi(rj, t), namely xi ¼ eijk
ovk oxj
ð1Þ
and thus xi is solenoidal, oxi ¼ 0: oxi
ð2Þ
By applying the curl operator to the vorticity definition and by using the continuity constraint for the incompressible flows the following vector elliptic Poisson equation for velocity vector is obtained: o2 v i oxk þ eijk ¼ 0: oxj oxj oxj
ð3Þ
Eq. (3) represents the kinematics of an incompressible fluid motion, expressing the compatibility and restriction conditions between velocity and vorticity field functions. The kinetic aspect is governed by the parabolic diffusion–convection vorticity transport equation obtained by applying the curl operator to the momentum equation, oxi oxi oxj vi o2 x i þ vj ¼ þm : ot oxj oxj oxj oxj
ð4Þ
The present investigation deals with constant material properties, kinematic viscosity m. One of the advantages of using velocity–vorticity formulation lies in removing pressure gradient terms from the solution process, resulting in a higher numerical stability of the computational scheme. The pressure does not appear in the solution procedure and has no influence on the velocity field, a fact that is of course valid only for incompressible fluid approximation. The pressure can be computed after velocity field is being solved. For the solution of the Poisson pressure equation the Neumann conditions are specified over the boundary of the domain, but at least in one point the pressure value has to be specified. This could well be the point on the inlet or outlet, thus the computed pressure field will correspond to the specified boundary condition and computed velocity field. In the case of compressible fluid the pressure is a thermodynamic variable and has direct influence on the density of the fluid, which further influences velocity and vorticity fields. The velocity–vorticity BEM method is currently being extended to the solution of this class of problems. For the case of the 2D flow, which is considered in the paper, the vorticity vector has just one component, perpendicular to the plane of the flow, simplifying Eqs. (3) and (4) to
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4199
o2 v i ox þ eij ¼ 0; oxj oxj oxj
ð5Þ
ox ox o2 x þ vj ¼m : ot oxj oxj oxj
ð6Þ
3. Integral representations The unique advantage of the BDIM originates from the application of GreenÕs fundamental solutions as particular weighting functions. Since they only consider the linear transport phenomenon, an appropriate selection of a linear differential operator L½ is of main importance in establishing stable and accurate singular integral representations corresponding to the original differential conservation equations. All different conservation models can be written in the following general form: L½u þ b ¼ 0;
ð7Þ
where the operator L½ can be either elliptic or parabolic, u(rj, t) is an arbitrary field function, and the nonhomogeneous term b(rj, t) is applied for non-linear transport effects or pseudo body forces. In the following, integral representations will be derived using different fundamental solutions based on particular type of governing equation under consideration, i.e. elliptic Poisson partial differential equation (PDE) for flow kinematics and elliptic diffusion–convection PDE for vorticity transport. 3.1. Integral representations for flow kinematics When transforming the elliptic Poisson PDE (3) into algebraic form the first step is the choice of the linear differential operator L½, which in the case of (3) reads as L½ ¼
o2 ðÞ ; oxj oxj
ð8Þ
allowing us to rewrite Eq. (3) as L ½ v i þ bi ¼
o2 v i þ bi ¼ 0; oxj oxj
ð9Þ
where bi ¼ eijk
oxk : oxj
ð10Þ
There are two possibilities in deriving the final integral statement of flow kinematics, suitable for the discretization and consequent use in the computational BDIM algorithm. The first one was presented in [5] and allowed the use of macro-element based discretization. Despite the low computational demands of this approach, resulting in sparse flow kinematics matrices, the formulation was limited to simple geometries with prescribed no-slip boundary conditions. In [1] an excellent explanation for this behaviour can be found noting that the solution of flow kinematics equation in [5] does not necessarily satisfy both continuity constraint and conservation equation (2), although vice versa is always true. In order to exclude non-conservative solutions additional enforcement of continuity constraint in equation has to be undertaken, as was shown recently in [11]. This formulation presents the basis for the flow kinematics integral equations, used with the presented BDIM computational algorithm. The derivation of
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4200
the integral form starts with the choice of uw as the Laplace fundamental solution, application of the GreenÕs theorems for scalar functions and additional implementation of the continuity constraint and vorticity definition, the final integral form of flow kinematics is obtained [11] which for the planar case reads as Z Z Z ouH ouH ouH dC ¼ dC eij x vj dX ð11Þ cðnÞvi ðnÞ þ vi on os oxj C C X and for the general 3D case Z Z Z ~ H ~ ~ H ~ ~ H Þ dX: cðnÞ~ vðnÞ þ ðru nÞ~ v dC ¼ ðru nÞ ~ v dC þ ð~ x ru C
C
ð12Þ
X
Eqs. (11) and (12) represent the starting point for discretization by means of boundary elements and domain cells. It can be used for evaluating unknown velocities in the interior of computational domain, whereas for determination of values on the boundary of computational domain (mainly vorticities) the discretization leads to the singular system matrix [2]. Therefore for the boundary, when the unknowns are the boundary vorticity values or the tangential velocity component values, one has to use the tangential component of the vector, Eq. (12), namely Z Z Z H H ~ ~ ~ H Þ dX cðnÞ~ nðnÞ ~ vðnÞ þ ~ nðnÞ ðru ~ nÞ~ v dC ¼ ~ nðnÞ ðru ~ nÞ ~ v dC þ ~ nðnÞ ð~ x ru C
C
X
ð13Þ in order to obtain the appropriate non-singular implicit system of equations. When the normal velocity components to the boundary are unknown, the normal form has to be employed: Z Z Z ~ H ~ ~ H ~ ~ H Þ dX: cðnÞ~ nðnÞ ~ vðnÞ þ ~ nðnÞ ðru nÞ~ v dC ¼ ~ nðnÞ ðru nÞ ~ v dC þ ~ nðnÞ ð~ x ru ð14Þ C
C
X
3.2. Integral representations for flow kinetics We discretize
du dt
in Eq. (4) by the backward Euler scheme approximation, given by
ou uF uF 1 ; ot Dt
ð15Þ
where Dt is the time step and F denotes the time step index, the linear differential operator can be written in the following diffusion-convective non-homogeneous PDE form, namely L½ ¼ m
o2 ðÞ oðÞ ðÞ vj ; oxj oxj oxj Dt
ð16Þ
where the velocity field vj was decomposed into an average constant vector vj and a perturbed one ^vj , i.e. vj ¼ vj þ ^vj . The following can now be stated: L½xi þ bi ¼ m
o2 xi oxi xi vj þ bi oxj oxj oxj Dt
ð17Þ
with bi ¼ ^vj
oxi oxj vi xi;F 1 : þ þ oxj oxj Dt
The second term on the right-hand side is zero in case of 2D flows considered in the present paper.
ð18Þ
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4201
As fundamental solution u* the diffusion–convection fundamental solution is implemented, being the solution of equation m
o2 u ou 1 u þ dðn; sÞ ¼ 0; þ vj oxj oxj oxj Dt
ð19Þ
given for the planar case as 1 ðvj rj Þ K 0 ðlrÞ exp ; 2p 2m 2 v 1 l2 ¼ ; þ 2m m Dt u ¼
ð20Þ
ð21Þ
where K0 is modified Bessel functions of the second kind, v is the magnitude of the average part of the velocity vector, s is the field point and n the source point. By applying the GreenÕs theorem [5] the following integral statement is obtained: Z Z Z Z Z ou ox 1 1 ou 1 u dC dC ¼ cðnÞxðnÞ þ x xvn u dC þ x^vj dX þ xF 1 u dX ð22Þ m C m X mDt X on oxj C C on with vn = vjnj. The velocity decomposition is necessary since the fundamental solution of Eq. (16) is known only for the linear diffusion-convective PDE with constant coefficients, i.e. known average velocity, the result of the evaluation of flow kinematics equations in the current BDIM iteration step.
4. Algebraic representations of integral equations The set of integral forms of governing equations, Eqs. (11) and (22), is transformed to its algebraic form first by dividing boundary C into boundary elements Ce and domain X into domain cells Xc. Second, geometry, the field functions vi, x and their derivatives are approximated by the use of quadratic interpolation functions [5]. Due to the use of macro-element approach in discretization of flow kinetics equations different positions of collocation nodes on boundary elements for flow kinematics and flow kinetics have to be selected, leading to two types of implemented boundary elements, Fig. 1: • Continuous quadratic elements: collocation nodes coincide with geometry nodes, this discretization is used for boundary elements in flow kinematics. • Discontinuous quadratic elements: collocation nodes are shifted from corners into the interior of the boundary element, thus avoiding the problem of obtaining an overdetermined system of equations in case of macro-element discretization approach [6] where these elements are used.
Fig. 1. Discrete model: continuous and discontinuous quadratic boundary elements (left), continuous quadratic domain cell (right).
4202
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
The detailed description of the implemented macro-element approach is given in Section 5. The use of continuous boundary elements in flow kinematics discretization ensures continuity of boundary vorticity values at sharp corners, what is of great advantage over the previous BDIM algorithm [5] where discontinuous boundary elements were used for both, flow kinematics and flow kinetics computations. After collocating over all nodes one obtains the final system of algebraic equations. In case of flow kinematics the matrix form of the final system of equations reads as ½H e fvx g ¼ ½H t fvy g ½Dy fxg;
ð23Þ
½H e fvy g ¼ ½H t fvx g þ ½Dx fxg:
ð24Þ
Forms (23) and (24) are discretized forms of Eq. (11) and are used for computation of velocity values for domain nodes. For computation of unknowns on the boundaries of computational domain equations (13) and (14) have to be used. After assembling equations written for all source points, the final matrix statements for tangential and normal forms read as ð½H 1 þ ½H t1 Þfvt g ¼ ð½H 2 ½H t2 Þfvn g þ ½Dn fxg;
ð25Þ
ð½H 1 þ ½H t1 Þfvn g ¼ ð½H t2 ½H 2 Þfvt g þ ½Dt fxg:
ð26Þ
In the case of vorticity kinetics the final matrix form of equations is obtained after discretization of Eq. (22) and reads as ox 1 1 ½DfxgF 1 : ½H fxg ¼ ½G ð27Þ þ ð½Dj ½C j Þ½^vj fxg þ on m mDt The elements of matrix [Cj] are computed using continuous boundary elements, which is necessary since both [Dj] and [Cj] are the result of transformation of the first term on the right-hand side bi of flow kinetics, Eq. (18), into the boundary and domain integral parts. In Eqs. (23)–(27), all values of field functions correspond to the actual time step (denoted also as F) and (F 1) denotes the previous time step values.
5. Subdomain and macro-element technique One of the main drawbacks of BEMs were always large matrices, that one had to store in-core or out-ofcore and corresponding large computing times, partly due to computationally expensive integration approaches and partly due to large computing times needed for solution of resulting sets of linear equations. In both cases, considerable saving in computational time and memory requirements are possible by applying subdomain technique to original set of integral equations. The original idea of the subdomain technique is to decompose computational domain into a few relatively large non-overlapping subdomains, interconnected by suitable interface conditions [8]. In comparison with other solution techniques the subdomain technique could be denoted as a direct Domain Decomposition Method (DDM), similar to Schur complement type DDM method [13]. In the context of BDIM this was successfully implemented in [4]. From the numerical approach point of view the technique could be applied to any kind of integral equation however it was shown [4] that in the case of flow kinematics computation it does seriously affect convergence of the coupled velocity–vorticity iteration scheme. This is attributed to the loss of interconnectivity between computational nodes, as the influence of fundamental solution reaches far beyond one domain cell or boundary element. Boundary unknowns, mainly unknown vorticity values on solid boundaries, must therefore be computed from the original one-domain computational mesh.
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4203
Since domain values of velocities are computed explicitly from Eqs. (23) and (24) it is possible to perform computation of these values in a segmentation mesh. The segmentation mesh is composed of several segments, which are geometrically identical to subdomains. However, segmentation technique has to be distinguished from the subdomain technique since segments are independent of each other and no interface conditions have to be imposed. In order to get accurate results the boundary values on segment boundaries have to be computed explicitly from the one-domain mesh. This is the main idea of segmentation technique, presented in [4]. In [4], the unknowns in vorticity transport equation were computed using standard subdomain technique, therefore the computational mesh used for segmentation technique could be used also for the assembly of vorticity transport matrices of integrals. The fundamental solution, used with this approach did not account for any field function values and the numerical scheme was not very stable for higher values of Reynolds number. When the subdomain technique is applied in its limit version, i.e. one subdomain is equal to one internal cell with corresponding boundary elements, the macro-element approach is derived. A macro-element is in general a small subdomain and requires application of implicit interface (I) conditions between neighboring macro-elements (for example 1 and 2), i.e. 1
2
xjI ¼ xjI ; 1 2 ox ox m ¼ m ; on I on I
ð28Þ
in the assembly of the final system of equations. The macro-element approach was used with the flow kinematics velocity–vector form in the BDIM algorithm [5] and tested for various test examples with Dirichlet boundary conditions. The approach was unfortunately limited to regular geometries of computational domains. On the other hand, application of cell or macro-element based diffusion–convection fundamental solution, developed in [5], showed that the numerical stability and accuracy of computational algorithm can be substantially increased with its implementation. Additionally, it is very important that the second main advantage of the macro-element discretization are sparse matrices of integrals and sparse system matrix, generally containing 25 entries in a row. In the computational scheme of flow kinetics (vorticity and heat transport) in [11] all the nodes (domain and boundary) had to be included into the system matrices in order to be able to compute higher Re number value flows, resulting in large fully populated system matrices. When compared with standard fully populated BEM system matrices, the macro-element discretization of flow kinetics is therefore clearly the best possible solution in the construction of a fast computational algorithm. The present paper presents a novel numerical approach to the BDIM which combines segmentation technique for flow kinematics with macro-element approach for flow kinetics i.e. vorticity transport. In this way the algorithm combines conservation properties of vector–potential formulation of flow kinematics and accurate computation of diffusion–convection transport phenomena, performed by macro-element based computation of flow kinetics. The new numerical approach is composed of the following steps: 1. Computation of boundary vorticity and velocity values from flow kinematics equation, discretized by the continuous full domain discretization, Fig. 2 (left), presented here in a final matrix form of system of equations obtained after implementation of boundary conditions and rearrangement of Eqs. (26) and (25): ½AC fxgC ¼ bCþX :
ð29Þ
4204
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
Fig. 2. Sequence of computational meshes.
2. Computation of velocity values on segment boundaries: fvx gS i ¼ ½H t fvy gC ½Dy fxgCþX ;
ð30Þ
fvy gS i ¼ ½H t fvx gC þ ½Dx fxgCþX :
ð31Þ
3. Computation of velocity values inside the segments, Fig. 2 (middle): fvx gXi ¼ ½H t fvy gCi ½Dy fxgCi þXi ;
ð32Þ
fvy gXi ¼ ½H t fvx gCi þ ½Dx fxgCi þXi :
ð33Þ
4. Transformation of velocity and vorticity values from segmentation mesh (Xi) to macro-element mesh (XME): p1 : fvi gXi ! fvi gXME ;
ð34Þ
p2 : fxgC ! fxgCME :
ð35Þ
5. Computation of domain vorticity values and vorticity fluxes from macro-element mesh, presented here in a final matrix form of system of equations obtained after implementation of boundary and interface conditions and rearrangement of Eqs. (27): ½ACME þXME fxgCME þXME ¼ bCME þXME :
ð36Þ
6. Transformation of vorticity values from macro-element mesh (XME) to segmentation mesh (Xi): r : fxgXME ! fxgXi :
ð37Þ
6. Computational aspects of the method 6.1. Inner iteration loop The kinematics and kinetics of the fluid flow motion can be seen as two interlaced problems in a computational procedure. For the general time dependent flows two iteration loops are necessary in order to obtain a solution. The first loop, the time advancing loop or the outer loop, is using known information about the motion at an instant time level for computation of values at subsequent time level, here achieved by the approximation (15). The other iteration loop, the inner loop, serves as a means of computing a divergence free velocity and vorticity field values. In the BDIM context, results of computation of flow kinetics (36) are values of vorticities inside the computational domain, and the result of computation of flow kinematics (29) are vorticity
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4205
values at the boundary of computational domain. Vorticity solution in the whole domain must therefore be searched in an iterative manner, where, with the use of under-relaxation parameter (H), new domain vorticities are computed: iþ1
fxgXME ¼ H iþ1 fxgXME þ ð1 HÞ i fxgXME :
ð38Þ
The inner loop is stopped after the following convergence check for a predefined error (usually 106, Nc = number of computational nodes) is satisfied: PN c iþ1 2 xj i xj Þ j¼1 ð ð39Þ C¼ PN c iþ1 2 : xj Þ j¼1 ð These two iteration loops constitute a computational loop which advances the solution from initial time level to a final time level. 6.2. Recomputation of matrices of integrals of flow kinetics Since the value of velocity vector is changing during each inner iteration formally this would also mean that the entries in matrices of integrals of kinetics must change, as fundamental solution (20) includes the average part of the velocity vector, vj , in a macro-element. This would strongly increase the computational cost of the method therefore the average part of velocity vector is frozen until the relative change in its magnitude has reached a certain level. In practice this means the average velocity in fundamental solution (20) and hence the integrals of a macro-element are updated after the criterion jvj jvj ð40Þ jvj P int is satisfied in the macro-element. Since recomputation of integrals of a macro-element leads to the need of updating the flow kinetics system matrix, this would lead to the increase of overall computation time and could to some extent cancel out positive effects of freezing the velocity in the fundamental solution. To overcome this problem a simple but effective strategy is to implement the ratio of number of macro-elements Nnew, where the criterion (40) was reached, to the overall number of macro-element NME, N new P int N ME
ð41Þ
as the criterion for recomputation of integrals and hence recomputation of system matrix entries. The effect of implementing criteria (40) and (41) can be seen in Fig. 3. Although the full update of integrals of flow kinetics in each inner loop iteration does result in the fastest convergence in the early stage of iterating, the overall convergence for different int is quite similar and leads to the conclusion, that freezing the integrals for several inner iterations does not significantly influence the overall stability of the BDIM iteration procedure. The value of int = 0.05 was therefore chosen and used in all computations in the test examples section. 6.3. Recomputation of system matrix of flow kinetics In order to construct a computationally favourable BDIM algorithm we have to deal also with the process of reformation of system matrix of flow kinetics. The need for the reformation of this matrix comes from the change in integrals but also from including the convection contributions in the entries of the system matrix of flow kinetics, corresponding to unknown vorticity values, i.e.
4206
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
0 -0.5 -1
log(C)
-1.5 epsint=0% epsint=5% epsint=10% epsint=25%
-2 -2.5 -3 -3.5 -4
10
20
30
40
50 60 iteration
70
80
90
100
Fig. 3. Convergence of BDIM scheme for different thresholds (epsint = int) for recomputation of integrals of kinetics. Test case of driven cavity at Re = 100 with non-equidistant 10 · 10 mesh and Dt = 1016.
1 1 ½ACME þXME ¼ ½H þ ½C j ½vj ½Dj ½^vj m m with the corresponding right-hand side ox 1 ½DfxgF 1 : fbg ¼ ½G þ on mDt
ð42Þ
ð43Þ
In this way the system matrix contains the majority of information and the right-hand side of system of equations consists mainly of the known boundary conditions and vorticity contribution from the previous time step. Unfortunately, this means that the flow kinetics system matrix has to be updated in each inner iteration due to the change of vj and hence ^vj . By freezing velocity values of terms in the system matrix, say vj,f, for several inner loop iterations, the system matrix remains constant, 1 1 ½ACME þXME ¼ ½H þ ½C j ½vj;f ½Dj ½vj;f vj m m and only a change in the vector of the right-hand side has to be made: ox 1 1 ½DfxgF 1 : fbg ¼ ½G þ ð½Dj ½C j Þ½vj vj;f fxg þ on m mDt
ð44Þ
ð45Þ
In the case of no freezing, the vj = vj,f and the second term on the right-hand side is zero. In the case of freezing velocities in the system matrix the update of integration is synchronized with the freezing process. An example for this can be seen from comparison of percentage of recomputed integrals of flow kinetics for different numbers of frozen iterations (nfroz), Fig. 4. Selecting larger values of frozen iterations do increase the percentage of newly integrated macro-elements, but the overall number of recomputed integrals is still significantly lower than for the case of update of velocity field for system matrix in each inner iteration (nfroz = 1).
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4207
% of integrated macro elements
100
80
60 nfroz=1 nfroz=10 nfroz=25
40
20
20
40
60
80
100 120 iteration
140
160
180
200
Fig. 4. Comparison of percentage of macro-elements where new integrals were computed based on relative change of 5% in velocity field. Test case as in Fig. 3.
The comparison of convergence of the C error for different numbers of frozen iterations is shown in Fig. 5. Differences in convergence of BDIM in the early stage of inner iteration loop result from the fact that the velocity field, used for the computation of initial integrals of kinetics, represents the solution to the potential flow in the computational domain, as x = 0 is the initial guess and hence the first solution of flow kinematics solves the equation o2 v i ¼ 0: oxj oxj
ð46Þ
0
log(C)
-1
-2
nfroz=1 nfroz=10 nfroz=25
-3
-4
-5 20
40
60
80 100 iteration
120
140
160
180
Fig. 5. Comparison of inner iteration loop convergence for different freezing steps. Test case as in Fig. 3.
4208
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
This is the reason for slower decrease of error C when increasing the number of frozen iterations, as shown in Fig. 5. After the update of velocity field in the integrals of flow kinetics the convergence significantly increases, although in the later stage of inner iteration loop the convergence tends to be similar regardless of the number of frozen iterations. 6.4. Use of diffusion–convection fundamental solution for large cell Re number values The argument of the exponential function in the diffusion–convection fundamental solution (20), pffiffiffiffiffiffiffiffiffi vME AME xe ; ð47Þ 2m where AME is the area of a macro-element, can present numerical problems due to limitations of computer processor in handling very large numbers (overflow). For a given fluid the xe depends on macro-element size and velocity magnitude. In order to avoid computational problems either mesh size or velocity magnitude should be decreased. With a given computational mesh the second possibility is applicable, therefore macro-element based maximum velocity magnitude is defined as 2mxe;max vME;max ¼ pffiffiffiffiffiffiffiffiffi ; AME
ð48Þ
where the xe,max is the largest argument that does not result in an overflow during computation. If the actual magnitude of the average velocity vector in a macro-element exceeds the value (48), local scaling is applied, i.e. vME;max ¼ /MEvME ;
ð49Þ
where /ME is a scaling factor. The scaling is performed simultaneously with recomputation of integrals of flow kinetics, 6.2, and the difference ð1 /ME ÞvME of average part of velocity is treated similar to frozen iteration approach, described in Section 6.2. An interesting observation can be made here. The argument (47) can be equated to half of a macroelement Re number, pffiffiffiffiffiffiffiffiffi vME AME 1 ¼ ReME ð50Þ 2 2m pffiffiffiffiffiffiffiffiffi with characteristic length L ¼ AME and characteristic velocity v ¼ vME applied in definition of ReME. For the AMD based 32-bit personal computer, on which all the computations were made, the maximum macroelement Re number was ReME = 2 · 308 = 616. 6.5. Iterative solvers for solution of resulting systems of equations To lower computational cost of the method, iterative Krylov subspace type solvers are used for the solution of resulting systems of equations: • diagonally preconditioned CGS or GMRES solver for flow kinematics system of equations [4], • level based ILU preconditioned CGS or GMRES solver for flow kinetics system of equations [7]. When compared with block-oriented ILU preconditioning, used in [4], the macro-element discretization allows more efficient use of level based ILU preconditioned, developed in [7]. Freezing the velocities in system matrix of kinetics has a consequence of freezing the preconditioned, so the same preconditioned is used
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4209
for computation of all frozen iterations. Thus, freezing of velocity field in integrals and system matrix of kinetics greatly decreases computational cost of the proposed algorithm. The developed BDIM algorithm has also lower memory demands as the original method [4], which is due to the macro-element discretization of flow kinetics instead of using large subdomains (i.e. segments) as done in [4]. There are further parts of BDIM that need improvements regarding memory demands, especially flow kinematics with its fully occupied system matrix. Wavelet transform techniques present here an effective possibility, a work that is already under development for the BDIM.
7. Test problems The presence of non-smooth boundaries can significantly affect the fluid flow conditions as compared with the flows in confined rectangular cavities. Abrupt changes in velocity profiles, caused by geometric irregularities, are connected with large values of vorticity, therefore an accurate determination of boundary vorticity values is a key step towards accurate solution of the governing set of equations. Except for the first test case, the driven cavity problem, the selected test cases have abrupt changes in geometry. In order to accurately resolve high gradients of velocities and vorticities, highly stretched computational meshes have to be used [19] and the following examples show the application of such meshes in BDIM. 7.1. Flow in a square driven cavity Fluid flow in a square cavity is one of the most difficult benchmark problems in flows of incompressible fluids, which is currently used for testing different numerical schemes for solution of Navier–Stokes equations [14–17,20]. The imposed velocity field at the top moving wall of the cavity drives a large recirculation region inside the cavity. With increasing moving wall velocity and hence the Re number value, additional smaller recirculation zones appear in the corners of the cavity. The complexity of the flow field in connection with a simple geometry makes this problem ideal for testing CFD codes. On the other side, the driven cavity problem can be found in a variety of important engineering applications, like flows in mixing vessels equipped with disk turbines, flows in polymer extruders and environmental flows. The test case was computed on two mesh densities, Fig. 6. The coarse mesh consisted of 30 · 30 macroelements with ratio 4 between the shortest and the longest element side, and 4 large segments for computation of flow kinematics. The fine mesh consisted of 40 · 40 macro-elements with ratio 20, and also 4 segments. The comparison of results between the two meshes showed that the coarse mesh gave accurate results up to Re = 5000, but was unable to give a converged solution for Re = 10,000. For Re = 10,000, the fine mesh was used. A check of mesh convergence with a 60 · 60 ratio 20 mesh with 4 subdomains
Fig. 6. Computational meshes for driven cavity: 30 · 30 with ratio 4 (left) and 40 · 40 macro-elements with ratio 20 (right).
4210
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
showed no improvement in the quality of results, therefore the 40 · 40 mesh was used for all computations of Re = 5000 and higher. All computations were performed by time stepping procedure, where the typical time steps were Dt = 1 for Re 1000 and Re = 5000 and Dt = 0.1 for Re = 7500 and Re = 10,000 (Figs. 7–13). For computation of
Fig. 7. Streamlines and vorticity isolines for Re = 1000.
Fig. 8. Streamlines and vorticity isolines for Re = 5000.
Fig. 9. Streamlines and vorticity isolines for Re = 7500.
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4211
Fig. 10. Streamlines and vorticity isolines for Re = 10,000.
higher value Re number flows starting values of lower Re number value computations were used. Underrelaxation parameter was 0.01 for all the computed cases. In Table 1 a comparison of primary vortex strength values with results of other authors is presented. 7.2. Flow in L shaped driven cavity In the flow in L shaped driven cavity case, the motion of fluid in the cavity is caused by imposed velocity field at the top of the cavity, which, due to the action of viscous forces, drives the recirculating flows inside the lower right and upper left cavity. In contrast to the standard driven cavity flow, the presence of a step inside the cavity causes the flow field to form three large recirculation zones, two in the upper part of the cavity and one in the lower, smaller part of the cavity (Fig. 14). The coarse computational mesh consisted of 300 macro-elements and the dense one consisted of 1200 (three times 20 · 20) macro-elements. In both cases 3 segments were used for segmentation of flow kinematics (Fig. 15).
1
Ghia Re=1000 BDIM Re=1000 Ghia Re=5000 BDIM Re=5000 Ghia Re=7500 BDIM Re=7500 Ghia Re=10,000 BDIM Re=10,000
0.8
y
0.6
0.4
0.2
0 -0.6
-0.4
-0.2
0
0.2 Vx(0.5,y)
0.4
0.6
Fig. 11. Comparison of vx velocity profiles at x = 0.5.
0.8
1
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4212
0.6
0.4
Vy(0.5,y)
0.2
0 Ghia Re=1000 BDIM Re=1000 Ghia Re=5000 BDIM Re=5000 Ghia Re=7500 BDIM Re=7500 Ghia Re=10,000 BDIM Re=10,000
-0.2
-0.4
-0.6 0
0.2
0.4
0.6
0.8
1
x
Fig. 12. Comparison of vy velocity profiles at y = 0.5.
200
0
-200
-400 Ghia Re=1000 BDIM Re=1000 Ghia Re=5000 BDIM Re=5000 Ghia Re=7500 BDIM Re=7500 Ghia Re=10,000 BDIM Re=10,000
y
-600
-800
-1000
-1200
-1400
-1600 0
0.2
0.4
0.6
0.8
1
w(x,1.0)
Fig. 13. Comparison of vorticity profile at the moving boundary.
The computations were performed for Re number values in the range 1000–7500 (Figs. 16–19). Here, the Reynolds number value was defined as (h stands for the total height of the cavity wall), Re ¼
vtop h : m
ð51Þ
As seen from the computation results the irregular geometry of the cavity causes formation of three large vortices (Table 1). The main vortex is at the right top part of the cavity and generally exhibits the characteristics of the main cavity formed in the driven cavity test case except for the lack of smaller corner recirculation regions, which is due to extension of the cavity to the left and to the bottom from the sharp corner.
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4213
Table 1 Comparison of primary vortex strength values Method
Re number value 1000
5000
7500
10,000
Wmin, present, BDIM Wmin [14] FDM Wmin [17] FEM Wmin [9] Poly-BEM Wmin [15] FDM Wmin [3] BEM–BDIM
0.1187 0.1179 0.1189 0.1193 0.1157 0.113
0.1199 0.1189 – 0.1221 0.1123 –
0.1245 0.1199 0.1224 – 0.1052 –
0.1307 0.1197 0.1224 – – 0.109
FDM—finite difference method, FEM—finite element method, BEM—boundary element method, BDIM—boundary domain integral method, Poly-BEM—poly-region boundary element method.
Fig. 14. Presentation of L shaped driven cavity problem.
Fig. 15. Computational mesh for L shaped driven cavity: 300 macro-elements with ratio 4 (left) and 1200 macro-elements with ratio 10 (right).
As a consequence of this large vortex two additional vortices appear in the computational domain, one left and the other at the bottom of the main vortex. They are rotating in the direction opposite to the main vortex.
4214
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
Fig. 16. Streamlines and vorticity isolines for Re = 1000.
Fig. 17. Streamlines and vorticity isolines for Re = 3200.
Fig. 18. Streamlines and vorticity isolines for Re = 5000.
Interesting flow behaviour can be observed on the top of the left vortex, where a thin stratified layer between the moving top wall and the left vortex is formed. It is becoming thinner as the value of Re number is increasing and more and more it exhibits characteristics of a semi closed recirculation region, identified also in the vorticity value increase in the middle of the moving top wall, Fig. 22.
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4215
Fig. 19. Streamlines and vorticity isolines for Re = 7500.
1
0.8
y
0.6
Re=1000 Re=3200 Re=5000 Re=7500
0.4
0.2
0 -0.6
-0.4
-0.2
0
0.2 Vx(0.75,y)
0.4
0.6
0.8
1
Fig. 20. Comparison of vx velocity profiles at x = 0.75.
Each of the additional large vortices has its own recirculation regions at the corners, the left vortex one and the bottom vortex two recirculation regions. They are becoming larger as the Re number value is increased. In the case of the bottom vortex results for Re = 7500 show that a new recirculation region is beginning to form on the right wall between the main and the bottom vortex. The comparison of velocity profiles along the lines x = 0.75 and y = 0.75, presented in Figs. 20 and 21 shows that the main characteristics of the driven cavity flow problem are replicated also in all of the three vortices: almost linear velocity profiles and constant vorticity values inside the vortex regions together with large velocity gradients and hence large vorticity values at vortex boundaries. 7.3. Flow in a Z shaped channel The last test example deals with the flow in a channel that has two bends and the total geometry has a Z-like form. The test case is interesting because of the sudden expansion over the first bend and sudden
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4216
0.6
0.4
Vy(x,0.75)
0.2
0
-0.2 Re=1000 Re=3200 Re=5000 Re=7500
-0.4
-0.6
-0.8 0
0.2
0.4
0.6
0.8
1
x
Fig. 21. Comparison of vy velocity profiles at y = 0.75.
0
-100
vorticity
-200
-300 Re=1000 Re=3200 Re=5000 Re=7500
-400
-500
-600 0
0.2
0.4
0.6
0.8
1
x
Fig. 22. Comparison of vorticity profiles at the moving boundary.
contraction over the second bend which influence the main flow path. Both abrupt changes in the geometry have ratio 2. The test case is interesting also because of mixed type boundary conditions, as the flow conditions at the outlet boundary are not known in advance and are a result of the numerical simulation. The geometry of the channel is presented in Fig. 23. The dimensions correspond to those of [18]. The Re number is defined as Re ¼
vh ; m
ð52Þ
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4217
Fig. 23. Z shaped channel geometry, boundary conditions.
where v is the average velocity at the inlet of the channel and h is the height of the channel at its inlet. In our case v ¼ 1 and h = 1 was chosen. Boundary conditions were prescribed as follows: at the inlet a fully developed laminar flow profile was prescribed, at the outlet the developed flow was assumed by vt = 0 and ox ¼ 0, and the rest of boundaries on was prescribed as no-slip solid wall, i.e. vx = 0, vy = 0. The result of applying the described outlet boundary conditions leads to mixed type boundary conditions in flow kinematics as well as in flow kinetics. In flow kinematics, vorticity is unknown all over the boundary C except for the outlet, where velocities normal to the outlet boundary are unknowns. In vorticity transport equation, the unknowns at the boundary are vorticity fluxes except for the outlet, where vorticities are unknowns. This combination of boundary conditions is another test for the numerical stability of the proposed computational method. The computational domain consisted of 896 macro-elements and seven segments (I–VII, Fig. 23). Computations were performed for Re = 100, 200 and 400. The time step and the under-relaxation values in computations were 0.25 and 0.01, respectively. Computation of higher value Re number flows was done in the same way as in all other test cases, i.e. by using results of the lower Re number value flow as initial conditions for the next higher Re number value. Computational results, presented in Figs. 24 and 25 show that three main recirculation regions are formed regardless of the value of Reynolds number values under consideration. The first one is in the upper right corner, the second is in the lower left corner and the third is formed behind the second bend on the upper wall of the contracted region. For the case of the highest computed Re number value, Re = 400, additional recirculation regions occur. The first additional region is formed along the left wall of the expanded region and is a consequence of the stronger recirculation in the lower left corner. The length of this region is clearly visible from Fig. 26, where vorticity values on the solid wall at x = 2.0 are plotted for various Re number values. The second rather small additional recirculation region appears behind the large recirculation region in the lower left corner. From Fig. 26 one can see that this region is starting to form already at Re = 200, although it becomes pronounced at Re = 400. The length of this region is 0.3 units for Re = 400 and 0.19 units for Re = 200. Also, high vorticity gradients at the sharp corner of the expansion can be observed, although because of the quadratic interpolation functions used these values show oscillations in the vicinity of the corner. The third additional recirculation region is a consequence of redirection of the mean flow after it hits the lower solid wall at the beginning of the contracted region. This region forms on the lower solid wall of the contracted region and is comparable in length with the recirculation region behind the contraction corner. The length of this region for Re = 400 is 1.73 units.
4218
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
Fig. 24. Z shaped channel flow at Re = 200: streamlines (top) and vorticity isolines (bottom).
Fig. 25. Z shaped channel flow at Re = 400: streamlines (top) and vorticity isolines (bottom).
Results, obtained for Re = 200 are in good agreement with results in [18], where QUICK finite differencing numerical scheme was used.
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
4219
4.5
4
3.5
3
Re=100 Re=200 Re=400
y
2.5
2
1.5
1
0.5
0 0
2
4
6
8
10
vorticity on solid wall at (2.0,y)
Fig. 26. Comparison of vorticity x values along the solid wall at x = 2.0.
8. Conclusions A new computational algorithm for the numerical solution of incompressible fluid flow on basis of BEMs has been presented. The algorithm is based on the BDIM and combines two different approximation techniques, approximation of flow kinematics by means of boundary elements and internal cells and approximation of flow kinetics by means of macro-elements. The first ensures divergence free numerical solution and the latter, in connection with the diffusion–convection fundamental solution, delivers accurate solutions to diffusion–convection transport phenomena at higher values of Re number. In order to decrease computational cost of the method, segmentation technique is applied in the flow kinematics computation. Computed test examples of flow in complex geometries show that the developed algorithm is capable of accurate solutions of high Re value flows. Since the BDIM presents an extension of the BEM it is straightforward to couple the method with the established BEM methods and to make use of the advantages BEM has to offer.
References [1] J.C. Wu, J.F. Thompson, Numerical solution of time dependent incompressible Navier–Stokes equations using an integrodifferential formulation, Comput. Fluids 1 (1973) 197–215. [2] L. Skerget, A. Alujevic, C.A. Brebbia, G. Kuhn, Natural and forced convection simulation using the velocity–vorticity approach, in: C.A. Brebbia (Ed.), Topics in Boundary Element Research, vol. 5, Springer-Verlag, Berlin, 1989 (Chapter 4). [3] Z. Rek, L. Skerget, Boundary element method for steady 2D high-Reynolds-number flow, Int. J. Numer. Methods Fluids 19 (4) (1994) 343–362. [4] M. Hribersˇek, L. Skerget, Iterative methods in solving Navier–Stokes equations by the boundary element method, Int. J. Numer. Methods Engrg. 39 (1996) 115–139. [5] L. Sˇkerget, M. Hribersˇek, G. Kuhn, Computational fluid dynamics by boundary domain integral method, Int. J. Numer. Methods Engrg. 46 (1999) 1291–1311. [6] M. Ramsak, L. Sˇkerget, Mixed boundary elements for laminar flows, Int. J. Numer. Methods Fluids 31 (1999) 861–877. [7] M. Hribersˇek, L. Sˇkerget, Fast boundary-domain integral algorithm for computation of incompressible fluid flow problems, Int. J. Numer. Methods Fluids 31 (1999) 891–907. [8] L.C. Wrobel, Applications in thermo-fluids and acoustics. The Boundary Element Method, vol. 1, Wiley, New York, 2002.
4220
M. Hribersˇek, L. Sˇkerget / Comput. Methods Appl. Mech. Engrg. 194 (2005) 4196–4220
[9] M.M. Grigoriev, G.F. Dargush, A poly-region boundary element method for incompressible viscous fluid flows, Int. J. Numer. Methods Engrg. 46 (1999) 1127–1158. [10] W.F. Florez, H. Power, Comparison between continuous and discontinuous boundary elements in the multidomain dual reciprocity method for the solution of the two-dimensional Navier–Stokes equations, Engrg. Anal. Bound. Elem. 25 (2001) 57–69. ˇ unicˇ, Natural convection flows in complex cavities by BEM, Int. J. Numer. Methods Heat Fluid [11] L. Sˇkerget, M. Hribersˇek, Z. Z Flow 13 (5/6) (2003) 720–736. [12] M.S. Ingber, A vorticity method for the solution of natural convection flows in enclosures, Int. J. Numer. Methods Heat Fluid Flow 13 (5/6) (2003) 655–671. [13] W. Hackbusch, Iterative Lo¨sung Grosser Schwachbesetzter Gleichungssysteme, B.G. Teubner, Stuttgart, 1991. [14] U. Ghia, K.N. Ghia, C.T. Shin, High-Re solutions for incompressible flow using the Navier–Stokes equations and a multigrid method, J. Comput. Phys. 48 (1982) 387–411. [15] O. Goyon, High-Reynolds number solutions of Navier–Stokes equations using incremental unknowns, Comput. Methods Appl. Mech. Engrg. 130 (1996) 319–335. [16] C.H. Liu, D.Y.C. Leung, Development of a finite element solution for the unsteady Navier–Stokes equations using projection method and fractional H-scheme, Comput. Methods Appl. Mech. Engrg. 190 (2001) 4301–4317. [17] E. Barragy, G.F. Carey, Stream function–vorticity driven cavity solution using p finite elements, Comput. Fluids 26 (5) (1997) 453–468. [18] C.Y. Perng, R.L. Street, A coupled multigrid-domain-splitting technique for simulating incompressible flows in geometrically complex domains, Int. J. Numer. Methods Fluids 13 (1991) 269–286. [19] P. Sockol, Multigrid solution of the Navier–Stokes equations on highly stretched grids, Int. J. Numer. Methods Engrg. 17 (1993) 543–566. [20] S. Turek, Tools for simulating non-stationary incompressible flow via discretely divergence free finite element models, Int. J. Numer. Methods Fluids 18 (1994) 71–105.