**2. Mesh Deformation Methods**

In this section, the RBF and IDW methods are introduced and compared. In addition, an improved RBF method to remedy its shortcomings is proposed. A mesh for a circular cylinder was deformed to the mesh of NACA 0012 using the RBF, IDW, and proposed methods. The red points in Figure 2 denote the control points on the surface of the circle and the blue ones indicate the displaced points on the surface of NACA 0012. The number of control points is 256 and the points move only in the vertical direction. The vertical movements make the normal vectors of the boundary faces rotate.

The initial grid for the circular cylinder is shown in Figure 3. The grid is a polyhedral mesh generated using snappyHexMesh, a standard built-in utility of OpenFOAM (The OpenFOAM Foundation Ltd., London, U.K.). The number of cells is 103,160 and the thickness of first boundary layer cell is 1% of the cylinder diameter. The top, bottom, right, and left boundaries were set as fixed boundary conditions. The front and back boundaries were set as symmetric boundary conditions. To satisfy the fixed boundary condition, the grid points of the fixed boundary were added to the control points. The displacements of the fixed boundary were set as zero.

**Figure 2.** Control points (red) and displaced points (blue).

#### *2.1. RBF Method*

The RBF method is an interpolation method that uses the distance between the grid point and control points. It was proposed by Sieger et al. [7], Boer et al. [11], Botch and Kobbelt [12], Jakobsson and Amoignon [13], and Michler [14]. The basic formula for displacement is presented in Equation (1).

$$\mathbf{D}(\mathbf{x}, \mathbf{y}, \mathbf{z}) = \mathbf{a}\_1 + \mathbf{a}\_2 \mathbf{x} + \mathbf{a}\_3 \mathbf{y} + \mathbf{a}\_4 \mathbf{z} + \sum \mathbf{w}\_i \mathbf{U}(\mathbf{R}\_i) \tag{1}$$

$$\mathbf{R}\_{\mathbf{i}} = \left| \left( \mathbf{x}\_{\mathbf{i}}, \mathbf{y}\_{\mathbf{i}}, \mathbf{z}\_{\mathbf{i}} \right) - \left( \mathbf{x}, \mathbf{y}, \mathbf{z} \right) \right| \tag{2}$$

Here, *Ri* is the distance between the grid point and the *i*-th control point, and *U* is a basis function. In this study, a thin plate spline (TPS) is applied as a basis function. The TPS provides a minimal and smooth displacement distribution. The details about basis functions can be found in [15].

$$\text{Thin plate spline}: \text{ : } \mathcal{U}(R) = |R|^2 \log R \tag{3}$$

The unknowns of Equation (1), *ai* and *wi*, are obtained by calculating Equation (4). The partitioned matrix *K* is determined by the distance between the control points, and *P* is composed of coordinates of the control points. Column vector <sup>→</sup> *v* contains the displacements of control points. Equation (4) is calculated with LU decomposition instead of an iterative method because the matrix is dense and small.

Because the size of the matrix in Equation (4) is proportional to the number of control points, every four points on the fixed boundary is added to the control points. Sheng and Allen [16] applied a greedy data reduction algorithm to reduce the matrix size and calculation time. Coulier and Darve [17] developed the inverse fast multipole method to reduce the computational time. After the mesh deformation, the normal component of displacement of the grid points on the symmetric plane is removed to satisfy the symmetric boundary condition.

$$
\begin{bmatrix}
\begin{array}{c}
K & P \\
P^T & O
\end{array}
\end{bmatrix}
\begin{bmatrix}
\begin{array}{c}
\stackrel{\rightarrow}{w} \\
\stackrel{\rightarrow}{a}
\end{array}
\end{bmatrix} = 
\begin{bmatrix}
\begin{array}{c}
\stackrel{\rightarrow}{v} \\
\stackrel{\rightarrow}{o}
\end{array}
\end{bmatrix}
\tag{4}
$$

$$K\_{ij} = \mathcal{U} \Big( \Big| \left( \mathbf{x}\_{i\prime} \, \mathbf{y}\_{j\prime} \, \mathbf{z}\_i \right) - \left( \mathbf{x}\_{j\prime} \, \mathbf{y}\_{j\prime} \, \mathbf{z}\_j \Big) \Big| \Big) \Big| \Big. \tag{5}$$

