Local Inverse Mapping Implicit Hole-Cutting Method for Structured Cartesian Overset Grid Assembly

An automatic hole-cutting method is proposed to search donor cells between a structured Cartesian mesh and an overlapping body-fitted mesh. The main flow is simulated on the structured Cartesian mesh and the viscous flow near the solid boundary is simulated on the body-fitted mesh. Through the spatial interpolation of flux, the surface boundary information on the body-fitted mesh is transferred to the Cartesian mesh nodes near the surface. Cartesian mesh box near a body-fitted mesh cell is selected as a local inverse map. The Cartesian nodes located inside the donor cells are marked by the relative coordinate transformation, so that all Cartesian nodes can be classified and the hole boundaries are implicitly cut. This hole-cutting process for overset grid assembly is called Local Inverse Mapping (LIM) method. In the LIM method, spatial interpolation of flux is carried out synchronously with the marking of Cartesian nodes. The LIM method is combined with the in-house finite-difference solver to simulate the unsteady flow field of moving bodies. The numerical results show that the LIM method can accurately mark the Cartesian hole boundary nodes, the search efficiency of donor cells is high, and the result of spatial interpolation is accurate. The calculation time of overset grid assembly (OGA) can be less than 3% of the total simulation time.


Introduction
Cartesian mesh is simple to generate, requires fewer computing resources, and requires no update for moving bodies. Peskin [1] proposed the immersed boundary method (IBM), which makes it feasible to calculate the flow field of objects with complex shapes in a Cartesian grid. The "Euler-Cartesian grid" is used to represent the fluid region and the "Lagrangian body surface grid" is used to represent the solid region. The effect of solids on fluids is realized through the construction and dispersion of "surface force sources". By using the "spatial discrete function", the Lagrangian force source term on the body surface grid nodes is extrapolated to a group of nearby Euler fluid Cartesian grid nodes, and the discrete Euler force source term is obtained. Then the Euclidean source term is put into the governing equations to solve the flow field to satisfy the boundary condition of the body surface.
Immersed boundary method has a simple force source structure and is easy to solve, but the mesh resolution near the wall surface is low, resulting in low solution accuracy of the flow field near the wall [2][3][4]. For viscous flow simulation with a high Reynolds number, a large number of grid nodes is needed to ensure adequate grid resolution near the wall. A hierarchical Cartesian grid is often used to improve the grid resolution and reduce grid node requirements. Chen Congcong et al. [5] states that a hierarchical Cartesian grid (as shown in Figure 1) can save 70% of computational grids compared to the simple stretch grid. The overset grid method [6], which solves the viscous flow on the body-fitted mesh around the wall and solves the main flow on the Cartesian mesh, can retain the advantages of simple generation, low memory requirements and high computational efficiency of Cartesian mesh. The overset grid method consists of two core technologies: Overset Grid Assembly (OGA) and Interpolated Data Transfer (IDT). The OGA process marks the Cartesian grid nodes that overlap the body-fitted mesh, and finds out the donor cells for information transfer. The IDT between two sets of grids is realized by spatial interpolation of physical flux.
The traversal method [7,8] is the most common OGA method used to identify the positional relationship between the body-fitted mesh and the Cartesian mesh. In the case of three dimensions of the traversal method, a very small Cartesian mesh spacing will result in a sharp increase in the computation consumption. Stencil Walk method or Stencil Jump method [9], and Neighbor-to-Neighbor Search method [10,11] are improved traversal methods, which can increase the computational efficiency of OGA to some extent but still require a huge amount of computation for 3D meshes. The Wall-Distance method [12,13] reduces the number of donor cells to be searched, thus reducing the amount of computation. When the number of Grid nodes exceeds 10 million, parallel computing method [14][15][16] should be adopted for Overset Grid Assembly to improve computing efficiency. The alternating digital tree (ADT) algorithm [17,18] applied to the adaptive Cartesian grid can effectively reduce the computation cost of OGA, because of its efficient data query structure. However, such methods cannot be applied to the structured Cartesian mesh.
In the structured overlapping Cartesian grid, WENO [19] scheme and other high-order finite difference schemes can be used to give full play to the advantages of the Cartesian grid. NASA carried out a series of studies in the 1980s and 1990s, and proposed hole mapping method [20], inverse map method [21] and X-Ray method [22], which significantly improved the OGA computing efficiency of the structured overset grid. At present, the research progress of the structured overset grid assembly method is not much, except Exact Inverse Map (EIM) method [23].
In order to improve the computational efficiency of structured overset grid assembly, we proposed an implicit hole cutting (IHC) method [24,25] named Local Inverse Mapping (LIM) method, for fluid points near the surface (FPNS, as shown in Figure 2) marking and flux interpolation between Cartesian mesh and overlapping body-fitted mesh. Structured Cartesian mesh is used as a reference scale, Newton iteration method [26] and local inverse map are used to calculate the relative position between a given point and a donor cell of body-fitted mesh. Thus, marking the Cartesian mesh nodes near the body surface. At the same time, spatial interpolation [4] is carried out on FPNS points.  Figure 2 shows a Cartesian mesh block consisting of nine points (yellow dotted line range) surrounding a hexahedral cell of the body-fitted cylinder mesh. This Cartesian grid block is defined as a "local inverse mapping box". The relative coordinates of those nine points of the local inverse mapping box are calculated, respectively. If the relative coordinate value of a point is within the range of [0, 1], this point is marked as Fluid Point Near Surface (FPNS). The flux at the FPNS point needs to be interpolated from the donor cell of the body-fitted mesh by Lagrange polynomial interpolation formula. Repeat the above steps to traverse all hexahedral donor cells of body-fitted mesh, thus all FPNS points on the Cartesian mesh are marked. Finally, through an X-Ray method, all cartesian mesh nodes located inside the solid surface are marked.

