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 finiteelement 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 ghostcontrol 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.


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(Peskin , 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. 2001Udaykumar et al. , 1999Udaykumar et al. , 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 fi t 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.Th is is because the boundary may cut the Cartesian grids anywhere on the cells and create new arbitrary shape.It would make it more diffi cult to discretize the equations and calculation of fl uxes 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.Th e present approach is diff erent from the cut-cell method such that boundary cell shapes remain unchanged.

nuMerIcAl AlGorItHM
Th e 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 defi ned in the middle of each element.Th is local coordinate system divides each element into 4 SCVs.Each SCV is associated with an element node at its vertex.Th erefore, 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.
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: 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).

( ) ∑ ( ) ( ) ∑ ( ) ∑
Immersed Boundary Method Based on the Implementation of Conservation Equations Along the Boundary using Control-Volume Finite-Element Scheme

SuB-COnTROl VOlumES AnD nODE TYPES
Discretization of governing equations and calculation of fl uxes are done on SCVs; hence, the classifi cation of diff erent sub-control volumes and nodes is described here.Th ere 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 classifi ed into 3 types in this paper.Th e fi rst type of SCV is the "ordinary" or "fl uid" 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.Th is 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 defi ned when the boundary is placed in the SCVs of fl uid 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 fl ow fi eld 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 fl uid 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 fl ow variables to the ghost nodes.Th ese values are mostly assigned by a kind of interpolation scheme (Mittal and Iaccarino 2005;Majumdar et al. 2001;Ghias et al. 2004Ghias et al. , 2007)).In the present method, however, fl ow 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.

GOVERninG EQuATiOnS AnD DiSCRETiZATiOn
In Eq. 3 there is a detail analysis of how Navier-Stokes equations were discretized.Th e integral form of the incompressible Navier-Stokes equations for 2-D fl ow is given by (3) 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: 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 implementation 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: 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).
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: 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 In Eq. 7 the convection term is represented in stream wise direction and q = (u 2 + v 2 ) 1/2 .Expression for convected velocity u is obtained on integration points which encompass all relevant variables related to flow condition.The convecting velocity u ˆ on ip is obtained from Eq. 8 as follows: (4) (5) 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 ) ( ) (4) ̂ ̂ ( 6) ) (4) boundary.As mentioned before, in the present method fl ow, variables on the ghost nodes are determined by implementation of conservation laws and the boundary condition on ghost SCVs.Th erefore, the key-point in the present method is to clearly implement conservation laws on ghost SCVs along the boundaries.Th is 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.Th us, mass and momentum conservation equations on ordinary SCV 1 would be: 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 2 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 2 1 is the normal surface vector of sub-surface1; ds y 2 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.where: ν 1 is the volume of ordinary SCV 1; the superscript " o " 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 x 2 1 , 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 ν 2 this part.Substituting these parameters in Eq. 3 for SCV 2 it results in: On SSb, fl ux vectors E, F, G, and H are evaluated on ip b .Th ese fl ux vectors are evaluated for SS2 on ip 2 .Th e discrete form of Eq. 12 is given by where is the convecting velocity vector; → ds b is the normal surface vector in the outward direction and (5) (4) (4) ) (5) (4) (4) ) ) (4) (4) (1) ) ) (4) ds xb and ds yb are the components of → ds 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 ( → pq ˆ b .→ ds b ) = 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 and are evaluated using bilinear interpolaion defi ned 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: Nodes 2 and 1, as well as scν 2 and scν 1 are considered 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 2 is volume of the grey area; v 1 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 1 1 is normal surface vector of sub-surface 1; ds y 1 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.

where:
→ ds b 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 fi nal solution of this method, the 2 following equations are combined: → ds b h ) = 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: Immersed Boundary Method Based on the Implementation of Conservation Equations Along the Boundary using Control-Volume Finite-Element Scheme Here ν 2 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: Since the boundary b is solid, then momentum flux ( → pq ˆ b . → ds b h )u b would be 0; p b is evaluated using bilinear interpolation, and velocity gradient terms of and 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. 2001Udaykumar et al. , 1999Udaykumar et al. , 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).
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.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 ( )u b is also calculated from known values of the exact solution; p b is evaluated from 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.
To obtain the difference between the results of the exact solution and the present method we introduce the following equations for P and KE: where N is the number of fluid nodes in the domain; subscripts n and e indicate numerical and the exact solutions, respectively.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.
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 ( → pq ˆ b .→ ds b h ) in Eq. 13 is calculated by known values of the exact solution.In x-momemtum conservation equation introduced by Eq. 14, ( → pq ˆ b .→ ds b )u 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.
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.
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.
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, 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     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.
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.
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 solution 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.
Boundary conditions are implemented through conservation 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, and 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.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.

Figure. 1 .
Figure. 1.(a) Defi nitions 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 fi gure is control volume made up of 4 surrounding SCVs and the light grey area is SCV.

Figure 3 .
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 6 .
Figure6.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 2 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 2 is the normal surface vector on SSb and is equal to -→ ds b 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 ( → pq ˆ b .

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

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

Figure 10 .P
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.

Figure 15 .
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 .
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.

Figure 17 .
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.
Based on the Implementation of Conservation Equations Along the Boundary using Control-Volume Finite-Element Scheme

Figure 18 .
Figure18.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 20 .Figure 21 .Figure 22 .
Figure 20.Decay of temporal KE for 2 Reynolds numbers of 100 and 1,000.Comparison with the exact solution.Time

Figure 23 .
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.

Figure 24 .
Figure24.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.
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.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.