Immersed Boundary Method Based on the Implementation of Conservation Equations along the Boundary using Control-Volume Finite-Element Scheme

ABSTRACT

In this study conservation equations were implemented along the boundaries via ghost control-volume immersed boundary method. The control-volume finite-element method was applied on a cartesian grid to simulate 2-D incompressible flow. In this approach, mass and momentum equations were conserved in the whole domain including boundary control volumes by introducing ghost-control volume concept. The Taylor problem was selected to validate the present method. Four different case studies of Taylor problem encompassing both inviscid and viscous flow conditions in ordinary and 45º rotated grid were used for more investigation. Comparisons were made between the results of the present method and those obtained from the exact solution. Results of the present method indicated accurate predictions of the velocity and pressure fields in midline, diagonal, and all boundaries. The agreement between the results of the present method and the exact solution was very good throughout the whole temporal domain. Furthermore, comparison of the rate of kinetic energy decay in viscous case showed same level of agreement between the results.

KEYWORDS:
Immersed boundary method; Control-volume-based finite element; Sub-control volumes; Conservation of mass and momentum equations; Ghost node; Ghost sub-control volume

INTRODUCTION

The immersed boundary method (IBM) is known as a powerful approach for simulating flows in moving boundary and complex geometry problems. In this method, discretization of equations is carried out on a Cartesian grid, which is simple to generate. However, the boundary does not conform to the grid lines, and therefore indirect methods are employed to apply the boundary conditions. This creates a range of different methods developed in the context of IBM which are applied to elastic (Peskin 1972, 1982; Beyer Jr 1992; Fauci and McDonald 1995; Zhu and Peskin 2003) and solid (Berger and Aftosmis 1998; Khadra et al. 2000; Tseng and Ferziger 2003; Saiki and Biringen 1996) boundaries. The conventional ghost-node method is currently used in problems with solid boundaries, where the value of ghost-node is set as to meet the boundary conditions. In ghost-node methods, finite difference scheme is usually used to simulate the flow field and the value of ghost-node is determined using a kind of interpolation schemes (Mittal and Iaccarino 2005; Majumdar et al. 2001; Ghias et al. 2004; Mittal et al. 2008). While these approaches are considered fairly fast in convergence and simple in application, mass and momentum equations are not conserved in applying boundary conditions. However, the so-called cut-cell method is a complicated approach based on Cartesian grid (Clarke et al. 1986; Udaykumar et al. 2001, 1999, 1996; Ye et al. 1999), which implements conservation laws in boundary cells. In this method the shape of Cartesian cells in the vicinity of the boundary is changed to fit the boundary. In cut-cell method, cells are divided by the boundary, and conservation laws are implemented in divided cells conforming to the boundary. Comparing to ghost node methods used in IBM, cut-cell method is extremely complicated. This is because the boundary may cut the Cartesian grids anywhere on the cells and create new arbitrary shape. It would make it more difficult to discretize the equations and calculation of fluxes particularly in 2- and 3-D and moving problems.

In the present study, an immersed boundary method based on CVFE scheme is proposed in the context of ghost node concept in which conservation of conserved quantities is enforced. Importantly, the present method has the capability to conserve mass and momentum equations along the boundary. The present approach is different from the cut-cell method such that boundary cell shapes remain unchanged.

NUMERICAL ALGORITHM

The governing equations in the present method are solved via CVFE scheme, which was presented by Minkowycz et al. (1988) to discretize governing equations. Sub-control-volume (SCV) and node types are further explained to implement boundary conditions.

CONTROL-VOLUME FINITE-ELEMENT METHOD

In this scheme, solution domain is always discretized into a number of Cartesian elements. As shown in Fig. 1a, a local coordinate system (s,t) is defined in the middle of each element. This local coordinate system divides each element into 4 SCVs. Each SCV is associated with an element node at its vertex. Therefore, as shown in Fig. 1c, the grey area represents a control volume made from surrounding 4 SCVs neighbour elements. All primitive variables are located at the vertices of the elements, placing them in the middle of each control volumes.

Figure. 1
(a) Definitions of the element. Local coordinate system of (s,t) is located in the middle of the element, sub-control surface is indicated, and integral points are shown via cross symbols in the middle of sub-control surfaces; (b) The grey area is the SCV, and surface normal vectors are indicated in its outward direction; (c) The dark grey area in the center of the figure is control volume made up of 4 surrounding SCVs and the light grey area is SCV.

Although governing equations are finally conserved on control volumes such as the one shown in Fig. 1b, their formation are done through the assembly of elemental equations (Minkowycz et al. 1988). Elemental equations of each element include conservation of governing equations on the 4 SCVs of that element. Variables and their gradients should be evaluated at the integration points (Fig. 1a) to determine the flux at each sub-control surface. Variables with elliptic nature or of diffusion type such as pressure and diffusion can be calculated using bilinear interpolation. Minkowycz et al. (1988) presented a bilinear shape function to determine the value of variables everywhere in the element (Fig. 1a). Accordingly the value of variable φ and its gradients can be determined by:

(1) $\mathit{\phi }\left(\mathit{s},\mathit{t}\right)=\sum _{\mathit{i}=1}^4{\mathit{N}}_{\mathit{i}}\left(\mathit{s},\mathit{t}\right){\Phi }_{\mathit{i}}$

(2) $\frac{\partial \mathit{\phi }\left(\mathit{s},\mathit{t}\right)}{\partial \mathit{x}}=\sum _{\mathit{i}=}^4\frac{\partial {\mathit{N}}_{\mathit{i}}\left(\mathit{s},\mathit{t}\right)}{\partial \mathit{x}}{\Phi }_{\mathit{i}}=\sum _{\mathit{i}=1}^4{\mathit{D}}_{\mathit{x},\mathit{i}}{\Phi }_{\mathit{i}}$

where ϕ_i is the value of φ at the vertices of each element; N_i is the i^th bilinear shape function.

Modelling of other variables without elliptic nature or diffusion type such as velocity components in mass fluxes and convection terms will be discussed in more details later. Details of the CVFE method and the formation of the system of governing equation were presented by Minkowycz et al. (1988).

SUB-CONTROL VOLUMES AND NODE TYPES

Discretization of governing equations and calculation of fluxes are done on SCVs; hence, the classification of different sub-control volumes and nodes is described here. There are 4 SCVs in each element as previously explained according to Fig. 1a. Depending on the location of elements in the domain, the SCVs and nodes are classified into 3 types in this paper. The first type of SCV is the "ordinary" or "fluid" one that is in the middle of the solution domain and it has no boundary in its SCV or in its related element (Fig. 2). An ordinary node is assigned to each related ordinary SCV. In the second type the boundary has crossed the SCV. This type of SCV and its pertaining node are called ghost SCV type I and ghost node type I, respectively (Fig. 3). Lastly, as shown in Fig. 4, the third type is defined when the boundary is placed in the SCVs of fluid nodes in the element. These SCVs are called ghost SCV type II and accordingly each related node is called ghost node type II (Fig. 4). To conclude, in this method, whenever the immersed boundary is placed within an element, nodes outside of the flow field are called ghost nodes (nodes 2 and 3 in Figs. 3 and 4) and their corresponding SCVs (SCVs 2 and 3 in Figs. 3 and 4) are called ghost SCVs. In this paper boundary conditions are applied via ghost SCVs (Fig. 5). Note that SCVs of both fluid and ghost nodes are always considered as ordinary SCVs or ghost SCVs, respectively, regardless of the boundary location. As noted earlier in conventional IBMs (sharp interface methods - Seo and Mittal 2011; Ghias et al. 2007) boundary conditions are applied via the assignment of appropriate values for the flow variables to the ghost nodes. These values are mostly assigned by a kind of interpolation scheme (Mittal and Iaccarino 2005; Majumdar et al. 2001; Ghias et al. 2004, 2007). In the present method, however, flow variables on the ghost nodes are determined by implementation of conservation laws and the boundary condition on ghost SCVs. Details of the method will be discussed in following section.

Figure 2
Ordinary SCV and node. all 4 SCVs and nodes are ordinary (fl uid).

