Treatment of solid objects in the Pencil Code using an immersed boundary method and overset grids
Abstract
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 bodyconformal 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
1 Introduction
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., bodyfitted structured meshes are commonly used to represent the object(s) in the flow. Bodyfitted 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 finitevolume or finiteelement 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 highorder 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 cutcell 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 cutcell 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 bodyforce, present due to nonconforming boundaries in the flow, is introduced in the NavierStokes 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 ghostpoints inside the solid and mirror/imagepoints 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 highorder accurate reconstructions of velocities near the surface can be overcome (Seo and Mittal, 2011; Xia et al., 2014). Further, finitedifference IBMs are, in general, not mass conserving. Finitevolume approaches with cutcell 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 bodyconformal 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 yinyang 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 finitevolume implementations of overset grids, see PärtEnander and Sjögreen (1994) and Zang and Street (1995)). The interpolation is done from donorpoints on one grid to fringepoints on another. Many different interpolation procedures have been explored for this purpose, and several studies have found that using highorder 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 highorder 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 finitedifference discretization, time integration, communication procedures, etc, are used. The purpose of this paper is to perform such a comparison for solids represented by a ghostpoints 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 semicircular and a semielliptical 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 Methodology
2.1 Governing equations
The governing equations of the flow are the continuity equation:
(1) 
and the momentum equation:
(2) 
where , , , and are the density, time, velocity vector, pressure, and dynamic viscosity (, with kinematic viscosity ), respectively , and
(3) 
is the material derivative operator. The compressible rate of strain tensor is given by:
(4) 
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 rewritten to:
(5) 
which is the form solved in the computations performed in this study.
2.2 Numerical methods
The governing equations (Eqs. (1) and (5)) are discretized with sixthorder finitedifferences in space and a thirdorder memory efficient RungeKutta scheme in time (Williamson, 1980).
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 NavierStokes 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 onedimensional 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 noslip 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 Ghostzone 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 fluidpoints, where the governing equations are solved, while the intersection points of dashed grid lines are grid points inside an immersed solid (solidpoints). At the solidpoints the governing equations are not solved. Rather, some points are used as ghostpoints 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 ghostpoint immersed boundary method. The IBM in the Pencil Code is such a ghostpoint method. An uncommon feature of the IBM implementation in the Pencil Code is that a several points deep ghostzone is used, rather than a single strip of ghostpoint inside the solid object. This ensures that the sixthorder finitedifference stencils can be used without any modifications in the vicinity of a solid object. The overhead related to computation of two additional layers of ghostpoints is negligible when compared to the computational cost of the fluid solver itself.
As sixthorder 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 fluidpoints 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 fluidpoints to the left). The points in the ghostzone, and are set using corresponding mirrorpoints and , respectively. With noslip and impermeability for velocity and zero gradient for density, the relationship between a ghost point and mirrorpoint is simply:
(6)  
(7) 
Note, that a secondorder 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 mirrorpoints do not coincide with grid points or grid lines, and need therefore to be interpolated from surrounding points. This is done by bilinear Lagrangian interpolation in 2D and trilinear Lagrangian interpolation in 3D. For mirrorpoints 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 fluidpoints, these mirrorpoints 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 fluidpoints ( 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 mirrorpoint near the surface the velocity is decomposed into cylindrical components, and the radial component is computed by:
(8) 
where is the radial velocity at the grid line interception point, and and are the distances from the mirrorpoint 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 fluidpoints very close to the surface of an object. This is done to avoid spurious errors due to delocalization dependencies in the finitedifference 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 fluidpoint close to the surface in the usual way, by using the finitedifference 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 mirrorpoints 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 mirrorpoint positions using surface normals is to use mirrorpoints along grid lines. This simplifies interpolation (making all interpolation onedimensional, along grid lines) and has proven promising in reducing the errors due to delocalization 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 finitedifference stencils (as would be the case for and in Fig. 2(a)). Flow variables with Neumann boundary conditions are set from mirrorpoints along surface normals as in the method described above. In this study, we will stick to the more mature method of using mirrorpoints from surface normals. Details on the alternative gridline ghostzone IBM can be found in Aarnes et al. (2018b).
2.2.2 Localtime restricted overset grids
Unlike the representation of solid bodies with most methods (IBM, bodyfitted 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). Yinyang grids are not considered (although a yinyang grid implementation exists in the Pencil Code).
For a flow with a solid object represented using a bodyconfined 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 bodyfitted grid, (3) solution of the governing equations on the bodyfitted grid, (4) communication of data from the bodyfitted grid to the background grid. The solution step on the bodyfitted grid requires implementation of a NavierStokes 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 NavierStokes 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 fringepoints. Fringepoints are fluidpoints set by interpolation from nearby donorpoints on the overlapping grid, rather than computed using discretization of the governing equations. Donorpoints are computed in the same way as an ordinary fluidpoint, unless the interpolation is implicit, meaning that a fringepoint may be used as a donorpoint (Chesshire and Henshaw, 1990), which is not the case here. A third class of points found on overset grids are holepoints. 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 holepoints on an overset grid generated to represent a circular object are seen in Fig. 3.
During the first stage of the intergrid communication, data is sent to the outermost grid points of the bodyfitted grid (Fig. 3(a)). At the second stage, data is sent back to Cartesian fringepoints (Fig. 3(b)). The fringepoints on the Cartesian grid are identified during preprocessing, where points within a distance from the solid, illustrated by red circular curves in Fig. 3(b), are selected. For moving objects, fringepoint locations on the background grid must be recalculated every time step. A several layers thick zone of fringepoints is used on both grids to enable the use of the sixthorder centered stencils at the outer edges of each grid, equivalent to the use of a ghostzone for the IBM.
In the illustration in Fig. 3 each fringepoint is surrounded by four donorpoints, as necessary in bilinear 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 highorder 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. Highorder Lagrangian interpolation and quadratic splines are implemented for overset grid interpolation in the Pencil Code. For simplicity, we will restrict ourselves to bilinear Lagrangian interpolation here.
At the solidfluid interface, we use summationbyparts (SBP) finitedifference operators to enhance stability of the solution. This means modifying the finitedifference stencils in the nine grid points closest to the surface (including the surface point) along each radial grid line, to asymmetric stencils (onesided at the surface). The order of accuracy for the SBPoperators are thirdorder for a sixthorder finitedifference 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 bodyfitted 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 bodyfitted 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 bodyfitted grid the reader is referred to Aarnes et al. (2018a).
2.2.3 A note on dissipation
The centered finitedifference schemes used for discretization of the governing equations are nondissipative. This can cause problems due to the potential growth of highfrequency modes, leading to numerical instability.
To some extent, the summationbyparts 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 highfrequency modes increases as the grid spacing decreases, and may lead to diverging solutions as the grid is refined. To suppress the highfrequency modes, a highorder lowpass 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:
(9) 
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 tridiagonal 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 sixthorder hyperdiffusion operator, which is already implemented in the PencilCode (see e.g. Haugen and Brandenburg, 2006), could have been used to filter the solution. The benefits of this approach is that the hyperdiffusion 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 fifthorder upwinding for the advection operators of the density rather than centraldifference 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 secondorder 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 fifthorder, while the accuracy of the tangential velocity and density was between second and thirdorder, when secondorder 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. Twodimensional 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 threedimensional effects in the latter flow are suppressed as we restrict ourselves to a twodimensional domain. A rectangular domain with domain size is used, with a bodyfitted 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 rootmeansquare 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 twodimensional 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 tradeoff 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, socalled forcepoints are used. The forcepoints 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 rootmeansquare lift coefficients seen in Fig. 5 may be due to a change in the position of forcepoints 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 fringepoints are interpolated with overset grids than mirrorpoints interpolated with the ghostpoint 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 ghostzone on the Cartesian grid in overset grids since no special handling for fluidpoints close to the fringepoint region is used. For the grid, the total number of fringepoints for the overset grids is approximately 1600 (% of total number of grid points). Note that the fraction of fringepoints to grid points decreases as the grid spacing decreases ( grid points are fringepoints when ). The number of mirrorpoints in an IBM simulation with the same grid spacing is approximately 200. A fringepoint is updated only once every RungeKutta time step, while mirrorpoints are updated every subtime 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  

Pan (2006)  1.32  0.16  
Haugen and Kragset (2010)  1.328  0.166  
Park et al. (1998)  ,  1.33  0.165  
Shi et al. (2004)  1.318  0.164  
Mittal (2005)  1.322  0.226  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, IBM  1.351  0.232  0.166  
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, rootmeansquare 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 rootmeansquare 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 sinusoidallike 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 finitevolume, finitedifference, finiteelement, spectral element and latticeBoltzmann 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 bodyfitted 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 rootmeansquare 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 forcepoints 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 finitedifferences 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 ghostzone 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 ghostcell 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 levelset 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 twoway 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 twodimensional immersed boundaries are represented by many small linesegments. Each linesegment 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 fluidpoints and mirrorpoints around the solid object.
The first difference lies in the identification of a given grid point to be a fluidpoint or a solidpoint. 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 linesegment’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 solidpoint 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/fluidpoints, three layers of ghost points are assigned to construct a sixorder central finitedifference 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 mirrorpoint can be obtained. The mirrorpoints are set either by symmetry over the solid’s line element (corresponding to the way mirrorpoints 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 3Dgeometries, to ensure that every mirrorpoint is surrounded by fluidpoints only. The same interpolation procedure as for the simple circular geometry (i.e., the bilinear interpolation method) can be adopted for the calculations of the parameters of the mirrorpoints. An optional way is the inverse distance weight interpolation method (Chaudhuri et al., 2011). Finally, the flow variables at the ghostpoints can be calculated with the aid of the mirrorpoints 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 threedimensional 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 bilinear interpolation procedure can obtain a local secondorder accuracy. Systematic validations have also been conducted through calculations of flow past an elliptical cylinder, square cylinder, semicylinder, 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 twodimensional flow, we have simulated flow past geometries constructed by combining a semicircle and semielliptical cylinder. The geometry is seen in Fig. 8. The radius of the semicircle () and the major axis of the semiellipse () 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 semiellipse 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, rootmeansquare 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.
2.0  1.55  0.34  0.190 
1.0  1.35  0.26  0.175 
0.5  1.17  0.11  0.168 
5 Concluding remarks
In this study we have described and compared the two solid body representations available in the highorder finitedifference 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 ghostpoint 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 mirrorpoints or nonconforming 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 twodimensional 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 45 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 opensource software. Perhaps, in time, overset grids can become more flexible, and the highorder 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, highorder interpolation of mirrorpoints 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 opensource software development.
Acknowledgements
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 pencilcode directory, the sample cases can be found in:
./samples/2dtests/cylinder_deposition
./samples/2dtests/cylinder_deposition_ogrid
The postfix ogrid denotes the overset grid sample case. To compile and run a sample case, use the commands pc_build
, pc_start
and pc_run
.
Both sample cases are simulations of a particleladen 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.
References
 Aarnes et al. (2018a) Aarnes, J.R., Andersson, H.I., and Haugen, N.E.L., 2018a. Highorder 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 freestream turbulence effects on the transitioninwake 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 3d 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 highaccuracy 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 higherorder boundary filters for the NavierStokes 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 cutcell 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 (36), 561–572.
 Kageyama and Sato (2004) Kageyama, A. and Sato, T., 2004. âyinyang 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 ImmersedBoundary finitevolume 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 LatticeBoltzmann 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 ghostcell immersed boundary method for the simulations of heat transfer in compressible flows under different boundary conditions partii: Complex geometries. International Journal of Heat and Mass Transfer, 104, 98–111.
 Luo et al. (2016) Luo, K., et al., 2016. A ghostcell 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 threedimensional 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ärtEnander and Sjögreen (1994) PärtEnander, E. and Sjögreen, B., 1994. Conservative and nonconservative 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 twodimensional 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 twodimensional bodies. Computers & Fluids, 23 (1), 125–142.
 Seo and Mittal (2011) Seo, J.H. and Mittal, R., 2011. A highorder immersed boundary method for acoustic wave scattering and lowmach number flowinduced sound in complex geometries. Journal of Computational Physics, 230 (4), 1000–1019.
 Sherer and Scott (2005) Sherer, S.E. and Scott, J.N., 2005. Highorder compact finitedifference 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 (13), 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 FED5. 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 (13), 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 preprocessor for oversetgrid cfd. In: 32nd AIAA Fluid Dynamics Conference and Exhibit. 3186.
 The Pencil Code (2018) The Pencil Code, 2018. (pencilcode.nordita.org [Internet]). Stockholm (SE): NORDITA; [updated February 2018]. Available from: https://github.com/pencilcode.
 Tseng and Ferziger (2003) Tseng, Y.H. and Ferziger, J.H., 2003. A ghostcell 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 cartesianbased embedded geometry technique with adaptive highorder 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. Highorderaccurate 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 higherorder finitedifference 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 nonconservative interpolation techniques in overset grid finitevolume methods. Computers & Fluids, 148, 39–55.
 White (2006) White, F.M., 2006. Viscous Fluid Flow. 3rd ed. McGrawHill Higher Education Boston.
 Williamson (1980) Williamson, J., 1980. Lowstorage rungekutta schemes. Journal of Computational Physics, 35 (1), 48–56.
 Xia et al. (2014) Xia, J., Luo, K., and Fan, J., 2014. A ghostcell based highorder immersed boundary method for interphase 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.