A Moving Boundary Flux Stabilization Method for Cartesian Cut-Cell Grids using Directional Operator Splitting

A Moving Boundary Flux Stabilization Method for Cartesian Cut-Cell Grids using Directional Operator Splitting

W.P. Bennett wpb22@cam.ac.uk    N. Nikiforakis nn10005@cam.ac.uk    R. Klein rupert.klein@math.fu-berlin.de University of Cambridge,
Laboratory for Scientific Computing, Cavendish Laboratory,
Cambridge CB3 0HE, UK
Freie Universität Berlin, Institut für Mathematik, Arnimallee 6, 14195 Berlin

An explicit moving boundary method for the numerical solution of time-dependent hyperbolic conservation laws on grids produced by the intersection of complex geometries with a regular Cartesian grid is presented. As it employs directional operator splitting, implementation of the scheme is rather straightforward. Extending the method for static walls from Klein et al., Phil. Trans. Roy. Soc., A 367, no. 1907, 4559–4575 (2009), the scheme calculates fluxes needed for a conservative update of the near-wall cut-cells as linear combinations of “standard fluxes” from a one-dimensional extended stencil. Here the standard fluxes are those obtained without regard to the small sub-cell problem, and the linear combination weights involve detailed information regarding the cut-cell geometry. This linear combination of standard fluxes stabilizes the updates such that the time-step yielding marginal stability for arbitrarily small cut-cells is of the same order as that for regular cells. Moreover, it renders the approach compatible with a wide range of existing numerical flux-approximation methods. The scheme is extended here to time dependent rigid boundaries by reformulating the linear combination weights of the stabilizing flux stencil to account for the time dependence of cut-cell volume and interface area fractions. The two-dimensional tests discussed include advection in a channel oriented at an oblique angle to the Cartesian computational mesh, cylinders with circular and triangular cross-section passing through a stationary shock wave, a piston moving through an open-ended shock tube, and the flow around an oscillating NACA 0012 aerofoil profile.

Keywords: cut-cell; embedded boundary; dimension splitting; moving boundary; compressible flow; flux stabilization

preprint: thanks: Corresponding Author

1 Introduction

Embedded boundary (or immersed boundary) methods comprise a wide range of techniques for discretising complex geometries ‘embedded’ into a simple Cartesian mesh. Cartesian mesh-based solvers can yield significant advantages over body-fitted and unstructured mesh methods in terms of mesh generation simplicity, reduction in cell skewness around complex geometries, automated mesh generation and improved parallelisation. This is particularly true when considering the movement of solid boundaries over time, where non-Cartesian mesh-based alternatives can require expensive regeneration of the mesh at each time step.

One class of embedded boundary methods in common use employs artificial volume forces distributed over a near-boundary zone of a few grid cells to effectively represent the pressure force exerted by a rigid wall, peskin1972flow (); goldstein1993modeling (); fadlun2000combined (); peskin2002immersed (); tyagi2005large (). The advantage of this approach is that the wall geometry appears only in the distributed momentum forces, while the surface geometry does not otherwise affect the scheme. In particular, intersections of the wall geometry with computational cells do not have to be accounted for explicitly.

“Ghost-fluid methods” fedkiw1999non (); fedkiw2002coupling () belong, in contrast, to the class of sharp interface schemes. That is, the sub-cell geometries of Cartesian grid cells cut by the boundary are explicitly accounted for, and no artificial smearing of the solution occurs near the boundaries. These schemes work by first extrapolating the solution at a given time step beyond the actual computational domain into a halo of “ghost cells” in the immediate vicinity of the boundary, and then employing standard flux balances for complete cells to advance the solution. The wall geometry is re-inserted into the updated full grid cells after a completed time step. A disadvantage of ghost fluid methods is that they are generally not conservative with respect to individual cut-cells, such that mass, momentum, and energy losses across a rigid wall may occur.

We limit our literature review to cut-cell methods. The field of Immersed Boundary Methods in general is large. An exhaustive review of this field is beyond the scope of the current work. For more information, see mittal2005immersed ().

Cut-cell methods provide a conservative alternative, where Cartesian cells underlying the solid boundary are each divided into constituent components. For flows involving a solid geometry in a surrounding fluid, the Cartesian cell underlying the boundary is divided into fluid and solid components. The boundary is therefore simply ‘cut out’ of the surrounding Cartesian mesh. Boundary conditions are subsequently imposed along the intersecting surface. Dividing the Cartesian cell in this manner does, however, generate small sub-cells at the boundary which may feature arbitrarily small volume and highly distorted geometry. These ‘small cells’, which may be orders of magnitude smaller than the surrounding Cartesian cells, create issues of stability and efficiency for explicit time integration methods. This is termed the “small cell problem”. A brief outline of existing conservative cut-cell approaches can now be described.

Flux redistribution schemes Pember95 (); Colella06 (); HuEtAl2006 (); schneiders2013accurate () are conservative and obtain stable updates for small cut-cells at time steps comparable to those achievable on a regular grid in two steps. First, cut-cell updates are determined from standard flux balances that are locally computed and disregard the small sub-cell problem. In the second step, a suitable fraction of these updates is redistributed to a stencil of neighbouring cells. Colella Colella1990 () argues that handling the time updates for cut-cells emerging from the intersection of a regular Cartesian grid with a complex geometry would necessitate “unsplit” numerical schemes, i.e., schemes that avoid directional operator splitting, and accordingly introduced an unsplit Godunov-type higher-order scheme in this paper. Yet, in Klein09 () we have demonstrated that stable individual updates for small cut-cells can be formulated in an operator splitting scheme, and we show here that the same methodology can be extended to moving walls.

Recent developments to the method of Schneiders et al. schneiders2013accurate () has extended the method to describe freely moving boundaries and multi cut-cells, alongside improvements in accuracy and efficiency through an improved interpolation stage. The improved method is described in Schneiders et al. schneiders2016efficient (), and is used to model particles in isotropic turbulence by Schneiders et al. schneiders2017accuracy (); schneiders2017direct ().

