Numerical simulation method for gas-liquid-solid multiphase coupling of solid with arbitrary complex shape entering water

Abstract

The application discloses a numerical simulation method for gas-liquid-solid multiphase coupling of a solid with an arbitrary complex shape entering water, which comprises the following steps: initializing a gas-liquid interface level set function, a solid boundary level set function and initial states of a fluid and a solid; solving a fluid momentum equation without considering a fluid-solid coupling forcing term in a current Runge-Kutta substep to obtain an intermediate velocity field; integrating pressure and viscous stress on a solid surface to obtain a resultant force and a resultant moment acting on the solid; inputting the resultant force and the resultant moment into a solid six-degree-of-freedom dynamics equation set based on Hamiltonian mechanics and Euler parameters, solving and updating a center-of-mass position, a center-of-mass velocity, an attitude, an angular momentum and an angular velocity of the solid; reconstructing the solid boundary level set function by using a ray casting algorithm; constructing a solid velocity field according to the updated center-of-mass velocity and the angular velocity of the solid; calculating an immersed boundary force source term; updating the velocity field, a pressure field and the gas-liquid interface level set function; repeating the above substep and advancing a time step to obtain a flow field distribution, gas-liquid interface evolution and a solid motion state in a water entering process. The application is suitable for simulation scenes of a free surface and large-amplitude motion of a solid with an arbitrary complex shape.

Description

Technical Field

[0001] This invention relates to multiphase flow numerical simulation technology in fluid mechanics, specifically to a numerical simulation method for multiphase coupling of gas, liquid, and solid in water with arbitrarily complex shapes. Background Technology

[0002] In marine engineering, ship and coastal structure design, and multiphase flow dynamics research, the fluid-structure interaction (FSI) is a crucial scientific and engineering problem. This is especially true when solid structures (such as ships, offshore wind turbine platforms, and floating aquaculture structures) move in waves, involving complex physical processes such as free surface deformation, six-degree-of-freedom motion of the solid, and strong nonlinear coupling of wave loads. Accurate simulation of such gas-liquid-solid three-phase coupling problems is of great significance for solid structure safety assessment, motion response prediction, and optimal design.

[0003] Currently, numerical methods for simulating fluid-solid interactions can be mainly classified into three categories: body-fitted mesh methods, overlapping mesh methods, and immersed boundary methods. Body-fitted mesh methods (such as the Arbitrary Lagrange-Eulerian method, ALE) use a dynamically adjusted mesh around the structure, which can effectively resolve boundary layer flows. However, they are prone to mesh distortion when handling large-scale motions, leading to numerical instability or even computational interruption. Overset mesh methods, by combining a fixed background mesh with a body-following mesh, can effectively simulate large-scale motions, but require interpolation in the overlapping regions to transfer information, introducing interpolation errors and significantly increasing computational cost and implementation complexity. Immersed boundary methods (IBM), as a non-body-fitted mesh method, embed the structure into a fixed Eulerian mesh in the form of force source terms, offering advantages such as simple mesh generation and ease of handling complex geometry and large-scale motions. Among them, the Direct Forcing IBM method further improves computational efficiency and stability by directly applying velocity boundary conditions.

[0004] However, existing immersion boundary methods still face problems such as significant spurious velocities, blurred interface transition layers, and easy instability in simulation of large rotations of solid structures when dealing with the coupling of free-moving solids and gas-liquid two-phase flows, which limit the accuracy and reliability of wave load and motion response calculations. Summary of the Invention

[0005] This invention provides a numerical simulation method for gas-liquid-solid multiphase coupling of solids with arbitrary complex shapes entering water, in order to solve the technical problems of significant false velocities at the fluid-solid interface, easy instability in simulation of large rotations of solid structures, and low accuracy and efficiency of coupling calculations in existing methods when simulating wave-solid interactions.

[0006] The technical solution to achieve the purpose of this invention is as follows:

[0007] A numerical simulation method for multiphase coupling of gas, liquid, and solid phases in water with solids of arbitrary complex shapes includes:

[0008] Step S1: Initialize the parameters for the iterative solution of the gas-liquid-solid coupling calculation model, including: initializing the flow field containing air, water, and free liquid surfaces on a fixed orthogonal staggered grid; defining the gas-liquid interface level set function; defining the solid boundary level set function to describe the solid surface; and initializing the fluid velocity, pressure, and the solid's center of mass position, center of mass velocity, attitude, angular velocity, and angular momentum.

[0009] Step S2: In the current Runge-Kutta substep of solving the gas-liquid-solid coupling calculation model, the fluid control equations are discretized using the finite difference method. Momentum prediction is performed without adding the submerged boundary force source term to obtain the intermediate velocity field.

[0010] Step S3: Based on the current flow field pressure and viscous stress tensor, perform area integration on the discrete triangular mesh representing the solid surface to calculate the resultant force and resultant moment acting on the solid.

[0011] Step S4: Input the resultant force, resultant torque and external force into the solid six-degree-of-freedom dynamic equation system based on Hamiltonian mechanics and Euler parameters, solve and update the solid's center of mass position, center of mass velocity, attitude, angular momentum and angular velocity;

