Developing an immersed boundary method for compressible flow

Development of a 3D sharp interface ghost node immersed boundary method (IBM) within a fluid solver based on the compressible Navier–Stokes equations is underway. The objective of the IBM is to accurately apply boundary conditions at fluid–solid interfaces. The Navier–Stokes solver is currently being verified using various test cases, including the classic cylinder in cross flow. As part of this verification process, particular attention was given to investigating the effects of different lateral boundary conditions. The results demonstrate that extrapolation boundary conditions exhibit better agreement with the literature compared to symmetry conditions, in cases with relatively narrow domains. These findings highlight the potential benefits of extrapolation boundary conditions in reducing confinement effects and removing nonphysical waves in external flow problems.

cylinder in cross flow, in Section 4, and present and discuss the verification results in Section 5.The work presented in this article builds on previous work [11], and parts of the text have been reused with adaptation.

NAVIER-STOKES EQUATIONS FOR COMPRESSIBLE FLOW
Compressible fluid flow is governed by the Navier-Stokes equations: Here,  represents time,  is the mass density, ⃗  = [, , ]  denotes the velocity,  =  +   represents the total specific enthalpy,  is the total specific energy,  represents pressure,  is the identity matrix, ∇ is the nabla operator,  denotes thermal conductivity,  represents temperature, and represents the components of the viscous stress tensor  for a Newtonian fluid, where ,  ∈ {, , } denote Cartesian coordinates.Here,  represents dynamic viscosity, and   denotes the Kronecker delta.The equations of state for a perfect gas yield:  = ( − 1) where is the ratio of specific heats at constant pressure and volume, c p and c v , respectively.The transport properties  and  are modeled as functions of temperature  only.The viscosity is given by Sutherland's law: where  0 and  0 are reference states for  and , respectively and  represents Sutherland's constant divided by  0 .Assuming constant Prandtl number , the thermal conductivity is given by To non-dimensionalize the governing equations, we introduce the following non-dimensional quantities, denoted by superscript * for non-dimensional values and subscript 0 for reference states: Here,  is a characteristic length scale, and denotes the Reynolds number of the reference state.By employing these non-dimensional formulations, the non-dimensionalized governing equations retain the same form as the dimensional ones, and hence, the superscript * will be dropped henceforth.The reference Mach number is defined as , where  0 represents the speed of sound in the reference state.To minimize cancellation errors for low Mach number flow, the perturbation formulation [12] is adopted.With this approach we solve for the flow variables' perturbations around a reference state rather than the variables themselves.We decompose the flow variables as: where  represents any flow variable,  0 denotes the reference value and  ′ is the perturbation.This perturbation formulation is applied to conserved and primitive variables but not to transport properties.By setting the reference state of velocity to zero, the governing equations retain their forms, except for the convective flux in the energy equations, which changes from  ⃗  to (() 0 + () ′ ) ⃗ .

Discretization in Space
The spatial discretization of the fluid domain is done by second order central difference operators.This approach was chosen for its compact numerical stencil.Thus, first derivatives, for example, in the continuity equation are approximated as Here, ,  are grid indices in the  and  directions respectively.Mixed derivatives, as found in the viscous terms of the Navier-Stokes equation, are discretized by applying the central difference operator twice, In the case of repeated derivatives, applying the operator twice would result in a wider stencil.Instead, the compact second derivative scheme is used, where the viscosity at the adjacent midpoints are computed as averages, The compact stencil Equation ( 10) is also used for the thermal diffusion terms in the energy equation and the heat equation.Besides being compact this discretization has the added benefit of better damping of high wave number oscillations.

Discretization in time
The spatial discretization of the Navier-Stokes equations for compressible flow Equation (1) in the interior of the domain leads to the system of ordinary differential equations The temporal discretization is done using the fourth order explicit classical Runge-Kutta method (RK4). where 2, = 4, = The time-step size Δ is chosen with regards to the stability criterion [13]: where and The constants  CFL = 0.9 and  4 = 0.68 are set smaller than their theoretical maximum values, 2 and 1.5, based on experience, to ensure stability of the solution.

