On the Application of a Mobile Grid Technology to Computational Fluid Dynamics

: On the basis of harmonic mapping theory, a mobile grid technology is applied to computational fluid dynamics (CFD). Starting from the observation that standard fixed-grid techniques often fail in addressing problems with large deformations, we elaborate a new algorithm relying on the software COMSOL Multiphysics 5.3a to solve the coupling of the mobile grid equation and the governing differential equations for fluid flow. The motion of water in a water tank when the tank waggles is simulated. We demonstrate that this technology can be implemented without a significant increase in the computational cost with respect to standard numerical methods.


Introduction
The phenomenon of fluid flow exists in the nature and various fields of engineering extensively, and all these processes are supported by the fundamental laws of physics such as conservation of mass, conservation of momentum, and conservation of energy, and namely, it needs to meet a certain governing differential equation. Computational fluid dynamics is the analysis and calculation conducted on the system of isophase physical phenomena including fluid flow and heat conduction by computer numerical calculation and image display. The basic thought of CFD can be summarized as replacing the physical fields (such as pressure field and velocity field) on the original time domain and spatial domain by a series of collection of variable values on a finite number of discrete points, and through a certain principle and method, the algebraic equation set about the relationship among the field variables on these discrete points is established, and then, the approximate value of the field variables gained by the algebraic equation set is solved [Zhou (1995); Huang and Russell (2011); Huang and Kamenski (2017)], and namely, the primary problem of CFD is to discretize these differential equations controlling fluid flow, dividing the solution area into the grid model which is convenient for solving. The conventional CFD discrete method includes finite difference method, finite element method and finite volume method, and among which, the finite volume method is applied most widely [Versteeg and Malalasekera (1995)]. They all have the methods for dividing the grid respectively, and it is collectively called as grid generation technique. The conventional discrete method is to divide the solution domain into uniform grid, while in the actual application, the problems appear frequently, such as grid large deformation, great changes happened by some physical quantities in a very small area, including fluid-solid coupling boundary, laminar flow, shear layer in turbulence, as well as explosion, and for these problems, the reasonable accuracy is hard to reach if adopting the simple uniform grid solution, and therefore, it is extremely important to generate the computational grid which matches with the problems.

Governing equations of fluid dynamics and its solutions 2.1 Governing differential equation of fluid dynamics
Any flow problem must meet the law of conservation of mass. This law can be described as that the increase of mass of fluid element per unit time equals to the net mass flowing into the element at the same time. According to this law, the mass conservation equation can be gained: In the equation, ρ is fluid density; u , v and w are the components of velocity vector in three directions of x, y and z; t is time.
Any fluid shall meet the law of conservation of momentum. This law can be described as that the rate of change of fluid momentum with respect to time in element equals to the sum of external forces acting on the element. According to this law, the famous Navier-Stokes equation can be gained: In the equation, u is velocity vector; I is unit matrix; τ is viscous stress tensor of fluid; F is the volume force of fluid, while ∇ is Hamilton operator.
Any fluid shall meet the law of conservation of energy. This law can be described as that the increase rate of element energy equals to the sum of net heat absorption of the element and external acting on the element. According to this law, the energy conservation equation can be gained: In the equation, T is the absolute temperature of fluid; cp is specific heat capacity of fluid; k is the heat transfer coefficient of fluid; ST is the internal heat source of fluid and the parts of fluid mechanical energy transferring into the heat energy because of action of clay, and sometimes, it is called as viscous dissipation term. The Eqs.
(1), (2) and (3) are the differential equations controlling fluid flow, and essentially, the process of solving computational fluid dynamics is to solve the above partial differential equation set.

Headings
No matter for the fluid flow problem or the fluid heat transfer problem, and no matter for the steady state problem or the transient problem, the solution process of CFD is as follows: a. Establish governing differential equation and confirm initial boundary value condition; b. Divide computational grid and generate the calculation node; c. Establish discrete equations, and substitute into the initial boundary value condition; d. Give solution control parameters and solve the discrete equation; e. Judge whether the solution is convergence. No, go to b; yes, go to f; f. Show the output calculation results. From the above process, it can be seen that essentially, CFD is the partial differential equation set of solution Eqs.
(1)- (3), while the emphasis and difficulty is the subdivision of solution domain grid. Under a certain situation, the traditional grid technology can well simulate the fluid flow issue (such as the aerodynamic characteristics of objects), however, for some specific issues (such as the large deformation of grid), no matter how highly accurate the numerical method is, if the grid condition does not match with the actual problem, the gained numerical solution will completely distorted. That is to say, the grid technology plays the same role for the accuracy of the numerical solution of partial differential equation as what the numerical method does.

