Numerical prediction method for three-phase flow of hypersonic vehicle entering water with strong coupling of structure

Abstract

The application discloses a hypersonic vehicle water-entry three-phase flow-structure strong coupling numerical prediction method. The method is based on the physical and causal relationship of 'pressure wave propagation-super-cavitation initial generation-pressure reloading-structure displacement feedback', takes the pressure prediction result as the driving source of the super-cavitation generation and extinction criterion and the fluid-solid interface load updating, realizes the synchronous prediction of pressure evolution, super-cavitation behavior and structure response in the same calculation framework, thereby avoiding the problems of pressure peak distortion, super-cavitation starting time deviation and structure load prediction lag under the conditions of weak coupling and incompressible assumption, and improving the prediction reliability of the pressure peak and its time sequence evolution under high-pressure impact working conditions.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation technology for the entry process of hypersonic moving bodies into water, specifically a numerical prediction method for three-phase flow-structure strong coupling of hypersonic moving bodies entering water. Background Technology

[0002] During the entry of a hypersonic vehicle into water, the vehicle crosses the free liquid surface at high speed from the air medium, accompanied by complex multi-physics coupling phenomena such as free liquid surface rupture, severe deformation of multiphase interfaces, cavitation generation and collapse, pressure wave propagation, water hammer impact loads, and structural dynamic response. Due to the characteristics of this process—high speed, short duration, large pressure peaks, rapid abrupt changes in interface morphology, and strong transient structural loads—real-world experiments require the simultaneous measurement of key data such as pressure field, cavitation evolution, fluid-solid interface loads, and structural deformation response on a microsecond timescale. The related measurement conditions are extremely demanding, resulting in high experimental costs. Therefore, establishing a numerical simulation method for the entry process of a hypersonic vehicle into water is of significant necessity and engineering importance. Data obtained through simulation, including pressure field, velocity field, cavitation morphology, vapor phase volume fraction, interface loads, and structural stress, strain, and displacement, can be used to analyze the formation mechanism of water impact loads, determine the location of cavitation generation and collapse, assess structural strength safety and deformation risk, and provide quantitative basis for vehicle shape design, material selection, protective structure design, and water entry attitude optimization.

[0003] Existing numerical simulation methods for hypersonic water ingress primarily control the velocity field or interface morphology, with the pressure field typically used as a result for subsequent analysis. These methods have not yet fully integrated the unified control of cavitation evolution criteria, structural load updates, and fluid-structure interaction responses. Therefore, current technologies struggle to achieve simultaneous, high-precision prediction of hypersonic water ingress impact loads, cavitation evolution, and structural responses. Summary of the Invention

[0004] This invention addresses the challenges of pressure transient changes dominating multiphysics evolution, significant liquid compressibility impacts, and difficulties in accurately and synchronously predicting structural responses during hypersonic body entry into water. It proposes a three-phase flow-cavitation-fluid-structure interaction numerical prediction method, specifically a three-phase flow-structure interaction numerical prediction method for hypersonic body entry into water, with pressure wave propagation prediction as the dominant physical mechanism. This method uses pressure wave propagation caused by liquid compressibility as the core physical thread and pressure prediction results as the unified control variable. By constructing a closed-loop calculation process involving pressure-cavitation evolution-structural response-mesh update, it achieves stable prediction of key physical quantities during high-speed water entry. This method is applicable not only to rigid body entry problems but also to fluid-structure interaction problems involving hyperelastic and viscoelastic bodies.

[0005] To achieve the above objectives, the present invention adopts the following technical solution: A numerical prediction method for strongly coupled three-phase flow-structure interaction of a hypersonic moving body entering water is designed, characterized in that the method includes at least the following steps: Step 1: Define the computational domain and boundary conditions; Step 2: Establish a three-layer overset mesh system; Step 3: Define the air-water-steam three-phase volume fraction model; Step 4: Define the three-phase volume fraction transport equation; Step 5: Establish a liquid-phase compressibility model based on the Tait equation of state; Step 6: Construct a cavitation evolution model based on local pressure criteria; Step 7: Set up the structural domains and structural material constitutive models; Step 8: Set the three-phase flow control equations; Step 9: Establish a fluid-structure interaction model; Step 10: Based on the initial entry velocity, initial entry angle and initial attitude information of the hypersonic motion body, and combined with Steps 1 to 9, perform data prediction through a two-way strongly coupled pressure-displacement iteration.

[0006] Compared with existing technologies, the beneficial effects of this invention are as follows: The numerical prediction method for three-phase flow-structure strong coupling of hypersonic moving bodies entering water provided by this invention differs from existing numerical methods that focus on velocity field or interface evolution. The method of this invention: (1) Starting from the physical causal relationship of "pressure wave propagation - cavitation generation - pressure reloading - structural displacement feedback", the pressure prediction result is used as the control quantity for cavitation evolution criterion and fluid-solid interface load update. The pressure evolution, cavitation behavior and structural response are predicted synchronously within the same calculation framework, thereby avoiding the problems of pressure peak distortion, cavitation initiation time offset and structural load prediction lag that are common under weak coupling and incompressible conditions. This improves the prediction reliability of pressure peak and its temporal evolution under high pressure impact conditions.

[0007] (2) The cavitation evolution criterion is directly coupled with the pressure prediction results, realizing continuous time series prediction of the cavitation generation, development and collapse process, reducing the problem of lag or position shift of cavitation in traditional models.

[0008] (3) The pressure-displacement bidirectional strong coupling strategy is adopted to synchronously update the interface pressure, structural displacement and foreground body-fitting follower mesh configuration within a single time step. This not only effectively reduces the common load lag and structural response distortion problems in coupling methods, but is also applicable to rigid bodies, flexible structures, deformable structures and large deformation fluid-structure coupling problems. It can more realistically reflect the bidirectional feedback mechanism of "impact load - structural deformation - flow field reconstruction - load redistribution".

[0009] (4) The three-layer overset grid system can maintain grid quality and variable transfer stability under conditions of hypersonic speed, large displacement, cross gas-liquid interface and structural deformation, and realize continuous numerical calculation of complex water entry process.

[0010] (5) A unified multiphysics prediction framework based on "pressure wave propagation - phase change cavitation - structural response" was constructed. Liquid compressibility is no longer used as an auxiliary correction factor, but as the core control factor that determines the formation of pressure peak, pressure wave propagation speed, cavitation evolution sequence and structural load response. Therefore, this method can not only be used for the assessment of water impact load and structural safety analysis of hypersonic moving bodies, but also for the impact response analysis and optimization design of cross-medium high-speed equipment, protective structures and energy-absorbing structures, and has clear engineering application value. Ignoring fluid compressibility will lead to the prediction distortion of pressure wave propagation, impact load peak and cavitation evolution process, while large structural deformation will significantly change the load transfer path and energy dissipation process of the fluid-solid interface.

[0011] The prediction method proposed in this invention is not a simple superposition of multiple physical models. Instead, it introduces a liquid compressibility model based on the Tait equation of state to dynamically calculate the formation of pressure peaks, local sound velocity changes, and pressure wave propagation characteristics. The pressure prediction results are then uniformly used as the control variables for cavitation evolution criteria, fluid-structure interface load updates, and structural response solutions, forming a unified "pressure-phase change-structural response" prediction framework. Under this mechanism, the evolution of the three-phase interface, changes in vapor phase volume fraction, structural displacement updates, and follow-up mesh updates are all advanced collaboratively within the same timescale, thereby achieving stable, continuous, and high-precision prediction of key physical quantities in hypersonic water entry processes. Furthermore, the bidirectional fluid-structure coupling framework in this method is applicable not only to rigid body motion prediction but also to the impact coupling analysis of flexible structures and structures with large deformations, thus expanding the engineering applicability of this invention. Therefore, this invention adopts a bidirectional strongly coupled modeling approach that simultaneously considers the high-speed compressibility of liquids and the large deformation response of structures to improve the completeness and consistency of the prediction of key physical quantities in hypersonic water entry processes. Attached Figure Description

[0012] To better understand the technical solution of the present invention, the present invention will be described in conjunction with the accompanying drawings. The drawings are used to illustrate the basic process and key modules of the present invention and do not constitute a limitation on the scope of protection of the present invention.

[0013] Figure 1 is a flowchart of one embodiment of the numerical prediction method for three-phase flow-structure strong coupling of hypersonic moving body entering water according to the present invention.

[0014] Figure 2 is a schematic diagram of the structure of a three-layer overset overlapping grid system in an embodiment of the numerical prediction method for three-phase flow-structure strong coupling of hypersonic moving body entering water according to the present invention; wherein: 1 is the extrapolation layer, 2 is the foreground body-fitting follower grid, 3 is the receiving layer, 4 is the interpolation layer, and 5 is the background grid.

[0015] Figure 3 is a schematic diagram of the evolution of cavitation structure over time in one embodiment of the numerical prediction method for three-phase flow-structure strong coupling of hypersonic moving body entering water according to the present invention, which is used to illustrate the ability of the method of the present invention to capture changes in cavitation morphology, volume fraction and morphological changes; wherein, the white block shape is indicated as a hypersonic moving body.

[0016] Figure 4 is a schematic diagram of the computational domain and boundary conditions of one embodiment.

[0017] Figure 5 is a schematic diagram of the fluid pressure distribution during the bidirectional strong coupling iteration of pressure-displacement at the fluid-solid interface in an embodiment of the numerical prediction method of three-phase flow-structure strong coupling for hypersonic moving body entering water according to the present invention (the high pressure zone is an arc shape extending towards the rear of the moving body, and the pressure decreases sequentially from the front end of the head to the rear).

[0018] Figure 6 Here are two diagrams showing the evolution of actual flow field characteristics (grayscale image) according to one embodiment and the evolution of flow field characteristics obtained using the method of this invention; wherein, Figure 6 In group (a), the grayscale image on the left is the actual flow field feature at 0.05 ms, and the image on the right is the flow field feature at 0.05 ms obtained using the method of this invention. Figure 6 In group (b), the grayscale image on the left is the actual flow field characteristic at 0.25ms, and the image on the right is the flow field characteristic at 0.25ms obtained using the method of this invention. Figure 6 In group (c), the grayscale image on the left is the actual flow field characteristic at 0.45ms, and the image on the right is the flow field characteristic at 0.45ms obtained using the method of this invention. Figure 6 The grayscale image on the left in group (d) is the actual flow field characteristics at 0.9 ms, and the image on the right is the flow field characteristics at 0.9 ms obtained using the method of this invention.

[0019] Figure 7 Here are two diagrams: one showing the actual shock wave and pressure before water entry, and the other showing the pressure cloud map in the air obtained through simulation calculation using the method of this invention. Figure 7 (a) is the actual diagram of shock wave and pressure before entering the water, and (b) is the pressure cloud diagram in the air medium obtained by the method of the present invention. Detailed Implementation

[0020] The technical solution of the present invention will be described below with reference to specific embodiments.

