Multiblock hybrid grid finite volume method to solve flow in irregular geometries

Comput. Methods Appl. Mech. Engrg. 196 (2006) 321–336 www.elsevier.com/locate/cma

Multiblock hybrid grid finite volume method to solve flow in irregular geometries Masoud Darbandi *, Alireza Naderi Department of Aerospace Engineering, Sharif University of Technology, P.O. Box 11365-8639, Tehran, Iran Received 25 August 2004; received in revised form 22 December 2005; accepted 11 April 2006

Abstract In this work, a finite-volume-based finite-element method is suitably developed for solving incompressible flow and heat transfer on collocated hybrid grid topologies. The method is generally applicable to arbitrarily shaped elements and orientations and, thus, challenges the potential to unify many of the different grid topologies into a single formulation. The key point in this formulation is the correct estimation of the convective and diffusive fluxes at the cell faces using a novel physical influence scheme. This scheme remarkably enhances the achieved solution accuracy. It is shown that the extended formulation is robust enough to treat any combination of multiblock meshes with dual element shape employments. In this regard, the solution domain is broken up into a number of different multiblock arrangements in which each block is filled with only one type of finite element shape. The combined grid can decrease the computational time and memory requirements and increase the numerical accuracy. The results of the extended formulation are validated against different benchmark and other solutions. The current results are in excellent agreement with the other solutions without exhibiting any disturbances around the block boundaries.  2006 Elsevier B.V. All rights reserved. Keywords: Finite volume; Finite element; Hybrid grid; Multiblock; Incompressible flow; Upwinding

1. Introduction The efficient numerical solution of the fluid flow governing equations is not achieved unless suitable grid arrangements are constructed and utilized. Different grid generation techniques have become very encouraging because of the demand on quick and efficient generation of quality grids in complex geometries [61]. The complexity in either solution domain geometry or fluid flow behavior promotes the research toward constructing and employing different grid topologies and orientations in a single flow problem. Subsequently, the adoption of this strategy may reduce the need for employing either very fine grid distributions or unstructured grid topologies. It is known that the unstructured grids normally require greater computational

*

Corresponding author. Tel.: +98 216 616 4944; fax: +98 216 602 2731. E-mail address: [email protected] (M. Darbandi).

0045-7825/$ - see front matter  2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2006.04.005

cost per grid point. Alternatively, one suggestion is to utilize mixed element shapes in complex geometries in order to achieve solutions with higher accuracies within reasonable solution periods. Indeed, the importance of mixed element shapes and meshes has caused motivations in side fields such as mesh generation activities [38,61]. There are a few papers in the relevant literature which report the benefits of using mixed grid topologies in different aspects of numerical methodologies and/or their applications. Davidson [17] eliminates the grid imperfection, which is caused by utilizing a single grid topology, by employing irregular combinations of multiblock meshes filled with either quadrilaterals or triangles. Beskok and Warburton [4] use modal expansion of Jacobi polynomials in mixed triangular and quadrilateral elements. Essentially, in a two-dimensional mesh refinement study, quadrilaterals are easily split into smaller quadrilaterals or even triangles, while triangles are only split into smaller triangles. However, the essence of the preceding work has been to

M. Darbandi, A. Naderi / Comput. Methods Appl. Mech. Engrg. 196 (2006) 321–336

distribute structured quadrilateral grids inside each triangle of an unstructured grid topology. Ainsworth and Coyle [1] construct a set of hierarchic base functions for the Galerkin discretization of the space H suitable for hybrid meshes containing both quadrilateral and triangular elements with arbitrary non-uniform order polynomial approximation. From the solution algorithm perspective, Toreja and Rizwan-uddin [56] develop the hybrid nodal-integral finite element (and analytic) method to solve the steadystate convection–diffusion equation in arbitrary geometries using both rectangular and triangular nodes with the latter restricted to boundaries that are not parallel to the utilized Cartesian coordinates. Moinier et al. [44] extend the employment of a multigrid approach to hybrid grid context. They construct a sequence of coarse hybrid grids on which the same residual operator can be applicable. Their hybrid grids involve a structured quadrilateral distribution with elements stretched in the vicinity of the airfoil and along it. The rest of the domain is filled with triangles. In a 3D attempt, Liou and Zheng [43] generate composite grids using an unstructured triangular-based grid close to the row of blades in a cascade and a structured quadrilateral-based distribution far from the blades and also in the channel between the vanes. Alternatively, they use a background H-type grid placed to cover the channel between the vanes, an O-type viscous grid used to resolve the region around the vane, and an unstructured grid located in the overlapping region between the H- and O-type grids. From the element shape perspective, Feistauer et al. [23] utilize a hybrid grid, which consists of triangular and dual meshes. A dual finite volume mesh can be constructed from a triangular mesh by the proper connection of the medians of the adjacent triangles in a solution domain filled only with triangular elements [24]. Ref. [23] uses a structured triangular mesh in the vicinity of the viscous walls while the rest of domain is filled with a dual mesh. From the grid adaption perspective, Struijs et al. [55] and Khawaja et al. [34] demonstrate the benefits of using dual cells for adaptive griding where triangular and quadrilateral elements can be easily combined to allow structured layers such as boundary layers near solid walls. Hwang and Wu [31] also present an adaptive cell-vertex finite volume approach on unstructured mixed quadrilateral–triangular meshes. In their approach, the quadrilaterals are approximately stretched in the flow regions having one-dimensional features, e.g., an oblique shock wave reflected from a wall. The advantages of employing hybrid grid topologies have promoted the current authors to extend a finitevolume-based finite-element method [9,13] to fluid flow and heat transfer applications on multiblock topologies using mixed triangular and quadrilateral elements. This work provides the required modifications which are necessary to be carried out in order to include the advantages of utilizing mixed element shapes in a finite volume context. To perform the general applicability of the current formulation on different grid topologies, the solution domain is suitably broken into different multiblock structures. It is

shown that the accuracy of the current numerical solution is excellent despite using coarse grid distributions. Additionally, in order to prove that the current methodology is not restricted to multiblock mesh employments, a scattered distribution of the two element shapes is also constructed and examined in a single-block solution domain. The numerical solution shows that the extended formulation is robust enough to treat any combination of multiblock meshes and/or single element shapes. The investigation also shows that the current solution does not show disturbances when passing through the blocks boundaries. 2. Domain discretization In finite element approaches, the solution domains are initially broken into infinitesimal elements. In two dimensions, there are two basic choices to select the element shape, i.e., triangular and quadrilateral shapes. In the literature, it is often observed that the finite element method normally utilizes only one type of element shape to discretize the solution domain. As was mentioned in Section 1, there are papers which alternatively address the use of hybrid grids. Generally speaking, the grid generation in irregular geometries can be remarkably simplified if the domain is properly divided into several smaller blocks. Consequently, each block can be filled with only one type of element which suitably fits the boundaries of that block [17,43,61]. However, the division of a domain into subdomains and the choice of suitable grid topologies in each subdomain are only one perspective of the problem in hand. The other perspective is to extend the numerical procedure in a manner which does not degrade its original robustness and superiorities. Fig. 1 shows a triangular cavity broken into four subdomains or blocks. As is seen, three

Grid=32 (on top face)

0.9 0.8 0.7 0.6 0.5

Y

322

0.4 0.3 0.2 0.1 0 -0.5

-0.25

0

0.25

0.5

X Fig. 1. A triangular cavity broken into four mesh blocks and filled with two types of element shapes.

M. Darbandi, A. Naderi / Comput. Methods Appl. Mech. Engrg. 196 (2006) 321–336

blocks are filled using triangles with different orientations and the fourth is filled with quadrilaterals. This type of grid distribution suppresses the adverse grid orientation effect [15]. Fig. 2 demonstrates an arbitrary distribution of quadrilateral and triangular elements in a part of solution domain filled with hybrid grids. The solid lines show the element boundaries. Fig. 3 shows two individual triangular and quadrilateral elements. After a suitable distribution of the elements within the solution domain, each element is broken into a number of sub-control volumes (SCV). A few of these SCVs are shown in Fig. 2. If these SCVs are properly assembled, the control volumes can be constructed. Fig. 3 shows how an element can be suitably broken into