Figure 3
Ghost SCV and node type I; grey area indicated in SCVs 2 and 3 are ghost SCVs type I, nodes 1 and 4 are fl uid nodes and nodes 2 and 3 are ghost node type I.

Figure 4
Ghost SCV and node type II. The grey area indicated in SCVs 1 and 4 are ghost SCV type II related to nodes 2 and 3 (ghost nodes type II). Nodes 1 and 4 are fl uid ones.

Figure 5
SCV and node types classification.

GOVERNING EQUATIONS AND DISCRETIZATION

In Eq. 3 there is a detail analysis of how Navier-Stokes equations were discretized. The integral form of the incompressible Navier-Stokes equations for 2-D flow is given by

(3) $\begin{array}{c}\hfill {\int }_{\mathit{v}}\frac{\partial \mathit{Q}}{\partial \mathit{t}}\mathit{dV}+{\int }_{\mathit{s}}\left({\mathit{Eds}}_{\mathit{x}}+{\mathit{Fds}}_{\mathit{y}}\right)=\\ \hfill ={\int }_{\mathit{s}}\left({\mathit{Gds}}_{\mathit{x}}+{\mathit{Hds}}_{\mathit{y}}\right)\end{array}$

where Q is the vector of conserved quantities; E and F are convection flux vectors; G and H are diffusion flux vectors; ν is volume.

The extended form of these vectors is:

(4) $\begin{array}{c}\mathit{Q}=\left(\begin{array}{c}\mathit{o}\\ \mathit{\rho }\mathit{U}\\ \mathit{\rho }\mathit{V}\end{array}\right);\mathit{E}=\left(\begin{array}{c}\left(\mathit{\rho }\mathit{u}\right)\\ \left(\mathit{\rho }\mathit{u}\right)\mathit{u}+\mathit{p}\\ \left(\mathit{\rho }\mathit{u}\right)\mathit{v}\end{array}\right);\hfill \\ \mathit{F}=\left(\begin{array}{c}\left(\mathit{pv}\right)\\ \left(\mathit{pv}\right)\mathit{u}\\ \left(\mathit{pv}\right)\mathit{v}+\mathit{p}\end{array}\right);\mathit{G}=\left(\begin{array}{c}0\\ {\mathit{\tau }}_{\mathit{xx}}\\ {\mathit{\tau }}_{\mathit{xy}}\end{array}\right)\mathit{\tau };\mathit{H}=\left(\begin{array}{c}0\\ {\mathit{\tau }}_{\mathit{yx}}\\ {\mathit{\tau }}_{\mathit{yy}}\end{array}\right)\hfill \\ {\mathit{\tau }}_{\mathit{xx}}=2\mathit{\mu }\frac{\partial \mathit{u}}{\partial \mathit{x}};{\mathit{\tau }}_{\mathit{yy}}=2\mathit{\mu }\frac{\partial \mathit{v}}{\partial \mathit{y}};\hfill \\ {\mathit{\tau }}_{\mathit{xy}}={\mathit{\tau }}_{\mathit{yx}}=\mathit{\mu }\left(\frac{\partial \mathit{u}}{\partial \mathit{y}}+\frac{\partial \mathit{v}}{\partial \mathit{x}}\right)\hfill \end{array}$

where ρ represents density; u and v are velocity in x and y directions, respectively; τ is shear stress; µ is viscosity; p means pressure.

Upper-case letters were used to indicate nodal values and the lower-case ones, to show the values of variables on integral points (ip). After substituting stress tensor within G and H, the simplified form would be as shown in Eq. (5).

(5) $\mathit{G}=\left(\begin{array}{c}0\\ \mathit{\mu }\frac{\partial \mathit{u}}{\partial \mathit{x}}\\ \mathit{\mu }\frac{\partial \mathit{v}}{\partial \mathit{x}}\end{array}\right);\mathit{H}=\left(\begin{array}{c}0\\ \mathit{\mu }\frac{\partial \mathit{u}}{\partial \mathit{y}}\\ \mathit{\mu }\frac{\partial \mathit{v}}{\partial \mathit{y}}\end{array}\right)$

Firstly the ordinary SCV is explained. Navier-Stokes equations should be discretized in all of the four SCVs of each element in order to form element-level equations. In a case of ordinary SCV, the process of discretizing is straightforward as described in Karimian and Schneider (1994a). This process is explained in more details as follows. SCV 1 in Fig. 1a is considered here where Eq. 3 is written for this SCV as follows:

(6) $\begin{array}{c}\hfill {\int }_{\mathit{SCV}1}\frac{\partial \mathit{Q}}{\partial \mathit{t}}\mathit{dV}+{\int }_{\mathit{SS}1\&\mathit{SS}4}\left({\mathit{EDS}}_{\mathit{x}}+{\mathit{Fds}}_{\mathit{y}}\right)=\\ \hfill {\int }_{\mathit{SS}1\&\mathit{SS}4}\left({\mathit{Gds}}_{\mathit{x}}+{\mathit{Hds}}_{\mathit{y}}\right)\end{array}$

where SS stands for the inner sub-control surface shown in Fig. 1a; ds_x and ds_y are the components of normal surface vector in the outward direction.

The volume integral of the transient term is estimated using a lumped approach. Surface-integrals of E, F, G, and H are calculated by their average values over SS at the midpoint location, which is the ip. For the diffusion flux vectors G and H, bi-linear interpolation (Eq. 2) is used to directly evaluate the components of stress tensor (Karimian and Schneider 1994b). In the convection flux vectors E and F, pressure is evaluated using bilinear interpolation (Eq. 1), and the momentum fluxes are linearized with respect to mass fluxes and. Velocity components u and v in mass fluxes are called integral-point convecting velocities and have been previously denoted by (ρu) and (ρv) (Karimian and Schneider 1994a). Other values of u and v in the momentum fluxes, which are convected by the mass fluxes through the control-volume surfaces, are called convected velocities. Convecting and convected velocities are cell-face, which are modelled in terms of nodal values of velocity and pressure.

Karimian and Schneider (1994a) reported the implementa- tion of the corresponding governing equations of flow to derive cell-face velocities (convected and convecting velocities) (Karimian and Schneider 1994a). In this method convected velocity is obtained from the following equation:

(7) $\mathit{\rho }\frac{\partial \mathit{u}}{\partial \mathit{t}}+\mathit{\rho }\mathit{q}\frac{\partial \mathit{u}}{\partial \mathit{s}}+\frac{\partial \mathit{p}}{\partial \mathit{x}}-\mathit{\mu }{\nabla }^2\mathit{u}=0$

In Eq. 7 the convection term is represented in stream wise direction and q = (u² + v²)^1/2. Expression for convected velocity u is obtained on integration points which encompass all relevant variables related to flow condition. The convecting velocity û on ip is obtained from Eq. 8 as follows:

(8) $\begin{array}{c}\left(\mathit{\rho }\frac{\partial \mathit{u}}{\partial \mathit{t}}+\mathit{pq}\frac{\partial \mathit{u}}{\partial \mathit{s}}+\frac{\partial \mathit{p}}{\partial \mathit{x}}-\mathit{\mu }{\nabla }^2\mathit{u}\right)-\hfill \\ -\mathit{u}\left(\mathit{\rho }\frac{\partial \mathit{u}}{\partial \mathit{x}}+\mathit{\rho }\frac{\partial \mathit{v}}{\partial \mathit{y}}\right)=0\hfill \end{array}$

For details about the modeling of cell-face velocities and their role in resolving pressure velocity decoupling in incompressible flow, see (Karimian and Schneider 1994a).

In the current research after completing the discretization of Navier-Stokes equations, a fully coupled algorithm is used to solve the resulted system of equations to obtain the flow variables (pressure and velocity components: p, u, and v). This system of equations is solved simultaneously using a band solver.

BOUNDARY CONDITIONS AND GHOST SUB-CONTROL VOLUMES

In IBM, flow variables are assigned so that their value guarantees satisfaction of boundary condition on the immersed boundary. As mentioned before, in the present method flow, variables on the ghost nodes are determined by implementation of conservation laws and the boundary condition on ghost SCVs. Therefore, the key-point in the present method is to clearly implement conservation laws on ghost SCVs along the boundaries. This process is explained here for the ghost SCVs types one and two.