[0021] The present invention provides a numerical prediction method for three-phase flow-structure strong coupling of hypersonic moving bodies entering water, which includes at least the following steps: Step 1: Define the computational domain and boundary conditions; Based on the solid structure and dimensions of a hypersonic vehicle with a streamlined head and axisymmetric shape, a structural domain model of the vehicle is constructed. In this invention, the structural domain, structural response, and structural material parameters all refer to the hypersonic vehicle itself. As the vehicle enters the water from the air, a computational domain including both the air and water domains is established, with the computational domain's scale larger than that of the vehicle. Based on actual monitoring of the entire process of the vehicle entering the water from the air, the longitudinal range of the water domain is generally taken as 20 to 30 times the length of the vehicle. The boundary conditions of the computational domain are as follows: (a) Solid wall (i.e., solid surface) no-slip boundary: , This represents the volume-average velocity vector of a three-phase mixture of air, water, and steam. Represents the local interface normal vector; (b) Free surface pressure boundary: , This represents the local pressure in the flow field; Indicates atmospheric pressure in the air domain; (c) The velocity derivative is zero (open boundary): , This represents the volume-average velocity vector of a three-phase mixture of air, water, and steam. Represents the local interface normal vector; (d) Pressure derivative is zero (commonly used at lateral boundaries): , Indicates the local pressure in the flow field. This is the unit normal vector outside the local boundary; As the moving body moves from the air domain into the water domain, the air phase, water phase, and vapor phase within the computational domain form a fluid domain, while the moving body is the structural domain. The fluid domain uses the following boundary conditions: (1) Free surface: Atmospheric pressure is applied And adopt local pressure of the flow field Normal velocity gradient condition Meanwhile, the gas-liquid interface shape is tracked using a three-phase volume fraction interface capture method; (2) Side walls and bottom solid walls: Apply no-slip conditions ; (3) Top air region: Apply pressure outlet conditions to allow pressure waves to propagate freely outward, satisfying the following conditions. , Indicates the local pressure in the flow field. Let be the unit normal vector outside the local boundary.

[0022] The boundary conditions of the structural domain are as follows: (1) The outer wall of the structure is subjected to motion or deformation response calculation based on the selected structural domain model; when a deformable structural model is used, the structural response is solved according to the structural dynamics equation; when a rigid body type is used, the overall motion response is updated according to the equivalent load of the fluid-structure interface. (2) If there is a fixed region in the structure, then apply displacement constraints to that region. ; (3) The fluid load on the fluid-solid interface is transferred to the structural domain in real time by the interface traction stress calculated by the fluid domain at the interface through a mapping algorithm; in the implementation method where the normal pressure is the main force, the fluid load can be approximately represented by the transient pressure field.

[0023] Step 2: Establish a three-layer overset mesh system; In step one, a background grid is set up within the computational domain to describe large-scale water flow, cavitation cloud structure, pressure wave propagation, and surface rupture behavior. The moving body is located within the background grid, and an overset moving grid is constructed around its surface. This overset moving grid is body-fitted according to the geometry of the moving body, with its grid boundaries consistent with the outer surface of the moving body. This grid is used to precisely describe the velocity gradient, pressure gradient, and three-phase interface changes in the near-wall region of the moving body.

[0024] The foreground body-fitted follower mesh consists of multiple element nodes, each representing a discrete spatial position within the mesh. The mesh is formed by connecting adjacent discrete nodes. It is generated initially based on the geometric model of the moving body and establishes a correspondence with the surface nodes of the structural domain through node mapping, enabling the foreground body-fitted follower mesh to update synchronously with the movement of the moving body structure.

[0025] In the numerical calculation process, after the structural dynamics solution yields the displacements of the moving body's structural nodes, these displacements are synchronously transferred to the corresponding nodes of the foreground body-fitting follower mesh through node mapping relationships. This updates the node coordinates of the foreground body-fitting follower mesh, enabling the body-fitting mesh to move synchronously with the moving body. The foreground body-fitting follower mesh consists of multiple mesh elements and their discrete nodes. Each node represents a spatial discrete point of the body-fitting mesh element, and the node number is... Its spatial coordinates are , is used to describe the position of the body mesh in the computational domain.

[0026] A three-layer overset mesh system is constructed between the background mesh domain and the foreground body-fitting mesh. This three-layer overset mesh system includes an extrapolation layer, an interpolation layer, and a receiving layer, used to achieve mesh crossing, variable transfer, and interface mapping coordination for the moving body during large displacement processes. The specific construction method is as follows: Within the overlapping region of the background mesh domain and the foreground body-fitted follower mesh, the layer properties of each mesh element are determined by their geometric positional relationship: When the geometric center of a mesh element in the background mesh domain falls within the envelope region of the foreground body-fitted follower mesh, this element is marked as an extrapolation layer element, i.e., a redundant element, and does not participate in the solution of fluid domain field variables and interpolation calculations. When the geometric center of a mesh element in the background mesh domain is located outside the envelope region of the foreground body-fitted follower mesh but the distance from the mesh element in the foreground body-fitted follower mesh is not greater than a set value, this element participates in the interpolation calculation as a candidate donor element, and the interpolation layer is composed of donor elements. Nodes or elements in the foreground body-fitted follower mesh that receive data from the background mesh domain are defined as the receiving layer.

[0027] The interpolation layer transmits field variables between the background grid domain and the foreground body-fitted servo grid using a distance-weighted interpolation method. These field variables include pressure, velocity components, and volume fractions of air, water, and vapor phases. This ensures continuous transmission and consistent updating of physical quantities between the background and foreground body-fitted servo grids, adapting to hypersonic, large-displacement water entry conditions. The receiving layer receives the interpolation results output by the interpolation layer and updates the fluid domain field variables on the foreground body-fitted servo grid.

[0028] To ensure interpolation stability, a transition interpolation region is set outside the foreground body-fitted follower mesh. Neighboring background mesh elements are searched based on the node positions of the foreground body-fitted follower mesh. Discrete nodes or receiving elements in the foreground body-fitted follower mesh receive the field variables interpolated by the background mesh donation elements, thus achieving continuous transfer of physical quantities such as pressure, velocity, and volume fraction from the overlapping region to the foreground body-fitted follower mesh. Several background mesh elements closest to the foreground body-fitted follower mesh nodes are selected as donation elements for calculating the interpolation weights of the field variables. This completes the continuous mapping of field quantities such as pressure, velocity, and volume fraction between the two mesh sets. The three-layer Overset overlapping mesh system can maintain mesh stability under hypersonic large-displacement water entry conditions, avoiding the distortion of traditional moving meshes.

[0029] The field variables obtained by the receiver layer from the foreground body-fitted mesh are obtained by the interpolation layer using the distance-weighted interpolation method. Let the field variable of a single mesh node in the foreground body-fitted mesh be... ,Should From the background grid domain, the distance to the grid node is not greater than a set threshold. The interpolation is obtained by weighted averaging of the field variables of each donation unit. The interpolation formula is:

[0030] in, For the first The field variable values ​​corresponding to each donation unit; The first kernel, determined based on the inverse distance ratio or the RBF kernel function, is... Weight of each donation unit:

[0031] For the first The distance between each donation unit and the foreground-fitted moving mesh node. For the first The distance between each donation unit and the foreground-fitted follow-mesh node.

[0032] This interpolation method ensures stable and continuous variable transfer of the foreground body-fitted moving mesh under hypersonic large displacement conditions. Among these, the field variables to be interpolated... The field variables transmitted within the overlapping grid region include, but are not limited to, pressure, velocity components, and volume fractions of each phase, to ensure the consistency of physical quantity distribution between the background grid domain and the foreground body-fitted follower grid.

[0033] Step 3: Define the air-water-steam three-phase volume fraction model; To describe the three-phase mixed flow of air, water, and steam within the aforementioned overlapping grid system, a three-phase volume fraction model is defined based on the grid structure. To describe the inlet interface rupture, cavitation generation, and steam phase evolution, this invention employs a three-phase volume fraction interface capture model. The volume fractions of the air, water, and steam phases are defined as follows: , , And satisfy the three-phase volume fraction constraint relationship: .

[0034] The density and viscosity of a mixture are defined as a volume-fraction-weighted average:

[0035]

[0036] in, , , These represent the densities of the air phase, water phase, and vapor phase, respectively. , , These represent the dynamic viscosity of the air phase, water phase, and vapor phase, respectively. and These represent the equivalent density and equivalent viscosity of the three-phase mixture, respectively. The mixing density serves as a key physical quantity input for subsequent pressure field solutions and the liquid-phase compressibility model.

[0037] It should be noted that in the method of this invention, the air phase density, vapor phase density, three-phase dynamic viscosity, Tait equation of state parameters, saturated vapor pressure, empirical coefficients of the cavitation model, vapor core volume fraction, vapor core equivalent radius, structural material parameters, damping parameters, boundary conditions, and initial conditions are all pre-given known input quantities. In addition to the above-mentioned known input quantities, the pressure field, velocity field, liquid phase density field, mixed density field, three-phase volume fraction distribution, fluid-solid interface equivalent external load, structural displacement, structural velocity, structural acceleration, and foreground body-fitted follower mesh node coordinates are all calculated step by step by the control equations and coupled iterative process, without the need to introduce new unknown test quantities during the calculation process.

[0038] Step 4: Define the three-phase volume fraction transport equation; The PLIC (Piecewise Linear Interface Construction) interface reconstruction algorithm is employed to ensure the geometric accuracy of the interface shape under high-speed impact. To describe the motion of the air-water-steam three-phase interface, this invention uses a three-phase volume fraction transport equation, as follows:

[0039] in, The corresponding markings for the air phase, water phase, and vapor phase; For the first Phase volume fraction; For gradient operators; This represents the volume-average velocity vector of a three-phase mixture of air, water, and steam. As a phase change source term, evaporation and condensation occur only between the liquid and vapor phases, satisfying the following conditions: .

[0040]

[0041] Let the phase change mass source term from liquid to vapor phase be denoted as ,but:

[0042] This three-phase volume fraction model can maintain interface sharpness under high-speed impact and strong interfacial stretching conditions, and supports complex morphological changes (such as interface rupture, cavitation region splitting and merging). Compared with the traditional two-phase model, it can significantly improve the ability to describe the transient multiphase flow structure in water.

[0043] Step 5: Establish a liquid-phase compressibility model based on the Tait equation of state; During high-speed water ingress, the liquid phase undergoes intense compression within a very short time, generating high-amplitude water hammer pressure. Its density increases considerably with pressure changes. To accurately characterize the compressibility response of the liquid phase under high-pressure impact, this invention employs the Tait equation of state to establish a functional relationship between liquid phase pressure and density.

[0044] The liquid phase density and pressure satisfy a unified Tail equation of state relationship:

[0045] In the formula, For the local pressure of the flow field, This represents the liquid phase density at the corresponding pressure. and These represent the liquid phase density and pressure under reference conditions, respectively. This is an empirical constant related to the bulk modulus of liquids. This is the exponential coefficient. All calculations involving the relationship between liquid phase density and pressure in this paper use the above unified expression. Among them, the empirical constants in the Tait equation of state... The preferred value range is (2.5~3.5)×10 8 Pa, exponential coefficient The preferred value range is 6 to 8.

[0046] Furthermore, the local sound velocity can be obtained from the Tait equation of state. expression:

[0047] in, This refers to the local sound velocity in the liquid phase. This refers to the local pressure in the flow field. The density of the liquid phase is calculated using the Tait equation of state. Reference liquid phase density; These are the empirical constants in the Tait equations of state; The above equation represents the pressure response to density changes under weakly compressible conditions in the liquid phase, obtained by differentiating the Tait equation of state with respect to density. This local sound velocity is used to describe the propagation speed of pressure waves in the liquid and is subsequently used for time step control based on CFL conditions.

[0048] Local sound speed This clearly demonstrates the compressibility behavior of a liquid under pressure, and its changes will directly affect the time step limit based on the CFL number:

[0049] in, The volume-average velocity vector of the three-phase mixture; This refers to the local sound velocity in the liquid phase. This refers to the time step for numerical computation. The CFL number is the characteristic length or equivalent element size of the smallest grid cell (excluding redundant cells) in the current three-layer Overset grid system; it is used to ensure the numerical stability of pressure wave propagation and velocity field calculations. To ensure the stability of pressure wave propagation calculations during hypersonic water entry, the CFL number is preferably controlled within the range of 0.2 to 0.6.

[0050] Step 6: Construct a cavitation evolution model based on local pressure criteria; Hypersonic entry into water is accompanied by significant pressure fluctuations, and the local pressure of the liquid may fall below the saturation pressure, inducing cavitation. Therefore, this invention employs a cavitation evolution model based on a local pressure criterion to describe the vapor phase generation and collapse behavior caused by local pressure reduction. When the local pressure of the flow field... Not higher than saturated vapor pressure When the liquid phase evaporates and generates vapor, the phase change source term manifests as an evaporation source term; when the local pressure of the flow field... Higher than saturated vapor pressure When the steam condenses and transforms into a liquid phase, the phase change source term manifests as a condensation source term.

[0051] vapor phase volume fraction In the evolution equation, the phase transition source term can be expressed as:

[0052] in, For evaporation source term, For condensation source terms, they are defined as follows: when hour, ; ; when hour, ; ; In the formula, , This is an empirical coefficient; Vapor core volume fraction; The equivalent radius of the vapor core; , These are the densities of the liquid phase and the vapor phase, respectively.

[0053] This cavitation evolution model based on local pressure criteria can respond to sudden pressure drops within microseconds. By updating the steam volume fraction, it simulates the evaporation and condensation processes, enabling numerical calculations of cavitation generation, development, and collapse behavior under high-pressure changing conditions.

[0054] Step 7: Set up the structural domains and structural material constitutive models; Based on the geometry, structural dimensions, material type, and service response characteristics of the hypersonic moving body, a structural domain model is established, and corresponding structural material constitutive models, mass parameters, damping parameters, initial states, and boundary constraints are set to describe the dynamic response behavior of the moving body under high-pressure impact loads.

[0055] The structural domain model can be used not only for rigid body type motion response prediction, but also for deformable structural model for structural dynamic response analysis; when a deformable structural model is used, the structural material constitutive model can be selected as a viscoelastic constitutive model or a hyperelastic constitutive model according to the material type and deformation characteristics.

[0056] in: (1) When structural materials exhibit stress-strain rate-related, time-related deformation, hysteresis response or energy dissipation characteristics under impact load, a viscoelastic constitutive model can be used to describe them. (2) When the structural material is under large deformation conditions, the hyperelastic constitutive model can be used to describe it in order to characterize the recoverable deformation characteristics of the material in the large strain range; Furthermore, the structural domain can be set as a rigid body, viscoelastic structure, or hyperelastic structure according to the analysis requirements; correspondingly, the structural response variables include structural displacement, structural velocity, structural acceleration, structural strain, and structural stress.

[0057] By setting the structural domain type and material constitutive relation, the method of the present invention is not only applicable to the problem of rigid body entering water, but also applicable to the prediction of fluid-structure interaction response of hyperelastic and viscoelastic structures under hypersonic water impact conditions, thereby improving the applicability of the method of the present invention in complex material systems and large deformation conditions.

[0058] Step 8: Set the three-phase flow control equations; The three-phase flow control equations specifically include the mixing continuity equation, the momentum equation, and the turbulent viscosity model.

[0059] Mixed continuity equations:

[0060] Momentum equation:

[0061] in: is the density of a three-phase mixture, used to characterize the equivalent density of the three phases of air, water and steam within a unit volume; For time; The term represents the change in the density of the mixture over time; For gradient operators; The volume-average velocity vector of the three-phase mixture; This represents the rate of change of the momentum of a mixture with respect to time. This refers to the local pressure in the flow field. The molecular dynamic viscosity of the three-phase mixture is obtained by weighting the viscosity of each phase by its volume fraction. Turbulent viscosity, used to describe the enhancing effect of turbulent fluctuations on momentum transport, is calculated using a turbulence model; turbulent viscosity is calculated using the k–ε model or an improved model. It represents turbulent kinetic energy and characterizes the intensity of fluid velocity fluctuations; It represents the turbulent kinetic energy dissipation rate, used to describe the rate at which turbulent kinetic energy is dissipated to smaller scales; By simultaneously solving the mixing continuity equation and the momentum equation, the transient pressure and velocity fields within the fluid domain can be obtained. The interfacial traction force distribution at the fluid-solid interface will serve as the load input for the fluid's action on the structure. In the implementation method where normal pressure is dominant, the interfacial traction can be approximately determined by the interfacial pressure distribution. The mixture density is dynamically updated by combining the liquid-phase compressibility model described in step five with the three-phase volume fraction relationship.

[0062] The multiphase volume fraction equation is used to update the three-phase distribution; the compressible equation of state is used to correct for local density; the momentum equation is used to update the velocity field; turbulent viscosity is provided by the turbulence model; the above equations are solved simultaneously in each time step, forming a coupled solution system of three-phase flow-compressible-turbulent flow. The turbulence model provides efficient turbulent viscosity calculation, ensuring stable solutions in high gradient regions.

[0063] Step 9: Establish a fluid-structure interaction model; To reflect the dynamic response of a structure under transient impact loads during hypersonic water entry and its reverse influence on the flow field, this invention establishes a two-way strong coupling model of the fluid-structure interface. This model is not only applicable to the prediction of the water entry response of a moving body as a rigid body, but also applicable to fluid-structure coupling problems of flexible structures, deformable structures and under large deformation conditions described by viscoelastic or hyperelastic constitutive models. It is used to describe the interface position changes, local flow field reconstruction and load redistribution processes caused by structural deformation.

[0064] The fluid-structure interface achieves mechanical and kinematic coordination through the following conditions: (1) Dynamic consistency: ; (2) Kinematic consistency: ; (3) Displacement continuity: ; in: , These are the fluid-side stress tensor and the structural-side stress tensor, respectively. The interface traction vector represents the contact force per unit interface area transmitted from the fluid to the structure or fed back from the structure to the fluid. , These represent the velocities at the fluid-solid interface on the fluid side and the structural side, respectively. , These represent the positions of the fluid side and the structure side at the fluid-structure interface, respectively. When the fluid-structure interface satisfies the coupling condition, the traction, velocity, and position on both sides of the interface are continuous, i.e., they satisfy... , , .

[0065] The external loads borne by the structure are based on the fluid pressure or interface stress at the fluid-solid interface. After interface integration and node mapping, the equivalent external load vectors on the discrete nodes of the structure are formed.

[0066] The interfacial traction stress calculated at the fluid-solid interface in the fluid domain, after interface integration and nodal mapping, is converted into an equivalent external load vector on the structural surface. Its general expression is:

[0067] in, Let T be the shape function matrix corresponding to the discrete nodes of the structure, and let T denote the transpose operation. Let be the interfacial drag vector of the fluid acting on the fluid-solid interface. This represents the interface where the fluid domain and the structural domain come into contact and loads are transferred. This represents the area integral element on the interface; in this invention, the fluid-solid interface is specifically the boundary surface of the hypersonic moving body that is in direct contact with the surrounding three-phase fluid of air, water or steam.

[0068] satisfy:

[0069] In the formula, For fluid-side stress tensor, The external normal unit vector of the fluid-solid interface.

[0070] Since this invention addresses the high-speed water impact condition of a hypersonic moving body, during the initial stages of contact between the moving body and the liquid surface and water entry, the fluid-solid interface experiences significant normal compressive impact. The peak interface pressure is much greater than the tangential viscous shear, and the external load on the structure is mainly determined by the impact pressure transmitted by the fluid along the interface normal. Therefore, the interface traction vector is represented by the normal action of the interface pressure. Correspondingly, the equivalent external load vector on the discrete nodes of the structure can be further expressed as:

[0071] in, The transient pressure at the fluid-solid interface. This is the unit vector normal to the fluid-structure interface. The equivalent external load serves as the external force input term in the structural dynamics equations, driving the structural response calculation.

[0072] The structural dynamics equations are:

[0073] in, The structural mass matrix, Here is the structural damping matrix. Let be the structural node displacement vector. The internal force vector of the structure is determined by the constitutive relation of the structural material, the geometric nonlinear relation, and the discrete form of the element. The equivalent external load vector of the structure is formed by the interface integration and nodal mapping of the traction stress at the fluid-structure interface. This represents the velocity vector of the structural nodes. This represents the acceleration vector of the structural nodes.

[0074] Specifically, the internal force vector of the structure The formation process is as follows: First, the strain field is obtained from the structural nodal displacement field through the strain-displacement relationship. Then, the stress tensor is calculated from the material constitutive relation. Subsequently, the stress tensor is integrated through the element virtual work equivalent integration and overall assembly to form the equivalent internal force vector on the discrete nodes of the structure. Therefore, the structural stress components are used to describe the local stress state at the material points, while the structural internal force vector is used to describe the overall equivalent internal force on the discrete nodes of the finite element. The two have different physical meanings, but their calculation processes correspond to each other. The above-mentioned method of calculating the strain field from the structural displacement field, obtaining the stress tensor from the material constitutive relation, and further forming the structural internal force vector through element integration and overall assembly is a commonly used method in the dynamic analysis of displacement-type finite element structures. See the relevant literature on finite element theory: Bathe K J. Finite element procedures [M]. Klaus-Jurgen Bathe, 2006.

[0075] In this invention, when the moving body structure exhibits significant geometric nonlinearity, material nonlinearity, or large deformation response, the internal force vector of the structure... This is used to describe the nonlinear structural response under conditions of flexible, deformable, and large deformation. For applications where structural deformation is small, the material response approximately satisfies a proportional relationship between stress and strain, and the effects of geometric nonlinearity are negligible, the internal force vector can be approximated as the product of the stiffness matrix and the nodal displacement vectors, expressed as: , The structural stiffness matrix is ​​obtained by finite element discretization and assembly of the structural material elastic parameters, geometric characteristics, element discretization form and boundary constraints. It is used to characterize the structure's resistance to deformation disturbances.

