Wedge-shaped body water entry slamming BEM-SPH coupling numerical method based on overlapping computational domain thought

Abstract

The invention relates to the technical field of computational fluid mechanics, and discloses a wedge-shaped body water entry slamming BEM-SPH coupling numerical method based on an overlapping computational domain idea. According to a specific simulation problem, a corresponding numerical model is established, a water entry slamming BEM model under a full nonlinear boundary condition is adopted to calculate a slamming load in the initial water entry process of the wedge, after high-amplitude transient slamming force capture is completed, internal flow field distribution and control point transmission speed are obtained according to speed potential and free liquid level elevation calculated by the BEM, and the initial water entry process of the wedge is completed. And an SPH calculation model is initialized by taking the result as input, so that wedge-shaped body motion attitude simulation and hydrodynamic force stress analysis after slamming occur are realized. According to the method, high-frequency pressure oscillation errors of a single SPH model in water-entry slamming peak capture are avoided, the defect that motion simulation after slamming cannot be achieved through a single BEM model is overcome, and effective simulation of the whole water-entry slamming process of the wedge is achieved.

Description

Technical Field

[0001] The present invention relates to the technical field of computational fluid dynamics, and in particular to a BEM-SPH coupled numerical method for wedge-shaped body water slamming based on the concept of overlapping computational domains. Background Art

[0002] Water slamming is the violent interaction between fluid and structure caused by the relative motion of fluid and structure within a short period of time. The slamming process is extremely brief and is characterized by dramatic changes in the wetted surface of the structure, strong nonlinear motion of the fluid, curling and even fragmentation of the free surface, and strong nonlinear coupling between the fluid and the structure. It is a transient, strongly nonlinear hydrodynamic problem. Intense slamming loads can easily cause local deformation and damage to the structure, affecting its overall motion performance. The water slamming problem of fluid and structure is very common in engineering. Examples include the water slamming of underwater vehicles during launch, the bow slamming of ships sailing in waves, and the slamming of seaplanes landing on the sea. These real-world slamming scenarios can be approximated by the physical model of the water slamming of wedge-shaped structures. Therefore, the study of simulation methods for the water slamming problem of wedge-shaped structures is of great significance in the field of marine engineering, providing theoretical guidance and technical support for the development of high-tech equipment.

[0003] The SPH (Smoothed Particle Hydrodynamics) method, a meshless technique that has garnered significant attention in recent years, can automatically track the free surface, exhibits excellent conservation properties, and avoids numerical dissipation caused by discretization errors in convection terms. For water slamming, a typical example of severe fluid-structure interaction, the SPH method demonstrates significant advantages. However, due to the extremely short duration of high-amplitude load responses during slamming, SPH numerical simulations suffer from severe pressure oscillations, resulting in inaccurate capture of the slamming load peak and hindering analysis of the cross-media load-bearing capacity of structures. The Boundary Element Method (BEM), a relatively well-developed method with a relatively sophisticated theoretical foundation and numerical techniques, has been widely recognized for its accuracy in simulating wedge-shaped body slamming. However, due to the continuity requirements of the boundary mesh, it cannot simulate structural motion and subsequent flow field after the slamming event. Furthermore, it lacks the ability to simulate the liquid surface curling and droplet splashing induced by the slamming event. Summary of the Invention

[0004] To solve the problems raised in the background technology, the present invention is implemented through the following technical solutions: a BEM-SPH coupled numerical method for wedge-shaped body water slamming based on the overlapping computational domain concept, including:

[0005] S1. Initialize the water impact BEM calculation model under full nonlinear boundary conditions in the extended coordinate system, set the initial state parameters of the BEM model, and set the BEM stage simulation duration according to the physical height of the wedge;

[0006] S2. At each computational time step, the boundary integral equations constructed with fully nonlinear boundary conditions are solved to obtain the velocity potential at the object surface boundary and the normal derivative of the velocity potential at the free surface. The hybrid Euler-Lagrangian method is used to update the liquid surface lift and velocity potential changes during the water entry process. With the help of auxiliary functions, the slamming pressure on the object surface is calculated. The slamming force is obtained by surface integral and substituted into the Newtonian equation of motion to obtain the motion response of the wedge. The Runge-Kutta method is used to complete time stepping, realizing the time-domain numerical simulation of the wedge entry stage of the BEM.

[0007] S3, coupled time step, calculate the source strength of each boundary element, limit the integration domain of the influence coefficient to [-1, 1] with the help of shape function, and use the six-node Gauss-Legendre formula for numerical integration;