“Cell merging” methods Clarke85 (); Quirk94 (); Berger03 (); Ingram03 (); Xu97 (); yang1997cartesian (); barton2011conservative (), avoid the small sub-cell problem by creating larger effective control volumes of size comparable to that of a regular cell through merging cut-cells with neighbouring larger cells and then computing common updates for the merged cells. There is some arbitrariness in these schemes as to which cells and sub-cells are to be merged to create the enlarged effective control volumes, and a potential loss of resolution of the near-wall flow states, both of which we intend to avoid. Notably, Falcovitz et al. falcovitz1997two () designed a technique based on directional operator splitting, and this is conceptually closest to the present approach. It uses a cell merging for cut-cells with volume fractions less than a user defined cell ‘partial length’ to maintain stability in small cells. The cell merging approach has been shown to evoke spurious oscillations when applied to moving boundaries, as shown in for example, Schneiders et al. schneiders2013accurate (). This motivates the use of alternative approaches, such as the one discussed in this article.

An approach that extends the merged-cell technique to obtain individual updates for small cut-cells has been proposed by Hartmann et al. HartmannEtAl2011 (). In their scheme, a stable update for a merged cell is computed first, and then a multi-dimensional reconstruction algorithm is invoked to calculate a conservative improved local update for the cut-cell portion of the merged cell. It does not appear obvious how their scheme could be formulated in the context of directional operator splitting as is the goal of the present work. The method of Hartmann et al. HartmannEtAl2011 () has been extended recently by Pogorelov et al. pogorelov2015cut (); pogorelov2016effects (); pogorelov2017adaptive () to study axial turbomachinery flows using Large Eddy Simulation (LES) based on the Monotone Implicit Large Eddy Simulation (MILES) approach.

The “H-Box method” by Berger, Helzel et al. Berger03 (); Helzel05 () is a locally second order accurate method that achieves stability for small cut-cells by resorting to LeVeque’s wave propagation scheme, LeVequeBook2002 (), which avoids flux balancing for control volumes but rather obtains updates by superimposing wave-type increments that emerge from (linearized) Riemann problems solved at grid cell interfaces. Using Roe’s approach for the formulation of an approximate Riemann problem, Roe1981 (); LeVequeBook2002 (), such wave superposition schemes can be formulated such that the resulting scheme conserves the primary conserved quantities, i.e., mass, momentum, and energy for gas dynamics. Through this construction, the proposed embedded boundary scheme is inherently tied to the wave propagation approach, however, and is not flexible with respect to the use of a user’s flux approximation of choice. Also, the distribution of wave portions in the vicinity of cut cells involves complex geometrical computations. Thus, even though a substantially simpler variant of the original scheme was proposed in BergerHelzel2012 (), implementation of this approach in three space dimensions remains rather tedious. Recently, Mudalidharan & Menon muralidharan2016high () describe a high order (third order inviscid fluxes and fourth order viscous fluxes) conservative finite volume cut-cell method that uses a high order reconstruction using a cell-centre piece-wise polynomial approximation. A cell merging technique is used for small cut-cells. Mudalidharan & Menon muralidharan2016high () attempt to overcome numerical oscillations at the cut-cell boundary caused by the irregular stencil used in the k-exact CENO reconstruction by applying the piece-wise polynomial reconstruction to clusters of cells. The method is not, however, demonstrated for high Mach number compressible flows such as those of interest in the current work and is cited as an avenue for further work. Other recent progress in the development of incompressible higher order cut-cell methods is Krause & Kummer krause2017incompressible () in which a Discontinuous Galerkin (DG) discretization approach is developed based on the extended Discontinuous Galerkin (XDG) approach of Kummer kummer2017extended (). This method is simply mentioned here for completeness, as we limit the scope of this literature review to compressible finite volume methods.

The challenges involved in introducing a cut-cell technique are exacerbated when a moving boundary is required. In addition to problems caused by continuous reduction in the volume fraction of a single cut-cell over time, the characteristics of a cut-cell (and the neighbouring regular cells) generally alter over time. This can result in an individual cell changing from a cut-cell to a solid cell over a single time step, a regular cell changing to a cut-cell, or a large cut-cell (volume fraction ) changing to a small cut-cell (or the opposite of all three scenarios). This has significant implications for any of the embedded boundary techniques summarized above. In addition, any cut-cell method development must satisfy the consistency condition, whereby a solid body moving at the same relative velocity as the surrounding fluid should produce an exact solution. Published literature on cut-cell methods for moving boundary problems include cell merging yang1997cartesian (); barton2011conservative () methods, implicit time-stepping aftosmis2000parallel (); murman2003implicit () and flux redistribution methods MeinkeEtAl2013 (); schneiders2013accurate (). A review of embedded boundary methods is given by Mittal et al. mittal2005immersed ().

We describe an alternative technique for the solution of moving boundary flows using a dimension-splitting cut-cell approach, derived from the static wall method introduced by Klein et al. Klein09 (). Our scheme calculates fluxes needed for a conservative update of the near-wall cut-cells as linear combinations of “standard fluxes” from a one-dimensional extended stencil. Here the standard fluxes are those obtained without regard to the small sub-cell problem, and the linear combination weights involve detailed information regarding the cut-cell geometry. This linear combination of standard fluxes stabilizes the updates such that the time-step yielding marginal stability for arbitrarily small cut-cells is of the same order as that for regular cells. Moreover, it renders the approach compatible with a wide range of existing numerical flux-approximation methods. The scheme is extended here to time dependent rigid boundaries by reformulating the linear combination weights of the stabilizing flux stencil to account for the time dependence of cut-cell volume and interface area fractions.

The basic scheme from Klein09 () and the extension to moving walls are described in section 2 below, while results for a number of smooth and discontinuous flow test cases are presented in section 3. Section 4 provides conclusions and an outlook to future work.

2 Numerical Method

2.1 Finite Volume Discretisation