[0076] When a moving structure exhibits a viscoelastic response, the structural stress depends not only on the current strain but also on the rate of strain change or the historical deformation state. In this case, the time-varying mechanical response of the material can be described using a viscoelastic constitutive relation, the general form of which is:

[0077] in, For structural stress tensor, For the structural strain tensor, Let be the structural strain rate tensor. For the elastic constitutive tensor, For viscous constitutive tensors, This represents the tensor double dot product operation. The stress tensor obtained from this viscoelastic constitutive relation, after element integration and global assembly, forms the corresponding structural internal force vector. It is used to describe the hysteretic deformation and energy dissipation response of viscoelastic structures under impact load.

[0078] When a moving structure exhibits large deformation or nonlinear material response, its recoverable large deformation behavior can be described using hyperelastic constitutive relations. In this case, the structural stress is obtained by differentiating the strain energy density function with respect to the strain variable, and its general form is:

[0079] in, For the second type of Piola-Kirchhoff stress tensor, Let be the strain energy density function. Let be the right Cauchy-Green deformation tensor. Depending on the specific material type, the strain energy density function can be selected from the Neo-Hookean model, the Mooney-Rivlin model, or the Ogden model. The stress tensor obtained from the hyperelastic constitutive relation in the above equation, after element integration, coordinate transformation, and overall assembly, forms the nonlinear structural internal force vector. It is used to describe the large deformation response of hyperelastic structures under hypersonic water impact loads.

[0080] Therefore, the structural dynamics solution in this invention is not limited to rigid or viscoelastic structures. For rigid structures, their overall motion response can be solved. For viscoelastic structures, time-dependent structural stresses and internal forces can be obtained through viscoelastic constitutive relations. For hyperelastic structures, nonlinear stresses and internal forces under large deformation conditions can be obtained through strain energy density functions. All of the above-mentioned structural models can be coupled with equivalent external loads at the fluid-structure interface, and the position of the fluid-structure interface and the foreground body-fitted follower mesh can be updated through structural displacement or overall motion feedback.

[0081] Obtain the structural nodal displacement vector Then, the foreground body-fitting mesh is updated to ensure that the structural surface is consistent with the foreground body-fitting mesh. The foreground body-fitting mesh is not a single node, but is composed of multiple discrete mesh nodes; during the mesh update process, the coordinates of each node in the foreground body-fitting mesh are updated separately.

[0082] When performing multiple fluid-structure interaction sub-iterations within a single physical time step, the foreground body-fitted follower mesh... The position of each node is updated using incremental displacement, and its expression is:

[0083] in, For the foreground body-fitted follower mesh after the s-th sub-iteration at the e-th physical time step. The coordinates of each node, The coordinates of the node after the previous sub-iteration. The node displacement increment corresponding to this sub-iteration satisfies:

[0084] in, and These are the structural node displacements corresponding to the s-th and s-1-th sub-iterations within the e-th physical time step, respectively.

[0085] When only one mesh update is performed within a physical time step, the above equation can be degenerated into:

[0086] By updating all discrete nodes in the foreground body-fitted follower mesh point by point, the geometric consistency between the boundary of the foreground body-fitted follower mesh and the surface of the moving body structure can be maintained, and the mesh in the near-wall region can be guaranteed to change synchronously with the structure movement during the subsequent flow field solution.

[0087] Step 10: Based on the initial water entry velocity, initial water entry angle and initial attitude information of the hypersonic moving body, and combined with Steps 1 to 9, perform data prediction through a two-way strongly coupled pressure-displacement iteration. Step 10.1: Based on the initial entry velocity, initial entry angle, and initial attitude information of the hypersonic vehicle, combined with the initial geometric distribution and boundary conditions of the computational domain from Step 1, as well as the initial position of the free surface and the physical properties of the hypersonic vehicle, determine the time. The initial discrete state. The initial discrete state includes the following three types of data: (1) Initial fluid domain field variables, including the initial pressure field Initial velocity field Initial liquid phase density field Initial mixing density field And the initial three-phase volume fraction field, i.e., the initial air phase volume fraction. Initial liquid volume fraction and initial vapor phase volume fraction ; (2) Initial structural domain state variables, including initial structural displacements Initial structural velocity and initial structural acceleration The variables mentioned above are used to characterize the dynamic state of the discrete nodes or degrees of freedom of the structure at the initial moment; (3) Initial geometric quantities of the foreground body-fitting follower mesh, including the initial coordinates of the nodes of the foreground body-fitting follower mesh. It is used to characterize the spatial position of the foreground body-fitted follower mesh at the initial moment, and does not belong to fluid field variables or structural field variables.

[0088] The above three types of data correspond to the initial input information required for fluid solution, structural response solution, and follower mesh update, respectively. They have different physical meanings and should not be used interchangeably.

[0089] The initial mixing density field satisfies:

[0090] in, and The densities of the air phase and the vapor phase are respectively, and are known constants given in advance. For the basic embodiment of the present invention, the air phase density... With vapor phase density All of them are input as pre-given known constants and are not treated as new unknowns to be solved during the calculation process.

[0091] Step 10.2: Set the time step label to e and the sub-iteration step label to s; when the first physical time step e=1, the time scale corresponds to The first sub-iteration step is s=1; according to the Tait state equation in step five, according to Local sound speed is obtained from the derivative of the state equation. ; ,in, The characteristic length or equivalent element size of the smallest grid cell (excluding redundant cells) in the current three-layer Overset mesh system is used; the CFL number is preferably controlled within the range of 0.2 to 0.6.

[0092] Based on the initial pressure field Velocity field and mixed density field Predicting the velocity field by solving the discrete momentum equation :

[0093] in, The time step size of the first physical time step. It is the gravitational acceleration vector. Given the velocity field The equivalent viscous stress tensor of the structure. The vector tensor product, i.e., the outer product of two velocity vectors, is represented by the equivalent viscous stress tensor, which is defined as:

[0094] in, For a unit second-order tensor, This represents the effective dynamic viscosity. In the initial sub-iteration at the first physical time step e=1 (i.e., s=1), the turbulent additional viscosity is set to zero. At this point, the effective dynamic viscosity is determined by the molecular dynamic viscosity, i.e. At other physical time steps, the effective dynamic viscosity is the sum of the molecular dynamic viscosity and the turbulent additional viscosity, which is calculated using a turbulent viscosity model.

[0095] In the initial transient phase of a hypersonic vehicle entering water, the flow field evolution is mainly controlled by the compressibility of the liquid, the rapid propagation of pressure waves, and fluid-structure interaction. The influence of turbulent shear diffusion on momentum transport within the current physical time step is relatively small. Therefore, the turbulent additional viscosity is taken as zero in this stage, and the effective dynamic viscosity is determined by the molecular dynamic viscosity. As the flow develops, when significant shear layer development, local separation and backflow, or enhanced turbulence occur in the flow field, the influence of turbulent fluctuations on momentum transport cannot be ignored. At this time, the effective dynamic viscosity is the sum of the molecular dynamic viscosity and the turbulent additional viscosity, which is calculated using a turbulent viscosity model.

[0096] Dynamic viscosity of mixture molecules satisfy:

[0097] in, , and These are the pre-defined dynamic viscosity constants for the air phase, water phase, and vapor phase, respectively.

[0098] In obtaining the predicted velocity field Instead of directly reversing the momentum equation to predict the pressure field, a pressure correction factor is introduced. And establish a pressure update relationship:

[0099] Meanwhile, the velocity correction relationship is written as:

[0100] The liquid phase density is updated according to the Tait equation of state described in step five, as it changes with pressure. The liquid phase density is then adjusted at a known pressure. Linear expansion yields the linearized formula for the liquid phase density:

[0101] In the formula .

[0102] According to the Tait state equation described in step five:

[0103] We can obtain:

[0104] In the formula, , , and The parameters of the liquid phase equation of state are given in advance, therefore The coefficient is obtained by direct calculation from the known quantity.

[0105] Furthermore, according to the definition of density of a three-phase mixture:

[0106] In the current inner iteration layer, the three-phase volume fraction is temporarily taken as the known value from the previous iteration layer. , , The liquid phase density is updated and corrected based on the current pressure. Then calculate the mixing density. Then we have:

[0107] Substituting the linearized liquid phase density formula into the above equation, we obtain the formula for calculating the mixing density field:

[0108] Therefore, in the pressure correction substep, the mixing density is not solved separately as an independent unknown, but is updated from the liquid phase density correction result under the condition that the volume fraction of the current sub-iteration is known.

[0109] Subsequently, the velocity correction relation and the mixing density update relation are substituted into the discrete mixing continuity equation:

[0110] It is possible to obtain only the pressure correction amount The pressure correction equation with the only master unknown is:

[0111] Solving the pressure correction equation yields the pressure correction amount for the first sub-iteration at the first physical time step. Then adjust according to the pressure. The pressure field under the first sub-iteration at the first physical time step is obtained. Velocity field Liquid phase density field and mixed density field Thus, the pressure-velocity-density update process of the first sub-iteration in the first physical time step is completed.

[0112] Step 10.3: Based on the pressure field obtained in Step 10.2 And the cavitation evolution model in step six, based on the local pressure of the flow field, determines the phase change source term, combined with the liquid phase density field obtained in step 10.2. The evaporation source term and condensation source term were calculated separately, and the vapor phase volume fraction was obtained. Then, based on the three-phase volume fraction transport equation in step four, the three-phase volume fraction constraint relationship in step three, and the velocity field obtained in step 10.2... The air phase volume fraction field under the first sub-iteration at the first physical time step was calculated. With liquid phase volume fraction field This completes the update of the three-phase volume fraction field within the fluid discrete domain under the current sub-iteration.

[0113] Step 10.4: Based on the fluid domain pressure field obtained in Step 10.2 Combining the fluid-structure interaction relationship from step nine, the interfacial traction vector acting on the structural surface is extracted at the fluid-structure interface. The fluid-side interfacial traction vector is determined by the fluid-side stress tensor and the outward normal unit vector of the fluid-structure interface, i.e. .

[0114] At the physical time step e=1, take (Since the fluid-solid interface load during the initial impact phase of a hypersonic vehicle entering water is mainly determined by the contribution of the normal impact pressure, the tangential viscous stress at the interface is relatively small.) The transient pressure at the fluid-solid interface in the current sub-iteration step. This is the unit vector normal to the fluid-structure interface. Then, the interface traction vector... Substituting the equivalent external load mapping relationship from step nine, and after fluid-structure interface integration and node mapping, the fluid-structure interface traction vector is converted into an equivalent external load vector on the discrete nodes of the structure. .