Ghost Sub-Control Volume Type I

In Fig. 6 an element with ordinary SCVs 1 and 4, and ghost SCVs 2 and 3 is presented. Implementation of Eq. 3 on ordinary SCVs 1 and 4 is done as described in previous section. Thus, mass and momentum conservation equations on ordinary SCV 1 would be:

Figure 6
Ghost SCV type I (grey area in SCV 2 is considered); SSl is the left part of sub-surface 2; SS2r is the right part of sub-surface 2 along the grey area; SSb is the boundary portion in SCV 2; v₂ is the volume of the grey area; ds_b is normal surface vector of boundary in SCV 2 in direction to outward of the grey area; ds_x²¹ is the normal surface vector of sub-surface1; ds_y²_2r is normal surface vector related to right part of sub-surface 2; ∆x and ∆y are grid dimensions; points 1, 2, 3, and 4 indicated with cross symbols are ip.

(9) $\mathit{\rho }{\hat{\mathit{u}}}_1{\mathit{ds}}_{{\mathit{x}}_1^1}+\mathit{p}{\hat{\mathit{v}}}_4{\mathit{ds}}_{{\mathit{y}}_4^1}=0$

(10) $\begin{array}{c}\hfill \frac{{\mathit{\rho }\mathit{\nu }}_1\left({\mathit{U}}_1-{\mathit{U}}_1^0\right)}{\Delta \mathit{t}}+\left(\mathit{p}{\hat{\mathit{u}}}_1{\mathit{ds}}_{{\mathit{x}}_1^1}\right){\mathit{u}}_1+\left(\mathit{\rho }{\hat{\mathit{v}}}_4{\mathit{ds}}_{{\mathit{y}}_4^1}\right){\mathit{u}}_4+\\ \hfill +{\mathit{p}}_1{\mathit{ds}}_{{\mathit{x}}_1^1}-\mathit{\mu }\frac{\partial \mathit{u}}{\partial \mathit{x}}{\mathit{ds}}_{{\mathit{x}}_1^1}-\mathit{\mu }\frac{\partial \mathit{u}}{\partial \mathit{y}}{\mathit{ds}}_{{\mathit{y}}_4^1}=0\end{array}$

(11) $\begin{array}{c}\frac{\mathit{p}{\mathit{\nu }}_1\left({\mathit{V}}_1-{\mathit{V}}_1^0\right)}{\Delta \mathit{t}}+\left(\mathit{p}{\hat{\mathit{u}}}_1{\mathit{ds}}_{{\mathit{x}}_1^1}\right){\mathit{v}}_1+\left(\mathit{\rho }{\hat{\mathit{v}}}_4{\mathit{ds}}_{{\mathit{y}}_4^1}\right){\mathit{v}}_4+\hfill \\ +{\mathit{p}}_4{\mathit{ds}}_{{\mathit{y}}_4^1}-\mathit{\mu }\frac{\partial \mathit{v}}{\partial \mathit{x}}{\mathit{ds}}_{{\mathit{x}}_1^1}-\mathit{\mu }\frac{\partial \mathit{v}}{\partial \mathit{y}}{\mathit{ds}}_{{\mathit{y}}_4^1}=0\hfill \end{array}$

where: ν₁ is the volume of ordinary SCV 1; the superscript " º " denotes value from the previous time step; ∆t is the time step. Lower and upper numeric indices in the normal surface vector components, for instance ds_x21, denote that ds_x is calculated on sub-surface 1 for the SCV 2. Similar equations can be obtained for other ordinary SCVs in the domain, e.g., SCV 4 in this element.

Next the implementation of Eq. 3 on ghost SCVs is explained. Ghost SCVs 2 and 3 are type I. The grey area in Fig. 6 represents the ghost SCV 2 in the flow field. This is an "effective" volume of ghost SCV 2 denoted by ν₂ this part. Substituting these parameters in Eq. 3 for SCV 2 it results in:

(12) $\begin{array}{c}\hfill {\int }_{\mathit{SCV}2}\frac{\partial \mathit{Q}}{\partial \mathit{t}}\mathit{dV}+{\int }_{\mathit{SSS}1\&\mathit{SS}{2}_{\mathit{r}}\&\mathit{SSb}}\left({\mathit{Eds}}_{\mathit{x}}+{\mathit{Fds}}_{\mathit{y}}\right)-\\ \hfill -{\int }_{\mathit{SS}!\&\mathit{SS}{2}_{\mathit{r}}\&\mathit{SSb}}\left({\mathit{Gds}}_{\mathit{x}}+{\mathit{Hds}}_{\mathit{y}}\right)=0\end{array}$

On SSb, flux vectors E, F, G, and H are evaluated on ip_b. These flux vectors are evaluated for SS2 on ip₂. The discrete form of Eq. 12 is given by

(13) $\begin{array}{c}\mathit{p}{\hat{\mathit{u}}}_1{\mathit{ds}}_{{\mathit{x}}_1^2}+\mathit{\rho }{\hat{\mathit{v}}}_2{\mathit{ds}}_{{\mathit{y}}_{2\mathit{r}}^2}+\left(\mathit{\rho }\overrightarrow{{\hat{\mathit{q}}}_{\mathit{b}}}.\overrightarrow{{\mathit{ds}}_{\mathit{b}}}\right)=0\hfill \\ \frac{\mathit{p}{\mathit{\nu }}_2\left({\mathit{U}}_2-{\mathit{U}}_2^0\right)}{\Delta \mathit{t}}+\left(\mathit{p}{\hat{\mathit{u}}}_1{\mathit{ds}}_{{\mathit{x}}_1^2}\right){\mathit{u}}_1+\left(\mathit{\rho }{\hat{\mathit{v}}}_2{\mathit{ds}}_{{\mathit{y}}_{2\mathit{r}}^2}\right){\mathit{u}}_2+\hfill \end{array}$

(14) $\begin{array}{c}+\left(\mathit{\rho }\overrightarrow{{\hat{\mathit{q}}}_{\mathit{b}}}.\overrightarrow{{\mathit{ds}}_{\mathit{b}}}\right){\mathit{u}}_{\mathit{b}}+{\mathit{p}}_1{\mathit{ds}}_{{\mathit{x}}_1^2}+{\mathit{p}}_{\mathit{b}}{\mathit{ds}}_{\mathit{xb}}-\hfill \\ {{\displaystyle -\mathit{\mu }\frac{\partial \mathit{u}}{\partial \mathit{x}}\mid }}_1{\mathit{ds}}_{{\mathit{x}}_1^2}-\mathit{\mu }{{\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{x}}\mid }}_{\mathit{b}}{\mathit{ds}}_{{\mathit{x}}_{\mathit{b}}}-\mathit{\mu }{{\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{y}}\mid }}_{\mathit{b}}{\mathit{ds}}_{\mathit{yb}}-\hfill \\ {{\displaystyle -\mathit{\mu }\frac{\partial \mathit{u}}{\partial \mathit{y}}\mid }}_2{\mathit{ds}}_{{\mathit{y}}_{2\mathit{r}}^2}=0\hfill \end{array}$

(15) $\begin{array}{c}\frac{{\mathit{\rho }\mathit{\nu }}_2\left({\mathit{V}}_2-{\mathit{V}}_2^0\right)}{\Delta \mathit{t}}+\left(\mathit{\rho }{\hat{\mathit{u}}}_1{\mathit{ds}}_{{\mathit{x}}_1^2}\right){\mathit{v}}_1+\left(\mathit{\rho }{\hat{\mathit{v}}}_2{\mathit{ds}}_{{\mathit{y}}_{2\mathit{r}}^2}\right){\mathit{v}}_2+\hfill \\ +\left(\mathit{\rho }\overrightarrow{{\hat{\mathit{q}}}_{\mathit{b}}}.\overrightarrow{{\mathit{ds}}_{\mathit{b}}}\right){\mathit{v}}_{\mathit{b}}+{\mathit{p}}_2{\mathit{ds}}_{{\mathit{y}}_2^2}+{\mathit{p}}_{\mathit{b}}{\mathit{ds}}_{{\mathit{y}}_{\mathit{b}}}-\hfill \\ {{\displaystyle -\mathit{\mu }\frac{\partial \mathit{v}}{\partial \mathit{y}}\mid }}_1{\mathit{ds}}_{{\mathit{x}}_1^2}-\mathit{\mu }{{\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{x}}\mid }}_{\mathit{b}}{\mathit{ds}}_{{\mathit{x}}_{\mathit{b}}}-\hfill \\ {{\displaystyle -\mathit{\mu }\frac{\partial \mathit{v}}{\partial \mathit{y}}\mid }}_{\mathit{b}}{\mathit{ds}}_{\mathit{yb}}-\mathit{\mu }{{\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{y}}\mid }}_2\mathit{ds}{}_{{\mathit{y}}_{2\mathit{r}}^2}=0\hfill \end{array}$