We consider a computational domain containing a fluid region, with one or more embedded solid regions of essentially arbitrary shape. Each solid region can move freely through the fluid in any direction with velocity, , and rotation rate, , as shown in Figure 1 (a). Following standard finite volume methodology, the computational domain is divided into discrete cells in the th spatial direction. The solid region is then cut out of this domain, leaving just the fluid domain discretised by the Cartesian mesh, as shown in Figure 1 (b). The current cut-cell implementation represents a curved solid interface using a planar representation in each computational cell. The resultant curvature is therefore mesh dependent. The solid-fluid interface is imposed within each cut-cell using an appropriate boundary condition.

(a) (b)

Figure 1: Demonstrative low resolution Cartesian mesh containing embedded solid regions moving with velocity and rotation rate . (a) Physical geometry with indicative signed distance function, . (b) Computational domain discretization using a Cartesian mesh.

In this paper, we consider an inviscid, single-phase gas for the fluid region. The compressible Euler equations comprise of the mass, momentum and energy conservation equations. These are defined in compact form as:


where the vector of conserved variables is given by and the flux vector is given by . These equations are closed using the perfect gas equation of state, , with a constant isentropic exponent, .

Following integration of Equation 1 over an arbitrary control volume, , and after applying the Gauss divergence theorem, the discrete form of the Euler equations can be obtained from the resultant integral form by replacing the flux integral around the control volume by a summation of the surface normal fluxes at each computational cell interface. For regular Cartesian cells, away from the cut-cell surface, this is given by:


where the shorthand cell index, and is the regular cell volume. Cut-cells are defined as having volume, , where the cut-cell volume fraction, , is defined as the ratio of the fluid volume in the cut-cell to the fluid volume of a regular Cartesian cell, at time step . In the current moving boundary formulation, the solid region is permitted to move over the course of each time step. This movement results in a change to the fraction of the Cartesian cell covered by the solid, altering the fluid volume fraction, , between time steps and . The discrete form of the Euler equations applied to a moving boundary cut-cell is then given by:


In this formulation, is an additional contribution due to the solid wall movement, denoted here as a wall ‘flux’ for convenience. The kronecker delta, , is therefore unity at solid boundaries and zero otherwise. The conservative flux is required normal to both inter-cell boundaries and solid-fluid boundaries, . The total number of cell faces in a Cartesian cut-cell may be more than, less than, or equal to the number of faces in a regular Cartesian cell. The total number of cut-cell faces is dependent on the local discretisation of the solid geometry and the geometry curvature. The moving boundary cut-cell approach described in this paper is a dimension-splitting approach, which isolates the flux component in each spatial direction, considering each as a separate Strang-splitting stage strang1968construction (); LeVequeBook2002 (), in which the operator sequence is chosen to maintain second-order accuracy.

2.2 Static Boundary Method Overview

The conservative, dimension-splitting cut-cell method of Klein et al. Klein09 () is extended in this paper to flows involving one or more moving solid boundaries. This section starts with a brief overview of the static boundary method of Klein et al. Klein09 () to provide context before describing the current extension to moving boundary flows. Figure 2 shows a two-dimensional solid-fluid interface extending across three cut-cells, with indices , and . Considering the -sweep i.e., the horizontal sweep, in isolation, the cut-cell face area fluid fractions at and are given by and respectively. These are labelled in Figure 2 (a). The difference between these cut-cell face area fractions, for cell , is given by . This represents the projection of the solid-fluid interface onto the cut-cell face, .


Figure 2: Demonstrative solid-fluid interface in three cut-cells. (a) Interface and volume fractions for the two-dimensional cut-cell method of Klein et al. Klein09 (). (b) Flux components for the two-dimensional cut-cell method of Klein et al. Klein09 ().

The flux stabilization method of Klein et al. Klein09 () permits a global time step to be imposed in both regular Cartesian cells and cut-cells. This is achieved by modifying the regular cell flux at inter-cell boundaries bordering a cut-cell. To describe this flux ‘stabilization’, we consider the cut-cell face in Figure 2 (a). The projection of the solid-fluid interface onto the cell face , given by , is termed in Klein et al. Klein09 () the ‘shielded’ region. Specifically, this region is directly influenced in the x-direction by the presence of the solid-fluid interface. The remaining fluid portion of the cut-cell face area, i.e., , is the interface portion not directly influenced by the solid-fluid interface. This is termed the ‘unshielded’ region.

Figure 2 (b) shows the shielded and unshielded portions of the inter-cell flux for the x-direction sweep along the interface . These fluxes are defined as and respectively. The total flux along this cell interface, , is defined by a weighted combination of the shielded and unshielded fluxes, as:


For the convex geometry shown in Figure 2, the unshielded flux component, , is simply the regular Cartesian cell flux:


The shielded flux creates a stable flux by combining the boundary flux, , with the regular cell flux, , as:


where is the average distance between the interface at and the solid-fluid interface, as shown in Figure 2 (b). The average distance, , is computed from the shielded volume fraction, , shown in Figure 2 (a), as:


For clarity, the notation of the average distance from the cut-cell interface to the stabilized cell interface, , is distinct from the cut-cell shielded volume fraction, , and the total cut-cell volume fraction, . The total cut-cell volume fraction is defined as the combination of the shielded and unshielded cut-cell volume fractions. The shielded flux in Equation 6 has the property that, as the cut-cell volume fraction is reduced, and therefore , then . As the cut-cell tends towards a regular cell, , the shielded flux tends towards the regular cell flux, .

The boundary flux, , is shown in Figure 2 (b) along the solid-fluid interface. The flux at this interface is computed using a reflective wall boundary condition. Caution is required here, as in order that the dimension splitting remains conservative, the advective part of the boundary flux in all directions must be computed for a common solid-fluid interface state. That is, we can separate the full flux at the boundary into an advective flux and a pressure flux:


where, the x-direction (horizontal) advective and pressure fluxes for the Euler equations in two-dimensions are respectively and based on the boundary state, . Similarly, the advective flux and pressure flux in the y-direction are given by and respectively. The advective flux is computed before the onset of the first dimensional sweep, at the start of each time-step. The pressure flux in the momentum equation is, however, computed at each dimensional splitting stage using the updated boundary state.