[0115] Then, by utilizing the equivalent external loads on the structural surface, the structural dynamics equations are solved to obtain the predicted structural nodal displacement vectors. Structural velocity vector and structural acceleration vector To mitigate interface oscillations caused by structural displacement updates during the fluid-structure interaction sub-iteration, the structural nodal displacement vectors predicted by the structural dynamics equations are relaxed and used as the structural nodal displacement vectors obtained in this sub-iteration. The expression for this relaxation is:

[0116] in, Let $\mathbf{e}$ be the structural node displacement vector obtained in the $s$-th sub-iteration at the $e$-th time step. The structural nodal displacement vectors are predicted by the structural dynamics equations at the e-th time step and the s-th sub-iteration. The displacement vector of the structural nodes is obtained in the (s-1)th sub-iteration at the e-th time step. The displacement relaxation coefficient is preferably in the range of 0.1 to 0.7.

[0117] The structural velocity vector will be predicted from the structural dynamics equations. Structural acceleration vector and the structural node displacement vector obtained through relaxation processing As a structural response state quantity.

[0118] Based on the structural node displacement vectors obtained in this sub-iteration, and following the methods in steps two and nine, the vectors are transferred to the foreground body-fitted servo mesh via node mapping or proximity interpolation. The coordinates of each mesh element node in the foreground body-fitted servo mesh are then updated point-by-point to obtain the updated foreground body-fitted servo mesh node coordinates after the first sub-iteration. That is, the new grid configuration state quantity.

[0119] It should be noted that the pressure field, velocity field, and volume fraction field in the foreground body-fitted follower mesh are discrete field variables of the fluid domain, used to describe the local flow state in the near-wall region of the moving body, and are not directly equivalent to the nodal loads actually borne by the moving body structure.

[0120] Step 10.5: Update the foreground body moving mesh node coordinates based on the first sub-iteration obtained in Step 10.4. Following the method in step two, a new three-layer overset mesh system is established. Using the distance-weighted interpolation method from step two, and taking the global field variables of the fluid domain updated in steps 10.2 to 10.4 as source data, the pressure field, velocity field, liquid density field, mixing density field, and volume fraction fields of air, liquid, and vapor phases are transferred and reconstructed to the updated foreground body-fitted follower mesh nodes or elements, thereby obtaining the local flow field state variables corresponding to the foreground body-fitted follower mesh. This completes one sub-iteration step within one time step.

[0121] Step 10.6: Keep the time step unchanged, i.e., e=1, and increment the sub-iteration step by 1, i.e., s=2. Use the local flow field state variables, structural response state variables, and mesh configuration state variables corresponding to the foreground body-fitted follower mesh obtained after completing a sub-iteration step in Step 10.5 as the initial data for this sub-iteration step. Keeping the process unchanged, repeat steps 10.2 to 10.5 to complete the execution of the second sub-iteration step. Use the local flow field state variables, structural response state variables, and mesh configuration state variables corresponding to the foreground body-fitted follower mesh obtained after the previous sub-iteration step as the initial data for the next sub-iteration step. Repeat the above sub-iteration execution process until the number of sub-iteration steps reaches the preset value or the current time step converges, then the iteration update of the time step is completed.

[0122] The convergence criterion for the current time step iteration is: when the interface pressure field and structural displacement field between two adjacent sub-iteration steps respectively satisfy the following convergence conditions:

[0123] in, and These are the convergence thresholds for structural node pressure and structural node displacement, respectively. Indicates the number of steps within the current time step. Discrete pressure vector on the fluid side of the fluid-structure interface during secondary fluid-structure strong coupling sub-iteration; Indicates the number of steps within the current time step. Discrete displacement vector of the fluid-structure interface on the structural side during the secondary fluid-structure strong coupling sub-iteration; This indicates the sequence number of the fluid-structure interaction (FSI) sub-iteration within the current time step. When the above conditions are met, the current time step is considered to have reached the FSI convergence state.

[0124] Step 10.7: At the second physical time step, e=2, and the time scale corresponds to... Its first sub-iteration step is s=1; according to the Tait state equation in step five, according to Local sound speed is obtained from the derivative of the state equation. , This represents the liquid phase density field obtained after the last sub-iteration converges in the first physical time step, and serves as the known input for the starting sub-iteration of the second physical time step.

[0125]

[0126] in, The characteristic length or equivalent cell size of the smallest grid cell (excluding redundant cells) in the three-layer Overset overlapping mesh system obtained from the last sub-iteration of the first physical time step; This represents the maximum value of the velocity vector magnitude on the effective fluid discrete element after removing redundant elements, obtained from the last sub-iteration of the first time step; E is the number of sub-iterations in the first physical time step; the CFL number is preferably controlled within the range of 0.2 to 0.6; the local flow field state variables, structural response state variables, and mesh configuration state variables corresponding to the foreground body-fitted follower mesh obtained from the last sub-iteration of the first physical time step are used as the initial data for the first sub-iteration step of the second physical time step. Replace with Repeat steps 10.2 through 10.6 to complete the iterative update of the second time step.

[0127] Step 10.8: Use the local flow field state variables, structural response state variables, and mesh configuration state variables corresponding to the foreground body-fitted follower mesh after the iterative update of the e-th time step as the initial data for the first sub-iteration step of the e+1-th physical time step; and determine the local sound velocity of the e+1-th physical time step based on the liquid phase density field obtained from the last sub-iteration of the e-th time step, and then calculate its time interval. Subsequently, repeat steps 10.2 to 10.6 using this time interval to complete the iterative update of the e+1-th time step; advance the calculation of subsequent time steps in the above manner until the time step reaches the preset value, and complete the data prediction.

[0128] It should be noted that in each sub-iteration step, after the foreground body-fitted servo mesh is updated with the structural response, the field variables need to be transferred and reconstructed under the new three-layer overset mesh system. This ensures that the updated foreground body-fitted servo mesh can obtain initial values ​​of local flow field state variables such as pressure field, velocity field, liquid density field, mixing density field, and volume fraction fields of air, liquid, and vapor phases from the background mesh domain, while maintaining the continuity and consistency of field variable distribution within the overlapping region of the background mesh domain and the foreground body-fitted servo mesh. Based on this, subsequent sub-iterations or the next physical time step are performed under the new mesh configuration. Combining the liquid-phase compressibility model described in step five, the cavitation evolution model described in step six, and the fluid-structure interaction model described in step nine, a unified iterative update of the pressure field, velocity field, density field, three-phase volume fraction field, structural response, and mesh configuration state variables is completed. This achieves unified time-domain prediction of three-phase flow, liquid-phase compressibility effect, cavitation evolution, fluid-structure interface load transfer, structural response, and foreground body-fitted servo mesh updates during the entry of a hypersonic moving body into water.

[0129] Fluid-structure interface traction extraction and equivalent external load on structure Mapping: When the pressure-velocity iteration converges, at the fluid-structure interface Upper extraction interface pressure Combined with the external normal unit vector Formation of fluid interface traction:

[0130] And mapped to discrete nodes of the structure Equivalent external load (in terms of structural shape function) Weighted integral):

[0131] in, For discrete nodes of the structure The corresponding shape function, For nodes The equivalent external loads of all nodes are assembled to form the equivalent external load vector of the fluid-structure interface, which is used as the external load input term on the right-hand side of the structural dynamics equation to solve for structural displacement and velocity. The structural displacement is then fed back to update the position of the fluid-structure interface and the following mesh, and the process proceeds to the next fluid-structure coupling iteration or the next time step.

[0132] During hypersonic water entry, the peak fluid pressure exhibits microsecond-level variations, while the structural response has a longer inertial timescale, resulting in significant time-domain differences between the two. This invention achieves both dynamic and kinematic consistency at the fluid-structure interface by performing multiple strongly coupled iterations of interface pressure-structural displacement within a single time step, thereby avoiding the hysteresis error and interface mismatch problems commonly found in weakly coupled methods.

[0133] Structural domain dynamics solution: Under the equivalent external load at the fluid-structure interface, the governing equations of structural dynamics are solved to obtain the structural displacement, velocity, and acceleration responses. Based on the structural displacement field, the structural strain field is calculated through the strain-displacement relationship, and then the structural stress tensor is obtained by combining the constitutive relations of the structural materials. The structural stress tensor can be decomposed into normal stress components. , , and shear stress components , , Based on this, the von Mises equivalent stress is used to characterize the structural strength response, and its expression is:

[0134] in, For von Mises equivalent stress; , , These are the normal stress components of the structural material in the x, y, and z directions, respectively. , , These represent the shear stress components of the structural material in the xy, yz, and zx planes, respectively. Each of these stress components is calculated from the displacement response obtained from structural dynamics, combined with strain relations and material constitutive relations.

[0135] The von Mises equivalent stress is used to evaluate the high-stress region and safety margin of the structure. It is a post-processing evaluation index of the structural response results. The von Mises equivalent stress can be compared with the allowable stress or yield strength of the structural material to determine whether the structure meets the strength safety requirements under water impact conditions.

[0136] Time-stepping and numerical stability control: To capture the microsecond-level pressure peaks and cavitation evolution behavior during hypersonic water entry, this invention employs an adaptive time-stepping control strategy based on CFL conditions. The time step is dynamically determined based on the characteristic length of the effective computational grid cells currently involved in the fluid solution, the local velocity modulus of the fluid, and the local sound velocity in the liquid phase. Under hypersonic water entry conditions, the adaptively determined time step is typically on the order of 10⁻⁷ to 10⁻⁶ s to ensure the numerical stability and temporal resolution of pressure wave propagation, phase transition evolution, and fluid-structure interaction calculations.

[0137] Example: This embodiment uses a hypersonic water-entering body with a characteristic diameter D = 0.006 m as an example, with an initial entry velocity of 645.5 m / s and vertical entry. The computational domain size of the background water area is 50D × 30D × 30D, the minimum size of the background mesh is 0.01D, and the minimum size of the foreground mesh is 0.005D. Under the corresponding working condition in this embodiment, after adaptive control according to CFL conditions, the time step stabilizes at... The order of magnitude, and the CFL number is controlled below 0.4. The structural material density is 7850 kg / m³. 3 The elastic modulus is 206 GPa. The actual flow field characteristic evolution process of this embodiment and the flow field characteristic evolution process obtained by the method of this invention are as follows: Figure 6 As shown, the actual flow field characteristics evolution process can be found in the following literature: G. Zhao, Z. Guo, T. Chen, W. Zhang, Experimental investigation on high-speed vertical water entry of spheres with different densities and velocities, Ocean Eng. 295 (2024)116775.

[0138] Figure 6 The results show a precise correspondence between the numerical simulation and experimental findings across four key time phases of the impact process: in the initial impact phase (t=0.05 ms), the calculated results accurately reproduce the shock wave-interface interaction mode observed in the experiment; in the transition phase (t=0.25 ms), the numerical results not only preserve the topological features of the cavitation necking morphology but also restore the details of the shock wave leading edge curvature; and in the deep penetration phase (t=0.45 ms to t=0.9 ms), both maintain a highly consistent cavitation profile and a coherent shock wave propagation trajectory. This multi-stage consistency across different time scales—including cavitation dynamics evolution and the transient flow field structure of the shock wave—strongly validates, from a qualitative perspective, the accuracy, reliability, and physical fidelity of the numerical model constructed in this study for transonic high-speed water entry processes.