SS

8 SS+

7

+

SS6

9 SS

1

+

+

SCV2

+ SS

+ SCV3

SS4

SCV1 +ip1 ip2+

+ S3 S

2 SS

ds

ds

Fig. 2. A control volume constructed from the assemblage of several subcontrol volumes (SCVs).

several SCVs. As is seen, each quadrilateral element is divided into four SCVs by drawing lines between the midpoints on the opposite edges, see Fig. 3(a). These lines are shown by dashed lines and labelled SS1 to SS4. On the other hand, each triangle can be divided into three SCVs by the help of its three medians, see Fig. 3(b). The medians are also shown by dashed lines and labelled SS1 to SS3. At this point, it is important to notify that the labelling of SCVs and SSs in Fig. 3 is different than that in Fig. 2. In fact, Fig. 2 labels them within a cell while Fig. 3 labels them within an element. Fig. 2 shows that the proper assemblage of the neighboring SCVs around one arbitrary node results in a control volume which surrounds that node. The shaded area in Fig. 2 presents a control volume consisted of five SCVs from five different elements. The dashed lines around the shaded area represent the cell faces. There are 10 cell faces around the shaded cell. They are labelled SS1–SS10. In the current formulations, node is treated as the location of all problem unknowns. This type of grid arrangement is called collocated grid because all problem unknowns are located at one sort of nodes. The geometric simplicity of the collocated grid arrangement is very attractive and it will be significant if the pressure checkerboard problem can be suppressed [45]. This problem is further discussed in Section 4. In Section 4, it is necessary to integrate the governing equations. The integrations are taken over the cell faces, i.e., SS1–SS10 faces around the shaded cell in Fig. 2. Such integrals are approximated by evaluating their arguments at the midpoint of cell faces. These midpoints are denoted as integration points ip and are shown in Figs. 2 and 3 by plus symbols. Two of them are labelled as ip1 and ip2 in Fig. 2. They surround SCV1. At this stage, it is interesting to notify that there are four integration points within each quadrilateral element, hence, a pure quadrilateral grid requires to evaluate 4Nel convective (or diffusive) fluxes where Nel denotes the number of

2

2 1 SS1

L up up

ip1+

SCV1 SCV2

vip1

dn

Ldn

b

SS 2

L up

up

ip1+

b

ip2 +

SS2

ip4 SS 4 +

1

SCV2

SS

SS10

+ SCV5

SS5

SCV4

323

ip2 +

vip1 L dn

SCV1

a

SCV3

SCV4

3

4

(a) Quadrilateral element.

SS3

SS3

+ip3

a

+ ip3

SCV3

3

(b)Triangular element.

Fig. 3. The nomenclature and velocity upwinding in utilized finite elements.

dn

1

324

M. Darbandi, A. Naderi / Comput. Methods Appl. Mech. Engrg. 196 (2006) 321–336

elements in the domain. Alternatively, if the same quadrilateral is broken up into two triangles, it requires 6Nel flux term evaluations. It is because there are two triangles in each quadrilateral and three integration points in each triangle. Therefore, the use of more quadrilaterals in the grid is computationally more desirable. In another words, the hybrid grid strategy can remarkably improve the cost per grid point in our formulation [13]. Before closing this section, it is necessary to mention that we use upper case letters such as U, V, and P to denote the magnitudes at nodes or cell centers, which are the location of all problem unknowns. On the other hand, the magnitudes at cell faces are denoted by the lower case letters, e.g., u, v, and p. The secondary unknown variables at the cell faces cannot be included in the list of primary unknown dependent variables because the number of unknowns overwhelms the number of equations and the problem becomes ill-posed. To make our implicit procedure wellposed, the unknown variables at cell faces need to be expressed in terms of the unknown variables at cell centers. The detail is further discussed in Section 4. 3. The governing equations The two-dimensional steady-state governing equations is given by oF oG oR oT þ ¼ þ þ B; ox oy ox oy

ð1Þ

where (F and G) and (R and T), respectively represent convection and diffusion flux terms. They are defined as T

F ¼ fu; qu2 þ p; quv; quhg ; G ¼ fv; qvu; qv2 þ p; qvhgT ;

ð2Þ T

R ¼ f0; sxx ; sxy ; usxx þ vsxy  qx g ; T ¼ f0; syx ; syy ; usyx þ vsyy  qy gT :

ð3Þ

Using the Boussinesq’s approximation, the buoyant force is given by T

B ¼ f0; q0 gx bðt  t0 Þ; q0 gy bðt  t0 Þ; 0g :

ð4Þ

In the above equations, q is the fluid density, V ¼ u^i þ v^j is the flow velocity, p is the static pressure, h is total enthalpy, t is temperature, and g ¼ gx^i þ gy^j is the gravitational acceleration. Additionally, q0, t0, and b represent the reference density, temperature, and thermal expansion coefficient magnitudes, respectively. The stress tensor and heat flux components are given by   ou ov ou ov þ sxx ¼ 2l ; syy ¼ 2l ; sxy ¼ syx ¼ l ; ð5Þ ox oy oy ox ot ot ð6Þ qx ¼ k ; qy ¼ k ; ox oy where l and k represent the fluid viscosity and thermal conductivity coefficient, respectively.

4. Computational modelling The discretization of the governing equations is carried out using the finite volume method. The integration of Eq. (1) over an arbitrary control volume, e.g., the shaded cell in Fig. 2, and the use of divergence theorem eventually yield Z Z Z Z ðF^i þ G^jÞ  dS ¼ ðR^i þ T^jÞ  dS cell face cell face Z Z Z þ B dt; ð7Þ cell area

where the two first integrals are taken over the cell faces (i.e., SS1–SS10 faces in Fig. 2) and the last integral is taken over the cell area (i.e., the shaded area in Fig. 2). Additionally, dt indicates the differential area. The vector dS is considered as an outward normal vector to the cell face. Fig. 2 shows this vector at two ip1 and ip2 midpoints. This vector is given by dS ¼ ðDS x Þ^i þ ðDS y Þ^j;

ð8Þ

where DSx and DSy represent the components of cell face length in ^j and ^i directions, respectively. Using the definition of dS, the cell face integrals in Eq. (7) can be approximated at the midpoint of cell faces. The last integral in Eq. (7) can be also approximated using a lumped-mass approach. These considerations finally yield imax imax X X ½FðDS x Þ þ GðDS y Þi ¼ ½RðDS x Þ þ TðDS y Þi þ Bt; i¼1

ð9Þ

i¼1

where i denotes the midpoint of cell faces around the chosen control volume. As is seen in Fig. 2, imax = 10 for the shaded cell. Additionally, t is the area of cell. Eq. (9) includes four conservative statements for the mass, momentums, and energy equations. As it is seen in Eq. (9), the buoyant term B is calculated at the cell center while the convection (F and G) and diffusion (R and T) flux terms are required to be calculated at the midpoint of each cell face. Considering U (or u), V (or v), P (or p), and T (or t) as the main dependent variables in our algorithm, the convective flux terms and the energy dissipation terms are non-linear and need to be linearized with respect to the main dependent variables. A simple linearization scheme yields T

F  fu; quu þ p; quv; quhg ; T

ð10Þ

G  fv; qvu; qvv þ p; qvhg ; T

R  f0; sxx ; sxy ; usxx þ vsxy  qx g ; T

T  f0; syx ; syy ; usyx þ vsyy  qy g :

ð11Þ

The bar over the velocity components means that the velocities are approximated from the preceding iteration. In the current work, we solve the governing equations in a steady form and the steady-state solutions are computed through non-linear iterations. The simple linearization given in Eqs. (10) and (11) would suffice to treat the incompressible flow

M. Darbandi, A. Naderi / Comput. Methods Appl. Mech. Engrg. 196 (2006) 321–336