$$P = \begin{bmatrix} 1 & \mathbf{x}\_1 & \mathbf{y}\_1 & \mathbf{z}\_1 \\ 1 & \mathbf{x}\_2 & \mathbf{y}\_2 & \mathbf{z}\_2 \\ \vdots & \vdots & \vdots & \vdots \\ 1 & \mathbf{x}\_p & \mathbf{y}\_p & \mathbf{z}\_p \end{bmatrix} \tag{6}$$
 
$$\begin{bmatrix} & 0 & 0 & 0 & 0 \end{bmatrix}$$

$$O = \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} \tag{7}$$

$$
\stackrel{\rightarrow}{w} = \begin{bmatrix} w\_1 \\ w\_2 \\ \vdots \\ w\_p \end{bmatrix} \stackrel{\rightarrow}{a} = \begin{bmatrix} a\_1 \\ a\_2 \\ a\_3 \\ a\_4 \end{bmatrix} \stackrel{\rightarrow}{v} = \begin{bmatrix} \stackrel{\rightarrow}{D\_1} \\ \stackrel{\rightarrow}{D\_2} \\ \vdots \\ \stackrel{\rightarrow}{D\_p} \end{bmatrix} \stackrel{\rightarrow}{o} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} \tag{8}
$$

**Figure 3.** Initial grid shape for circular cylinder. (**a**) whole domain (side view), (**b**) surface mesh (perspective view), (**c**) around front part of the body (leading edge), (**d**) around rear part of the body (trailing edge).

The deformed mesh using the RBF method is illustrated in Figure 4. The overall domain is deformed smoothly, but the boundary layer is thinner than the initial grid and the non-orthogonality is worse than in the initial grid. The maximum skewness and maximum non-orthogonality are compared in Table 1. The mesh deformation takes approximately 45 s. The results are similar to those of He et al. [9], who conducted similar mesh deformation using the RBF and IDW methods with a 2D structured grid.

**Figure 4.** Deformed mesh using the radial basis function (RBF) method. (**a**) whole domain (side view), (**b**) near body (side view), (**c**) around leading edge (side view), (**d**) around trailing edge (side view).

**Table 1.** Quality of deformed mesh by the RBF method.


#### *2.2. IDW Method*

The IDW method is an interpolation method. The displacement of the grid points is calculated using Equation (9), where *Ri* denotes the distance as defined in Equation (3). *Mi* and *Ti* in Equation (10) are the rotation matrix and translation displacement of the boundary face, respectively. The weighting function is given by Equation (11). Luke et al. [8] suggested these values: *a* = 3, *b* = 5, and α = 0.25. *Lde f* is recommended as the maximum distance of mesh points, and *Ai* is the area of the boundary face. In this study, *a*, *b*, and α are set as 3, 0, and 0, respectively. Because of the irregular face area distribution of the boundaries, the thickness of the boundary layer becomes uneven. The calculation time is also reduced since the second term is not calculated.

$$D(x, y, z) = \frac{\sum w(R\_i)s(R\_i)}{\sum w(R\_i)}\tag{9}$$

$$s(R\_i) = M\_i R\_i + T\_i - R\_i \tag{10}$$

$$\mathbf{w}(R\_i) = A\_i \times \left[ \left( \frac{L\_{dcf}}{R\_i} \right)^a + \left( \frac{aL\_{dcf}}{R\_i} \right)^b \right] \tag{11}$$

Figure 5 depicts the deformed mesh using the IDW method. The non-orthogonality of the boundary cell is much better than that by the RBF method. The thickness of the boundary cell is almost equal to that of the initial grid except for the leading and trailing edges. The thickness around the leading and trailing edges is slightly larger than that of the initial grid. The displacements of grid points far from the deformed boundary are smaller than those by the RBF method. The grid quality of deformed mesh is compared with that of the initial grid in Table 2. The time for mesh deformation is approximately 115 s, which is ≈2.6 times larger than that by the RBF method.


**Table 2.** Quality of deformed mesh using the IDW method.
