Efficient mesh motion using radial basis functions with volume grid points reduction algorithm

Journal of Computational Physics 348 (2017) 401–415

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

Efficient mesh motion using radial basis functions with volume grid points reduction algorithm Liang Xie ∗ , Hong Liu School of Aeronautics and Astronautics, Shanghai Jiao Tong University, Shanghai, 200240, China

a r t i c l e

i n f o

Article history: Received 15 August 2016 Received in revised form 27 May 2017 Accepted 24 July 2017 Available online 27 July 2017 Keywords: Radial basis function Grid deforming Data reduction

a b s t r a c t As one of the most robust mesh deformation technique available, the radial basis function (RBF) mesh deformation has been accepted widely. However, for volume mesh deformation driven by surface motion, the RBF system may become impractical for large meshes due to the large number of both surface (control) points and volume points. Surface points selection procedure based on the greedy algorithm results in an efficient implementation of the RBF-based mesh deformation procedure. The greedy algorithm could reduce the number of surface points involved in the RBF interpolation while acquire an acceptable accuracy as shown in literature. To improve the efficiency of the RBF method furthermore, an issue that how to reduce the number of the volume points needed to be moved is addressed. In this paper, we propose an algorithm for volume points reduction based on a wall distance based restricting function which is added to the formulation of the RBF based interpolation. This restricting function is firstly introduced by the current article. To support large deformation, a multi-level subspace interpolation is essentially needed, although this technique was used to improve the efficiency of the surface points selection procedure in the existed literature. The key point of this technique is setting the error of previous interpolation step as the object of current step, and restricting interpolation region gradually. Because the tolerance of the error is decreased hierarchically, the number of the surface points is increased but the number of the volume points needed to be moved is reduced gradually. Therefore, the CPU cost of updating the mesh motion could be reduced eventually since it scales with the product of these two numbers. The computational requirement of the proposed procedure is reduced evidently compared with the standard procedure as proved by some examples. © 2017 Elsevier Inc. All rights reserved.

1. Introduction

Grid deformation technique is very important for the simulation of flows with a deforming boundary, such as fluidstructure interaction (FSI) problems [1–3], aerodynamic shape optimization [4] and numerical icing simulation [5]. Once the boundary is deformed, the interior points should be adjusted to avoid negative volume. Several algorithms have been developed to accomplish this task, such as the spring analogy [6,7], the linear elasticity analogy [8], Delaunay graph mapping [9], the interpolation method based on radial basis functions (RBF) [10–20], etc.

*

Corresponding author. E-mail addresses: [email protected] (L. Xie), [email protected] (H. Liu).

http://dx.doi.org/10.1016/j.jcp.2017.07.042 0021-9991/© 2017 Elsevier Inc. All rights reserved.

402

L. Xie, H. Liu / Journal of Computational Physics 348 (2017) 401–415

Among these methods, grid deforming method based on the RBF has been accepted widely in the past decade due to its simplicity, robustness, flexibility and achieved mesh quality. RBF was firstly used as a technique of grid deformation by Beckert and Wendland [10]. Their work has been extended to large-scale spatial coupling problems by Ahrem and his co-workers [3]. This technique has recently used by Allen and Rendall [1,2] to develop a unified meshless approach for fluid-structure coupling and mesh motion. The RBF-based mesh deformation technique has been employed in the field of aerodynamic shape optimization by Jacobsson and Amoignon [12]. To enhance the efficiency of the RBF method, Rendall and Allen [13,14] proposed a surface points selection procedure based on a greedy algorithm. To improve the efficiency of the greedy algorithm, an incremental approach has been proposed recently by Selim and his co-workers [16]. Their method decreases the computational complexity of the greedy algorithm from O ( N s4 ) to O ( N s3 ), where N s is the number of selected surface points. Although generally the RBF mesh deformation technique is good in preserving orthogonality, the displacements in different directions are uncoupled. Gillebaart and his co-authors [19] proposed a solution of incorporating orthogonality explicitly in the formulation by adding additional points near the boundary. Although this approach is particularly attractive for high quality mesh deformation, it is quite expensive for large mesh size. The computational cost of the RBF method for the mesh motion is of size N s × N v , where N s is the number of the surface points defining the motion, and N v is the number of the volume mesh points to be moved. Although the efficiency of the RBF technique used in the mesh deforming algorithm has been improved greatly by the greedy algorithm proposed by Allen and Rendall [13,14], the greedy algorithm itself requires very more CPU cost to select surface points. To alleviate the CPU requirement of the selection procedure itself, a multi-level subspace points selection algorithm has been developed by Wang and his collaborators [17]. It should be pointed out that the multi-level idea was employed by Wendland [21] to obtain a convergence results of RBF system before the research of Wang [17]. The works of Allen, Rendall [13,14] and Wang [17] only focused on the reduction of the number of surface points, namely N s . To reduce time requirement of mesh motion step, a new hybrid way has been developed by Wang and his co-workers [20], which combines the advantages of Delaunay graph and RBF method. Their algorithm provides perfect performance on the save of CPU cost of mesh updating procedure. Although the multi-level route [17] offers improved surface accuracy and reduces the CPU cost of selection procedure, it involves much more surface points. Therefore, the efficiency of volume grid deformation process is harmed by this algorithm. To solve this problem, a research [22] has been carried out by the first author of the current article. That article developed a new procedure to limit the deforming region based on the multi-level subspace points selection algorithm provided by Wang and co-workers [17], although they only use their technique to improve the efficiency of the surface points select procedure. The algorithm hierarchically restricts the volume region need to be moved at each level so that the number of volume points influenced by each level is reduced in a stepped manner. However, the previous work [22] is inconvenient and the improvement is not obvious because the restricting methodology is too complicated. This drawback motivates the present work. The principle of the algorithm of the previous and the current work [22] is that the deformation region required by the previous step is relatively large, but the number of the surface nodes chosen is small; the deformation region of next step can be set to be relatively small because the interpolation target is the error of the previous steps, although much more surface nodes are involved. The result of this algorithm is that when N s is large, N v is relatively small; when N v is large, N s is relatively small. Therefore, the computational cost of the grid motion could be saved significantly. Except the multi-level interpolation described above, another important key to construct efficient volume data reduction algorithm is how to limit the region where volume nodes need to be moved. This issue is not solved well in the previous paper [22]. Compared with the former research [22], a wall distance based restricting function has been introduced by the current article to fix up this problem by a simple way. The remainder of this paper is organized as follows: the RBF method used in the mesh motion algorithm is summarized briefly in Section 2; the greedy algorithm and multi-level subspace technique used in the surface points selection procedure are described in Section 3; a new method to reduce the number of the volume points is developed in Section 4, which is the central contribution of the current work; numerical results are presented in Section 5 to illustrate the efficiency of the new algorithm provided by this paper; finally, conclusions are given in Section 6. 2. Radial basis function mesh motion A radial basis function is a real-valued function whose value depends only on the distance to a support point. The general expression of the RBF is

