A linear MARS method for interface tracking in three-dimensional multiphase flow simulation

Abstract

This paper discloses a linear MARS method for interface tracking in three-dimensional multiphase flow simulation. Given initial liquid tracer information and the flow field in which it resides, this method utilizes Yin set representation and time integration to obtain the liquid surface tracer information at a target moment during the flow process under the influence of the flow field, ultimately determining the interface morphology. This method achieves consistent fourth-order accuracy in both time and space, maximizes the geometric topological properties of the flow phase, and is economical, fast, and highly reusable.

Description

Technical Field

[0001] The present invention belongs to the technical field of liquid sloshing simulation, and in particular relates to a linear MARS (an analytic framework of interface tracking via mapping and adjusting regular semi-algebraic sets) method for the interface tracking problem in three-dimensional multiphase flow simulation. Background Art

[0002] Aircraft airworthiness certification requires a reasonable assessment of the sloshing response of wing fuel tanks during flight. Liquid sloshing simulation, based on mechanical theory and mathematical background, simulates the sloshing behavior of liquid in a fuel tank. This method ensures that the wing tank structure meets dynamic and control requirements when the aircraft is put into service, and prevents joints from breaking or loosening during sloshing, which could cause accidents. Sloshing, short for tank sloshing, generally refers to the phenomenon in which liquid sloshing occurs in a vehicle with a liquid tank (such as a ship, vehicle, or aircraft) under the influence of environmental loads (such as acceleration and external forces). Liquid sloshing is highly random and fluid, and the impact loads it generates can have a significant destructive effect on bulkheads, causing structural deformation at the least and even damage at the worst. Furthermore, liquid sloshing can significantly affect the gas-liquid interface, thereby impacting the operating environment of the tank intake. From a mathematical perspective, the sloshing problem with a free surface can be reduced to solving the Navier-Stokes equations (or Euler equations) for liquid motion. These equations are generally nonlinear, and the boundary conditions of the free surface often manifest as nonlinear equations as well, with the liquid surface position varying with time. Liquid sloshing is a complex dynamic phenomenon that depends on parameters such as the container's geometry, liquid depth, excitation, and wall properties. Solving the general sloshing equation for liquids is quite difficult. Actual industrial products (such as aircraft wing tanks, liquid fuel tanks, and liquid-carrying vehicles) are not suitable for purely theoretical research and analysis due to their complex three-dimensional models and variable operating conditions. Therefore, computer-aided techniques are often used to simulate and analyze liquid sloshing on real products.

[0003] MARS is an analytic framework for interface tracking via mapping and adjusting regular semialgebraic sets. The linear MARS interface tracking method was developed as a numerical computational method to better track interfaces in multiphase flow simulations. It provides high-precision numerical simulations of the interaction of liquids with different phases during tank sloshing. The goal is to use numerical simulations to adjust and optimize product designs, reduce the impact of sloshing, and effectively control the flow behavior of liquids within the container. Figure 1(a) shows an example of sloshing within an aircraft wing tank. When the tank contains relatively little oil, sloshing is likely to occur, impacting the tank walls. Using specialized physical structures can effectively mitigate this sloshing, as shown in Figure 1(b), but this also increases the complexity of the interface. This requires high-precision simulations of the multiphase flow during the numerical simulation phase to estimate the potential flow behavior.

[0004] When capturing liquid flow behavior, a very small mesh size is often required at the gas-liquid interface. However, since liquid slosh simulation is transient, the gas-liquid interface is constantly changing. This necessitates the use of a three-dimensional interface tracking method to accurately capture the gas-liquid interface at any given moment. In traditional liquid slosh simulation software (such as FLUENT), the VOF method [Hirt and Nichols, 1981 J. Comput. Phys.] and the LS method (level set method) are commonly used as interface tracking solvers.