Advantages of mobile grid technology and its basic thoughts
As one of the adaptive grid, mobile grid technology was proposed by Liao et al. [Liao and Anderson (1992)] in 1992, and later, it was widely applied by some scholars into the field of partial differential equations, achieving a certain effect [Song and Quan (2004); Li, Liu and Ma (2004)]. In the calculation of fluid of large deformation of grid, the traditional Lagrange method and Euler method have deficiencies: although Lagrange method possesses the advantages of being capable of accurately distinguishing the interface of matter, the grid intersecting caused by large deformation will cause the calculation interrupting; in the Euler method, the geometrical shape of grid is good, but the problems such as interface tracking at boundaries, hybrid grid and reconstruction are hard to process. Mobile grid method is a method for solving the evolution equation, and it timely adjusts the density and shape of grid according to changes of physical solution.
Mobile grid based on Arbitrary Lagrange-Euler method (for short: ALE) can both accurately distinguish material interface and keep geometrical features of grid, and therefore, when solving the large deformation of grid, it has the incomparable advantages. The features of ALE are to determine the characteristics of grid with the characteristics of solution in the solution area, and to constantly update the grid in the solution process, thus to make it automatically match with the solution. For realizing this process, a kind of mobile grid method based on harmonic mapping needs to be applied [Li, Tang and Zhang (2001); Li, Tang and Zhang (2002)], and the basic thought of this method is to totally separate the grid deformation with the solution of governing differential equation, thus, it only needs to construct a control function which suits for our problems. In this process, it needs to introduce a logic area, the motion of grid is realized by changing the physical area to logical area, and the physical area is the actual fluid area which needs to be solved.

Harmonic mapping
The most critical and difficult point of mobile grid method is to find grid mapping or grid transformation which meets a certain condition. Harmonic mapping is put forward by Fuller in the 1940s, and later, it was widely applied into the field of mathematical physics and adaptive grid field [Gowrisankar and Natesan (2013)] and the definition on harmonic mapping is as follows: For two n dimension calculation fields DP and DL, a smooth mapping is given: : p L D D ξ  → , and its energy density is: The energy of mapping ξ is the integration of energy density: In the equation, Dp and DL are the physical area and logical area respectively of the calculation area. From the above derivation, it can be seen that the essence of harmonic mapping is the promotion of harmonic function and minimal manifolds [Li (2010)]. Hamilton, Schoen, et. al established relatively complete theorem of existence and uniqueness of harmonic mapping [Hamilton (1975); Schoen and Yau (1978)], and presented the corresponding conditions. The conclusion of existence uniqueness is beneficial for avoiding the calculation interrupting caused by the interlacing of grids.

Control function
Control function is a given positive definite function matrix on solution domain, and it is the key factor for generating mobile grid. Control function can be used to control the quality of the grid in variational form and couple the grid with the physical solution of the solution domain, and at that time, control function can be used to measure the physical quantity on the physical region, while the metric matrix on the physical region can be taken as the inverse element of control function. The general principle for the construction of control function is that the gradient of numerical solution is used for judging which parts' solution in solution domain changes fast, thus the grids can be centralized on those parts. The general form of the control function constructed by gradient is: In the equation,α , β and γ are positive parameters; I is unit matrix; ∇ is Hamilton operator.
It also constructs control function by connecting control function with error estimation [Chen and Liu (2016) Their principles are the same with the control function about gradient.

Movement of grid
For the specific realization method of changes from physical region Dp of the solution region to DL of logic region, it is to confirm the mapping situation of fluid on boundaries first after having selected logic region DL and given a uniform partition T (suppose its node as x i ) on flow field DP, then the mapping on this boundary is used for the limitation showed that [Cao, Huang and Russell (2001)] for the physical region of the same solution domain, selecting different logic region will not generate significant influences on the calculation results.
If it wants to complete the movement of grid, it still needs to be realized through the moving vector fields. It assumes the corresponding physical grid of t moment is t T (its knot is i t X ), and through calculating control function M, it solves elliptic equation set on grid t T .
And namely, it can gain the image i t ξ of the node i t X of flow field grid under the harmonic mapping, and furthermore, its difference with the logic grid * i i i t δξ ξ ξ = − is gained. In virtue of the changed first derivative from logic grid to flow field grid, it is interpolated as a shard vector field on flow field grid, thus, the movement direction ,* i X δ of each node on the flow field is gained, and therefore, the node of flow field grid is updated as: In the equation, η is the step length of the grid movement. Therefore, the grid of the solution domain reaches the goal of updating. The procedure of applying mobile grid method into CFD numerical simulation is as follows: Solving control differential equations at t=t n Solving equation (11) to obtain ξ *

.1 The model background
During the process of water tank tilting and shaking, the water in the water tank rocks back and forth under the action of gravitational vector, because there are no constraints on liquid level, its deformation is huge and irregular, it is almost unable to solve with the conventional uniformed fixed grid, while this kind of problem can be well solved if we apply mobile grid technology. In this case, computational fluid dynamics module (CFD module) and mobile grid module (ALE module) in multiple physical field coupling software COMSOL Multiphysics 5.3a are applied to simulate the movement rule of water in the water tank. As the first global truly COMSOL Multiphysics, it adopts the mode of modeling based on partial differential equation, and therefore, it has its unique advantages when handling with the problems of mobile grid and governing equation coupling. The model parameters are: the rectangular tank, its height is enough, there is a certain amount of water inside, the solution domain length is 1m, width is 0.4 m, the density of water is 1000 kg/m 3 , the viscosity is 1.01×10 3 pa * s, and the maximum angle of inclination of tank swing is 4°. The fluid is Newtonian fluid, the model adopts one-way laminar flow model, and it does not give consideration to heat conduction process. According to grid's own powerful grid division function, it is freely divided into triangular gird, and the initial grid of its solution domain is shown in the following figure:

Boundary condition of the CFD and mobile grid
After the solution domain of the control equation is determined, it needs to set up the boundary condition of governing differential equation and mobile grid. The confirmation of fluid governing differential equation boundary condition: as shown in the Fig. 2, the left boundary, right boundary and lower boundary of rectangle are the sliding wall liquid level's boundary condition; the upper boundary is the free liquid level, and for opening the boundary, the normal stress is zero. The definition of mobile grid's boundary condition: as shown in Fig. 2, for the left boundary and right boundary of the rectangular region, the specified grid displacement is zero in horizontal direction, and in vertical direction, it is free deformation; for the lower boundary, the specified grid is zero on both horizontal direction and vertical direction; for the free liquid level, it adopts the boundary coordinate system, and the speed at the specified grid normal is "u*nx+v*ny", and namely, the horizontal speed u and the vertical speed v of the grid are decomposed to normal, and the tangential velocity is not restricted. Through the above setting, it is actually couple the mobile grid with the solution region, when conducting the transient analysis, the grid will generate a deformation after every time substep passes, and after updating, the grid model will become the initial condition of the next time substep, and time and again, the deformation condition of grid within a certain period can be gained.

Simulation results and analysis
In this case, the transient solver provided by software is adopted for solving, the simulation time is 5 s, and the time step is 0.1 s. Among which, the velocity vector and grid deformation figure of the time t being 1.1 s, and 1.7 s respectively is shown in the following figure: (a) t=1.1 s (b) t=1.7 s  Fig. 3, it can be seen that when the water tank waggles, the liquid in the tank also rocks back and forth under the action of gravitational vector. In this case, the COMSOL Multiphysics is applied for coupling mobile grid with the governing differential equation of fluid flow, and the deformation condition of grid is used for truly simulating the movement condition of water in water tank. For the movement of grid, it is just the adjustment of grid node position, while the node number of the data structure of grid node will not change, and therefore, the memory overhead in the computer can be saved largely, the calculation time is shorten, and the debugging of the program also becomes much simpler. Through the mobile grid technology，we can also track the movement condition of the free liquid level, and as shown in Fig. 4, the relation of elevation of the left-end liquid level and the elevation right-end liquid level of rectangle changing as time is shown as follows:  Fig. 4, it can be seen that when the elevation of the left-end free liquid level reaches the wave crest, the elevation of the right-end free liquid level reaches the trough, and vice versa; the elevation of the left-end liquid level is exactly one phase different from the elevation of the right-end liquid level, which conforms to the actual situation. Through this case, it explains that it is feasible and reasonable to apply the mobile grid technology into computational fluid dynamics.

Conclusion
CFD issue can ultimately summarized as the issue of solving partial differential equations, and therefore, we can introduce a powerful tool for solving the partial differential equation: mobile grid technology. mobile grid has the advantages that the conventional fixed grid cannot compare with, and applying it into the numerical simulation field of CFD, it can well solve the large deformation problem of the flow region widely existed in the nature, and moreover, in the simulation process, it does not change the node and grid number, and there is no extra expense for the computer resources. Therefore, mobile grid technology has the wide application prospect in the research field of numerical simulation of CFD.