A Meshless Solution Method for Unsteady Flow with Moving Boundary

Abstract

Using the concept of overlapping mesh method for reference, a new method called as Overlapping Clouds of Points Method (OCPM) is firstly proposed to simulate unsteady flow with moving boundary problems based on meshless method. Firstly, a set of static background discrete points is generated in the whole calculation zone. Secondly, moving discrete points are created around moving body. According to the initial position of moving object in the flow field, the two sets of discrete points can be overlapped. With the motion of moving objects in the calculation field, moving discrete points around the moving body will inherently move. The exchange of flow field information between static points and moving points is realized by the solution of the clouds of points made up of static and moving discrete points using weighted meshless method nearby overlapping boundary. Four cases including piston problem, NACA0012 airfoil vibration flow around a moving sphere in supersonic and multibody separation are given to verify accuracy and practicability of OCPM. The numerical results agree well with exact solution and experimental results, which shows that the proposed OCPM can be applied to the simulation of unsteady flow problem.

1. Introduction

The unsteady flow with moving objects is complex in reality. Simulating complex flow with moving boundary, such as free surface problem, multibody separation, aeroelastic problem, forced vibration, and structural coupling, is a difficulty in the computational fluid dynamics. In mesh methods, the difficulty in disposing the unsteady, real, and complex numerical simulation with moving boundary is how to form computational meshes legitimately and quickly and how to make the mesh match well with the movement of objects so as not to induce large deformation. At the same time, suitable boundary conditions and solver need to be chosen carefully [1–3]. Under the mesh framework, such problems can be disposed by mesh deformation method, mesh deformation plus partial, all flow field remeshing method, and overlapping mesh method. The first two methods both ask mesh for strong deform ability and good mesh quality after deformation, which is necessary to maintain numerical accuracy of solutions and expand the use scope of distortion method. Therefore, it needs a great deal of computational effort to guarantee that serious distortion would not happen in the mesh. The last method needs to exchange complex flow field information among the mesh.

Different from mesh algorithm, meshless method [4] is a numerical algorithm that only needs nodes information and does not need to connect nodes into grid cell. In the history of meshless method, it is mainly employed to solve the fracture problems of solid mechanics. By the 1990s [5], it began to be applied to computational fluid dynamics. Large numbers of researches are made by Liu and Su [6], Morinishi [7–9], Anandhanarayanan [10, 11], and Katz and Jameson [12–14] later. The numerical simulation of unsteady flow with moving boundary in meshless method only involves the movement of the discrete points, but not restriction of stretch or distortion of mesh. Compared with mesh methods, meshless method has more advantages of geometric flexibility in dealing with complex aerodynamic calculation for complex configuration and complex flow process. Numerical simulation of unsteady flow based on meshless method has been studied by Morinishi [15], Anandhanarayanan [10, 11], Wang et al. [16], Zhou and Xu [17, 18], Wang et al. [19, 20], and Pu and Chen [21].

In this paper, a new method, called Overlapping Clouds of Points Method (OCPM), is proposed to dispose moving boundary problems based on Batina's meshless method and using the concept of overlapping grid method under mesh framework. The method is verified by piston problem with analytical solution and applied to numerical simulation of NACA0012 airfoil with pitching vibration and flow around a sphere in supersonic, multibody separation.

2. The Basic Principle of Least Squares Meshless Method in CFD

Meshless method in CFD mainly involves generation of discrete points, searching satellite points, and fitting spatial derivative. Using the concept of the clouds of points (Figure 1), the governing partial differential equations (PDEs) are solved directly by performing local least squares curve fits in each cloud of points. After that, the resulting curve-fit equations are analytically differentiated to approximate the derivatives of the PDEs. In order to strengthen the flexibility to choose satellite points so as to process flow problem with complex configuration, the weighted function [22] is introduced in most meshless methods. The values of weight function range from 0 to 1 in compacted support domain. And the function value decreased with increasing of distance between satellite points and central point.

Figure 1:

Cloud of Points.

2.1. Searching Satellite Points