equations. However, Darbandi and Mokarizadeh [10] present a conceptual linearization and linearize the convection terms with respect to both the velocity and mass flux components. For example, the non-linear (quv) term in Eq. (10) becomes ðquv ¼ quv þ vqu  quvÞ. This linearization is suitable for the flow fields with high density variations, i.e., compressible flow fields. As was mentioned in Section 2, the unknown dependent variables at the cell faces (the variables shown by lower case letters such as u and v) cannot be treated as the major dependent variables in our algorithm. In finite volume methods, it is customary to derive expressions which relate the unknown variables at the cell faces to them at the cell centers. Therefore, to connect Fi, Gi, Ri, and Ti in Eq. (9) (which are treated at the midpoints of cell faces) to the main dependent variables at nodes, one may use the finite element shape functions. For example, the pressure at ip3 midpoint in Fig. 3 can be related to the pressure magnitudes at the element vertices using pip3 ¼

N nd X

ðN j Þip3 P j ;

ð12Þ

j¼1

where Nnd counts the number of vertices in the element where ip3 is located. Additionally, the subscripts of (Nj)ip3 indicates that the shape functions are calculated at ip3. Eqs. (5) and (6) indicate that all the terms in R and T, i.e., Eq. (3), are in gradient form. Since these terms are diffusive type, one straightforward suggestion is to treat them using the gradients of the finite element shape functions. For example, the approximation of ou/ox term at an arbitrary cell face i in Fig. 3 yields    N nd  X ou oN j  U j; ð13Þ ox i ox i j¼1 where Nnd indicates the number of element vertices. Choosing a bilinear interpolation, Nnd is three for the triangular and four for the quadrilateral elements. To simplify our notations, we use i instead of ips in the rest of this section. It is important to notify that although there is no restriction in the choice of element order, we have chosen linear shape functions to fulfill this study, see Figs. 2 and 3. Contrary to the terms in R and T, the terms in F and G are not diffusive type and cannot be simply expressed by the finite element shape functions and their gradients. The past efforts show that there have been serious challenges in accurate modelling of convection–diffusion problem. As is known, the standard finite-element methods are unable to simulate the correct behavior of flow when the advection is the dominant part in a mixed advection–diffusion problem [32]. For the past couple of decades, different classes of methods have been developed to overcome the encountered disabilities. One class of methods known as Eulerian method utilizes the primitive modified techniques such as upstream weighting to generate more accurate approximations [46]. Such methods generally suffer from

325

excessive truncation errors. To surmount the inaccuracies, more advanced techniques were constructed, e.g., the Petrov–Galerkin FEM [3], optimal test function method [6] and the streamline diffusion FEM [5,28]. These methods add the numerical diffusion only in the direction of streamline. In fact, they enforce a free parameter which determines the amount of diffusion applied and therefore has a great effect on the achieved accuracy. A second class of methods known as characteristic method makes the use of hyperbolic nature of the convection–diffusion equation. They use a combination of Eulerian fixed grids to treat the diffusion component and Lagrangian coordinates by tracking particles along the characteristics to treat the advective component [41,46,58]. These methods allow the advantage of using large time steps to march in time. However, these methods basically have difficulty in mass conservation and general boundary condition treatment. The Eulerian– Lagrangian localized adjoint methods ELLAM are considered as the extension of the characteristic methods which conserve mass and treat general boundary conditions properly [7,30,54]. Similar to the finite element context, there have been numerous efforts in the finite-difference and finite volume methods in order to incorporate the correct physics of flow in the formulations. The upstream upwinding [50] and the skewed upwinding [51] schemes are two basic remedies which produce less and more false diffusion. To increase the accuracy, higher order interpolators such as quadratic upstream [40,47] have been suggested. Indeed, the history in progress of the cell face flux interpolators has been long started and its research is still open to be worked on [29,37]. Unfortunately, few of the past interpolators consider both the mathematical and physical point of views. Their bases are mainly either physical or mathematical. Alternatively, Darbandi and Schneider [11] suggest new statements which would respect both views. Because of developing an allspeed algorithm, they concentrated on a momentum-based formulation instead of a velocity-based one to treat compressible flow as well. Their formulations include the contribution of all convection, diffusion, and pressure terms in calculating the convective fluxes at the cell faces. In this work, the developed formulation is suitably extended in order to consider the correct transportation of advective fluxes using a velocity-based formulation instead of a momentum-based one. In this regard, the pressure terms in F and G are firstly calculated using the finite element shape functions because of the elliptic nature of pressure field. The approximation at the ith cell face in Fig. 3 yields pi 

N nd X ðN j Þi P j :

ð14Þ

j¼1

As was mentioned, the approximation of the active velocity components in F and G needs some more sophisticated schemes because of their non-elliptic physics. To achieve this, we derive the cell face statement for u velocity from

326

M. Darbandi, A. Naderi / Comput. Methods Appl. Mech. Engrg. 196 (2006) 321–336

the x-momentum governing equation using both finitedifference and finite-element contributions. In this regard, the x-momentum equation can be written in the streamwise direction as ou 1 op þ ¼ mr2 u; ð15Þ os q ox pffiffiffiffiffiffiffiffiffiffiffiffiffiffi where V tot ¼  u2 þ v2 . We do not need including the buoyant term in this equation. It is because Darbandi and Schneider [12] have already shown that the inclusion of this term in Eq. (15) can effectively improve the accuracy of the solution in solving high thermobuoyant flow fields. Such flow fields are out of concern in this study. On the other hand, the inclusion or exclusion of the buoyant term (and even the other terms) in Eq. (15) is just a choice not a force. At this stage, we want to discretize Eq. (15) term-by-term. Indeed, the key point in the current physical-based formulation is the prominent role of convection terms which have been written in the streamwise direction. This form provides the correct direction of upwinding in the streamwise direction, see Fig. 3. Considering an arbitrary cell face i, the employment of a first-order upwinding in the streamwise direction yields     ou u  uup V tot  ðV tot Þi os i Lup i XN nd N nd ðN j Þi U j  ðuup Þi X j¼1  ½ðN j Þi ðV tot Þj  ; ð16Þ ðLup Þi j¼1 V tot

where Lup and uup are illustrated in Fig. 3 for i = ip1 in both quadrilateral and triangular element cases. The upstream locations are found by intersecting the extension of the streamline direction at the ith integration point with the edges of the same element. Then, the magnitude of (uup)i is linearly interpolated between the neighboring nodes at the same edge. For example, the magnitude of uup at ip1 in Fig. 3 is obtained using (u(up))ip1 = (a/b)U2 + (1  a/b)U3. Considering the elliptic nature of pressure field, see Eq. (14), and using a treatment similar to Eq. (13), the pressure term in Eq. (15) is discretized to    N nd  1 op 1X oN j  P j: ð17Þ q ox i q j¼1 ox i The elliptic nature of the diffusion term in Eq. (15) also promotes the use of finite element shape functions. We explicitly calculate the diffusion term from the known magnitudes of the velocity components in the preceding iteration. In another words, the diffusion term is treated as a source term. Nevertheless, we use the divergence theorem and a mass-lumped approach to approximate the Laplacian operator in Eq. (15). It yields Z Z 1 ðr2 uÞi  ðruÞi  dS; ð18Þ te element face where te indicates the area of the element.

The modelling of all terms in Eq. (15) is complete now. If these models, including the models in Eqs. (16)–(18), are plugged into Eq. (15) and similar terms are combined together, an algebraic expression for the u component at the ith integration point within an element is obtained. As is seen in Fig. 3, there are three and four ip’s in triangular and quadrilateral elements, respectively. We can derive the algebraic expressions for the u component at all ips. The resulting expressions show that the cell face u velocity components of an element are in terms of the major dependent variables (i.e., U, V, and P) of the same element. For the sake of brevity, the compact form of all the algebraic expressions can be written as fug ¼ ½C uu fU g þ ½C up fP g þ fC u g:

ð19Þ

Considering an element with a quadrilateral (or triangular) shape, the Cuu and Cup are two 4 · 4 (or 3 · 3) matrices which indicate how the u velocity components at the cell faces within an element can be expressed in terms of U and P variables at the vertices of that element. The 4 · 1 (or 3 · 1) array of Cu includes all known parts of the assembled terms, e.g., the source term given in Eq. (18). To distinguish the matrices from each other, the first superscript of C means that it belongs to u velocity component. The second superscript indicates which array of dependent variables it is multiplied by. For example, Cup indicates that the matrix includes coefficients which are multiplied by the pressure unknown in u velocity expression. Similar to the steps taken to derive an expression for u, an expression can be readily obtained for v using the same procedure; however, starting the procedure from the ymomentum equation. Similar to Eq. (19), the final algebraic expression for the v velocity components can be written as fvg ¼ ½C vv fV g þ ½C vp fP g þ fC v g:

ð20Þ

The u and v cell face expressions given by Eqs. (19) and (20) are subsequently substituted only in the active velocity components of the convective terms in the x and y momentum statements in Eq. (10). By this substitution and the substitutions from Eqs. (13) and (14), the momentum conservative statements become free from u, v, and p which were not considered as our major dependent variables in our algorithm. Therefore, it remains to close the conservative statements for the continuity and energy equations. To solve the pressure checkerboard problem [45] in our pressure-based algorithm, it requires some remedy. One idea is to use the staggered grid arrangement which is out of concern in our collocated formulations. An alternative choice is to use dual velocity definitions at cell faces. This is why we did not substitute Eqs. (19) and (20) in the conservative statement for the continuity equation, see Eq. (10). In this work, collocation of variables at cell centers is achieved using the pioneer idea of Rhie and Chow [52]. However, the current alternative velocities are derived in a completely different manner. In another words, we once more use our proposed physical-based strategy to derive

M. Darbandi, A. Naderi / Comput. Methods Appl. Mech. Engrg. 196 (2006) 321–336

the second set of velocity component expressions at the cell faces. To achieve our purpose, the new expressions can be obtained by a combination of both momentum and continuity equation errors. Darbandi and Bostandoost [8] suggest a combination of both momentum and continuity equation errors, i.e., [(the discretized momentum equation error)  u(the discretized continuity equation error) = 0]. They show that the velocity cell face which is obtained from this expression suffices to suppress completely spurious pressure oscillations in the flow field. Following their one-dimensional investigation, we suggest a suitable expression to derive our new u velocity component. It is written as   ou 1 op ou ov V tot þ u þ ð21Þ ¼ mr2 u: os q ox ox oy

327

tum equations are treated using the bilinear interpolation given in Eq. (14). Back to the conservative statements given in Eq. (9) and the above treatments and substitutions, the final conservative statements for each cell in the domain can be written as 9 9 8 28 3 ^v ^u > > > > imax = = < < X 6 7 q^vu  syx ðDS y Þ5 4 q^uu þ p  sxx ðDS x Þ þ > > > > ; ; : :  i¼1 q^uv  sxy q^vv þ p  syy i 9 8 0 > > = < ð24Þ ¼ t q0 gx bðt  t0 Þ ; > > ; :  q0 gy bðt  t0 Þ

The hat on this velocity component indicates that it differs from that given in Eq. (19). The constructed matrices and vectors and the definitions of their superscripts are as before except using a hat over u. Contrary to Eq. (19), both U and V dependent variables as well as P variable appear in the second cell face velocity expression. Similar to Eq. (22), a second v velocity expression can be derived starting from a suitable combination of the y-momentum and continuity equations. The treatment eventually yields

where t is the volume of cell. As was mentioned before, the buoyant term is approximated using a mass-lumped approach. This is why this term appears in the right-hand-side of Eq. (24). If the conservative statements in Eq. (24) are carefully inspected it is observed that U, V, and P are the only unknown dependent variables which actively remain in the equations. The achieved statements are linear with respect to our dependent variables. Fortunately, all the dependent variables appear in all the conservative statements. This is considered as a big superiority of our formulation because it results in a strong coupling among U, V and P dependent variables in our system of linear algebraic equations. If we obtain Eq. (24) for all the cells in the domain and implement the required boundary conditions, a system of linear algebraic equations is obtained. We can readily solve this system for the pressure and velocity fields. The last stage is to close the conservative statement for the energy equation. As is indicated by Eq. (24), the energy equation is not simultaneously solved with the flow equations. To close the energy equation, the non-linear convection terms in the energy equation, see Eq. (10), are further linearized with respect to the temperature variable. This results in

f^vg ¼ ½C^vu fU g þ ½C^vv fV g þ ½C^vp fP g þ fC^v g:

quh  qcp ut þ quV 2tot =2;

This equation is a combination of Eq. (15) and the continuity equation. All of the terms in the momentum part of this equation are treated similar to those performed for treating Eq. (15). However, the two other continuity terms in the brackets are discretized using the approximation given in Eq. (13). These treatments finally result in a new u velocity expression at the cell face. Similar to Eq. (19), we can compact all the u cell face velocities located in an element in vector form, i.e., f^ ug ¼ ½C ^uu fU g þ ½C ^uv fV g þ ½C ^up fP g þ fC ^u g:

ð22Þ

ð23Þ

The velocity statements given in Eqs. (22) and (23) are plugged into the conservative statement of the continuity equation, see Eq. (10). The experience shows that these substitutions would fully suppress the non-physical oscillations in the pressure and velocity fields [8]. At this stage, we are ready to construct the system of linear algebraic equations. In this regard, the velocity components given in Eqs. (22) and (23) are plugged in the continuity equation. Since the pressure variable appears in both ^ u and ^v expressions, it results in achieving one major advancement toward coupling of pressure and velocity fields in our fully implicit algorithm. The foregoing velocities are also plugged in the lagged velocities in the momentum equations. On the other hand, the velocity components given in Eqs. (19) and (20) are plugged in the convective terms of the x and y momentum equations, respectively. All the stress term components in the momentum equations are treated suitably using the tactic given in Eq. (13). Additionally, the pressure terms in the momen-

qvh  qcp vt þ qvV 2tot =2:

ð25Þ

To obtain suitable cell face temperature expression, we use the differential form of the energy equation which has been written in the streamwise direction. If the dissipation terms of the energy equation are ignored, the equation yields V tot

ot k 2 1 ~  V~Þ: ¼ r tþ ðpr os qcv qcv

ð26Þ

Comparing Eq. (26) with Eq. (15), it shows that their convection and diffusion terms are the same if m is changed with k/(qcv) and u with t. Therefore, these two terms can be treated in manners similar to those performed in Eq. (15). For example, the approximation of the diffusion term is done similar to the one in Eq. (18). Additionally, the pressure term in Eq. (26) is approximated by N nd X ~  V~Þ  ðr ~  V~Þ ðpr ðN j Þi P j ; ð27Þ i i j¼1

where all the velocity gradients are calculated using the strategy given in Eq. (13). Substituting the above

328

M. Darbandi, A. Naderi / Comput. Methods Appl. Mech. Engrg. 196 (2006) 321–336

approximations into Eq. (26) and writing Eq. (26) for all the midpoints located in an element. They finally result in ftg ¼ ½C tt fT g þ fC t g:

ð28Þ

The definitions of C and its superscripts are similar to those given in Eq. (19) except that the current matrix and vector are 1 · 3 (or 1 · 4) and 1 · 1, respectively. This cell face temperature expression can be plugged in Eq. (25) and the results can be plugged in the conservative statement for the energy equation, Eq. (9). Similar to Eq. (24), the energy conservative statement for an arbitrary cell in the domain is given by imax X ^vtgðDS y Þ ^ ½fqcp  utgðDS x Þ þ fqcp  i¼1

  sxx þ vsxy  qx  q uV 2tot =2 ðDS x Þ ¼ u   þ  usyx þ vsyy  qy  qvV 2tot =2 ðDS y Þi :

ð29Þ