[0139] Example 2 This embodiment takes a hyperelastic rotating body at 1000 m / s in air as an example, before it enters the water surface. Figure 7The diagrams show the shock wave and pressure before water entry. The left image shows a photograph of the head shock wave and parameters before and after the wave in a planar cylinder at Mach number 2. Subscript 1 represents parameters upstream of the shock wave, and subscript 2 represents downstream parameters after passing through the shock wave. The right image shows the pressure cloud map in the air medium before entering the water surface, obtained from simulation calculations of a hyperelastic rotating body at 1000 m / s using the method of this invention.

[0140] The calculation results show that a typical bow-shaped shock wave is formed in front of the rotating body's head, and its overall shape is consistent with the experimental phenomenon in the left figure. A relatively short normal shock wave structure can also be observed directly in front of the head's axis; as the bow-shaped shock wave extends to the sides and rear, the tilt of the shock wave gradually increases, corresponding to different shock wave angles. The pressure decreases continuously along the direction of extension. Regarding pressure distribution, a significant high-pressure concentration zone appears in front of the head of the rotating body, with the peak value located in the stagnant region directly in front of the head, gradually decreasing backward along the shock wave's shape. The main reason for this phenomenon is that the shock wave compression causes a sudden change in the fluid state downstream of the shock wave. According to the normal shock wave relationship, the pressure, density, and temperature downstream of the shock wave are all higher than upstream, thus exhibiting a higher pressure level near the head of the rotating body. Conversely, low-pressure zones mainly appear at the head and tail edges of the rotating body. These locations are prone to boundary layer separation, leading to localized overexpansion of the flow, resulting in pressure reduction and the formation of a low-pressure distribution.

[0141] The simulation process includes the following stages: interface fracture stage, cavitation generation stage, high-speed penetration stage, cavitation collapse stage, and structural dynamic response stage. Within the strongly coupled framework of this invention, all of these stages evolve synchronously.

[0142] This invention can output the following key physical quantities: vapor phase volume fraction evolution, pressure field distribution, pressure peak value, water hammer wave propagation process, velocity decay and attitude change of moving body, interface rupture and cavitation morphology evolution characteristics, structural stress, strain and dynamic response, and secondary impact characteristics caused by cavitation collapse.

[0143] The inputs of the method of the present invention include: (1) a three-dimensional geometric model of the moving body; (2) background and foreground meshes; (3) air, water and steam physical property parameters; (4) Tait equation of state parameters; (5) saturated vapor pressure and cavitation model parameters; (6) structural material properties and boundary conditions; (7) working conditions such as water entry velocity and angle.

[0144] The outputs of the method of the present invention include: (1) the evolution of the volume fraction of the vapor phase over time; (2) the pressure field distribution and the peak value of the water hammer pressure; (3) the velocity decay and the change of the attitude of the moving body; (4) the structural stress, strain and dynamic response; (5) the interface rupture and cavitation morphology; and (6) the pressure wave propagation process and the location of cavitation collapse.

[0145] This invention proposes a three-phase flow-structure strongly coupled numerical prediction method for the entry process of hypersonic vehicles into water. This method uses pressure wave propagation as the core prediction driver, synchronously coupling the liquid compressibility effect, cavitation evolution process, and structural dynamic response within a unified computational framework, achieving continuous prediction of key physical quantities throughout the entire hypersonic entry process.

[0146] By introducing a compressible equation of state, the model internally calculates pressure, density, and sound velocity changes self-consistently, and uses the pressure prediction results as a unified criterion for cavitation evolution and a control variable for updating the fluid-solid interface load, thus constructing a closed-loop prediction mechanism for pressure-phase transition-structural response. This mechanism enables the three-phase interface evolution, cavitation generation and collapse, and structural response to advance synergistically within the same timescale, mitigating the problems of pressure peak distortion, cavitation initiation shift, and load prediction lag in traditional non-coupling and weakly coupled methods.

[0147] By combining a three-layer overset grid with a strong pressure-displacement coupling iteration strategy within a single time step, this invention can maintain numerical calculation stability under hypersonic and large displacement water entry conditions, and achieve continuous prediction of moving bodies crossing the gas-liquid interface and deep water regions.

[0148] In summary, this invention is not a simple combination of multiple numerical models, but rather a numerical prediction method for the entire hypersonic water entry process, with pressure wave propagation prediction induced by liquid compressibility as the dominant physical mechanism. It possesses a clear physical framework and good engineering applicability. This method is applicable not only to water entry impact prediction of rigid moving bodies, but also to fluid-structure interaction response analysis of flexible structures and structures with large deformations. It can be used for water entry impact prediction analysis in the engineering design and feasibility study stages, as well as for impact load assessment and structural safety analysis of cross-medium high-speed equipment, protective structures, and energy-absorbing structures.

[0149] The embodiments described are used to illustrate the principles and methods of the present invention, and are not intended to limit the scope of protection of the present invention. Those skilled in the art can make various equivalent modifications to the three-layer Overset mesh structure, three-phase flow model, compressible liquid phase model, cavitation evolution model, structural material constitutive model, coupled iteration form, and time-progression strategy without departing from the core idea of ​​the present invention, and all such modifications should fall within the scope of protection of the present invention.

[0150] Any aspects not covered in this invention are applicable to existing technologies.

Claims

1. A high-speed vehicle water-entry three-phase flow - structure strongly coupled numerical prediction method, characterized in that, The method includes at least the following steps: Step 1: Define the computational domain and boundary conditions; Based on the solid structure and dimensions of the hypersonic vehicle with a streamlined head and axisymmetric shape, a structural model of the vehicle is constructed. The vehicle enters the water from the air, and a computational domain including both the air and water domains is established. The scale of the computational domain is larger than that of the vehicle, with the longitudinal extent of the water domain being 20–30 times the length of the vehicle. The boundary conditions of the computational domain are as follows: (a) Solid wall no-slip boundary: , represents the volume-averaged velocity vector of the three-phase mixture of air, water, and vapor phases, represents the local interfacial normal vector; (b) free surface pressure boundary: , represents the local pressure of the flow field; represents the atmospheric pressure of the air domain; (c) The velocity derivative is zero: , This represents the volume-average velocity vector of a three-phase mixture of air, water, and steam. Represents the local interface normal vector; (d) The pressure derivative is zero: , Indicates the local pressure in the flow field. This is the unit normal vector outside the local boundary; As the moving body moves from the air domain into the water domain, the air phase, water phase, and vapor phase within the computational domain form a fluid domain, while the moving body is the structural domain. The fluid domain uses the following boundary conditions: (1) Free surface: Atmospheric pressure is applied And adopt local pressure of the flow field Normal velocity gradient condition Meanwhile, the gas-liquid interface shape is tracked using a three-phase volume fraction interface capture method; (2) Side walls and bottom solid walls: Apply no-slip conditions ; (3) Top air region: Apply pressure outlet conditions to allow pressure waves to propagate freely outward, satisfying the following conditions. , Indicates the local pressure in the flow field. This is the unit normal vector outside the local boundary; The boundary conditions of the structural domain are as follows: (1) The outer wall of the structure is subjected to motion or deformation response calculation based on the selected structural domain model; when a deformable structural model is used, the structural response is solved according to the structural dynamics equation; when a rigid body type is used, the overall motion response is updated according to the equivalent load of the fluid-structure interface. (2) If there is a fixed region in the structure, then apply displacement constraints to that region. ; (3) The fluid load on the fluid-solid interface is transferred to the structural domain in real time by the interface traction stress calculated by the fluid domain at the interface through a mapping algorithm; Step 2: Establish a three-layer overset mesh system; In step one, a background grid domain is set in the computational domain to describe large-scale water flow, cavitation cloud structure, pressure wave propagation and liquid surface rupture behavior; the moving body is located in the background grid domain, and a foreground body-fitting follower grid is constructed around the surface of the moving body. The foreground body-fitting follower grid is body-fitted according to the geometric shape of the moving body, and its grid boundary is consistent with the outer surface of the moving body. The foreground body-fitting follower mesh consists of multiple element nodes. Each node is used to represent the spatial discrete position in the body-fitting follower mesh. The mesh is formed by connecting adjacent discrete nodes. It is generated according to the geometric model of the moving body at the initial moment and establishes a correspondence with the surface nodes of the structural domain through node mapping relationship, so that the foreground body-fitting follower mesh is updated synchronously with the movement of the moving body structure. A three-layer overset mesh system is constructed between the background mesh domain and the foreground body-fitting mesh. This three-layer overset mesh system includes an extrapolation layer, an interpolation layer, and a receiving layer, used to achieve mesh crossing, variable transfer, and interface mapping coordination for the moving body during large displacement processes. The specific construction method is as follows: Within the overlapping region of the background mesh domain and the foreground body-fitted follower mesh, the layer attributes of each mesh element are determined by their geometric positional relationship: when the geometric center of a mesh element in the background mesh domain falls within the envelope region of the foreground body-fitted follower mesh, the element is marked as an extrapolation layer element, i.e., a redundant element, and does not participate in the solution of fluid domain field variables and interpolation calculation; when the geometric center of a mesh element in the background mesh domain is located outside the envelope region of the foreground body-fitted follower mesh but the distance from the mesh element of the foreground body-fitted follower mesh is not greater than a set value, the element participates in the interpolation calculation as a candidate donor element, and the interpolation layer is formed by the donor elements; nodes or elements in the foreground body-fitted follower mesh that receive the background mesh domain are defined as the receiving layer; The interpolation layer transmits field variables between the background grid domain and the foreground body-fitted follower grid using a distance-weighted interpolation method. These field variables include pressure, velocity components, and volume fractions of air, water, and vapor phases. This enables continuous transmission and consistent updating of physical quantities between the background grid and the foreground body-fitted follower grid, adapting to hypersonic large-displacement water entry conditions. The receiving layer receives the interpolation results output by the interpolation layer and updates the fluid domain field variables on the foreground body-fitted follower grid. Step 3: Define the air-water-steam three-phase volume fraction model; To describe the three-phase mixed flow of air, water, and steam within the aforementioned three-layer overset grid system, a three-phase volume fraction model is defined based on the grid structure. To describe the inlet interface rupture, cavitation generation, and steam phase evolution, a three-phase volume fraction interface capture model is employed. The volume fractions of the air, water, and steam phases are defined as follows: , , And satisfy the three-phase volume fraction constraint relationship: ; The density and viscosity of a mixture are defined as a volume-fraction-weighted average: ， ， in, , , These represent the densities of the air phase, water phase, and vapor phase, respectively. , , These represent the dynamic viscosity of the air phase, water phase, and vapor phase, respectively. and These represent the equivalent density and equivalent viscosity of the three-phase mixture, respectively. Step 4: Define the three-phase volume fraction transport equation; The PLIC interface reconstruction algorithm is employed to ensure the geometric accuracy of the interface shape under high-speed impact. To describe the motion of the air-water-steam three-phase interface, a three-phase volume fraction transport equation is used: ， in, The corresponding markings for the air phase, water phase, and vapor phase; For the first Phase volume fraction; For gradient operators; This represents the volume-average velocity vector of a three-phase mixture of air, water, and steam. As a phase change source term, evaporation and condensation occur only between the liquid and vapor phases, satisfying the following conditions: ; ， Let the phase change mass source term from liquid to vapor phase be denoted as ,but: ， Step 5: Establish a liquid-phase compressibility model based on the Tait equation of state; The Tait equation of state is used to establish a functional relationship between liquid phase pressure and density. The liquid phase density and pressure satisfy a unified Tait equation of state relationship: ， In the formula, For the local pressure of the flow field, This represents the liquid phase density at the corresponding pressure. and These represent the liquid phase density and pressure under reference conditions, respectively. This is an empirical constant related to the bulk modulus of liquids. For exponential coefficients; Furthermore, the local sound velocity can be obtained from the Tait equation of state. expression: ， in, This refers to the local sound velocity in the liquid phase. This refers to the local pressure in the flow field. The density of the liquid phase is calculated using the Tait equation of state. Reference liquid phase density; These are the empirical constants in the Tait equations of state; The above equation represents the pressure response to density change under weakly compressible conditions in the liquid phase, obtained by differentiating the density from the Tait equation of state; this local sound velocity is used to describe the propagation speed of pressure waves in the liquid and is used for subsequent time step control based on CFL conditions. Local sound speed This clearly demonstrates the compressibility behavior of a liquid under pressure, and its changes will directly affect the time step limit based on the CFL number: ， in, The volume-average velocity vector of the three-phase mixture; This refers to the local sound velocity in the liquid phase. This refers to the time step for numerical computation. The characteristic length or equivalent element size of the smallest grid cell (excluding redundant cells) in the current three-layer Overset grid system is used to ensure the numerical stability of pressure wave propagation and velocity field calculation. To ensure the stability of pressure wave propagation calculation during hypersonic water entry, the CFL number is preferably controlled within the range of 0.2 to 0.