[0008] S4, coupled time step, convert the BEM boundary geometric coordinates to the physical coordinate system, with the free liquid surface elevation and the current position of the wedge as the geometric boundary. The SPH tank boundary is the same as the BEM model, and the internal space is filled with SPH fluid particles to form an SPH particle space distribution geometric model;

[0009] S5, coupled time step, using the SPH fluid particles converted back to the extended coordinate system as control points, calculating the induced velocity of each boundary element on each control point, and further combining the wedge entry velocity to obtain the transfer velocity vector field from the BEM stage to the SPH stage;

[0010] S6. Initialize the SPH calculation model, assign the transfer velocity vector field calculated by the BEM model to the corresponding fluid particles, set the water density, sound velocity coefficient, density diffusion coefficient and water kinematic viscosity required for SPH simulation, and select the smoothing function type and smoothing length;

[0011] S7. Start the SPH model, calculate the particle density, pressure gradient force, and viscous force, obtain the velocity increment and density increment of all fluid particles, update the particle density and velocity according to the prediction-correction time format, complete the iteration of the particle position, and enter the solution of the next time step until the set calculation time is reached.

[0012] Furthermore, the initial state parameters of the BEM model include the mass of the wedge, the bottom lift angle, the water entry speed, the water entry angle, the water tank width, the wall height, the number of surface grids, the number of free liquid surface grids, the number of side wall grids, the number of water tank bottom grids, and the calculation time.

[0013] Furthermore, the geometric dimensions and velocity variables are established in the extended coordinate system to ensure that the BEM simulation maintains a reasonable and stable grid scale from the initial water contact to the liquid surface rise. The conversion relationship between the relevant physical quantities between the extended coordinate system and the physical coordinate system is:

[0014]

[0015] Among them, x, z and v are the coordinate values ​​and velocities of the BEM model in the physical coordinate system. and To expand the grid node coordinates and velocities in the coordinate system, s is the vertical distance of the wedge entering the water, s = ∫Wdt, and W is the water entry velocity.

[0016] Furthermore, a local coordinate system is introduced on the free surface The liquid level is updated during the BEM model calculation process, where: The axis is on the tangent plane of the free surface, Axis perpendicular to Axis, and pointing outside the domain, local coordinate system The conversion relationship with the Cartesian coordinate system is as follows:

[0017]

[0018] in, and for The components of the unit vector of the axis in the Cartesian coordinate system, and for The components of the unit vector of the axis in the Cartesian coordinate system; in the local coordinate system, the node motion is updated using material derivatives, and the free surface position and velocity potential are updated using the following formulas:

[0019]

[0020] in, is the velocity potential, is the wave surface elevation.

[0021] Furthermore, in said S2, As a harmonic function, the auxiliary function method is used to solve it to ensure the accuracy of the impact force calculation during the dynamic change of the grid. Split into:

[0022]

[0023] Among them, the auxiliary functions χ1 and χ2 satisfy the Laplace equation in the watershed and meet the following boundary conditions:

[0024]

[0025] Where n is the normal direction of the boundary element mesh pointing to the fluid domain, It is along The axis weight, Yes and is the velocity potential along Axis and The axis's weight;

[0026] Through Green's third formula, the boundary condition of the auxiliary function is converted into a boundary integral equation. After grid discretization, numerical integration, collocation transformation and LU decomposition, the solutions of χ1 and χ2 are obtained. Then, the vertical slamming force on the wedge is obtained. The specific formula is:

[0027]

[0028] Among them, F z is the vertical slamming force, l B is the boundary area of ​​the object surface, p is the impact pressure on the object surface, and ρ is the water density.

[0029] Furthermore, in S3, the linear algebraic equation is solved on the boundary surface element. Calculate the surface element source strength σ j , where the influence coefficient The isoparametric method is used to solve c using a two-node bar element. ij , the specific calculation formula is:

[0030]

[0031] Among them, h k is the shape function. The specific form of the one-dimensional shape function is J j (ξ n ) is the Jacobian matrix, ξ n is the coordinate of the Gaussian point, w n is the corresponding integral weight, is the horizontal coordinate value of the source point q, is the vertical coordinate value of the source point q, is the horizontal coordinate value of the node j of the unit, is the vertical coordinate value of the node of unit j, k is the node number of unit j, is the field point p i The horizontal coordinate value of is the field point p i The vertical coordinate value of The source point q and the field point p i The distance between them.

[0032] Furthermore, in the above S4, the BEM free surface boundary is cut off from the jet reversal point, and the physical coordinates of the cut off free surface boundary are used as the upper boundary of the SPH fluid domain, and the internal fluid particles are filled.

[0033] Furthermore, in S5, by accumulating the surface element source strength σ j With each boundary grid j control point P i Unit induced speed Get the induced velocity at the control point Unit induced speed The specific calculation formula is:

[0034]

[0035] in, is the control point P i The horizontal coordinate value of is the control point P i The vertical coordinate value of The source point q and the control point P i The distance between them.

[0036] Furthermore, in S6, the velocity initialization of the SPH fluid particles is completed with the help of the particle number Idp, thereby realizing the transfer of the BEM calculation results to the SPH calculation model.

[0037] Furthermore, in S7, the SPH format continuity equation used contains a density diffusion term represented by the total density to improve the pressure field distribution near the wall boundary. The SPH fluid particle density increment is calculated as follows:

[0038]

[0039] Among them, ρ i , p i , c i , v i , r i is the density, pressure, sound speed, velocity vector and position vector of particle i; m j ,ρ j , p j , v j , r j is the mass, density, pressure, velocity vector and position vector of particle j; v ij , r ij , z ij are the velocity vector difference, position vector difference, total density difference, static density difference, hydrostatic pressure difference, and vertical distance difference between particles i and j, respectively; N is the total number of particles in the support domain of particle i, W ijis the smooth function of the influence of particle j on particle i, δ is the density diffusion coefficient, ψ ij is the density diffusion term, h is the smoothing length, γ=7, ρ0 is the standard density, c0 is the speed of sound corresponding to the standard density, h swl is the maximum still water level, coef sound is the speed of sound coefficient.

[0040] Furthermore, in S7, the SPH format momentum equation used contains pressure terms, laminar viscous stress terms, and Reynolds stress terms. The SPH fluid particle velocity increment is calculated as follows:

[0041]

[0042] Among them, υ0 is the kinematic viscosity of water, η 2 =0.01h 2 , τ ij is the sub-grid stress tensor.

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

[0044] In the numerical study of the slamming problem of a wedge entering water, a BEM-SPH coupling method based on the overlapping computational domain concept proposed in this paper achieves effective numerical simulation of the entire wedge entry process. This method fully utilizes the BEM method's ability to capture the peak slamming load of the wedge entering water with high precision, avoids the numerical error caused by the high-frequency pressure oscillation of the SPH method, and realizes the motion simulation and load simulation of the structure after the slamming peak is captured. This provides a means for evaluating the dynamic performance of the wedge after entering water. At the same time, the SPH method can also effectively capture the liquid surface curling and droplet splashing after the development of the slamming jet, realizing the evolution simulation of the strong turbulent flow field caused by the slamming. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 Schematic diagram of the BEM numerical model of the wedge-shaped body slamming into water according to the present invention;

[0046] Figure 2 This is a schematic diagram of the arrangement of pressure measurement points on the surface of the wedge body of the present invention;

[0047] Figure 3 This is the time history curve of the impact pressure coefficient at each measuring point in the BEM stage of the present invention;

[0048] Figure 4 The calculation results of the influence coefficient, velocity potential and source strength during the coupling time step of the present invention;

[0049] Figure 5 It is a schematic diagram of the free liquid surface cutoff of the present invention;

[0050] Figure 6 Schematic diagram of the initial spatial distribution of SPH particles in the present invention;

[0051] Figure 7 It is the vector diagram of fluid particle transfer velocity during the coupling time step of the present invention;

[0052] Figure 8 The particle evolution diagram of the SPH stage at each typical moment of the present invention;

[0053] Figure 9 Schematic diagram of the simulation of the wedge-shaped body entering water and slamming. DETAILED DESCRIPTION

[0054] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0055] The embodiment of the BEM-SPH coupled numerical method for wedge-shaped body water slamming based on the overlapping computational domain concept is as follows:

[0056] See also Figures 1-9 The BEM-SPH coupled numerical method for wedge-shaped body slamming into water based on the overlapping computational domain concept includes:

[0057] Step 1: Initialize the water impact BEM calculation model under fully nonlinear boundary conditions in the extended coordinate system. Set the initial state parameters of the BEM model, including the wedge mass, bottom lift angle, water entry velocity, water entry angle, flume width, wall height, number of surface grids, number of free surface grids, number of side wall grids, and number of flume bottom grids. The BEM simulation duration is set by the physical height of the wedge.

[0058] The geometric dimensions and velocity variables are established in the extended coordinate system to ensure that the BEM simulation maintains a reasonable and stable grid scale from the initial water contact to the liquid surface rise. The conversion relationship between the relevant physical quantities between the extended coordinate system and the physical coordinate system is:

[0059]

[0060] Among them, x, z and v are the coordinate values ​​and velocities of the BEM model in the physical coordinate system. and To expand the grid node coordinates and velocities in the coordinate system, s is the vertical distance of the wedge entering the water, s = ∫Wdt, and W is the water entry velocity.