Considering the minor role of conduction in our study, qx and qy terms have been lagged and transferred to the right-hand-side. If Eq. (29) is written for all the cells in the domain and the required boundary conditions are implemented, a system of linear algebraic equations is obtained. This system of equations can be solved for the temperature field T. After updating the temperature field, the flow field is updated. This sequence is continued until converging to a solution with sufficient accuracy. A careful investigation in the extended formulation shows that the pressure field has a second-order accuracy. Alternatively, the velocity field has a first-order accuracy due to streamwise upwinding utilized in Eqs. (15) and (21). Despite low order of the velocity field, the current formulation performs a rate of convergence higher than the first-order upwinding because of the adopted reconstructions [16]. This can be clearly seen in the Results section where very accurate solutions are obtained despite using coarse grid resolutions. 5. Results There are three important aspects need to be respected in developing a multiblock-based scheme. The scheme is expected to be stable, conservative, and accurate around the block boundaries irrespective of the type of grids in the neighboring blocks and their connections at the joint boundaries. We present a few of works concentrated on these aspects. From the stability perspective, de Nicola, et al. modelled quasi-one-dimensional [18] and twodimensional convection diffusion [19] Euler equations and showed that a block interface could cause serious instability in the numerical solution of multiple-zone schemes [2]. Wood and Kleb [59] studied the effect of grid change and its orientation in the oscillatory behavior of the numerical solution. They showed that the change in the shape of grid elements could result in oscillatory behavior in the finitevolume solutions. They also reported that the convergence

time could increase by a factor of 2.7 when the oscillatory solution appeared in the domain. Since we are using finite element volume method, the role of finite element shape functions is very important in the solution. At the block boundaries, the shape functions vary from one block to the adjacent block. The variation may cause discontinuities in the gradients and discrepancies in the solution at the vicinity of block boundaries. Of course, the discrepancies become more serious when a high gradient region passes through the block boundaries [22]. From the conservation perspective, Rai [49] developed a conservative finite-difference-based zonal-boundary approach using the integration scheme to update the points on the block boundaries. Kassies and Tognaccini [33] used local grid refinement in a finite volume context to allow cell-length and skewness discontinuities of the grid normal to the block boundaries. However, there are variety of mesh refinements which do not respect the conservation laws. For example, the chimera method is a major stepstone toward a hybrid grid. While its geometrical flexibility is widely attractive, interpolation of data in the overlapped regions is done in a non-conservative fashion. There are remedies to such imperfectness [42,43,61]. Lange et al. [39] investigated the effect of local block refinement and used patched structured grids for the blockwise refinement to solve incompressible flows with heat transfer. They demonstrated the sensitivity of the solution to the refinement rate in their collocated finite volume method. From the accuracy perspective, we focus on several benchmark solutions. As is known, any newly extended formulation in CFD should be initially validated against the standard benchmark problems which provide accurate and reliable solutions for comparison. We also validate our novel extended formulations by solving standard test cases using hybrid grid topologies. For the sake of brevity, we do not demonstrate mesh study for all of them. The lid-driven square cavity is chosen as the first test case. The problem is tested at Re = 5000. Fig. 4 illustrates two types of grid distributions in the cavity. In Fig. 4(a), the solution domain is broken into two interior and exterior blocks. This hybrid structured topology has been similarly suggested and studied by Kim and Choi [35]. Alternatively, Fig. 4(b) demonstrates a scattered distribution of triangles among quadrilaterals. These two grids are identified as C1 and C2. There are 24% and 38% element reductions in C1 and C2 comparing with a domain filled fully with triangles. The cavity problem is tested using C1 and C2 grid topologies as well as single either triangular (T) or quadrilateral (Q) grids. The grid resolutions are 64 · 64 and 32 · 32 of which the latter one is shown in Fig. 4. It should be notified that the flow at Re = 5000 is a complex one with a primary circulation and several secondary circulations in the domain. Indeed, the correct detection of small vorticities in the domain needs utilizing sufficient fine grid resolutions. Fig. 5 illustrates the centerline velocity profiles in the cavity using fine and coarse grid resolutions. The current

M. Darbandi, A. Naderi / Comput. Methods Appl. Mech. Engrg. 196 (2006) 321–336 C1 Type

C2 Type 1

Y

Y

1

0.5

0

329

0

0.5

0.5

0

1

0

0.5

X

1

X

(a) Two interior–exterior blocks C1.

(b) Irregular scattered blocks C2.

Fig. 4. Two types of multiblock topologies filled with dual element considerations.

U -0.4

-0.2

0

U

0.2

0.4

0.6

0.8

1

1

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

1

0.4

0.4

0.3

Re=5000, Grid=64 × 64 0.75 C1 C2 Q T Ghia

Y 0.5

0.3

0.2

Re=5000, Grid=32×32

0.75

C1 C2 Q T Ghia

0.1 0

V

Y 0.5

0.1 0

-0.1

-0.3

-0.2 0.25

-0.3

-0.4 -0.5

0 0

0.25

0.5

0.75

1

X

(a) Fine grid.

V

-0.1

-0.2 0.25

0.2

-0.4 -0.5

0 0

0.25

0.5

0.75

1

X

(b) Coarse grid.

Fig. 5. Centerline velocity profiles in cavity using four different grid topologies with fine and coarse grid resolutions, Re = 5000.

results are compared with those of benchmark [26]. The benchmark utilizes a grid resolution of 257 · 257. Fig. 5 shows that there is no apparent distinction between the results of different grid topologies. However, Fig. 5(b) shows that there are slight differences among the solutions when very coarse grid is utilized. This is evident because a very coarse grid generally generates false diffusion in the regions with high flow gradients. Refs. [15,48] solve similar test problem using different triangular mesh orientations. They report that the proper choice of triangular mesh orientation can significantly improve the accuracy of the coarse grid solution. Ref. [14] has further investigated the regions beside the four corners of the cavity in order to ensure that the solution accuracy is not deteriorated even

in the regions with higher sensitivities and complexities. Considering theses points, it can be said that the current results show excellent accuracy despite using a relatively coarse grid distribution. This important outcome originates from the fact that the convective and diffusive flux terms are accurately estimated at the cell faces, see Eqs. (19), (20), (22) and (23). The results in Fig. 5 show that there is no discrepancy at or around the block boundaries. This ensures that the extension to hybrid grid employments has been fulfilled successfully. The next is to examine the natural convection heat transfer problem in the square cavity. In this regard, the top and bottom edges are insulated. The left and right walls are maintained at constant cold and hot temperatures,

330

M. Darbandi, A. Naderi / Comput. Methods Appl. Mech. Engrg. 196 (2006) 321–336

respectively. The problem is tested at a Rayleigh number of 105. Different grid topologies are similarly chosen to treat the buoyant flow. We only present the results of using grid type C1 here. Fig. 6 shows the streamlines and isotherms inside the cavity. The qualitative and quantitative comparisons with the benchmark solutions [12,20] indicate great accuracy of the current solutions. For the sake of brevity, the details of quantitative study performed on the Nusselt number, the maximum magnitudes in the velocity profiles, etc. are not reported here. Similar to the fluid flow part, no disturbances are observed at or around the block boundaries. Furthermore, to validate the extended formulations, the second test case is the lid-driven flow in standard skewed cavity. The problem is solved for two skew angles of 30 and 45 of which the latter one is seen in Fig. 7. Following the benchmark solution of Demirdzic et al. [21], the problem is tested at Re = 1000. This cavity is broken into suitable blocks. Two types of multiblock meshes are illustrated in Fig. 7 and labelled as SC1 and SC2. Fig. 7(a) breaks the domain into two triangular and one rectangular blocks. Following Davidson [17], Fig. 7(b) demonstrates six blocks involving four triangles and two rectangles. Comparing with Fig. 7(a), Fig. 7(b) replaces more triangles with rectangles. Therefore, the total number of elements reduces from

2400 in Fig. 7(a) to 1920 in Fig. 7(b). However, the number of quadrilaterals increases to 896. Consequently, one side advantage of multiblock employment is to reduce the number of elements within a specified solution domain, i.e., it is computationally cost effective. For example, the defined skewed cavity needs a total of 2816 triangles to cover a domain with 44 · 32 grid nodes. However, the multiblock approach reduces the required elements to 2400 in Fig. 7(a) and to 1920 in Fig. 7(b). In another words, there are 15% and 32% reductions in the number of elements using SC1 and SC2 topologies, respectively. The centerline u and v velocity profiles are plotted in Fig. 8 using a grid resolution of 88 · 64 and two skew angles of 45 and 30. The latter angle has been also investigated by Davidson [17] using different combinations of multiblock meshes. The vertical centerline is now skewed. The employed grid sizes are twice finer than the grids shown in Fig. 7 while the grid topologies are the same. The current results are compared with each other and those of the benchmark solutions [21]. The present predictions are in excellent agreement with the benchmark computations. As is observed, no discrepancies are observed around the block boundaries. In Fig. 8(b), the current results are also compared with the best results of Davidson [17]. Davidson utilizes a total of 2707 cells in his fine grid.

