The governing equations used in the finite volume solver for laminar flows on structured grids are presented here. The mass, momentum, and energy conservation laws for the compressible fluid flow in conservative differential form are ∂ q ⃗ c ∂ t + ∇ ⋅ F ⃗ I − F ⃗ V = 0 . (1) In the aforementioned equation, F ⃗ I constitutes the inviscid fluxes, and F ⃗ V constitutes viscous fluxes. The vector q ⃗ c consists of the conservative variables which includes mass, momentum, and energy per unit volume. q ⃗ c = ρ ρ u ρ v ρ w ρ E t , (2) F ⃗ I = ρ u ρ u 2 + p ρ u v ρ u w ρ u H t n x + ρ v ρ u v ρ v 2 + p ρ v w ρ v H t n y + ρ w ρ u w ρ v w ρ w 2 + p ρ w H t n z , (3) F ⃗ V = 0 t x x t x y t x z u t x x + v t x y + w t x z − q x n x + 0 t x y t y y t y z u t y x + v t y y + w t y z − q y n y + 0 t z x − τ z x t z y − τ z y t z z − τ z z u t z x + v t z y + w t z z − q z n z . (4) In the aforementioned equations, u, v, and w are the Cartesian components of velocity along the X, Y, and Z directions, respectively, ρ is the fluid density, p is the fluid pressure, and H _t is the total specific enthalpy. Also, t ̃ ̃ is the laminar stress tensor, and q ⃗ is the laminar heat flux. The quantities n _x , n _y , and n _z are the projections of the cell face area (vector) along the X, Y, and Z axes, respectively. In order to close the equations, the equation of state for an ideal gas is used as follows: P = ρ R T , (5) where R is the specific gas constant for air.

FEST-3D (Sandhu et al., 2020), an in-house developed parallel code for structured grids based on a finite volume framework, is used. The capabilities of this solver range from laminar to turbulent regimes comprising various spatial and time-discretization schemes. In the present work, the solver is run with the advection upstream splitting method (AUSM) (Liou and Steffen, 1993), which is an upwind scheme, for inviscid flux calculation, MUSCL for second-order reconstruction (spatial) as presented in Hirsch (2007), and preconditioned LUSGS (Kitamura et al., 2011)/RK4 for time integration. The use of upwind schemes for inviscid flux formulation with MUSCL reconstruction can be seen in the works of Hartmann et al. (2009), Tamaki and Imamura (2018), Qu et al. (2018), Anand Bharadwaj and Ghosh (2020), and Luo et al. (2006). All the simulations presented in this work are 2D computations. Molecular viscosity is modeled with Sutherland’s law, and a laminar Prandtl number of 0.72 is used for air. The scope of the work presented in this study is specific to laminar and two-dimensional flows.

As mentioned in Ramakrishnan et al. (2016), the immersed boundary is represented by a set of line segments, wherein an outward normal is also specified along with the coordinates of its endpoints. The cells in the entire domain are classified as field cells (white), band cells (blue), and interior cells (red) based on the distance of their centers with respect to the immersed surface determined using a signed distance approach, as shown for a NACA0012 airfoil in Figure 1. In addition, the cell faces shared by a band and field cell are designated as band faces.

The solution at the band faces is then reconstructed (solution forcing) to implicitly obey the desired boundary conditions at the immersed surface. Specifically, the boundary conditions considered for velocity forcing are no-slip and no-penetration conditions for viscous flows and only no-penetration conditions for inviscid flows. The variable reconstruction at the band face center is performed by first constructing an interpolation point normal to the immersed surface, as discussed as follows. This is different from the approach adopted in the work of Ramakrishnan et al. (2016), wherein only the solution at the nearest field cell to the band face was used for the solution reconstruction at the band face. The modification in the solution forcing approach was carried out to avoid the dependence of the forcing method on the relative location of the band face center with respect to the immersed surface, as required in the previous method. The velocity reconstruction strategy adopted in Ramakrishnan et al. (2016) is described as follows for the benefit of the reader.