The neighbouring fluxes, at and , are also shown in Figure 2 (b). These fluxes are, in turn, a combination of the local shielded and unshielded flux contributions from the solid-fluid interface in cells and respectively for this example.

The complete two-dimensional update for cell is therefore given by:


Klein et al. Klein09 () show that, using the flux stabilization from Equation 6, this update can be rewritten for cut-cells in such a way that the division by is eliminated and singular behaviour for vanishing cut-cell volume fractions is avoided.

2.3 Moving Boundary Extension

2.3.1 Moving Boundary Issues

A number of additional challenges arise in the extension of the static cut-cell method of Klein et al. Klein09 () to accommodate moving solids. Whereas the status of the cut-cells and regular Cartesian fluid cells are unchanging over the course of a static boundary computation, a transient cut-cell status now exists due to the motion of the solid-fluid interface across each Cartesian cell. As the solid-fluid interface first crosses one of the regular cell boundaries, the regular Cartesian cell changes to a cut-cell which requires stabilization of the boundary fluxes, Figure 3 (a). As the solid-fluid interface moves further into the cut-cell, the cut-cell volume fraction becomes a smaller fraction of the regular cell fluid volume, Figure 3 (b), tending towards a zero volume fraction when the solid-fluid interface passes through the final cut-cell fluid boundary. At this point the cut-cell changes to a solid cell, which is no longer used in the fluid computation, Figure 3 (c).

(a) (b) (c)

Figure 3: Interaction of a moving solid-fluid interface with a cut-cell. (a) The solid first appears in the subject Cartesian cell, turning the regular cell into a cut-cell. (b) The cut-cell volume fraction gradually decreases as the solid covers more of the cut-cell. (c) The solid covers the Cartesian cell completely, turning the cut-cell into a solid cell. Solid line: Solid-fluid interface at the current time step. Dashed line: Solid-fluid interface at the previous time step. The arrow indicates the direction of the solid-fluid interface movement.

Conversely, a cut-cell can be created from a solid cell when the solid-fluid interface passes through one of the inter-cell boundaries of a solid cell. As the solid-fluid interface moves across the new cut-cell, the fluid volume fraction becomes an increasing fraction of the regular cell volume, until the solid-fluid interface passes out of the cut-cell and the cut-cell becomes a regular Cartesian fluid cell.

The problem caused by the decreasing cut-cell volume fraction, and the related changes of the cut-cell face area fractions, can be identified by rearranging the discrete Euler equation, Equation 3, to make the state variable the subject, i.e:


where the shorthand index, . As the cut-cell volume fraction tends towards zero over time, the reciprocal of the cut-cell volume fraction at the updated time, , causes the state variables at the new time, , to increase to large values. Merging cut-cell methods, as used for moving boundaries by Xu et al. Xu97 (), Yang et al. yang1997cartesian () and Barton et al. barton2011conservative (), for example, deal with this issue by grouping cut-cells of small volume fraction with an adjacent regular cell (or large cut-cell) to form one cell with a combined volume. In contrast, Schneiders et al. schneiders2013accurate (); schneiders2016efficient () deal with this issue through a redistribution of the numerical error from the small volume cut-cell to the surrounding cells, following a time integration of all cells using regular cell fluxes. The vanishing volume fraction issue is dealt with in the current method through an augmentation of the cut-cell shielded flux.

For a solid-fluid interface moving at a velocity that matches the local velocity, the work done on the fluid by the solid-fluid interface reduces to zero and a constant wall pressure is recovered, . This is described in, for example, Toro Toro (). In terms of the present cut-cell method, this constraint dictates that the mass, momentum or energy displaced by the solid-fluid boundary movement within a cut-cell should equal the total advective mass, momentum or energy flux through the computational cell boundaries, for each spatial direction. This provides a sound test of any new moving boundary method, as any resultant pressure fluctuations at the solid-fluid interface are therefore due to errors in matching the boundary movement to the advective flux, or due to fluctuations in the solid geometry through inadequate geometry discretisation and movement across the Cartesian mesh over time.

2.3.2 Moving Boundary Geometry Approximation

The moving boundary cut-cell method described in this paper is a straightforward extension to the static boundary method of Klein et al. Klein09 () described in Section 2.2. The moving boundary extension maintains the principle of dividing each cut-cell interface flux into a ‘shielded’ and ‘unshielded’ contribution, based on whether the inter-cell boundary directly faces the solid-fluid interface during each dimensional sweep.

For the purpose of illustration, we can consider the same solid-fluid interface as in the static boundary section, intersecting three cut-cells, , and , this time with movement of the solid-fluid interface over a typical time step, , as shown in Figure 4. The direction of the solid-fluid interface movement is indicated by arrows. Over the time interval, , the cut-cell interface in cell moves such that the intersection points with the Cartesian cell boundaries, and , at time become and at time . The cut-cell volume correspondingly changes from to .

Figure 4: Moving solid interface across three cut-cells. The solid-fluid interface in cell moves from A-B to A’-B’ between time level (dashed line) and (solid line).

Murman et al. murman2003implicit () classify the numerical representation of the space-time geometry development in increasing order of fidelity as either ‘sequential-static’, ‘staircase-in-time’, ‘linear-in-time’ or ‘full space-time approximation’. In terms of this classification, the current method would be considered a ‘linear-in-time’ approach, in which the geometry movement, and hence the cut-cell area and volume fractions, are assumed to change linearly between time steps and . The current moving boundary method is introduced by initially considering a simpler ‘moving wall’ boundary condition approach applied to the method of Klein et al. Klein09 () in Section 2.3.3.

2.3.3 The Moving Wall Boundary Condition Approach

The simplest extension to the method of Klein et al. Klein09 () involves fixing the geometry at each time level and changing the boundary flux to a ‘moving wall’ boundary flux, as described for a compressible Riemann-based approach in Toro Toro (). At the moving solid-fluid interface, the boundary flux incorporates the local boundary velocity by imposing a ghost-cell velocity, , according to:


