A dimensionally split Cartesian cut cell method for the compressible Navier-Stokes equations
AA dimensionally split Cartesian cut cell method for the compressibleNavier-Stokes equations
Nandan Gokhale a, ∗ , Nikos Nikiforakis a , Rupert Klein b a Laboratory for Scientific Computing, Cavendish Laboratory, University of Cambridge, Cambridge, CB3 0HE, UK b Institut f¨ur Mathematik, FB Mathematik und Informatik, Freie Universit¨at Berlin, Arnimallee 6, 14195 Berlin, Germany
Abstract
We present a dimensionally split method for computing solutions to the compressible Navier-Stokes equa-tions on Cartesian cut cell meshes. The method is globally second order accurate in the L norm, fullyconservative, and allows the use of time steps determined by the regular grid spacing. We provide a de-scription of the three-dimensional implementation of the method and evaluate its numerical performanceby computing solutions to a number of test problems ranging from the nearly incompressible to the highlycompressible flow regimes. All the computed results show good agreement with reference results from theory,experiment and previous numerical studies. To the best of our knowledge, this is the first presentation of adimensionally split cut cell method for the compressible Navier-Stokes equations in the literature. Keywords:
Cartesian grid, Cut cell, Dimensional splitting, Navier-Stokes, Adaptive Mesh Refinement,Immersed boundary method
1. Introduction
The Cartesian cut cell mesh generation procedure involves computationally ‘cutting out’ the geometryfrom a background Cartesian grid to produce a resulting mesh with a sharp representation of the interface.The procedure is attractive since it allows for rapid, automatic mesh generation for complex geometrieswhile retaining the computational conveniences offered by the use of Cartesian grids. Furthermore, cut cellmethods, unlike other types of immersed boundary methods [1], are designed to be fully conservative (seeBerger [2] for a comprehensive recent overview of cut cell methods). However, the resulting mesh can containcut cells of arbitrarily small volume fraction at the interface which impose severe constraints on the timestep for explicit numerical schemes.Since the early 1980s, a number of ways have been presented in the literature to deal with this ‘small cellproblem’ in the context of the Euler equations [3, 4, 5, 6, 7, 8]. Over the past decade, most of these techniqueshave been extended for solving the compressible Navier-Stokes equations. Hartmann et al. [9] have usedcell linking (which is related to the intuitive concept of cell merging) to develop a three-dimensional cut cellmethod that is implemented in an adaptive octree grid. Using the static boundary cut cell formulation ofHartmann et al. [9], Schneiders et al. [10, 11] have successfully developed a method to compute movingboundary problems in 3D by introducing an interpolation routine and flux redistribution step. Berger etal. [12] use pseudo-time stepping and a multigrid approach to compute steady state solutions for the 2DReynolds-averaged Navier-Stokes (RANS) equations. In a subsequent publication [13], they describe a cutcell implementation of a novel ODE-based wall model which, unlike conventional equilibrium wall functions,has the advantage that it can be applied further away from the interface in the wake region of a turbulentboundary layer. Graves et al. [14] extend the ‘flux redistribution’ technique of Colella et al. [7] for small cell ∗ Corresponding author
Email addresses: [email protected] (Nandan Gokhale), [email protected] (Nikos Nikiforakis), [email protected] (Rupert Klein)
Preprint submitted to Journal of Computational Physics September 18, 2018 a r X i v : . [ phy s i c s . c o m p - ph ] S e p tability to develop a second order accurate method. Another high order discretisation was demonstratedby Muralidharan and Menon [15], who extended the flux mixing technique of Hu et al. [8] to develop a thirdorder accurate scheme.The aforementioned techniques are all implemented in an unsplit fashion. We are particularly interestedin adopting a dimensionally split approach which is a convenient way to extend one-dimensional methods tosolve multi-dimensional problems. In that context, Gokhale, Nikiforakis and Klein [16] recently presented asimple dimensionally split cut cell method for hyperbolic conservation laws using a ‘Local Proportional FluxStabilisation’ (LPFS) approach and demonstrated its performance through the computation of solutionsto a number of challenging problems for the Euler equations. The LPFS method is an improvement onthe original split method of Klein, Bates and Nikiforakis (KBN) [17]. Although both methods are firstorder accurate at the interface, LPFS was shown to produce more accurate solutions near boundaries forhyperbolic problems, and it allows the use of larger Courant numbers [16]. In this paper, we combine andextend the LPFS and KBN methods to solve compressible Navier-Stokes problems involving rigid embeddedboundaries. To the best of our knowledge, this is the first presentation of a dimensionally split cut cellmethod for the compressible Navier-Stokes equations in the literature.The rest of this paper is organised as follows. In Section 2, we outline the governing equations andsolution framework that we use. A description of the cut cell geometric parameters required by the methodis provided in Section 3. In Section 4, we describe the numerical method in detail. In Section 5, we presentnumerical solutions for a number of multi-dimensional test problems to demonstrate the performance of themethod. Finally, conclusions and areas for future work are provided in Section 6.
2. Governing equations and solution framework
The compressible Navier-Stokes equations are ∂ t ρ + ∇ · ( ρ u ) = 0 ,∂ t ( ρ u ) + ∇ · ( ρ u ⊗ u + pI ) = ∇ · σ, (1) ∂ t E + ∇ · [( E + p ) u ] = ∇ · ( σ u ) + ∇ · ( ξ ( ∇ T )) , where ρ is density, u is velocity, p is pressure, I is the identity matrix, ξ is thermal conductivity and T istemperature. We assume the fluid under consideration is Newtonian such that σ , the stress tensor, is σ = µ ( ∇ u + ∇ u T ) + λ ( ∇ · u ) I, (2)where µ is the dynamic viscosity and we use the Stokes’ hypothesis λ = − µ . E is the total energy per unitvolume, given by E = ρ (cid:18) | u | + e (cid:19) , (3)where e is the specific internal energy. To close the system of equations Eq. 1-Eq. 3 we use the ideal gasequation of state e = pρ ( γ − , (4)where γ , the heat capacity ratio, is assumed to be 1.4.We defer our description of the procedure for computing the explicit inviscid (Godunov-based) andviscous fluxes till Section 4.1, although it may be noted that the flux stabilisation approach is independentof the particular choice of flux methods used. Hierarchical Adaptive Mesh Refinement (AMR) [18] is usedto refine areas of interest such as the cut cell interface, or shock waves, while allowing the use of coarserresolutions elsewhere for the sake of computational efficiency. In this work, we manually specify the extentof the refined region for all subsonic problems. For supersonic problems, dynamic refinement of the regioncovering the shock is carried out by comparing a density gradient based metric to a user-defined threshold.Multi-dimensional updates are performed using Strang splitting [19] in 2D, and straightforward Godunov2plitting [20] in 3D, although Strang splitting could also be used in 3D if time order of accuracy was importantfor the problem at hand.As in Dragojlovic et al. [21], we assume that the global time step, ∆ t , is restricted by the minimum ofthe hyperbolic and diffusive time steps∆ t = min [∆ t hyp , ∆ t diff ] , = min (cid:34) C cfl min d,i (cid:32) ∆ x d,i W max d,i (cid:33) , min d,i (cid:32) ∆ x d,i µ i ρ i , ξ i ( ρc p ) i ) (cid:33)(cid:35) , (5)where d is the index of the coordinate direction, i is the index of a computational cell, and c p is the specificheat at constant pressure. ∆ x d,i and W max d,i are the spatial resolution and max wave speed for cell i inthe d direction respectively. As in similar previous works [21, 9, 10, 15], we operate the algorithm in theregion where the hyperbolic condition is most restrictive, and all simulations in this paper were run using aCourant number of 0.5.The wave speed for cell i in the d direction, W d,i , is computed using the following estimate suggested byToro [22]: W d,i = | u d,i | + a i , (6)where u d,i is the component of the velocity in cell i in the d direction. a i is the speed of sound in cell i ,given by a i = (cid:114) γp i ρ i . (7)
3. Mesh generation I I I I A BC DE FG H
FluidSolid xyz
Figure 1: Illustration of a cut cell in 3D.
Consider Fig. 1, which shows a 3D cut cell where the solid-fluid interface intersects the cell at four points I , I , I and I . The geometric parameters required by the method are the following:(i) The face fraction, β ∈ [0 , A b , of the reconstructed interface in the cell and ˆn b , the interface unit normal. The super-script b is short for ‘boundary’. 3iii) The volume fraction, α ∈ [0 , x c , of the fluid part of the cell.(v) The interface centroid, x bc of the reconstructed solid interface.We proceed by treating the interface implicitly as the zero level-set of a signed distance function. Thetechnique of Mauch [23] is used to compute the signed distance function, φ ( x ), at the vertices of the cell. Asdescribed in [16], this information can be used to work out the coordinates of the intersection points I - I and parameters (i)-(iii) defined above. We avoid repeating the details here for the sake of brevity. It maybe noted that our mesh generation procedure is only valid when the intersection of the cell and geometrycan be described by a single interface. Sufficient resolution must therefore be used to ensure that all cutcells are ‘singly cut’. The signed distance function can be used to deal with ‘split cells’ created by multipleintersections of the geometry with the cell by using a ‘Marching Cubes’ [24] based approach as in Guntheret al. [25], but this is beyond the scope of this work. Next, we describe the calculation of parameters (iv)and (v) in two and three dimensions.In 2D, the fluid part of a cut cell is a polygon. x c = [ x c,x , x c,y ] T can therefore be worked out using theformula for the centroid of a non-self-intersecting polygon with n ordered vertices ( x , y ) , ..., ( x n , y n ) [26] x c,x = 16 A n (cid:88) i =1 ( x i + x i +1 )( x i y i +1 − x i +1 y i ) , (8) x c,y = 16 A n (cid:88) i =1 ( y i + y i +1 )( x i y i +1 − x i +1 y i ) , (9)where A is the area of the polygon [26] A = 12 n (cid:88) i =1 ( x i y i +1 − x i +1 y i ) . (10)Note that in Eq. 8-Eq. 10, the summation index is periodic so that ( n + 1) = 1.Since we assume that cut cells are singly cut, x bc in 2D is worked out trivially as the average of thepositions of the two intersection points.In 3D, the fluid part of the cell is a polyhedron. x c can therefore be calculated using the formula for thecentroid of an arbitrary polyhedron with N F faces [27] x c = 34 (cid:104)(cid:80) N F i =1 ( x c , i · ˆn i ) x c , i A i (cid:105)(cid:104)(cid:80) N F i =1 ( x c , i · ˆn i ) A i (cid:105) , (11)where x c , i is the centroid of face i (computed by the appropriate use of Eq. 8-Eq. 10), and ˆn i and A i arethe outward unit normal and area of face i respectively.To compute x bc in 3D, we follow the following procedure:1. Rotate the interface plane with n vertices ( x , y , z ) , ..., ( x n , y n , z n ) to give a polygon aligned withthe x - y Cartesian plane. This polygon has vertices ( x rot1 , y rot1 , z rot ) , ..., ( x rot n , y rot n , z rot ). Note thatbecause the interface reconstructed from the signed distance function may not be perfectly planar, the z coordinates of the rotated interface may differ slightly from one another. In practice, we set z rot tobe the numerical average of these varying z coordinates.2. Use Eq. 8-Eq. 10 on the polygon with ordered vertices ( x rot1 , y rot1 ) , ..., ( x rot n , y rot n ) to get x b, rot c,x and x b, rot c,y , which are the x and y coordinates of the interface centroid in the rotated frame, x b, rot c . Notethat x b, rot c = [ x b, rot c,x , x b, rot c,y , z rot ] T .3. Rotate x b, rot c back to the original frame to give x bc as required.4 . Numerical method In this subsection, we describe the procedure we use to compute the explicit inviscid and viscous fluxes.Their stabilisation at the cut cells is described in Section 4.2.
In this work, we compute the inviscid intercell fluxes that appear in the divergence terms to the left ofthe equal signs in Eq. 1 using an exact Riemann solver and the MUSCL-Hancock scheme in conjunctionwith the van-Leer limiter [22]. This scheme is second order accurate in smooth regions.The calculation of the viscous fluxes that appear in the divergence terms to the right of the equal signsin Eq. 1 requires the computation of ∇ u and ∇ T at the cell faces. In Fig. 2, we illustrate our procedurefor calculating ( ∇ φ ) i − / ,j for a general scalar φ at the left face of cell ( i, j ) in an x dimensional sweep intwo dimensions. i − i − / ij − jj + 1 (a) Regular region away fromthe interface. i − i − i − / ij − jj + 1 xy (b) Irregular region near the interface.Figure 2: Illustration of the calculation of viscous face derivatives in irregular and regular regions. In a regular region away from the interface, the required face derivatives can be worked out using centraldifferencing as illustrated in Fig. 2a. The derivative in the current coordinate direction is calculated tosecond order accuracy as (cid:18) ∂φ∂x (cid:19) i − / ,j = φ i − φ i − ∆ x . (12)The transverse derivative is computed to second order accuracy as the average of the neighbouring transversederivatives (cid:18) ∂φ∂y (cid:19) i − / ,j = 12 (cid:18) φ i − ,j +1 − φ i − ,j − y + φ i,j +1 − φ i,j − y (cid:19) . (13)It may be noted that in 3D, the face derivatives can be calculated analogously on a 2 × × ∇ φ ) LS i − ,j and ( ∇ φ ) LS i,j at the volumetric centroids of cells ( i − , j ) and ( i, j )respectively. Note that the stencil for the weighted least squares calculations contains all the regular andcut cells from the 8 (26 in 3D) neighbours of the cell. This is illustrated using dotted lines for cell ( i − , j )in Fig. 2b. ( ∇ φ ) i − / ,j is calculated as( ∇ φ ) i − / ,j = ( ∇ φ ) LS i − ,j + ( ∇ φ ) LS i,j . (14)5s with other cut cell discretisations [9, 12], our approximation of the face derivatives in irregular regions isonly first order accurate. This does not impact the overall accuracy of the method, however, since the fluxstabilisation is also first order accurate at the boundary. To ensure conservation in a dimensionally split scheme, advective boundary fluxes have to be treateddifferently to the pressure and diffusive fluxes. We describe the treatment for the inviscid Euler equationsin [16]. Essentially, the advective boundary flux in a cell needs to be evaluated using a ‘reference’ boundarystate computed by solving a wall normal Riemann problem. The reference state should be computed beforethe start of a time step and kept constant for the duration of the time step. The velocity resulting fromthe wall Riemann problem is tangential to the wall. As a consequence, the sum of advective fluxes acrossthe interface accummulated over a full cycle of split steps cancels exactly, thus ensuring zero net mass fluxacross the wall. The boundary pressure, on the other hand, is to be updated in between sweeps for the sakeof accuracy. For the Navier-Stokes equations, we need to compute also the stress tensor at the boundaryand update it in between sweeps in order to compute accurate viscous momentum boundary fluxes. Sincewe consider only static (no-slip), adiabatic boundaries in this work, note that there are no viscous heatingor thermal conductivity boundary fluxes to consider for the energy equation.To compute ∇ u at the boundary, we follow the process illustrated for cut cell ( i, j ) in Fig. 3. Note that inthe rest of this section, subscripts specified with greek letters are assumed to range from 1 to n d , the numberof dimensions, and repeated indices of that kind imply the use of the Einstein summation convention. Let x µ represent a Cartesian coordinate system, and ˆ x µ represent an orthonormal coordinate system with a unitvector pointing normal to the cell interface, and unit vector(s) in the interface tangential plane. i − ijj + 1 x interp h ˆ x ˆ x xy Figure 3: Illustration of the procedure used to compute ∇ u at the boundary of cell ( i, j ). u ∗ µ at x interp is reconstructed fromthe weighted least squares gradient computed at the nearest cell volumetric centroid (filled white circle). We start by computing the normal boundary derivatives for each velocity component u µ as (cid:18) ∂u µ ∂ ˆ x (cid:19) = u ∗ µ − u bµ h = u ∗ µ h , (15)where u bµ , the value of the velocity component at the boundary is 0 as required by the no-slip boundarycondition. u ∗ µ is the value of the velocity component at the interpolation point x interp which is located ata distance of h from the interface centroid in the normal direction. We reconstruct the value of u ∗ µ using aweighted least squares gradient computed at the cell volumetric centroid closest to x interp . h is calculatedas in Meyer et al. [28] using h = 0 . (cid:113) ˆn bν ∆ x ν , (16)where ∆ x ν is the spatial resolution in the ν coordinate direction. Note that when using a uniform meshspacing of ∆ x in all coordinate directions, h = 0 . x .6ith the normal derivatives calculated, and the tangential boundary derivatives known to be 0 becauseof no-slip and no space dependent wall motion, we can construct the tensor ∂u µ /∂ ˆ x ξ at the boundary.However, we want to compute the tensor ∇ u which has Cartesian components ∂u µ /∂x ν that are given by ∂u µ ∂x ν = ∂ ˆ x ξ ∂x ν ∂u µ ∂ ˆ x ξ . (17) ∂ ˆ x ξ /∂x ν can be represented by a matrix whose columns are the unit vectors of the ˆ x µ coordinate system.The normal vector ˆn b is already known from the information provided by the signed distance function (seeSection 3), while the vectors spanning the tangential plane are calculated from ˆn b using a Gram-Schmidtorthogonalisation process. Eq. 17 can then be used to compute ∇ u , and hence the stress tensor σ at theboundary as required to compute the momentum diffusion boundary fluxes.Finally, although beyond the scope of this work, it may be noted that if needed, it is possible to introducethe use of a simple algebraic or ODE-based turbulent wall function into the above approach [12, 29, 30]. Thereconstructed state at x interp and the known state at the wall interface centroid can be used as boundaryconditions to calculate the wall shear stress. ∆ x α ∆ x ∆ x F ni − / F ni − / F b,n i − i (a) KBN. i − i ∆ x α ∆ xα ∆ x ∆ x ∆ t ∆ t hypcc t n t n +1 t x (b) LPFS.Figure 4: Illustration of the KBN and LPFS flux stabilisation procedures for a boundary cut cell neighbouring a regular cellin 1D. In this subsection, we describe the combination and extension of the KBN and LPFS methods for theNavier-Stokes equations. Fig. 4a and Fig. 4b illustrate the one-dimensional flux stabilisation proceduresused in both approaches for a boundary cut cell i neighbouring a regular cell.Let U ni represent the conserved variable state vector for cell i at time level n , and let F ni ± i/ represent theexplicit numerical fluxes (inviscid and viscous) computed at its ends. As described in [16, 17], the stabilisedKBN flux can be derived to be F KBN ,ni − / = F b,n + α ( F ni − i/ − F b,n ) , (18)where we note that the expression is consistent with respect to the natural limits of the grid so that as α → F KBN ,ni − / → F ni − i/ , and as α → F KBN ,ni − / → F b,n .The KBN flux Eq. 18 uses only the geometric parameter α to derive a stable flux. The LPFS approachcombines both geometric and wave speed information to derive a flux that was shown in [16] to producemore accurate boundary solutions for hyperbolic problems. Consider Fig. 4b, which shows the boundary7ut cell neighbouring the regular cell in the x - t plane for one time step. ∆ t is the global stable hyperbolictime step which is determined in part by the fastest wave speed in the domain, W max (see Eq. 5). Forthe configuration of Fig. 4b, we illustrate the ‘small cell problem’ at the cut cell as being caused by theleft-going wave from the solution of the boundary Riemann problem. Stability would therefore require theuse of the smaller ∆ t hypcc ∆ t hypcc = C cfl α ∆ xW i , (19)where W i is the wave speed for the cut cell.In the LPFS approach, the basic idea is to use the explicit flux F ni − i/ for the part of the time stepfor which it is stable, ∆ t hypcc , and to employ a different flux which can maintain stability for the duration(∆ t hypcc , ∆ t ]. This gives the LPFS method an inherent advantage in regions of low velocity - although α < t hypcc relative to ∆ t , part of the reduction is offset if W i < W max , an effect which is mostpronounced in regions of low velocity near a stagnation point or in the boundary layer. For larger cut cells,the ratio ∆ t hypcc / ∆ t can be greater than 1 so that no flux stabilisation is required.As described in [16], a suitable LPFS flux is F LPFS ,ni − / = ∆ t hypcc ∆ t F ni − / + (cid:18) − ∆ t hypcc ∆ t (cid:19) F KBN,mod ,ni − / , (20)where F KBN,mod ,ni − / = F b,n + αα ( F ni − i/ − F b,n ) , (21)and αα = ∆ t hypcc ∆ t = (cid:15) αW max W i . (22)The ‘wave speeds uncertainty’ parameter (cid:15) ∈ [0 ,
1] is introduced to account for any errors arising from theuse of Eq. 6 to estimate the cut cell wave speeds. We found setting (cid:15) to 0.8 to be a robust choice for thewide range of problems tackled in this paper.The derivation of the LPFS flux Eq. 20 implicitly assumes that the cut cell time step is limited by thelocal hyperbolic time step ∆ t hypcc . Although the global time step for all problems in this work is indeedlimited by the hyperbolic time step, it is possible for the local time step at some cut cells to be limited bythe local diffusion time step ∆ t diffcc = ( α ∆ x ) µ i ρ i , ξ i ( ρc p ) i ) , (23)in which case one can use the KBN stabilisation Eq. 18 but not the LPFS stabilisation.For a completely stable discretisation, then, we compare the relative magnitudes of ∆ t hypcc and ∆ t diffcc before deciding how to stabilise the flux. For the usual case ∆ t hypcc ≤ ∆ t diffcc , we use the LPFS flux Eq. 20 atthe cut cell cell. If ∆ t diffcc < ∆ t hypcc , on the other hand, we use the KBN flux Eq. 18.8 .3. Multi-dimensional extension cell L cell R β C = β US + β SS,R β US β SS,R β L β R α shielded R (a) Geometric parameters. cell L cell R F US ,n F SS,R ,n F B ,n R (b) Intercell flux components.Figure 5: Illustration of the parameters used in the flux stabilisation process for multi-dimensional simulations. As described in detail in [17, 16], the extension of the 1D LPFS and KBN flux stabilisation approachesto multiple dimensions requires attention to be given to the irregular nature of the cut cells. With referenceto Fig. 5, we briefly summarise the procedure for computing the intercell flux between neighbouring cellsfor an arbitrary dimensional sweep in 2D. The process is exactly the same in 3D.The basic idea is to divide the interface into ‘unshielded’ (US) regions which do not ‘face’ the boundary,and ‘shielded’ regions which are ‘covered’ by the boundary in the current sweep direction. To minimisethe diffusive impact of the first order flux stabilisation, the aim is to deviate as little as possible from thestandard explicit update for regular cells. The flux on the unshielded region does not call for stabilisation sothe regular explicit flux (calculated as per Section 4.1.1) is applied there. The stabilised flux is applied onlyon the shielded regions. In Fig. 5, part of the interface is ‘singly-shielded from the right’ (SS,R) and so theflux there would be stabilised, with the place of α in Eq. 18 or Eq. 20 being taken by α shielded R , which is theaverage distance from the cell face to the boundary in the current coordinate direction, non-dimensionalisedby the corresponding regular cell spacing. The boundary flux F B ,n R is calculated as per Section 4.1.2. Notethat all the geometric parameters of Fig. 5a can be computed from the geometric information described inSection 3. The details are left out for the sake of brevity.The single modified flux F modified ,n C to be applied at the cell face, then, is taken as an area-weighted sumof the individual components: F modified ,n C = 1 β C (cid:2) β US F US ,n + β SS,R F SS,R ,n (cid:3) . (24)With the modified fluxes computed, the multi-dimensional update for a cell ( i, j ) using straightforwardGodunov splitting, for example, would be given by U n +1 / i,j = U ni,j + ∆ tα i,j ∆ x (cid:104) β i − / ,j F modified ,ni − / ,j − β i +1 / ,j F modified ,ni +1 / ,j − (cid:0) β i − / ,j − β i +1 / ,j (cid:1) F B ,ni,j (cid:105) , U n +1 i,j = U n +1 / i,j + ∆ tα i,j ∆ y (cid:104) β i,j − / G modified ,ni,j − / − β i,j +1 / G modified ,ni,j +1 / − (cid:0) β i,j − / − β i,j +1 / (cid:1) G B ,ni,j (cid:105) , (25)where we use F and G to denote fluxes acting in the x and y directions respectively. A 3D simulation wouldinvolve one more sweep using fluxes H acting in the z direction.For concave geometries, it is quite possible for part of the interface to be shielded by the boundary fromboth the left and right sides. Devising a completely stable flux to be applied in these ‘doubly-shielded’ regionsis a challenging research problem. A simple practical solution is to use the ‘mixing flux’ and ‘conservativecorrection’ procedures as described in [16]. 9 . Results Re = 20 lid-driven cavity problem The lid-driven cavity problem from Kirkpatrick et al. [31] was used to verify the accuracy of the numericaldiscretisation. Fig. 6a illustrates the simulation set-up. The left, right and bottom domain edges are staticno-slip boundaries. The top boundary is also no-slip, but with a parabolic x velocity profile u lid which variesfrom 0 at the boundary edges to u maxlid at the centre. u maxlid was set to correspond to M = 0 .
1. The Reynoldsnumber of the flow based on the lid length L and u maxlid is 20.Fig. 6a shows the contours of normalised x velocity for the ‘reference’ solution computed on a fine400 ×
400 grid, on which the minimum encountered cut cell volume fraction was 2 . × − . The simulationwas performed at five coarser resolutions, with the errors for these runs computed relative to the referencesolution. The dashed box in Fig. 6a shows the region used for the error computations. In [16], we showedthat for a scheme which is first order accurate at the cut cells and second order accurate elsewhere, the L p norm of the global solution error converges as O (∆ x p +1 p ). As seen in Fig. 6b, the computed solutionconverges with first order at the cut cells, while the global error measured by the L norm converges withsecond order accuracy as expected. LL L/ u lid − . . (a) u/u maxlid contours for the 400 ×
400 cells referencesolution. . . .
11 10 100 L e rr o r N First orderCut cell errorGlobal errorSecond order (b) Convergence of cut cell and global L error norms for velocitymagnitude. N is the number of cells along one coordinate direction.Figure 6: Velocity results for the Re = 20 lid-driven cavity problem. This test involves the computation of a two-dimensional flat plate boundary layer for M ∞ = 0 . Re L = 30000. We run the test in configurations with the plate both coordinate aligned and non-aligned asillustrated in Fig. 7. AMR is used to refine the no-slip region of the plate. Note that a uniform velocity10rofile is specified at the inflow boundary so that the boundary layer develops on the no-slip part of theplate. LRe L (a) Grid-aligned plate. LRe L (b) Plate aligned at 5 ◦ to the horizontal.Figure 7: Grid aligned and non-aligned configurations for the laminar flat plate boundary layer problem. The shaded regionof length L is the no-slip part of the plate. The remainder of the plate is specified as a slip boundary. For the grid-aligned configuration, the first row of finite volumes adjacent to the plate are all cut cellswith the same volume fraction. Fig. 8a shows a comparison of the computed boundary layer profile atthe centre of the plate for a series of 4 simulations with progressively smaller cut cell volume fractions.Four levels of AMR are employed such that on the finest level, there are roughly 30 cells resolving the 99%boundary layer thickness δ . All the computed profiles show good agreement with the theoretical Blasiussolution, which is a similarity solution for a steady, two-dimensional laminar boundary layer forming oversemi-infinite plate in an incompressible flow [32].We use twice the resolution for the non-aligned plate configuration such that there are roughly 60 cellsresolving the δ thickness at the centre of the plate. Like Graves et al. [14], we compare the computedboundary layer velocity profiles along ‘wall normal rays’ emanating from every cut cell in the range 5000 ≤ Re x ≤ . . . .
81 0 2 4 6 8 u / u ∞ η No cut cells α = 0 . α = 0 . α = 0 . (a) Boundary layer velocity profiles at Re x = 15000 for vary-ing cut cell volume fractions of the lowest row for the grid-aligned configuration. . . . .
81 0 1 2 3 4 5 6 7 u / u ∞ η SimulationBlasius (b) Velocity profiles along wall-normal rays emanating fromevery cut cell in the range 5000 ≤ Re x ≤ ◦ .Figure 8: Computed boundary layer profiles for the (a) co-ordinate aligned, and (b) non-aligned, laminar flat plate boundarylayer problem. Note that η = y (cid:113) uνx . The next problem considered is that of two-dimensional flow over a circular cylinder at Re = 40 and Re = 100. Large amounts of experimental and numerical results exist in the literature for both cases. At Re = 40, the eddies behind the cylinder are attached and we can evaluate the accuracy and convergencerate of the computed surface solutions. At Re = 100, the wake is unstable and we can also compute thefrequency of the vortex shedding process.The simulation set-up is illustrated in Fig. 9a. The left boundary is subsonic inflow while subsonicoutflow boundary conditions are specified on the other domain edges. As illustrated in Fig. 9b, AMR isused to resolve the near wall solution and cylinder wake. M ∞ is set to 0.1. D D D D (a) Domain dimensions (cylinder not to scale). (b) AMR mesh in the vicinity of the cylinder.Figure 9: Domain dimensions and an illustration of the Adaptive Mesh Refinement used for the simulations of the flow over acircular cylinder. D on the finest AMR level. At thefinest resolution, the minimum encountered cut cell volume fraction is 6 . × − . The results from the finestresolution simulations are treated as the reference solutions to which errors from the coarser simulations arecompared.As shown in Table 1, the computed drag coefficient from the reference solution for the Re = 40 case isin very good agreement with previous experimental and numerical studies. Table 1: Comparison of the computed drag coefficient with previous experimental and numerical studies for the Re = 40cylinder flow problem. Study C D Present work (resolution: D/ ∆ x = 160) 1.57Tritton [33] (experiment) 1.57Tseng and Ferziger [34] (simulation, resolution: πD/ ∆ x ≈
72) 1.53Meyer et al. [28] (simulation, resolution: D/ ∆ x = 72) 1.56Fig. 10a shows the computed pressure distribution over the cylinder for the Re = 40 case. The resultscompare well with experimental measurements from Grove et al. [35] over the whole of the cylinder. Fur-thermore, as shown in Fig. 10b the computed drag coefficient shows the expected first order convergencerate to the reference solution. − − . . . C p θ SimulationExperiment (a) Comparison of the numerical and experimental [35] surface pres-sure distributions. . . E rr o r D/ ∆ xC D errorFirst order (b) Convergence plot of C D computation.Figure 10: Results of surface pressure distribution and drag coefficient convergence for the Re = 40 cylinder flow problem. At Re = 100, the simulation produces a periodic vortex shedding flow pattern (Fig. 11a shows in-stantaneous computed vorticity contours in the wake of the cylinder). Table 2 compares the computedtime-averaged drag coefficient C D and Strouhal number St with previous experimental and numerical stud-ies. The frequency of the vortex shedding is calculated from the Fourier Transform of the solution for y velocity at a point in the middle of the wake located a distance of 1 . D behind the cylinder. The simula-13ion results are in good agreement with the previous studies. As shown in Fig. 11b, the computed dragcoefficients also converge with first order as expected. Table 2: Comparison of the computed drag coefficient and Strouhal number with previous experimental and numerical studiesfor the Re = 100 cylinder flow problem. Study C D St Present work ( D/ ∆ x = 160) 1.41 0.165Relf [33] (experiment) 1.39 -Wieselsberger [36] (experiment) 1.40 -Williamson [37] (experiment) - 0.165Tseng and Ferziger [34] (simulation, resolution: πD/ ∆ x ≈
72) 1.42 0.165Lai and Peskin [37] (simulation, resolution: D/ ∆ x ≈
38) 1.45 0.165 (a) Instantaneous vorticity contours over-laid on the mesh in the vicinity of the cylin-der. . . E rr o r D/ ∆ xC D errorFirst order (b) Convergence plot of C D computation.Figure 11: Results of vorticity and drag coefficient convergence for the Re = 100 cylinder flow problem. For the compressible flow problem of the shock reflection from a wedge, we use the simulation parametersfrom Graves et al. [14] shown in Table 3. As illustrated in Fig. 12, the initial conditions are zero velocity witha discontinuity in pressure and density. The horizontal ground is also discretised with cut cells of volumefraction 10 − . A base resolution of 1024 ×
512 cells is employed with two AMR levels of refinement factor2 each to resolve the waves and cut cell interfaces. The minimum encountered cut cell volume fraction onthe wedge surface is 4 . × − . The ground and wedge are no-slip boundaries, while transmissive boundaryconditions are specified at the left, right and top domain edges.14 able 3: Initial conditions (with reference to Fig. 12) for the shock reflection from wedge problem. Variable Value p (Pa) 1 . × p (Pa) 7 . × ρ (kg / m ) 3 . × − ρ (kg / m ) 3 . × − µ (kg/(m s)) 1 . × − C v (J/(kg K)) 3 . × ξ (W/(m K)) 1 . × − γ / ◦ .
09 m0 . .
15 m0 .
075 m p , ρ p , ρ Figure 12: Simulation set-up for the shock reflection from wedge problem.
The solution evolves to form a left-propagating rarefaction and a right-propagating M = 7 . t = 9 . µ s. Note that the extra wave visible behind the shock in Fig. 13a is thenumerical ‘start-up error’ [38] caused by starting the simulation with sharply discontinuous initial conditions.15 S (a) L S is the distance travelled by the Mach stem up thewedge. L M (b) Close-up of the Mach stem of length L M .Figure 13: Density contour plots at t = 9 . µ s of the solution for the shock reflection from wedge problem. The parameter of interest is R M = L M L S , (26)the ratio of the Mach stem length to the distance travelled by the Mach stem up the wedge. L S is determinedfrom a line-out of the solution taken along the wedge surface. To avoid ambiguity in the determination of theposition of the triple point, we compute L M via an algebraic approach. Consider Fig. 14, which illustratesthe position of the incident shock (ET) and Mach stem (TD) at the end of the simulation. Once distancesBD ( L S ) and AC are known, TD ( L M ) and hence, R M can be calculated via trignometry. AC is determinedfrom another line-out of the solution taken horizontally through the incident shock. Assuming an error of ± ∆ x in the determination of the average positions of the incident shock and Mach stem from the lineouts,we calculate R M = 0 . ± . R M with previousnumerical studies. The results from all three studies agree to 1 significant figure - differences in measurementprocedures used in the studies are the most likely explanation for why closer agreement is not obtained. Inour results, for example, ∆ x/L M ≈ . R M is very sensitive to themethod used to measure L M . 49 ◦ ET DA B C
Figure 14: Illustration of the positions of the incident shock (ET) and Mach stem (TD) at the end of the shock reflection fromwedge simulation. able 4: Comparison of the computed value of R M with previous numerical studies for the shock reflection from wedge problem. Study R M Present work 0 . ± . The final test problem is that of three-dimensional, supersonic flow over a sphere for M ∞ = 2 . Re D = 6 . × . The sphere of diameter D = 0 . . D from the inlet in adomain of size 10 D × D × D . The free stream density ρ ∞ , velocity u ∞ and temperature T ∞ are set tobe 0.2199 kg / m , 527.2 m/s and 173 K respectively as in Uddin et al. [40]. A supersonic inflow boundarycondition is specified at the inlet, while transmissive boundary conditions are used at all the other domainboundaries. A base resolution of 100 × ×
50 cells is employed and three levels of AMR of refinement factor2 each are used to resolve the interface and sphere bow shock. The minimum encountered cut cell volumefraction is 2 × − . The simulation is run until a time of t = 100 D/u ∞ . Fig. 15 shows a plot of the AMRpatches overlaid on contours of normalised velocity magnitude ( | v | /u ∞ ) along the x - y symmetry plane ofthe sphere at the end of the simulation. Clearly visible in the figure are the bow shock ahead of the sphereand the viscous wake behind it.The two results parameters of interest to us to validate the solution are the shock standoff distance δ/D , and the drag coefficient C D . δ is the shortest distance from the shock to the sphere leading edge. Asseen in Table 5, our computed values for δ/D and C D show good agreement with previous experimentaland numerical studies. Note that using one less AMR refinement factor produced no change to the dragcoefficient computed to three significant figures, suggesting that our final resolution is sufficient to producea numerically converged solution. Table 5: Comparison of the computed shock standoff distance and drag coefficient with previous experimental and numericalstudies for the three-dimensional supersonic flow over a sphere problem.
Study δ/D C D Present work 0.175 1.01Krasil ' shchikov and Podobin [41] (experiment) 0.17 1.00Bailey and Hiatt [42] (experiment) - 1.00Uddin et al. [40] (simulation) 0.18 1.07Al-Marouf and Samtaney [39] (simulation) 0.183 0.9617.00 1.23 Figure 15: Plot of AMR patches on the second and third refinement level overlaid on a contour plot of | v | /u ∞ .
6. Conclusions
In this paper, we presented a dimensionally split Cartesian cut cell method for the compressible Navier-Stokes equations. The numerical performance of the scheme was investigated through a number of testproblems ranging from the nearly incompressible to the highly compressible flow regimes.The computation of a nearly incompressible laminar boundary layer over a flat plate in both horizontaland inclined configurations was used to demonstrate that the method can handle coordinate aligned and non-aligned interfaces. The numerical convergence of the method was verified explicitly by computing solutionsfor a Re = 20 lid-driven cavity problem, and for flow over a circular cylinder at Re = 40 and Re = 100.For the highly compressible problem of a M = 7 . Re = 6 . × . The computed shock standoff distance anddrag coefficient showed good agreement with previous experimental and numerical studies.It is relatively straightforward to extend this work to be applicable for problems in the turbulent flowregime by adopting an explicit or implicit Large Eddy Simulation approach combined with a wall functionto compute the wall shear stresses. We have made progress on implementing such a ‘wall-modelled LES’scheme in the context of the split method and hope to present this as part of a future publication. Forthe future, modifying the method to achieve second-order accuracy at the boundary for both inviscid andviscous flows would also be a very useful contribution. Acknowledgements
Nandan Gokhale acknowledges the Cambridge Commonwealth, European & International Trust for fi-nancially supporting his research, and thanks Lukas Wutschitz for useful discussions. Rupert Klein acknowl-edges support by Deutsche Forschungsgemeinschaft through the Collaborative Research Centers CRC 102918project C01) and CRC 1114 (project C01). The authors would like to thank Philip Blakely for providingsupport with the parallel AMR code.
ReferencesReferences [1] R. Mittal, G. Iaccarino, Immersed boundary methods, Annu. Rev. Fluid Mech. 37 (2005) 239–261.[2] M. Berger, Cut cells: Meshes and solvers, in: Handbook of Numerical Analysis, Vol. 18, Elsevier, 2017, pp. 1–22.[3] D. K. Clarke, H. Hassan, M. Salas, Euler calculations for multielement airfoils using cartesian grids, AIAA journal 24 (3)(1986) 353–358.[4] J. J. Quirk, An alternative to unstructured grids for computing gas dynamic flows around arbitrarily complex two-dimensional bodies, Computers & fluids 23 (1) (1994) 125–142.[5] C. Helzel, M. J. Berger, R. J. LeVeque, A high-resolution rotated grid method for conservation laws with embeddedgeometries, SIAM Journal on Scientific Computing 26 (3) (2005) 785–809.[6] M. Berger, C. Helzel, A simplified h-box method for embedded boundary grids, SIAM Journal on Scientific Computing34 (2) (2012) A861–A888.[7] P. Colella, D. T. Graves, B. J. Keen, D. Modiano, A Cartesian grid embedded boundary method for hyperbolic conservationlaws, Journal of Computational Physics 211 (1) (2006) 347–366.[8] X. Hu, B. Khoo, N. A. Adams, F. Huang, A conservative interface method for compressible flows, Journal of ComputationalPhysics 219 (2) (2006) 553–578.[9] D. Hartmann, M. Meinke, W. Schr¨oder, A strictly conservative Cartesian cut-cell method for compressible viscous flowson adaptive grids, Computer Methods in Applied Mechanics and Engineering 200 (9) (2011) 1038–1052.[10] L. Schneiders, D. Hartmann, M. Meinke, W. Schr¨oder, An accurate moving boundary formulation in cut-cell methods,Journal of Computational Physics 235 (2013) 786–809.[11] L. Schneiders, C. G¨unther, M. Meinke, W. Schr¨oder, An efficient conservative cut-cell method for rigid bodies interactingwith viscous compressible flows, Journal of Computational Physics 311 (2016) 62–86.[12] M. Berger, M. J. Aftosmis, S. R. Allmaras, Progress Towards a Cartesian Cut-Cell Method for Viscous CompressibleFlow, Aerospace Sciences Meetings, American Institute of Aeronautics and Astronautics, 2012.[13] M. J. Berger, M. J. Aftosmis, An ODE-based Wall Model for Turbulent Flow Simulations, AIAA Journal (2017) 1–15.[14] D. Graves, P. Colella, D. Modiano, J. Johnson, B. Sjogreen, X. Gao, A cartesian grid embedded boundary method for thecompressible Navier–Stokes equations, Communications in Applied Mathematics and Computational Science 8 (1) (2013)99–122.[15] B. Muralidharan, S. Menon, A high-order adaptive Cartesian cut-cell method for simulation of compressible viscous flowover immersed bodies, Journal of Computational Physics 321 (2016) 342–368.[16] N. Gokhale, N. Nikiforakis, R. Klein, A dimensionally split Cartesian cut cell method for hyperbolic conservation laws,Journal of Computational Physics 364 (2018) 186–208.[17] R. Klein, K. Bates, N. Nikiforakis, Well-balanced compressible cut-cell simulation of atmospheric flow, PhilosophicalTransactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 367 (1907) (2009)4559–4575.[18] M. J. Berger, J. Oliger, Adaptive mesh refinement for hyperbolic partial differential equations, Journal of ComputationalPhysics 53 (3) (1984) 484–512.[19] G. Strang, On the construction and comparison of difference schemes, SIAM Journal on Numerical Analysis 5 (3) (1968)506–517.[20] R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Vol. 31, Cambridge University Press, 2002.[21] Z. Dragojlovic, F. Najmabadi, M. Day, An embedded boundary method for viscous, conducting compressible flow, Journalof Computational Physics 216 (1) (2006) 37–51.[22] E. F. Toro, Riemann solvers and numerical methods for fluid dynamics: a practical introduction, Springer Science &Business Media, 2013.[23] S. Mauch, A fast algorithm for computing the closest point and distance transform, Tech. Rep. caltechASCI/2000.077,California Institute of Technology (2000).[24] W. E. Lorensen, H. E. Cline, Marching cubes: A high resolution 3d surface construction algorithm, in: ACM siggraphcomputer graphics, Vol. 21, ACM, 1987, pp. 163–169.[25] C. G¨unther, D. Hartmann, L. Schneiders, M. Meinke, W. Schr¨oder, A cartesian cut-cell method for sharp moving bound-aries, AIAA Paper 3387 (2011) 27–30.[26] B. Paul, Calculating The Area And Centroid Of A Polygon), URL , Accessed: 16-02-2018.[27] Z. Wang, Improved formulation for geometric properties of arbitrary polyhedra, AIAA journal 37 (10) (1999) 1326–1327.[28] M. Meyer, A. Devesa, S. Hickel, X. Hu, N. Adams, A conservative immersed interface method for large-eddy simulationof incompressible flows, Journal of Computational Physics 229 (18) (2010) 6300–6317.[29] F. Capizzano, Turbulent Wall Model for Immersed Boundary Methods, AIAA journal 49 (11) (2011) 2367–2381.[30] Z. L. Chen, S. Hickel, A. Devesa, J. Berland, N. A. Adams, Wall modeling for implicit large-eddy simulation and immersed-interface methods, Theoretical and Computational Fluid Dynamics 28 (1) (2014) 1–21.
31] M. Kirkpatrick, S. Armfield, J. Kent, A representation of curved boundaries for the solution of the Navier–Stokes equationson a staggered three-dimensional Cartesian grid, Journal of Computational Physics 184 (1) (2003) 1–36.[32] P. Kundu, L. Cohen, Fluid mechanics, 5th Edition, Academic Press, 2014.[33] D. Tritton, Experiments on the flow past a circular cylinder at low Reynolds numbers, Journal of Fluid Mechanics 6 (4)(1959) 547–567.[34] Y.-H. Tseng, J. H. Ferziger, A ghost-cell immersed boundary method for flow in complex geometry, Journal of Computa-tional Physics 192 (2) (2003) 593–623.[35] A. S. Grove, F. Shair, E. Petersen, An experimental investigation of the steady separated flow past a circular cylinder,Journal of Fluid Mechanics 19 (1) (1964) 60–80.[36] C. Wieselsberger, New data on the laws of fluid resistance, Physikalische Zeitschrift 22 (1921).[37] M.-C. Lai, C. S. Peskin, An immersed boundary method with formal second-order accuracy and reduced numericalviscosity, Journal of Computational Physics 160 (2) (2000) 705–719.[38] H. Glaz, P. Colella, I. Glass, R. Deschambault, A numerical study of oblique shock-wave reflections with experimentalcomparisons, in: Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, Vol.398, The Royal Society, 1985, pp. 117–140.[39] M. Al-Marouf, R. Samtaney, A versatile embedded boundary adaptive mesh method for compressible flow in complexgeometry, Journal of Computational Physics 337 (2017) 339–378.[40] H. Uddin, R. Kramer, C. Pantano, A Cartesian-based embedded geometry technique with adaptive high-order finitedifferences for compressible flow around complex geometries, Journal of Computational Physics 262 (2014) 379–407.[41] A. P. Krasil ' shchikov, V. P. Podobin, Experimental study of sphere aerodynamic characteristics in free flight up to m ˜15, Fluid Dynamics 3 (4) (1972) 132–134.[42] A. Bailey, J. Hiatt, Sphere Drag Coefficients for a Broad Range of Mach and Reynolds Numbers, AIAA Journal 10 (11)(1972) 1436–1440.shchikov, V. P. Podobin, Experimental study of sphere aerodynamic characteristics in free flight up to m ˜15, Fluid Dynamics 3 (4) (1972) 132–134.[42] A. Bailey, J. Hiatt, Sphere Drag Coefficients for a Broad Range of Mach and Reynolds Numbers, AIAA Journal 10 (11)(1972) 1436–1440.