This section describes the velocity reconstruction methodology proposed by Ramakrishnan et al. (2016), which is henceforward referred to as IBM(A). This velocity forcing at the band face is loosely based on the boundary conditions applied at ghost cells adjacent to walls in finite volume solvers on body-fitted grids. The reconstruction for the slip- and no-slip wall is outlined as follows.

In this work, we construct an interpolation point on the line passing through the band-face center and normal to the nearest immersed surface (line segment). The distance of the interpolation point from the solution forcing point (which can be a face center or cell center) is often fixed for the entire domain as a function of local grid spacing(s) (Chi et al., 2017) or a fixed value as in the study by Capizzano (2011). In this work, however, to determine the position of the interpolation point, an interpolation stencil of nine cells is constructed, which consists of the field cell that has the band face and its eight neighbors. This ensures that more number of field cells are used in interpolating the primitive variables at the interpolation point. An inverse distance approach proposed by Ghosh et al. (2010) is then used to first determine the location of the interpolation point.

To do this, the perpendicular distances of the field cell centers from the line passing through the band face center and normal to the nearest IB line segment (for instance, d p i ), as shown in Figure 4, are determined. Using the inverse of these distances as weights and the distance along the line (for instance, d n i ) as the variable, the distance of the interpolation point from the band face center along the line is determined. Finally, the primitive variables are also interpolated at the interpolation point based on the same algorithm by using the nearby fluid cell data.

The formulation is as shown as follows, wherein k is the total number of field cells in the 3 × 3 stencil considered. Here, d _n,IP is the distance of the interpolation point from the band face center along the normal line; ϕ _IP and ϕ _i are scalars (velocity components, pressure, and density) at the interpolation point and the field cell centers in the stencil, respectively. d n , I P = ∑ i = 1 k w i d n i ∑ i = 1 k w i , (11) ϕ I P = ∑ i = 1 k w i ϕ i ∑ i = 1 k w i , (12) w i = 1 d p i , i = 1 : k .

The basic idea here is to reconstruct the state at the band face such that the presence of the immersed boundary is “felt” by the rest of the flow. This is carried out by using the state (or some gradient information) of the flow variable at the immersed boundary and the data reconstructed at the interpolation point.

This test case simulates the inviscid transonic flow past a bump in a channel. The channel is 3.0 m long and 1.0 m in height. The bump is located halfway along the length of the lower wall. The thickness-to-chord ratio of the bump is 10%.

This internal flow test case was performed on a structured body-fitted grid by Favini et al. (1996), and the solution is compared with the same work. This simulation is carried out on grids with dimensions 130 × 42 and 180 × 64, for comparison of the present method with that presented in Ramakrishnan et al. (2016), and a sequence of successively refined grids 48 × 16, 96 × 32, 192 × 64, and 384 × 128 in the x and y directions, respectively, with clustering near the ends of the bump and minimum grid spacing—Δx _min = Δy _min = 0.005 m for the 192 × 64 grid—same as in Kumar et al. (2020). The inlet Mach number, pressure, and temperature are equal to 0.675, 1.0e5 Pa, and 300 K, respectively.

Boundary conditions applied are shown in Figure 7. The non-uniform mesh is shown in Figure 8 with the IB rendered as a yellow curve. The pressure contours obtained using the IBM from Ramakrishnan et al. (2016), referred to henceforward as IBM(A), and the present IBM are compared in Figure 9.

The contours indicate that both the IBM(A) and the proposed IBM produce similar results with the 130 × 42 grid. However, as the number of grid points is increased to 180 × 64, the solution using IBM(A) becomes non-physical whereas the present IBM generates physically correct results. It is to be noted that this incorrect behavior of the IBM(A) is present in spite of using the limiting function Eq. 7. Without the limiting function, the solution diverges.

The solution convergence with grid refinement can be observed from the surface plots of the Mach number in Figure 10A with the 192 × 64 grid and 384 × 128 grids having virtually identical surface Mach number distribution. Also, as shown in Figure 10B, the Mach number plots along the lower and upper surfaces of the domain (192 × 64) show overall excellent agreement with the results obtained by Favini et al. (1996).