[0012] Step S5: Based on the updated position and orientation of the solid centroid in step S4, recalculate the solid boundary level set function corresponding to the current position of the solid using the ray projection algorithm to achieve dynamic updating of the fluid-solid interface.

[0013] Step S6: Based on the updated solid centroid velocity and angular velocity, construct a solid velocity field compatible with the solid motion on the Eulerian mesh nodes;

[0014] Step S7: Based on the solid boundary level set function, the intermediate velocity field obtained in step S2, and the solid velocity field obtained in step S6, calculate the volume force source term in the immersion boundary method;

[0015] Step S8: Add the volume force source term calculated in step S7 to the predicted momentum or intermediate velocity field, and obtain the velocity field and pressure field that satisfy the mass conservation and fluid-solid interface no-slip conditions through pressure correction, and update the gas-liquid interface horizontal set function according to the updated velocity field.

[0016] Step S9: Based on the initialization parameters, repeat steps S2 to S8 until all Runge-Kutta substeps within the current physical time step are completed, and then proceed to the next time step. Repeat steps S1 to S9 to obtain the flow field velocity, pressure distribution, gas-liquid interface position, and solid motion state during the solid entry process.

[0017] Furthermore, in step S1, the gas-liquid interface level set function is defined as follows:

[0018] The gas-liquid interface level set function Defined as the signed shortest distance from any point in the computational domain to the gas-liquid interface; it takes a positive value in the liquid region, a negative value in the gas region, and a zero value at the gas-liquid interface, based on the gas-liquid interface level set function. Determine the density and kinematic viscosity distributions of the gas and liquid phases.

[0019] Further, in step S1, defining the solid boundary level set function includes:

[0020] Based on a solid STL geometry file and a ray casting algorithm, the shortest distance from the background mesh nodes to the solid surface is calculated, and this shortest distance is converted into a solid boundary sign distance function Φ. s , as a level set function of the solid boundary, where positive values ​​represent the fluid region, negative values ​​represent the solid interior, and zero values ​​represent the fluid-solid interface.

[0021] Further, step S3 includes:

[0022] The fluid force and torque are obtained by integrating the flow field pressure and viscous stress along the discrete triangular surface of the solid. The resultant force and resultant torque are obtained by superimposing the gravity, mooring force or other external forces with the fluid force and torque.

[0023] Step S4 includes:

[0024] Euler parameters are used to represent the attitude of the solid. Hamiltonian mechanics is used to establish a first-order dynamic equation system for the translation and rotation of the solid. The solid dynamic equation system is then integrated over time to obtain the instantaneous position of the center of mass, the velocity of the center of mass, the attitude, the angular momentum and the angular velocity of the solid.

[0025] Furthermore, in step S5, the spatial coordinates of each triangular facet in the STL geometry file are first transformed, and then the solid boundary level set function Φ is recalculated on a fixed Eulerian background mesh using a ray casting algorithm. s The solid boundary level set function Φ s Determined by the current position and orientation of the solid, not by the horizontal set function of the gas-liquid interface. renew.

[0026] Furthermore, in step S7, the volumetric force source term is calculated using density and kinematic viscosity interpolation methods and a discrete form; wherein, the activation range of the volumetric force source term is determined by the solid boundary level set function Φ. s It is determined that the continuous interpolation of density and kinematic viscosity is determined by the gas-liquid interface level set function. This method eliminates the artificial density gradient caused by abrupt density changes at the fluid-structure interface. This improvement effectively suppresses spurious velocities near the interface caused by non-physical density differences, allowing for the use of thinner, smoother transition layers (e.g., ε-column). s =0.6Δx), thereby significantly improving the flow resolution and calculation accuracy of fluid load near the fluid-solid interface while ensuring numerical stability.

[0027] Furthermore, the volumetric force source term is:

[0028] f (k+1) ≈ H s (Φ s * ) · [P(u * ) - u * ] / (2α k Δt)