6. Step 6: Construct a cavitation evolution model based on local pressure criteria; Hypersonic entry into water is accompanied by significant pressure fluctuations, and the local pressure of the liquid may fall below the saturation pressure, inducing cavitation. Therefore, this invention employs a cavitation evolution model based on a local pressure criterion to describe the vapor phase generation and collapse behavior caused by the reduction in local pressure. Not higher than saturated vapor pressure When the liquid phase evaporates and generates vapor, the phase change source term manifests as an evaporation source term; when the local pressure of the flow field... Higher than saturated vapor pressure At this time, steam condenses and transforms into the liquid phase, and the phase change source term manifests as a condensation source term; vapor phase volume fraction In the evolution equation, the phase transition source term can be expressed as: ， in, For evaporation source term, For condensation source terms, they are defined as follows: when hour, ; ; when hour, ; ; In the formula, , This is an empirical coefficient; Vapor core volume fraction; The equivalent radius of the vapor core; , These are the densities of the liquid phase and the vapor phase, respectively. Step 7: Set up the structural domains and structural material constitutive models; Based on the geometry, structural dimensions, material type, and service response characteristics of the hypersonic moving body, a structural domain model is established, and corresponding structural material constitutive models, mass parameters, damping parameters, initial states, and boundary constraints are set to describe the dynamic response behavior of the moving body under high-pressure impact loads. Step 8: Set the three-phase flow control equations; The three-phase flow control equations specifically include the mixing continuity equation, the momentum equation, and the turbulent viscosity model; Mixed continuity equations: ， Momentum equation: ， in: is the density of a three-phase mixture, used to characterize the equivalent density of the three phases of air, water and steam within a unit volume; For time; The term represents the change in the density of the mixture over time; For gradient operators; The volume-average velocity vector of the three-phase mixture; This represents the rate of change of the momentum of a mixture with respect to time. This refers to the local pressure in the flow field. The molecular dynamic viscosity of the three-phase mixture is obtained by weighting the viscosity of each phase by its volume fraction. Turbulent viscosity, used to describe the enhancing effect of turbulent fluctuations on momentum transport, is calculated using a turbulence model; turbulent viscosity is calculated using the k–ε model or an improved model. It represents turbulent kinetic energy and characterizes the intensity of fluid velocity fluctuations; It represents the turbulent kinetic energy dissipation rate, used to describe the rate at which turbulent kinetic energy is dissipated to smaller scales; Step 9: Establish a fluid-structure interaction model; To reflect the dynamic response of the structure under transient impact loads during hypersonic water entry and its reverse influence on the flow field, a two-way strongly coupled fluid-structure interface model is established; the fluid-structure interface achieves mechanical and kinematic coordination through the following conditions: (1) Kinetic consistency: ; (2) Kinematic consistency: ; (3) Displacement continuity: ; in: , These are the fluid-side stress tensor and the structural-side stress tensor, respectively. The interface traction vector represents the contact force per unit interface area transmitted from the fluid to the structure or fed back from the structure to the fluid. , These represent the velocities at the fluid-solid interface on the fluid side and the structural side, respectively. , These represent the positions of the fluid side and the structural side at the fluid-solid interface, respectively. The interfacial traction stress calculated at the fluid-solid interface in the fluid domain, after interface integration and nodal mapping, is converted into an equivalent external load vector on the structural surface. Its general expression is: ， in, Let T be the shape function matrix corresponding to the discrete nodes of the structure, and let T denote the transpose operation. Let be the interfacial drag vector of the fluid acting on the fluid-solid interface. This represents the interface where the fluid domain and the structural domain come into contact and loads are transferred. Denote the area integral infinitesimal element on the interface; satisfying: ， In the formula, For fluid-side stress tensor, The external normal unit vector of the fluid-structure interface; The structural dynamics equations are: ， in, The structural mass matrix, Here is the structural damping matrix. Let be the structural node displacement vector. The internal force vector of the structure. This is the structural equivalent external load vector. Represents the velocity vector of the structural nodes; Represents the acceleration vector of the structural nodes; Obtain the structural nodal displacement vector Then, the foreground body-fitting follower mesh is updated to make the structure surface consistent with the foreground body-fitting follower mesh; the foreground body-fitting follower mesh is not a single node, but is composed of multiple discrete mesh nodes; during the mesh update process, the coordinates of each node in the foreground body-fitting follower mesh are updated separately. When performing multiple fluid-structure interaction sub-iterations within a single physical time step, the foreground body-fitted follower mesh... The position of each node is updated using incremental displacement, and its expression is: ， in, For the foreground body-fitted follower mesh after the s-th sub-iteration at the e-th physical time step. The coordinates of each node, The coordinates of the node after the previous sub-iteration. The node displacement increment corresponding to this sub-iteration satisfies: ， in, and These are the structural node displacements corresponding to the s-th and s-1-th sub-iterations within the e-th physical time step, respectively; Step 10: Based on the initial water entry velocity, initial water entry angle and initial attitude information of the hypersonic moving body, and combined with Steps 1 to 9, perform data prediction through a two-way strongly coupled pressure-displacement iteration. Step 10.1: Based on the initial entry velocity, initial entry angle, and initial attitude information of the hypersonic vehicle, combined with the initial geometric distribution and boundary conditions of the computational domain from Step 1, as well as the initial position of the free surface and the physical properties of the hypersonic vehicle, determine the time. The initial discrete state; the initial discrete state includes the following three types of data: (1) Initial fluid domain field variables, including the initial pressure field Initial velocity field Initial liquid phase density field Initial mixing density field And the initial three-phase volume fraction field, i.e., the initial air phase volume fraction. Initial liquid volume fraction and initial vapor phase volume fraction ; (2) Initial structural domain state variables, including initial structural displacements Initial structural velocity and initial structural acceleration The above variables are used to characterize the dynamic state of the discrete nodes or degrees of freedom of the structure at the initial moment. (3) Initial geometric quantities of the foreground body-fitting follower mesh, including the initial coordinates of the nodes of the foreground body-fitting follower mesh. It is used to characterize the spatial position of the foreground body-fitted follower mesh at the initial moment; The initial mixing density field satisfies: ， in, and The densities of the air phase and the vapor phase are given in advance as known constants. Step 10.2: Set the time step label to e and the sub-iteration step label to s; in the first time step, e=1, and the time scale corresponds to... The first sub-iteration step is s=1; according to the Tait state equation in step five, according to Local sound speed is obtained from the derivative of the state equation. ; ,in, The characteristic length or equivalent element size of the smallest mesh element excluding redundant elements in the current three-layer Overset mesh system; the CFL number is controlled within the range of 0.2 to 0.6; Based on the initial pressure field Velocity field and mixed density field Predicting the velocity field by solving the discrete momentum equation : ， in, The time step size of the first physical time step. It is the gravitational acceleration vector. Given the velocity field The equivalent viscous stress tensor of the structure; This represents the vector tensor product, which is the cross product of two velocity vectors. The equivalent viscous stress tensor is defined as follows: ， in, For a unit second-order tensor, The effective dynamic viscosity is determined by the molecular dynamic viscosity in the initial sub-iteration at the first physical time step e=1. The turbulent additional viscosity is set to zero, and the effective dynamic viscosity is then determined by the molecular dynamic viscosity. At other time steps, the effective dynamic viscosity is the sum of the molecular dynamic viscosity and the turbulent additional viscosity, which is calculated using a turbulent viscosity model. Dynamic viscosity of mixture molecules satisfy: ， in, , and The pre-defined dynamic viscosity constants are for the air phase, water phase, and vapor phase, respectively. In obtaining the predicted velocity field Instead of directly reversing the momentum equation to predict the pressure field, a pressure correction factor is introduced. And establish a pressure update relationship: ， Meanwhile, the velocity correction relationship is written as: ， The liquid phase density is updated according to the Tait equation of state described in step five, as it changes with pressure; the liquid phase density is then adjusted at a known pressure. Linear expansion yields the linearized formula for the liquid phase density: ， In the formula, ; According to the Tait state equation described in step five: ， We can obtain: ， In the formula, , , and The parameters of the liquid phase equation of state are given in advance, therefore The coefficients are obtained by direct calculation from known quantities; Furthermore, according to the definition of density of a three-phase mixture: ， In the current inner iteration layer, the three-phase volume fraction is temporarily taken as the known value from the previous iteration layer. , , The liquid phase density is updated and corrected based on the current pressure. Then calculate the mixing density. Then we have: ， Substituting the linearized liquid phase density formula into the above equation, we obtain the formula for calculating the mixing density field: ， Therefore, in the pressure correction substep, the mixing density is not solved separately as an independent unknown, but is updated from the liquid phase density correction result under the condition that the volume fraction of the current sub-iteration is known; Subsequently, the velocity correction relation and the mixing density update relation are substituted into the discrete mixing continuity equation: ， It is possible to obtain only the pressure correction amount The pressure correction equation with the only master unknown is: ， Solving the pressure correction equation yields the pressure correction amount for the first sub-iteration at the first physical time step. Then adjust according to the pressure. The pressure field under the first sub-iteration at the first physical time step is obtained. Velocity field Liquid phase density field and mixed density field Thus, the pressure-velocity-density update process of the first sub-iteration in the first physical time step is completed. Step 10.3: Based on the pressure field obtained in Step 10.2 And the cavitation evolution model in step six, based on the local pressure of the flow field, determines the phase change source term, combined with the liquid phase density field obtained in step 10.