1.
On one time-step "t", a suitable "donor cell" of the Cartesian grid (the red quadrilateral in Figure 2), surrounding an extra boundary point of the body-fitted grid is found, and the flow field information on that "donor cell" of the Cartesian grid is interpolated into the boundary point of the body-fitted grid as the "additional boundary condition"; 2.
By solving the boundary layer flow field near the surface on the body-fitted grid, and then the information of the surface force at time-step "t + dt/2" can be obtained; 3.
The fluid points on the Cartesian gird overlapping with the body-fitted grid are defined as fluid points near the surface (FPNS, as shown in Figure 2). Searching the "donor cell" of the body-fitted grid (the yellow quadrilateral in Figure 2) surrounding an FPNS point. Then, the surface flow information at the "t + dt/2" time-step on the body-fitted grid is interpolated to the FPNS point through the spatial interpolation formula. Repeat the above step to obtain surface flow information at all FPNS points; 4.
Finally, the updated surface flow information on FPNS points is used as an "additional boundary condition" to solve the flow field on the whole Cartesian grid, and the solution of the time-step "t + dt" is completed.

Flux Interpolation Methods
Taking the donor cell hexahedron shown in Figure 3 as an example, the interpolated physical flux φ P at point P can be calculated according to the interpolation formula [4]: where, φ (0∼1,0∼1,0∼1) is the physical flux at a vertice of the donor cell, which is under the coordinate system (ξ, η, ζ). According to (ξ P , η P , ζ P ), the coordinates of point P, the interpolation coefficients C i,j,k (ξ, η, ζ) at eight vertices can be calculated as follows:

Strategies for Donor Cells Searching and Cartesian Nodes Classification
For a Cartesian mesh point P in the local inverse mapping box as shown in Figure 2, its spatial position relationship between itself and the donor cell can be determined by solving the nonlinear equations and obtaining the relative coordinates. If the relative coordinate value of point P, (ξ P , η P , ζ P ) is within the range of [0, 1], it indicates that this point is inside the donor cell.
The nonlinear equations are established to solve the relative coordinates (ξ P , η P , ζ P ) of point P, by given the Cartesian coordinates (x P , y P , z P ) of point P, and the Cartesian coordinates x (i,j,k) , y (i,j,k) , z (i,j,k) , i, j, k ∈ {0, 1} of 8 vertex points of the donor cell, and substituting these coordinates of nine points into Equation (1). The form of this nonlinear system is as follows:

X-ray Method and Cartesian Hole-Cutting Results
The FPNS Cartesian nodes marked by the LIM method are shown as the lower left vertex of the red square in Figure 4a. The Cartesian solid nodes inside the cylinder are marked by the X-ray method. Take the leftmost boundary of the Cartesian mesh as a starting point and move the unmarked point to the right along the ray. If the adjacent point to the left of the moving point is marked and the point to the right is unmarked, then this point is recorded as the starting point. Continue to move the unmarked point to the right along the ray until the adjacent point to the left of the moving point is unmarked and the point to the right is marked, then record the point as the end point. Mark all points between the start point and the end point as solid points. Both the start and end points are marked as ghost-cells. The last end point is then used as the new ray starting point, and the preceding steps are repeated until all the inner solid points are marked.
The solid points marked by the X-Ray method are shown in the lower left vertex of the blue square in Figure 4b, while the ghost cell points are shown in the lower left vertex of the green square in Figure 4. The donor cells searching method and flux interpolation method adopted by LIM are simple and efficient. The computation cost of spatial interpolation and information exchange could be less than 3% of the flow field simulation.