Meshless method only has node information. The space derivation of central point is fitted by flow field values of its adjacent point, which is called satellite point of central point. Satellite points are determined in an automatic searching way—compacted support domain radius r₀ is firstly determined (the size varies with the number of needed satellite points). Suppose there are some points in discrete area and the number of a central point is i (only with several points’ number in the figure) as shown in Figure 1. If the distance r between a certain point and the point i is less than or equal to r₀, this point would be regarded as a satellite point of the central point.

The points in discrete area are arranged out of order. In order to determine the satellite points of the central point i, the central point i needs to compare the distance from all other points to the central point with r₀. If we regard the comparison between the distance from a point to the central point and r₀ as a calculation, the determination of all satellite points in discrete space requires calculation n × (n − 1) times. As a result, when discrete space is a three-dimensional space, the number of discrete points of the region will be up to more than hundred thousands, which is time consuming in searching satellite points. Furthermore, if we double the number of points, the number of calculations will increase by about three times. In order to speed up searching satellite points, we will use the method of zoning search. Suppose the points in computational domain are divided into four zones whereas the middle is the four subareas’ boundary as shown in Figure 2. Based on the need of clouds of points structure and searching satellite point, before searching satellite points, coordinate information of hollow points near the overlapping boundary is transferred to adjacent subzones. Then searching satellite points will be done in each subzone and hollow points belonging to adjacent subzones. If computational domain is divided into four zones, the number of computations to determine clouds of point's structure in the whole discrete area is about 4 × n/4 × (n/4 − 1), making it reduced by three-quarters.

Figure 2:

Zoning sketch.

2.2. Governing Equations and Numerical Method

Axisymmetrical Euler equations in differential type are shown in a Cartesian coordinate system for

$\frac{\partial \mathbf{U}}{\partial t}+\frac{\partial {\mathbf{F}}_{\mathbf{1}}\left(\mathbf{U}\right)}{\partial x}+\frac{\partial {\mathbf{F}}_{\mathbf{2}}\left(\mathbf{U}\right)}{\partial y}+\mathbf{G}=\mathrm{0}.$

(1)

In (1),

$\begin{array}{c}\mathbf{U}=\left(\begin{array}{c}\rho \\ \rho u\\ \rho v\\ e\end{array}\right),\hspace{1em}\hspace{1em}{\mathbf{F}}_{\mathbf{1}}\left(\mathbf{U}\right)=\left(\begin{array}{c}\rho u\\ \rho u^{\mathrm{2}}+p\\ \rho uv\\ u\left(e+p\right)\end{array}\right),\\ {\mathbf{F}}_{\mathbf{2}}\left(\mathbf{U}\right)=\left(\begin{array}{c}\rho v\\ \rho vu\\ \rho v^{\mathrm{2}}+p\\ v\left(e+p\right)\end{array}\right),\hspace{1em}\hspace{1em}\mathbf{G}=\frac{\mathrm{1}}{y}\left(\begin{array}{c}\rho v\\ \rho uv\\ \rho u{v}^{\mathrm{2}}\\ \left(e+p\right)v\end{array}\right).\end{array}$

(2)

In (2), ρ is fluid density, p is pressure, e is unit volume total internal energy, and u, v are the velocity component of x, y direction, respectively. If computational space is two-dimensional, then G = 0.

When y → 0,

$\begin{array}{c}\frac{v}{y}=\mathrm{0.5}\hspace{0.17em}\left(X_m-\mathrm{1}\right)\frac{\partial u}{\partial x},\\ X_m=\frac{u}{c},\\ c=\sqrt{\frac{\gamma p}{\rho }}.\end{array}$

(3)

Equation of state of gas is

$p=\left(\gamma -\mathrm{1}\right)\rho \left[e-\frac{\mathrm{1}}{\mathrm{2}}\left(u^{\mathrm{2}}+v^{\mathrm{2}}\right)\right].$

(4)

In (4), γ is ratio of specific heat; here it is 1.4.