Figure 11 shows the pressure contour plots for the present computations and the CUS-IBM (Kumar et al., 2020) solution on the 192 × 64 grid. The shock structure and location show good agreement between the solvers, with the present solver producing a somewhat smeared shock foot.

The exit Mach number is probed at the lower and upper surfaces (listed in Table 2) and is compared against the body-fitted grid values obtained by Favini et al. (1996). The percentage errors for the lower and upper surfaces are 2.3 % and 2.54%, respectively.

This is an external flow computation of flow over a cylinder. A supersonic flow of M _inlet = 3, P _inlet = 103320 Pa, and T = 300 K over a half-cylinder is simulated. The computational domain is [ −1 m, 0 m] × [ −2 m, 2 m]. The cylinder is centered at (0 m, 0 m) with a radius of 0.5 m.

Uniform Cartesian grids are chosen with grid sizes 50 × 200, 100 × 400, and 200 × 800 along the x and y directions, respectively, which are the same as reported in Kumar et al. (2020). A symmetry or slip-wall boundary condition is used for the top and bottom boundaries as the flow is expected to be parallel near these boundaries, with the bow shock exiting the right boundary. Boundary conditions used are as shown in Figure 12.

Figure 13A plots the pressure coefficient along the surface of the upper-right quarter of the cylinder on different grids. It can be observed that the pressure distribution on the fine and medium grids is close to each other compared to that on the coarse grid, which suggests that the solution on the fine grid can be considered to be grid-converged. A comparison is also made between the grid-converged solution on the fine grid (200 × 800) using the present IBM and the results of the CUS-IBM (Kumar et al., 2020) on the same grid. Good agreement between the predictions can be seen, as shown in Figure 13B.

Furthermore, post-shock stagnation pressure in the computational domain is also probed at y = 0 and compared to the calculated analytical value (listed in Table 2). The percentage error of the computed value is determined to be about 0.57%. The contour plot of the Mach number is shown in Figure 14, and a comparison is also made with the CUS-IBM (Kumar et al., 2020). The contour plots are in good agreement. It can be seen from the contours that there is a detached bow shock formed in front of the cylinder, and the solutions on the two different IBM solvers compare well.

This is a laminar flow simulation with a large separation vortex near the trailing edge of the airfoil. Free stream flow parameters are M _∞ = 0.8, p _∞ = 103320 Pa, T _∞ = 300 K, and Re _∞ = 500, and the angle of attack is 10°.

This test case is chosen as the flow is transonic and involves flow separation with a recirculation bubble forming on the suction side of the airfoil. This tests whether the IBM can accurately predict flow separation, which is very important for the reliable prediction of viscous flows.

The computations for this case are carried out on the non-uniform grids of sizes 162 × 256, 324 × 512, and 648 × 1,024 along x and y directions, respectively. The size of the grid 648 × 1,024 is based on the grid used in Kumar et al. (2020), with a minimum grid spacing of Δx _min = Δy _min = 0.005 m at the leading edge. Boundary conditions considered and the zoomed-in view of the grid used are shown in Figure 15. The solution of the present IB method is compared with the body-fitted grid simulation by Swanson and Langer (2016).

Figure 16A shows the grid convergence results obtained using the presented IB method. It can be observed that the surface pressure coefficient predictions on both the medium and fine grids overlap, and the solution on the 648 × 1,024 grid can be considered as grid-converged. In Figure 16B, pressure coefficient distribution on the airfoil with the present IBM (648 × 1,024) is further compared with the 1,280 × 512 C grid computations by Swanson and Langer (2016). The C _p distributions for the two methods match closely, indicating that the computations with the present IBM produce accurate surface pressure with a grid having a similar size as used for a body-fitted simulation.