where q̂_b = (û_{2 b} + v̂_2b)^1/2 is the convecting velocity vector; $\overrightarrow{\mathit{d}}{\mathit{s}}_{\mathit{b}}$ is the normal surface vector in the outward direction and ds_xb and ds_yb are the components of $\overrightarrow{\mathit{d}}{\mathit{s}}_{\mathit{b}}$ in x and y directions, respectively.

Depending on the boundary condition, appropriate constraints can be forced in Eqs. 13 to 15. For instance, if the boundary b is solid, then ($\left(\overrightarrow{\mathit{p}}{\hat{\mathit{q}}}_{\mathit{b}}·\overrightarrow{\mathit{d}}{\mathit{s}}_{\mathit{b}}\right)$ ) = 0, u_b = 0 and v_b = 0; p_b is described based on the nodal pressures of element using bi-linear interpolation. Moreover, velocity gradients of ${{\displaystyle \frac{\partial u}{\partial x}\mid }}_b,{{\displaystyle \frac{\partial u}{\partial y}\mid }}_b,{{\displaystyle \frac{\partial v}{\partial x}\mid }}_b,\mathrm{and}{{\displaystyle \frac{\partial v}{\partial y}\mid }}_b,$ are evaluated using bilinear interpola- ion defined in Eq. 2.

Ghost Sub-Control Volume Type II

In Fig. 7 an element with 2 fluid nodes and two ghost nodes is shown. As mentioned in the section "Sub-Control Volumes and Node Types", SCVs 1 and 4 are considered ordinary SCVs, and SCVs 2 and 3 are ghost SCVs type II. Since it is important to remain in the IBM general framework, any point within the flow field (i.e. inside the boundary) and its SCV are considered ordinary. Here Eq. 3 is applied to the whole area of SCV 1, i.e. the area between SS1, SS4, and node 1. Conservation laws for an ordinary SCV were introduced by Eqs. 9-11 in section ghost SCV type I. The actual area within the flow field is the dotted area between SSb, SS4, and node 1 which is shown by grey area in Fig. 6. This area is assigned to ghost node 2 and is called ghost SCV 2. Conservation laws (Eq. 3) are written for this ghost SCV, and boundary condition is applied in these equations. In the present study, boundary condition is applied via ghost SCVs, and not necessarily via the SCVs containing the boundary. Combination of conservation laws for ordinary SCV 1 and ghost SCV 2 will result in the conservation of conserved quantities for the dotted area in SCV 1, which is actually within the flow filed. Implementation of Eq. 3 on ghost SCVs is explained next. Mass conservation equation for the grey area is written as:

Figure 7
Ghost SCV type II (SCV 1 and SCV 2 are considered); the grey area is ghost-SCV type II assigned to node 2; the dotted area is the difference area between complete SCV 1 area and the grey area which contains fl uid; 4r is the right part of sub-surface 4; 4l is the left part of sub-surface 4 along the grey area; SSb is the boundary portion in SCV 1; v₂ is volume of the grey area; v₁ is complete volume of SCV 1; ds_b^h is normal surface vector of boundary in SCV 1 in the direction outward from the grey area; ds_b^d is the normal surface vector on SSb in the direction outward from the dotted area; ds_x¹₁ is normal surface vector of sub-surface 1; ds_y¹_4l is normal surface vector related to the left part of sub-surface 4; 1,2,3, and 4 points indicated with cross symbol are ip.

(16) $\mathit{p}{\hat{\mathit{u}}}_1{\mathit{ds}}_{{\mathit{x}}_1^1}+\mathit{p}{\hat{\mathit{v}}}_4{\mathit{ds}}_{{\mathit{y}}_{4\mathit{l}}^1}+\left(\mathit{p}\overrightarrow{{\hat{\mathit{q}}}_{\mathit{b}}}.{\overrightarrow{\mathit{ds}}}_{\mathit{b}}\mathit{h}\right)=0$

where: $\overrightarrow{\mathit{d}}{\mathit{s}}_{{\mathit{b}}^{\mathit{h}}}$ is the normal surface vector on SSb. As shown in Fig. 6 surface 4 is divided into 2 parts where the left side is denoted by 4l and the right side denoted by 4r. Mass conservation equations of SCV 1 and SCV 2 are written in the system of equations, and solved simultaneously. In order to obtain the final solution of this method, the 2 following equations are combined:

(17) $\begin{array}{c}\mathit{p}{\hat{\mathit{v}}}_4{\mathit{ds}}_{{\mathit{y}}_{4\mathit{r}}^1}-\left(\mathit{p}\overrightarrow{{\hat{\mathit{q}}}_{\mathit{b}}}.{\overrightarrow{\mathit{ds}}}_{{\mathit{b}}^{\mathit{h}}}\right)=0\hfill \\ \mathit{p}{\hat{\mathit{v}}}_4{\mathit{ds}}_{{\mathit{y}}_{4\mathit{r}}^1}-\left(\mathit{p}\overrightarrow{{\hat{\mathit{q}}}_{\mathit{b}}}.{\overrightarrow{\mathit{ds}}}_{{\mathit{b}}^{\mathit{d}}}\right)=0\hfill \end{array}$

where $\overrightarrow{\mathit{d}}{\mathit{s}}_{{\mathit{b}}^{\mathit{d}}}$ is the normal surface vector on SSb and is equal to $-\overrightarrow{\mathit{d}}{\mathit{s}}_{{\mathit{b}}^{\mathit{h}}}$.

Equation 17 is in fact the mass conservation equation for the dotted area with the actual SCV related to node 1. Depending on the boundary condition of the problem, appropriate constraints can be forced in Eq. 16. For instance, if boundary b is solid, then ($\left(\overrightarrow{\mathit{p}}{\hat{\mathit{q}}}_{\mathit{d}}·\overrightarrow{{\mathit{ds}}_{{\mathit{b}}^{\mathit{h}}}}\right)$) = 0.

Similar procedure is applied for momentum conservation equations. The x-momentum conservation equation for the grey area in Fig. 6 is written as follows:

(18) $\begin{array}{c}\frac{\mathit{p}{\mathit{\nu }}_2\left({\mathit{U}}_2-{\mathit{U}}_2^0\right)}{\Delta \mathit{t}}+\left(\mathit{p}{\hat{\mathit{u}}}_1{\mathit{ds}}_{{\mathit{x}}_1^1}\right){\mathit{u}}_1+\left(\mathit{p}{\hat{\mathit{v}}}_4{\mathit{ds}}_{{\mathit{y}}_{4\mathit{l}}^1}\right){\mathit{u}}_4+\hfill \\ +\left(\mathit{p}\overrightarrow{{\hat{\mathit{q}}}_{\mathit{b}}}.{\overrightarrow{\mathit{ds}}}_{{\mathit{b}}^{\mathit{h}}}\right)\mathit{u}{}_{\mathit{b}}+{\mathit{p}}_1{\mathit{ds}}_{{\mathit{x}}_1^1}+{\mathit{p}}_{\mathit{b}}{\mathit{ds}}_{{\mathit{x}}_{\mathit{b}}^{\mathit{h}}}-\hfill \\ {{\displaystyle -\mathit{\mu }\frac{\partial \mathit{u}}{\partial \mathit{x}}\mid }}_1{\mathit{ds}}_{{\mathit{x}}_1^1}-\mathit{\mu }{{\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{x}}\mid }}_{{\mathit{b}}^{\mathit{h}}}{\mathit{ds}}_{{\mathit{x}}_{\mathit{b}}^{\mathit{h}}}-\hfill \\ {{\displaystyle -\mathit{\mu }\frac{\partial \mathit{u}}{\partial \mathit{y}}\mid }}_{{\mathit{b}}^{\mathit{h}}}{\mathit{ds}}_{\mathit{yb}}^{\mathit{h}}{{\displaystyle -\mathit{\mu }\frac{\partial \mathit{u}}{\partial \mathit{y}}\mid }}_2{\mathit{ds}}_{{\mathit{y}}_{4\mathit{l}}^1}=0\hfill \end{array}$