[0005] In 1981, Hirts and Nichols first proposed the VOF method in a paper that numerically simulated the interface tracking problem of dam break and surge. Since its introduction, the basic idea of ​​the VOF method has been widely recognized, especially for its groundbreaking approach to the reconstruction of the transport interface. VOF defines a variable labeling function for the fluid phase M being tracked, and then for each control volume C i,j The VOF method can be further divided into two steps to define the volume ratio at each time step: 1. The reconstruction step determines a tangent or polynomial curve in the boundary cell based on the volume ratio and velocity normal in the boundary control volume. This approximates the location of the explicit boundary of the flow phase and uses the resulting boundary as the boundary condition. 2. The convection step recalculates the change in volume ratio of each control volume within that time step based on the boundary conditions obtained in the reconstruction step and the given flow field conditions. VOF maintains volume conservation, facilitating the calculation of complex interface changes such as wave rollover, breakup, and merging, and has been used in many practical simulations.

[0006] The LS method, also known as the isosurface function method, is a popular method for numerical simulation of multiphase flow interface tracking and has been applied in many fields. The level set method [Osher and Sethian 1988 J. Comput. Phys.] was originally proposed in 1988 by American computational mathematician Stanley Osher and others. The core idea of ​​the level set method is to set an isosurface function φ(x, t) in the flow field that satisfies the governing equation:

[0007]

[0008] The isosurface function The value at the interface of the flow phases is 0; at any time, by solving the control equation to obtain its zero isosurface, the interface of the flow phases can be located. The level set method does not need to explicitly represent the interface. This feature improves its ability to track complex interfaces and can more easily handle situations where topological changes occur. In order to optimize the LS algorithm, Osher et al. conducted in-depth research on the level set reinitialization method. Considering that the boundary only changes position in a certain local area, the control equation can be solved only in this smaller area. This method is called the fast local LS method. The fast local level set method greatly speeds up the computational efficiency of interface tracking, while also reducing the storage capacity of the calculation and improving its practicality.

[0009] However, the above existing technologies also have certain limitations. The limitations of the VOF method and the level set method are mainly reflected in the following aspects:

[0010] 1. The calculation accuracy is not high enough. When there is a sharp point on the motion interface (C 1 When there are three or more phases in the multiphase flow, such cusps are bound to exist, which means that the accuracy of the level set method and the VOF method will further decrease with the appearance of cusps.

[0011] 2. Inability to handle topological changes of the fluid. It is difficult for the level set method and the VOF method to accurately reflect the physical properties of the fluid phase into the interface tracking algorithm [Zhang and Fogelson 2014SIAM J.Sci.Comput.]; from a physical point of view, during the interface tracking process, when the fluid phase undergoes topological changes, the specific details of the change should be considered based on the physical properties of the fluid phase, rather than simply "automating" to handle it uniformly. In addition, when the velocity field is sufficiently continuous, the corresponding stream function should be homeomorphic, while the level set method and the VOF method may change the topology of the fluid phase during the tracking process, which is contrary to the fact that the stream function is a homeomorphic mapping. Finally, the existing methods may produce large errors when dealing with topological changes of the fluid phase. For example, when the VOF method is processing an interface merge with a scale smaller than the side length h of the control body, each control body that produces a topological change will cause O(h D ) error (D is the dimension of the computational domain), because the number of control volumes on the boundary is O(1 / h D-1 ) and the time step is O(1 / h). Therefore, when the topology of the multiphase flow changes too frequently, the accumulated errors in time and space may cause the interface tracking results of the VOF method to not converge. Summary of the Invention

[0012] In view of the above, the present invention provides a linear MARS method for the interface tracking problem in three-dimensional multiphase flow simulation, which has consistent fourth-order accuracy in time and space and maximizes the geometric topological characteristics of the flow phase, and is economical, fast and highly reusable.

[0013] A linear MARS method for interface tracking in three-dimensional multiphase flow simulation is used to determine the liquid interface morphology at a target time during the flow of a tank liquid under the action of a flow field. The linear MARS method comprises the following steps:

[0014] (1) Obtain the shape of the liquid at the initial moment and use the triangulation algorithm to record the liquid surface position information at the initial moment t0, that is, the position information of each tracking point on the liquid interface;

[0015] (2) Divide the entire time period from the initial time t0 to the target time into multiple small segments. For the first small segment from t0 to t1, calculate the liquid level position information at time t1;

[0016] (3) The liquid surface position information at time t1 is used as the liquid surface position information at the start time of the next small segment, thereby calculating the liquid surface position information at the end time of the next small segment, and so on to obtain the liquid surface position information at the target time, and finally determine the shape of the liquid interface.

[0017] Furthermore, the linear MARS method uses Yin Ji as a mathematical model to describe the area enclosed by a closed interface in three-dimensional space, and uses a linear triangulation surface to approximate the tank liquid interface, and tracks the dynamic changes of the tank liquid level over time by calculating the dynamic position of the interface.

[0018] Furthermore, the triangulation algorithm in step (1) needs to select points that are approximately evenly distributed on the liquid interface as tracking points, and provide a triangulation of the tracking points, so as to use the tracking of the tracking points on these triangulation surfaces to replace the tracking simulation of the entire liquid interface.

[0019] Furthermore, the specific implementation method of step (2) is as follows: first, the position information of each tracing point on the liquid interface at time t1 is calculated using a time integration algorithm based on the flow field information and the liquid surface position information at time t0, and the tracing point is the node in the triangular mesh; then, the side length and area of ​​each triangle in the triangular mesh at this moment are judged to increase or decrease the nodes in the mesh and update the mesh, and finally the liquid surface position information at time t1 is determined.

[0020] Furthermore, each triangle in the triangular mesh is traversed. If the side length of a triangle is greater than a threshold, the midpoint of the side is added as a tracer point, and the original triangle is then split into two triangles using the tracer point, and the entire triangular mesh is updated.

[0021] Furthermore, each triangle in the triangular mesh is traversed. If the area of ​​a triangle is smaller than the threshold and the cosine value of the maximum angle of the triangle is smaller than the threshold, no adjustment is made.

[0022] If the area of ​​a triangle is less than the threshold and the cosine value of the maximum angle of the triangle is greater than or equal to the threshold, further judgment is made: if the triangle has at least two sides less than the threshold, the triangle is shrunk into a tracer point and the entire triangle mesh is updated; if only one side of the triangle is less than the threshold, the side is shrunk into a tracer point and the entire triangle mesh is updated; if each side of the triangle is greater than the threshold, the triangle is tried to be optimized by side swapping and the entire triangle mesh is updated; if all three sides cannot be optimized by side swapping, the shortest side is shrunk into a tracer point and the entire triangle mesh is updated.

[0023] Furthermore, the shrinking of the triangle into a tracing point is to delete the three vertices of the triangle and add the centroid of the triangle as the tracing point.

[0024] Furthermore, shrinking the edge into a tracer point means deleting an endpoint of the edge whose angle of maximum normal vectors corresponding to the plurality of triangular planes is relatively small.

[0025] Furthermore, the side-exchange optimization is that for any side of a triangle, if two triangles with the side as a common side are combined into a convex quadrilateral, the diagonals of the convex quadrilateral are exchanged.

[0026] Based on the above technical solution, the present invention has the following beneficial technical effects:

[0027] 1. Second- or third-order accuracy. The linear MARS method of the present invention uses linear interpolation to add tracking points and a polyhedron composed of triangular planes to approximate the Yin set in three-dimensional space, thereby achieving an accuracy of 2α for the three-dimensional linear MARS method. However, in practical applications in three-dimensional space, a large value of α results in low algorithm efficiency. To ensure a balance with flow field calculations, α is often set to 1.5. Third-order accuracy ensures a balance in computational complexity with flow field simulation, ensuring that the computational complexity of interface tracking does not exceed the complexity of flow field calculations.

[0028] 2. Topological structure guaranteed. In three dimensions, during the adjustment step of linear MARS, the present invention limits the adjustment operation by analyzing the local flatness, thereby further preserving the interface information of the flow phase. In addition, when the stream function is a homeomorphic mapping, linear MARS can maintain its topological structure even if the flow phase deformation is large.

[0029] 3. Flexibility. The linear MARS of the present invention allows users to define the relationship between the interface scale and the control volume (i.e., the internal scale of the fluid phase) after weighing the calculation accuracy and efficiency, providing users with a more flexible and free interface. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1(a) shows an example of liquid sloshing in a wing tank.

[0031] Figure 1(b) is an example diagram of slowing down liquid sloshing under the baffle structure.

[0032] Figure 2 Schematic diagram of the linear MARS method of the present invention.

[0033] Figure 3 Schematic diagram of the tracer points on the forward tracking boundary.

[0034] Figure 4 Schematic diagram of adding tracer points to make the length of each triangle smaller than the upper bound.

[0035] Figure 5 To delete three edges smaller than r tiny h L Schematic diagram of a triangle.

[0036] Figure 6 To delete an edge smaller than r tiny h L Schematic diagram of a triangle.

[0037] Figure 7 Schematic diagram of trying to optimize the triangle using the side-swapping method for a convex quadrilateral.

[0038] Figure 8 Schematic diagram of a standard triangle with three sides deleted. DETAILED DESCRIPTION

[0039] In order to describe the present invention more specifically, the technical solution of the present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0040] In order to numerically simulate the liquid sloshing process in the fuel tank during vehicle movement, instead of physical tests (for example, replacing full-scale tests in aircraft airworthiness certification), thereby achieving the goal of saving industrial costs and improving industrial efficiency, the present invention proposes an interface tracking method that can simulate the changes of liquid interfaces with time-space consistent fourth-order accuracy, maximizes the geometric topological characteristics of the flow phase, and is economical, fast, and highly reusable; Figure 2 As shown, the main steps of the linear MARS method are as follows:

[0041] (1) Extract the tracer points on the liquid interface, obtain the position information of these tracer points at the initial moment, and use the triangulation surface to fit the approximate liquid interface surface.

[0042] We use Yin Ji [Zhang and Li 2020Math.Comput.] as a mathematical model to describe the area enclosed by a closed crystal interface in a three-dimensional plane, and use a linear triangulation surface to approximate the tank liquid interface. By calculating the dynamic position of this interface, we track the dynamic changes of the tank liquid level over time.

[0043] In the actual calculation process, we need to select approximately uniformly distributed points on the interface as tracing points and give a triangulation of the tracing points, so as to replace the tracking simulation of the entire interface with the tracking of the tracing points on these triangulation surfaces.

[0044] (2) Figure 3 As shown, the tracer points after flow mapping at the next moment are calculated, that is, the ordinary differential equation is solved by the velocity field u:

[0045]

[0046] (3) Increase the number of tracer points so that the side length of all triangles in the triangulation is no greater than the upper bound h L If the side length is greater than h L , trace back the positions of the two end tracer points of this side at the previous moment q1q2, and map the midpoint q0 of these two points to the next moment, such as Figure 4As shown, a new tracer point is obtained and added to the tracer point set {p i Repeat this process until all side lengths are less than the upper bound. The specific implementation process is as follows:

[0047] 3.1 In Upper positioning The vertex preimage q i ,q j ,q k , will be greater than the upper bound h L The sides of are bisected by their midpoints and these midpoints are marked as new tracer points;

[0048] 3.2 Add new edges, use q i ,q j ,q k Construct a new local triangulation with the new tracer points;

[0049] 3.3 Replace with a new local triangulation

[0050] 3.4 For these new tracer points, solve the ordinary differential equation of the velocity field u to obtain the new tracer points at t n +k moment position and add it to the vertex sequence {p j}middle.

[0051] (4) In the new triangulation obtained in the previous step, traverse each triangle to process triangles with too small area and too small side length to make the mesh approximately uniform.

[0052] 4.1 Delete three edges that are smaller than the lower bound r tiny h L and check the resulting change in side length, such as Figure 5 As shown, if it is greater than the upper bound, add a point to this edge according to the method in the previous step and add the newly added triangle to the unchecked triangle stack.

[0053] 4.2 Delete area smaller than And only one side is smaller than The triangle is changed to the midpoint of the side and the length of the side is checked. Figure 6 As shown, if it is greater than the upper bound, add a point to this edge according to the method in the previous step and add the newly added triangle to the unchecked triangle stack, r tiny and r a All are given coefficients.

[0054] 4.3 For a triangle with three standard sides, try to optimize the triangle by changing the sides, such as Figure 7 If that doesn't work, delete the triangle and check the resulting change in side length, as shown in Figure 8As shown; if it is greater than the upper bound, add a point to this edge according to the method in the previous step and add the newly added triangle to the unchecked triangle stack.

[0055] (5) Get a new triangulation information to approximate t n+1 The above steps are repeated until the termination time, and the final triangulation is used to approximate the liquid interface.

[0056] The above description of the embodiments is intended to facilitate understanding and application of the present invention by those skilled in the art. It is apparent that those skilled in the art can readily make various modifications to the above embodiments and apply the general principles described herein to other embodiments without requiring creative effort. Therefore, the present invention is not limited to the above embodiments. Any improvements or modifications made by those skilled in the art based on the disclosure of the present invention should fall within the scope of protection of the present invention.

Claims

1. A linear MARS method for the interface tracking problem in three-dimensional multiphase flow simulation is used to determine the liquid interface morphology at a target moment during the flow process of the tank liquid under the action of the flow field. It is characterized by: The linear MARS method uses Yin Ji as a mathematical model to describe the area enclosed by closed interfaces in three-dimensional space, and uses a linear triangulation surface to approximate the tank liquid interface. The dynamic position of the tank liquid interface is calculated to track the dynamic changes of the tank liquid level over time. Specifically, the method includes the following steps: (1) Obtain the shape of the liquid at the initial moment and use the triangulation algorithm to record the liquid surface position information at the initial moment t0, that is, the position information of each tracking point on the liquid interface; The triangulation algorithm needs to select points that are approximately evenly distributed on the liquid interface as tracking points, and give a triangulation of the tracking points, so as to use the tracking of the tracking points on these triangulation surfaces to replace the tracking simulation of the entire liquid interface; (2) Divide the entire time period from the initial time t0 to the target time into multiple small segments. For the first small segment from t0 to t1, calculate the liquid surface position information at time t1. The specific implementation method is as follows: first, calculate the position information of each tracer point on the liquid interface at time t1 using a time integration algorithm based on the flow field information and the liquid surface position information at time t0. The tracer point is the node in the triangular mesh; then, the side length and area of ​​each triangle in the triangular mesh at this moment are judged to increase or decrease the nodes in the mesh and update the mesh, and finally determine the liquid surface position information at time t1; (3) The liquid surface position information at time t1 is used as the liquid surface position information at the start time of the next small segment, thereby calculating the liquid surface position information at the end time of the next small segment, and so on to obtain the liquid surface position information at the target time, and finally determine the shape of the liquid interface.

2. The linear MARS method according to claim 1, characterized in that: Traverse each triangle in the triangular mesh. If the side length of a triangle is greater than the threshold, add the midpoint of the side as a tracer point, and then use the tracer point to split the original triangle into two triangles, and update the entire triangular mesh.

3. The linear MARS method according to claim 1, characterized in that: Traverse each triangle in the triangular mesh. If the area of ​​a triangle is less than the threshold and the cosine value of the maximum angle of the triangle is less than the threshold, no adjustment is made. If the area of ​​a triangle is less than the threshold and the cosine value of the maximum angle of the triangle is greater than or equal to the threshold, further judgment is made: if the triangle has at least two sides less than the threshold, the triangle is shrunk into a tracer point and the entire triangle mesh is updated; if only one side of the triangle is less than the threshold, the side is shrunk into a tracer point and the entire triangle mesh is updated; if each side of the triangle is greater than the threshold, the triangle is tried to be optimized by side swapping and the entire triangle mesh is updated; if all three sides cannot be optimized by side swapping, the shortest side is shrunk into a tracer point and the entire triangle mesh is updated.

4. The linear MARS method according to claim 3, characterized in that: The shrinking of the triangle into a tracing point is to delete the three vertices of the triangle and add the center of gravity of the triangle as the tracing point.

5. The linear MARS method according to claim 3, characterized in that: The shrinking of the edge into a tracer point is to delete an endpoint with a relatively small angle between the maximum normal vectors of the two end points of the edge and the corresponding multiple triangular planes.

6. The linear MARS method according to claim 3, characterized in that: The side-exchange optimization is that for any side of a triangle, if two triangles with the side as a common side are combined into a convex quadrilateral, the diagonals of the convex quadrilateral are exchanged.

Images