Treatment of solid objects in the Pencil Code using an immersed boundary method and overset grids
Two methods for solid body representation in flow simulations available in the Pencil Code are the immersed boundary method and overset grids. These methods are quite different in terms of computational cost, flexibility and numerical accuracy. We present here an investigation of the use of the different methods with the purpose of assessing their strengths and weaknesses. At present, the overset grid method in the Pencil Code can only be used for representing cylinders in the flow. For this task it surpasses the immersed boundary method in yielding highly accurate solutions at moderate computational costs. This is partly due to local grid stretching and a body-conformal grid, and partly due to the possibility of working with local time step restrictions on different grids. The immersed boundary method makes up the lack of computational efficiency with flexibility in regards to application to complex geometries, due to a recent extension of the method that allows our implementation of it to represent arbitrarily shaped objects in the flow.
encil Code; immersed boundary method; overset grids; compressible fluid dynamics; complex geometries
Fluid flow in a domain that contains an immersed solid object is a common case in computational fluid dynamics. Obstructions in the flow include (but are not limited to) cylinders, spheres, flat plates, rectangular or elliptical cylinders and spheroids, triangles, and complex geometries made out of a combination of these. Finding a method to represent such objects in the best possible way in simulations is not a trivial task, and the method used is often chosen specifically to the problem at hand.
For many generic shapes, such as cylinders, spheres, plates, etc., body-fitted structured meshes are commonly used to represent the object(s) in the flow. Body-fitted structured meshes conform to the object(s) in the flow domain and to the domain’s other physical boundaries (inlet, outlet, walls, etc.). Depending on the flow domain and object in the flow, this may require a deformation of the grid to conform to domain boundaries, in addition to the mapping procedures to map the grid in the flow domain to a simple computational domain. This may result in a grid with unnecessary local variations of the grid (e.g., a grid that is denser than necessary in certain areas of the domain) and time consuming grid generation (Versteeg and Malalasekera, 2007). A popular alternative to such meshes, particularly when the shape of the flow domain or objects in the flow domain is more complex, is unstructured meshes. Unstructured meshes provide the highest flexibility in grid adaptation to a particular flow geometry, and is a good alternative for complex geometries when finite-volume or finite-element formulations of the governing equations are used (Mavriplis, 1997). Disadvantages of such grids are much larger storage requirements than for structured grids (Pletcher et al., 2012), the need for intricate mesh generation techniques (Owen, 1998), and the difficulty in achieving high-order of accuracy.
By other choices of grid methods, the object(s) in the flow and the flow domain can be represented without the grid conforming to the object(s). Typically this is done by using a Cartesian grid, with a modification in either the flow equations or the grid cells in the immediate vicinity of the solid object(s). Popular methods of this type include immersed boundary methods (IBMs) (Peskin, 1972, 2002; Mittal and Iaccarino, 2005) and cut-cell methods (Quirk, 1994; Causon et al., 2000; Günther et al., 2011). These methods differ in that the IBM uses a Cartesian grid in the entire flow domain, while in cut-cell methods grid cells are ‘cut’ near the objects or domain boundaries that do not conform to the grids, and the flow equations are solved on the new, modified cells (Ingram et al., 2003). Due to this cell cutting, care must be taken such that the cut cells do not become too small, since this may be a potential source of numerical instabilities.
For the IBM, rather than modifying the grid cells near the solid object, the boundary conditions on the solid are imposed directly in the flow equations. This is done either by a continuous or a discrete forcing technique. In both cases a body-force, present due to non-conforming boundaries in the flow, is introduced in the Navier-Stokes equations. This is done either before discretization (continuous forcing) or after (direct forcing) (Mittal and Iaccarino, 2005). The latter of these is the preferred method for IBM used to represent rigid boundaries. One development of the discrete forcing method is to treat the immersed boundary as a sharp interface, and to impose the boundary conditions directly. This is done by using a combination of ghost-points inside the solid and mirror/image-points in the flow domain (set by interpolation) to reconstruct the solid (Tseng and Ferziger, 2003; Berthelsen and Faltinsen, 2008). An advantage with this approach is that the boundary conditions are handled without any added force in the flow equations, hence, the method can easily be implemented in an existing flow solver. The disadvantage is the accuracy reduction in the vicinity of the surface, although recent developments show that some of the challenges related to high-order accurate reconstructions of velocities near the surface can be overcome (Seo and Mittal, 2011; Xia et al., 2014). Further, finite-difference IBMs are, in general, not mass conserving. Finite-volume approaches with cut-cell methodology is appropriate if mass and momentum conservation needs to be guaranteed (Mittal and Iaccarino, 2005).
About ten years after the emergence of the first IBM, a new method was proposed to represent solids in the flow by using several grids overset one another (Steger et al., 1983; Benek et al., 1985). Such overset grid methods (often called Chimera methods) employ body-conformal grids at the boundaries of objects in the flow, but the grids do not extend to the physical boundaries of the domain. Rather, a background grid (typically uniform Cartesian) is used, and updated flow information of overlapping grid regions is communicated between grids at every time step. Note that special overset grid methods without background grids exist, like yin-yang grids where two identical component grids are used to cover a spherical surface, thus avoiding very small grid cells close to poles of the spherical geometry (Kageyama and Sato, 2004).
The flow domain resolved with overset grids may contain a single grid overlapping another, or several grids overlapping necessitating a priority of communication and computation of solutions of the different grids (Steger and Benek, 1987; Chesshire and Henshaw, 1990). For complex configurations, this may require extensive preprocessing for fixed objects (Suhs et al., 2002) or intricate grid handling at run time for moving bodies (Noack, 2005). Overset grid methods are, in general, not mass conserving, since interpolation is necessary between grids overset one another (although exceptions do exist, for finite-volume implementations of overset grids, see Pärt-Enander and Sjögreen (1994) and Zang and Street (1995)). The interpolation is done from donor-points on one grid to fringe-points on another. Many different interpolation procedures have been explored for this purpose, and several studies have found that using high-order interpolation between grids is beneficial in regards to the overall accuracy and stability of flow computations (Sherer and Scott, 2005; Chicheportiche and Gloerfelt, 2012; Völkner et al., 2017).
In the high-order compressible flow solver known as the Pencil Code (The Pencil Code, 2018), solid objects in the flow can be represented by different schemes. This makes it possible to compare different surface representations not only for the same flow problems, but for simulations where the same finite-difference discretization, time integration, communication procedures, etc, are used. The purpose of this paper is to perform such a comparison for solids represented by a ghost-points IBM and overset grids. The performance of these surface representations is assessed in terms of computational cost and accuracy for a common benchmarking case. The flow case used is the frequently appearing fluid mechanics problem of flow past a circular cylinder. Further, we wish to shed light on an advantage of the IBM implementation in the Pencil Code by simulating flow past a complex geometry. The complex geometry used as an example case is a combination of a semi-circular and a semi-elliptical cylinder.
The structure of the paper is the following: In Section 2 the governing flow equations and the two methods for solid object representation are described, with details on their implementation in the Pencil Code included. Performance of the different methods for the flow past a cylinder in both the steady regime and the unsteady vortex shedding regime is compared in Section 3. Following this, we elaborate on the extension of the IBM to complex geometries in Section 4, before concluding remarks are made in Section 5.
2.1 Governing equations
The governing equations of the flow are the continuity equation:
and the momentum equation:
where , , , and are the density, time, velocity vector, pressure, and dynamic viscosity (, with kinematic viscosity ), respectively , and
is the material derivative operator. The compressible rate of strain tensor is given by:
where is the identity matrix. The pressure is computed by the isothermal ideal gas law, where is the speed of sound. With a constant speed of sound (for the isothermal case) and a constant kinematic viscosity, the momentum equation (Eq. (2)) can be re-written to:
which is the form solved in the computations performed in this study.
2.2 Numerical methods
Many different types of domain and enforcements of boundary conditions are available in the Pencil Code. For simplicity we consider a domain with a uniform mean flow, using Navier-Stokes characteristic boundary conditions (NSCBC) on both the inlet and the outlet of the flow domain, and periodic boundary conditions in all other directions. The NSCBC is a formulation developed by Poinsot and Lele (1992) that makes use of one-dimensional characteristic wave relations to allow acoustic waves to pass through the boundaries.
We place an object in the flow domain by representing it with one of the two available methods in the Pencil Code. The grids used in each of the methods are quite different. An illustration of the grids for IBM and overset grid representation of a circular solid can be seen in Fig. 1. For both methods of solid body representation, we use boundary conditions of no-slip and impermeability for velocity components, and zero gradient for density in the direction normal to the surface, on the solid’s surface. The latter condition can be derived from the ideal gas law and the boundary layer approximation for pressure normal to the boundary ( (White, 2006)). In the remainder of this section, details of how the different boundary representations are implemented in the Pencil Code are given.
2.2.1 Ghost-zone immersed boundary method
In the illustration of a circular object in a flow domain represented by an IBM (Fig. 1(a)) the intersections of solid grid lines represent fluid-points, where the governing equations are solved, while the intersection points of dashed grid lines are grid points inside an immersed solid (solid-points). At the solid-points the governing equations are not solved. Rather, some points are used as ghost-points for the fluid solver and some are unused points. As mentioned, the boundary conditions of a solid object may be imposed directly on the flow variables by the ghost-point immersed boundary method. The IBM in the Pencil Code is such a ghost-point method. An uncommon feature of the IBM implementation in the Pencil Code is that a several points deep ghost-zone is used, rather than a single strip of ghost-point inside the solid object. This ensures that the sixth-order finite-difference stencils can be used without any modifications in the vicinity of a solid object. The overhead related to computation of two additional layers of ghost-points is negligible when compared to the computational cost of the fluid solver itself.
As sixth-order central differencing is used, three points on each side of a grid point are necessary to update the solution. This is illustrated in Fig. 2(a), where stencils of fluid-points and will include grid points within the solid object to update the horizontal gradients of the velocity components and density. Point will need information from , and , will need information from , and , etc. (in addition to information from fluid-points to the left). The points in the ghost-zone, and are set using corresponding mirror-points and , respectively. With no-slip and impermeability for velocity and zero gradient for density, the relationship between a ghost point and mirror-point is simply:
Note, that a second-order accurate method to set the Neumann boundary condition has recently been implemented, but will not be described here. Details, and testing of different boundary conditions, can be found in Luo et al. (2016).
In general, the mirror-points do not coincide with grid points or grid lines, and need therefore to be interpolated from surrounding points. This is done by bi-linear Lagrangian interpolation in 2D and tri-linear Lagrangian interpolation in 3D. For mirror-points close to the surface, one or more of the surrounding grid points may be inside the solid object (see Fig. 2(b)). Rather than using four surrounding fluid-points, these mirror-points are set by interpolation using the boundary intercept point (that is, the point where the surface normal through the ghost point intercepts the solid’s boundary) and the point where the surface normal intercepts the first grid line outside the solid. Data is first interpolated to the interception point of the surface normal with the first grid line from neighboring fluid-points ( in Fig. 2(b)) by linear Lagrangian interpolation. The velocity components normal to the surface is expected to scale as when approaching the surface, where is the distance from the boundary. To interpolate velocity in a mirror-point near the surface the velocity is decomposed into cylindrical components, and the radial component is computed by:
where is the radial velocity at the grid line interception point, and and are the distances from the mirror-point and grid line interception point to the boundary interception point, respectively. Remaining velocity components are obtained by linear interpolation. No special handling is used for density.
Special handling is used for fluid-points very close to the surface of an object. This is done to avoid spurious errors due to de-localization dependencies in the finite-difference stencils, a detrimental effect that occurs when flow variables quite far from a grid point is indirectly used in the update of said grid point. To see this, consider in Fig. 2(a) that the horizontal velocity component of grid points surrounding will affect , since is set by . Rather than computing flow variables in a fluid-point close to the surface in the usual way, by using the finite-difference stencils, they are set directly by interpolation, as seen in Fig. 2(c). The interpolation procedure for these grid points is the same as for mirror-points very close to the boundary, as described above. Note that this type of special handling is only possible for variables with a Dirichlet boundary condition on the surface.
An alternative to setting mirror-point positions using surface normals is to use mirror-points along grid lines. This simplifies interpolation (making all interpolation one-dimensional, along grid lines) and has proven promising in reducing the errors due to de-localization dependencies. In such an approach, a ghost point can have several values for each flow variable that has a Dirichlet boundary condition, one used in horizontal and one in vertical finite-difference stencils (as would be the case for and in Fig. 2(a)). Flow variables with Neumann boundary conditions are set from mirror-points along surface normals as in the method described above. In this study, we will stick to the more mature method of using mirror-points from surface normals. Details on the alternative grid-line ghost-zone IBM can be found in Aarnes et al. (2018b).
2.2.2 Local-time restricted overset grids
Unlike the representation of solid bodies with most methods (IBM, body-fitted structured or unstructured grids), codes using overset grids require splitting of the flow solver, as one solver is needed for each grid. We limit this study to a single grid on top of a background grid, for a more general discussion see Chesshire and Henshaw (1990) or Meakin (1995). Yin-yang grids are not considered (although a yin-yang grid implementation exists in the Pencil Code).
For a flow with a solid object represented using a body-confined grid over a Cartesian background grid, the governing equations are, in principle, solved for two different flow domains, one with and one without a solid object present in the flow. At least one boundary in each domain is set by interpolating flow variables from another grid, and in this way the presence of the solid affects the flow on all grids. Admittedly, this makes overset grids somewhat more unwieldy than IBMs.
To consider this more systematically, let us regard a fluid time step as split into four parts: (1) solution of the governing equations on the background grid, (2) communication of data from the background grid to the body-fitted grid, (3) solution of the governing equations on the body-fitted grid, (4) communication of data from the body-fitted grid to the background grid. The solution step on the body-fitted grid requires implementation of a Navier-Stokes solver applicable to the type of grid that is used to resolve the bluff body’s boundary. For our case of a cylinder in a cross flow, a Navier-Stokes solver applicable to cylindrical coordinates is necessary.
In the illustration of the mesh used in the (cylindrical) overset grid method (Fig. 1(b)) there are no grid points inside the circular object. Strictly speaking, the background Cartesian grid is present in the entire domain (also inside the limits of the curvilinear mesh and inside the solid object), but only few points inside the curvilinear mesh are used, and not a single Cartesian grid point inside the solid is (or should ever be) used.
The points of the background grid that are in use inside the curvilinear mesh are mostly fringe-points. Fringe-points are fluid-points set by interpolation from nearby donor-points on the overlapping grid, rather than computed using discretization of the governing equations. Donor-points are computed in the same way as an ordinary fluid-point, unless the interpolation is implicit, meaning that a fringe-point may be used as a donor-point (Chesshire and Henshaw, 1990), which is not the case here. A third class of points found on overset grids are hole-points. These are unused grid points, typically found outside of the fringe of a grid (like Cartesian grid points far inside the curvilinear grid). Fringe-, donor-, and hole-points on an overset grid generated to represent a circular object are seen in Fig. 3.
During the first stage of the inter-grid communication, data is sent to the outermost grid points of the body-fitted grid (Fig. 3(a)). At the second stage, data is sent back to Cartesian fringe-points (Fig. 3(b)). The fringe-points on the Cartesian grid are identified during pre-processing, where points within a distance from the solid, illustrated by red circular curves in Fig. 3(b), are selected. For moving objects, fringe-point locations on the background grid must be re-calculated every time step. A several layers thick zone of fringe-points is used on both grids to enable the use of the sixth-order centered stencils at the outer edges of each grid, equivalent to the use of a ghost-zone for the IBM.
In the illustration in Fig. 3 each fringe-point is surrounded by four donor-points, as necessary in bi-linear Lagrangian interpolation. If higher order interpolation is desired the amount of donor points for each fringe points must be increased accordingly. Such an increase is straightforward for overset grids, unlike for immersed boundary methods where more intricate interpolation stencils are needed for high-order interpolation to avoid using grid points that are inside the bluff body. Note, however, that a straightforward extension from second to third order interpolation (or higher) does not guarantee a better solution. This is due to possible overshoots in the interpolation polynomials. High-order Lagrangian interpolation and quadratic splines are implemented for overset grid interpolation in the Pencil Code. For simplicity, we will restrict ourselves to bi-linear Lagrangian interpolation here.
At the solid-fluid interface, we use summation-by-parts (SBP) finite-difference operators to enhance stability of the solution. This means modifying the finite-difference stencils in the nine grid points closest to the surface (including the surface point) along each radial grid line, to asymmetric stencils (one-sided at the surface). The order of accuracy for the SBP-operators are third-order for a sixth-order finite-difference method. Details on these operators can be found in Strand (1994) (first derivatives) and Mattsson and Nordström (2004) (second derivatives).
A peculiar feature of the overset grid implementation in the Pencil Code is how the restrictions on the time step are handled. The advective and diffusive time step restrictions are and , respectively, where is the time step, the smallest grid spacing in any direction, and and are the diffusive and advective Courant numbers, respectively. For a weakly compressible flow we typically require a very short time step, increasingly so if grid stretching is used in order to have a fine grid in the vicinity of the solid object. However, when overset grids are used these restrictions are no longer global restrictions on the time step, but local. Hence, by performing several time steps on the body-fitted grid for each time step on the background grid, the efficiency of the code may be greatly improved. In particular, this allows for a very fine resolution close to the surface (on the body-fitted grid) without using small time steps for flow far from the surface (on the background grid).
For all the overset grid computations in the present study the diameter of the cylindrical grid is three times that of the solid cylinder it is fitted to. For a consideration of the extent of the domain covered by the body-fitted grid the reader is referred to Aarnes et al. (2018a).
2.2.3 A note on dissipation
The centered finite-difference schemes used for discretization of the governing equations are non-dissipative. This can cause problems due to the potential growth of high-frequency modes, leading to numerical instability.
To some extent, the summation-by-parts boundary conditions suppress such instabilities that are due to boundary conditions when overset grids are used, but these boundary stencils are not sufficient to suppress all oscillations in the solution when grid stretching is used on the curvilinear grid. Such oscillations are most prominent in the density field. The detrimental effect of the high-frequency modes increases as the grid spacing decreases, and may lead to diverging solutions as the grid is refined. To suppress the high-frequency modes, a high-order low-pass filter is used on the curvilinear part of the overset grid. The filter is a 10th order Padé filter, with boundary stencils of 8th and 6th order. On the interior of the domain, the filter is given by:
where and are components of the filtered and unfiltered solution vectors, respectively, is a free parameter () and are fixed parameters dependent only on (details in Visbal and Gaitonde, 1999). Boundary stencils can be found in Gaitonde and Visbal (2000). The Padé filter is implicit, and requires us to solve a tri-diagonal linear system at each grid point, in the radial direction, and a cyclic tridiagonal system in the direction tangential to the surface. The free parameter is set to 0.1. With such a small value for , filtering the solution once per Cartesian time step is found sufficient to get a stable and accurate solution.
Alternatively, a sixth-order hyper-diffusion operator, which is already implemented in the Pencil-Code (see e.g. Haugen and Brandenburg, 2006), could have been used to filter the solution. The benefits of this approach is that the hyper-diffusion operator is explicit and fast, and does not require extra communication between processors. It is, however, expected to be less sharp than the 10th order Padé filter, as Padé filters are known to outperform explicit filtering schemes (Visbal and Gaitonde, 1999, 2002).
When IBM is used rather than overset grids, some dissipation can be turned on by using fifth-order upwinding for the advection operators of the density rather than central-difference stencils (details in Dobler et al., 2006). The problem of oscillations in the density field is, however, much less prominent when a uniform Cartesian mesh is used, so while Padé filtering is on by default when overset grids are used, a dissipative solution by upwinding is optional for other simulations.
3 Simple geometry
In previous studies, the order of accuracy of the solid object representations described above has been assessed for steady flow computation. Using a slightly modified handling of Neumann boundary conditions, Luo et al. (2016) showed that, regardless of boundary condition, the IBM implementation in the Pencil Code is second-order accurate in the vicinity of a resolved circular boundary. For the same geometry, using overset grids, Aarnes et al. (2018a) showed that the order of accuracy in the vicinity of the solid differed for different flow variables. The radial velocity component was computed with order of accuracy between third and fifth-order, while the accuracy of the tangential velocity and density was between second and third-order, when second-order Lagrangian interpolation was used for communication between grids. Both the mentioned studies also demonstrated that characteristic flow parameters, like drag, lift and shedding frequency (for unsteady flow), could be reproduced to good agreement with previous studies, with the respective boundary representation in use.
We will not repeat an assessment of accuracy for steady flow computations here. Rather, we investigate the boundary representations by a direct comparison for the case of flow past a circular cylinder in different shedding regimes. Two-dimensional flow past a circular cylinder in the vortex shedding regime is a classical benchmarking case for fluid dynamic simulations with a solid object present in the flow domain. We simulate such a flow with Reynolds number 100, a Reynolds number where the von Kármán vortex street can be observed in the cylinder’s wake. The Reynolds number is defined as , where is the incoming flow velocity and is the diameter of the cylinder obstructing the flow. In addition, we test each boundary representation for a steady flow () and an unsteady flow with more chaotic tendencies (). Note that three-dimensional effects in the latter flow are suppressed as we restrict ourselves to a two-dimensional domain. A rectangular domain with domain size is used, with a body-fitted grid with diameter used in the overset grid simulations. For overset grids, grid stretching is used in the radial direction to obtain approximately quadratic cells at the cylinder surface and in the region where the interpolation between the grids is performed. In this region the cells of both grids are similar in size. The inflow Mach number () is set to 0.1, and the Reynolds number is varied by adjusting the value of the kinematic viscosity. The vorticity component normal to the -plane for the three different Reynolds number flows with mean flow in the -direction can be seen in Fig. 4.
Consider Fig. 5, depicting normalized deviations of mean drag coefficient and root-mean-square lift coefficient computed at different resolutions with the two methods for . For both cases, the results from the finest grid is used for normalization. Strouhal number (dimensionless shedding frequency) is not included in the figure, since it is barely affected by the grid spacing and is therefore not a good measure of grid independence. The number of grid points per diameter, on the Cartesian grid, is given on the horizontal axis. Note that this might be somewhat misleading, as the overset grid uses two grids to cover the flow domain. For the flow domain and grid sizes used in this study, an overset grid simulation uses 10% more grid points than a corresponding IBM simulation when of the two simulations is the same. No matter this difference, it is clear that the overset grid method greatly outperforms the IBM with respect to the necessary grid required to reach grid independency under the conditions of these simulations. Using a background grid with the deviation from results on a grid is less than 0.16% for the overset grid. The results from the IBM calculations converge much slower. To reach a comparable level of grid independence to that of the overset grid with , a grid using IBM requires . Such a fine grid yields a deviation of less than in drag and in lift, from the results at the finest grid level (). If these grids ( for overset grids, for IBM) are deemed sufficiently accurate resolutions for grid independent solutions for the different solid object representations, the IBM requires 18.1 times as many grid points as the overset grid method on the two-dimensional domain used in these simulations. In these simulations, the advective restriction is more strict than the viscous restriction, hence, the time step is proportional to the grid spacing. This means that there is a factor 4.7 difference in time step between the and , that has an additional large impact on the (in)efficiency of the IBM as compared to the overset grids.
For practical application, it may perhaps be excessive to require deviation for results to be deemed grid independent. With a resolution in the simulations performed with IBM, there is less than 0.5% deviation in drag and less than 1% deviation in lift. Choosing such a resolution will, in many cases, be an acceptable trade-off between accuracy and efficiency. This reduces the difference between the overset grid and IBM somewhat, although it warrants the use of a somewhat coarser grid for overset grid computations as well. Note that in computing drag and lift forces on the cylinder represented by the IBM, so-called force-points are used. The force-points are distributed uniformly around the cylinder, and viscous and pressure forces are approximated at these points, using data from surrounding grid points. Some of the oscillations in the computed mean drag and root-mean-square lift coefficients seen in Fig. 5 may be due to a change in the position of force-points when the grid is refined.
A relevant consideration when the different costs associated with overset grids and IBM are compared is the computational cost of interpolation in the two different methods. With equally spaced Cartesian grids, more fringe-points are interpolated with overset grids than mirror-points interpolated with the ghost-point IBM method. This is due to overset grids having two zones of interpolation (one for interpolation from Cartesian to cylindrical, and one for interpolation back to Cartesian), a larger circumference of the interpolation regions (interpolation farther from the cylinder) and the need for a deeper ghost-zone on the Cartesian grid in overset grids since no special handling for fluid-points close to the fringe-point region is used. For the grid, the total number of fringe-points for the overset grids is approximately 1600 (% of total number of grid points). Note that the fraction of fringe-points to grid points decreases as the grid spacing decreases ( grid points are fringe-points when ). The number of mirror-points in an IBM simulation with the same grid spacing is approximately 200. A fringe-point is updated only once every Runge-Kutta time step, while mirror-points are updated every sub-time step. Hence, approximately 2.5 times as much interpolation is performed when the solid is represented by overset grids rather than by IBM, if the same grid spacing is used in the different solid object representations. As much finer grids are required with IBM than with overset grids, the advantage of a smaller interpolation cost with IBM is lost.
|Kim et al. (2001)||1.33||0.165|
|Haugen and Kragset (2010)||1.328||0.166|
|Park et al. (1998)||,||1.33||0.165|
|Shi et al. (2004)||1.318||0.164|
|Stålberg et al. (2006)||,||1.32|
|Li et al. (2009)||1.336||0.164|
|Posdziech and Grundmann (2007)||1.350||0.167|
|Posdziech and Grundmann (2007)||1.312||0.163|
|Qu et al. (2013)||1.326||0.2191||0.166|
|Qu et al. (2013)||1.310||0.2151||0.165|
|Present, overset grid.||1.347||0.234||0.166|
To verify that the flow is computed accurately, the resolutions from the discussion above ( for overset grids, for IBM) are used in a simulation on a large domain for each of the solid body representations. The domain size is set to . The resulting mean drag, root-mean-square lift, and Strouhal frequency are compared to results reported from other studies, in Tab. 1. Domain sizes and types are listed in the table, along with the most relevant flow coefficients. Some of the listed values for root-mean-square lift coefficients are scaled values, as only amplitude of the lift coefficient was reported from these particular studies. A scaling factor of has been used (since the lift coefficient is a smooth sinusoidal-like function with zero mean value), and the scaled results are marked with a superscript (*). The studies in Tab. 1 use a wide range of numerical methods to compute the flow, including finite-volume, finite-difference, finite-element, spectral element and lattice-Boltzmann methods. The top three studies in Tab. 1 use immersed boundary methods to represent the solid cylinder, while the remaining studies (present IBM simulations excluded) use body-fitted methods. Only Haugen and Kragset (2010) and Li et al. (2009) simulate compressible flows (where the former of these uses the Pencil Code with the IBM described here, but with different domain size and resolution). Table 1 includes two results from each of the studies by Posdziech and Grundmann (2007) and Qu et al. (2013), to include results from both comparable domain sizes to the present study, and highly accurate results from very large domains. The present results, both those computed with the IBM and the results from overset grid simulations, agree well with the results found in the literature.
Grid independence results for and are depicted in Fig. 6. The results are similar to those obtained for : there is a much more rapid convergence to grid independent solutions with overset grids than with IBM. With overset grids, a background grid with yields less than 0.2% deviation in the drag and lift coefficients, at . This means using overset grids . With IBM, a grid with is necessary to get comparable grid independence (deviation less than 0.25% in drag and lift coefficients). That means using a grid for this specific domain size, and a factor four larger time step than on the background grid in the overset grid simulation. For the steady flow, the lift coefficient is not defined, and for this reason we include only the drag coefficient in the grid independence comparison. The results in Fig. 6(a) show that is sufficient to obtain drag with less than 0.4% deviation from the finest grid result for the overset grid computation. With IBM, is needed to get the deviation down to the same level.
From these tests it is clear that in representing the simple geometry of a circular cylinder with the Pencil Code, the method of overset grids is far superior to the IBM in terms of efficiency and accuracy. Not only is far less grid points required (and, consequently a larger time step allowed) to reach a grid independent solution, there is also far less variation in the solution before grid independence is reached. With the IBM we may have to accept a deviation of, say, 1.0% in mean drag and root-mean-square lift coefficients from one grid to a finer one. In the case of overset grids, on the other hand, a deviation one order of magnitude smaller than this can be achieved at reasonable computational costs. As mentioned, however, some of the variation seen in the IBM results may be attributed to the way the coefficients themselves, and not the flow, are computed. The positioning of force-points where drag and lift are computed is affected by the choice of resolution, but has no influence on the solution of the flow equations.
That being said, we should now address the limitations of the overset grid method. In short, the problem of adaptability to different geometries is a major drawback of this method. Even an extension from a 3D cylinder to a sphere would require a completely new grid handling, with updates needed all the way down to the level of finite-differences in the code. For more complex shapes, where an analytic transformation from Cartesian space to the fitted grid coordinates is not available, this will become increasingly difficult, if not impossible. It is in this respect that the full potential of the IBM can be achieved. The simple handling of boundaries and lack of any modification needed in the treatment of the governing equations, make the ghost-zone IBM ideal for complex geometries of all kinds. How this is done in the Pencil Code is the topic of the remainder of this paper.
4 Complex geometries
One of the difficulties in extension to irregular geometries of a ghost-cell immersed boundary method’s lies in how to track the boundaries correctly. To the best of our knowledge, two ways to overcome this exist: the unstructured triangle surface mesh (Gilmanov et al., 2003; Mittal et al., 2008; Nagendra et al., 2014) and the combination with level-set signed distance functions (Liu and Hu, 2014; Uddin et al., 2014). The first method can be used to represent arbitrary geometries and has gained its popularity in biological fluid mechanics. For example, interactions between a very complex body, such as a bluegill sunfish pectoral fin or a false vocal fold, and its surrounding flows have been studied in a two-way coupled manner via the first method (Zheng et al., 2009; Dong et al., 2010). This method, with an unstructured surface mesh for the complex boundary, was introduced and implemented in the Pencil Code by Luo et al. (2016). Arbitrary two-dimensional immersed boundaries are represented by many small line-segments. Each line-segment is identified by two vertices as shown in Fig. 7. The general procedure is still the same as illustrated in Section 2.2.1, other than some special handling of fluid-points and mirror-points around the solid object.
The first difference lies in the identification of a given grid point to be a fluid-point or a solid-point. For a circular object, the distance from a given grid point to the center of the circle is calculated, and compared with the radius of the circular object to identify if this grid point is a solid grid point or not, i.e. if it is inside the object or not. This becomes an ineffective method for a complex geometry, where no single radius can be found for the object. In this case, for a given grid point, the closest surface element is detected first. Secondly, the dot product between the closest line-segment’s normal vector and the direction vector pointing from the centroid of the closest facet to the given grid point is calculated. The sign of the dot product determines the identification of the grid point. Generally, a negative result indicates a solid-point for the convex boundary shown in Fig. 7(a). The treatment of some special cases that may occur during this process is described in Luo et al. (2017).
After the identification of solid/fluid-points, three layers of ghost points are assigned to construct a six-order central finite-difference stencil as shown in Fig. 7(b). Following this, a corresponding boundary intercept point is determined for each one of them. The method to detect the boundary intercept points is different from that of the simple circular geometry. First, the vertex closest to a given ghost point is determined. Then, the set of line/surface elements sharing that vertex can be identified and a search is carried out among these elements to find the boundary intercept point (which should lie within the line/surface elements) as shown in Fig. 7(b). While conceptually simple, the implementation can be very complicated and special attention is needed to find the correct intercept point. Here, we adopt a method based on the robust procedure proposed in Mittal et al. (2008). For details, see Luo et al. (2017).
Once boundary intercept points are determined for every ghost point, a corresponding mirror-point can be obtained. The mirror-points are set either by symmetry over the solid’s line element (corresponding to the way mirror-points are set with a simple geometry, see Fig. 2), or at a constant distance away from the boundary intercept point. This distance, in Fig. 7(b), is typically set to for 2D geometries and for 3D-geometries, to ensure that every mirror-point is surrounded by fluid-points only. The same interpolation procedure as for the simple circular geometry (i.e., the bi-linear interpolation method) can be adopted for the calculations of the parameters of the mirror-points. An optional way is the inverse distance weight interpolation method (Chaudhuri et al., 2011). Finally, the flow variables at the ghost-points can be calculated with the aid of the mirror-points and the given boundary conditions at the boundary intercept point by linear interpolation. Three types of boundary conditions, the Dirichlet, Neumann or Robin boundary condition have been implemented, similarly as for the simple circular object discussed in Section 2.2.1. For more details and test of the boundary conditions, the reader is referred to Luo et al. (2016).
This new method can be straightforwardly extended to complex three-dimensional geometries, where triangular surface elements can be adopted to represent the surfaces. Details about the method and its implementation can be found in Luo et al. (2017) and will not be repeated here. In Luo et al. (2017), a spatial convergence test indicates that only the bi-linear interpolation procedure can obtain a local second-order accuracy. Systematic validations have also been conducted through calculations of flow past an elliptical cylinder, square cylinder, semi-cylinder, as well as a NACA0012 airfoil. Quantitative comparisons with reported results in literature show that the present method can accurately reproduce the main features of the fluid flow past solid objects with complex geometry, quantified by coefficients such as drag and lift coefficients, Nusselt number and Strouhal number.
To demonstrate the IBM capabilities for a two-dimensional flow, we have simulated flow past geometries constructed by combining a semi-circle and semi-elliptical cylinder. The geometry is seen in Fig. 8. The radius of the semi-circle () and the major axis of the semi-ellipse () can be varied to construct different geometries. Three cases, with 2.0, 1.0 (circle), and 0.5, respectively, are considered. For each case, 360 line segments are used to resolve the immersed geometry. Other parameters related to the computational domain are kept consistent with the case in the grid refinement part of Section 3, except that the solid body is placed at a distance from the inlet, in the streamwise direction, rather than in the center of the flow domain ( from the inlet).
The vorticity component normal to the view plane for flow past the three different geometries can be seen in Fig. 9. It can be seen that as the length of the semi-ellipse is increased, the length of the bound vortex increases accordingly. This results in different patterns of von Kármán vortex streets for each of the three cases. The corresponding mean drag coefficient, root-mean-square lift coefficient, and Strouhal frequency number are listed in Tab. 2. It is obvious that even though the immersed boundary of the third geometry is the longest, the drag force it experiences is the least. This is perhaps not surprising, as the shape of the object is closer to a streamlined body for the largest value of , a shape that is known for low drag.
5 Concluding remarks
In this study we have described and compared the two solid body representations available in the high-order finite-difference code known as the Pencil Code. The two methods, the immersed boundary method and overset grids, are fundamentally different in many aspects. These differences can be summed up as:
The ghost-point IBM can be implemented straightforwardly in an existing flow solver, by extending the code without requiring major modifications to the existing solver. Using overset grids requires a more generalized flow solver, able to handle all grids that are overset one another (Cartesian, cylindrical, etc.). This may require a modification of the flow solver itself, when overset grids are first implemented in an existing fluid dynamics code.
Neither the IBM nor overset grids are mass conserving, as they both rely on interpolation of flow quantities to either mirror-points or non-conforming grid points, respectively. As the interpolation is moved away from the solid surface when overset grids are used, the accuracy loss following from interpolation is expected to have a reduced impact on flow properties directly related to the solid object (boundary layer properties, etc.).
For a circular cylinder, using the overset grid method in the Pencil Code is far superior to using IBM. Reaching grid independent solutions with IBM required 4.7 and 4 times as many grid points in each direction for and 400, respectively, as compared to the background grid used in the overset grids method. In total, this meant using 18.1 and 14.5 times as many grid points in our two-dimensional simulations with IBM as compared to that with overset grids at these Reynolds numbers. In addition there comes a much stricter limitation on the time step for the fine grid used in the IBM. Such a limitation is only present on the curvilinear grid in the overset grid method, while a 4-5 times as large time step can be used on the coarse background grid.
The IBM is highly flexible. The implementation in the Pencil Code can handle complex geometries, that is, surfaces where an analytic surface representation is not available. This opens up a large area of research, that cannot be studied with overset grids.
The development of both IBM and overset grids is far from over, and the evolution of these methods in the Pencil Code is destined to continue as long as researchers use the code and implement their own improvements into this open-source software. Perhaps, in time, overset grids can become more flexible, and the high-order accuracy achieved for the simple geometry can be available for more complex shapes (e.g., by using many grids overset one another). Alternatively (or, perhaps, in addition), the accuracy of the IBM implementation may be improved through implementation of stable, high-order interpolation of mirror-points in the vicinity of the solid object. We hope that with this paper, more researchers will be attracted to use the Pencil Code for simulations of flow past solid objects. Only in this way can the advancements of the solid object representations in the software continue, in the spirit of open-source software development.
We would like to acknowledge that this research is funded by the Research Council of Norway (Norges Forskingsråd) under the FRINATEK Grant [grant number 231444] and by the National Natural Science Foundation of China [No. 51576176]. The research is supported in part by The GrateCFD project [grant 267957/E20], which is funded by: LOGE AB, Statkraft Varme AS, EGE Oslo, Vattenfall AB, Hitachi Zosen Inova AG and Returkraft AS together with the Research Council of Norway through the ENERGIX program. Computational resources were provided by UNINETT Sigma2 AS [project numbers NN9405K, NN2649K].
Appendix A Sample cases
To get started with simulations of flow past a cylindrical geometry using the Pencil Code, sample cases that are available with the download of the code (The Pencil Code, 2018) may be useful. From the pencil-code directory, the sample cases can be found in:
The postfix ogrid denotes the overset grid sample case. To compile and run a sample case, use the commands
Both sample cases are simulations of a particle-laden flow past a cylinder at , in which particles may impact and deposit on the cylinder surface. For documentation on the handling of particles and particle deposition with the cylinder represented by IBM and overset grids, the reader is referred to Haugen and Kragset (2010) and Aarnes et al. (2018a), respectively.
- Aarnes et al. (2018a) Aarnes, J.R., Andersson, H.I., and Haugen, N.E.L., 2018a. High-order overset grid method for detecting particle impaction on a cylinder in a cross flow. (Manuscript submitted for publication).
- Aarnes et al. (2018b) Aarnes, J.R., Andersson, H.I., and Haugen, N.E.L., 2018b. Numerical investigation of free-stream turbulence effects on the transition-in-wake state of flow past a circular cylinder. Journal of Turbulence, 19 (3), 252–273.
- Benek et al. (1985) Benek, J.A., Buningt, P.G., and Steger, J.L., 1985. A 3-d chimera grid embedding technique. In: AIAA 7th Computational Fluid Dynamics Conference, AIAA Paper No. 1523.
- Berthelsen and Faltinsen (2008) Berthelsen, P.A. and Faltinsen, O.M., 2008. A local directional ghost cell approach for incompressible viscous flow problems with irregular boundaries. Journal of Computational Physics, 227 (9), 4354–4397.
- Causon et al. (2000) Causon, D.M., et al., 2000. Calculation of shallow water flows using a cartesian cut cell approach. Advances in Water Resources, 23 (5), 545–562.
- Chaudhuri et al. (2011) Chaudhuri, A., Hadjadj, A., and Chinnayya, A., 2011. On the use of immersed boundary methods for shock/obstacle interactions. Journal of Computational Physics, 230 (5), 1731–1748.
- Chesshire and Henshaw (1990) Chesshire, G. and Henshaw, W.D., 1990. Composite overlapping meshes for the solution of partial differential equations. Journal of Computational Physics, 90 (1), 1–64.
- Chicheportiche and Gloerfelt (2012) Chicheportiche, J. and Gloerfelt, X., 2012. Study of interpolation methods for high-accuracy computations on overlapping grids. Computers & Fluids, 68, 112–133.
- Dobler et al. (2006) Dobler, W., Stix, M., and Brandenburg, A., 2006. Magnetic field generation in fully convective rotating spheres. The Astrophysical Journal, 638 (1), 336–347.
- Dong et al. (2010) Dong, H., et al., 2010. Computational modelling and analysis of the hydrodynamics of a highly deformable fish pectoral fin. Journal of Fluid Mechanics, 645, 345–373.
- Gaitonde and Visbal (2000) Gaitonde, D.V. and Visbal, M.R., 2000. Padé-type higher-order boundary filters for the Navier-Stokes equations. AIAA Journal, 38 (11), 2103–2112.
- Gilmanov et al. (2003) Gilmanov, A., Sotiropoulos, F., and Balaras, E., 2003. A general reconstruction algorithm for simulating flows with complex 3D immersed boundaries on cartesian grids. Journal of Computational Physics, 191 (2), 660–669.
- Günther et al. (2011) Günther, C., et al., 2011. A cartesian cut-cell method for sharp moving boundaries. In: 20th AIAA Computational Fluid Dynamics Conference. 3387.
- Haugen and Brandenburg (2006) Haugen, N.E.L. and Brandenburg, A., 2006. Hydrodynamic and hydromagnetic energy spectra from large eddy simulations. Physics of Fluids, 18 (7), 075106.
- Haugen and Kragset (2010) Haugen, N.E.L. and Kragset, S., 2010. Particle impaction on a cylinder in a crossflow as function of Stokes and Reynolds numbers. Journal of Fluid Mechanics, 661, 239–261.
- Ingram et al. (2003) Ingram, D.M., Causon, D.M., and Mingham, C.G., 2003. Developments in cartesian cut cell methods. Mathematics and Computers in Simulation, 61 (3-6), 561–572.
- Kageyama and Sato (2004) Kageyama, A. and Sato, T., 2004. âyin-yang gridâ: An overset grid in spherical geometry. Geochemistry, Geophysics, Geosystems, 5 (9).
- Kim et al. (2001) Kim, J., Kim, D., and Choi, H., 2001. An Immersed-Boundary finite-volume method for simulations of flow in complex geometries. Journal of Computational Physics, 171 (1), 132–150.
- Li et al. (2009) Li, Y., et al., 2009. Prediction of vortex shedding from a circular cylinder using a volumetric Lattice-Boltzmann boundary approach. European Physical Journal: Special Topics, 171 (1), 91–97.
- Liu and Hu (2014) Liu, C. and Hu, C., 2014. An efficient immersed boundary treatment for complex moving object. Journal of Computational Physics, 274, 654–680.
- Luo et al. (2017) Luo, K., et al., 2017. A ghost-cell immersed boundary method for the simulations of heat transfer in compressible flows under different boundary conditions part-ii: Complex geometries. International Journal of Heat and Mass Transfer, 104, 98–111.
- Luo et al. (2016) Luo, K., et al., 2016. A ghost-cell immersed boundary method for simulations of heat transfer in compressible flows under different boundary conditions. International Journal of Heat and Mass Transfer, 92, 708–717.
- Mattsson and Nordström (2004) Mattsson, K. and Nordström, J., 2004. Summation by parts operators for finite difference approximations of second derivatives. Journal of Computational Physics, 199 (2), 503–540.
- Mavriplis (1997) Mavriplis, D., 1997. Unstructured grid techniques. Annual Review of Fluid Mechanics, 29 (1), 473–514.
- Meakin (1995) Meakin, R., 1995. The chimera method of simulations for unsteady three-dimensional viscous flow. CFD Review, 70–86.
- Mittal et al. (2008) Mittal, R., et al., 2008. A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries. Journal of Computational Physics, 227 (10), 4825–4852.
- Mittal and Iaccarino (2005) Mittal, R. and Iaccarino, G., 2005. Immersed Boundary Methods. Annual Review of Fluid Mechanics, 37 (1), 239–261.
- Mittal (2005) Mittal, S., 2005. Excitation of shear layer instability in flow past a cylinder at low reynolds number. International Journal for Numerical Methods in Fluids, 49 (10), 1147–1167.
- Nagendra et al. (2014) Nagendra, K., Tafti, D.K., and Viswanath, K., 2014. A new approach for conjugate heat transfer problems using immersed boundary method for curvilinear grid based solvers. Journal of Computational Physics, 267, 225–246.
- Noack (2005) Noack, R., 2005. Suggar: a general capability for moving body overset grid assembly. In: 17th AIAA Computational Fluid Dynamics Conference. 5117.
- Owen (1998) Owen, S.J., 1998. A survey of unstructured mesh generation technology. In: Proceedings of 7th international meshing roundtable, Dearborn, MI, USA. Sandia National Laboratories, 239–267.
- Pan (2006) Pan, D., 2006. An immersed boundary method on unstructured cartesian meshes for incompressible flows with heat transfer. Numerical Heat Transfer, Part B: Fundamentals, 49 (3), 277–297.
- Park et al. (1998) Park, J., Kwon, K., and Choi, H., 1998. Numerical solutions of flow past a circular cylinder at Reynolds numbers up to 160. KSME International Journal, 12 (6), 1200–1205.
- Pärt-Enander and Sjögreen (1994) Pärt-Enander, E. and Sjögreen, B., 1994. Conservative and non-conservative interpolation between overlapping grids for finite volume solutions of hyperbolic problems. Computers & Fluids, 23 (3), 551–574.
- Peskin (1972) Peskin, C.S., 1972. Flow patterns around heart valves: A numerical method. Journal of Computational Physics, 10 (2), 252–271.
- Peskin (2002) Peskin, C.S., 2002. The immersed boundary method. Acta Numerica, 11, 479–517.
- Pletcher et al. (2012) Pletcher, R.H., Tannehill, J.C., and Anderson, D., 2012. Computational fluid mechanics and heat transfer. CRC Press.
- Poinsot and Lele (1992) Poinsot, T.J. and Lele, S., 1992. Boundary conditions for direct simulations of compressible viscous flows. Journal of Computational Physics, 101 (1), 104–129.
- Posdziech and Grundmann (2007) Posdziech, O. and Grundmann, R., 2007. A systematic approach to the numerical calculation of fundamental quantities of the two-dimensional flow over a circular cylinder. Journal of Fluids and Structures, 23 (3), 479–499.
- Qu et al. (2013) Qu, L., et al., 2013. Quantitative numerical analysis of flow past a circular cylinder at Reynolds number between 50 and 200. Journal of Fluids and Structures, 39, 347–370.
- Quirk (1994) Quirk, J.J., 1994. An alternative to unstructured grids for computing gas dynamic flows around arbitrarily complex two-dimensional bodies. Computers & Fluids, 23 (1), 125–142.
- Seo and Mittal (2011) Seo, J.H. and Mittal, R., 2011. A high-order immersed boundary method for acoustic wave scattering and low-mach number flow-induced sound in complex geometries. Journal of Computational Physics, 230 (4), 1000–1019.
- Sherer and Scott (2005) Sherer, S.E. and Scott, J.N., 2005. High-order compact finite-difference methods on general overset grids. Journal of Computational Physics, 210 (2), 459–496.
- Shi et al. (2004) Shi, J.M., et al., 2004. Heating effect on steady and unsteady horizontal laminar flow of air past a circular cylinder. Physics of Fluids, 16 (12), 4331–4345.
- Stålberg et al. (2006) Stålberg, E., et al., 2006. High order accurate solution of flow past a circular cylinder. Journal of Scientific Computing, 27 (1-3), 431–441.
- Steger et al. (1983) Steger, J.L., Dougherty, F.C., and Benek, J.A., 1983. A chimera grid scheme. In: K.N. Ghia and U. Ghia, eds., Advances in Grid Generation, ASME FED-5. 59–69.
- Steger and Benek (1987) Steger, J.L. and Benek, J.A., 1987. On the use of composite grid schemes in computational aerodynamics. Computer Methods in Applied Mechanics and Engineering, 64 (1-3), 301–320.
- Strand (1994) Strand, B., 1994. Summation by parts for finite difference approximations for d/dx. Journal of Computational Physics, 110 (1), 47–67.
- Suhs et al. (2002) Suhs, N., Rogers, S., and Dietz, W., 2002. Pegasus 5: An automated pre-processor for overset-grid cfd. In: 32nd AIAA Fluid Dynamics Conference and Exhibit. 3186.
- The Pencil Code (2018) The Pencil Code, 2018. (pencil-code.nordita.org [Internet]). Stockholm (SE): NORDITA; [updated February 2018]. Available from: https://github.com/pencil-code.
- Tseng and Ferziger (2003) Tseng, Y.H. and Ferziger, J.H., 2003. A ghost-cell immersed boundary method for flow in complex geometry. Journal of Computational Physics, 192 (2), 593–623.
- Uddin et al. (2014) Uddin, H., Kramer, R., and Pantano, C., 2014. A cartesian-based embedded geometry technique with adaptive high-order finite differences for compressible flow around complex geometries. Journal of Computational Physics, 262, 379–407.
- Versteeg and Malalasekera (2007) Versteeg, H.K. and Malalasekera, W., 2007. An Introduction to Computational Fluid Dynamics: the Finite Volume Method. Pearson Education.
- Visbal and Gaitonde (1999) Visbal, M.R. and Gaitonde, D.V., 1999. High-order-accurate methods for complex unsteady subsonic flows. AIAA Journal, 37 (10), 1231–1239.
- Visbal and Gaitonde (2002) Visbal, M.R. and Gaitonde, D.V., 2002. On the use of higher-order finite-difference schemes on curvilinear and deforming meshes. Journal of Computational Physics, 181 (1), 155–185.
- Völkner et al. (2017) Völkner, S., Brunswig, J., and Rung, T., 2017. Analysis of non-conservative interpolation techniques in overset grid finite-volume methods. Computers & Fluids, 148, 39–55.
- White (2006) White, F.M., 2006. Viscous Fluid Flow. 3rd ed. McGraw-Hill Higher Education Boston.
- Williamson (1980) Williamson, J., 1980. Low-storage runge-kutta schemes. Journal of Computational Physics, 35 (1), 48–56.
- Xia et al. (2014) Xia, J., Luo, K., and Fan, J., 2014. A ghost-cell based high-order immersed boundary method for inter-phase heat transfer simulation. International Journal of Heat and Mass Transfer, 75, 302–312.
- Zang and Street (1995) Zang, Y. and Street, R.L., 1995. A composite multigrid method for calculating unsteady incompressible flows in geometrically complex domains. International Journal for Numerical Methods in Fluids, 20 (5), 341–361.
- Zheng et al. (2009) Zheng, X., et al., 2009. A computational study of the effect of false vocal folds on glottal flow and vocal fold vibration during phonation. Annals of Biomedical Engineering, 37 (3), 625–642.