1

1

T Q

0.5

0

Y

Y

T Q

0

0.5

1

0.5

0

0

0.5

X

(a) Streamlines.

X

(b) Isotherms.

Fig. 6. The natural convection in convecting cavity using grid type C1, Ra = 105.

(a) Three mesh blocks; SC1.

(b) Six mesh blocks; SC2.

Fig. 7. A skewed cavity broken into two different multiblock topologies.

1

M. Darbandi, A. Naderi / Comput. Methods Appl. Mech. Engrg. 196 (2006) 321–336

U

U 0

0.25

331

0.5

0.75

1

-0.25 0.03

1

0

0.25

0.5

0.75

1

0.01

0.02 0.75

0.01

+ Y

0 0.5

0.005

0.75

V

+ Y

0

V

0.5

-0.01

Re=1000, Grid=88× 64

-0.005

Re=1000, Grid=88× 64 -0.02

SC1 SC2 T Demirdzic

0.25

-0.01

SC1 SC2 T Davidson Demirdzic

0.25 -0.03 -0.04

0

-0.015

-0.02

0 0

0.25

0.5

0.75

0.015

1

1

0

+ X

0.25

0.5

0.75

1

+ X

(a) Skew angle of 45.

(b) Skew angle of 30.

Fig. 8. The velocity profiles along the centerlines of two different skewed cavities, Re = 1000.

10

-1

10-2

Res

Similar to the current work, Davidson validates his newly pressure-correction multiblock method against a number of benchmark solutions including the skewed cavity problem. The comparison shows that the accuracy in Ref. [17] is not comparable with that of this work. To evaluate the performance of the developed algorithm in solving different multiblock structures, we study the history of pressure residuals using SC1 and SC2 topologies. The residuals are computed using the root-mean-squared of the unknown velocity components and pressure variable. Fig. 9 depicts the pressure residual history with iterations. As is observed, the two histories perform almost identical trends. Such behavior was similarly observed in studying the convergence histories of the other residuals as well. The skewed cavity is also investigated for solving the natural convection problem using the boundary conditions implemented in Fig. 6. The problem was studied for the two multiblock structures. There is no standard benchmark solution to validate the current results. Therefore, Fig. 10 only demonstrates the qualitative behavior of the flow and heat in the cavity. The streamlines and isotherms are presented at Ra = 105. Comparing with Fig. 6, there is only one stretched vortex in the skewed cavity instead of two round vorticities which appear in the squared cavity. The next test case is the steady recirculating viscous flow in an equilateral triangular cavity. The top horizontal wall is sliding with a constant velocity. A primary eddy and several secondary eddies at different regions indicate the complexity of the flow field. Following Refs. [27,53], the problem is tested at two Reynolds of 100 and 500. Fig. 1 shows the solution domain which is discretized considering a hybrid distribution of elements in four different blocks. The grid size is 32 on lid and 22 on side walls in this figure. However, the problem is solved for a grid size of 128 and

10

-3

10

-4

SC1 SC2

10-5 10

-6

10-7 10

-8

10

-9

50

100

150

200

Iteration Fig. 9. The history of pressure residual in solving the flow inside the skewed cavity.

88 in our computations. This grid topology provides 25% reduction in the number of elements with respect to a domain fully covered with the triangular elements. Fig. 11 shows the streamlines and pressure contours. As is observed, the recirculation on the left wall and several recirculations on the centerline emphasize the complexity of flow at higher Reynolds. The depicted streamlines are in excellent agreement with those of Refs. [27,53]. Unfortunately, the references do not provide the details of solution inside the cavity. Ref. [16] solves the problem using a multiblock strategy with combined structured and unstructured distribution of triangles. Comparing the results of the current work with those of Ref. [16], it indicates that

332

M. Darbandi, A. Naderi / Comput. Methods Appl. Mech. Engrg. 196 (2006) 321–336

the multiblock decomposition incorporated with a hybrid grid has absolute superiority with respect to a pure triangular structured–unstructured grid. The study has proven that the current hybrid strategy needs much less grid sizes

to achieve the same accuracy which can be reached by using a structured–unstructured grid topology. It should be notified that this point is numerously important in solving the fluid flow problems utilizing implicit algorithms.

(b) Isotherms.

(a) Streamlines.

Fig. 10. The natural convection in a convecting skewed cavity using grid type SC2, Ra = 105.

0.9

Re=100, Grid=128, Stream

0.9 0.8

0.7

0.7

0.6

0.6

0.5

0.5

Y

Y

0.8

Re=500, Grid=128, Stream

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0 -0.5

-0.25

0

0.25

0 -0.5

0.5

-0.25

0

X (a) Streamlines at Re = 100.

0.9

0.5

(b) Streamlines at Re = 500.

0.9

Re=100, Grid=128, P Contour

0.0

0.8

-0.0075

Re=500, Grid=128 , P Contour 03 0.

0.8

0.25

X

-0

5 07 3 .0 -0 -0.0

.03

3

3 0.0 75 0 -0.0

0.7

0.6

0.6

Y

-0.0

075

-0.12

0.5

Y

0.5

-0.03

-0.03

0.4

-0.0075

5

0.4

0.0

0

07

0.7

0.3

0.3

0.2

0.2

0.1

0.1

0

0

0

0 -0.5

-0.25

0

X (c) Isobars at Re = 100.

0.25

0.5

0 -0.5

-0.25

0

0.25

X (d) Isobars at Re = 500.

Fig. 11. The streamlines and isobars in the triangular cavity using two Reynolds of 100 and 500 and a grid resolution of 128 · 88.

0.5

M. Darbandi, A. Naderi / Comput. Methods Appl. Mech. Engrg. 196 (2006) 321–336

The last test case is a cylinder placed in an unbounded uniform flow. If the diameter of the cylinder d and the free stream velocity U1 are considered as the length and velocity scales, respectively, the Reynolds number is given by Re = U1d/m. In order to achieve a steady solution, the computation is carried out at Re = 20 and 40. We use two types of convex and semi-convex grids to discretize the solution domain. The past investigation has shown that the use of a quadrilateral structured grid around the immersed bodies such as airfoils and cylinders can effectively improve the numerical solution close to those bodies, e.g., see Refs. [34–36,43]. In another words, the uniform structured grid provides better efficiency in treating the boundary layer adjacent to the immersed bodies. In the current work, the simulations were performed on a domain with 20d · 30d dimensions, see Fig. 12. Two types of hybrid grid topologies G1 and G2 were utilized in this work. Fig. 12(a) demonstrates a two-block hybrid grid consisted of structured and unstructured grids. This grid was employed using two different fine and coarse mesh resolutions. Fig. 12(a) presents the coarse grid with 17 · 17 grid points in the streamwise and transverse directions. However, to ensure the mesh independent solution, the computations have been also carried out using a fine grid of 33 · 33 resolutions. Fig. 12(b) demonstrates a semi-convex grid topology G2 consisted of three structured and unstructured blocks in the domain. We have tested this topology using either fine or coarse grid resolutions. Contrary to G1 topology, the number of nodes at the four boundaries is not necessarily the same. As is known, it is peculiar to many finite volume approaches based upon collocated meshes to have poor performance when the mesh topology is such that the shaded area in Fig. 2 is non-convex. In this instance, the pressure reconstruction is particularly problematic even if the scheme is fully conservative. To check the performance of current formulation on non-convex grids, we have discretized the solution domain using very stretched, non-uniform meshes. It would be interesting to evaluate the results