[0061] Step 2: At each calculation time step, solve the boundary integral equation constructed by the full nonlinear boundary condition, calculate the velocity potential of the object surface boundary and the unknown term of the velocity potential normal derivative of the free liquid surface; use the hybrid Euler-Lagrangian method to update the free liquid surface, solve the boundary value problem with the Euler method in each time step, and use the Lagrangian method to track the free liquid surface in the subsequent time step to achieve effective simulation of the large deformation free liquid surface; calculate with the help of auxiliary functions Solve the time difference problem of dynamic grids and complete the coupled analysis of fluid flow and structural motion. Use the Runge-Kutta method to time step to the specified coupling time step to complete the BEM calculation and analysis of the initial stage of the wedge entering the water.

[0062] Velocity potential under fully nonlinear boundary conditions The solution conditions are:

[0063]

[0064] in, and To expand the grid node coordinates and velocity potential in the coordinate system, is the wave surface elevation, s is the vertical distance of the wedge into the water, s=∫Wdt, W is the water entry velocity, g is the acceleration of gravity, n is the normal direction of the object surface pointing to the fluid domain, It is along Axis component; S B is the object surface boundary, S F is the free liquid surface boundary, S C The far control surface boundary.

[0065] According to Green's third formula, the value of the velocity potential satisfying the Laplace equation at any point in the domain can be expressed using the value on the boundary and the normal derivative. Therefore, the boundary condition of the velocity potential can be rewritten as the following boundary integral equation:

[0066]

[0067] l q =l B +l F +l C ;

[0068] Where A is the solid angle, R pq is the distance from the source point q to the field point p on the integral line, l B is the integration interval of the object surface boundary, l F is the integration interval of the free liquid surface boundary, l C is the integration interval of the far control surface boundary, l q is the integral interval of the entire boundary, l q =l B+l F +l C After the grid is discretized, the equation is organized into a matrix form:

[0069]

[0070] Among them, H corresponds to Item, G corresponds to Item, N d is the total number of boundary grids in the BEM calculation model, Δ q Refers to a grid. The normal derivative of the velocity potential on the object surface at each moment is known through boundary conditions. and the velocity potential on the free surface Using the collocation method, we move the unknown terms to the left side of the equation and the known terms to the right side of the equation, and we get:

[0071]

[0072] The subscript b indicates that the node is on the object surface or the boundary of the remote control surface, and the subscript f indicates that the node is on the free surface. The solution to this linear system of equations is obtained by LU decomposition and back substitution, yielding the velocity potential at the object surface boundary and the normal derivative of the velocity potential on the free surface.

[0073] Considering the liquid surface rise and jet deformation during the water impact process, a local coordinate system is introduced on the free surface. To achieve the liquid level update during the calculation process, where The axis is on the tangent plane of the free surface, Axis perpendicular to Axis, and points outside the domain. In the local coordinate system, the node motion is updated using the material derivative method. The transformation relationship between the local coordinate system and the Cartesian coordinate system is as follows:

[0074]

[0075] in, and for The components of the unit vector of the axis in the Cartesian coordinate system, and for The components of the unit vector of the axis in the Cartesian coordinate system. In the local coordinate system, the free surface position and velocity potential are updated using the following formulas:

[0076]

[0077] in, is the wave surface elevation in the local coordinate system.

[0078] In order to realize the coupled analysis of fluid flow and structural motion, it is necessary to integrate the surface pressure to calculate the external force, specifically the spatial derivative and time partial derivative of the velocity potential. The present invention calculates the spatial derivative of the velocity potential by direct interpolation of the grid and uses the auxiliary function method to calculate Specifically:

[0079]

[0080] The auxiliary functions χ1 and χ2 satisfy the Laplace equation in the flow domain and meet the following boundary conditions:

[0081]

[0082] Adoption and speed potential The same solution strategy is used to solve the boundary value problem composed of auxiliary functions χ1 and χ2, and we get The solution to , and then, the vertical hydrodynamic force (slamming force) on the wedge is obtained:

[0083]

[0084] Among them, F z is the vertical slamming force, p is the slamming pressure on the surface of the object, and ρ is the water density.

[0085] Substituting into Newton's equation of motion, we get the motion response of the wedge:

[0086]

[0087] Where m is the mass of the wedge, F g is the weight of the wedge.

[0088] Step 3: Couple time steps and solve linear algebraic equations on the BEM boundary Calculate the surface element source strength σ j , where the influence coefficient With the help of shape functions, the integration domain is limited to [-1,1], and the six-node Gauss-Legendre formula is used to perform numerical integration of the two-dimensional Green's function on the integration line. The integral operation is converted into a multinomial summation to complete the solution of the influence coefficient. ij The matrix is ​​decomposed by LU, and the source strength σ is obtained by back substitution j .