[0029] Among them, f (k+1) For the submerged boundary force source term of the (k+1)th Runge-Kutta substep, u * For the intermediate velocity field obtained in the prediction step, P(u * α is the projection operator that projects the intermediate velocity field onto the solid velocity field. k Φ is the coefficient of the k-th Runge-Kutta substep, Δt is the physical time step, and Φ s * H is the solid boundary level set function updated in the current substep. s (Φ s * ) is the solid boundary Heaviside function, used to smoothly activate the submerged boundary force source term;

[0030] The continuous interpolation of density and kinematic viscosity is:

[0031] ρ = ρ w H gl ( ) + ρ a [1 - H gl ( )];

[0032] ν = ν w H gl ( ) + ν a [1 - H gl ( )];

[0033] Where ρ is the density of the mixed fluid, and ν is the kinematic viscosity of the mixed fluid; ρ w and ρ a Densities of the aqueous phase and gas phase, respectively, ν w and ν a These are the kinematic viscosities of the aqueous phase and the gas phase, respectively; H gl ( ) is a smooth Heaviside function for gas-liquid property interpolation, which approaches 1 in the liquid phase region and approaches 0 in the gas phase region.

[0034] Furthermore, in step S2, the convection term in the fluid momentum equation is discretized using a fifth-order weighted essentially non-oscillatory scheme, and the diffusion term in the fluid momentum equation is discretized using a second-order central difference scheme; in step S4, a third-order low-storage Runge-Kutta method is used to perform time integration on the six-degree-of-freedom dynamic equations of the solid; in step S9, the fluid control equations and the solid dynamic equations are synchronously advanced according to the third-order low-storage Runge-Kutta method.

[0035] Compared with the prior art, the present invention has the following beneficial effects:

[0036] In numerical simulation studies of fluid-structure interaction problems of freely floating structures in waves, this invention proposes an improved direct force immersion boundary method. This method achieves high-precision and high-efficiency full-process simulation of the large-scale motion and strongly nonlinear fluid interaction of solids with complex geometries under wave action. By introducing a novel density interpolation method, this method effectively eliminates artificial density gradients near the fluid-structure interface and significantly suppresses the generation of spurious velocities, thereby greatly improving the accuracy and numerical stability of wave load calculation and motion response prediction. Furthermore, by using Euler parameters (quaternions) combined with Hamiltonian mechanics to describe the six-degree-of-freedom motion of the solid, it completely avoids… The gimbal lock problem in the traditional Euler angle method can simulate arbitrary large-amplitude rotations of solids without singularity. In addition, the sharp interface capture technique based on a fixed background mesh and ray casting algorithm avoids the complexity and computational overhead caused by dynamic mesh deformation or overlapping mesh interpolation, simplifies preprocessing, and enhances the robustness of the method. The strongly coupled solution framework established in this invention realizes the synchronous advancement of fluid and structure in time, and can accurately simulate the complete dynamic process from wave impact, large-scale structural motion to mooring system response, providing an efficient and reliable numerical tool for the design optimization and safety assessment of floating structures in marine engineering. Attached Figure Description

[0037] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0038] Figure 1 This is a schematic diagram of the calculation process of an interface capture method based on IBM and the Level Set method, provided for an embodiment of the present invention.

[0039] Figure 2 (a) shows the STL geometric model, mesh, and level set function graph. Figure 2 In the middle (b), the symbolic distance field and zero level set of the fluid-structure interface are shown. The white circles represent STL geometry.

[0040] Figure 3 A schematic diagram illustrating the coordinate system and related conventions used in the processing method for three-dimensional solid rotation. Figure 3 (a) is a schematic diagram of the inertial coordinate system (x,y,z) and the body coordinate system (ξ,η,ζ); Figure 3 (b) is a diagram of Euler angles. Figure 3 (c) is a schematic diagram of Euler parameters.

[0041] Figure 4 This is a schematic diagram of the solid density processing strategy in an embodiment of the present invention. Figure 4 (a) shows the density distribution of discontinuous solids. Figure 4 (b) is a density distribution diagram of a continuous fluid. Figure 4 (c) is a schematic diagram of the solid Heaviside function, which represents the distribution of the transition layer. The white circles represent the solid outline.

[0042] Figure 5 This is a schematic diagram of a three-dimensional simulation model of a non-rotating body freely falling into water, provided as an embodiment of the present invention.

[0043] Figure 6 This is a comparison chart of experimental and numerical simulation results of the evolution of non-rotating bodies entering water at different times in an embodiment of the present invention. Detailed Implementation

[0044] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of them. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0045] In order to at least solve the technical problems proposed in this invention, embodiments of this invention provide a coupled numerical simulation method for gas-liquid-solid multiphase flow based on the coupling of the immersion boundary method and the Level Set method.

[0046] Figure 1 A flowchart illustrating a numerical simulation method for coupled gas-liquid-solid multiphase flow based on the submerged boundary method and the Level Set method, as provided in this embodiment of the invention, is shown below. Figure 1 As shown, the method provided in this embodiment of the invention includes the following steps:

[0047] S1: Initialize all variables, or obtain the values ​​of fluid (fluid velocity and pressure) and solid (center of mass position, center of mass velocity, attitude, angular velocity, and angular momentum) from the previous time step for the next calculation. First, a fixed staggered grid needs to be established to cover the entire computational domain, including the air domain, water domain, and regions where solids may move.

[0048] S101: Initialization of flow field variables:

[0049] Velocity field initialization: Initialize the fluid velocity vector u across the entire computational domain. Typically, at the initial moment (t=0), the fluid is considered stationary, i.e., u = 0. For water entry problems, the initial velocity of the solid impacting the water surface can be calculated based on the set release height, but the flow field velocity is still set to zero.

[0050] S102: Pressure field initialization:

[0051] The pressure field p is initialized based on the hydrostatic pressure distribution. In the air domain, the pressure is typically set to atmospheric pressure (or a relative pressure of 0); in the water domain, the pressure increases linearly with water depth, i.e. , where h is the water depth. In the formula, ρ is the fluid density, and g is the acceleration due to gravity.

[0052] S103: Initialization of the gas-liquid interface level set function:

[0053] Define the gas-liquid interface level set function Γ(t) is considered as the signed distance function of the gas-liquid interface. The zero isosurface of (x,t). For any point in the computational domain, It takes a positive value in the liquid region, a negative value in the gas region, and a zero value at the gas-liquid interface. Its mathematical expression is:

[0054] (x,t)>0: Liquid phase region;

[0055] (x,t)=0: Gas-liquid interface Γ(t);

[0056] (x,t)<0: Gas phase region;

[0057] Based on the initial position of the free surface, calculate the signed shortest distance from each grid node to that surface, which is... The initial value.

[0058] S104: Initialization of solid boundary level set function:

[0059] Define the solid boundary level set function Φ s Used to describe the surface of a solid. For example... Figure 2 As shown in a), the solid boundary is represented by a discrete STL (Stereolithography) geometric file, which consists of a series of triangular facets capable of accurately describing arbitrarily complex solid shapes. To identify solid regions within a fixed Eulerian mesh, a ray-casting algorithm is used to calculate the shortest distance from each background mesh node to the surface of the STL triangular mesh. Depending on whether the point is located inside or outside the solid, it is assigned a negative or positive value, thus constructing the signed distance function Φ of the solid boundary throughout the entire Eulerian domain. s .

[0060] S105: Solid state initialization:

[0061] Set the initial state of the solid;

[0062] Center of mass position c: Define the initial coordinates of the solid's center of mass in the inertial coordinate system as c = (c0, c1, c2). T [m]

[0063] Center of mass velocity U c Define the initial translational velocity U of the solid's center of mass. c For the problem of free fall into water, the initial velocity is usually 0 or the impact velocity calculated based on the release height.

[0064] Attitude initialization: such as Figure 3 As shown, Figure 3 Figure (a) shows a schematic diagram of the inertial coordinate system (x, y, z) and the body coordinate system (ξ, η, ζ); where CM represents the origin of the body coordinate system, c represents the coordinates of the solid's center of mass, and x... b Let x represent the coordinates of a point in a solid in the body coordinate system, and let x represent the coordinates of that point in the inertial frame. Figure 3 (b) is a schematic diagram of Euler angles, where the roll angle Φ, pitch angle Θ, and yaw angle Ψ represent rotations about the x, y, and z axes, respectively. Figure 3(c) is a schematic diagram of Euler parameters, where e = (e0, e1, e2, e3). T Let be a quaternion, n be the unit vector along the rotation axis, and θ be the rotation angle. Using Euler parameters (quaternions) e to describe the initial spatial attitude of the solid avoids the "gimbal lock" problem that may occur with traditional Euler angle methods. Given Tait-Bryan Euler angles (roll Φ, pitch Θ, yaw Ψ), the initial values ​​of the quaternions can be calculated using the following formula:

[0065]

[0066] In the formula, c(·) and s(·) represent the cosine function and the sine function, respectively. Angular velocity ω b Initialization with angular momentum h: Set the initial angular velocity ω in the body coordinate system. b,0 [rad / s] (usually 0 at the initial moment of entry into the water), and calculate the initial value h0=Iω of the angular momentum vector in the body coordinate system based on the moment of inertia tensor I [kg·m²]. b,0 .

[0067] Through the above steps, all necessary variables have been initialized, preparing for the subsequent Runge-Kutta time-advancement loop.

[0068] S2: Flow field prediction step: In the current Runge-Kutta substep, the fluid momentum equation is solved first without considering the fluid-structure interaction forcing term to obtain the intermediate velocity field.

[0069] Begin the iteration of the k-th Runge-Kutta substep. In the explicit third-order low-storage Runge-Kutta scheme of this method, the fluid momentum equation is first solved without considering the submerged boundary force source terms of fluid-structure interaction to obtain an intermediate velocity field. The formula for calculating this prediction step is:

[0070]

[0071] Where the superscript (k-1) represents the value of the previous substep, u * For the intermediate velocity field, u (k-1) For the velocity field of the previous substep, p (k-1) Let be the pressure field of the previous substep, where ρ and ν are the fluid density and kinematic viscosity, respectively, and g is the acceleration due to gravity; α1=4 / 15, α2=1 / 15, α3=1 / 6, γ1=8 / 15, γ2=5 / 12, γ3=3 / 4, ζ1=0, ζ2=−17 / 60, ζ3=−5 / 12, k=1, 2, 3. Corresponding to step S2, the convection term in the fluid momentum equation is discretized using a fifth-order WENO scheme, and the diffusion term is discretized using a second-order central difference scheme.

[0072] S3: Fluid load integral: Based on the current fluid pressure field and viscous stress tensor, perform surface integral on a discrete triangular mesh representing the solid surface to calculate the total fluid force and torque acting on the solid.

[0073] Based on the current flow field pressure p and viscous stress tensor τ, the total force exerted by the fluid on the solid is accurately calculated by performing surface integration on a discrete triangular mesh representing the solid surface. and torque The calculation formula is as follows:

[0074]

[0075]