Here ν₂ is the volume of grey area in Fig. 6. As mentioned for the mass conservation equation, Eqs. 10 and 18 are written in the system of equations and solved simultaneously. In order to obtain the final solution of this method, the 2 following equations are combined:

(19) $\begin{array}{c}\frac{{\mathit{pv}}_1\left({\mathit{U}}_1-{\mathit{U}}_1^0\right)}{\Delta \mathit{t}}-\frac{\mathit{\rho }{\mathit{v}}_2\left({\mathit{U}}_2-{\mathit{U}}_2^0\right)}{\Delta \mathit{t}}+\left(\mathit{\rho }{\hat{\mathit{v}}}_4{\mathit{ds}}_{{\mathit{y}}_{4\mathit{r}}^1}\right){\mathit{u}}_4-\hfill \\ -\left(\mathit{\rho }\overrightarrow{{\hat{\mathit{q}}}_{\mathit{b}}}.{\overrightarrow{\mathit{ds}}}_{{\mathit{b}}^{\mathit{h}}}\right){\mathit{u}}_{\mathit{b}}-{\mathit{p}}_{\mathit{b}}{\mathit{ds}}_{{\mathit{x}}_{\mathit{b}}^{\mathit{h}}}-{{\displaystyle \mathit{\mu }\frac{\partial \mathit{u}}{\partial \mathit{x}}\mid }}_{{\mathit{b}}^{\mathit{h}}}{\mathit{ds}}_{{\mathit{x}}_{\mathit{b}}^{\mathit{h}}}-\hfill \\ -{{\displaystyle \mathit{\mu }\frac{\partial \mathit{u}}{\partial \mathit{y}}\mid }}_{{\mathit{b}}^{\mathit{h}}}{\mathit{ds}}_{{\mathit{y}}_{\mathit{b}}^{\mathit{h}}}-{{\displaystyle \mathit{\mu }\frac{\partial \mathit{u}}{\partial \mathit{y}}\mid }}_2{\mathit{ds}}_{{\mathit{y}}_{4\mathit{r}}^1}=0\hfill \end{array}$

After substitution, Eq. 19 will become:

(20) $\begin{array}{c}\left\{\frac{{\mathit{\rho }\mathit{\nu }}_1\left({\mathit{U}}_1-{\mathit{U}}_1^0\right)}{\Delta \mathit{t}}-\frac{{\mathit{\rho }\mathit{\nu }}_2\left({\mathit{U}}_2-{\mathit{U}}_2^0\right)}{\Delta \mathit{t}}\right\}+\left(\mathit{\rho }{\hat{\mathit{v}}}_4{\mathit{ds}}_{{\mathit{y}}_{4\mathit{r}}^1}\right){\mathit{u}}_4+\hfill \\ +\left(\mathit{\rho }\overrightarrow{{\hat{\mathit{q}}}_{\mathit{b}}}.{\overrightarrow{\mathit{ds}}}_{{\mathit{b}}^{\mathit{d}}}\right){\mathit{u}}_{\mathit{b}}+{\mathit{p}}_{\mathit{b}}{\mathit{ds}}_{{\mathit{x}}_{\mathit{b}}^{\mathit{d}}}+{{\displaystyle \mathit{\mu }\frac{\partial \mathit{u}}{\partial \mathit{x}}\mid }}_{{\mathit{b}}^{\mathit{h}}}{\mathit{ds}}_{{\mathit{x}}_{\mathit{b}}^{\mathit{h}}}+\hfill \\ +{{\displaystyle \mathit{\mu }\frac{\partial \mathit{u}}{\partial \mathit{y}}\mid }}_{{\mathit{b}}^{\mathit{h}}}{\mathit{ds}}_{{\mathit{y}}_{\mathit{b}}^{\mathit{h}}}-{{\displaystyle \mathit{\mu }\frac{\partial \mathit{u}}{\partial \mathit{y}}\mid }}_2{\mathit{ds}}_{{\mathit{y}}_{4\mathit{r}}^1}=0\hfill \end{array}$

This is the equation of x-momentum conservation for the dotted area which is the actual SCV of node 1, if U₁ and U₂ are assumed to be the same in transient term. Here (ν₁ - ν₂) is equal to the volume of the dotted area. Conservation equation of y-momentum for the grey area in Fig. 6 can be obtained from the following equation:

(21) $\begin{array}{c}\frac{{\mathit{\rho }\mathit{\nu }}_2\left({\mathit{V}}_2-{\mathit{V}}_2^0\right)}{\Delta \mathit{t}}+\left(\mathit{\rho }{\hat{\mathit{u}}}_1{\mathit{ds}}_{{\mathit{x}}_1^1}\right){\mathit{v}}_1+\left(\mathit{\rho }{\hat{\mathit{v}}}_4{\mathit{ds}}_{{\mathit{y}}_{4\mathit{l}}^1}\right){\mathit{v}}_4+\hfill \\ +\left(\mathit{\rho }\overrightarrow{{\hat{\mathit{q}}}_{\mathit{b}}}.{\overrightarrow{\mathit{ds}}}_{{\mathit{b}}^{\mathit{h}}}\right){\mathit{v}}_{\mathit{b}}+{\mathit{p}}_{´4}{\mathit{ds}}_{{\mathit{y}}_4^1}+{\mathit{p}}_{\mathit{b}}{\mathit{ds}}_{{\mathit{y}}_{\mathit{b}}^{\mathit{h}}}-\hfill \\ -{{\displaystyle \mathit{\mu }\frac{\partial \mathit{v}}{\partial \mathit{x}}\mid }}_{{\mathit{b}}^{\mathit{h}}}{\mathit{ds}}_{{\mathit{x}}_{\mathit{b}}^{\mathit{h}}}-{{\displaystyle \mathit{\mu }\frac{\partial \mathit{v}}{\partial \mathit{y}}\mid }}_{{\mathit{b}}^{\mathit{h}}}{\mathit{ds}}_{{\mathit{y}}_{\mathit{b}}^{\mathit{h}}}-{{\displaystyle \mathit{\mu }\frac{\partial \mathit{v}}{\partial \mathit{y}}\mid }}_2{\mathit{ds}}_{{\mathit{y}}_{4\mathit{l}}^1}=0\hfill \end{array}$

Since the boundary b is solid, then momentum flux $\left(\overrightarrow{\mathit{p}}{\hat{\mathit{q}}}_{\mathit{b}}·\overrightarrow{\mathit{d}}{\mathit{s}}_{{\mathit{b}}^{\mathit{h}}}\right){\mathit{u}}_{\mathit{b}}$ would be 0; p_b is evaluated using bilinear interpolation, and velocity gradient terms of ${{\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{x}}\mid }}_{{\mathit{b}}^{\mathit{h}}}\mathit{and}{{\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{y}}\mid }}_{{\mathit{b}}^{\mathit{h}}}$ are evaluated using bilinear interpolation. Identical procedure can be applied to obtain a similar equation to Eq. 20 for y-momentum conservation equation.

The key-point of the present method is to remain in the IBM context while implementing conservation laws all over the domain including boundary control volumes. Therefore control volumes of nodes within the flow field are always considered ordinary and complete. In this method, boundary conditions are implemented via the ghost control volumes where conservation laws are applied to determine variables values on ghost nodes. This is in contrast to other IBM methods in which interpolation functions (Mittal and Iaccarino 2005; Majumdar et al. 2001; Ghias et al. 2004), or cut cell methods (Clarke et al. 1986; Udaykumar et al. 2001, 1999, 1996; Ye et al. 1999) are used to determine variables values of ghost nodes to satisfy boundary conditions. The application of this method can be extended to moving boundary problems in IBM context to reduce the spurious pressure oscillations. This is due to local mass conservation errors observed in simulations of moving boundary problems with typical immersed boundary methods (Seo and Mittal 2011).

