Method, system and medium for predicting air, liquid and solid three-phase coupled elastic thin shell slamming load

Abstract

The application discloses a kind of elastic thin shell slamming load prediction method, system, medium considering gas, liquid, solid three-phase coupling, steps are as follows: S1: construction PLIC-VOF compressible two-phase flow model, capture sub-millimeter level air cushion free surface deformation, roll and break, solve compressible flow control equation;S2: introduce virtual fluid GFM method, obtain high-precision slamming pressure field;S3: build elastic thin shell finite element dynamics model, represent the transient large deformation of the shell of vehicle body;S4: using partition implicit strong coupling algorithm, fluid and structure model are solved synchronously iterated, and output slamming pressure, structure dynamic response and air cushion evolution law.The application discards traditional incompressible fluid, rigid structure assumption, solves compressible air cushion and elastic deformation strong coupling solving problem, can represent the three-phase coupling effect of cross-medium water entry, provides high-precision numerical simulation means for the anti-slamming structure design of vehicle body.

Description

Technical Field

[0001] This invention relates to the fields of shipbuilding and marine engineering, and cross-medium aircraft design technology, specifically to a high-precision numerical simulation method for the water impact problem of a vehicle, and in particular to a method, system, and medium for predicting the impact load of a gas-liquid-solid three-phase coupled elastic thin shell. Background Technology

[0002] Cross-medium vehicles are subjected to enormous impact loads upon entering the water during their traversal of the gas-liquid interface, posing a serious threat to their structural safety and the stability of their entry trajectory. This impact process involves strong nonlinear coupling between the gas, liquid, and solid phases, accompanied by complex physical phenomena such as large deformation, tumbling, and fragmentation of the free liquid surface, as well as the evolution of the compressible air cushion at the vehicle's bottom. Accurately predicting the impact load and the structural dynamic response is a key challenge in the design of cross-medium vehicles. Existing numerical methods suffer from the following shortcomings: most fluid models are based on the incompressibility assumption, neglecting the compressibility of air and failing to accurately simulate the buffering effect of the bottom air cushion during high-speed water entry and its impact on slamming loads; the capture of the gas-liquid interface often employs algebraic VOF methods (such as HRIC and MULES), which suffer from severe numerical dissipation when dealing with sub-millimeter-scale thin air cushions, leading to interface ambiguity and an inability to accurately simulate the compression and escape processes of the air cushion; the vehicle is often simplified as a rigid structure, ignoring the transient large deformation of the elastic thin shell under slamming loads and failing to reveal the essence of fluid-structure interaction; traditional fluid-structure interaction methods are mostly weakly or explicitly coupled, exhibiting poor numerical stability and insufficient synchronization in the millisecond-level response time domain where air cushion pressure and structural deformation are strongly coupled. Therefore, there is an urgent need to develop a cross-medium slamming load prediction method that can simultaneously consider gas-liquid compressibility, high-precision free liquid surface tracking, elastic thin shell structural deformation, and strong fluid-structure interaction effects. Summary of the Invention

[0003] The technical problem to be solved by the present invention is to address the shortcomings of the prior art by providing a method, system and medium for predicting the impact load of an elastic thin shell that takes into account the coupling of gas, liquid and solid phases, so as to achieve accurate prediction of the spatiotemporal distribution of impact pressure, the dynamic response of the elastic thin shell and the evolution process of the compressible air cushion.

[0004] The technical problem to be solved by the present invention is achieved through the following technical solution: a method for predicting the impact load of an elastic thin shell considering three-phase coupling of gas, liquid, and solid phases, the steps of which are as follows:

[0005] S1: Construct two-phase flow control equations that consider gas-liquid compressibility, and use the geometric fluid volume (VOF) method based on piecewise linear interface reconstruction (PLIC) to establish a high-precision free surface tracking model for compressible two-phase flow, which is used to accurately capture the evolution of gas-liquid interface including sub-millimeter-scale air cushions in cross-medium processes.

[0006] S2: Introducing the Virtual Fluid Method (GFM) to precisely maintain the sharp discontinuity of density and pressure at the gas-liquid interface, suppressing the spurious velocity introduced by the dynamic pressure term, and obtaining an accurate slamming pressure field;

[0007] S3: Establish a finite element dynamic model of the elastic thin shell structure to solve the transient large deformation behavior of the thin shell of the aircraft body under impact load;

[0008] S4: Using a partitioned coupling framework and an implicit strong coupling algorithm, the compressible two-phase flow model established in steps S1-S2 and the structural dynamics model established in step S3 are solved synchronously and iteratively to achieve strong coupling calculation of the air cushion compression process and the elastic thin shell deformation process in the millisecond time domain; output the spatiotemporal distribution of impact pressure, structural dynamic response, and evolution characteristics of the compressible air cushion of the elastic thin shell structure during the cross-medium impact process.

[0009] The two-phase flow control equations considering the compressibility of water vapor in S1 include the mass conservation equation, momentum conservation equation, total energy conservation equation, and volume fraction transport equation. The gas-liquid interface is reconstructed and captured using the PLIC geometric VOF method, specifically expressed as follows:

[0010] (1) Mass conservation equation:

[0011] ;

[0012] in, For the density of the mixed fluid, For time, It is a velocity vector. Represents the divergence operator. This is the local partial derivative of density with respect to time;

[0013] (2) Momentum conservation equation:

[0014] ;

[0015] in, For pressure, For viscous stress tensor, It is the gravitational acceleration vector. This is the surface tension volume force term. This represents the gradient operator;

[0016] (3) Total energy conservation equation:

[0017] ;

[0018] in, Total energy per unit mass This is the heat flux density vector;

[0019] The specified unit mass always satisfies:

[0020] ;

[0021] in, Internal energy per unit mass Kinetic energy per unit mass;

[0022] (4) Volume fraction transport equation:

[0023] ;

[0024] in, It represents the volume fraction of the aqueous phase. This refers to the gas phase volume fraction. The interface compression velocity is used to suppress numerical diffusion at the interface and maintain the clarity of the gas-liquid interface; the gas-liquid interface is reconstructed into a piecewise linear interface within the interface unit according to the volume fraction distribution using the PLIC geometric VOF method.

[0025] Furthermore, the physical properties of the mixed fluid satisfy:

[0026] ;

[0027] ;

[0028] in, The density of the aqueous phase, For gas phase density, For mixed fluid dynamic viscosity, For the dynamic viscosity of the aqueous phase, This refers to the dynamic viscosity of the gas phase.

[0029] The aqueous and gas phases satisfy their respective equations of state to establish the coupling relationship between pressure, density, and energy.

[0030] Preferably, the PLIC-based geometric VOF method in S1 includes the following sub-steps:

[0031] S11: Define the volume fraction parameter for each grid cell within the computational domain; define the proportion of the target phase volume in the two-phase flow field to the total volume of the grid cells as the volume fraction. ,in, This indicates that the mesh cell is completely filled with the target phase. This indicates that there is no target phase within this mesh cell. This indicates that the mesh element is an interface element. By initializing the volume fraction of all mesh elements within the computational domain, a discrete volume fraction field is established to characterize the distribution of the two-phase interface.

[0032] S12: Identify interface units and extract local interface information; based on the volume fraction values ​​of each grid unit, filter those that meet the requirements. The mesh elements are used as interface elements; at the same time, the volume fraction distribution information of the interface elements and their adjacent mesh elements is obtained to provide a data basis for subsequent interface normal vector solving and interface geometry reconstruction.

[0033] S13: Calculate the interface normal vector within the interface unit; calculate the interface normal vector based on the spatial gradient of the volume fraction field surrounding the interface unit, wherein the interface normal vector can be expressed as:

[0034] ;

[0035] in, For the interface normal vector, It is the volume fraction. For volume fraction gradient;

[0036] In specific implementation, the volume fraction gradient can be obtained through central difference, least squares method or weighted reconstruction method to improve the accuracy of interface orientation recognition;

[0037] S14: Construct a piecewise linear interface within the interface unit; using the PLIC method, approximate the real phase interface within each interface unit using a linear interface, wherein the linear interface satisfies the following expression:

[0038] ;

[0039] in, Let be the position vector of any point on the interface. The interface intercept constant;

[0040] The orientation of the linear interface is determined based on the interface normal vector obtained in step S13, and the intercept constant is then solved by combining the volume fraction constraints within the interface element. This ensures that the geometric volume of the target phase segmented by the linear interface within the current grid cell is consistent with the target phase volume corresponding to the volume fraction of that grid cell.

[0041] S15: Solve for the interface intercept constant based on the volume fraction constraint; let the total volume of the current interface element be... The target phase volume fraction is Then the linear interface must satisfy the following volume conservation relationship:

[0042] ;

[0043] in, Represents the normal vector and intercept constant The target phase volume intercepted within the current mesh cell by the determined linear interface;

[0044] In practice, the intercept constant can be obtained by analytical solution, table lookup, or iterative solution. To complete the interface geometry reconstruction;

[0045] S16: Construct a swept volume on a grid surface based on the velocity field and calculate the phase volume flux; at a given time step... Within this time step, a sweeping region is constructed based on the normal velocity and area of ​​each grid surface, and the transport volume of the target phase through the grid surface is determined. For any grid surface, its sweeping volume can be expressed as:

[0046] ;

[0047] in, For the normal velocity of the mesh surface, For the area of ​​the grid surface, For time step;

[0048] Furthermore, combining the linear interface positions obtained in S14 and S15, the geometric overlap volume between the swept region and the target phase region in the upstream unit is calculated, which is used as the geometric flux of the target phase through the grid surface.

[0049] S17: Update the volume fraction of each grid cell based on geometric flux; according to the target inflow and outflow volumes of each grid surface, perform a conservation update on the volume fraction of each grid cell, using the following formula:

[0050] ;

[0051] in, and They represent the first The volume fraction of each grid cell at the current time and the next time step. For the first Volume of each grid cell This represents the total outflow volume of the target phase within the current time step for this mesh cell. The total volume of the target phase inflow within the current time step for this mesh cell;

[0052] S18: Apply boundedness correction to the updated volume fraction; to ensure that the volume fraction satisfies the physical constraints, apply boundedness restrictions to the updated volume fraction to ensure that it satisfies:

[0053] ;

[0054] When the updated volume fraction of a certain grid cell is less than 0, it is corrected to 0; when the updated volume fraction is greater than 1, it is corrected to 1, in order to avoid non-physical interface diffusion or spurious phase volume during numerical calculation.

[0055] S19: Repeat the interface reconstruction and volume fraction transport process; in each time step, repeat steps S12 to S18 to update the position and shape of the two-phase interface in real time, thereby realizing continuous capture and evolution calculation of the gas-liquid interface.

[0056] Preferably, the virtual fluid method (GFM) introduced in S2 includes the following sub-steps:

[0057] S21: Identify the gas-liquid interface location; determine the interface unit based on the distribution results of the volume fraction function in step S1, and combine it with the piecewise linear interface reconstructed by the PLIC geometric VOF method to obtain the local position, normal direction, and medium properties on both sides of the interface, which serve as the geometric basis for constructing the virtual fluid state. GFM requires first clarifying the interface location, and then processing the fluids on both sides in the interface neighborhood according to the single-medium format.

[0058] S22: Extract the real fluid state on both sides of the interface; in the adjacent cells of the interface, extract the real fluid state quantities on the liquid side and the gas side respectively. The real fluid state quantities include at least density, pressure, velocity, total energy and the speed of sound determined by the equation of state, for subsequent calculation of local interface interactions.

[0059] S23: Establish a local Riemann problem along the interface normal; using the interface unit normal as the solution direction, project the actual fluid velocities on both sides obtained in step S22 onto the normal direction to construct a local one-dimensional Riemann problem along the interface normal direction to characterize the propagation and interaction process of pressure waves, velocity waves and contact discontinuities at the gas-liquid interface.

[0060] S24: Solve for the intermediate state of the interface; Solve the local Riemann problem established in step S23 to obtain the intermediate interface pressure and the intermediate interface normal velocity; where the intermediate interface pressure is used to characterize the pressure equilibrium state after the interaction of fluids on both sides of the interface, and the intermediate interface normal velocity is used to characterize the propagation velocity or interface velocity along the normal direction at the interface; For compressible multi-medium problems, linearized interface interaction methods or approximate Riemann solvers are often used in the literature to calculate the interface state; In the case of surface tension, a pressure transition term caused by curvature can also be added to the interface pressure condition;

[0061] S25: Construct virtual fluid states on both sides of the interface; based on the intermediate interface pressure and intermediate interface normal velocity obtained in step S24, construct corresponding virtual fluid states on both sides of the interface respectively; the virtual fluid states include at least virtual pressure, virtual normal velocity, virtual density, and virtual total energy. The virtual pressure and virtual normal velocity are determined by the intermediate interface state, and the virtual density can be obtained based on the corresponding medium state equation and combined with isentropic extrapolation or characteristic extrapolation, thereby avoiding non-physical mixing and pressure oscillations caused by direct interpolation across the interface;

[0062] S26: Complete the single-phase discrete template near the interface; using the virtual fluid state constructed in step S25, complete the single-phase calculation template of the unit near the interface in the differential, finite volume or finite reconstruction process, so that the gas phase region and the liquid phase region can be solved independently according to the control equation of their respective medium during numerical discretization, without having to directly call the heterogeneous unit data across the real gas-liquid interface during the discretization process.

[0063] S27: Correct the numerical flux near the interface; Based on the single-phase template completed in step S26, calculate the mass flux, momentum flux and energy flux of the unit near the interface, and substitute the numerical flux into the two-phase flow control equation for discrete update, so as to reduce the non-physical oscillation, interface diffusion and numerical overheating caused by abrupt changes in physical properties near the interface.

[0064] S28: Update the interface position and repeat the virtual fluid construction process. After completing the flow field update at the current time step, continue to update the volume fraction field and the gas-liquid interface position, and repeat steps S21 to S27 for the new interface neighborhood, thereby realizing continuous simulation of interface propagation, wave system reflection and transmission, and strong discontinuous interaction in the transient compressible water-gas two-phase flow process.

[0065] Preferably, the finite element dynamic equation of the elastic thin shell structure in S3 is:

[0066] ;

[0067] in, The overall mass matrix of the elastic thin-shell structure. For the overall damping matrix, For the overall stiffness matrix, Let be the nodal displacement vector. For the node velocity vector, For the nodal acceleration vector, The vector of external loads acting on the structure. This is the fluid-structure interaction load vector acting on the surface of the thin-shell structure.

[0068] Preferably, the implicitly strongly coupled algorithm in S4 includes the following sub-steps:

[0069] S41: At the beginning of each physical time step, the structural displacement converged in the previous time step is passed to the fluid solver as the initial boundary condition;

[0070] S42: The fluid solver solves a compressible two-phase flow model based on the current structural boundary, calculating the fluid dynamic load vectors acting on the fluid-structure interaction interface. ;

[0071] S43: The radial basis function interpolation method is used to interpolate the hydrodynamic load. And conservatively map from the fluid mesh to the structural mesh, while simultaneously transferring structural displacements. A consistent mapping from the structural mesh to the fluid mesh is used to maintain energy balance at the coupling interface;

[0072] S44: The structural solver receives the mapped load, solves the finite element dynamic equations of the elastic thin-shell structure, and obtains the new structural displacements. ;

[0073] S45: Calculate the... Force residual vector after the second iteration and displacement residual vector And calculate their discrete values ​​respectively. Norm:

[0074] ;

[0075] S46: Use the relative convergence criterion to determine whether the iteration has converged:

[0076] ;

[0077] in For the force convergence threshold, This is the displacement convergence threshold; if the above conditions are met, the sub-iteration is exited and the next time step is initiated.

[0078] S47: If convergence fails, update the structural boundary displacements using the interface quasi-Newtonian relaxation method:

[0079] ;

[0080] in The adaptive relaxation factor is then used, and the process returns to step S42 to continue iterating until the convergence criterion is met or the preset maximum number of iterations is reached.

[0081] Preferably, the force convergence threshold and displacement convergence threshold The average value is .

[0082] Preferably, the conservative mapping ensures the conservation of total force through the coupling interface, and the uniform mapping ensures the continuity of the displacement field.

[0083] Preferably, a system for predicting the impact load of a gas-liquid-solid three-phase coupled elastic thin shell is provided. This system is used to execute any of the above methods for predicting the impact load of a gas-liquid-solid three-phase coupled elastic thin shell, including a compressed two-phase flow calculation module for executing steps S1-S2.

[0084] The elastic thin-shell structure calculation module is used to execute the steps in S3.

[0085] An implicitly strongly coupled iterative module is used to execute the steps of S4;

[0086] Within each time step, the compressible two-phase flow calculation module and the elastic thin-shell structure calculation module are driven to perform synchronous iterative solutions; and the spatiotemporal distribution of impact pressure, structural dynamic response, and compressible air cushion evolution characteristics are output.

[0087] Preferably, a computer-readable storage medium stores a computer program that, when executed by a processor, implements the above-mentioned method for predicting impact loads on a gas-liquid-solid three-phase coupled elastic thin shell.

[0088] Compared with the prior art, the beneficial technical effects of the present invention are:

[0089] (1) By combining the geometric VOF method (PLIC) and the virtual fluid method (GFM), a sharp and spurious-free capture of sub-millimeter-level compressible air cushions was achieved in the slamming problem of elastic thin shell structures, solving the key problems of severe numerical dissipation and pressure distortion in traditional methods.

[0090] (2) It breaks through the two traditional assumptions of fluid incompressibility and structural rigidity, and establishes the most complete physical model that includes gas and liquid compressibility and elastic thin shell deformation, which more realistically reflects the physical nature of cross-medium impact of elastic thin shell structure.

[0091] (3) An implicit strong coupling algorithm is adopted to ensure that the numerical calculation has excellent stability and convergence in the strong coupling scenario where air cushion pressure and structural deformation are interdependent, and can accurately reveal the gas-liquid-solid three-phase coupling mechanism in the millisecond time domain.

[0092] (4) It can provide accurate numerical analysis tools for the anti-slamming optimization design and load reduction method research of cross-medium navigation bodies, ship structures, seaplanes and other structures, shorten the research and development cycle and reduce the test cost. Attached Figure Description

[0093] Figure 1 A flowchart of the overall method for predicting the impact load of elastic thin-shell structures that takes into account the coupling of gas, liquid, and solid phases;

[0094] Figure 2 The diagram shows the geometric model and monitoring point layout of the elastic thin-shell curved wedge structure, where (a) is the geometric model and (b) is the monitoring point layout.

[0095] Figure 3 The diagram shows the computational domain meshing and boundary conditions in the embodiment, where (a) is the fluid domain mesh and (b) is the structural domain mesh.

[0096] Figure 4 The diagram shows the interface reconstruction of the geometric VOF method based on PLIC, where (a) is a typical mesh diagram of a single plane cutting, and (b) is a flowchart of the surface interpolation scheme.

[0097] Figure 5 A schematic diagram illustrating the principle of interface state construction in the Virtual Fluid Method (GFM);

[0098] Figure 6 This is an iterative flowchart for an implicitly strongly coupled algorithm.

[0099] Figure 7 The calculated spatiotemporal distribution cloud map of the slamming pressure;

[0100] Figure 8 The calculated displacement distribution cloud map of the elastic thin shell structure;

[0101] Figure 9 The calculated stress time history curve of the elastic thin-shell structure;

[0102] Figure 10 The diagram shows a comparison of the calculated free surface evolution process, where (a) to (d) represent the free surface morphology at different times. Detailed Implementation

[0103] The specific technical solutions of the present invention will be further described below with reference to the accompanying drawings, so as to enable those skilled in the art to further understand the present invention, without constituting a limitation on its rights.

[0104] Example 1: A method, system, and medium for predicting impact loads of an elastic thin shell considering three-phase coupling of gas, liquid, and solid phases. This example uses the vertical entry of an elastic thin shell curved wedge structure into water as an example, but the scope of protection of this invention is not limited to this example.

[0105] This invention addresses the strong transient impact problem occurring during high-speed crossing of the gas-liquid interface in transmedium vehicles, partially thin-shell ship structures, and elastic profiles of seaplanes. It constructs a numerical calculation method that simultaneously considers the compressibility of the gas phase, the compressibility of the liquid phase, and the large deformation response of the elastic thin-shell structure. This method achieves unified calculation of air cushion compression, interface evolution, and structural response during transmedium impact through a strongly coupled solution of the compressible two-phase flow governing equations, the PLIC geometric VOF interface tracking method, the virtual fluid method, and the finite element dynamics model of the elastic thin shell. Figure 1 This is an overall flowchart of the impact load method for an elastic thin-shell structure that takes into account the three-phase coupling of gas, liquid, and solid phases in this invention.

[0106] Example model and parameter settings:

[0107] A calculation method was implemented for the problem of a uniformly vibrating surface wedge entering water vertically and slamming into it. Figure 2 This is a diagram showing the geometric model of the computational model structure and the layout of monitoring points used in this embodiment of the invention. The geometry of the wedge-shaped body is defined by the distance from the keel to the side of the hull. Geometric half-width and average bottom rise angle To define it. The structure is considered as a beam with simply supported boundary conditions from the keel to the side, where , The impact velocity was The density of water Kinematic viscosity of water The structural material is steel, and the elastic modulus is... mass density , plate thickness Because the geometric VOF method is used to capture the free surface, the air phase is also included in the current hydrodynamic calculations. Air density. air kinematic viscosity .

[0108] The profile of a two-dimensional curved surface wedge can be represented by the following analytical formula:

[0109] ;

[0110] ;

[0111] Where β is the average heave angle along the surface of the wedge, and B is the width of the curved wedge. χ=1 represents a linear wedge, χ>1 corresponds to a concave wedge with negative curvature, and χ<1 corresponds to a convex wedge with positive curvature. The embodiment selected in this patent is a convex wedge with an average heave angle of β=10° and a curvature factor of χ=0. To capture the mechanical response of the curved wedge during water entry, its arc length is divided into six equal parts, and five of these division points are selected as characteristic monitoring points for impact pressure, deflection, and stress. The measuring points are numbered sequentially according to the order of water entry: the position that first contacts the water is defined as measuring point 1, and the rest are numbered sequentially as measuring points 2, 3, 4, and 5.

[0112] Grid generation:

[0113] Considering structural symmetry, only the finite element model of the left elastic wedge is established, where the hinged beams at both ends correspond to... Figure 2 The length of the annotation is The wedge-shaped base plate. Figure 3 This is a schematic diagram of the computational domain meshing and boundary conditions in an embodiment of the present invention, where (a) is the fluid domain mesh and (b) is the structural domain mesh. The fluid domain uses an unstructured hexahedral mesh, with local refinement in the wedge's entry into the water region. The mesh size on the wedge surface is 0.0025m to capture the evolution of the free surface. The total number of meshes in the fluid domain is approximately 94,000. The elastic wedge is discretized using linear reduced integral solid elements (C3D8R), and the structural domain mesh size is uniformly set to 0.005m. The structural domain contains 200 elements.

[0114] In two-way fluid-structure interaction calculations, the initial and boundary conditions of the structural domain are set as follows: an initial velocity field is defined for all nodes, and the keel and side nodes are constrained to ensure that they move at a velocity of 100 rpm. The wedge falls vertically at a constant velocity while its horizontal displacement is restricted, ensuring that the entire wedge enters the water vertically at a constant velocity. Furthermore, based on the fundamental assumptions of the two-dimensional numerical model, corresponding constraints are imposed on the out-of-plane translational degrees of freedom of all nodes.

[0115] S1: Establish a high-precision free surface tracking model for compressible two-phase flow. Figure 4 This is a schematic diagram of interface reconstruction using the PLIC-based geometric VOF method in this invention, where (a) is the actual gas-liquid interface and (b) is the piecewise linear interface reconstruction result within the mesh cell. For specific implementation details of this step, please refer to... Figure 2 The principle shown.

[0116] S11: Volume fraction initialization: Define the volume fraction of the aqueous phase within the computational domain. Initial air domain waters The gas-liquid interface is located on a horizontal plane. The initial position of the elastic thin-shell curved wedge structure is set at the critical position above the free surface.

[0117] S12: Identify Interface Units: At each time step, scan all grid units and filter out... The unit is used as the interface unit.

[0118] S13: Calculate the interface normal vector: For each interface element, calculate the volume fraction gradient using the least squares method based on the volume fraction field of its neighboring elements. This leads to the unit normal vector. The normal vector points towards the gas phase (air) side.

[0119] S14-S15: Constructing Piecewise Linear Interfaces: Using the PLIC method, a linear interface is constructed within each interface unit. The intercept is solved iteratively. This makes the target phase volume separated by the linear interface equal to .

[0120] S16-S17: Calculate geometric flux and update volume fraction: Construct a swept volume for each grid surface based on the velocity field, calculate the geometric overlap volume between the swept volume and the target phase region in the upstream cell, and obtain the volume flux. Then, according to the formula... Update the volume fraction of each unit.

[0121] S18: Boundedness Correction: For the updated Limiting the amplitude: To prevent numerical overflow.

[0122] S19: Cyclic execution: Repeat S12~S18 at each time step to track the evolution of the gas-liquid interface in real time.

[0123] In this embodiment, an air layer (air cushion) is formed during the pressing down of the bottom of the elastic thin-shell curved wedge. The PLIC-VOF method clearly captures the evolution of this air cushion, with a sharp interface and no visible numerical diffusion.

[0124] S2: Introducing the Virtual Fluid Method (GFM) to precisely maintain sharp discontinuities in density and pressure at the gas-liquid interface and suppress spurious velocities introduced by the dynamic pressure term. Figure 5 This is a schematic diagram illustrating the virtual fluid (GFM) method used in this invention to construct a virtual fluid along the interface normal. Figure 5 As shown, by constructing virtual fluid states on both sides of the interface, the density gradient appears only on the boundary between the real and virtual values, thereby suppressing spurious velocities.

[0125] S21: Identify the interface location: Based on the volume fraction field and PLIC interface reconstructed by S1, determine the precise location and normal of each interface element.

[0126] S22: Extract the actual fluid state on both sides of the interface: Extract the liquid phase pressure in the adjacent cells on both sides of the interface. Gas phase side pressure Liquid phase side normal velocity Gas phase side normal velocity Liquid phase density gas phase density Liquid phase sound velocity Harmony and the speed of sound Furthermore, the acoustic impedances of the liquid phase and the gas phase are defined as follows: , .

[0127] S23-S24: Establishing the local Riemann problem and solving for the intermediate state of the interface: Calculating the intermediate pressure of the interface using a linearized Riemann solver. and intermediate normal velocity .

[0128] ;

[0129] ;

[0130] S25: Construct virtual fluid states: Construct virtual fluid states on both sides of the interface. , .

[0131] S26-S27: Complete the discrete template and correct the flux: Fill the single-phase discrete template of the adjacent cells of the interface with virtual fluid state, and then calculate the numerical flux of the liquid phase and gas phase respectively, and substitute it into the momentum equation and energy equation to solve.

[0132] S28: Update the interface and repeat: Repeat the above process at each time step to achieve stable calculation of the compressible two-phase flow.

[0133] S3: Establish a dynamic model of the elastic thin-shell structure. The dynamic equations of the elastic wedge are established using the finite element method:

[0134] in, The overall mass matrix of the elastic thin-shell structure. For the overall damping matrix, For the overall stiffness matrix, Let be the nodal displacement vector. For the node velocity vector, For the nodal acceleration vector, The pressure distribution at the fluid-structure interaction interface is calculated in steps S1-S2.

[0135] S4: Implicitly strongly coupled iterative solution. Figure 6This is an iterative flowchart of the implicitly strongly coupled algorithm in this invention. For example... Figure 6 As shown, within each physical time step, the fluid solver and structure solver iterate repeatedly until the force and displacement residuals satisfy the convergence criteria. A partitioned coupling framework is adopted, connecting the fluid solver and the structure solver CalculiX, which are based on OpenFOAM secondary development, through the preCICE library. The time step for the fluid-structure interaction calculation is... The following sub-iterations are executed within each physical time step ( ):

[0136] S41: Initial Prediction: The structural displacement converged in the previous time step As a fluid boundary.

[0137] S42: Fluid Solution: Call the solver from steps S1-S2 to calculate the pressure distribution at the fluid-solid interface. .

[0138] S43: Data Mapping: Radial basis function (RBF) interpolation is used to conservatively map fluid nodal forces to structural nodes and uniformly map structural nodal displacements to fluid nodes.

[0139] S44: Structural Solution: The structural solver receives the load and solves for the new displacement. .

[0140] S45: Calculate residuals: Calculate force residuals and displacement residual and calculate Norm.

[0141] S46: Convergence Criterion: If and If the iteration converges, the sub-iteration exits.

[0142] S47: Relaxation Update: If convergence fails, update the relaxation factor using the interface quasi-Newton method. ,make Return to S42.

[0143] In this embodiment, on average, 15 to 20 sub-iterations are required to achieve convergence at each time step in the early stage of water entry; in the later stage of water entry, the number of sub-iterations is reduced to 5 to 10.

[0144] III. Analysis of Calculation Results and Summary of Implementation Effects

[0145] To verify the effectiveness of the method proposed in this invention, this embodiment conducted a numerical simulation of the uniform vertical water entry process of a convex wedge with χ=0. The spatiotemporal distribution of the impact pressure, the structural displacement contour map, the stress time history curve, and the evolution process of the free surface were calculated, as shown below. Figures 7 to 10 As shown.

[0146] Figure 7 The spatiotemporal distribution cloud map of the slamming pressure calculated for the example is shown in the figure. As shown, in the initial stage of water entry ( A significant high-pressure zone appears in the central region of the wedge's bottom, with peak pressures exceeding 100 kPa, and this high-pressure zone exhibits a locally concentrated distribution. As the depth of entry into the water increases ( The high-pressure zone gradually diffuses towards the edges of the wedge, and the peak pressure drops below 50 kPa. At that moment, the pressure distribution tends to be uniform, and the peak value further decreases. The pressure pulsation obtained by the prediction method of this invention is smoother, thanks to the Virtual Fluid Method (GFM) introduced in step S2 of this invention, which effectively suppresses spurious velocities at the gas-liquid interface and avoids non-physical pressure oscillations, thereby more realistically reflecting the buffering effect of the air cushion under the structure on the slam load.

[0147] Figure 8 The figure shows the displacement distribution cloud map of the elastic thin-shell structure calculated for this embodiment. It illustrates the deflection distribution of the wedge-shaped base plate at different times. Initially, the maximum displacement occurs in the forward part of the base plate during the initial water entry phase. With the continued action of the impact load, the maximum displacement region concentrates towards the center of the base plate at the peak of the vertical impact force. Subsequently, the structure undergoes elastic vibration, and the displacement gradually decreases. Compared to the rigid structure assumption, this invention, considering elastic deformation, allows the structure to absorb some of the impact energy through deformation, thereby effectively reducing the local peak pressure.

[0148] Figure 9 The stress-time history curves of the elastic thin-shell structure calculated for this example are shown, where P1-5 represent measuring points 1-5. The time history curve for P1 is a dashed line, for P2 it's a dotted line, for P3 it's a solid line, for P4 it's a dotted line with a different spacing, and for P5 it's a double-dotted line. As shown, the stress at each measuring point exhibits a trend of initially rising rapidly, reaching a peak, and then decaying and oscillating. The peak stress at measuring point 1 (the point that first comes into contact with water) appears at... The peak stress at measuring point 3 (midpoint) is approximately 129 MPa; the peak stress at measuring point 3 (midpoint) is slightly later. The peak stress at measuring point 5 (edge ​​point) was 92 MPa; the minimum peak stress at measuring point 5 (edge ​​point) was 81 MPa. It is worth noting that the stress-time history curve showed the lowest peak stress during the slamming stage (…). The stress exhibits a certain degree of high-frequency fluctuation, which is due to the millisecond-level interaction between the air cushion pressure and structural deformation captured by the implicit strong coupling algorithm in step S4; after entering the transition stage, the stress exhibits low-frequency decay oscillation, corresponding to the structure's natural vibration characteristics.

[0149] Figure 10The figure shows a comparison of the free surface evolution process calculated in the example, where (a) to (d) correspond to the free surface morphology at different times. The figure compares the calculation results of the PLIC free surface capture method (left side) and the MULES free surface capture method (right side) of this invention. At time , the free surface lifting height obtained by the two methods is basically the same, and the jet morphology is not significantly different. The above comparison fully demonstrates that this invention, by combining the PLIC-based geometric VOF method in step S1 with the GFM in step S2, achieves high-precision capture of the complex evolution of sub-millimeter-level air cushions and free surfaces (including jet separation and flipping), while the MULES method based on traditional algebraic VOF cannot accurately simulate such fine free surface features.

[0150] Through the above specific embodiments, this invention can achieve high-precision numerical simulation of the gas-liquid-solid three-phase coupling mechanism during cross-medium impact, and is particularly suitable for the analysis of local strong nonlinear impact on thin-shell structures and the evolution of air cushions during water entry. The main technical effects are summarized as follows:

[0151] (1) Accurate pressure field: By combining the compressible two-phase flow model with GFM, the false velocity at the interface is eliminated, the distortion of the slam pressure peak is avoided, and the air cushion buffering effect is reflected more realistically.

[0152] (2) Reliable structural response: The elastic thin shell dynamic model and implicit strong coupling algorithm are adopted to realize the synchronous solution of structural deformation and air cushion pressure, which provides a more reliable basis for structural safety assessment.

[0153] (3) Fine capture of free surface: The geometric VOF method based on PLIC realizes high-resolution simulation of complex evolution forms such as free surface jet separation and flipping, breaking through the numerical diffusion limitation of the traditional algebraic VOF method.

[0154] (4) Strong engineering applicability: The wedge structure, material parameters, water entry velocity and other parameters used in this embodiment are typical and representative. The calculation is stable and has good convergence, indicating that the present invention can be extended to the prediction of slam load and impact resistance design of various cross-medium navigation bodies, seaplane hulls and other structures.

Claims

1. A method for predicting the impact load of an elastic thin-shell considering three-phase coupling of gas, liquid, and solid phases, characterized in that, Includes the following steps: S1: Construct two-phase flow control equations that consider gas-liquid compressibility, and use the geometric fluid volume (VOF) method based on piecewise linear interface reconstruction (PLIC) to establish a high-precision free surface tracking model for compressible two-phase flow, which is used to accurately capture the evolution of gas-liquid interface including sub-millimeter-scale air cushions in cross-medium processes. S2: Introducing the Virtual Fluid Method (GFM) to precisely maintain the sharp discontinuity of density and pressure at the gas-liquid interface, suppressing the spurious velocity introduced by the dynamic pressure term, and obtaining an accurate slamming pressure field; S3: Establish a finite element dynamic model of the elastic thin shell structure to solve the transient large deformation behavior of the thin shell of the aircraft body under impact load; S4: Using a partitioned coupling framework, the compressible two-phase flow model established in steps S1-S2 and the structural dynamics model established in step S3 are solved synchronously and iteratively through an implicit strong coupling algorithm, so as to realize the strong coupling calculation of the air cushion compression process and the elastic thin shell deformation process in the millisecond time domain. The spatiotemporal distribution of impact pressure, structural dynamic response, and evolution characteristics of compressible air cushion of the output elastic thin-shell structure during cross-medium impact process are studied.

2. The method for predicting the impact load of an elastic thin shell considering three-phase coupling of gas, liquid, and solid phases as described in claim 1, characterized in that: The two-phase flow control equations considering the compressibility of water vapor in S1 include the mass conservation equation, momentum conservation equation, total energy conservation equation, and volume fraction transport equation. The gas-liquid interface is reconstructed and captured using the PLIC geometric VOF method, specifically expressed as follows: (1) Mass conservation equation: ； in, For the density of the mixed fluid, For time, It is a velocity vector. Represents the divergence operator; (2) Momentum conservation equation: ； in, For pressure, For viscous stress tensor, It is the gravitational acceleration vector. This is the surface tension volume force term. This represents the gradient operator; (3) Total energy conservation equation: ； in, Total energy per unit mass This is the heat flux density vector; The specified unit mass always satisfies: ； in, Internal energy per unit mass Kinetic energy per unit mass; (4) Volume fraction transport equation: ； in, It represents the volume fraction of the aqueous phase. This refers to the gas phase volume fraction. The interface compression velocity is used to suppress numerical diffusion at the interface and maintain the clarity of the gas-liquid interface; the gas-liquid interface is reconstructed into a piecewise linear interface within the interface unit according to the volume fraction distribution using the PLIC geometric VOF method. The physical properties of the mixed fluid satisfy: ； ； in, The density of the aqueous phase, The density is the gas phase density. For mixed fluid dynamic viscosity, For the dynamic viscosity of the aqueous phase, This refers to the dynamic viscosity of the gas phase. The aqueous and gas phases satisfy their respective equations of state to establish the coupling relationship between pressure, density, and energy.

3. The method for predicting the impact load of an elastic thin shell considering three-phase coupling of gas, liquid, and solid phases as described in claim 1, characterized in that: The PLIC-based geometric VOF method in S1 includes the following sub-steps: S11: Define the volume fraction parameter for each grid cell within the computational domain; define the proportion of the target phase volume in the two-phase flow field to the total volume of the grid cells as the volume fraction. ,in, This indicates that the mesh cell is completely filled with the target phase. This indicates that there is no target phase within this mesh cell. This indicates that the grid cell is an interface cell. By initializing the volume fraction of all grid cells in the computational domain, a discrete volume fraction field is established to characterize the distribution of the two-phase interface. S12: Identify interface units and extract local interface information; based on the volume fraction values ​​of each grid unit, filter those that meet the requirements. The mesh elements are used as interface elements; at the same time, the volume fraction distribution information of the interface elements and their adjacent mesh elements is obtained to provide a data basis for subsequent interface normal vector solving and interface geometry reconstruction. S13: Calculate the interface normal vector within the interface unit; calculate the interface normal vector based on the spatial gradient of the volume fraction field surrounding the interface unit, wherein the interface normal vector can be expressed as: ； in, For the interface normal vector, It is the volume fraction. For volume fraction gradient; The volume fraction gradient can be obtained through central difference, least squares method or weighted reconstruction method to improve the accuracy of interface orientation recognition; S14: Construct a piecewise linear interface within the interface unit; using the PLIC method, approximate the real phase interface within each interface unit using a linear interface, wherein the linear interface satisfies the following expression: ； in, Let be the position vector of any point on the interface. The interface intercept constant; The orientation of the linear interface is determined based on the interface normal vector obtained in step S13, and the intercept constant is then solved by combining the volume fraction constraints within the interface element. This ensures that the geometric volume of the target phase segmented by the linear interface within the current grid cell is consistent with the target phase volume corresponding to the volume fraction of that grid cell. S15: Solve for the interface intercept constant based on the volume fraction constraint; let the total volume of the current interface element be... The target phase volume fraction is Then the linear interface must satisfy the following volume conservation relationship: ； in, Represents the normal vector and intercept constant The target phase volume intercepted within the current mesh cell by the determined linear interface; The intercept constant can be obtained using analytical methods, table lookup methods, or iterative methods. To complete the interface geometry reconstruction; S16: Construct a swept volume on a grid surface based on the velocity field and calculate the phase volume flux; at a given time step... Within this time step, a sweeping region is constructed based on the normal velocity and area of ​​each grid surface, and the transport volume of the target phase through the grid surface is determined. For any grid surface, its sweeping volume can be expressed as: ； in, For the normal velocity of the mesh surface, For the area of ​​the grid surface, For time step; Furthermore, combining the linear interface positions obtained in S14 and S15, the geometric overlap volume between the swept region and the target phase region in the upstream unit is calculated, which is used as the geometric flux of the target phase through the grid surface. S17: Update the volume fraction of each grid cell based on geometric flux; according to the target inflow and outflow volumes of each grid surface, perform a conservation update on the volume fraction of each grid cell, using the following formula: ； in, and They represent the first The volume fraction of each grid cell at the current time and the next time step. For the first Volume of each grid cell This represents the total outflow volume of the target phase within the current time step for this mesh cell. The total volume of the target phase inflow within the current time step for this mesh cell; S18: Apply boundedness correction to the updated volume fraction; to ensure that the volume fraction satisfies the physical constraints, apply boundedness restrictions to the updated volume fraction to ensure that it satisfies: ； When the updated volume fraction of a certain grid cell is less than 0, it is corrected to 0; when the updated volume fraction is greater than 1, it is corrected to 1, in order to avoid non-physical interface diffusion or spurious phase volume during numerical calculation. S19: Repeat the interface reconstruction and volume fraction transport process; in each time step, repeat steps S12 to S18 to update the position and shape of the two-phase interface in real time, thereby realizing continuous capture and evolution calculation of the gas-liquid interface.

4. The method for predicting the impact load of an elastic thin shell considering three-phase coupling of gas, liquid, and solid phases as described in claim 1, characterized in that: The virtual fluid method (GFM) introduced in S2 includes the following sub-steps: S21: Identify the gas-liquid interface location; determine the interface unit based on the distribution results of the volume fraction function in step S1, and combine the piecewise linear interface reconstructed by the PLIC geometric VOF method to obtain the local position, normal direction and medium properties on both sides of the interface, as the geometric basis for constructing the virtual fluid state; GFM needs to first clarify the interface location, and then process the fluids on both sides in the interface neighborhood according to the single medium format. S22: Extract the actual fluid state on both sides of the interface; In the interface adjacent cell, the real fluid state variables of the liquid phase side and the gas phase side are extracted respectively. The real fluid state variables include at least density, pressure, velocity, total energy and sound speed determined by the equation of state, which are used for subsequent calculation of local interface interactions. S23: Establish a local Riemann problem along the interface normal; using the interface unit normal as the solution direction, project the actual fluid velocities on both sides obtained in step S22 onto the normal direction to construct a local one-dimensional Riemann problem along the interface normal direction to characterize the propagation and interaction process of pressure waves, velocity waves and contact discontinuities at the gas-liquid interface. S24: Solve for the intermediate state of the interface; Solve the local Riemann problem established in step S23 to obtain the intermediate interface pressure and the intermediate interface normal velocity; where the intermediate interface pressure is used to characterize the pressure equilibrium state after the interaction of fluids on both sides of the interface, and the intermediate interface normal velocity is used to characterize the propagation velocity or interface velocity along the normal direction at the interface; For compressible multi-medium problems, linearized interface interaction methods or approximate Riemann solvers are often used in the literature to calculate the interface state; In the case of surface tension, a pressure transition term caused by curvature can also be added to the interface pressure condition; S25: Construct virtual fluid states on both sides of the interface; based on the intermediate interface pressure and intermediate interface normal velocity obtained in step S24, construct corresponding virtual fluid states on both sides of the interface respectively; the virtual fluid states include at least virtual pressure, virtual normal velocity, virtual density, and virtual total energy; wherein, the virtual pressure and virtual normal velocity are determined by the intermediate interface state, and the virtual density can be obtained based on the corresponding medium state equation and combined with isentropic extrapolation or characteristic extrapolation, thereby avoiding non-physical mixing and pressure oscillation caused by direct interpolation across the interface; S26: Complete the single-phase discrete template near the interface; using the virtual fluid state constructed in step S25, complete the single-phase calculation template of the unit near the interface in the differential, finite volume or finite reconstruction process, so that the gas phase region and the liquid phase region can be solved independently according to the control equation of their respective medium during numerical discretization, without having to directly call the heterogeneous unit data across the real gas-liquid interface during the discretization process. S27: Correct the numerical flux near the interface; Based on the single-phase template completed in step S26, calculate the mass flux, momentum flux and energy flux of the unit near the interface, and substitute the numerical flux into the two-phase flow control equation for discrete update, so as to reduce the non-physical oscillation, interface diffusion and numerical overheating caused by abrupt changes in physical properties near the interface. S28: Update the interface position and repeat the virtual fluid construction process; after completing the flow field update at the current time step, continue to update the volume fraction field and the gas-liquid interface position, and repeat steps S21 to S27 for the new interface neighborhood, thereby realizing continuous simulation of interface propagation, wave system reflection and transmission and strong discontinuous interaction in the transient compressible water-gas two-phase flow process.

5. The method for predicting the impact load of an elastic thin shell considering three-phase coupling of gas, liquid, and solid phases as described in claim 1, characterized in that: The finite element dynamic equations for the elastic thin-shell structure in step S3 are as follows: ； in, The overall mass matrix of the elastic thin-shell structure. For the overall damping matrix, For the overall stiffness matrix, Let be the nodal displacement vector. For the node velocity vector, For the nodal acceleration vector, The vector of external loads acting on the structure. This is the fluid-structure interaction load vector acting on the surface of the thin-shell structure.

6. The method for predicting the impact load of an elastic thin shell considering three-phase coupling of gas, liquid, and solid phases as described in claim 1, characterized in that: The implicitly strongly coupled algorithm in S4 includes the following sub-steps: S41: At the beginning of each physical time step, the structural displacement converged in the previous time step is passed to the fluid solver as the initial boundary condition; S42: The fluid solver solves a compressible two-phase flow model based on the current structural boundary, calculating the fluid dynamic load vectors acting on the fluid-structure interaction interface. ; S43: The radial basis function interpolation method is used to interpolate the hydrodynamic load. And conservatively map from the fluid mesh to the structural mesh, while simultaneously transferring structural displacements. A consistent mapping from the structural mesh to the fluid mesh is used to maintain energy balance at the coupling interface; S44: The structural solver receives the mapped load, solves the finite element dynamic equations of the elastic thin-shell structure, and obtains the new structural displacements. ; S45: Calculate the... Force residual vector after the second iteration and displacement residual vector And calculate their discrete values ​​respectively. Norm: ； S46: Use the relative convergence criterion to determine whether the iteration has converged: ； in For the force convergence threshold, This is the displacement convergence threshold; if the above conditions are met, the sub-iteration is exited and the next time step is initiated. S47: If convergence fails, update the structural boundary displacements using the interface quasi-Newtonian relaxation method: ； in The adaptive relaxation factor is then used, and the process returns to step S42 to continue iterating until the convergence criterion is met or the preset maximum number of iterations is reached.

7. The method for predicting the impact load of an elastic thin shell considering three-phase coupling of gas, liquid, and solid phases as described in claim 6, characterized in that: The force convergence threshold and displacement convergence threshold The average value is .

8. The method for predicting the impact load of an elastic thin shell considering three-phase coupling of gas, liquid, and solid phases as described in claim 6, characterized in that: Conservative mapping ensures the conservation of total force through the coupling interface, while uniform mapping ensures the continuity of the displacement field.

9. A system for predicting impact loads on a gas-liquid-solid three-phase coupled elastic thin-shell shell, characterized in that: This system is used to execute any one of the methods for predicting the impact load of a gas-liquid-solid three-phase coupled elastic thin shell according to claims 1-8, including: The compressed two-phase flow calculation module is used to execute steps S1-S2; The elastic thin-shell structure calculation module is used to execute the steps in S3. An implicitly strongly coupled iterative module is used to execute the steps of S4; Within each time step, the compressible two-phase flow calculation module and the elastic thin-shell structure calculation module are driven to perform synchronous iterative solutions. It also outputs the spatiotemporal distribution of impact pressure, structural dynamic response, and the evolution characteristics of compressible air cushions.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that: When executed by the processor, the program implements the method for predicting the impact load of an elastic thin shell considering three-phase coupling of gas, liquid, and solid phases, as described in any one of claims 1 to 8.

Images