[0089] Among them, the isoparametric method is used to solve c using a two-node bar element. ij , the specific calculation formula is:

[0090]

[0091] Among them, h k is the shape function. The specific form of the one-dimensional shape function is J j (ξ n ) is the Jacobian matrix, ξ n are the coordinates of the Gaussian point, is the horizontal coordinate value of the source point q, is the vertical coordinate value of the source point q, is the horizontal coordinate value of the node j of the unit, is the vertical coordinate value of the node of unit j, k is the node number of unit j, is the field point p i The horizontal coordinate value of is the field point p i The vertical coordinate value of The source point q and the field point p i The distance between n is the corresponding integral weight, and the specific values ​​are shown in the following table:

[0092] n <![CDATA[ξ n ]]> <![CDATA[w n ]]> 1 -0.9324695142031520 0.1713244923791700 2 -0.6612093864662640 0.3607615730481380 3 -0.2386191860831960 0.4679139345726910 4 0.2386191860831960 0.4679139345726910 5 0.6612093864662640 0.3607615730481380 6 0.9324695142031520 0.1713244923791700

[0093] Step 4: During the coupled time step, the BEM boundary geometry in the extended coordinate system is converted to the physical coordinate system. The free surface is truncated at the jet reversal point, and the truncated surface elevation data for the BEM model is output. The elevation data is used as the free surface geometric boundary. The wedge position and the water tank boundary are the same as in the BEM model. The internal space is filled with SPH fluid particles, forming an SPH particle spatial distribution geometric model.

[0094] Step 5: Couple the time steps and convert the spatial coordinates of the SPH fluid particles back to the extended coordinate system, using the SPH fluid particles in the extended coordinate system as the control point P i , calculate the unit induced velocity of each boundary grid to the control point and the source strength σ of the corresponding boundary grid j The induced velocity at the control point is obtained by multiplication, and then further converted to the physical coordinate system to obtain the transfer velocity vector field v(P i ), the specific calculation formula is:

[0095]

[0096] in, is the control point P i The horizontal coordinate value of is the control point P i The vertical coordinate value of is the source point q and the control point P i The distance between them.

[0097] Step 6: Initialize the SPH calculation model, assign the transfer velocity vector field calculated by the BEM model to the corresponding fluid particles, set the water density, sound velocity coefficient, density diffusion coefficient and water kinematic viscosity required for SPH simulation, and select the smoothing function type and smoothing length;

[0098] Among them, the velocity initialization of the SPH fluid particles is completed with the help of the particle number Idp, and the transfer of the BEM calculation results to the SPH calculation model is realized.

[0099] Step 7: Start the SPH model, calculate the particle density, pressure gradient force, and viscous force, obtain the velocity increment and density increment of all fluid particles, update the particle density and velocity according to the prediction-correction time format, complete the iteration of the particle position, and enter the next time step solution until the set calculation time is reached.

[0100] The Tait-type state equation is used to calculate the particle pressure from the particle's own density and internal energy using the nonlinear relationship between density and pressure of a compressible liquid under isothermal conditions:

[0101]

[0102] Among them, γ=7, ρ0 is the standard density, c0 is the speed of sound corresponding to the standard density, h swl is the maximum still water level, coef sound is the speed of sound coefficient.

[0103] The SPH format continuity equation used contains a density diffusion term represented by the total density to improve the pressure field distribution near the wall boundary. The SPH fluid particle density increment is calculated as follows:

[0104]

[0105] Among them, ρ i , p i , c i , v i , r i is the density, pressure, sound speed, velocity vector and position vector of particle i; m j ,ρ j , p j , v j , r j is the mass, density, pressure, velocity vector and position vector of particle j; v ij , r ij , z ijare the velocity vector difference, position vector difference, total density difference, static density difference, hydrostatic pressure difference, and vertical distance difference between particles i and j, respectively. N is the total number of particles in the support domain of particle i, and W ij is the smooth function of the influence of particle j on particle i, δ is the density diffusion coefficient, h is the smooth length, ψ ij is the density diffusion term.

[0106] The Navier-Stokes equation discretized using the SPH format contains pressure terms, laminar viscous stress terms, and Reynolds stress terms. The velocity increment of SPH fluid particle i is calculated as follows:

[0107]

[0108] Among them, υ0 is the kinematic viscosity of water, η 2 =0.01h 2 , τ ij is the sub-grid stress tensor.