Information Entropy Analysis of OGA Computing Process
Information entropy can be used to characterize the complexity of OGA computation. The more complex the OGA process is, the more redundant computation is involved and the higher the information entropy is. Information theory provides a theoretical foundation to quantify the information content, or the uncertainty, of a random variable represented as a distribution. Formally, let X be a discrete random variable with alphabet χ and probability mass function p(x), x ∈ χ. The Shannon entropy of X is defined as [27] The key to applying the concept of entropy to OGA problems lies in how to specify a suitable random variable X and define the probability function p(x). In most cases, these probability functions can be defined heuristically to suit the needs of individual applications [28]. The probability p(x) of the OGA calculation process can be defined as the probability that any cell on a body-fitted grid is a donor cell of a background grid node. Assuming that the scale of the background grid is the same as that of the body-fitted grid, the total number of nodes in the background grid is n b , and the total number of cells in the body-fitted grid is n a , then the probability mass function p(x) corresponding to the traversal method is as follows p(x) = 1 n a n b Shannon entropy of the traversal method is 1 n a n b log 1 n a n b = log(n a n b ) (6) Under the same assumption, since a body-fitted grid cell only needs to be judged by the 27 Cartesian background grid nodes surrounding it, the probability mass function p(x) corresponding to the LIM method is as follows Shannon entropy of the LIM method is The information entropy difference between the traversal method and the LIM method is as follows When n b is much larger than 27, the information entropy of the traversal method will be much larger than that of the LIM method. Therefore, through information analysis, it can be seen that the LIM method has lower information entropy, less redundant calculation and higher computational efficiency than the traversal method. Especially for the 3D case with a large number of background grid nodes, the difference in OGA computational efficiency is more obvious.

Numerical Solver
The LIM method is added to the in-house code developed by the authors for computational fluid dynamics simulation. The governing equations, spatial and temporal discrete methods, boundary conditions and turbulence models used in the CFD program are briefly introduced in this chapter.

Governing Equations
Governing equations on the fixed orthogonal Cartesian grid and the movable bodyfitted grid is the compressible Navier-Stokes equations, which can be expressed in the Cartesian coordinate system (x, y, z, t) as follows [29]: where, the vector U is the conserved flux, the vector F, G, H are the inviscid flux terms (convection terms), and vector F v , G v , H v are the viscous flux terms (dissipative terms). The expressions of the above 7 vectors are as follows: where, ρ is the density, p is the pressure, T is the temperature, e is the total energy, τ is the viscosity stress and κ is the thermal conductivity. (u, v, w) T is the velocity vector under the coordinate system (x, y, z), (u r , v r , w r ) T is the relative velocity vector, and (u e , v e , w e ) T is the convected velocity vector of the movable body-fitted grid.

Spatial Discretization
Because of the strong hyperbolic property of convection terms in N-S equations, upwind schemes, such as the WENO scheme, which has the characteristics of high discontinuity resolution, high difference accuracy, slight numerical oscillation and small numerical dissipation, should be selected. Therefore, our in-house CFD code adopts flux-vector splitting schemes [6] combined with 5th-order WENO scheme [19,30] to carry out the spatial dispersion of convection terms. The 6th-order central difference scheme is used for spatial discretization of dissipative terms due to its dissipative property.

Temporal Discretization
The explicit three-step and third-order R-K method [31] is adopted. Since two sets of nonuniform meshes, background meshes and surface meshes, are, respectively, used to solve the mainstream and boundary layer flows, the two sets of meshes can adopt different global time steps. A small local time step is used on the surface body-fitted grid and a large time step is used on the Cartesian grid. Every step of the flow field solution is pushed forwards in time on the Cartesian grid, and several steps on the body-fitted grid, so as to improve the computational efficiency. This time-step strategy works in a similar way as Jameson's dual time-step approach [32].

Boundary Conditions and Turbulence Models
Whitfield and Janus's "non-reflection boundary condition" [33] derived from Euler equations based on the characteristic theory, is applied to all boundaries of the Cartesian grid. While "Non-slip wall boundary condition" is adopted on the solid surface. Spalart-Allmaras turbulence model [34] is used for viscous flow simulation.