Boundary conditions at the boundary of the domain
At the boundary of the domain we use combinations of the following different boundary conditions: inlet, outlet, periodic, symmetry, and extrapolation.At the inlet, we set a fixed inlet velocity, along with a temperature fluctuation of  = 0, and the pressure is linearly extrapolated from the interior.The conserved variables , , , ,  are derived from the primitive variables , , , , , using the ideal gas relation.The outlet boundary fixes only the pressure fluctuation to be  = 0, while the momentum components and density are linearly extrapolated from the interior.At the periodic boundaries, the values for all flow variables are copied from the boundary-adjacent nodes at the opposite side of the domain.The symmetry condition represents a full-slip wall, providing confinement without friction.Numerically, this is done using a layer of ghost nodes.The ghost nodes copy the values of the primitive variables from the boundary adjacent nodes in the interior, except for the normal velocity component, which changes sign.The conserved variables are derived from the primitive ones.The last boundary condition involves extrapolating all conserved variables from the interior.This is also implemented using a layer of ghost nodes.The flow variables in these nodes are linearly extrapolated from the interior nodes that coincide with the boundary and the interior boundary adjacent nodes:  ghost = 2 boundary −  adjacent , where  is any conserved variable.This boundary condition allows flow over the boundary, but should not be used if considerable inflow is expected, as the extrapolation would be based on downstream values.To approximate the flow variables at the ghost nodes, the following procedure is followed.A normal probe is extended from each ghost node to the closest point on the interface, referred to as the body intercept point.This normal probe is further extended into the fluid domain to determine the position of an image point.By performing trilinear interpolation of the surrounding fluid nodes, we can compute the flow variables at the image points.As depicted in the figure, values from the boundary conditions at the body intercept point can also be utilized in the interpolation process.The values at the ghost nodes are determined based on the boundary conditions at the body intercept points and the values at the image points.

Boundary conditions at immersed boundaries
To enforce no-slip conditions at the immersed surfaces, Dirichlet conditions are applied to the velocity components.The velocity ⃗  BI in the body intercept point will be 0 for a stationary surface, and the surface velocity otherwise.We employ linear extrapolation from the body intercept and image points to the ghost node: The immersed surfaces are modeled as adiabatic walls, resulting in a zero normal gradient for the temperature at these surfaces.Additionally, we assume the boundary layer approximation for the pressure, implying zero normal pressure gradient.These homogeneous Neumann conditions for  and  are implemented as: Here, ⃗  BI − ⃗  GN represents the vector from a ghost node to its body intercept point, as illustrated in Figure 1, and

TEST CASE
One of the test cases used to verify the method is the cross flow around a cylinder.The test case is set up as a rectangular domain with artificial boundaries, as seen in Figure 2. The flow wise direction aligns with the positive x-axis, while the axial direction of the cylinder corresponds to the z-axis, in which the flow does not change.The length of the domain in the flowwise direction is set to 40 times the diameter of the cylinder, and the width in the lateral direction along the y-axis is 20 times the cylinder diameter.The implementation of the boundary conditions at the domain boundary are outlined in Section 3.3.In this case the inlet is at  = 0, on the left side, and on the opposite side we have the outlet at  = 40.At the lateral boundaries  ∈ {0, 20} we simulate two different configurations: symmetry conditions and extrapolation.The cylinder is modeled as an adiabatic surface with no-slip conditions and the boundary layer approximation for the pressure, as discussed in Section 3.4.
The simulations presented here were conducted at Reynolds number  0 = 40 and Mach number  0 = 0.25, on a grid of 1600 × 802 × 3 nodes.The initial condition is uniform stagnation flow at the reference condition, and the inlet velocity linearly builds up to the target value 0.25 over the first 50 units of time. is the drag coefficient.The other variables are explained in Figure 3.