on such mesh topologies for which many finite-volume techniques do not perform well. There are 4844 triangular and 160 quadrilateral elements in G1 fine grid resolution. The number of triangles reduces to 2990 in semi-convex grid whereas the number of quadrilaterals does not change. This indicates that the semi-convex grid with fine resolution is coarser than the convex grid with fine resolution. Subsequently, this may cause deterioration in the accuracy of the results of the latter grid in comparison with that of the former grid. Fig. 13 shows the streamlines around the cylinder using either G1 or G2 grid topologies. The recirculation regions behind the cylinder are clearly observed in all test cases. We have concentrated on the bubbles generated behind the cylinder in Fig. 13(b) and (c). However, we have illustrated the streamlines up to the outlet section in Fig. 13(a). It is because we want to show that there are no disturbances at the downstream of the cylinder because of passing through either block boundaries or non-convex volumes. The figure shows that the streamlines are quite smooth. Additionally, comparing with the plots depicted by other computational workers, the present results compare qualitatively excellent. Fig. 14 demonstrates the residuals of pressure variable in terms of iterations at Re = 40. The residuals are plotted using both types of G1 and G2 grid topologies introduced in Fig. 12. As is observed, the non-convex grid does not degrade the convergence rate. This outcome can be considered as a great benefit of employing the current robust formulation. However, to document the rate of convergence at least from numerical observation, Fig. 14 also presents the number of iterations required to achieve a residual criterion of 109 using both G1 and G2 topologies with either fine or coarse grid resolutions. The figure shows that the convergence rates are very fast if coarse grid resolutions are utilized. Comparing the rate of convergence for the fine grids with those of coarse grids, it indicates that the order of convergence is higher than a first-order method. This conclusion is consistent with the discussion provided at

5

5

Y

10

Y

10

0

0

-5

-5

-10 -10

0

X

10

(a) Regular convex grid; G1.

20

333

-10 -10

0

X

10

(b) Irregular semi-convex grid; G2.

Fig. 12. Two types of hybrid grids around a cylinder in a channel using multiblock structured and unstructured topologies.

20

334

M. Darbandi, A. Naderi / Comput. Methods Appl. Mech. Engrg. 196 (2006) 321–336

(a) Using semi-convex grid G2, see Fig. 12(b); Re = 40.

(b) Using regular convex grid G1; Re = 20.

(c) Using regular convex grid G1; Re = 40.

Fig. 13. The streamlines in the domain using two types of hybrid grid topology with and without convex volumes.

Pressure Residual

10

10

-1

10

-3

10

-5

10

-7

10

-9

10

the end of Section 4 about the order of magnitude in the current work. To quantify the current results, they are compared with those of Refs. [25,57,60] in Table 1. The comparisons are performed at Re = 20 and 40 using two types of G1 and G2 grid topologies. The current results are in excellent agreement with those of others. It should be notified that the maximum differences in the values of Xc for both Reynolds are confined to less than 3% despite using a coarse grid distribution. The slight differences are mainly due to the use of coarse grid distribution in the domain and the use of linear interpolation to locate the bubbles’ reattachment point.

1

Grid 1 Grid 2 converged at G1 coarse at 134th itr. G1 fine at 163th itr. G2 coarse at 130th itr. G2 fine at 176th itr.

-11

0

50

100

150

200

250

300

350

400

450

Iterations Fig. 14. A comparative study of pressure residual histories using G1 and G2 topologies and the number of iterations is required to achieve a residual criterion of 109 using different fine and coarse grid choices, Re = 40.

Table 1 The quantification of the bubble’s length calculated in the current work with those reported in other references

Re = 20 Re = 40

Grid 1

Grid 2

Fornberg [25]

Tseng and Ferziger [57]

Ye et al. [60]

0.94 2.21

0.89 2.21

0.91 2.24

– 2.21

0.92 2.27

M. Darbandi, A. Naderi / Comput. Methods Appl. Mech. Engrg. 196 (2006) 321–336

6. Concluding remarks An efficient finite-volume-based finite-element algorithm was developed for solving fluid flow and heat transfer problems on multiblock hybrid mesh topologies. The new algorithm uses a combined triangular/quadrilateral grid, which offers great flexibility in mesh generation. This extension permits converting triangulations to quadrangulations and vice versa for either the whole domain or a part of the domain as needed. Thus, the approach is very amenable to quickly create quality grids for complex geometries with various individual components. Additionally, the approach fully preserves the advantageous features of both the triangular and quadrilateral grids while eliminating their deficiencies in grid orientation problem and the improper applications to arbitrary irregular geometries. A physical influence scheme is utilized to predict the conservative statements at the cell faces using dual element shape features. The flow solutions indicate the satisfaction of conservation properties through the block boundaries. Despite using coarse grid distributions, the current results are in good agreement with the benchmark solutions and measured data, which is an indication of the reliability of the method. No deteriorations and disturbances are observed in the solutions even on cells which are located at the block boundaries. Indeed, the solutions through joint boundaries are fully smooth and free of noises emanating from using different interpolators. Since the current method is fully implicit, the new multiblock hybrid grid can drastically reduce the computer memory which is normally required for the unstructured grid solvers and make more efficient implicit solvers accessible. Acknowledgement The authors would like to acknowledge the grant received from the Research Center of Sharif University under contraction number SUT-RC-1382-4105. References [1] M. Ainsworth, J. Coyle, Hierarchic hp-edge element families for Maxwell’s equations on hybrid quadrilateral/triangular meshes, Comput. Methods Appl. Mech. Engrg. 190 (2001) 6709–6733. [2] M.M. Alishahi, M. Darbandi, Multiple-zone potential solution around wing-body configurations, J. Aerospace Engrg. 6 (1993) 329–346. [3] J.W. Barrett, K.W. Morton, Approximate symmetrization and Petrov–Galerkin methods for diffusion–convection problems, Comput. Methods Appl. Mech. Engrg. 45 (1984) 97–122. [4] A. Beskok, T.C. Warburton, An unstructured hp finite-element scheme for fluid flow and heat transfer in moving domains, J. Thermophys. 174 (2001) 492–509. [5] A. Brooks, T.J.R. Hughes, Streamline upwind Petrov–Galerkin formulations of convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 32 (1982) 199–259. [6] M.A. Celia, I. Herrera, E. Bouloutas, J.S. Kindred, A new numerical approach for the advective–diffusive transport equation, Numer. Methods Partial Differential Equations 5 (1989) 203–226.

335