As can be seen from Figure 17, which presents the Mach number contour plots using the present IBM and body-fitted grid simulations (Swanson and Langer, 2016), there is good agreement between the predicted solutions. In Figure 18, streamlines are shown in a zoomed-in view of the region near the trailing edge to compare the predictions of the recirculation zone between the present IBM and body-fitted simulation by Swanson and Langer (2016). It can be observed that the location and size of the recirculation regions look very similar. Furthermore, the errors in the lift coefficient c _l and the drag coefficient c _d , between the present computation and value reported by Swanson and Langer (2016) using body-fitted grid simulations, are determined to be about 1.07% and 2.6%, respectively. The exact values of the lift and drag coefficients are listed in Table 2.

Flow conditions for this external flow simulation are M _∞ = 2.0, p _∞ = 103320 Pa, T _∞ = 300 K, and Re = 300. The domain extent is 60 D × 40 D (D = 1 m), the same as in Kumar et al. (2020) and Qiu et al. (2016), and the cylinder is centered at (24,20).

Non-uniform meshes with grid sizes 158 × 108, 316 × 216, and 632 × 432 along x and y directions, respectively, were used for the grid convergence study. There is a uniform mesh region around the cylinder of size 1.7 D × 1.7 D with a grid spacing of 0.1 D, 0.05 D, and 0.025 D for three levels of grids as in Takahashi et al. (2014). This grid spacing of the finest grid used is considered sufficient to capture the boundary layer effects for the specific Reynolds number. The domain with the boundary conditions employed and the grid are shown in Figure 19.

As shown in Figure 20A, it can be observed that the predicted pressures on the medium and fine grids are virtually overlapping, and thus, the solution on the 632 × 432 grid can be considered as grid-converged.

The pressure coefficient on the cylinder surface plotted in Figure 20B compares the converged solution on the finer grid with the published results by Takahashi et al. (2014) and shows excellent agreement.

Furthermore, the error in the drag coefficient, c _d , between the present computation and the value reported by Takahashi et al. (2014) for Cartesian grid simulations is determined to be about 0.4%. The exact values of the lift and drag coefficients are listed in Table 2. Contour plots of the density ratio in Figure 21 show that the bow shock is captured well with the present IBM and the flow field is qualitatively similar to that predicted by Takahashi et al. (2014).

This section discusses the performance of the proposed IB solver by performing grid convergence tests on three grids of different sizes for two of the test cases presented earlier in this work: inviscid Mach 3 flow past a circular cylinder and viscous Mach 2.0 flow past a circular cylinder at Re of 300. For the inviscid flow test case, the three grids considered are 50 × 200, 100 × 400, and 200 × 800 on the computational domain of [−1,0]×[−2,2] with uniform grid spacing. For the viscous flow simulations, the three grid sizes chosen are 158 × 108, 316 × 216, and 632 × 432 with a uniform grid spacing of 0.1 D, 0.05 D, and 0.025 D, respectively, in the rectangular region 1.7 × 1.7 around the cylinder.

For the inviscid test case, post-shock stagnation pressure along the line y = 0 is considered to determine the error in the computation. The reference value in this case is the theoretical post normal-shock value at Mach 3.0 for air. For the viscous test case, the error in the solution is calculated by computing the drag coefficient C _d . An error function is computed as ϵ h = | ϕ − ϕ ref | ϕ r e f , where ϕ = [P ₀, C _D ] and ϕ _ref is the zero-grid-spacing reference value, obtained using Richardson’s extrapolation (Slater, 2022, Roache, 1994). The plot shown in Figure 22 shows the variation of the error function (ϵ _h ) epsilon versus the grid spacing parameter (r). The grid spacing parameter r is defined as the ratio of the average grid spacing of a grid with respect to the average grid spacing for the finest grid. The observed order of convergence obtained is in the range of 1.4–1.6, which is determined using that of Slater (2022). This suggests that the solution obtained from the proposed IBM is less than second-order accurate, although the finite volume method used here is second-order accurate. This can be attributed to the reduced order of solution reconstruction at cell faces in the neighborhood of the immersed surface, as discussed in Section 3.4.5. Furthermore, although linear interpolation is used for forcing the velocity at the band face, density and pressure are extrapolated, which can also contribute to the overall reduction in the accuracy of the solver.