φ = φ(r − ri ),

(1)

where ri is coordinate of the support point, r is the position where the function is evaluated. Usually, the norm · is Euclidean distance. There are three categories of the RBF: global, local and compact. Global bases are always non-zero and grow moving away from the support point. Local bases are also always non-zero but decay moving from the origin. Compact bases decay moving from the origin like the local functions but reach zero at a finite distance. As a compact basis, the Wendland’s C 2 function, which is given by

 φ(η) =

(1 − η)4 (4η + 1), 0 ≤ η < 1 , 0, η≥1

(2)

L. Xie, H. Liu / Journal of Computational Physics 348 (2017) 401–415

403

is chosen as the RBF in the current work due to its satisfactory performances [13]. A weighted sum of N s radial basis functions is used to approximate a given function, which is expressed by

I (r) =

Ns 

αi φ(

r − r i 

i =1

R

(3)

),

where I (r) is the function to be evaluated at location r, N s is the number of the support points, αi is the weight corresponding to the i-th support point, and R is the support radius. For mesh deformation problems, I (r) indicates the deformation of the  volume node,the support  points are the  mesh nodes located on the moving surface. In the current paper, we use V = v 1 , v 2 , . . . , v N v and S = s1 , s2 , . . . , s N s represent the set of the volume and surface nodes respectively. The weights αi are solved by requiring the following condition must be satisfied:

Xs = Max , Ys = Ma y ,

(4)

Zs = Maz , where

T  Xs = xs1 , xs2 , · · · , xs N s

(5)

is a vector contains the displacements of the x component of the surface mesh nodes (Ys , Zs could be defined analogously),

ax =



αx,s1 , αx,s2 , · · · , αx,s N s

T

(6)

is a vector contains the coefficients corresponding to the displacements of the x direction (a y and az are analogous). The elements of the matrix M are defined by

M j ,i = φ(

  rsj − rsi  R

),

(7)

which indicates the basis function evaluated on the distance between surface points s j and si with a support radius of R. Once the coefficients have been computed, the displacements of the volume nodes could be calculated by Eqn. (3). With the definition of a vector X v = [x v1 , x v2 , . . . , x v N v ] T , which represents displacements of all volume nodes, the updating of the volume nodes could be written as

X v = Hax ,

(8)

where the element of the matrix H is defined by

   r v j − rsi  ),  R