[7] M.A. Celia, T.F. Russell, I. Herrera, R.E. Ewing, An Eulerian– Lagrangian localized adjoint method for the advection–diffusion equation, Adv. Water Resour. 13 (1990) 187–206. [8] M. Darbandi, S.M. Bostandoost, A new formulation toward unifying the velocity role in collocated variable arrangement, Numer. Heat Transfer B 46 (2004) 361–382. [9] M. Darbandi, K. Mazaheri-Body, S. Vakilipour, A pressureweighted upwind scheme in unstructured finite-element grids, in: M. Feistauer, V. Dolejsi, P. Knobloch, K. Najzar (Eds.), Numerical Mathematics and Advanced Applications, Springer-Verlag, 2004, pp. 250–259. [10] M. Darbandi, V. Mokarizadeh, A modified pressure-based algorithm to solve the flow fields with shock and expansion waves, Numer. Heat Transfer B 46 (2004) 497–504. [11] M. Darbandi, G.E. Schneider, Momentum variable procedure for solving compressible and incompressible flows, AIAA J. 35 (1997) 1801–1805. [12] M. Darbandi, G.E. Schneider, Thermobuoyancy treatment for electronic packaging using an improved advection scheme, ASME J. Elect. Packag. 125 (2003) 244–250. [13] M. Darbandi, G.E. Schneider, K. Javadi, N. Solhpour, The performance of a physical influence scheme in structured triangular grids, AIAA Paper 2003-0436, The 41st AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, January 2003. [14] M. Darbandi, G.E. Schneider, A.R Naderi, A finite element volume method to simulate flow on mixed element shapes, AIAA Paper 20033638, The 36th AIAA Thermophysics Conference, Orlando, FL, June 2003. [15] M. Darbandi, G.E. Schneider, A.R. Naderi, The mesh orientation impact in performance of physical-based upwinding in structured triangular grids, in: Proc. the 11th Annual CFDSC Conference 2003, CFDSC, Ottawa, Canada, 2003, pp. 688–695. [16] M. Darbandi, G.E. Schneider, S. Vakilipour, A modified upwindbiased strategy to calculate flow on structured-unstructured grid topologies, AIAA Paper 2004-0435, The 42nd AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, January 2004. [17] L. Davidson, A pressure correction method for unstructured meshes with arbitrary control volumes, Int. J. Numer. Meth. Fluids 22 (1996) 265–281. [18] C. de Nicola, G. Pinto, R. Tognaccini, On the numerical stability of block structured algorithms with applications to 1-D advection– diffusion problems, Comput. Fluids 24 (1995) 41–54. [19] C. de Nicola, G. Pinto, R. Tognaccini, Stability of two-dimensional model problems for multiblock structured fluid-dynamics calculations, Comput. Fluids 26 (1997) 43–58. [20] G. De Vahl Davis, Natural convection of air in a square cavity: a bench mark numerical solution, Int. J. Numer. Meth. Fluids 3 (1983) 249–264. [21] I. Demirdzic, Z. Lilek, M. Peric, Fluid flow and heat transfer test problems for non-orthogonal grids: bench-mark solutions, Int. J. Numer. Meth. Fluids 15 (1992) 329–354. [22] R.E. Ewing, Z. Li, T. Lin, Y. Lin, The immersed finite volume element methods for the elliptic interface problems, Math. Comput. Simul. 50 (1999) 63–76. [23] M. Feistauer, J. Felcman, M. Lukacova-Medvid’ova, Combined finite element–finite volume solution of compressible flow, J. Comput. Appl. Math. 63 (1995) 179–199. [24] J. Fezoui, B. Stoufflet, A class of implicit schemes for Euler simulations with unstructured meshes, J. Comput. Phys. 84 (1989) 174–206. [25] B. Fornberg, A numerical study of steady viscous flow past a circular cylinder, J. Fluid Mech. 98 (1980) 819–855. [26] 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. [27] J.L. Guermond, L. Quartapelle, Calculation of incompressible viscous flows by an unconditionally stable projection FEM, J. Comput. Phys. 132 (1997) 12–33.

336

M. Darbandi, A. Naderi / Comput. Methods Appl. Mech. Engrg. 196 (2006) 321–336

[28] P. Hansbo, A. Szepessy, A velocity–pressure streamline diffusion finite element method for the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 84 (1990) 107–129. [29] T. Hayase, J.A.C. Humphrey, R. Greif, A consistently formulated QUICK scheme for fast and stable convergence using finite-volume iterative calculation procedures, J. Comput. Phys. 98 (1992) 108–118. [30] R.W. Healy, T.F. Russel, A finite-volume Eulerian–Lagrangian localized adjoint method for solution of the advection–dispersion equation, Water Resour. Res. 29 (1993) 2399–2413. [31] C.J. Hwang, S.J. Wu, Adaptive finite volume upwind approach on mixed quadrilateral–triangular meshes, AIAA J. 31 (1993) 61–67. [32] C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method, Cambridge University Press, Cambridge, 1987. [33] A. Kassies, R. Tognaccini, Boundary conditions for Euler equations at internal block faces of multi-block domains using local grid refinement, AIAA paper no. 1990-1590, 1990. [34] A. Khawaja, T. Minyard, Y. Kallinderis, Adaptive hybrid grid methods, Comput. Methods Appl. Mech. Engrg. 189 (2000) 1231– 1245. [35] D. Kim, H. Choi, A second-order time-accurate finite volume method for unsteady incompressible flow on hybrid unstructured grids, J. Comput. Phys. 162 (2000) 411–428. [36] R.M. Kirby, T.C. Warburton, I. Lomtev, G.E. Karniadakis, A discontinuous Galerkin spectral/hp method on hybrid grids, Appl. Numer. Math. 33 (2000) 393–405. [37] K.B. Kuan, C.A. Lin, Adaptive QUICK-based scheme to approximate convective transport, AIAA J. 38 (2000) 2233–2237. [38] L. Lammer, M. Burghardt, Parallel generation of triangular and quadrilateral meshes, Adv. Engrg. Softwares 31 (2000) 929–936. [39] C.F. Lange, M. Schafer, F. Durst, Local block refinement with a multigrid flow solver, Int. J. Numer. Methods Fluids 38 (2002) 21–41. [40] B.P. Leonard, A stable and accurate convective modelling procedure based on quadratic upstream interpolation, Comput. Methods Appl. Mech. Engrg. 19 (1979) 59–98. [41] X. Li, W. Wu, O.C. Zienkiewicz, Implicit characteristic Galerkin method for convection–diffusion equations, Int. J. Numer. Methods Engrg. 47 (2000) 1689–1708. [42] M.S. Liou, K.H. Kao, Progress in grid generation: from Chimera to DRAGON grids, NASA TM 106709, NASA, August 1994. [43] M.S. Liuo, Y. Zheng, A novel approach of three-dimensional hybrid grid methodology: Part 2. Flow solution, Comput. Methods Appl. Mech. Engrg. 192 (2003) 4173–4193. [44] P. Moinier, J.D. Muller, M.B. Giles, Edge-based multigrid and preconditioning for hybrid grids, AIAA J. 40 (2002) 1954–1960.

[45] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, 2nd ed., Hemisphere Publishing Co., New York, 1996. [46] O. Pironneau, Finite Element Methods for Fluids, John Wiley and Sons Ltd., Chichester, 1989. [47] A. Pollard, A. Siu, The calculation of some laminar flows using various discretization schemes, Comput. Methods Appl. Mech. Engrg. 35 (1982) 293–313. [48] C. Prakash, Examination of the upwind (Donor-Cell) formulation in control volume finite-element methods for fluid flow and heat transfer, Numer. Heat Transfer 11 (1987) 401–416. [49] M.M. Rai, An implicit, conservative, zonal-boundary scheme for Euler equation calculations, Comput. Fluids 14 (1986) 295–319. [50] G.D. Raithby, K.E. Torrance, Upstream-weighted differencing schemes and their application to elliptic problems involving fluid flow, Comput. Fluids 2 (1974) 191–206. [51] G.D. Raithby, Skew upstream differencing schemes for problems involving fluid flow, Comput. Methods Appl. Mech. Engrg. 9 (1976) 153–164. [52] C.M. Rhie, W.L. Chow, Numerical study of the turbulent flow past an airfoil with trailing edge separation, AIAA J. 21 (1983) 1525–1532. [53] C.J. Ribbens, L.T. Watson, C.Y. Wang, Steady viscous flow in a triangular cavity, J. Comput. Phys. 112 (1994) 173–181. [54] T.F. Russell, M.A. Celia, An overview of research on Eulerian– Lagrangian localized adjoint methods (ELLAM), Adv. Water Resour. 25 (2002) 1215–1231. [55] R. Struijs, P. Vankeirsbilck, H. Deconinck, An adaptive grid polygonal finite volume method for the compressible flow equations, AIAA Paper 1989–1957, 1989. [56] A.J. Toreja, Rizwan-uddin, Hybrid numerical methods for convection–diffusion problems in arbitrary geometries, Comput. Fluids 32 (2003) 835–872. [57] Y.H. Tseng, J.H. Ferziger, A ghost-cell immersed boundary method for flow in complex geometry, J. Comput. Phys. 192 (2003) 593–623. [58] M.F. Wheeler, C.N. Dawson, An Operator-Splitting Method for Advection–Diffusion-Reaction Problems, The Mathematics of Finite Elements and Applications (Uxbridge, 1987), vol. VI, Academic Press, London–New York, 1988, pp. 463–482. [59] W.A. Wood, W.L. Kleb, Diffusion characteristic of finite volume and fluctuation splitting schemes, J. Comput. Phys. 153 (1999) 353–377. [60] T. Ye, R. Mittal, H.S. Udaykumar, W. Shyy, An accurate Cartesian grid method for viscous incompressible flows with complex immersed boundaries, J. Comput. Phys. 156 (1999) 209–240. [61] Y. Zheng, M.S. Liuo, A novel approach of three-dimensional hybrid grid methodology: Part 1. Grid generation, Comput. Methods Appl. Mech. Engrg. 192 (2003) 4147–4171.