RESULTS

Taylor problem (Alisadeghi 2012; Alisadeghi and Karimian 2011; Mahesh et al. 2004; Darbandi and Vakilipour 2008) is selected to evaluate the present method. This problem corresponds to periodic and counter-rotating vortices whose strength decays in time at a rate determined by the viscosity. The Exact solution for Taylor problem (Alisadeghi 2012; Alisadeghi and Karimian 2011; Mahesh et al. 2004; Darbandi and Vakilipour 2008) for velocity components, pressure, and kinetic energy at square domain of unit size is given as:

(22) $\mathit{u}\left(\mathit{x},\mathit{y}\right)=\sin \left(\mathit{\pi }\mathit{x}\right)\cos \left(\mathit{\pi }\mathit{y}\right)\mathit{e}\left(-2{\mathit{\pi }}^2\raisebox{1ex}{$\mathit{t}$}\!\left/ \!\raisebox{-1ex}{$\mathit{Re}$}\right.\right)$

(23) $\mathit{v}\left(\mathit{x},\mathit{y}\right)=-\cos \left(\mathit{\pi }\mathit{x}\right)\sin \left(\mathit{\pi }\mathit{y}\right)\mathit{e}\left(-2{\mathit{\pi }}^2\raisebox{1ex}{$\mathit{t}$}\!\left/ \!\raisebox{-1ex}{$\mathit{Re}$}\right.\right)$

(24) $\begin{array}{c}\hfill \mathit{P}\left(\mathit{x},\mathit{y}\right)=0.25\left[\cos \left(2\mathit{\pi }\mathit{x}\right)\right.+\\ \hfill \left.+\cos \left(2\mathit{\pi }\mathit{y}\right)\right]\mathit{e}\left(\raisebox{1ex}{$-4{\mathit{\pi }}^2\mathit{t}$}\!\left/ \!\raisebox{-1ex}{$\mathit{Re}$}\right.\right)\end{array}$

(25) $\begin{array}{c}\mathit{KE}\left(\mathit{x},\mathit{y}\right)=0.5\left[{\sin }^2\right.\left(\mathit{\pi }\mathit{x}\right){\cos }^2\left(\mathit{\pi }\mathit{y}\right)+\hfill \\ \left.+{\sin }^2\left(\mathit{\pi }\mathit{y}\right){\cos }^2\left(\mathit{\pi }\mathit{x}\right)\right]\mathit{e}\left(\raisebox{1ex}{$-4{\mathit{\pi }}^2\mathit{t}$}\!\left/ \!\raisebox{-1ex}{$\mathit{Re}$}\right.\right)\hfill \end{array}$

In inviscid case where Reynolds number is infinite, the exponential terms in Eqs. 22-25 can be ignored. Pressure contours and streamlines of the exact solution for inviscid Taylor problem are plottet in Fig. 8. The solution of the flowfield will be constant with time if numerical solution is used in inviscid case. For the first test problem, the present method is applied to investigate this fact. Figure 9 shows grid structure and domain boundary used for the present method. Grid spacing in both directions are uniform and equal to 0.05.

Figure 8
Pressure contours and streamlines of the exact solution for inviscid Taylor problem.

Figure 9
Grid structure, domain boundary (blackline), and ghost nodes (black nodes) for the solution of Taylor problem.

Domain boundary, denoted by blackline, is immerssed within the elements close to boundary. Thus, domain size will be 0.93 by 0.93, which is less than unity. In this case, the outer grid nodes, depicted by the black squares, are ghost nodes. SCVs of these ghost nodes are type II (as Fig. 4). Boundary conditions are implemented in conservation equations of mass and momentums for ghost control volumes. In the mass conservation equation (Eq. 16), boundary mass flux ($\left(\overrightarrow{\mathit{p}}{\hat{\mathit{q}}}_{\mathit{d}}·\overrightarrow{{\mathit{ds}}_{{\mathit{b}}^{\mathit{h}}}}\right)$) is known. In x-momemtum conservation equation (Eq. 18), $\left(\overrightarrow{\mathit{p}}{\hat{\mathit{q}}}_{\mathit{d}}·\overrightarrow{{\mathit{ds}}_{{\mathit{b}}^{\mathit{h}}}}\right){\mathit{u}}_{\mathit{b}}$ is also calculated from known values of the exact solution; p_b is evaluated from bilinear interpolation from the nodal values of pressure. The same method is used for y-momemtum conservation equation. In the inviscid Taylor problem, solution is initially proceeded for the fist hundered time using Δt = 1e−6 s. Then velocity components calculated by the present method are compared with those of the exact solution along the midlines and diagonal of the domain shown in Fig. 10. As can be observed the present results have remained constant throughout the time and exactly the same as the results of the exact solution.

Figure 10
Velocity diagrams of Taylor problem. Comparison with the exact solution: (a) u and v along the horizontal and the vertical midlines; (b) u and v along the diagonal from south west to northeast.

To obtain the difference between the results of the exact solution and the present method we introduce the following equations for P and KE:

(26) ${\mathit{E}}_{\mathit{p}}=\sum _1^{\mathit{N}}\left(\right.\frac{\sqrt[2]{{\left({\mathit{P}}_{\mathit{n}}-{\mathit{P}}_{\mathit{e}}\right)}^2}}{{\mathit{P}}_{\mathit{e}}}\left./\mathit{N}\times 100\right)$

(27) ${\mathit{E}}_{\mathit{KE}}=\sum _1^{\mathit{N}}\left(\right.\frac{\sqrt[2]{{\left(\mathit{KEn}-{\mathit{KE}}_{\mathit{e}}\right)}^2}}{{\mathit{KE}}_{\mathit{e}}}\left./\mathit{N}\times 100\right)$

where N is the number of fluid nodes in the domain; subscripts n and e indicate numerical and the exact solutions, respectively. P and KE errors reached 0.004 and 9.7e−7%, respectively, and their values did not change during the first 100 time iterations. Error between the results of numerical and the exact solutions shows the high strenght of the present method.

To present the capability of the present method in handling immersed boundaries, Taylor problem is solved on a grid rotated 45º with respect to the solution domain as depicted in Fig. 11. Uniform grid spacing in both directions equal to 0.0707 m. Domain boundary, denoted by the black line in Fig. 11, is immerssed within the elements. Therefore, domain size will be 0.85 by 0.85, which is less than 1 unity. In this case, the outer grid nodes, denoted by the black color in Fig. 12, are ghost nodes. SCVs of these ghost nodes (Fig. 13) are type I and II, as previously defined.

Figure 11
Schematic of Taylor problem (blackline) on a grid rotated 45º.

Figure 12
Positions of ghost nodes (blacknodes) with respect to solution domain (blackline).

Figure 13
SCVs of ghost nodes of Fig. 12. (a) Type I; (b) Type II.

Boundary conditions are implemented using conservation equations of mass and momentums of ghost control volumes. Implementation of conservation equations for ghost SCV type II (Fig. 13b) is similar to the one in the previous test case. For the ghost subcontrol volume type I (Fig. 13a) a similar procedure is followed. This means that boundary mass flux ($\left(\overrightarrow{\mathit{p}}{\hat{\mathit{q}}}_{\mathit{d}}·\overrightarrow{{\mathit{ds}}_{{\mathit{b}}^{\mathit{h}}}}\right)$) in Eq. 13 is calculated by known values of the exact solution. In x-momemtum conservation equation introduced by Eq. 14, $\left(\overrightarrow{\mathit{p}}{\hat{\mathit{q}}}_{\mathit{d}}·\overrightarrow{{\mathit{ds}}_{{\mathit{b}}^{\mathit{h}}}}\right){\mathit{u}}_{\mathit{b}}$ is calculated using known values of the exact solution; p_b is calculated by bilinear enterpolation using the element nodal values; 1e−6 s is used as time step for the first 100 time steps.