[0109] After obtaining the velocity increment and density increment of all fluid particles, the density and velocity of the particles are updated according to the prediction-correction time format, and the iteration of the particle position is completed, and the solution of the next time step is entered until the set calculation time is reached. The calculation step size used is:

[0110] Δt=CFL·min(Δt f ,Δt cv )

[0111]

[0112] Among them, f i is the magnitude of the force acting on unit mass (i.e., acceleration), 1≤i,j≤N sph , N sph is the total number of particles in the SPH domain, and CFL is the safety factor.

[0113] The technical solution of the present invention is now illustrated by taking the problem of a wedge-shaped body slamming into water at a uniform velocity vertically as an example. The simulation includes the following steps:

[0114] Step 1: At the initial moment, establish the following in the extended coordinate system: Figure 1 The numerical model shown (in the extended coordinate system, the apex of the wedge Direction coordinates The coordinates of the wedge's apex do not represent its physical height or immersion depth. Initial values ​​for the required parameters are assigned to the model. In the simulation, the wedge's base rise angle β = 45°, height h = 1 m, bottom width l = 2h = 2 m, and entry velocity W = 10 m / s. The water depth d = 10 m, the tank's bottom width L = 2d = 20 m, and the left and right wall heights H = d = 10 m are used. The object surface grid count n1 = 40, the free surface grid count n2 = 120, the side wall grid count n3 = 50, and the bottom grid count n4 = 20.

[0115] Step 2: After the BEM calculation model is initialized, the response calculation begins, and the simulation duration is 0 <t<t c , t c is the specified coupling time step. In this example, t c = h / W = 0.1s, at which point the wedge's depth into the water is equal to its physical height, the initial slamming phase has been completed, and the free surface has not yet curled, deformed, and splashed, and the post-slamming fluid dynamic evolution phase has not yet begun. Four pressure measuring points are arranged on the wedge surface (such as Figure 2 As shown), BEM stage (ie 0 <t<t c ) The time history curve of the impact pressure coefficient at each measuring point is shown in the attached figure. Figure 3 As shown, the calculation formula of the slamming pressure coefficient is C p =P / (ρW 2 ), where P is the pressure value at the measuring point, unit is Pa, and the water density is ρ = 1000 kg / m 3 As the water depth increases, pressure responses appear at each measuring point one by one, and all the measuring points reach the stable section of slamming pressure before the coupling time step.

[0116] Step 3: At the coupling time step (t = t c =0.1s), considering the contribution of the boundary element, calculate the influence coefficient c ij (1≤i,j≤380). In the fully nonlinear theoretical model solved by BEM, the number of boundary grids will increase as the free liquid surface rises during the water entry process, t=t c When n1+n2+n3+n4=380. The influence coefficient c at the coupling time step in this example is ij , velocity potential and solving linear algebraic equations The obtained source strength σ j The calculation results are as follows Figure 4 shown.

[0117] Step 4: At the coupling time step (t = t c = 0.1s), s = 1, the point coordinate values ​​are consistent in the extended coordinate system and the physical coordinate system, and the free liquid surface is cut off from the jet reversal point (such as Figure 5(as shown), obtain the truncated free surface elevation calculation data of the BEM model at the current time step. Taking the truncated free surface considering the influence of the liquid surface rise during the water entry process as the geometric boundary, the wedge boundary and the water tank boundary are consistent with the BEM model. Among them, the wedge boundary moves to the position where the water entry depth s = 1m. Set the particle spacing Δx = 0.05m, and initialize the SPH particle spatial distribution model. The specific results are as Figure 6 shown. In this example, the number of fluid particles n f = 79297, the number of wedge particles n w = 441, and the number of water tank boundary particles n b = 799.

[0118] Step Five: At the coupled time step (t = t c = 0.1s), set the SPH fluid particles as the control points P i (1 ≤ i ≤ 79297), and calculate the unit induced velocity of the distributed sources on the boundary surface element j (1 ≤ j ≤ 380) to the control point P i . By traversing and multiplying with the source strength σ j , obtain the induced velocity at the control point P i . Furthermore, obtain the fluid particle velocity transferred by the BEM model In this example, the vector diagram of the fluid particle velocity transferred at the coupled time step is as Figure 7 shown. Due to the impact effect during the water entry of the wedge, the transfer velocity is generally greater in the near-wedge region than in the near-control surface region.

[0119] Step Six: Taking the particle number Idp as the medium, transfer the fluid particle velocity calculated by the BEM model to the SPH numerical model, and assign the values of the required physical quantities to the SPH calculation model. In this example, ρ0 = 1000 kg / m 3 , coef sound = 20, δ = 0.1, υ0 = 1.0×10 -6 m 2 / s, select the Wendland-type kernel function, and the smoothing length

[0120] Step Seven: Start the SPH numerical model, calculate the particle density, pressure gradient force and viscous force, calculate the velocity change rate and update the particle information, and step forward according to the prediction-correction time format until the set simulation duration T = 0.5s and the time step coefficient CFL = 0.3. The particle evolution results at typical moments in the SPH stage (t c < t < T) are as Figure 8 shown. Attached Figure 8 (a) The distribution of the velocity cloud diagram is the same as that of Attached Figure 7As the wedge further enters the water, the free surface begins to separate and gradually curls and deforms, realizing the simulation of the flow field evolution after the wedge is fully immersed.

[0121] Through the above steps, the entire structural dynamic evolution of a wedge-shaped object during uniform water entry, from impact with the liquid surface to complete immersion, was calculated. The numerical results of this example demonstrate the rationality and effectiveness of the simulation. The proposed BEM and SPH coupled numerical model effectively simulates the entire wedge-shaped object entry process.

[0122] While the embodiments of the present invention have been shown and described, it will be apparent to those skilled in the art that various changes, modifications, substitutions, and alterations can be made to the embodiments without departing from the principles and spirit of the invention.

Claims

1. A BEM-SPH coupled numerical method for wedge-shaped body slamming into water based on the concept of overlapping computational domains, including: S1. Initialize the water impact BEM calculation model under full nonlinear boundary conditions in the extended coordinate system, set the initial state parameters of the BEM model, and set the BEM stage simulation duration according to the physical height of the wedge; S2. At each computational time step, the boundary integral equations constructed with fully nonlinear boundary conditions are solved to obtain the velocity potential at the object surface boundary and the normal derivative of the velocity potential at the free surface. The hybrid Euler-Lagrangian method is used to update the liquid surface lift and velocity potential changes during the water entry process. With the help of auxiliary functions, the slamming pressure on the object surface is calculated. The slamming force is obtained by surface integral and substituted into the Newtonian equation of motion to obtain the motion response of the wedge. The Runge-Kutta method is used to complete time stepping, realizing the time-domain numerical simulation of the wedge entry stage of the BEM. S3, coupled time step, calculate the source strength of each boundary element, limit the integration domain of the influence coefficient to [-1, 1] with the help of shape function, and use the six-node Gauss-Legendre formula for numerical integration; S4, coupled time step, convert the BEM boundary geometric coordinates to the physical coordinate system, with the free liquid surface elevation and the current position of the wedge as the geometric boundary. The SPH tank boundary is the same as the BEM model, and the internal space is filled with SPH fluid particles to form an SPH particle space distribution geometric model; S5, coupled time step, using the SPH fluid particles converted back to the extended coordinate system as control points, calculating the induced velocity of each boundary element on each control point, and further combining the wedge entry velocity to obtain the transfer velocity vector field from the BEM stage to the SPH stage; S6. Initialize the SPH calculation model, assign the transfer velocity vector field calculated by the BEM model to the corresponding fluid particles, set the water density, sound velocity coefficient, density diffusion coefficient and water kinematic viscosity required for SPH simulation, and select the smoothing function type and smoothing length; S7. Start the SPH model, calculate the particle density, pressure gradient force, and viscous force, obtain the velocity increment and density increment of all fluid particles, update the particle density and velocity according to the prediction-correction time format, complete the iteration of the particle position, and enter the solution of the next time step until the set calculation time is reached.

2. The BEM-SPH coupled numerical method for wedge-shaped body water slamming based on the overlapping computational domain concept according to claim 1 is characterized by: The geometric dimensions and velocity variables are established in the extended coordinate system to ensure that the BEM simulation maintains a reasonable and stable grid scale from the initial water contact to the liquid surface rise. The conversion relationship between the relevant physical quantities between the extended coordinate system and the physical coordinate system is: Among them, x, z and v are the coordinate values ​​and velocities of the BEM model in the physical coordinate system. and To expand the grid node coordinates and velocities in the coordinate system, s is the vertical distance of the wedge entering the water, s = ∫Wdt, and W is the water entry velocity.

3. The BEM-SPH coupled numerical method for wedge-shaped body water slamming based on the overlapping computational domain concept according to claim 1 is characterized by: Introducing a local coordinate system at the free surface The liquid level is updated during the BEM model calculation process, where: The axis is on the tangent plane of the free surface, Axis perpendicular to Axis, pointing outside the domain, local coordinate system The conversion relationship with the Cartesian coordinate system is as follows: in, and for The components of the unit vector of the axis in the Cartesian coordinate system, and for The components of the unit vector of the axis in the Cartesian coordinate system; in the local coordinate system, the node motion is updated using material derivatives, and the free surface position and velocity potential are updated using the following formulas: in, is the velocity potential, is the wave surface elevation.

4. The BEM-SPH coupled numerical method for wedge-shaped body water slamming based on the overlapping computational domain concept according to claim 1 is characterized by: In said S2, As a harmonic function, the auxiliary function method is used to solve it to ensure the accuracy of the impact force calculation during the dynamic change of the grid. Split into: Among them, the auxiliary functions χ1 and χ2 satisfy the Laplace equation in the watershed and meet the following boundary conditions: Where n is the normal direction of the boundary element mesh pointing to the fluid domain, It is along The axis weight, Yes and is the velocity potential along Axis and The weight of the axis; Through Green's third formula, the boundary condition of the auxiliary function is converted into a boundary integral equation. After grid discretization, numerical integration, collocation transformation and LU decomposition, the solutions of χ1 and χ2 are obtained. Then, the vertical slamming force on the wedge is obtained. The specific formula is: Among them, F z is the vertical slamming force, l B is the boundary area of ​​the object surface, p is the impact pressure on the object surface, and ρ is the water density.

5. The BEM-SPH coupled numerical method for wedge-shaped body water slamming based on the overlapping computational domain concept according to claim 1 is characterized by: In S3, the linear algebraic equation is solved on the boundary surface element. Calculate the surface element source strength σ j , where the influence coefficient The isoparametric method is used to solve c using a two-node bar element. ij , the specific calculation formula is: Among them, h k is the shape function. The specific form of the one-dimensional shape function is J j (ξ n ) is the Jacobian matrix, ξ n is the coordinate of the Gaussian point, w n is the corresponding integral weight, is the horizontal coordinate value of the source point q, is the vertical coordinate value of the source point q, is the horizontal coordinate value of the node j of the unit, is the vertical coordinate value of the node of unit j, k is the node number of unit j, is the field point p i The horizontal coordinate value of is the field point p i The vertical coordinate value of The source point q and the field point p i The distance between them.

6. The BEM-SPH coupled numerical method for wedge-shaped body water slamming based on the overlapping computational domain concept according to claim 1 is characterized by: In the above S4, the BEM free surface boundary is cut off from the jet reversal point, and the physical coordinates of the cut off free surface boundary are used as the upper boundary of the SPH fluid domain, and the internal fluid particles are filled.

7. The BEM-SPH coupled numerical method for wedge-shaped body water slamming based on the overlapping computational domain concept according to claim 1 is characterized by: In S5, by multiplying the surface element source strength σ j With each boundary grid j control point P i Unit induced speed Get the induced velocity at the control point Unit induced speed The specific calculation formula is: in, is the control point P i The horizontal coordinate value of is the control point P i The vertical coordinate value of The source point q and the control point P i The distance between them.

8. The BEM-SPH coupled numerical method for wedge-shaped body water slamming based on the overlapping computational domain concept according to claim 1 is characterized by: In S6, the velocity initialization of the SPH fluid particles is completed with the help of the particle number Idp, so as to realize the transfer of the BEM calculation results to the SPH calculation model.

9. The BEM-SPH coupled numerical method for wedge-shaped body water slamming based on the overlapping computational domain concept according to claim 1 is characterized by: In S7, the SPH format continuity equation used contains a density diffusion term represented by the total density to improve the pressure field distribution near the wall boundary. The SPH fluid particle density increment is calculated as follows: Among them, ρ i , p i , c i , v i , r i is the density, pressure, sound speed, velocity vector and position vector of particle i; m j ,ρ j , p j , v j , r j is the mass, density, pressure, velocity vector and position vector of particle j; v ij , r ij , zij are the velocity vector difference, position vector difference, total density difference, static density difference, hydrostatic pressure difference and vertical distance difference between particles i and j respectively; N is the total number of particles in the support domain of particle i, W ij is the smooth function of the influence of particle j on particle i, δ is the density diffusion coefficient, ψ ij is the density diffusion term, h is the smoothing length, γ=7, ρ0 is the standard density, c0 is the speed of sound corresponding to the standard density, h swl is the maximum still water level, coef sound is the speed of sound coefficient.

10. The BEM-SPH coupled numerical method for wedge-shaped body water slamming based on the overlapping computational domain concept according to claim 1 is characterized by: In S7, the SPH format momentum equation used contains pressure terms, laminar viscous stress terms, and Reynolds stress terms. The SPH fluid particle velocity increment is calculated as follows: Among them, υ0 is the kinematic viscosity of water, η 2 =0.01h 2 , τ ij is the sub-grid stress tensor.

Images