2. The evaporation source term and condensation source term were calculated separately, and the vapor phase volume fraction was obtained. Then, based on the three-phase volume fraction transport equation in step four, the three-phase volume fraction constraint relationship in step three, and the velocity field obtained in step 10.2... The air phase volume fraction field under the first sub-iteration at the first physical time step was calculated. With liquid phase volume fraction field This completes the update of the three-phase volume fraction field within the fluid discrete domain under the current sub-iteration; Step 10.4: Based on the fluid domain pressure field obtained in Step 10.2 Combining the fluid-structure interaction relationship from step nine, the interface traction vector acting on the structural surface is extracted at the fluid-structure interface; the fluid-side interface traction vector is determined by the fluid-side stress tensor and the outward normal unit vector of the fluid-structure interface, i.e. ; At the physical time step e=1, take ,in, The transient pressure at the fluid-solid interface in the current sub-iteration step. The external normal unit vector of the fluid-structure interface; then the interface traction vector. Substituting the equivalent external load mapping relationship from step nine, and after fluid-structure interface integration and node mapping, the fluid-structure interface traction vector is converted into an equivalent external load vector on the discrete nodes of the structure. ; Then, by utilizing the equivalent external loads on the structural surface, the structural dynamics equations are solved to obtain the predicted structural nodal displacement vectors. Structural velocity vector and structural acceleration vector To mitigate interface oscillations caused by structural displacement updates during the fluid-structure interaction sub-iteration, the structural nodal displacement vectors predicted by the structural dynamics equations are relaxed and used as the structural nodal displacement vectors obtained in this sub-iteration. The expression for this relaxation is: ， in, Let $\mathbf{e}$ be the structural node displacement vector obtained in the $s$-th sub-iteration at the $e$-th time step. The structural nodal displacement vectors are predicted by the structural dynamics equations at the e-th time step and the s-th sub-iteration. The structural node displacement vector is obtained in the (s-1)th sub-iteration at the e-th time step. The structural velocity vector will be predicted from the structural dynamics equations. Structural acceleration vector and the structural node displacement vector obtained through relaxation processing As a structural response state quantity; Based on the structural node displacement vectors obtained in this sub-iteration, and following the methods in steps two and nine, the vectors are transferred to the foreground body-fitted servo mesh via node mapping or proximity interpolation. The coordinates of each mesh element node in the foreground body-fitted servo mesh are then updated point-by-point to obtain the updated foreground body-fitted servo mesh node coordinates after the first sub-iteration. , that is, the new grid configuration state quantity; Step 10.5: Update the foreground body moving mesh node coordinates based on the first sub-iteration obtained in Step 10.

4. Following the method in step two, a new three-layer overset mesh system is established. Using the distance weight interpolation method in step two, the global field variables of the fluid domain updated in steps 10.2 to 10.4 are used as source data. The pressure field, velocity field, liquid density field, mixing density field, and volume fraction fields of air, liquid, and vapor phases are transferred and reconstructed to the updated foreground body-fitted follower mesh nodes or elements, thereby obtaining the local flow field state variables corresponding to the foreground body-fitted follower mesh. This completes one sub-iteration step within one time step. Step 10.6: Keep the time step unchanged, i.e., e=1, and increment the sub-iteration step by 1, i.e., s=2. Use the local flow field state variables, structural response state variables, and mesh configuration state variables corresponding to the foreground body-fitted follower mesh obtained after completing a sub-iteration step in Step 10.5 as the initial data for this sub-iteration step. Keeping the process unchanged, repeat steps 10.2 to 10.5 to complete the execution of the second sub-iteration step; take the local flow field state, structural response state, and mesh configuration state of the foreground body-fitting follower mesh obtained after the previous sub-iteration step as the initial data for the next sub-iteration step, and repeat the above sub-iteration execution process until the number of sub-iteration steps reaches the preset value or the current time step converges, then the iteration update of the time step is completed; The convergence criterion for the current time step iteration is: when the interface pressure field and structural displacement field between two adjacent sub-iteration steps respectively satisfy the following convergence conditions: ， in, and These are the convergence thresholds for structural node pressure and structural node displacement, respectively. Indicates the number of steps within the current time step. Discrete pressure vector on the fluid side of the fluid-structure interface during secondary fluid-structure strong coupling sub-iteration; Indicates the number of steps within the current time step. Discrete displacement vector of the fluid-structure interface on the structural side during the secondary fluid-structure strong coupling sub-iteration; This indicates the sequence number of the fluid-structure interaction (FSI) sub-iteration within the current time step; when the above conditions are met, the current time step is considered to have reached the FSI convergence state. Step 10.7: At the second physical time step, e=2, and the time scale corresponds to... Its first sub-iteration step is s=1; according to the Tait state equation in step five, according to Local sound speed is obtained from the derivative of the state equation. , This represents the liquid phase density field obtained after the last sub-iteration of the first physical time step converges, and serves as the known input for the starting sub-iteration of the second physical time step. ， in, The characteristic length or equivalent cell size of the smallest grid cell (excluding redundant cells) in the three-layer Overset overlapping mesh system obtained from the last sub-iteration of the first physical time step; This represents the maximum value of the velocity vector magnitude on the effective fluid discrete element after removing redundant elements, obtained from the last sub-iteration of the first time step; E is the number of sub-iterations in the first physical time step; the CFL number is preferably controlled within the range of 0.2 to 0.6; the local flow field state variables, structural response state variables, and mesh configuration state variables corresponding to the foreground body-fitted follower mesh obtained from the last sub-iteration of the first physical time step are used as the initial data for the first sub-iteration step of the second physical time step. Replace with Repeat steps 10.2-10.6 to complete the iterative update of the second time step; Step 10.8: Take the local flow field state variables, structural response state variables, and mesh configuration state variables corresponding to the foreground body-fitted follower mesh after the iterative update of the e-th time step as the initial data for the first sub-iteration step of the e+1-th physical time step; and determine the local sound velocity of the e+1-th physical time step based on the liquid phase density field obtained in the last sub-iteration of the e-th time step, and then calculate its time interval; subsequently, repeat steps 10.2 to 10.6 using this time interval to complete the iterative update of the e+1-th time step; advance the calculation of subsequent time steps in the above manner until the time step reaches the preset value and the data prediction is completed.

2. The numerical prediction method for strong coupling of three-phase flow and structure for hypersonic moving bodies entering water according to claim 1, characterized in that, In step two, the field variables obtained by the foreground body-fitted follower mesh through the receiver layer are obtained by the interpolation layer using the distance-weighted interpolation method: the field variable of a mesh node of the foreground body-fitted follower mesh is set as follows. ,Should From the background grid domain, the distance to the grid node is not greater than a set threshold. The interpolation formula is obtained by weighted averaging of the field variables of each donation unit; the interpolation formula is: ， in, For the first The field variable values ​​corresponding to each donation unit; The first kernel, determined based on the inverse distance ratio or the RBF kernel function, is... Weight of each donation unit: ， For the first The distance between each donation unit and the foreground-fitted moving mesh node. For the first The distance between each donation unit and the foreground-fitted follow-mesh node.

3. The numerical prediction method for strong coupling of three-phase flow and structure for hypersonic moving bodies entering water according to claim 1, characterized in that, Empirical constants in Tait's equations of state The preferred value range is (2.5~3.5)×10 8 Pa, exponential coefficient The value range is 6 to 8.

4. The numerical prediction method for strong coupling of three-phase flow and structure for hypersonic moving bodies entering water according to claim 1, characterized in that, The structural domain model uses a rigid body type for motion response prediction or a deformable structural model for structural dynamic response analysis.

5. The numerical prediction method for strong coupling of three-phase flow and structure for hypersonic moving bodies entering water according to claim 4, characterized in that, When the structural domain model uses a rigid body type for motion response prediction, in step nine... , This is the structural stiffness matrix.

6. The numerical prediction method for strong coupling of three-phase flow and structure for hypersonic moving bodies entering water according to claim 4, characterized in that, When a deformable structural model is used, the structural material constitutive model can be selected as a viscoelastic constitutive model or a hyperelastic constitutive model according to the material type and deformation characteristics; wherein: (1) When structural materials exhibit stress-strain rate-related, time-related deformation, hysteresis response or energy dissipation characteristics under impact load, a viscoelastic constitutive model is used for description. (2) When the structural material is under large deformation conditions, the hyperelastic constitutive model is used to describe it in order to characterize the recoverable deformation characteristics of the material in the large strain range.

7. The numerical prediction method for strong coupling of three-phase flow and structure for hypersonic moving bodies entering water according to claim 6, characterized in that, When a moving structure exhibits a viscoelastic response, the structural stress is related not only to the current strain but also to the rate of strain change or the historical deformation state. In this case, the time-varying mechanical response of the material is described using a viscoelastic constitutive relation, the general form of which is: ， in, For structural stress tensor, For the structural strain tensor, Let be the structural strain rate tensor. For the elastic constitutive tensor, For viscous constitutive tensors, This represents the tensor double dot product operation; the stress tensor obtained from this viscoelastic constitutive relation, after element integration and global assembly, forms the corresponding structural internal force vector. It is used to describe the hysteretic deformation and energy dissipation response of viscoelastic structures under impact load.

8. The numerical prediction method for strong coupling of three-phase flow and structure for hypersonic moving bodies entering water according to claim 6, characterized in that, When a moving structure exhibits large deformation or nonlinear material response, its recoverable large deformation behavior is described using a hyperelastic constitutive relation. In this case, the structural stress is obtained by differentiating the strain energy density function with respect to the strain variable, and its general form is: ， in, For the second type of Piola-Kirchhoff stress tensor, Let be the strain energy density function. Let be the right Cauchy-Green deformation tensor; depending on the specific material type, the strain energy density function can be selected from the Neo-Hookean model, the Mooney-Rivlin model, or the Ogden model; the stress tensor obtained from the hyperelastic constitutive relation in the above equation, after element integration, coordinate transformation, and overall assembly, forms a nonlinear structural internal force vector. It is used to describe the large deformation response of hyperelastic structures under hypersonic water impact loads.

9. The numerical prediction method for strong coupling of three-phase flow and structure for hypersonic moving bodies entering water according to claim 1, characterized in that, In step 10.4, is the displacement relaxation coefficient. The value range is 0.1 to 0.7.

Images