[0076] Where n is the unit normal vector pointing outward from the solid surface, and r i Let ΔΩ be the distance vector from the solid's center of mass to the center of the i-th triangular element. i Let ρ be the area of ​​the i-th triangular element, N be the total number of triangular elements on the solid surface; p be the fluid pressure, ρ be the fluid density, ν be the fluid kinematic viscosity, and τ be the velocity gradient symmetric tensor, ρντ i This is the corresponding viscous stress tensor. This step accurately calculates the forces exerted by the flow field on the solid, providing input for solving the solid's equations of motion. If other external forces, such as mooring forces, exist, they are also calculated and superimposed in this step.

[0077] S4: Solving the solid dynamic equations and updating the solid state: Input the fluid load and external forces such as gravity and mooring force into the solid six-degree-of-freedom dynamic equations based on Hamiltonian mechanics and Euler parameters, solve and update the solid's center of mass position, center of mass velocity, attitude, angular momentum and angular velocity.

[0078] The fluid loads (as well as external forces such as gravity and mooring force) calculated by S3 are substituted into the six-degree-of-freedom dynamic equations of the solid established based on Hamiltonian mechanics and Euler parameters to solve and update the motion state of the solid. This dynamic equation set consists of the following first-order ordinary differential equations:

[0079]

[0080]

[0081]

[0082]

[0083] Where c is the position of the centroid, p c Let m be the momentum along the center of mass. s For the mass of a solid, F iThe resultant force acting on the solid is e = (e0, e1, e2, e3). T Let E be the Euler parameter (unit quaternion), and G be 3×4 matrices composed of quaternion components. Let Ġ be the time derivative of G, and R = EG. T Let M be the rotation matrix from the body coordinate system to the inertial coordinate system, I be the moment of inertia tensor in the body coordinate system, h be the angular momentum in the body coordinate system, and M be the rotation matrix from the body coordinate system to the inertial coordinate system. b Let ω be the resultant torque in the body coordinate system. b =I -1 h represents the angular velocity in the body coordinate system. This formulation avoids the gimbal lock singularity problem that occurs with traditional Euler angles at pitch angles of ±90°, and can simulate the rotational motion of solids with arbitrary large amplitudes without singularity. The translational and rotational equations of the solid are constructed based on the Hamiltonian mechanics framework and are integrated over time using the same explicit third-order low-storage Runge-Kutta scheme as the fluid control equations. This achieves complete synchronization and strong coupling between the flow field evolution and the solid motion in time. These equations are integrated over time using the same explicit third-order low-storage Runge-Kutta scheme as the fluid solver, thus achieving complete synchronization and strong coupling between the flow field evolution and the solid motion in time.

[0084] S5: Dynamic Update of Fluid-Solid Interface: Based on the updated position and orientation of the solid centroid, the solid boundary level set function Φ, used to identify the solid region, is recalculated and updated on a fixed Eulerian background mesh using a ray casting algorithm. s This ensures that the fluid-solid interface is consistent with the actual position of the solid.

[0085] In S5, accurate characterization of the fluid-structure interface is crucial. A ray casting algorithm is employed to directly calculate and signify the shortest distance from fixed Eulerian mesh nodes to the triangular mesh on the solid STL surface (negative inside the solid, positive outside). This constructs a sharp, sign-distance function independent of the background mesh size, i.e., the solid boundary level set function Φ used to identify the solid. s The solid boundary level set function is used to define the smooth solid boundary Heaviside function H. s (Φ s The function is used to continuously distinguish between fluid and solid regions and control the application range of the submerged boundary force source term; it is not used for interpolation of gas and liquid phase properties. The notation for the solid boundary level set function is: Φ s <0 indicates the interior of a solid, Φ s =0 indicates a fluid-structure interface, Φ s >0 indicates the fluid region; H s (Φ s )≈1 represents a solid region, H s (Φ s )≈0 represents the fluid region. Where, H sLet ε be the solid boundary heaviside function. s This is the thickness parameter of the solid boundary transition layer.

[0086] The updated quaternion calculated based on S4 and the position of the center of mass Update the spatial position and orientation of the solid. For each point t (represented as (ξ,η,ζ)) on a triangular facet in the STL geometry file, its new coordinates in the inertial coordinate system are... It can be calculated using the following formula:

[0087]

[0088] In the formula, r t Let ξ be the position vector of a point on the triangular surface element in the inertial coordinate system, c be the position vector of the solid's center of mass, R be the rotation matrix, and t = (ξ, η, ζ). T Let be the position vector of this point in the body coordinate system. After the coordinates of all triangular elements have been updated, the ray casting algorithm is used again to recalculate the signed distance from the background mesh nodes to the solid surface based on the new solid centroid position and orientation, thereby updating the solid boundary level set function Φ used to identify the solid region. s This ensures that the fluid-solid interface always precisely matches the actual position of the solid.

[0089] S6: Construction of solid-induced velocity field: Based on the updated solid centroid velocity and angular velocity, a solid velocity field compatible with the solid's kinematic state is constructed on the Eulerian grid nodes.

[0090] Based on the updated solid center of mass velocity U c and the angular velocity ω in the inertial frame (given by the angular velocity ω in the body coordinate system) b (Obtained through rotation matrix transformation), a kinematically compatible solid velocity field u is constructed on the Euler mesh nodes. s For any point in space, the corresponding solid velocity is the superposition of the translational velocity and rotational velocity of its center of mass, calculated as follows, where x is the position vector of the Eulerian grid point, c is the position vector of the solid's center of mass, and P(u * The projection operator is used to project the intermediate velocity field onto the solid velocity field.

[0091]

[0092] S7: Calculation of submerged boundary force source term: based on the updated solid boundary level set function Φ s The intermediate velocity field and the solid velocity field are used to calculate the volumetric force source term in the direct force immersion boundary method.

[0093] The force source term f is calculated in a discrete form, and its activation range is determined by the solid boundary level set function Φ. s Determined; density and kinematic viscosity are calculated using the current gas-liquid interface level set function. Perform continuous interpolation. Based on the updated solid boundary level set function Φ s The intermediate flow field velocity obtained from S2 and the solid velocity field u obtained from S6 s Calculate the volumetric force source term f in the direct force immersion boundary method. The activation of the force source term is determined by the solid boundary level set function Φ. s Control, in its discrete form, is approximately as follows:

[0094] f (k+1) ≈ H s (Φ s * ) · [P(u * ) - u * ] / (2α k Δt)

[0095] Among them, f (k+1) For the submerged boundary force source term of the (k+1)th Runge-Kutta substep, u * For the intermediate velocity field obtained in the prediction step, P(u * α is the projection operator that projects the intermediate velocity field onto the solid velocity field. k Φ is the coefficient of the k-th Runge-Kutta substep, Δt is the physical time step, and Φ s * H is the solid boundary level set function updated in the current substep. s (Φ s * ) is the solid boundary Heaviside function, used to smoothly activate the force source term within the transition layer near the fluid-solid interface;

[0096] H s (Φ s )≈1 represents a solid region, H s (Φ s )≈0 indicates the fluid region.

[0097] The key improvement of this method lies in its assumption that, when dealing with the physical properties near the fluid-solid interface, the density and kinematic viscosity of the virtual fluid inside the solid remain consistent with and smoothly transition to the properties of the adjacent real fluid (water or air). Throughout the computational domain, density and kinematic viscosity are determined solely by the gas-liquid interface level set function. Calculations, rather than introducing the density of the solid itself; gas-liquid interface level set function With the solid boundary level set function Φ sThey are independent of each other; the former is used for gas / liquid property interpolation, while the latter is used for fluid-solid interface identification and force source term activation.

[0098] ρ = ρ w H gl ( ) + ρ a [1 - H gl ( )]

[0099] ν = ν w H gl ( ) + ν a [1 - H gl ( )]

[0100] Where ρ is the density of the mixed fluid, and ν is the kinematic viscosity of the mixed fluid; ρ w and ρ a Densities of the aqueous phase and gas phase, respectively, ν w and ν a These are the kinematic viscosities of the aqueous phase and the gas phase, respectively; H gl ( ) is a gas-liquid property interpolation type smoothed Heaviside function. The above processing method eliminates the traditional method (ρ( ,Φ s )=ρ s H s (Φ s )+[1-H s (Φ s )]{ρ w H gl ( )+ρ a [1-H gl ( At the fluid-solid interface, the density ρ of the solid itself is introduced. s The resulting artificial density gradient effectively suppresses spurious velocities and allows for the use of thinner transition layers (e.g., ε). s =0.6Δx), which significantly improves the flow resolution and the calculation accuracy of fluid load near the fluid-solid interface.

[0101] S8: Add the volumetric force source term to the intermediate velocity field or predicted momentum, and solve the pressure Poisson equation for pressure correction to obtain the velocity and pressure fields that satisfy mass conservation and fluid-solid boundary conditions. Then, based on the updated velocity field, apply the pressure correction to the gas-liquid interface horizontal set function. Perform transport updates to complete the strongly coupled computation of the current substep.

[0102] Add the force source term f obtained from S7 to the intermediate velocity field or predicted momentum obtained from S2 to obtain a corrected velocity field that takes into account the effects of solid boundaries.

[0103] Subsequently, the pressure Poisson equation is solved to obtain a new pressure correction term and update the velocity field; and the gas-liquid interface level set function is then adjusted based on the updated velocity field. Perform transport updates.

[0104]

[0105] The Poisson equation is solved using an efficient algebraic multigrid preprocessor and the BiCGStab iterative solver to obtain the pressure correction term. Subsequently, using Update the pressure field and velocity field:

[0106]

[0107]

[0108] Where, p (k) and p (k+1) These are the pressure fields before and after pressure correction, u (k+1) For the corrected velocity field, p corr For pressure correction terms. Updated velocity field, pressure field, and gas-liquid interface level set function. It simultaneously satisfies the conservation of momentum, conservation of mass, kinematic conditions of the gas-liquid interface, and no-slip boundary conditions at the fluid-solid interface.

[0109] S9: Iterative Loop and Time Advancement: Iterative loop and time advancement based on initialization parameters: Repeat S2 to S8 until all Runge-Kutta substeps within the current physical time step are calculated; then advance to the next time step, and repeat the entire process until the set total simulation duration is reached.

[0110] After completing the current Runge-Kutta substep, if not all substeps of the third-order Runge-Kutta scheme have been completed, use the updated flow field and gas-liquid interface level set functions. The process continues to iterate over the solid state.

[0111] Once all Runge-Kutta substeps are completed, the current physical time step ends; the calculated flow field velocity, pressure, and gas-liquid interface level set function are then calculated. The position, velocity, attitude, angular velocity, and angular momentum of the solid's center of mass are used as the initial values ​​for the next time step.

[0112] Then proceed to the next time step, repeating S1 to S9 until the set total simulation duration is reached.

[0113] The above S1 to S9 constitute the same calculation process, where S2 to S8 are the fluid-structure interaction solutions within a single Runge-Kutta substep, and S9 is used to control the substep cycle and the physical time step advancement.

[0114] Therefore, the step numbers and execution order of this method are consistent in the claims, the invention description and the specific embodiments, all representing the same stage of the same method flow.

[0115] This embodiment also provides a numerical simulation device, characterized in that it includes: a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps of the numerical simulation method.

[0116] A computer storage medium storing an executable program, the executable program being executed by a processor to implement the steps of the numerical simulation method.

[0117] In one embodiment, to illustrate the effect of the gas-liquid-solid multiphase flow coupling numerical simulation method based on the horizontal set and immersion boundary method provided by the present invention, a non-rotating body vertical water immersion verification experiment is used as an example.

[0118] Figure 5 This is a schematic diagram of a three-dimensional simulation model of a non-rotating, free-falling body entering water, provided as an embodiment of the present invention. Figure 5 As shown, taking the vertical entry of a nonrotating solid into water as an example, with a major-minor axis ratio AR=3, the nonrotating solid is released from rest at a certain height above the liquid surface and undergoes free fall motion under the action of gravity. This simulation refines and captures the gas-liquid-solid multiphase flow coupling characteristics and cavitation evolution process during the nonrotating solid's entry into water. In this embodiment, a self-developed CFD solver based on the finite difference method is used. A cuboid water tank fluid domain is established, with dimensions of 0.5 m × 0.5 m × 0.5 m and a water depth of 0.4 m, enclosed by six faces: four side walls, a bottom surface, and a top surface. The six degrees of freedom, weight, center of gravity position, and moment of inertia of the three-dimensional solid model of the nonrotating solid are set, and the initial position of the nonrotating solid in the flow field is set at a certain height above the free liquid surface. To improve computational efficiency, a fixed orthogonal staggered mesh is used to discretize the computational domain, and local mesh refinement is performed around the free liquid surface and the nonrotating solid to accurately capture the drastic flow changes near the interface. The bottom and four sidewalls of the computational domain were set as no-slip wall boundary conditions, and the top was set as a pressure outlet boundary condition. The non-rotating solid surface was treated using the submerged boundary method to achieve no-slip conditions. The solution time was set to 0.5 s, and the time step was dynamically controlled by the CFL condition to ensure numerical stability.

[0119] Figure 6This document compares experimental and numerical simulation results of cavitation evolution at different times during water immersion in this embodiment of the invention. In this embodiment, the non-rotating body is calculated using the submerged boundary method and the Level Set method for gas-liquid-solid multiphase flow coupling. The experimental and numerical results of cavitation evolution at different times along the major and minor axes are compared as follows: Figure 6 As shown.

[0120] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A numerical simulation method for multiphase coupling of gas, liquid, and solid phases in water with solids of arbitrary complex shapes, characterized in that, include: Step S1: Initialize the parameters for the iterative solution of the gas-liquid-solid coupling calculation model, including: initializing the flow field containing air, water, and free liquid surfaces on a fixed orthogonal staggered grid; defining the gas-liquid interface level set function; defining the solid boundary level set function to describe the solid surface; and initializing the fluid velocity, pressure, and the solid's center of mass position, center of mass velocity, attitude, angular velocity, and angular momentum. Step S2: In the current Runge-Kutta substep of solving the gas-liquid-solid coupling calculation model, the fluid control equations are discretized using the finite difference method. Momentum prediction is performed without adding the submerged boundary force source term to obtain the intermediate velocity field. Step S3: Based on the current flow field pressure and viscous stress tensor, perform area integration on the discrete triangular mesh representing the solid surface to calculate the resultant force and resultant moment acting on the solid. Step S4: Input the resultant force, resultant torque and external force into the solid six-degree-of-freedom dynamic equation system based on Hamiltonian mechanics and Euler parameters, solve and update the solid's center of mass position, center of mass velocity, attitude, angular momentum and angular velocity; Step S5: Based on the updated position and orientation of the solid centroid in step S4, recalculate the solid boundary level set function corresponding to the current position of the solid using the ray projection algorithm to achieve dynamic updating of the fluid-solid interface. Step S6: Based on the updated solid centroid velocity and angular velocity, construct a solid velocity field compatible with the solid motion on the Eulerian mesh nodes; Step S7: Based on the solid boundary level set function, the intermediate velocity field obtained in step S2, and the solid velocity field obtained in step S6, calculate the volume force source term in the immersion boundary method; Step S8: Add the volume force source term calculated in step S7 to the predicted momentum or intermediate velocity field, and obtain the velocity field and pressure field that satisfy the mass conservation and fluid-solid interface no-slip conditions through pressure correction, and update the gas-liquid interface horizontal set function according to the updated velocity field. Step S9: Based on the initialization parameters, repeat steps S2 to S8 until all Runge-Kutta substeps within the current physical time step are completed, and then proceed to the next time step. Repeat steps S1 to S9 to obtain the flow field velocity, pressure distribution, gas-liquid interface position, and solid motion state during the solid entry process.

2. The numerical simulation method according to claim 1, characterized in that, In step S1, the gas-liquid interface level set function is defined as follows: The gas-liquid interface level set function Defined as the signed shortest distance from any point in the computational domain to the gas-liquid interface; it takes a positive value in the liquid region, a negative value in the gas region, and a zero value at the gas-liquid interface, based on the gas-liquid interface level set function. Determine the density and kinematic viscosity distributions of the gas and liquid phases.

3. The numerical simulation method according to claim 1, characterized in that, The step S1 defines the solid boundary level set function, including: Based on a solid STL geometry file and a ray casting algorithm, the shortest distance from the background mesh nodes to the solid surface is calculated, and this shortest distance is converted into a solid boundary sign distance function Φ. s , as a level set function of the solid boundary, where positive values ​​represent the fluid region, negative values ​​represent the solid interior, and zero values ​​represent the fluid-solid interface.

4. The numerical simulation method according to claim 1, characterized in that, Step S3 includes: The fluid force and torque are obtained by integrating the flow field pressure and viscous stress along the discrete triangular surface of the solid. The resultant force and resultant torque are obtained by superimposing the gravity, mooring force or other external forces with the fluid force and torque. Step S4 includes: Euler parameters are used to represent the attitude of the solid. Hamiltonian mechanics is used to establish a first-order dynamic equation system for the translation and rotation of the solid. The solid dynamic equation system is then integrated over time to obtain the instantaneous position of the center of mass, the velocity of the center of mass, the attitude, the angular momentum and the angular velocity of the solid.

5. The numerical simulation method according to claim 1, characterized in that, Step S5 first transforms the spatial coordinates of each triangular facet in the STL geometry file, and then recalculates the solid boundary level set function Φ on a fixed Euler background mesh using a ray casting algorithm. s The solid boundary level set function Φ s Determined by the current position and orientation of the solid, not by the horizontal set function of the gas-liquid interface. renew.

6. The numerical simulation method according to claim 1, characterized in that, In step S7, the volumetric force source term is calculated using density and kinematic viscosity interpolation and a discrete form; wherein, the activation range of the volumetric force source term is determined by the solid boundary level set function Φ. s It is determined that the continuous interpolation of density and kinematic viscosity is determined by the gas-liquid interface level set function. Sure.

7. The numerical simulation method according to claim 6, characterized in that, The volumetric force source term is: f (k+1) ≈ H s (Φ s * ) · [P(u * ) - u * ] / (2α k Δt) Among them, f (k+1) For the submerged boundary force source term of the (k+1)th Runge-Kutta substep, u * For the intermediate velocity field obtained in the prediction step, P(u * α is the projection operator that projects the intermediate velocity field onto the solid velocity field. k Φ is the coefficient of the k-th Runge-Kutta substep, Δt is the physical time step, and Φ s * H is the solid boundary level set function updated in the current substep. s (Φ s * ) is the solid boundary Heaviside function, used to smoothly activate the submerged boundary force source term; The continuous interpolation of density and kinematic viscosity is: p = p w H gl ( ) + p a [1 - H gl ( )]; n = n w H gl ( ) + n a [1 - H gl ( )]; Where ρ is the density of the mixed fluid, and ν is the kinematic viscosity of the mixed fluid; ρ w and ρ a Densities of the aqueous phase and gas phase, respectively, ν w and ν a These are the kinematic viscosities of the aqueous phase and the gas phase, respectively; H gl ( ) is a smooth Heaviside function for gas-liquid property interpolation, which approaches 1 in the liquid phase region and approaches 0 in the gas phase region.

8. The numerical simulation method according to claim 1, characterized in that, In step S2, the convection term in the fluid momentum equation is discretized using a fifth-order weighted essentially non-oscillatory scheme, and the diffusion term in the fluid momentum equation is discretized using a second-order central difference scheme. In step S4, the solid six-degree-of-freedom dynamic equations are integrated over time using a third-order low-storage Runge-Kutta scheme. In step S9, the fluid control equations and the solid dynamic equations are synchronously advanced according to the third-order low-storage Runge-Kutta scheme.

9. A numerical simulation device, characterized in that, include: Memory, processor, and storage on said memory A computer program, wherein the processor executes the computer program to implement the steps of the numerical simulation method according to any one of claims 1-8.

10. A computer storage medium, characterized in that, The computer storage medium stores an executable program, which is executed by a processor to implement the steps of the numerical simulation method according to any one of claims 1-8.

Images