where, is the local solid-fluid interface velocity vector and is the fluid velocity vector in the cut-cell , as shown in Figure 4. The ghost-cell density and pressure can be imposed for an inviscid reflective boundary condition according to and . Finding the resultant boundary state in terms of pressure, density and velocity results in a boundary velocity equal to the wall velocity, . For an advancing wall, the wall pressure, , and density, . Conversely, for a retreating wall, and . For a simple one-dimensional problem, approximating the non-penetration reflective boundary condition using this ‘moving wall’ boundary state results in a mass flux ‘through’ the solid-fluid interface which matches the fluid mass displaced by the movement of the solid-fluid interface over this time step. Using an example of a one-dimensional cut-cell reducing in length from to over the time interval , as shown in Figure 5, the cut-cell update is given by:

Figure 5: One-dimensional moving boundary example. The cut-cell reduces in size from to between times and , where . As the total cut-cell volume fraction, , the shielded cut-cell volume fraction, , and the average distance from the stabilized cell interface to the cut-cell interface, , are equivalent in a one-dimensional problem, we use the notation here.


Figure 6: Moving solid interface in three cut-cells. (a) Interface and volume fractions for the two-dimensional cut-cell method of Klein et al. Klein09 () based on the solid-fluid interface at time, . (b) Interface fluxes based on the solid fluid interface at time .

where the boundary flux, , and boundary state, , are non zero in mass, momentum and energy. We explicitly note that, in the one-dimensional case, the total cut-cell volume fraction, , the shielded cut-cell volume fraction, , and the average distance from the stabilized cell interface to the cut-cell interface, , are equivalent. For clarity, we therefore use the total cut-cell volume fraction notation, , in Equation 13 and Figure 5. It is noted that the moving boundary time integration/accuracy is explained further in Section 2.3.4.

The corresponding two-dimensional cut-cell update, for the cell shown in Figure 6, would then be given by:


where and are the wall velocity and direction components. is the boundary state. As shown in Ben-Artzi & Falcovitz ben2003generalized () and Falcovitz et al. falcovitz1997two (), the resultant contribution of the moving boundary in the and directions is then:


Thus, the contribution from the movement of the non-penetrative solid wall is through the boundary pressure in the momentum and energy equations. The boundary interface area that this pressure acts upon is the time averaged boundary area, as outlined in Section 2.3.4.

The necessity of precisely accounting for the movement of the solid boundary, in this case through a time averaged boundary in the conservation equation, is highlighted by Gokhale Gokhale14 () through consideration of the loss in conservation for a dimension-splitting method applied to a diagonal piston problem. This is presented here by considering the movement of the solid-fluid interface in the and directions independently, as shown in Figure 7. The resultant fluid volume displacement arising from the independent x-direction and y-direction movement does not necessarily equal the two-dimensional volume displacement.

(a) (b)

Figure 7: Comparison of the volume fraction change as a diagonal solid-fluid interface diagonally traverses a Cartesian cell. (a) Volume change after the x-direction sweep. (b) Volume change after the y-direction sweep. Dashed line: Solid-fluid interface at time level . Solid line: Solid-fluid interface at time level . Striped region: Volume change in each direction between time level and . Dash-dot line: Solid-fluid interface projection in the x-direction or y-direction.

As the rigid body interface moves across the Cartesian cell in Figure 7, the overall loss of conservation is accompanied by an under-prediction of the boundary work done by the solid-fluid interface. We attempt to preserve conservation for a moving rigid body by defining time averaged boundary properties in the conservation equations, together with a conservation check for newly covered cells. These are outlined in Section 2.3.4 and Section 2.3.5 respectively.

2.3.4 Time Averaged Boundary Modification

The moving boundary method developed here uses the time-averaged solid-fluid interface for defining the shielded and unshielded cell area fractions, , and the average shielded distance from the solid, . These can be formally defined as:


These are used for calculating , and . By using the time averaged cell areas, together with the conservative state after each dimensional sweep, the resultant flux through each cut-cell interface provides a better approximation than shown in Figure 7. Similarly, the solid-fluid interface unit normal is also computed as the time averaged surface normal.

(a) (b)

Figure 8: x-sweep and y-sweep for the diagonal moving piston. (a) x-sweep, starting from the dashed line, , and ending at the solid line, . Bold solid line: Time averaged boundary for the x-sweep. (b) y-sweep, starting from the dashed line, , and ending at the solid line, . Bold solid line: Time averaged boundary for the y-sweep.

The x-sweep and y-sweep for a single time step, between times and , is shown in Figure 8 for the same diagonal piston problem of a solid moving diagonally across a static fluid. Figure 8 (a) shows the motion of the solid during the x-sweep. The initial position of the solid-fluid interface is given by a dashed line. The final position of the solid-fluid interface, at the end of the x-sweep, is given by the solid line, with the solid shaded in grey. The final position of the solid-fluid interface, at the end of the time step , is given as a dash-dot line for reference. The time-averaged location of the solid-fluid interface during the x-sweep is given as a bold solid line, located between the dashed line and the thin solid line. Similarly, the y-sweep is shown in Figure 8 (b). The spatial location of the solid-fluid interface at the end of the x-sweep, , is given as a dashed line in this figure. This is the start position for the y-sweep. The location of the solid-fluid interface at the end of the y-sweep, at time , is given as a solid line with the solid at this time shaded in grey. The time-averaged location of the solid-fluid interface during the y-sweep is given by a bold solid line, located between the dashed line and the thin solid line. Following Ben-Artzi & Falcovitz ben2003generalized (), splitting the motion into individual x and y-sweeps (i.e. horizontal and vertical motion) simplifies the determination of the time averaged cut-cell face fraction, , the cut-cell volume fraction and hence the average cut-cell distance, compared to non-dimensional-splitting alternatives. For example, in Figure 8 (a) during an x-sweep, the boundary movement is isolated to the horizontal direction, maintaining an invariant projection of the solid surface (between adjacent geometric vertices) onto the vertical cell interfaces, i.e., . The same argument holds for the projection of the same solid-fluid interface section on the cell interface during the y-sweep in Figure 8 (b). The change in the cut-cell face area fractions, , and the average cut-cell distance, , over a time step then becomes a function of the angle between the cell interface and the local solid-fluid interface (i.e., the surface slope in the current piece-wise planar surface representation), together with the local wall velocity component in the direction of the current sweep.