Numerical Results
In order to verify the LIM method, an in-house CFD solver was used to simulate the flow around a 2-D cylinder and contra-rotating open rotors. The boundary conditions of this computational domain are:

1.
Nonreflection boundary condition are applied at the top and bottom; 2.
Nonslip wall boundary condition is implemented on the cylinder surface; 3.
Velocity is a uniform flow at the left boundary and zero-gradient at the right boundary.
The time-evolution of the flow field can be quantitatively represented by the variation in the drag and lift coefficients at the surface of the body, which are defined as: where D is the diameter of cylinder, ρ ∞ is the free-stream density, and U ∞ is the free-stream velocity, F d is the drag force and F L is the lift force acting on the circular cylinder. The dynamic force F acting on the surface of the cylinder can be computed by integrating the local pressure and stress distributions along the cylinder surface, which can be written as: where p b is the pressure, I is the identity matrix, n = (n x , n y ) is the surface unit normal vector, and the subscripts x and y denote the x-direction and y-direction, respectively. τ is the stress tensor and S b is the surface of the solid body.

Steady Flow Past a Fixed Cylinder at Re = 40
This case is steady laminar flow past a circular cylinder at Re = 40 and the freestream Mach number Ma = 0.1. Figure 6 is the result of velocity magnitude and streamlines of flow past the fixed cylinder. L w /D is the non-dimensional wake bubble length. Table 1 shows the good agreement in the drag coefficientC d and L w /D.  To investigate the sensitivity of the results to grid resolution, unsteady flow past a stationary circular cylinder at Re = 185 and Ma = 0.1 is simulated with three different grids spacings, which is ∆x = ∆y = 0.01D, 0.02D, 0.04D, 0.08D, respectively. Coefficients ofC d and St with different ∆x min are shown in Table 2. According to the comparison of the results in Table 2, the grid spacing ∆x of the subsequent simulations in this paper is equal to 0.02D. The Strouhal number St is defined as follows:

Velocity magnitude and streamlines
where f 0 is the vortex natural shedding frequency.  Figure 7 presents the drag and lift coefficients for this case. The average drag coefficient, C d , the root-mean-square value of the lift coefficient, C Lr.m.s. , and the Strouhal number St are reported in Table 3. As shown in Table 3, St is in good agreement with Williamson's universal law [38], and the flow parameters agree with other experimental and numerical results very well. From this, the value of f o can be calculated and is used in the following simulation for the transversely oscillating cylinder.

Oscillating 2-D Cylinder
The problem of a cylinder transversely oscillating in the free stream is a classic flow problem with moving boundaries, which can be used to verify the accuracy of the LIM method applied to OGA and flow field interpolation.
There are various experimental and numerical investigations available (Khalili et al. [2], Guilmineau and Queutey [39], Lu and Dalton [40], Wu [42], Gu et al. [43], Schneiders et al. [44]). The above studies showed that the vortex shedding frequency can synchronize with the frequency of the oscillatory forcing. This so-called lock-in of the vortex shedding is usually observed when the forcing frequency approaches the natural shedding frequency of the flow past a fixed cylinder.
The Reynolds number based on the cylinder diameter D and the freestream velocity u is Re = 185. The vertical cylinder position is imposed as y(t) = A e sin(2π f e t) with A e and f e being the oscillation amplitude and the excitation frequency, respectively. The quantity f 0 denotes the natural vortex-shedding frequency obtained for the corresponding flow past a fixed cylinder. f e is very important for the flow structure. LIM method is used to perform calculations for Re = 185, A e = 0.2D. The range of f e / f 0 = 0.8, 0.9, 1.0, 1.1, 1.12, 1.2 is considered. All the simulations start from the flow past a non-oscillating cylinder until the vortex shedding reaches to the stable state. Figure 8 shows the drag and lift coefficients for different values of f e / f 0 obtained by the LIM method and WENO CFD solver. The primary changes in the drag and lift coefficient patterns can be seen as f e / f 0 increases. The impact of the higher harmonic oscillations can be easily seen in this figure for excitation frequencies higher than f e / f 0 ≥ 1.