Space derivative is determined in local cloud of points by quadratic minimal surface approximation as shown in Figure 1. Discrete points k _n (n = 1, 2, 3, 4, 5) are satellite 9 points of ith central point in space. Every point can be expressed as (x, y, f); f is the point's flow function value in local coordinates. That is to say, i is the origin of coordinate and k _n 's coordinate is (Δx _kn , Δy _kn , Δf _kn ). Consider

$\begin{array}{c}\Delta x_{k_n}=x_{k_n}-x_i,\\ \Delta y_{k_n}=y_{k_n}-y_i,\\ \Delta f_{k_n}=f_{k_n}-f_i.\end{array}$

(5)

Suppose function f can be expressed by linear relation as follows:

$f=ax+by.$

(6)

The number of points in the cloud of point n is more than 2 leading to the contradictory equations:

$A\stackrel{⃑}{x}=\stackrel{⃑}{f}.$

(7)

Solving (7) in least square method, we can obtain a and b and then space derivative. The mathematical expression is

${\frac{\partial f}{\partial x}|}_i=a,\hspace{1em}\hspace{1em}{\frac{\partial f}{\partial y}|}_i=b.$

(8)

Referring dissipationmode [19, 20] of unstructured grid method, the second-order and fourth-order artificial dissipations are introduced for stability.

Supposing D(U) is additional dissipation the dissipation in ith point is calculated as follows:

$D_i\left(U\right)={\displaystyle \sum }_{k=\mathrm{1}}^K{d}_{i,k},$

(9)

where

$\begin{array}{c}{d}_{i,k}={\alpha }_{i,k}\hspace{0.17em}{\epsilon }_{i,k}^{\left(\mathrm{2}\right)}\left(U_k-U_i\right)-{\alpha }_{i,k}\hspace{0.17em}{\epsilon }_{i,k}^{\left(\mathrm{4}\right)}\left({\nabla }^{\mathrm{2}}{U}_k-{\nabla }^{\mathrm{2}}{U}_i\right),\\ {\nabla }^{\mathrm{2}}{U}_i={\displaystyle \sum }_{k=\mathrm{1}}^n\left(U_k-U_i\right),\\ {\epsilon }_{i,k}^{\left(\mathrm{2}\right)}=K^{\left(\mathrm{2}\right)}{\upsilon }_{i,k},\\ {\epsilon }_{i,k}^{\left(\mathrm{4}\right)}=\max \left(\mathrm{0},K^{\left(\mathrm{4}\right)}-{\epsilon }_{i,k}^{\left(\mathrm{2}\right)}\right),\\ {\upsilon }_{i,k}=\frac{\left|p_k-p_i\right|}{\left|p_k+p_i\right|},\\ \Delta t_i=\frac{\mathrm{1}}{\left|u\right|+c+\left|v\right|+c},\\ {\alpha }_{i,k}=\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\mathrm{1}}{\Delta t_i}+\frac{\mathrm{1}}{\Delta t_k}\right),\end{array}$

(10)

where K⁽²⁾ and K⁽⁴⁾ are specified by user. Their values are taken as 0.5 and 1/64, respectively, in present work. Semidiscrete Euler equations are shown as follows:

${\frac{dU}{dt}|}_i+a_i+b_i+G_i-D_i\left(U\right)=\mathrm{0}.$

(11)

Fourth-order Runge-Kutta format is applied in the process of time advance.

In order to reduce the effect on the physical solution's accuracy of space derivative because of distance, it is needed to improve the flexibility in choosing satellite points. When fitted by quadratic minimization surface in weight function, all of them have the features of nonnegativity and monotony. Making w(k) weight function, we get function about distance:

$w\left(k\right)=\{\begin{array}{cc}{\left(\frac{d_{ik}}{\stackrel{-}{d}}\right)}^{\alpha },\hfill & d_{ik}<\stackrel{-}{d},\hfill \\ \mathrm{1.0},\hfill & d_{ik}>\stackrel{-}{d},\hfill \end{array}$

(12)

where α numerical range is 1.0 to 2.0. d is reference distance and d _ik is the distance from point k to point i. With the introduced weight function, contradictory equations mentioned above can be expressed in

$WA\stackrel{⃑}{x}=W\stackrel{⃑}{f}.$

(13)

We can get space derivative after weight in least square method.

3. The Overlapping Clouds of Points Method

The fundamental idea of Overlapping Clouds of Points Method is to use two or more sets of independently formed discrete points to cover the whole flow field and form overlapping clouds of points in some specified region. The exchange of flow filed information nearby the boundary of each overlapped region is realized by using a special cell structure in the meshless method, which is called “clouds of points.”

Assume that a certain body A is moving in flow field, as shown in Figure 3. Firstly, background discrete points are generated in zone 1. Secondly, moving discrete points are generated around moving body A in zone 2. According to the initial position of object A in the flow field, the two sets of discrete points can be overlapped. As a result, a part of points in zone 1 are covered by those in zone 2. The attribute of the covered points in zone 1 and the points very close to the overlapping boundary is marked as “nonactive.” The attribute of the rest of the points in zone 1 and the points in zone 2 is marked as “active.” Before the movement of body A, the clouds of points structure in the whole flow field (including zone 1 and zone 2) should be determined. However, the satellite points are searched in the points of which attribute is active. Near the overlapping boundary, some satellite points of the central point which is in zone 1 come from zone 2. Similarly, some satellite points of central point which is in zone 2 come from zone 1. From Figure 4, the points marked with filled square in zone 1 are active points and the points marked with outlined square are the nonactive points covered by zone 2 or very close to the overlapping boundary, and the points marked with hollow circle come from zone 2. The fourth and fifth satellite points of the ith point of zone 1 come from zone 2, while the third and fourth satellite points of jth point of zone 2 come from zone 1. Once the clouds of points structure near overlapping boundary is determined, the exchange of flow field information will be realized between zone 1 and zone 2.

Figure 3:

Overlapping clouds of points.

Figure 4:

Structure of clouds of points near the overlapping boundary.

4. Numerical Example

In order to verify accuracy and practicability of OCPM, four cases including piston problem, NACA0012 airfoil vibration flow around a moving sphere in supersonic and multibody separation, are given.

4.1. The Problems of Piston

Supposing there is a hollow cylinder with a piston on the left, when a piston moves to the right at a constant speed, right-row compression wave will be produced in the right flow field of the piston. But when a piston does constant motion to the left, expansion wave will be produced in the right piston. The two waves’ forms both have analytic solution.

Choosing nondimensional computational area as 1 × 0.2, initial condition of the air in the right cylinder's cavity of the piston is as follows: pressure p = 1.0, density ρ = 1.4, x-velocity component u = 0.0, y-velocity component v = 0.0, and ratio of specific heat γ = 1.4. In Figure 5, it shows the comparison between analytic solution and numerical results that piston does constant stretching movement to the left and compression movement to the right. From the figure, we can see that numerical results match well with analytic solution, which shows that it is available to solve the unsteady flow field problems with moving objects.

Figure 5:

Pressure distribution induced by movement of a piston with different velocities.

4.2. NACA0012 Airfoil Vibration

Airfoil vibration problems are typical moving boundary problems. In this passage, NACA0012 airfoil pitching vibration problems around a quarter of chord are solved by Overlapping Clouds of Points Method based on meshless method. The inner boundary uses solid wall boundary condition, while the other exterior boundary uses nonreflecting boundary condition. Initial condition of the air i is as follows: pressure p = 1.0, density ρ = 1.4, Mach number Ma = 0.755, and the evolution of angle of attack α(t) is

$\alpha \left(t\right)={\alpha }_m+{\alpha }_{\mathrm{0}}\sin \left(\omega t\right).$

(14)

In formula (14), average value α _m of angle of attack is 0.016° and amplitude of vibration α₀ is 2.51°. Nondimensional angular frequency k is ωc/2U_∞, among which c is chord length and U_∞ is the velocity of free stream.

In this case, angular frequency k is 0.0814, the corresponding vibration period T is 51.1185, and aerofoil around a quarter of chord begins to vibrate clockwise from the initial position α _m = 0.016° for two periods. In Figure 6, lift coefficient and pitching moment coefficient change with angle of attack. The numerical result shows a slight difference when compared with the experimental result introduced in [23], which shows that numerical simulation results with moving boundary field by meshless method are acceptable.

Figure 6:

Changes in lift coefficient and pitching moment coefficient with angle of attack.

4.3. Flow around a Moving Sphere

The motion of a sphere in a static flow field is a widely known case of supersonic flow with moving boundary. It is convenient to choose a reference frame in which the sphere is stationary in a uniform coming flow. The initial conditions are velocity Ma = − 3 the pressure p₀ = 1.0, and the density ρ₀ = 1.4. When the sphere moves at velocity Ma = 3 in a static flow, the initial conditions of the air are pressure p₀ = 1.0 and density ρ₀ = 1.4. Figure 7 shows the variation of drag coefficients with time under two cases including the moving sphere and static sphere. The numerical result agrees well with the experimental results given by document [24]. It can be observed from Figure 8 that the evolutions of drag coefficient are consistent with time between two different situations of the static sphere and moving sphere. When t > 0.8, the drag coefficients in two different situations are very close. Figure 8 shows the pressure flood contours at different time for the stationary sphere and moving ball, which reveals the formation process of a shock wave out of body. Comparing isobaric flood contours at different time, we can see that the shape of shock wave out of body is similar as well as the pressure flood contours at the same time in two different situations.

Figure 7:

Comparison of drag coefficients for moving sphere and static sphere.

Figure 8:

The pressure distribution when the sphere moves.

4.4. Multibody Separation

The model in document [25] is employed for this case. The computational geometry with sudden enlargement in its cross section is a shock tube connected to a big chamber, as schematically shown in Figure 10. The shock tube is 33 long and 5 in diameter. The big chamber is 53 long and 20 in diameter. The projectile is 5 long and 5 in diameter. A projectile moving at supersonic speed in a tube, like a piston, drives a precursor shock wave ahead of it. When the precursor shock wave is discharged from the open-end of the shock tube, shock wave diffraction with an associated starting vortex and a jet flow will develop. Later, when the projectile itself moves out of the shock tube, it will interact with this diffracting shock system and the jet flow. Besides this, a bow shock wave will develop in front of the projectile as soon as it moves out of the first shock cell in the jet and will catch up and pass through the precursor shock wave some time later. Shocked flows induced by a supersonic projectile moving from a shock tube into a big chamber are closely related to the ram accelerator research and muzzle flow field.

The initial conditions for the numerical simulations are set as follows: consider a projectile moving down the shock tube and ahead of it conditions being ambient air of p = 1 and ρ₀ = 1.4 behind of it being ambient air of p = 100.0 and ρ₀ = 140.0. At the initial frame, the left of projectile is located at x = 0.0. In this study we aim at studying the separation process, so Newton's law is not taken into consideration. The projectile always moves at the velocity Ma = 2.0. Control of discrete points is shown in Figure 9 in the calculation process. Figure 10 shows the pressure flood contours at different time.

Figure 9:

Change of calculation domain while projectile separates from shock tube.

Figure 10:

The pressure flood contours at different time.

5. Conclusion

In this paper, a new method, Overlapping Clouds of Points Method, has been proposed to dispose moving boundary problems. The method is based on Batina's meshless method and using the concept of overlapping grid method. Because the meshless method only involves the movement of discrete points and not restriction of stretch and distortion of mesh to dispose moving boundary problems, it can adjust boundary deformation of any range. When an object moves, only the attribute of partial points again and search satellite points in partial clouds of points structure are judged, leading to higher computation efficiency. Its flexible disposing method is suitable to be improved and used in flow field simulation. The method has been successfully applied to numerically simulate the unsteady flow field with moving objects. The method can be applied not only to vibration problem with small deformation but also to multibody separation problem with large deformation. The simulation results show good agreement with the experiment result and mesh results available.