For the numerical examples shown in this paper, in which a non-deforming rigid body is defined, we have made the assumption that the time step, , is small enough such that the time averaged quantities can be approximated using a linear movement of in both and sweeps. It is recognised that this is not a good approximation in general, and an improvement to this assumption to define an efficient general method for describing the time averaged cut-cell area fractions and length fractions is under development. Falcovitz et al. falcovitz1997two () describes one such method.

Using the time averaged and , the modified version of the shielded flux is therefore given by:




The total flux at the boundary, , is then given by:


where the unshielded flux, is the regular Cartesian cell flux given by Equation 5.

In the current dimensional-split formulation, the time-averaged boundary contribution in the , and directions are applied independently. The final form can therefore be divided into a boundary movement in the direction, followed by a boundary movement in the direction, followed by a boundary movement in the direction, with the corresponding component of the time-averaged boundary movement applied at each stage. The three-dimensional form of the moving boundary conservation equations are therefore:


where, and are used for compactness. The three-dimensional computational cell orientation and solid region used in this definition are shown in Figure 9.

Figure 9: Three-dimensional computational cell orientation, with the solid region highlighted in grey. The remaining (white) space contains a single phase fluid.

To maintain a second-order accurate scheme, the x-y-z sequence is reversed from Equations 23-25 in the next sweep via Strang-splitting, away from the solid surface. A detailed introduction to Strang-splitting for second-order accurate multi-dimensional schemes is given by Leveque LeVequeBook2002 (). As explained in Leveque LeVequeBook2002 (), choosing to reverse the x-y-z sequence in the Strang-splitting approach can be a more efficient alternative than taking half time-steps.

Finally, it is recognised that in the limit as , i.e., for a stationary solid, these equations reduce to the static boundary version of Klein et al. Klein09 (), i.e., the two-dimensional form of Equations 14-15 reduce to Equations 9-10 respectively.

2.3.5 Covering and Uncovering Cut-Cells

During the transit of the solid-fluid interface through the computational domain, computational cells can change over the period of a time step from regular Cartesian cells to cut-cells, from cut-cells to solid cells, or the reverse of these two. The change from a regular cell to a cut-cell is dealt with implicitly by the moving boundary method outlined in the preceding sections, as summarised in Equations 20-25. The uncovering of a cut-cell (transition from a solid cell to a cut-cell) and the ‘covering’ of a cut-cell (transition from a cut-cell to a solid cell), however, require additional measures to maintain strict conservation.

The change from a solid cell to a cut-cell (uncovering) is shown diagrammatically in Figure 10 (a). The reverse situation, of a cut-cell changing to a solid cell, is shown in Figure 10 (b).

(a) (b)

Figure 10: The covering and uncovering of a cut-cell due to the solid-fluid interface moving from the dashed line at time level to the solid line at time level . (a) An uncovered cut-cell at . (b) A covered cut-cell at .

Solid cells in the context of the current method are not updated and therefore have no valid conservative variable state. The change over time from a solid cell to a cut-cell therefore requires initialisation of the cut-cell before continuing the cut-cell update procedure.

To maintain strict global conservation, an initial state is computed from a volume weighted average of the neighbouring cells, with the contribution from each neighbouring cell being reduced from its own conserved state. Following the method outlined in Schneiders et al schneiders2013accurate (), we define a volume weighting, , for a cell with shorthand index, , as:


where is the total number of neighbouring cells used as donors for the newly uncovered cell. This includes the uncovered cell itself. In this work, donor cells are defined as those neighbouring cells horizontal or vertical to the uncovered cell. For a two-dimensional computation this may include donor cells at , , , and the uncovered cell, . The uncovered cell, in Figure 10, is then initialised using a conserved variable state of:


where is the total number of donor cells (excluding the uncovered cell). Note that for a solid cell and for fluid cells. This ensures only cut-cells or regular fluid cells contribute to the new cell initialisation. Each of the donor cells is then updated with a modified state according to:


where subscripts and distinguish the donor and uncovered cell respectively for the general case. The penalty of maintaining conservation in this way is a local reduction in the primitive variables for those donor cells surrounding the uncovered cell. The influence of this penalty can be reduced through increased use of Adaptive Mesh Refinement (AMR), a method of locally refining the computational mesh around flow features of interest using a hierarchical mesh configuration based on a user defined mesh refinement criterion. The density is used as the refinement criterion in all the test cases described in this paper that use AMR, employing a user-defined refinement threshold which is individually tailored for each test case to maintain a tight confinement of the refined mesh around the compressible flow features, while maintaining efficiency in preventing unnecessary sub-division of the refined region into excessively small AMR blocks.

To attempt to maintain global conservation, the reverse situation of a cut-cell being covered completely by the solid, requires the conserved variables in the cut-cell, , to be completely expelled through the fluid portion of the inter-cell boundary. The cut-cell conserved variables are gradually passed to the neighbouring cells through the shielded flux as the cut-cell volume fraction . To ensure strict conservation is maintained after the cut-cell is completely covered, any residual fluid remaining in the cut-cell is distributed to the neighbouring cells. The residual conserved state remaining in the covered cut-cell is computed as:


This is distributed to the surrounding fluid cells as:


where, is the total number of recipient cells surrounding the covered cut-cell and subscript denotes the recipient cell index.

In the next section, the moving boundary cut-cell method described here is applied to a range of compressible flows involving the time dependent interaction of a moving rigid body with a single phase fluid. To compute the regular cell fluxes, and those of the cut-cell unshielded flux component, the MUSCL-Hancock solver (see, for example, Toro Toro ()) is used in combination with the HLLC Riemann solver and van-Leer slope limiter. Currently, complex geometries can be imported from the standard STL file format, allowing simulation of Computer Aided Design (CAD) models. However, as the test cases geometries in this paper are simple geometries, each geometry was explicitly specified using simple polygons, or the NACA profile.

3 Numerical Results

3.1 The Sloped Channel Flow - Moving Boundary Validation