H j ,i = φ( 

(9)

which represents the RBF function evaluated on the distance between the volume point v j and the surface point si . Y v and Z v could be defined and computed in a similar way. To simplify formulations, Eqn. (4) and Eqn. (8) could respectively be expressed in the following general form:

 S = Ma,

(10)

 V = Ha.

(11)

3. Surface point selection procedure 3.1. Greedy algorithm Although the RBF algorithm produces a deformation mesh with high quality, the computational cost could be too high to prevent the implementation of the procedure on a larger mesh. In the description of last section, M has a size of N s × N s , while H has a size of N v × N s . Namely, the CPU cost associated with solving the linear algebra equations (Eqn. (4)) scales with N s3 , while which of the mesh motion (Eqn. (8)) scales with N v × N s . Therefore, the reduction of the number of the surface points N s is essential to tackle both the computational challenges of solving and updating the interpolation. In order to achieve this goal, Rendall and Allen [13,14] provided a point selection procedure based on the greedy algorithm. Their method works by starting with a RBF interpolation system with only one single point and then selecting surface point where the error of the current interpolation system is the largest. When a new point is added to the selected active points list, a new coefficient is introduced to the RBF system which will correct the position of that point exactly. Then errors on all surface points will be evaluated to find the point with the maximal error. This point will be added to the selected active

404

L. Xie, H. Liu / Journal of Computational Physics 348 (2017) 401–415

points list. The above greedy loop will be executed some times until the errors of all surface points are sufficiently small. This algorithm reduces the number of the surface points involved into the RBF interpolation method as far as possible while an acceptable requirement of accuracy is meet. The Gauss elimination method could be used to solve the linear system of equations Eqn. (4). For a linear system of equations Ax = b, the algorithm could be described by the following pseudocode (before the procedure, the vector of b is stored in the vector of x) Algorithm 1 Forward elimination. 1: for k = 1 → N − 1 do 2: for j = k + 1 → N do 3: f ← A j ,k / A k,k 4: for i = k → N do 5: A j ,i ← A j ,i − f · A k ,i 6: end for 7: x j ← x j − f · xk 8: end for 9: end for

After the elimination procedure, the matrix A is an upper triangular matrix. Then the unknowns could be computed by the following back substitution process: Algorithm 2 Back substitution. 1: for j = N → 1 do 2: for i = N → j + 1 do 3: x j ← x j − xi · A j , i 4: end for 5: x j ← x j / A j, j 6: end for

The elimination and substitution procedure have an order of N 3 and N 2 respectively. Since N s points have been selected finally, the equations Eqn. (4) will be solved N s times, the total computational cost has an order of N s4 . However, when a new point is added to the collection, the matrix only add one new row and column, and the original matrix becomes a sub matrix of the new matrix. Therefore, it is not need to start the elimination process from the first elements of the matrix. On the contrary, only the newly added row and column need to be eliminated. Then, the totally computational cost of the elimination has an order of N s3 instead of N s4 . Although the back substitution procedure should be repeated throughout, since the cost of it has an order of N s2 , the total CPU cost is also an order of N s3 . This technique has been reported recently by Selim and his co-authors [16]. More details could be found in their article [16]. 3.2. Multi-level subspace point selection procedure Although the above points select procedure improves the efficiency of the interpolation especially for the updating step, this procedure itself requires much CPU cost because when a point is added to the active points list, the coefficients of these points should be solved repeatedly. Since the solving of the linear algebra equations requires N s3 computational cost, the total CPU requirement of the select procedure has an order of N s4 (for the standard version). To improve the efficiency of the point selection procedure described above, Wang and his co-workers [17] developed a multi-level subspace radial basis function interpolation algorithm. The key point of the method is setting the error of last interpolation step as the object of current step, and restricting the number of the points selected by each level. Since the computational cost of calculating the weights is N s4 , this subspace point selection procedure could save the CPU requirement. The details of the algorithm is described below. In the first step, the interpolation target is the displacements of all surface points, namely

 S ( 0) =  S .

(12)

(0)

N s points are selected from all surface points to construct the RBF interpolation system. The vector of coefficients is defined by ( 0)

( 0)

( 0)

a(0) = [αs1 , αs2 , . . . , αs

(0)

Ns

]T .

(13)

The residual of the interpolation using this set of points is

 S (1) =  S (0) − M(0) a(0) .

(14)

L. Xie, H. Liu / Journal of Computational Physics 348 (2017) 401–415

405

(0)

Theoretically, the interpolation result is exact on these N s selected points but not exact on those points which are not (0) selected. Namely, the N s elements in  S (1) corresponding to these selected points are zero. The influence of these surface points on the motion of volume nodes could be described by

 V (0) = H(0) a(0) ,

(15) (0)

where H(0) has a size of N s × N v . In the next step, the interpolation target is switch to the residual of the previous step  S (1) , then a new interpolation problem is given by

 S (1) = M(1) a(1) .

(16) (1)

With the same greedy algorithm described above, N s points effected by these points is described by

points will be selected for this step. Then the updating of volume

 V (1) = H(1) a(1) ,

(17) (1)

where H(1) has a size of N s

× Nv .

(k)

The above algorithm could be carried some times to make the error to be small enough. For the k-th level, N s could be set to be a small value M. If this loop has been executed N l times, the computational cost for constructing a RBF interpolation with Nl × M supporting points is on the order of N l × M 4 , whereas the order of computational costs for the traditional data reduction algorithm to obtain the same supporting points is ( N l × M )4 [17]. Therefore, the multi-level method could save the CPU requirement of the points selection procedure. Finally, the approximation of function values of surface points and volume points is respectively given by i = N l −1

S =



i = N l −1

 S (i ) =

i =0

 i =0

M(i ) a(i ) ,

(18)

i =0

i = N l −1

V =



i = N l −1

 V (i ) =



H(i ) a(i ) .

(19)

i =0

i =N −1

(i )

i =N −1

(i )

l Therefore, the CPU cost of mesh updating step scales with N v × i =0 l N s . Since N s is often larger than N s of i =0 the single-level greedy method, the multi-level process is less efficient than the single-level method on the mesh motion procedure.

4. Volume points reduction algorithm As described before, although the multi-level method improves the efficiency of surface points selection procedure, it requires much more CPU cost on the motion of volume nodes. We believe that the reduction of volume nodes need to be moved on each level is helpful to alleviate this problem. Although the reduction of N v only influences the CPU cost of grid updating, this idea is also appealing because for many cases the surface points selection procedure could be placed outside the simulation loop, though this would implying choosing the points before the simulation started using a benchmark deformation [13]. For FSI problems where a model analysis is employed for modeling the structure, the weights even could be solved just one time before the iteration loop and used for all step during the simulation. Consider how to achieve this task. One straightforward choice is reducing the support radius R when a compact basis is used. However, the practical restriction on the support radius is that it should be larger than the biggest motion of any surface point, usually by a factor of several times [13]. Large value of the supporting radius yields a good approximation order while small value leads to much more error. Therefore, limiting the deformation region by the reduction of the support radius is not practical. Another way [18] is adding some support points around the deformation surface and setting the displacements of these points to be zero, as depicted by Fig. 1. Then the deformation of mesh nodes outside these points is zero. This strategy was first used by Andreas [18] to locally restrict mesh deformation to the vicinity of the moving component and leave the surfaces of other components unaffected. This technique was adopted by the previous research of the first author [22]. However, it is too complicated to implement for complex configuration because the arrangement of the additional points could not be easily coded automatically. It also increases the number of the support points of the RBF system which could hamper the performance of the volume points reduction algorithm. Another drawback is that the additional points could not be placed very closely to the moving surface, otherwise singularity could be introduced into the interpolation matrix. Since most of grid nodes are located near the wall boundary for viscous mesh, the reduction of the CPU cost is not obvious for large scale viscous mesh as illustrated by our previous work [22]. Therefore, we do not adopt this technique in the present paper.

406

L. Xie, H. Liu / Journal of Computational Physics 348 (2017) 401–415

Fig. 1. Illustration of restricting technique by adding some support points around the deformation surface [18].

Instead, the current article introduces a wall distance based restricting function which has following expression:

ψ = ψ(

d(r) D

(20)

),

where d(r) indicates the distance from the position r to the motion surface, D is the support distance. The function should decay from ψ = 1 for the motion surface to ψ = 0 when the distance is larger than the support distance. The present work d(r) chooses the following expression (with the definition of ξ = D )



ψ(ξ ) =

1 − ξ, 0,

0≤ξ <1 . ξ ≥1

(21)

The above function could be added to the RBF interpolation, then the following restricted formulation of the RBF system could be obtained:

I (r) = ψ(

Ns d(r) 

D

)

i =1

αi φ(

r − r i  R

).

(22)

When the above formulation is employed, the interpolation value is zero in the region where d(r) > D. Therefore, only the volume nodes whose wall distance is smaller than D are needed to be moved. Namely, the number of the volume nodes needed to be moved N v could be decreased. In other word, the reduction of the volume nodes is independent of the choice of the support radius with the use of the above wall distance based restricting function. Since ψ(ξ ) = 1 for ξ = 0, Eqn. (22) recovers to Eqn. (3) on the moving surface. Therefore, the surface points selection procedure could be used without any difference for the restricted and un-restricted RBF interpolation procedure. But the volume deformation region could be restricted to the area where d(r) < D with the above expression. For small deformation, D could be set smaller enough to reduce the number of the volume nodes need to be moved. However, for larger deformation, the motion surface would cross the deformation region if D in Eqn. (22) is set to be smaller than the largest displacement occurred on the surface. To reduce the number of the space nodes to be moved for larger deformation, the current paper provides an algorithm which couples the above restricting process and the multi-level subspace interpolation procedure. The details of this method are given below. At first, the interpolation target is set to be the displacements of the surface points, namely

 S ( 0) =  S .

(23) (0)

The tolerance of the interpolation is set to be a value e (0) which is relatively large. The greedy method is used to select N s surface points to make the interpolation error smaller than the tolerance e (0) . Since the tolerance is relatively large, only a few surface points will be chosen. To support large deformation, a large value of support distance D (0) is needed here to prevent the deforming mesh to be overlap. In the present work, D (0) = k · ( S (0) )max is used, where ( S (0) )max is the maximal deformation on the motion surface, the free parameter k is a factor to provide an acceptable volume deformation region. Larger value of k admits larger deformation but leads to bigger N v , while the smaller value has adverse influence. The influence of k will be illustrated by example. In this step, only the volume nodes with a positive ψ(d(r)/ D (0) ) will be (0) influenced by the radial basis functions. The number of these volume points is denoted as N v . The vector of deformation ( 0) ˜ on these nodes is denoted as V which is computed by

˜ (0) a(0) ,  V˜ (0) = H ˜ (0) has a size of N (v0) × N s(0) . where H

(24)

L. Xie, H. Liu / Journal of Computational Physics 348 (2017) 401–415

407

Fig. 2. Original (left) and deformed (right) meshes around the NACA0012 airfoil.

Once the mesh has been updated by the first level RBF-based interpolation, the error of the surface points is set to be the interpolation target of the second level, namely

 S (1) =  S (0) − M(0) a(0) .

(25)

The tolerance e (1) of this step is set to be smaller than that of the previous step e (0) , which indicates that more surface (1) (0) points will be selected by this step, namely N s > N s . However, the support distance D (1) at this step could be set smaller ( 1) than the previous one because we choose D = k · ( S (1) )max .  S (1) < e (0) implies  S (1)   S (0) , namely D (1)  D (0) . (1) (1) Only the nodes with a positive ψ(d(r)/ D (1) ) need to be moved, the number of which is N v . Since D (1)  D (0) , N v is (0) relatively smaller than N v . The above multi-level restricting process could be executed more times to obtain high accuracy of interpolation. When (k) (k) (k) (k) more steps involved, N s will become bigger but N v will become smaller. Therefore, the product N v × N s could be preserved to be smaller than the single-level interpolation algorithm. Although multiple steps will be carried, we believe that the total computational cost will be still saved significantly by our work, while much more accuracy will be achieved. Finally, the deformation of volume points could still be expressed by Eqn. (19), but the matrix H(i ) is non-longer density ˜ Therefore, the CPU requirement of mesh updating step matrix. The size of non-zero elements in H is the same as that of H. i =Nl −1 (i ) (i ) (i ) by this restricting procedure scales with N × N . Since N v < N v , the current algorithm is more efficiency than v s i =0 Wang’s approach. We also believe that the CPU cost of the present algorithm is smaller than that of the single-level method. This declaration will be demonstrated by some numerical examples. It should be pointed out that the nearest point correction process, as defined by Rendall and Allen [15] and adopted by Wang [17], might be thought of as the very finest correction level used in the multi-level method here. The multi-level approach here bridges the gap from the coarsest to the finest level in a smooth manner. The provided algorithm involves the wall distance d(r), which is just computed one time before the iteration loop. Therefore, the calculation of the wall distance does not harm the efficiency of the current method. An algorithm based on the alternating digital tree (ADT) [23] is employed by the current work to compute d(r), the computational cost of which has an order of O ( N v × log( N s )). Some numerical examples will be given in the next section to verify the efficiency of the provided algorithm. 5. Numerical results 5.1. Evaluation of the impact of k The influence of the parameter k is illustrated by the deformation of a 2D NACA0012 airfoil. The mesh of the foil is described by Fig. 2, which contains 23020 volume cells and 240 surface elements. The deformation of this airfoil could be described by the following formula:

 y = 0.03 sin(4π x).

(26)

The original and deformed meshes are shown in Fig. 2. To evaluate the influence of the restricting procedure on the quality of the mesh, the following formula is adopted

D Q = Q d − Q o,

(27)

408

L. Xie, H. Liu / Journal of Computational Physics 348 (2017) 401–415

Fig. 3. Impact of the parameter k on the deformation of NACA0012 airfoil. Left: deformed mesh; middle: contour of D Q ; right: histogram of elements with D Q > 1◦ . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

where Q d and Q o are respectively the internal angles of elements of deformed and original meshes. The value of D Q indicates the change of the internal angle of deformed element relative to the un-deformed element. The smaller the value, the better the quality of deformed mesh is. The results of the un-restricted and limited deformed meshes with different value of k are displayed in Fig. 3. Pictures on the left, middle and right show respectively the distorted meshes, the contour of D Q and the stochastic histogram of the deformed elements (only cells with D Q > 1◦ are counted). From top to bottom, results of un-limited mesh, limited mesh with k = 5.0, 2.0, 1.2 are shown respectively. It is obvious that only elements near the foil are deformed when the restricting strategy is employed. The quality of restricted meshes with k = 5.0 and k = 2.0 is very close to that of the un-restricted mesh, because both the maximum values of D Q of these three meshes are about 22◦ . However, the quality of restricted mesh is not accepted for k = 1.2 since the maximum of D Q is near 50◦ . The histogram and the contour of D Q state clearly that the number of deformed elements of restricted mesh is smaller than that of the un-restricted deformed mesh. These results illustrate that smaller value of k, the less the deformed cells, and the worse the quality. Experience shows that k = 5 is acceptable for both efficiency and quality, even for a complex configuration. 5.2. Evaluation the new procedure on the deformation of 2D NACA0012 airfoil The same mesh and deformation formula of the previous case is employed to test and compare the RBF mesh deformation method with and without the volume points reduction algorithm. At first, we illustrate the performance of the different implementation of greedy algorithm on the surface points selection procedure by Fig. 4. This figure shows the number of selected surface points and the CPU requirement of the different greedy algorithm. It is obvious that Wang’s method selects more surface points although it requires less computational time compared with the standard greedy algorithm provided by Rendall and Allen. This result agrees with the original literature [17]. However, the method proposed recently by Selim [16] requires the least computational time and selects the same number of points as the original method. Therefore, Selim’s algorithm is adopted in our multi-level restricting method to selected the surface points at each level. Six level interpolations have been used for the multi-level restricting procedure for this problem. The interpolation tolerance of each level is given in Table 1. The convergence rate of displacement error in terms of the number of the surface points selected is shown in Fig. 5. It is clear that much more surface points needed for lower tolerance. However, as depicted by Fig. 6, the region of the volume points needed to be moved is restricted hierarchically. Namely, fewer volume grid nodes are needed to be moved for smaller tolerance. Fig. 7 shows the D Q -based quality of deformed mesh for the

L. Xie, H. Liu / Journal of Computational Physics 348 (2017) 401–415

409

Fig. 4. The number of selected surface points and CPU cost of the greedy algorithm for the NACA0012 case.

Table 1 Comparison of efficiency of single- and multi-level restricting interpolation algorithm on the mesh updating procedure for the NACA0012 case. Method

N v × Ns

Level

k

Relative tolerance

Ns

Nv

1 2 3 4 5 6 total

Non-restricting 5.0 5.0 5.0 5.0 5.0 –

1e-1 1e-2 1e-3 1e-4 1e-5 1e-7 1e-7

11 22 39 76 136 182 –

12986 6504 3559 1686 481 240 –

142846 143088 138801 128136 65416 43680 661967

3.394e-3 3.389e-3 3.435e-3 3.063e-3 1.571e-3 1.065e-3 1.5917e-2

Wang’s method ( M = 5) Wang’s method ( M = 10)

Total Total

– –

1e-7 1e-7

221 219

20289 20289

4483869 4443291

9.5920e-2 9.5115e-2

Single-level

–

–

1e-7

205

20289

4159245

8.7163e-2

Current

CPU time (s)

Fig. 5. Convergence history of the multi-level surface points select procedure for the NACA0012 case.

new and old algorithm. The maximum values of D Q of restricted and un-restricted deformed mesh are respectively 24o and 210 . Evidently, the performance of two methods on the mesh quality is very close. The comparison of efficiency of the proposed algorithm and the original procedure is given in Table 1. The total CPU cost of the multi-level restricting interpolation is about 1.59e-2s, whereas which of the standard single-level interpolation is 8.72e-2s. Namely, the cost of the current algorithm is 18 percentage of that of the standard procedure. However, Wang’s method requires 9.59e-2s and 9.51e-2s respectively for M = 5 and M = 10 which are little larger than that of the standard greedy approach. To illustrate the influence of the empirical parameter k in the multi-level restricting algorithm, Table 2 shows the comparison of efficiency and mesh quality with difference k on the mesh updating procedure for the NACA0012 case. These calculations use the same parameters as Table 1, except that the k values are different. The mesh quality is evaluated by the maximum value of D Q . It is clear that larger k leads to better mesh quality but more CPU cost while smaller k leads

410

L. Xie, H. Liu / Journal of Computational Physics 348 (2017) 401–415

Fig. 6. Deformed mesh around the NACA0012 airfoil by multi-level restricting method of 1st step (top left), 2nd step (top right), 4th step (bottom left) and 6th step (bottom right). Blue nodes: selected surface points; red region: volume nodes needed to be moved. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 7. Quality of deformed mesh of the NACA0012 case. Left: the un-restricted mesh; right: the multi-level restricted mesh (k = 5). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

to less computational time and worse results. Comprehensively considering the computational efficiency and mesh quality, our experiences show that k = 5 is appropriate for most numerical examples. Although for this problem, the CPU resume is small due the small size of the mesh, the proposed procedure has shown great potential to save the CPU requirements of the mesh updating process. We will use some 3D problems to illustrate the advantage of the present technique in next sub-sections.

L. Xie, H. Liu / Journal of Computational Physics 348 (2017) 401–415

411

Table 2 Comparison of efficiency and mesh quality with difference K on the mesh updating procedure for the NACA0012 case. k

CPU time (s)

D Q max

2 3 4 5 6 7 8 9 10

1.332e-2 1.432e-2 1.449e-2 1.592e-2 1.672e-2 1.789e-2 1.784e-2 1.811e-2 1.841e-2

36.481 26.667 25.264 24.418 23.856 23.460 23.128 22.905 22.694

Non-restricting

8.716e-2

21.150

Fig. 8. Original and deformed mesh around the ONERA-M6 wing.

5.3. Deformation of 3D ONERA-M6 wing Here we choose deformation of a 3D ONERA-M6 wing to show the behavior of the present method on the large deformation problem. The deformation is determined by the following expression

z·π   y = 0.25bz 1 − cos( ) , b

(28)

where b is the band width of the ONERA-M6 wing. The mesh used in this problem is a hybrid mesh which contains 627446 grid nodes. These nodes build 1295207 cells including 797936 tetrahedrons, 474720 hexahedrons and 22551 pyramids. The original and deformed meshes are drawn in Fig. 8, from which we could observe the large deformation of the wing clearly. To validate the present algorithm for viscous simulation, the height of the first mesh cell off the wall is set to be 1.5e-6(m) to achieve a target y+ value of 1.0 for viscous simulation (Re = 1.86e7 based on the chord). Both the multi-level restricting and single-level standard interpolation are adopted to perform the volume grid deformation to compare the efficiency. Seven layers are adopted for the multi-level restricting method. The final tolerances of the two methods are both 1e-8. The convergence history versus the number of selected points of the multi-level algorithm is depicted in Fig. 9. The selected surface points as well as the interpolation error on the surface points are displayed in Fig. 10. From these figures, the performance of the greedy procedure is shown clearly for each level. For each step of multilevel restricting method as well as the standard single-level algorithm and Wang’s multi-level approach, the tolerance, the number of the selected surface points, the number of the moved volume nodes and the CPU cost are given in Table 3. The total computational cost of the multi-level restricting method is about 3.273s, whereas which of the single-level procedure is about 20.57s. The CPU requirement of the former is only about 16 percentage of the later for this problem. However, Wang’s method requires 21s which is also little larger than that of the standard greedy approach. This case indicates the good performance of the current algorithm, namely, the CPU cost of the volume mesh updating step is reduced dramatically. Fig. 11 shows a comparison of the quality of the deformed meshes obtained respectively by the standard and the current algorithm. The maximum values of D Q on this problem are respectively 53.1o and 53.4o for un-restricted and restricted mesh. It is clearly that the performance of the proposed algorithm on the quality is very close to full RBF motion.

412

L. Xie, H. Liu / Journal of Computational Physics 348 (2017) 401–415

Fig. 9. Convergence history of the multi-level surface points select procedure for the M6 case.

Fig. 10. Interpolation error on the surface of the ONERA M6 wing by multi-level restricting method for 1st step (top left), 2nd step (top right), 4th step (bottom left) and 7th step (bottom right). Red nodes: selected surface points. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

L. Xie, H. Liu / Journal of Computational Physics 348 (2017) 401–415

413

Table 3 Comparison of efficiency of single- and multi-level restricting interpolation algorithm for the M6 case on the mesh updating procedure. Method

N v × Ns

Level

k

Relative tolerance

Ns

Nv

1 2 3 4 5 6 7 total

Non-restricting 5.0 5.0 5.0 5.0 5.0 5.0 –

1e-2 1e-3 1e-4 1e-5 1e-6 1e-7 1e-8 1e-8

11 39 149 541 1402 2271 1874 –

627446 326076 198907 121895 65329 23205 11167 –

6901906 12716964 29637143 65945195 91591258 52698555 20926958 280417979

Wang’s method ( M = 5) Wang’s method ( M = 10)

total Total

– –

1e-8 1e-8

2837 2841

627446 627446

1780064302 1782574086

20.87160 20.91104

Single-level

–

–

1e-8

2795

627446

1753711570

20.56931

Current

CPU time (s) 0.080975 0.141059 0.327682 0.778720 1.077351 0.620005 0.246853 3.272645

Fig. 11. Quality of deformed mesh of the M6 case. Left: the un-restricted mesh; right: the multi-level restricted mesh (k = 5). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

5.4. Deformation of 3D DLR-F6 configuration The DLF-F6 model is adopted to evaluate the performance of the current method on complex configuration. As described by the following equations, the deformation is a blend of a bend and a twist: z  y = δ0 sin(− 2b π ); z φ = φ0 (exp(− 2b ln(2)) − 1);

δ0 = 0.5c ; φ0 = 30◦ ;

(29)

where b and c are respectively the span size and mean aerodynamic chord length. The hybrid unstructured mesh is generated for this configuration, which is depicted in Fig. 12. The mesh contains 2257483 nodes and 3564918 cells including 1732295 hexahedrons, 1074650 tetrahedrons, 475568 pyramids and 282405 prisms. The normal height of the first element off the wall is taken 2.6e-3(mm) so that the target y+ equals to 1.0. The efficiency of the new multi-level method, the standard single-level approach and Wang’s multi-level procedure is compared by Table 4. Parameters of the computation are also given in the table. Obviously, the efficiency of the RBF method has been improved dramatically by the proposed method. The CPU cost of the present approach is only about 21.4 percentage of that of the single-level algorithm. The performance of Wang’s multi-level method is very close to that of the standard approach. The comparison of the quality based on D Q is shown by Fig. 13. The maximum values of D Q on this problem are respectively 39o and 34o for restricted and un-restricted mesh. It is evident that the quality of the two deformed meshes is very close. 6. Conclusions To save the computational cost of the mesh updating of the RBF based mesh deformation method, the current paper provides a multi-level based volume nodes restricting algorithm. To efficiently limit the volume region driven by the motion surface, a wall distance based restricting function is introduced. For each level, the target of the interpolation problem is the error produced by the previous step. Since the error is decreased, the updating volume region could be reduced hierarchically.

414

L. Xie, H. Liu / Journal of Computational Physics 348 (2017) 401–415

Fig. 12. The original hybrid unstructured mesh of the DLR-F6 configuration for the viscous simulation.

Table 4 Comparison of efficiency of single- and multi-level restricting interpolation algorithm for the DLF-F6 configuration on the mesh updating procedure. Level

k

Relative tolerance

Ns

Nv

N v × Ns

CPU time (s)

1 2 3 4 5 Total

Non-restricting 5.0 5.0 5.0 5.0 –

1e-3 1e-4 1e-5 1e-6 1e-8 1e-8

57 186 468 946 1966 –

2257483 1427610 655566 164039 82008 –

128676531 265535460 306804888 155180894 161227728 1017425501

1.402 2.916 3.539 2.017 1.887 11.761

Wang’s method ( M = 5) Wang’s method ( M = 10)

Total Total

– –

1e-8 1e-8

2030 2053

2257483 2257483

4582690490 4634612599

55.90 56.78

Single-level

–

–

1e-8

1986

2257483

4483361238

54.83

Method

Current

Fig. 13. Contour of D Q of deformed mesh of the F6 case (unit: degree). Left: the un-restricted mesh; right: the multi-level restricted mesh (k = 5). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

To demonstrate the efficiency of the proposed method for unstructured hybrid mesh, typical deformation cases of the 2D NACA0012 airfoil, 3D ONERA-M6 wing and DLF-F6 configuration are presented. The large deformation is shown to be supported by the current algorithm, surface accuracy has been improved while the computational cost is reduced remarkably for the mesh updating procedure. Acknowledgements This article is supported by Center for HPC Shanghai Jiao Tong University. References [1] T.C.S. Rendall, C.B. Allen, Fluid-structure interpolation and mesh motion using radial basis functions, Int. J. Numer. Methods Eng. 74 (10) (2008) 1519–1559. [2] C.B. Allen, T.C.S. Rendall, A unified approach to CFD-CSD interpolation and mesh motion using radial basis functions, in: 25th Applied Aerodynamics Conference, AIAA Paper No. AIAA-2007-3804, Miami, FL, 2007. [3] R. Ahrem, A. Beckert, H. Wendland, A new multivariate interpolation method for large-scale spatial coupling problems in aeroelasticity, in: Proceedings to Int Forum on Aeroelasticity & Structural Dynamics, 2001.

L. Xie, H. Liu / Journal of Computational Physics 348 (2017) 401–415

415

[4] A.M. Morris, C.B. Allen, T.C.S. Rendall, CFD-based optimization of airfoils using radial basis functions for domain element parameterization and mesh deformation, Int. J. Numer. Methods Fluids 58 (8) (2008) 827–860. [5] W.L. Kong, H. Liu, Development and theoretical analysis of an aircraft supercooled icing model, J. Aircr. 51 (3) (2014) 975–986. [6] J.T. Batina, Unsteady Euler algorithm with unstructured dynamic mesh for complex-aircraft aerodynamic analysis, AIAA J. 29 (3) (1991) 327–333. [7] C. Farhat, C. Degand, B. Koobus, M. Lesoinne, Torsional springs for two-dimensional dynamic unstructured fluid meshes, Comput. Methods Appl. Mech. Eng. 163 (1–4) (1998) 231–245. [8] S.H. Huo, F.S. Wang, W.Z. Yan, Z.F. Yue, Layered elastic solid method for the generation of unstructured dynamic mesh, Finite Elem. Anal. Des. 46 (10) (2010) 949–955. [9] X. Liu, N. Qin, H. Xia, Fast dynamic grid deformation based on Delaunay graph mapping, J. Comput. Phys. 211 (2) (2006) 405–423. [10] A. Beckert, H. Wendland, Multivariate interpolation for fluid-structure-interaction problems using radial basis functions, Aerosp. Sci. Technol. 5 (2) (2001) 125–134. [11] A.D. Boer, M.S. van der Schoot, H. Bijl, Mesh deformation based on radial basis function interpolation, Comput. Struct. 85 (11) (2007) 784–795. [12] S. Jakobsson, O. Amoignon, Mesh deformation using radial basis functions for gradient-based aerodynamic shape optimization, Comput. Fluids 36 (6) (2007) 1119–1136. [13] T.C.S. Rendall, C.B. Allen, Efficient mesh motion using radial basis functions with data reduction algorithms, J. Comput. Phys. 228 (17) (2009) 6231–6249. [14] T.C.S. Rendall, C.B. Allen, Reduced surface point selection options for efficient mesh deformation using radial basis functions, J. Comput. Phys. 229 (8) (2010) 2810–2820. [15] T.C.S. Rendall, C.B. Allen, Parallel efficient mesh motion using radial basis functions with application to multi-bladed rotors, Int. J. Numer. Methods Eng. 81 (1) (2010) 89–105. [16] M.M. Selim, R.P. Koomullil, A.S. Shehata, Incremental approach for radial basis functions mesh deformation with greedy algorithm, J. Comput. Phys. 340 (2017) 556–574. [17] G. Wang, H.M. Haris, Z.Y. Ye, J.D. Lee, Improved point selection method for hybrid-unstructured mesh deformation using radial basis functions, AIAA J. 53 (4) (2015) 1016–1025. [18] A.K. Michler, Aircraft control surface deflection using RBF-based mesh deformation, Int. J. Numer. Methods Eng. 88 (10) (2011) 986–1007. [19] T. Gillebaart, A.H. van Zuijlen, H. Bijl, Radial basis function mesh deformation including surface orthogonality, in: 54th AIAA Aerospace Sciences Meeting, AIAA SciTech, San Diego, California, USA, 4–8 January 2016. [20] Y. Wang, N. Qin, N. Zhao, Delaunay graph and radial basis function for fast quality mesh deformation, J. Comput. Phys. 294 (2015) 149–172. [21] H. Wendland, Scattered Data Approximation, Cambridge Monogr. Appl. Comput. Math., Cambridge University Press, Cambridge, UK, 2005. [22] L. Xie, M. Xu, B. Zhang, et al., Space points reduction in grid deforming method based on radial basis functions, J. Vib. Shock 32 (10) (2013) 141–146 (in Chinese, with English abstract). [23] D.A. Boger, Efficient method for calculating wall proximity, AIAA J. 39 (12) (2001) 2404–2406.