Streamlines obtained from the present method are plotted in Fig. 14. Note that 4 triangles around the central square domain in Fig. 14 are outside the flowfield. In general, streamlines in Fig. 14 are similar to those of Fig. 8. Errors for pressure and kinetic energy, E_p and E_KE, were 2.24 and 4.5e−6%, respectively. Similar to the previous test case, KE keeps constant in time as the problem is non-viscous.

Figure 14
Pressure contours and streamlines for the proposed IBM (Taylor problem with 45º grid rotation).

Velocity profiles derived from the present method are compared with the results of the exact solution in the middle of the flow field and along the diagonal from south west to north east in Fig. 15, and along left, right, up, and down boundaries in Fig. 16. The results show very good agreement with the results of the exact solution. Therefore it is seen that, in cases where grid structure is not aligned with domain boundary, the present method is still able to solve flowfield accurately. The differences between profiles in first and second test cases in Figs. 10 and 14 are due to the difference between sizes of solution domains in the 2 cases.

Figure 15
Velocity diagrams of Taylor problem with 45º grid rotation. Comparison with the exact solution: (a) u and v along the horizontal and the vertical midlines; (b) Velocity components along the diagonal from southwest to northeast.

Figure 16
Velocity diagrams of Taylor problem with 45º grid rotation. Comparison with the exact solution: (a) u and v along the left and right boundaries; (b) u and v along the up and down boundaries.

As our third test case, Taylor problem with Reynolds number of 1,000 is solved using the present method on regular grid domain. This case is chosen to investigate the capability of the present method to solve viscous transient problems. Schematic of the solution domain is similar to Fig. 9, but the location of the boundary is different as shown in Fig. 17. In the viscous case, pressure and velocity fields decay in time at a rate determined by the viscosity. As verified in Eq. 25, KE is a function of t/Re. So the higher the Reynolds number, the lower the rate of KE decay.

Figure 17
Grid structure for the solution of viscous Taylor problem on regular grid, domain boundary by black line, and ghost nodes by black nodes.

This test is solved using the present method. Ghost control volumes in Fig. 17 is type I and is considered as described before. In this case, viscosity terms should be calculated at the boundary. Therefore, ${{\displaystyle \mathit{\mu }\frac{\partial u}{\partial x}\mid }}_{b^h}{\mathrm{ds}}_{x_b^h}\mathrm{and}{{\displaystyle \mathit{\mu }\frac{\partial u}{\partial y}\mid }}_{b^h}{\mathrm{ds}}_{y_b^h}$ and can be found using known values of the exact solution. All the other terms in mass and momentum conservetion equations are calculated in a similar way as explained in the previous test case, although grid configuration is different from the previous test case. Solution proceeds until the pressure and velocity fields are completely decayed by viscosity. In this case, Δt = 0.005 s and the solution is itereted in every time step until the solution converges. Decay of the temporal kinetic energy is compared with the exact solution in Fig. 18. The result of the present method has excellent agreement with that of the exact solution. Grid independency study was carried out in present research. Viscous Taylor problem on regular grid was solved on 21 × 21, 31 × 31, and 41 × 41 grids. As a sample chosen from the results, pressure variations along the diagonal of the solution domain on 3 grids are compared with each other in Fig. 19. As can be seen, after a small change between results of grid 21 × 21 and grid 31 × 31, no significant changes can be notified between the results of grids 31 × 31 and 41 × 41. The same trend was observed on other results of these 3 grids. So all results are obtained on grid 31 × 31.

Figure 18
Decay of temporal KE: comparison between results of the present method and the exact solution for Reynolds number of 1,000 on a regular grid.

Figure 19
Pressure variation along the diameter for the purpose of grid convergence study. Grids with 21 × 21, 31 × 31 and 41 × 41 nodes.

KE decays faster as Reynolds number decreases. Decay of temporal KE are compared with the exact solution for 2 Reynolds numbers of 100 and 1,000 in Fig. 20. The results show very good agreement with the exact solution. As observed, fast decay of KE clearly occurs for Reynolds number of 100.

Figure 20
Decay of temporal KE for 2 Reynolds numbers of 100 and 1,000. Comparison with the exact solution.

Velocity profiles resulted from the present method are compared with the results of the exact solution in the middle of the flow field, and along the diagonal from southwest to northeast in Fig. 21. Results are presented at different times of of 0.0, 25, and 50 s. Again, as observed, all of the results match excellently with the exact solution obtained from Eqs. 22 and 23.

Figure 21
Velocity profiles of viscous Taylor problem in 3 instances of 0.0, 25, and 50.0 s. Comparison with the exact solution: (a) u and v along the horizontal and vertical midlines; (b) u and v along the diagonal from southwest to northeast.

In addition to the above comparisons, rates of pressure at the center of the solution domain and average of pressure over the solution domain with time are very similar to results of the exact solution as shown in Fig. 22.

Figure 22
Time variation of pressure at the center of the solution domain and average of pressure over the solution domain: comparison of the present method with the exact solution.

Here again, viscous Taylor problem is solved on a 45º rotated grid to present capability of present immersed boundary method on grids skewed with respect to solu- tion boundary. Solution domain and grid structure of this fourth test case are shown in Fig. 23. The only difference between Figs. 12 and 23 is the location of the boundary. In Fig. 23 the domain size is exactly 1 × 1 square, which is smaller in Fig. 12. In both cases the type of the ghost control volmes are the same as shown in Fig. 13, and a 31 × 31 grid is used in both cases. The outer grid nodes, denoted by black color in Fig. 23, are ghost nodes.

Figure 23
Grid structure for the solution of viscous Taylor problem on 45º rotated grid: domain boundary by black line and ghost nodes by black nodes.

Boundary conditions are implemented through conserva- tion equations of mass and momentums of ghost control volumes. Implementation of conservation equations for ghost subcontrol volumes in this case is exactly the same as the one applied in the second case with its grid in Fig. 13. The only difference in this case is the fact that viscosity terms are calculated at the boundary. Therefore, ${{\displaystyle \mathit{\mu }\frac{\partial \mathit{u}}{\partial \mathit{x}}\mid }}_{{\mathit{b}}^{\mathit{h}}}{\mathit{ds}}_{{\mathit{x}}_{\mathit{b}}^{\mathit{h}}}\mathit{and}{{\displaystyle \mathit{\mu }\frac{\partial \mathit{u}}{\partial \mathit{y}}\mid }}_{{\mathit{b}}^{\mathit{h}}}{\mathit{ds}}_{{\mathit{y}}_{\mathit{b}}^{\mathit{h}}}$ are calculated from known values of the exact solution. Solution is proceeded until velocity and pressure fields completely decay with time and is proceeded in time using time steps of 1e−3 s.

Decay of temporal KE and its comparison with the exact solution are illustrated in Fig. 24. Result of the present method follows the exact solution of KE decay. Velocity profiles derived from the present method are compared with the results of the exact solution in the middle of the flow field and along the diagonal from southwest to northeast in Fig. 25. Results are presented at various times of 0.0, 25, and 50 s. In all profiles excellent agreement is verified between results of the present method and those of the exact solution. Here again, rates of pressure at the center of the solution domain and the pressure average variations over the solution domain with time also agree very well with the results of the exact solution as shown in Fig. 26.

Figure 24
Decay of temporal KE: comparison between results of the present method and the exact solution for Reynolds number of 1,000 on a 45º rotated grid.

Figure 25
Velocity profiles of viscous Taylor problem in 3 instances of 0.0, 25, and 50 s. Comparison with the exact solution: (a) u and v along the horizontal and vertical midlines; (b) u and v along the diagonal from southwest to northeast.

Figure 26
Time variation of pressure at the center of the solution domain and average of pressure over the solution domain: comparison of the present method with the exact solution.

Now it is possible to claim with confidence that the present method can solve the flow filed in both viscous and non-viscous cases with high accuracy in comparison with the exact solution. The accuracy of the present method is not challenged even on grids not aligned with the boundary domain.

CONCLUSION