RESULTS
The results of the cylinder flow simulations with two different lateral boundary conditions are presented and discussed.Table 1 compares experimental and simulated values from the literature with those obtained in the present study.Overall, the results align well with the values reported in the literature.However, it is notable that the drag coefficient  D is significantly higher in the present study.We attribute this difference to the confinement effect of the boundaries.The width of the computational domain relative to the cylinder diameter has a substantial impact on the flow structures near the cylinder, particularly on the drag coefficient.For instance, when using symmetry conditions at the lateral boundaries, the calculated  D values are 2.58, 2.05, and 1.75 for domain widths of 5D, 10D, and 20D, respectively.The significance of the domain width has also been emphasized by Fornberg [17] and Linnick & Fasel [18].It should be noted that, in the other listed studies, the domain width is typically 40D or larger, whereas our verification process has been limited to 20D thus far.
Another observation is that the results obtained with extrapolation boundary conditions are closer to the results from other studies compared to those with symmetry conditions.Our hypothesis is that since the extrapolation allows flow over the boundary and thus imposes less confinement on the flow, it acts more as far-field boundary conditions, unlike the symmetry conditions.Additionally, the extrapolation boundary conditions allow density and pressure waves to pass through without reflection, rather than nonphysically bouncing back and forth between the boundaries.These two features are advantageous for simulating external flows.The reduced confinement may alleviate the demand for a larger domain width, while the non-reflective nature of the extrapolation boundary conditions may lead to faster convergence or reduced noise in transient simulations.
As for the temporal development, when the inlet velocity was increased from zero to the target value, this velocity propagated through the domain.Even though the velocity was increased gradually, this still created a planar density wave that traveled back and forth between the inlet and outlet boundaries.The wave was smooth and did not have a clear wavefront, but it still caused the flow to pulsate a bit.Consequently the separation bubble was elongated and compressed slightly, as the wave traveled with and against the flow direction.Naturally, this caused the output parameters in Table 1 to sway as well.At the end of the present simulations, the pulsation was still there, though much weaker than at the start.This means that the flow did not converge fully.The values in the table were obtained by computing means over 40 sample points in time close before the end of the simulation.Despite the pulsation, the flow stayed fully symmetric about the middle plane  = 10 throughout the simulation.There was no meandering of the wake, and no alternating vortex shedding.

CONCLUSIONS
The Navier-Stokes solver for compressible flow, utilizing a sharp interface ghost point immersed boundary method, is being verified using the cylinder cross-flow test case.The results show good agreement with experimental and simulated values reported in the literature, except for the drag coefficient.The higher drag coefficient observed in this study can be attributed to the narrower domain width compared to other investigations.It is worth noting that the extrapolation boundary conditions seem to reduce confinement effects and allow waves to exit the domain without reflection.Consequently, they may reduce the domain size requirement and convergence time for simulations of external flows.To draw definitive conclusions, further tests on external flow problems with different boundary condition types, e.g.periodic boundary conditions, are recommended.

A C K N O W L E D G M E N T S
The authors acknowledge the funding of the current research by the Research Council of Norway in the research project no.303218 entitled "Virtual Surgery in the Upper Airways-New Solutions to Obstructive Sleep Apnea Treatment (VirtuOSA)".

Figure 1 F I G U R E 1
Figure 1 illustrates the points close to an immersed boundary, along with their classifications.Ghost nodes are nodes located outside the fluid domain but included in the numerical stencil of at least one fluid node.

F I G U R E 2
Illustration of the cylinder test case.D is the cylinder diameter.TA B L E 1 Properties from the cylinder test case at Reynolds number 40.Experimental and simulated values from the literature and the present study.

F I G U R E 3
Schematic of the notation for parameters of interest. s is the separation length,  and  give the positions of the cores of the counter vortices relative to the cylinder's trailing edge, and  is the separation angle.The cylinder diameter is  = 1.The figure is taken from[10].