We consider as our first test case the propagation of a simple perturbation in density within an inclined channel. The simulation configuration for the simple wave problem is shown in Figure 11. Two parallel interfaces at an angle of to the horizontal plane form an embedded boundary, creating cut-cells with a large range of volume fractions. The initial condition consists of a simple density wave, centred at and travelling in a direction parallel to the channel walls. The channel itself moves in the positive x-direction over time. This test case is intended to assess how the moving boundary approach handles reducing and increasing cut-cell volume fractions, as well as the covering and uncovering of cut-cells. By eliminating the force of the cut-cell boundary movement on the fluid, we therefore effectively simulate a static channel test case resolved by a mesh moving in the negative horizontal direction. Conservation is maintained through an identical movement of the two channel walls. This is confirmed during the computation. The solid wall boundary condition is assessed for the moving boundary cut-cell method in subsequent test cases. The CFL number used in this test case is 0.5. Figure 12 shows the relative orientation of the Cartesian mesh, the channel boundaries and the initial density perturbation, which is shown as a grey scale pseudo-colour plot.

Figure 11: Sloping channel flow configuration. The test geometry consists of a channel oriented an angle of to the horizontal boundary. The channel moves in the positive horizontal direction over time.
Figure 12: Sloping channel case mesh showing the relative orientation of the Cartesian mesh, the channel boundaries and the initial density perturbation. The density perturbation is shown using a grey scale pseudo-colour plot.

The simulation is performed on a domain m by m, with the width of the channel set to m. The final time for each simulation is chosen as ms. This allows the density pulse to propagate across the majority of the domain, and the entire cut-cell surface to experience the full range of cut-cell volume fractions at each point.

The density along the lower channel surface at a time of ms is given in figure 13 for mesh resolutions of , and cells. This figure shows a comparison of the moving boundary results against the exact analytical density profile and against an equivalently resolved static boundary simulation using the current solver, applied to a static channel.

Figure 13: Density profiles along the lower channel surface at ms. Comparison of the moving boundary method with a static boundary method for a static channel. Comparison with an exact analytical density profile. (triangular symbols), (square symbols), (circles). Static cut-cell case (filled symbols), Moving cut-cell case (hollow symbols). The black line shows the exact analytical density profile at ms.

A satisfactory comparison between the moving boundary and static boundary results is shown for the three mesh resolutions given. Numerical diffusion of the cut-cell method is shown at the outer edges of the density pulse. This, of course, varies with mesh resolution, with the lowest diffusion evident in the results. Numerical dissipation is also evident at the centre of the pulse resulting in a reduction in the density pulse peak value.

The cut-cell error with increasing mesh resolution is shown for the static boundary results at mesh resolutions of , , , and in table 1. This is compared with the moving boundary cut-cell error in table 2.

Resolution Cell size norm order norm order norm order
2.5(-3)m 4.74(-2) - 1.03(-1) - 5.99(-1) -
1.25(-3)m 1.79(-2) 1.40 5.32(-2) 9.53(-1) 4.19(-1) 5.16(-1)
6.25(-4)m 5.61(-3) 1.67 2.64(-2) 1.01 2.53(-1) 7.28(-1)
3.12(-4)m 1.69(-3) 1.73 1.01(-2) 1.39 1.21(-1) 1.06
1.56(-4)m 4.61(-4) 1.87 3.82(-3) 1.40 5.87(-2) 1.07
Table 1: Static boundary results. Global , and error norms for the propagation of a density perturbation in a static boundary sloped channel at ms. Values given in compact form as to represent .
Resolution Cell size norm order norm order norm order
2.5(-3)m 6.84(-2) - 2.11(-1) - 8.82(-1) -
1.25(-3)m 3.01(-2) 1.18 1.07(-1) 9.80(-1) 7.11(-1) 3.11(-1)
6.25(-4)m 1.07(-2) 1.49 5.37(-2) 9.95(-1) 5.38(-1) 4.02(-1)
3.12(-4)m 3.32(-3) 1.69 2.51(-2) 1.10 2.99(-1) 8.47(-1)
1.56(-4)m 9.50(-4) 1.81 9.91(-3) 1.34 1.51(-1) 9.86(-1)
Table 2: Moving boundary results. Global , and error norms for the propagation of a density perturbation in a moving boundary sloped channel at ms. Values given in compact form as to represent .

The and norm error order indicate that the global order of accuracy increases towards second order as the mesh resolution increases. This is evident for the static boundary test case in table 1 and the moving boundary test case in table 2. The cut-cell method is locally first order accurate. The global error increases above unity as the influence of the first order boundary region diminishes with an increasing proportion of second order accurate regular cells at higher mesh resolutions. The norm error shows a smaller order of accuracy, increasing to around unity for both the static and moving boundary cases. The norm error uses the maximum error in the domain. This is therefore limited by the first order accurate cut-cell boundary region. The increase in the global order of accuracy with increasing mesh resolution justifies the introduction of Adaptive Mesh Refinement (AMR) to resolve flow features of large density gradient and regions close to the cut-cell boundary, restricting the spatial influence of the first order accurate cut-cell method with the current second-order accurate numerical solver. Refinement of the mesh using AMR is test case specific and can be chosen to refine across, for example, shock waves, contact surfaces or vortices through adjustment of a refinement criterion. Computational cells are flagged for refinement based on a user defined state variable threshold value. Adjacent computational cells with a difference in the chosen state variable are flagged for refinement. In the current computation, the refinement variable is density. If desired, in the current implementation, the rigid body surface can also be independently flagged for refinement.
The analysis presented so far in this section studies the convergence properties of the current explicit scheme at fixed CFL number, Tables 1 and 2, yielding the global accuracy of the method in both time and space as the mesh resolution increases. In an extension to this analysis, we isolate the space and time resolution contributions by adjusting the CFL number as the mesh is refined. Suppose we have a scheme that is formally second order accurate in both space and time. A measure of the spatial accuracy in isolation can then be computed by allowing the CFL number to vary between mesh resolutions, according to:


where, is a positive integer representing the current refinement level. This value is defined as for the least refined mesh case and increases by one as the mesh resolution is doubled in all spatial directions. in Equation 31 therefore denotes the CFL number used in the least refined case. With the CFL number linked to the mesh refinement in this manner, the time step approaches zero more rapidly than the mesh resolution, therefore the spatial error will dominate in this limit. Using this strategy, the global error for the moving channel case is given in Table 3. The most refined three mesh resolutions from Table 2 have been tested here. Table 3 documents the CFL number as the mesh resolution increases.

Resolution Cell size CFL norm order norm order norm order
6.25(-4)m 0.5 9.27(-3) - 3.92(-2) - 4.62(-1) -
3.12(-4)m 0.25 3.08(-3) 1.59 1.92(-2) 1.03 2.64(-1) 8.07(-1)
1.56(-4)m 0.125 8.70(-4) 1.82 7.45(-3) 1.37 1.38(-1) 9.36(-1)
Table 3: Moving boundary results. Isolation of the spatial component of the global , and error norms for the propagation of a density perturbation in a moving boundary sloped channel at ms. Values given in compact form as to represent .

Table 3 indicates that the influence of the first order cut-cell method is restricted to a small region close to the surface in terms of the global order of accuracy when the mesh is sufficiently well resolved. When considered together with Table 2, which demonstrated a similar convergence to second order accuracy, the implication is that the accuracy of the second-order method is globally maintained in both time and space for the most refined mesh resolution cases.

3.2 Circular Cylinder Interaction with a Stationary Shock Wave

An assessment of the current moving boundary cut-cell method for computing the interaction of a curved surface with a moderately strong shock wave is presented in the second test case. A stationary shock wave initially separates two regions of fluid with a discontinuous rise in pressure and density, accompanied by a decrease in fluid velocity, in the downstream direction. The circular cylinder moves downstream with the same velocity as the initial low pressure region, eventually passing through the stationary shock wave. This models the experimental situation of a stationary cylinder and moving shock wave with the same rise in pressure and density across the shock wave. An ideal gas is used, with specific heat capacity ratio, , and specific gas constant, J/kgK. The CFL number is 0.5.

Figure 14 shows the interaction of a stationary Mach 2.0 shock wave with a moving cylinder of diameter 0.04m. Specifically, numerical Schlieren plots are used to highlight compressible flow features at times of 6s, 35s, 106s and 126s after the initial contact of the shock wave with the circular cylinder surface.

(a) (b)

(c) (d)

Figure 14: Numerical Schlieren plot series showing a circular cylinder passing through a stationary shock wave. (a) 6s, (b) 35s, (c)106s, (d) 126s. Times taken with respect to the initial contact of the cylinder and shock wave.

This Mach 2.0 shock wave test case is initialised by globally imposing the stationary shock wave flow properties throughout the computational domain. Specifically, the flow properties to the left and right-hand sides of the shock wave are 0.595, 459.54, 50000 and 0.944, 289.86, 96410. The circular cylinder velocity is 459.54, matching the fluid velocity on the left-hand side of the shock wave. The cylinder is initially located one diameter upstream of the shock wave. Before contact with the shock wave, therefore, the circular cylinder initially moves with the same velocity as the surrounding fluid. The absence of pressure waves in the low pressure fluid at this time demonstrates that the consistency condition is satisfied by the moving cut-cell method.

The location of the initial shock wave and the development of additional compressible flow features, such as a bow shock and shock wave triple points are evident in Figure 14. As outlined in Zdravkovich zdravkovich1997flow (), the location, strength and development of these compressible flow features are dependent on the magnitude of the pressure and density rise across the shock wave. Even moderate increases in the shock wave Mach number result in an altered shock wave pattern, providing a sound basis with which to assess the moving boundary cut-cell surface condition. Figure 15 shows a comparison between the current moving boundary cut-cell method and an experiment by Bryson & Gross bryson1961diffraction () at a higher shock wave Mach number of . At this Mach number, the pressure and density ratio across the shock wave is , and the velocity ratio is . The Cartesian mesh resolution used for Figure 15 is , where is the circular cylinder diameter.

Figure 15: Circular cylinder passing through a Mach 2.81 shock wave. Top: Schlieren image from the experiment of Bryson & Gross bryson1961diffraction () in which a moving shock wave passes over a fixed cylinder. Bottom: Numerical Schlieren image using the current moving boundary cut-cell method and modelling a stationary shock wave with a moving cylinder.

The distance from the circular cylinder to the bow shock, and the location of the shock triple point, in the numerical Schlieren compare well with the experiment of Bryson & Gross bryson1961diffraction (). The contact discontinuity, which runs between the circular cylinder and the shock triple point is also well located, indicating the moving cut-cell method is able to capture the correct advection dominated flow behaviour regarding the shock wave interaction with the moving circular cylinder.

3.3 Schardin Experiment

The interaction of a solid body with an incident shock wave is extended in this test case to the classic experiment of Schardin Schardin57 (), in which a planar shock interacts with a triangular wedge. This interaction results in the generation of a relatively complex flow pattern, with transonic tip vortices extending downstream of the wedge, as well as Mach stems, triple points, reflected and diffracted shocks, slip lines and acoustic waves. Schardin’s original shadowgraph images provide an initial basis for comparison, with further experimental imaging by Chang and Chang chang2000shock () allowing for more detailed verification. As with the circular cylinder test case, we model a stationary shock wave, matching the experimental pressure and density rise across the shock wave. A triangular wedge, moving with the same velocity as the surrounding fluid on one side of the shock wave, is then allowed to pass through the stationary shock wave. This models the experiment, in which a fixed stationary wedge in static fluid is subjected to the impingement of a moving shock wave. An ideal gas is used, with specific heat capacity ratio, , and specific gas constant, J/kgK. The CFL number used in this test case is 0.5.

The solid body in this test case takes the form of an equilateral triangle, with base length 20mm. The triangular wedge is placed in a fluid with a density and pressure of kg m and MPa. The incident shock wave has Mach number 1.34, giving a post-shock state of kg m and MPa.

The numerical results are computed on the domain of size