In this paper a new immersed boundary method using control volume finite element scheme was introduced for discretization of governing equations. The advantage of this method is that boundary conditions are implemented by conservation of conserved quantities along the boundaries. In typical immersed boundary methods conservation equations are only satisfied within the boundary but not necessarily along the boundary. However, in this study, a new approach for the implementation of boundary conditions was presented in which mass and momentum conservation laws are fully conserved along the boundary as well as inside the domain. In the present method, in addition to the use of ghost node value applied in typical IBM, a new concept of ghost SCV was introduced. This new concept makes the implementation of the conservation laws in the vicinity of the boundary possible. The present method is validated by solving Taylor problem in both non-viscous and viscous cases with Reynolds numbers of 100 and 1,000. Steady and unsteady cases of Taylor problem on regular and 45º rotated grids were also solved for further investigation. Results both in pressure and velocity diagrams show an excellent agreement between the present method and of the exact solution in all cases even at the sharp corners. Based on these results, accurate simulation of the flow fields in physically complex problems can be expected from the present IBM method.

REFERENCES

• Alisadeghi H (2012) Modification and improvement of a pressure-velocity coupling algorithm for the numerical solution of two dimensional incompressible flows (PhD thesis). Tehran: Amirkabir University of Technology. In Persian.
• Alisadeghi H, Karimian SMH (2011) Different modelings of cell-face velocities and their effects on the pressure-velocity coupling, accuracy and convergence of solution. Int J Numer Meth Fluid 65(8):969-988. doi: 10.1002/fld.2224
  » https://doi.org/10.1002/fld.2224
• Berger M, Aftosmis M (1998) Aspects (and aspect ratios) of Cartesian mesh methods. In: Bruneau CH, editor. Sixteenth International Conference on Numerical Methods in Fluid Dynamics. Vol. 515. Berlin: Springer.
• Beyer Jr RP (1992) A computational model of the cochlea using the immersed boundary method. J Comput Phys 98(1):145-162. doi: 10.1016/0021-9991(92)90180-7
  » https://doi.org/10.1016/0021-9991(92)90180-7
• Clarke DK, Hassan HA, Salas MD (1986) Euler calculations for multielement airfoils using Cartesian grids. AIAA J 24(3):353-358. doi: 10.2514/3.9273
  » https://doi.org/10.2514/3.9273
• Darbandi M, Vakilipour S (2008) Developing implicit pressure-weighted upwinding scheme to calculate steady and unsteady flows on unstructured grids. Int J Numer Meth Fluid 56(2):115-141. doi: 10.1002/fld.1451
  » https://doi.org/10.1002/fld.1451
• Fauci LJ, McDonald A (1995) Sperm motility in the presence of boundaries. Bull Math Biol 57(5):679-699. doi: 10.1016/0092-8240(95)00022-I
  » https://doi.org/10.1016/0092-8240(95)00022-I
• Ghias R, Mittal R, Dong H (2007) A sharp interface immersed boundary method for compressible viscous flows. J Comput Phys 225(1):528-553. doi: 10.1016/j.jcp.2006.12.007
  » https://doi.org/10.1016/j.jcp.2006.12.007
• Ghias R, Mittal R, Lund TS (2004) A non-body conformal grid method for simulation of compressible flows with complex immersed boundaries. Proceedings of the 42nd AIAA Aerospace Sciences Meeting and Exhibit; Reno, USA.
• Karimian S, Schneider G (1994a) Numerical solution of two-dimensional incompressible Navier-Stokes Equations: treatment of velocity-pressure coupling. Proceedings of the 25th AIAA Fluid Dynamics Conference; Colorado Springs, USA.
• Karimian S, Schneider G (1994b) Pressure-based computational method for compressible and incompressible flows. Journal of Thermodynamics and Heat Transfer 8(2):267-274. doi: 10.2514/3.533
  » https://doi.org/10.2514/3.533
• Khadra K, Angot P, Parneix S, Caltagirone JP (2000) Fictitious domain approach for numerical modeling of Navier-Stokes equations. Int J Numer Meth Fluid 34(8):651-684. doi: 10.1002/1097-0363(20001230)34:8<651::AID-FLD61>3.0.CO;2-D
  » https://doi.org/10.1002/1097-0363(20001230)34:8<651::AID-FLD61>3.0.CO;2-D
• Mahesh K, Constantinescu G, Moin P (2004) A numerical method for large-eddy simulation in complex geometries. J Comput Phys 197(1):215-240. doi: 10.1016/j.jcp.2003.11.031
  » https://doi.org/10.1016/j.jcp.2003.11.031
• Majumdar S, Iaccarino G, Durbin PA (2001) RANS solvers with adaptive structured boundary non-conforming grids. Annual Research Briefs. Stanford: Center for Turbulence Research
• Minkowycz WJ, Sparrow EM, Murthy JY (1988) Handbook of numerical heat transfer. Vols. 379-421. New York: Wiley.
• Mittal R, Dong H, Bozkurttas M, Najjar FM, Vargas A, von Loebbecke A (2008) A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries. J Comput Phys 227(10):4825-4852. doi: 10.1016/j.jcp.2008.01.028
  » https://doi.org/10.1016/j.jcp.2008.01.028
• Mittal R, Iaccarino G (2005) Immersed boundary methods. Annu Rev Fluid Mech 37:1-487. doi: 10.1146/annurev.fluid.37.061903.175743
  » https://doi.org/10.1146/annurev.fluid.37.061903.175743
• Peskin CS (1972) Flow patterns around heart valves: a numerical method. J Comput Phys 10(2):252-271. doi: 10.1016/0021-9991(72)90065-4
  » https://doi.org/10.1016/0021-9991(72)90065-4
• Peskin CS (1982) The fluid dynamics of heart valves: experimental, theoretical, and computational methods. Annu Rev Fluid Mech 14:1-442. doi: 10.1146/annurev.fl.14.010182.001315
  » https://doi.org/10.1146/annurev.fl.14.010182.001315
• Saiki EM, Biringen S (1996) Numerical simulation of a cylinder in uniform flow: application of a virtual boundary method. J Comput Phys 123(2):450-465. doi: 10.1006/jcph.1996.0036
  » https://doi.org/10.1006/jcph.1996.0036
• Seo JH, Mittal R (2011) A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations. J Comput Phys 230(19):7347-7363. doi: 10.1016/j.jcp.2011.06.003
  » https://doi.org/10.1016/j.jcp.2011.06.003
• Tseng YH, Ferziger JH (2003) A ghost-cell immersed boundary method for flow in complex geometry. J Comput Phys 192(2):593-623. doi: 10.1016/j.jcp.2003.07.024
  » https://doi.org/10.1016/j.jcp.2003.07.024
• Udaykumar HS, Mittal R, Rampunggoon P, Khanna A (2001) A sharp interface Cartesian grid method for simulating flows with complex moving boundaries. J Comput Phys 174(1):345-380. doi: 10.1006/jcph.2001.6916
  » https://doi.org/10.1006/jcph.2001.6916
• Udaykumar HS, Mittal R, Shyy W (1999) Computation of solid-liquid phase fronts in the sharp interface limit on fixed grids. J Comput Phys 153(2):535-574. doi: 10.1006/jcph.1999.6294
  » https://doi.org/10.1006/jcph.1999.6294
• Udaykumar HS, Shyy W, Rao MM (1996) Elafint: a mixed Eulerian-Lagrangian method for fluid flows with complex and moving boundaries. Int J Numer Meth Fluid 22(8):691-712. doi: 10.1002/(SICI)1097-0363(19960430)22:8<691::AID-FLD371>3.0.CO;2-U
  » https://doi.org/10.1002/(SICI)1097-0363(19960430)22:8<691::AID-FLD371>3.0.CO;2-U
• Ye T, Mittal R, Udaykumar HS, Shyy W (1999) An accurate Cartesian grid method for viscous incompressible flows with complex immersed boundaries. J Comput Phys 156(2):209-240. doi: 10.1006/jcph.1999.6356
  » https://doi.org/10.1006/jcph.1999.6356
• Zhu L, Peskin CS (2003) Interaction of two flapping filaments in a flowing soap film. Phys Fluid 15(7):128-136. doi: 10.1063/1.1582476
  » https://doi.org/10.1063/1.1582476

Publication Dates

• Publication in this collection
  Jul-Sep 2017

History

• Received
  24 Apr 2016
• Accepted
  27 Sept 2016