Drag and Lift Coefficients Simulation Results
When f e / f 0 ≤ 1.0, the force behaves periodically in time because stable vortex shedding is established. However, periodicity is lost for f e / f 0 ≥ 1.1. Increasing f e ( f e / f 0 ≥ 0.9), the amplitudes for the lift coefficients are increasing. Since an opposite sign vortex is formed at the lower side of the cylinder as f e / f 0 increases ( f e / f 0 > 1), as can be seen in Figure 9, this phenomenon causes changes in the vorticity at the top and bottom of the shear layers. The above results show good agreement with reference [39,43].   Figure 9 presents the instantaneous streamlines with velocity in the left column and instantaneous spanwise vorticity contours in the right column. Three cases are shown where f e / f 0 is gradually increased from 0.8 to 1.12. It is observed that the wake pattern is suddenly changed from f e / f 0 = 1.0 to f e / f 0 = 1.12. Figures 10 and 11 are the comparison of the surface pressure coefficient for the cylinder located at the top-dead-center with the results of Guilmineau and Queutey [39]. Comparison results show good agreement of the distribution of the coefficient of pressure. (deg) Figure 11. Comparison of the surface pressure coefficient for the cylinder located at the top-dead center against reference [39], f e = 1.1.
The above simulation results show that the LIM method can accurately implement the assembly of the moving overlapping mesh, and ensure the correct flow field interpolation and simulation results.

Contra-Rotating Open Rotors
In order to further verify the accuracy and reliability of LIM method, the overset grid assembly calculation and unsteady flow simulation of the Contra-rotating open rotors are carried out. Figure 12 shows the surface mesh and body-fitted mesh of the Contra-rotating open rotor. Each blade block is an independent structured body-fitted grid with a node number of 0.204 million. The front and rear rotors have a total of 14 blades and 2.8 million grid nodes. Table 6 shows the number of nodes in the two sets of grids. The percentage of CPU time-consuming ratio of OGA is shown in Table 7. OGA computing is implemented on a personal computer with AMD Ryzen 3100 CPU. The traversal method and LIM method were used, respectively, to mark the FPNS Cartesian mesh points. The comparison of OGA calculation time is shown in Table 8. The comparison results show that the LIM method can reduce the computing time by 99% compared with the traversal method. Compared with the 2-D case, the computational efficiency of the LIM method is significantly improved in the 3-D case. Cartesian nodes overlapped with blade mesh were marked green by the LIM method as shown in Figure 13.

OGA Calculation Efficiency Comparison
This example demonstrates the applicability of the LIM method to the multibody case. The front and rear rotors rotate in opposite directions. Each blade grid block can move independently, and there are 14 moving blocks in total. There is no exchange of flow information between the blade grid blocks, and the flow information can be transmitted only through flux interpolation with the Cartesian mesh.

Flux Interpolation Strategy
High-Order Interpolation Method [47] and Conservative Interpolation Method [48][49][50] are very helpful to ensure the accuracy of simulation results of the overlapping grid. For structured grids, non-conservation interpolation [51] can also obtain satisfactory simulation results. The Lagrange polynomial spatial interpolation method used in this paper belongs to non-conservation interpolation. The physical quantity transferred by spatial interpolation is the vector U in Equation (11), rather than the original physical variables, such as density ρ and velocity u, v, w. Figure 14 shows the vorticity distribution unsteady simulation results of Contra-rotating open rotors. Figure 15 shows the simulation results of shock waves in an open rotor blade passage. The continuous vorticity and shock wave distribution verify that the LIM method is accurate and the spatial interpolation results are correct.

Parallel Computing Strategy of LIM
LIM method can support both Coarse-grained and Fine-grained parallel computing strategies. In the Coarse-grained parallelism strategy, each moving grid block, such as the blade grid block, occupies one thread and performs OGA calculation. For 14 blade multi-body grid blocks, 14 threads can be used for parallel computation. In the Fine-grained parallelism strategy, the structured body-fitted blade mesh can be divided into np subblocks (np = number of threads) according to the number of threads. LIM method is used to mark the FPNS nodes on each sub-block in parallel. After the OGA parallel calculation on all sub-blocks is completed, the nodes in the Cartesian mesh that overlap with the body-fitted mesh are marked, and the flux is interpolated from the blade mesh donor cells. OGA parallel computing code is under development, and the efficiency improvements will be compared in subsequent studies.

Discussion and Conclusions
Compared with the traversal OGA method, the local inverse mapping implicit holecutting method proposed in this paper can automatically perform the overset grid assembly calculation without manual intervention. Verification results show that this method is efficient and accurate. The calculation cost of points marking and flux interpolation accounts for 3% of the total calculation cost for the 2-D case and 14% for the 3-D case. This method could effectively reduce the computing consumption of donor cells searching and save memory.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: