A numerical simulation method for fluid motion based on the Euler-Lagrange coupled field

Abstract

The present invention discloses a numerical simulation method for fluid motion based on the Euler-Lagrange coupled field, including: obtaining a set of fixed grids and a set of virtual grids; each virtual grid moves a corresponding distance in the flow field according to its current velocity and acceleration vectors within dt time to obtain a new velocity and then reaches a new position; the moved virtual grid distributes the fluid volume and momentum it carries to four fixed grids that overlap with it in position nearby, and the virtual grid disappears automatically after distribution; according to the movement of the virtual grid and the distributed fluid volume and momentum information, the mass and momentum of the fluid in each fixed grid are updated to obtain the fluid thickness and velocity at the next calculation time step; the above process is repeated until the maximum calculation time set by the program to achieve the numerical analysis of the complete flow field dynamics. The present invention solves the problems of high time cost in the process of solving fluid mechanics by conventional CFD and the inability to be applicable to the simulation of any fluid dynamics scenario in a large scale range.

Description

Technical Field

[0001] The invention belongs to the technical field of numerical simulation methods of fluid mechanics processes, and the invention relates to a numerical simulation method of fluid motion based on Euler-Lagrangian coupling field. Background Art

[0002] Fluid mechanics is an eternal topic that cannot be avoided in the process of human beings getting along with nature. Phenomena such as turbulence have always been hot and difficult issues in the world. The Navier-Stokes (NS) equations derived from the laws of conservation of mass, momentum and energy are a set of nonlinear partial differential equations. They are universal equations that describe the motion of viscous fluids while considering nonlinear boundary conditions under the assumptions of fluid continuum and isotropic homogeneity. In addition to the potential flow theory under the assumptions of inviscidity, incompressibility and irrotation, which can reduce the accuracy and obtain analytical solutions after simplifying the NS equations, it is extremely difficult to obtain analytical solutions for the unsimplified NS equations, making it one of the seven millennium mathematical problems of the Clay Institute in the United States.

[0003] With the development of modern fluid mechanics, numerical mathematics and computer science, computational fluid dynamics (CFD) approximates the integral and differential terms in the fluid dynamics control equations as discrete algebraic equations, and obtains numerical solutions at discrete time / space points through numerical methods such as FDM, FVM, LBM and SPH. CFD numerical methods have the characteristics of low computational cost, high accuracy and realistic modeling scenarios, which have partially replaced expensive experimental schemes in the engineering field. However, the convection term with second-order nonlinearity derived from the conversion of the NS equations from Lagrangian description to Euler description makes the computational time cost of CFD numerical analysis huge, making it difficult to apply to a wide range of fluid scenarios.

[0004] Some numerical schemes that simplify the NS equations, such as Long-wave, Boussinesq, GN, RANS and other equations, can improve computational efficiency to a certain extent, but they seriously reduce the accuracy of the solution. Therefore, the problems of simple and efficient or complex and inefficient existing numerical methods seriously limit the modeling and engineering application needs of large-scale physical process scenes in fluid mechanics. Since the NS depth-averaged equation can more accurately reflect the fluid dynamics process under the depth-averaged assumption, the simplified equation is easy to solve without losing accuracy, so that the example can break through the size effect and realize the simulation of a large range of fluid scenes.

[0005] Therefore, how to adopt reasonable technical means to obtain numerical solutions of fluid mechanics in time and space, while achieving high accuracy and computational efficiency, is a technical problem that needs to be solved urgently. Summary of the invention

[0006] To solve the above problems existing in the prior art, the present invention provides a numerical simulation method for fluid motion based on an Euler-Lagrange coupled field, which solves the problems of high time cost in the conventional CFD for solving fluid mechanics and the inability to be applicable to the simulation of any fluid dynamics scenario in a large scale range. The technical problems to be solved by the present invention are realized through the following technical solutions:

[0007] A numerical simulation method for fluid motion based on an Euler-Lagrange coupled field, the fluid motion numerical simulation method comprising:

[0008] Step 1, obtain a set of fixed grids and a set of virtual grids, both the fixed grids and the virtual grids being rectangular grids composed of N*M square grids, and the fixed grids and the virtual grids coinciding in position at the initial moment of the time step;

[0009] Step 2, at the time step t+dt, obtain the coordinate position where the i-th square grid of the virtual grid moves into the fixed grid;

[0010] Step 3, obtain the new average velocity of the i-th square grid of the virtual grid after moving to the new position;

[0011] Step 4, after the position of the i-th square grid of the virtual grid is updated within dt time, update the index of the fixed grid closest to the center position of the i-th square grid of the moved virtual grid;

[0012] Step 5, based on the index of the fixed grid, use D c to normalize the vector between the center position of the i-th square grid of the moved virtual grid and the center position of the square grid of the fixed grid closest to it, where D c is the side length of the square grid;

[0013] Step 6, based on the vector, distribute the volume and / or momentum of the virtual grid with updated position to the four nearby fixed grids overlapping with it according to the distribution formula;

[0014] Step 7, update the volume and / or momentum of the fixed grid according to the distributed volume and / or momentum and the volume and / or momentum of the fixed grid itself;

[0015] Step 8, obtain the fluid thickness and velocity at the next calculation time step according to the mass and momentum of the fluid in the updated fixed grid;

[0016] Step 9, judge whether the time reaches the maximum calculation time. If not, repeat Steps 2 to 8. If so, complete the full numerical analysis of the flow field dynamics.

[0017] In an embodiment of the present invention, step 2 includes:

[0018] Step 2.1: At time step t + dt, obtain the new position of the i-th square grid of the virtual grid after movement.

[0019] Step 2.2: Based on the new position, driven by the average velocity and average acceleration, obtain the coordinate position of the i-th square grid of the virtual grid moving into the fixed grid.

[0020] In an embodiment of the present invention, the expression of the new position is:

[0021]

[0022] The expression of the coordinate position is:

[0023]

[0024] Where r i (x i , y i ) is the position of the i-th square grid of the virtual grid at the t-th moment, v i (t) is the velocity of the i-th square grid of the virtual grid at the t-th moment, a i (t) is the acceleration of the i-th square grid of the virtual grid at the t-th moment, is the new position of the i-th square grid of the virtual grid at the t + dt moment, is the x-axis coordinate value in the new position, x i is the x-axis coordinate value in the original position, v x (r i , t) is the x-axis velocity in the original position, a x (r i , t) is the x-axis acceleration in the original position, is the y-axis coordinate value in the new position, y i is the y-axis coordinate value in the original position, v y (r i , t) is the y-axis velocity in the original position, a y (r i , t) is the y-axis acceleration in the original position.

[0025] In an embodiment of the present invention, the average acceleration of the virtual grid is the sum of the gravitational acceleration and the dynamic friction acceleration.

[0026] In an embodiment of the present invention, the expression of the gravitational acceleration is:

[0027]

[0028] The expression for the dynamic friction acceleration is as follows:

[0029]

[0030] where H(r i , t) is the fluid thickness, S(r i , t) is the surface elevation, T(r i ) is the topographic boundary, μ d is the dynamic friction coefficient, and ρ w is the fluid density.

[0031] In an embodiment of the present invention, the expression for the updated index of the fixed grid is:

[0032]

[0033] The expression for the vector vector is:

[0034]

[0035] where NINT is the nearest integer function, and (n, m) is the index of the i-th square grid in the virtual network before movement.

[0036] In an embodiment of the present invention, the expression for the allocation formula is:

[0037]

[0038] where f1 represents that the volume or momentum is allocated to the nearest square grid in the fixed grid; if dx > 0, f2 represents that the volume or momentum is allocated to the square grid on the right side of the nearest square grid in the fixed grid, if dx < 0, f2 represents that the volume or momentum is allocated to the square grid on the left side of the nearest square grid in the fixed grid; if dy > 0, f3 represents that the volume or momentum is allocated to the square grid above the nearest square grid in the fixed grid, if dy < 0, f3 represents that the volume or momentum is allocated to the square grid below the nearest square grid in the fixed grid; if dx > 0 and dy > 0, f4 represents that the volume or momentum is allocated to the square grid in the upper right corner of the nearest square grid in the fixed grid, if dx < 0 and dy < 0, f4 represents that the volume or momentum is allocated to the square grid in the lower left corner of the nearest square grid in the fixed grid.

[0039] In an embodiment of the present invention, the expression for the volume of the updated fixed grid is:

[0040]

[0041] Among them, V[n, m] represents the volume of the virtual grid to be allocated, and the

[0042] and respectively represent the volumes of the fixed grids at the corresponding positions that may overlap with the virtual grid after it moves to the new position;

[0043] The expression for the momentum of the updated fixed grid is:

[0044]

[0045] Among them, M[n, m] represents the momentum of the virtual grid to be allocated, and the

[0046] and respectively represent the momenta of the fixed grids at the corresponding positions that may overlap with the virtual grid after it moves to the new position.

[0047] In an embodiment of the present invention, before the step 9, it further includes:

[0048] Performing smoothing processing on the velocity, thickness, and surface elevation.

[0049] In an embodiment of the present invention, the expression for the smoothing processing process of the momentum is:

[0050]

[0051] ΔM = fact·[v t0 (np, mp)·hc t0 (np, mp) - v t0 (n, m)·hc t0 (n, m)]

[0052] The expression for the smoothing processing process of the velocity is:

[0053] v t1 (n, m) = v t0 (n, m) + Δv

[0054] v t1 (np, mp) = v t0 (np, mp) - Δv

[0055] Δv = fact·[v t0 (np, mp) - v t0 (n, m)]

[0056] Among them, v t0 (n,m) and v t1 (n,m) respectively represent the velocities of the fluid at the center position (n,m) at time t0 before smoothing and at time t1 after smoothing, v t0 (np,mp) and v t1 (np,mp) respectively represent the velocities of the fluid at the position to be traversed (np,mp) at time t0 before smoothing and at time t1 after smoothing, hc t0 (n,m) and hc t0 (np,mp) respectively represent the fluid thicknesses at the positions (n,m) and (np,mp) at time t0 before smoothing, fact represents the correction coefficient related to the fluid thickness, h min represents the minimum limit value of the fluid thickness; max{} represents the maximum value among the specified values;

[0057] The expression for the smoothing process of the thickness is:

[0058] hc t1 (n,m) = hc t0 (n,m) + Δh

[0059] hc t1 (np,mp) = hc t0 (np,mp) - Δh

[0060] Δh = fact·[hc t0 (np,mp)·topo t0 (np,mp) - hc t0 (n,m)·topo t0 (n,m)]

[0061] Among them, when Δh > 0 and Δh > hc t0 (np,mp), Δh = hc t0 (np,mp); when Δh < 0 and Δh < -hc t0 (n,m), Δh = -hc t0 (n,m); topo t0 (n,m) and topo t0 (np,mp) respectively represent the bottom surface elevations of the fluid at the positions (n,m) and (np,mp) at time t0 before smoothing, hc t0 (n,m) and hc t0 (np,mp) respectively represent the fluid thicknesses at the positions (n,m) and (np,mp) at time t0 before smoothing, hc t1 (n,m) and hc t1(np, mp) respectively represent the fluid thickness at positions (n, m) and (np, mp) at time t1 after smoothing;

[0062] The expression for the smoothing process of the surface elevation is:

[0063]

[0064] where S t0 (n, m), S t0 (np, mp) respectively represent the surface elevation at positions (n, m) and (np, mp) at time t0 before smoothing, S t1 (n, m) represents the surface elevation at position (n, m) at time t1 after smoothing, w1 represents the elevation weight at position (n, m), and w2 represents the elevation weight at position (np, mp).

[0065] Advantages of the present invention:

[0066] (1) The most essential feature of the numerical simulation method provided by the present invention compared with the prior art methods is that: the starting point of thinking of the method of the present invention is not the control equations that macroscopically describe the hydrodynamic characteristics, but the interaction rules that strictly abide by the mass and momentum conservation laws described by Newton's second law between fluid micro-elements (meshes) behind the control equations. Therefore, the method of the present invention can more truly express the characteristics and essence of the macroscopic motion of complex fluid systems from the first principles that are more fundamental than solving the control equations. Based on the above, the method of the present invention can also more accurately reflect the hydrodynamic process under the assumption of depth-averaging, and can also save the computational time cost exponentially, improve the numerical calculation efficiency, and is applicable to simulating any hydrodynamic scenarios in a large scale range, including tsunamis, debris flows, sediment-laden floods, storm surges, glacier movements, and various two-phase flow movements, etc.

[0067] (2) The numerical simulation method provided by the present invention records the initial condition characteristics such as the volume, velocity, and acceleration of fluid micro-elements by setting up a square grid in the flow field. The method of the present invention moves the grid according to the fluid velocity and acceleration while distributing the grid mass and momentum to satisfy the mass and momentum conservation of the flow field without solving any control equations; essentially, it is to equivalently solve the linear difference equations through infinitesimal displacements, but it can avoid the computational stability problems brought by the infinitesimal motions in this process. Since the method of the present invention does not need to solve complex differential equations and avoids the calculation of individual numerical derivatives, complex boundary conditions of the flow field such as free surface and terrain boundary conditions can be adaptively satisfied during the movement of the numerical grid.

[0068] (3) The numerical simulation method provided by the present invention sets up two sets of numerical square grids: a fixed grid and a virtual grid. The fixed grid uses Eulerian to describe the kinematic information of the entire flow field, and the position of each grid remains fixed; the virtual grid uses Lagrangian to track the micro-element information of the flow field and only exists within one computational time step. After each virtual grid completes the volume and momentum distribution according to the calculation results, it disappears. The method of the present invention realizes the Euler-Lagrangian coupled field through two sets of grids, avoiding problems such as complex grid division and grid reconstruction in the traditional technology, and solving the problem of huge computational amount caused by the second-order non-linear convection term derived when the NS equation is transformed from Lagrangian description to Eulerian description. The numerical computational amount is reduced from N 2 to the order of 4N, exponentially reducing the computational time cost and greatly improving the computational efficiency.

[0069] (4) The numerical simulation method provided by the present invention makes assumptions of continuous medium, homogeneous in each item, incompressible and depth-averaged for the fluid. When calculating the fluid motion process according to the method of the present invention, only the source of the average acceleration of the micro-element grid needs to be focused on, and by expanding different external accelerations (derived from gravity, wind, air pressure, viscous resistance or boundary friction, etc.), various single-phase or multi-phase flow dynamics processes can be described; by expanding the acceleration under the mechanism of fluid-structure interaction, various fluid-structure coupling dynamics processes can be described. For compressible fluids, the fluid compressibility coefficient can be considered during the grid volume and momentum distribution process to ensure mass and momentum conservation. In summary, compared with the prior art, the method of the present invention has good scalability in both fluid physics and engineering application scenarios. The method of the present invention is relatively easy to discard the depth-averaged assumption and expand to a three-dimensional fluid numerical simulation method. The present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Description of the Drawings

[0070] Figure 1 is a schematic flow chart of a numerical simulation method for fluid motion based on the Euler-Lagrangian coupled field provided by an embodiment of the present invention;

[0071] Figure 2 is a schematic flow chart of another numerical simulation method for fluid motion based on the Euler-Lagrangian coupled field provided by an embodiment of the present invention;

[0072] Figure 3 is a schematic diagram of the specific content of the acceleration term and the corresponding flow field type in a numerical simulation method for fluid motion based on the Euler-Lagrangian coupled field provided by an embodiment of the present invention;

[0073] Figure 4 is a schematic diagram of the smoothing processing technology type in a numerical simulation method for fluid motion based on the Euler-Lagrangian coupled field provided by an embodiment of the present invention;

[0074] Figure 5 Schematic diagram of the smoothing range of a numerical simulation method for fluid motion based on the Euler-Lagrange coupled field provided by an embodiment of the present invention. Detailed implementation manners

[0075] The present invention will be further described in detail below with reference to specific embodiments, but the implementation manners of the present invention are not limited thereto.

[0076] Embodiment 1

[0077] Please refer to Figure 1 and Figure 2 , Figure 1 which are the schematic flowcharts of a numerical simulation method for fluid motion based on the Euler-Lagrange coupled field provided by an embodiment of the present invention, Figure 2 and

[0078] Step 1: Obtain a set of fixed grids and a set of virtual grids. Both the fixed grids and the virtual grids are rectangular grids composed of N*M square grids. The fixed grids and the virtual grids coincide in position at the initial moment of calculation (time step).

[0079] Specifically, the method of this embodiment simplifies the flow field to be homogeneous in all directions, incompressible, and depth-averaged, and divides the research object into an N*M rectangular grid composed of a number of square grids. The side length of each square grid is set to D c , and the center point records the coordinate information r i (t) = (x i , y i ). The fixed grids describe the kinematic information of the flow field in the Eulerian manner and remain fixed in position throughout the numerical process; the virtual grids track the microelement information of the flow field in the Lagrangian manner and can only change in position within dt time. The Euler-Lagrange coupled numerical field is realized through the fixed grids and the virtual grids. To ensure the stability of numerical calculation, the value of time dt must be a tiny time interval to ensure that the distance moved by the fluid within dt time in each grid does not exceed the range of an adjacent single grid. For example, the value range of dt is 0 < dt < (Dc / v i ). At time t, the fixed grids and the virtual grids coincide in position at the initial moment of the time step and both store the fluid thickness information H i (t) = H(r i , t) and the average velocity v i (t) = v(r i , t) and the average acceleration a i(t) = a(r i , t) and other information.

[0080] Step 2: Within the time step t + dt, obtain the coordinate position where the i-th square grid of the virtual grid moves to the fixed grid.

[0081] In a specific embodiment, Step 2 may specifically include:

[0082] Step 2.1: Within the time step t + dt, obtain the new position of the i-th square grid of the virtual grid after movement.

[0083] Specifically, at the time step t + dt, taking the i = 20-th square grid of the fixed grid as an example, assuming the index of this square grid in the fixed network is (n, m), this square grid is initialized with the fluid position r i =(x i , y i ), thickness H(r i , t), velocity v(r i , t), acceleration a(r i , t) and terrain boundary T[r i (x, y)] and other information. At this time, the virtual grid r i＝20 moves according to Newton's laws of motion (F = ma), then within the time step t + dt, the expression for the new position of the i-th square grid of the virtual grid after movement is:

[0084]

[0085] where r i (x i , y i ) is the position of the i-th square grid of the virtual grid at the t-th moment, v i (t) is the velocity of the i-th square grid of the virtual grid at the t-th moment, a i (t) is the acceleration of the i-th square grid of the virtual grid at the t-th moment, is the new position of the i-th square grid of the virtual grid at the t + dt moment.

[0086] Step 2.2: Based on the new position, driven by the average velocity and average acceleration, obtain the coordinate position where the i-th square grid of the virtual grid moves to the fixed grid.

[0087] Specifically, as can be seen from Step 2.1, each square grid of the virtual grid moves a corresponding distance within the dt time according to its current velocity and acceleration vectors in the flow field grid, and thus reaches a new position in the Euler fixed network after obtaining a new velocity. Then the expression for the coordinate position of the new position is:

[0088]

[0089] Among them, is the coordinate value of the x-axis in the new position, x i is the coordinate value of the x-axis in the original position, v x (r i , t) is the velocity of the x-axis in the original position, a x (r i , t) is the acceleration of the x-axis in the original position, is the coordinate value of the y-axis in the new position, y i is the coordinate value of the y-axis in the original position, v y (r i , t) is the velocity of the y-axis in the original position, a y (r i , t) is the acceleration of the y-axis in the original position.

[0090] In this embodiment, the average acceleration of the virtual grid is the sum of the gravitational acceleration and the dynamic friction acceleration.

[0091] The numerical process of this embodiment is equivalent to solving the flow field continuity control equation to ensure the conservation of flow field mass and momentum. The key to this numerical process is: accurately giving the average acceleration information of a single square grid within each time step. Please refer to Figure 3 , and the sources of this average acceleration include two aspects: gravitational acceleration and dynamic friction acceleration.

[0092] Regarding the gravitational acceleration of the fluid, first assume that the elevation of the free surface of the fluid within the virtual grid is S[r i (x, y)], and the terrain boundary is T[r i (x, y)], then the thickness of the fluid within the virtual grid is H[r i (x, y)] = S[r i (x, y)] - T[r i (x, y)], so the gravity and gravitational acceleration exerted on the fluid virtual grid are respectively:

[0093]

[0094] Among them, is the unit vector in the x-axis direction, is the unit vector in the y-axis direction, and M is the mass of the fluid.

[0095]

[0096] Among them, g is the gravitational acceleration.

[0097] The expression for the acceleration (i.e., dynamic friction acceleration) generated by the dynamic friction resistance of the fluid virtual grid considered in this embodiment is as follows:

[0098]

[0099] In the formula, μ d is the dynamic friction coefficient, and ρ w is the fluid density.

[0100] For the numerical analysis of fluid mechanics in complex situations, corresponding acceleration terms need to be added according to the situation. For example, considering the coupling of wind and waves, air pressure, viscous resistance, or boundary friction, etc., various single-phase or multiphase flow dynamic processes can be described; by extending the acceleration under the interaction mechanism between the fluid and the structure, various fluid-structure interaction dynamic processes can be described.

[0101] Step 3: Obtain the new average velocity of the i-th square grid of the virtual grid after moving to a new position.

[0102] Specifically, the new position of the i-th square grid of the virtual grid at the (t + dt) moment can also be expressed as:

[0103]

[0104] The expression for the new average velocity of the i-th square grid of the virtual grid after moving to a new position is:

[0105]

[0106] Step 4: After the position of the i-th square grid of the virtual grid is updated within dt time, update the index of the fixed grid closest to the center position of the i-th square grid of the moved virtual grid.

[0107] Specifically, after the position of the i-th square grid of the virtual grid is updated within dt time, the expression for the index of the fixed grid closest to the center position of the moved virtual grid can be obtained as:

[0108]

[0109] Among them, NINT is the nearest integer function, and (n, m) is the index of the i-th square grid in the virtual network before movement. Therefore, the center position coordinates of the fixed grid closest to the virtual grid should be

[0110] Step 5: Based on the index of the fixed grid, use D c to normalize the vector between the center position of the i-th square grid of the moved virtual grid and the center position of the square grid of the fixed grid closest to it. The expression for this vector is:

[0111]

[0112] Step 6: Based on the vector, distribute the volume and / or momentum of the virtual grid with updated position to the four nearby fixed grids that overlap with it according to the distribution formula. After the distribution, the virtual grid automatically disappears. The expression of the distribution formula is:

[0113]

[0114] where f1 represents that the volume or momentum is distributed to the nearest square grid in the fixed grid; if dx > 0, f2 represents that the volume or momentum is distributed to the square grid on the right of the nearest square grid in the fixed grid, if dx < 0, f2 represents that the volume or momentum is distributed to the square grid on the left of the nearest square grid in the fixed grid; if dy > 0, f3 represents that the volume or momentum is distributed to the square grid above the nearest square grid in the fixed grid, if dy < 0, f3 represents that the volume or momentum is distributed to the square grid below the nearest square grid in the fixed grid; if dx > 0 and dy > 0, f4 represents that the volume or momentum is distributed to the square grid in the upper right of the nearest square grid in the fixed grid, if dx < 0 and dy < 0, f4 represents that the volume or momentum is distributed to the square grid in the lower left of the nearest square grid in the fixed grid.

[0115] Step 7: Update the volume and / or momentum of the fixed grid according to the distributed volume and / or momentum and the volume and / or momentum of the fixed grid itself.

[0116] That is to say, the updated volume of the four nearby fixed grids should be the sum of their own volume and the distributed volume, and the updated momentum of the four nearby fixed grids should be the sum of their own momentum and the distributed momentum.

[0117] In this embodiment, the expression of the volume of the updated fixed grid is:

[0118]

[0119] where V[n,m] represents the volume of the virtual grid to be distributed, and the

[0120] and respectively represent the volumes of the fixed grids at the corresponding positions that the virtual grid may overlap with after moving to the new position.

[0121] In this embodiment, the expression of the momentum of the updated fixed grid is:

[0122]

[0123] Among them, M[n,m] represents the momentum of the virtual grid to be allocated. The

[0124] and respectively represent the momentum of the fixed grid at the corresponding position that may overlap with the virtual grid after it moves to the new position.

[0125] Step 8: Obtain the fluid thickness and velocity at the next calculation time step based on the updated mass and momentum of the fluid in the fixed grid.

[0126] Specifically, after updating the fluid mass and momentum of each square grid in the fixed grid, the fluid thickness and velocity at the next calculation time step can be calculated based on the updated mass and momentum. That is, calculate the thickness according to the mass (under the assumption of incompressibility, mass gives volume, and volume divided by the bottom area Dc*Dc of each grid gives the thickness); similarly, momentum divided by mass gives the velocity.

[0127] Step 9: Please refer to Figure 4 , and perform smoothing processing on the velocity, thickness, and surface elevation.

[0128] (1) Smoothing of the average velocity of the fixed grid fluid: As Figure 5 shown, centered on (n,m), traverse the surrounding square grid area (the traversal range is, for example, 3×3) for smoothing processing, and update the fluid velocity of all grids in the square grid area. When the smaller value of the fluid thickness at the center point (n,m) and the traversed position (np,mp) exceeds the lowest limit value, perform smoothing processing based on momentum, otherwise perform smoothing processing based on velocity.

[0129] In this embodiment, the momentum smoothing processing process is as shown in formulas (1) to (3):

[0130]

[0131] ΔM = fact·[v t0 (np,mp)·hc t0 (np,mp) - v t0 (n,m)·hc t0 (n,m)] (3)

[0132] In this embodiment, the velocity smoothing processing process is as shown in formulas (4) to (6):

[0133] v t1 (n,m) = v t0 (n,m) + Δv (4)

[0134] vt1 (np, mp) = v t0 (np, mp) - Δv (5)

[0135] Δv = fact·[v t0 (np, mp) - v t0 (n, m)] (6)

[0136] Wherein, v t0 (n, m), v t1 (n, m) respectively represent the average velocity of fluid movement at the central position (n, m) at time t0 before smoothing and at time t1 after smoothing; v t0 (np, mp), v t1 (np, mp) respectively represent the average velocity of fluid movement at the position to be traversed (np, mp) at time t0 before smoothing and at time t1 after smoothing; hc t0 (n, m), hc t0 (np, mp) respectively represent the fluid thickness at positions (n, m) and (np, mp) at time t0 before smoothing, fact represents the correction coefficient related to the fluid thickness, h min represents the minimum limit value of the fluid thickness; max{} represents the maximum value among the specified values.

[0137] (2) Grid fluid thickness smoothing: As Figure 5 shown, with (n, m) as the center, traverse the surrounding 3×3 square grid area for smoothing and update the fluid thickness of all grids within the square area. First, determine the correction amount Δh of the fluid thickness:

[0138] Δh = fact·[hc t0 (np, mp)·topo t0 (np, mp) - hc t0 (n, m)·topo t0 (n, m)]

[0139] Wherein, when Δh > 0 and Δh > hc t0 (np, mp), Δh = hc t0 (np, mp); when Δh < 0 and Δh < -hc t0 (n, m), Δh = -hc t0 (n, m); topo t0 (n, m), topo t0 (np, mp) respectively represent the bottom elevation of the fluid at positions (n, m) and (np, mp) at time t0 before smoothing.

[0140] Then, update the fluid thickness:

[0141] hc t1 (n,m) = hc t0 (n,m) + Δh

[0142] hc t1 (np,mp) = hc t0 (np,mp) - Δh

[0143] Among them, hc t0 (n,m), hc t0 (np,mp) respectively represent the fluid thickness at the positions of (n,m) and (np,mp) at time t0 before smoothing, hc t1 (n,m), hc t1 (np,mp) respectively represent the fluid thickness at the positions of (n,m) and (np,mp) at time t1 after smoothing.

[0144] (3) Mesh surface elevation smoothing: As Figure 5 shown, with (n,m) as the center, traverse the surrounding 3×3 square grid area for smoothing, and update the surface elevation at the center position of the square area. The surface elevation at the center position of the grid after smoothing is:

[0145]

[0146] Among them, S t0 (n,m), S t0 (np,mp) respectively represent the surface elevation at the positions of (n,m) and (np,mp) at time t0 before smoothing; S t1 (n,m) represents the surface elevation at the position of (n,m) at time t1 after smoothing; w1 represents the elevation weight at the position of (n,m), generally taking 1.0; w2 represents the elevation weight at the position of (np,mp), generally taking 0.2.

[0147] Step 10: Determine whether the time has reached the maximum calculation time set by the program. If not, repeat Steps 2 to 9. If so, complete the full numerical analysis of the flow field dynamics.

[0148] In summary, the above numerical process is equivalent to solving the flow field continuity control equation to ensure the conservation of mass and momentum in the flow field. The key to this numerical process is: accurately giving the average acceleration information of a single grid within each time step. The sources of this average acceleration include two aspects: gravitational acceleration and dynamic friction acceleration, etc. For the numerical analysis of fluid mechanics in complex situations, corresponding acceleration terms need to be added according to the situation. Uncertainties in the selection of calculation time steps, sudden changes in boundary conditions such as terrain, etc. will generate numerical noise, leading to instability in numerical calculations. The method of the present invention specifically proposes a technical solution for numerical smoothing.

[0149] (1) The most essential feature of the numerical simulation method provided by the present invention compared with the prior art methods is that the starting point of the method of the present invention is not the control equations that macroscopically describe the hydrodynamic characteristics, but the interaction rules that strictly abide by the mass and momentum conservation laws between fluid micro-elements (meshes) described by Newton's second law behind the control equations. Therefore, the method of the present invention can more truly express the characteristics and essence of the macroscopic motion of complex fluid systems from the first principles that are more fundamental than solving the control equations. Based on the above, the method of the present invention can also more accurately reflect the hydrodynamic process under the assumption of depth-averaging, and can exponentially save the computational time cost, improve the numerical calculation efficiency, and is applicable to simulating any hydrodynamic scenarios within a large scale, including tsunamis, debris flows, sand-carrying floods, storm surges, glacier movements, and various two-phase flow movements, etc.

[0150] (2) The numerical simulation method provided by the present invention records the initial condition characteristics such as the volume, velocity, and acceleration of fluid micro-elements by setting square meshes in the flow field. The method of the present invention distributes the mesh mass and momentum while moving the meshes according to the fluid velocity and acceleration, so as to satisfy the mass and momentum conservation of the flow field without solving any control equations; essentially, it is to equivalently solve the linear difference equation set through infinitesimal displacements, but the computational stability problem caused by the infinitesimal motion in this process can be avoided. Since the method of the present invention does not need to solve complex differential equation sets and avoids individual numerical derivative calculations, complex boundary conditions of the flow field such as free surface and terrain boundary conditions can be adaptively satisfied during the numerical mesh movement process.

[0151] (3) The numerical simulation method provided by the present invention sets two sets of numerical square meshes: fixed meshes and virtual meshes. The fixed meshes use Euler to describe the kinematic information of the entire flow field, and the position of each mesh remains fixed; the virtual meshes track the information of fluid micro-elements through Lagrangian description and only exist within one computational time step. After each virtual mesh completes the volume and momentum distribution according to the calculation results, it disappears. The method of the present invention realizes the Euler-Lagrangian coupled field through the two sets of meshes, avoids problems such as complex mesh division and mesh reconstruction required by the traditional technology, solves the problem of huge computational amount caused by the second-order nonlinear convection term derived when the NS equation is transformed from Lagrangian description to Eulerian description, reduces the numerical computational amount from N2 to the 4N order of magnitude, exponentially reduces the computational time cost, and greatly improves the computational efficiency.

[0152] (4) The numerical simulation method provided by the present invention makes assumptions of continuous medium, homogeneous in all aspects, incompressible, and depth-averaged for the fluid. According to the method of the present invention, when calculating the fluid motion process, only the source of the average acceleration of the micro-element grid needs to be focused on, and by expanding different external accelerations (derived from gravity, wind, air pressure, viscous resistance, or boundary friction, etc.), various single-phase or multi-phase flow dynamics processes can be described; by expanding the acceleration under the mechanism of fluid-structure interaction, various fluid-structure interaction dynamics processes can be described. For compressible fluids, the fluid compressibility coefficient can be considered during the grid allocation of volume and momentum to ensure mass and momentum conservation. In summary, compared with the prior art, the method of the present invention has good scalability both in terms of fluid physics and engineering application scenarios. The method of the present invention can be easily extended to a three-dimensional fluid numerical simulation method by abandoning the depth-averaged assumption.

[0153] In the description of the present invention, the terms "first" and "second" are only used for descriptive purposes and cannot be construed as indicating or implying relative importance or implicitly specifying the quantity of the indicated technical features. Thus, the features defined with "first" and "second" may explicitly or implicitly include one or more of such features. In the description of the present invention, "a plurality" means two or more, unless otherwise specifically defined.

[0154] In the description of this specification, the description with reference to terms such as "one embodiment", "some embodiments", "example", "specific example", or "some examples", etc. means that the specific features, structures, materials, or specific data points described in connection with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic representations of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials, or specific data points described can be combined in a suitable manner in any one or more embodiments or examples. In addition, those skilled in the art can combine and combine the different embodiments or examples described in this specification.

[0155] The above content is a further detailed description of the present invention in combination with specific preferred embodiments, and it cannot be determined that the specific implementation of the present invention is only limited to these descriptions. For those of ordinary skill in the technical field to which the present invention pertains, without departing from the concept of the present invention, several simple deductions or substitutions can still be made, and all should be regarded as belonging to the protection scope of the present invention.

Claims

1. A numerical simulation method for fluid motion based on the Euler-Lagrange coupled field, characterized in that The numerical simulation method of fluid motion is used to simulate the dynamic scenario of fluid motion; The numerical simulation method of fluid motion includes: Step 1: Obtain a set of fixed grids and a set of virtual grids. Both the fixed grid and the virtual grid are rectangular grids composed of N*M square grids, and the positions of the fixed grid and the virtual grid coincide at the initial moment of the time step. Among them, the fixed grid uses the Euler method to describe the kinematic information of the flow field, and the virtual grid tracks the microelement information of the flow field through the Lagrangian method. The Euler-Lagrange coupled field is realized through the fixed grid and the virtual grid; Step 2: At the time step t+dt, obtain the coordinate position where the i-th square grid of the virtual grid moves into the fixed grid; Step 3: Obtain the new average velocity of the i-th square grid of the virtual grid after moving to the new position; Step 4: After the position of the i-th square grid of the virtual grid is updated within dt time, update the index of the fixed grid closest to the center position of the i-th square grid of the moved virtual grid; Step 5. Based on the index of the fixed grid, use D c to normalize the vector between the central position of the i-th square grid of the virtual grid after movement and the central position of the square grid of the fixed grid closest to it, where D c is the side length of the square grid; Step 6: Based on the vector vector, distribute the volume and / or momentum of the virtual grid with updated position to the four nearby fixed grids that overlap with it according to the distribution formula; Step 7: Update the volume and / or momentum of the fixed grid according to the distributed volume and / or momentum and the volume and / or momentum of the fixed grid itself; Step 8: Obtain the fluid thickness and velocity at the next calculation time step according to the mass and momentum of the fluid in the updated fixed grid; Step 9: Judge whether the time reaches the maximum calculation time. If not, repeat Steps 2 to 8. If so, complete the numerical analysis of the full flow field dynamics, and obtain the fluid thickness and velocity corresponding to the maximum calculation time to realize the numerical simulation of fluid motion in the dynamic scenario; Step 2 includes: Step 2.1: At the time step t+dt, obtain the new position of the i-th square grid of the virtual grid after movement. The new position is calculated using the position of the i-th square grid of the virtual grid at the t-th moment, the velocity of the i-th square grid of the virtual grid at the t-th moment, and the acceleration of the i-th square grid of the virtual grid at the t-th moment. Step 2.2: Based on the new position, driven by the average velocity and average acceleration, obtain the coordinate position where the i-th square grid of the virtual grid moves into the fixed grid. The average acceleration of the virtual grid is the sum of the gravitational acceleration and the dynamic friction acceleration, and the dynamic friction acceleration is calculated through the fluid thickness, surface elevation, terrain boundary, dynamic friction coefficient, and fluid density.

2. The numerical simulation method of fluid motion based on the Euler-Lagrange coupled field according to claim 1, wherein The expression of the new position is: The expression of the coordinate position is: where r i (x i , y i ) is the position of the i-th square grid of the virtual grid at the t-th moment, v i (t) is the velocity of the i-th square grid of the virtual grid at the t-th moment, a i (t) is the acceleration of the i-th square grid of the virtual grid at the t-th moment, is the new position of the i-th square grid of the virtual grid at the (t + dt)-th moment, is the x-axis coordinate value in the new position, x i is the x-axis coordinate value in the original position, v x (r i , t) is the x-axis velocity in the original position, a x (r i , t) is the x-axis acceleration in the original position, is the y-axis coordinate value in the new position, y i is the y-axis coordinate value in the original position, v y (r i , t) is the y-axis velocity in the original position, a y (r i , t) is the y-axis acceleration in the original position.

3. The numerical simulation method for fluid motion based on the Euler-Lagrange coupled field according to claim 1, characterized in that The expression of the gravitational acceleration is: The expression of the dynamic friction acceleration is: Among them, is the unit vector in the x-axis direction, is the unit vector in the y-axis direction, H(r i , t) is the fluid thickness, S(r i , t) is the surface elevation, T(r i ) is the topographic boundary, μ d is the dynamic friction coefficient, ρ w is the fluid density.

4. The numerical simulation method for fluid motion based on the Euler-Lagrange coupled field according to claim 2, wherein The expression of the index of the updated fixed grid is: The expression of the vector vector is: Among them, NINT is the nearest integer function, and (n, m) is the index of the i-th square grid in the virtual network before movement.

5. The numerical simulation method for fluid motion based on the Euler-Lagrange coupled field according to claim 4, wherein The expression of the distribution formula is as follows: Among them, f1 represents that the volume or momentum is distributed to the nearest square grid in the fixed grid; if dx>0, f2 represents that the volume or momentum is distributed to the square grid on the right side of the nearest square grid in the fixed grid, if dx<0, f2 represents that the volume or momentum is distributed to the square grid on the left side of the nearest square grid in the fixed grid; if dy>0, f3 represents that the volume or momentum is distributed to the square grid above the nearest square grid in the fixed grid, if dy<0, f3 represents that the volume or momentum is distributed to the square grid below the nearest square grid in the fixed grid; if dx>0 and dy>0, f4 represents that the volume or momentum is distributed to the square grid in the upper right corner of the nearest square grid in the fixed grid, if dx<0 and dy<0, f4 represents that the volume or momentum is distributed to the square grid in the lower left corner of the nearest square grid in the fixed grid.

6. The numerical simulation method for fluid motion based on the Euler-Lagrange coupled field according to claim 5, wherein The expression of the volume of the updated fixed grid is as follows: Among them, V[n, m] represents the volume of the virtual grid to be allocated, and the and respectively represent the volumes of the fixed grids at the corresponding positions that may overlap with the virtual grid after it moves to the new position; The expression of the momentum of the updated fixed grid is as follows: Among them, M[n,m] represents the momentum of the virtual grid to be allocated, and the and respectively represent the momenta of the fixed grids at the corresponding positions that may overlap with the virtual grid after it moves to the new position.

7. The numerical simulation method for fluid motion based on the Euler-Lagrange coupled field according to claim 1, characterized in that Before the step 9, it further includes: Smoothing the velocity, thickness and surface elevation.

8. The numerical simulation method for fluid motion based on the Euler-Lagrange coupled field according to claim 7, characterized in that The expression of the smoothing process of the momentum is as follows: ΔM = fact·[v t0 (np,mp)·hc t0 (np,mp)-v t0 (n,m)·hc t0 (n,m)] The expression of the smoothing process of the velocity is as follows: v t1 (n,m) = v t0 (n,m) + Δv v t1 (np, mp) = v t0 (np, mp) - Δv Δv = fact·[v t0 (np,mp)-v t0 (n,m)] Among them, v t0 (n, m), v t1 (n, m) respectively represent the velocities of the fluid at the central position (n, m) at time t0 before smoothing and at time t1 after smoothing, v t0 (np, mp), v t1 (np, mp) respectively represent the velocities of the fluid at the position to be traversed (np, mp) at time t0 before smoothing and at time t1 after smoothing, hc t0 (n, m), hc t0 (np, mp) respectively represent the fluid thicknesses at the positions (n, m) and (np, mp) at time t0 before smoothing, fact represents the correction coefficient related to the fluid thickness, h min represents the minimum limit value of the fluid thickness; max{} represents the maximum value among the specified values; The expression of the smoothing process of the thickness is as follows: hc t1 (n,m) = hc t0 (n,m) + Δh hc t1 (np, mp) = hc t0 (np, mp) - Δh Δh = fact·[hc t0 (np, mp)·topo t0 (np, mp)-hc t0 (n, m)·topo t0 (n, m)] where, when Δh > 0 and Δh > hc t0 (np, mp), Δh = hc t0 (np, mp); when Δh < 0 and Δh < -hc t0 (n, m), Δh = -hc t0 (n, m); topo t0 (n, m), topo t0 (np, mp) respectively represent the bottom elevation of the fluid at the positions (n, m) and (np, mp) at time t0 before smoothing, hc t0 (n, m), hc t0 (np, mp) respectively represent the fluid thickness at the positions (n, m) and (np, mp) at time t0 before smoothing, hc t1 (n, m), hc t1 (np, mp) respectively represent the fluid thickness at the positions (n, m) and (np, mp) at time t1 after smoothing; The expression of the smoothing process of the surface elevation is as follows: Among them, S t0 (n,m), S t0 (np,mp) respectively represent the surface elevations at the positions of (n,m) and (np,mp) at time t0 before smoothing, S t1 (n,m) represents the surface elevation at the position of (n,m) at time t1 after smoothing, w1 represents the elevation weight at the position of (n,m), and w2 represents the elevation weight at the position of (np,mp).

Images