METHOD FOR ANALYZING HYDROELASTIC EFFECT OF A MARINE STRUCTURE

Abstract

Disclosed is a method for analyzing hydroelastic effect of a marine structure based on moving particle semi-implicit method and modal superposition method. The method is based on a fluid dynamic solver of a moving particle semi-implicit method (MPS) and a structural response solver of a mode superposition method, and the fluid dynamic solver and the structural response solver achieve strong coupling in a time domain in an iterative mode. According to the present disclosure, the deficiency of the traditional potential flow-based hydroelesticity calculation (i.e. cannot directly simulate the severe deformation of the free surface and the slamming phenomenon of the ship body) is improved, and the problem of mutual influence between the rigidity of the structure and the flexible motion modal is considered.

Description

BACKGROUND

As ships and marine structures become increasingly larger and faster, the natural frequency of their structures and the encountered wave frequency tend to be in the same range. This means that the hydroelastic effect is becoming more and more important for the overall fatigue analysis of marine structures.

The analysis of structural integrity in this case requires time-domain simulation. Since the traditional potential-flow-based solvers cannot directly simulate the violent deformation of the free surface during slamming, it requires a CFD-based solver to perform flow calculations. Partial combination or full application of the CFD model based on viscous flow is the ultimate goal of hydroelastic calculation.

Compared with the traditional grid-based CFD methods, particle methods such as SPH and MPS do not have grid distortion problems, and the flow field is updated by the Lagrangian method, which is very suitable for simulating large free surface deformation. Moving Particle Semi-implicit (MPS), as one of the mainstream particle methods, has been improved and applied to various severe fluid and structural interaction problems. In addition, the typical motion of the ship is large rigid body motion plus relatively small elastic deformation. This means that the modal superposition theory is sufficient to accurately calculate the elastic deformation. However, in the traditional hydroelastic analysis, the interaction between the rigid body and the elastic modal is not considered, and the rigid body and the elastic part are only related to each other in the sense of the hydrodynamic coefficient of the cross term.

SUMMARY

According to the technical problems mentioned above, a method for analyzing hydroelastic effect of a marine structure based on moving particle semi-implicit method and modal superposition method is provided. The marine structure can be a ship, a semi-submersible platform, or models thereof. The technical means adopted by the present invention are as follows:

Disclosed is a method for solving a strong nonlinear time domain hydroelasticity problem based on moving particle semi-implicit method and modal superposition method, and in the time domain, the solutions are decomposed into a number of time intervals with small intervals. Assuming that the information of the flow field and structural motion in the previous time step (or initial moment) is known, the following calculation is performed on the current time step. The specific steps are as follows:

S1. spatially discretizing a boundary of the flow field with mesh-less discrete points to avoid the difficulty of re-meshing caused by strong nonlinear free surface motion and large rigid body displacement of the structure, wherein the boundary is a surface of a marine structure, e.g., a ship, a semi-submersible platform, or models thereof, wherein if a previous time step is initial time, velocity and intensity of pressure of fluid particles need to be initialized;

S2. determining an estimated value of a structural boundary position of a current time step according to motion information of a solid boundary of the previous time step, wherein if the previous time step is the initial time, an initial value of a structural motion is adopted;

S3. updating a position and a velocity of fluid particles by Navier-Stokes equation according to the flow field information of the previous time step, excluding fluid pressure;

S4. deriving a pressure Poisson equation and solving the pressure Poisson equation;

S5. updating the position and the velocity of the fluid particles using the pressure value obtained in S4;

S6. solving the governing equation of solid response which couples the rigid body and elastic mode using the pressure value obtained in step S4;

S7. updating the structure motion and deformation by using dynamic information obtained in step S6, which in turn provides a new fluid boundary position;

S8. comparing structural position differences in step S2 and step S7, wherein if a convergence condition is satisfied, a next time step calculation is performed, otherwise, repeating steps S2 to S8 at a current time step until the next time step calculation is performed after the convergence condition is satisfied.

Specific steps of the step S1 are as follows:

In the MPS method, the gradient, divergence and Laplace operator of a physical quantity use the weighted average method, discretizing based on uniformly distributed particles, as shown in equation (1.1.1).

$\begin{array}{cc}\nabla \phi \left(r_i\right)=\frac{d_m}{n_0}\sum _{j\ne i}^N\frac{\phi \left(r_j\right)-\phi \left(r_i\right)}{r_{\mathrm{ij}}^2}\left(r_j-r_i\right)W\left(r_{\mathrm{ij}}\right)& \left(1\right)\\ {\nabla }^2\phi \left(r_i\right)=\frac{2d_m}{n_0\lambda }\sum _{j\ne i}^N\left[\phi \left(r_j\right)-\phi \left(r_i\right)\right]W\left(r_{\mathrm{ij}}\right)& \left(2\right)\left(1.1\mathrm{.1}\right)\\ \nabla ·\Phi \left(r_i\right)=\frac{d_m}{n_0}\sum _{j\ne i}^N\frac{\left(\Phi \left(r_j\right)-\Phi \left(r_i\right)\right)·\left(r_j-r_i\right)}{r_{\mathrm{ij}}^2}W\left(r_{\mathrm{ij}}\right)& \left(3\right)\end{array}$

wherein, ϕ and Φ represent arbitrary scalars and vectors respectively, d_m represents a number of dimensions of the problem under study: two dimensions or three dimensions, N is a number of particles in an affected area of the relevant particles; n₀ is a particle density at initial time, r represents a particle spacing, and r_i and r_j represent the coordinates of particle i and particle j respectively r_ij=|r_i−r_j|.

A kernel function W (r_ij) and a parameter λ are defined as equation (1.1.2) and equation (1.1.3):

$\begin{array}{cc}W\left(r_{\mathrm{ij}}\right)=\{\begin{array}{cc}\frac{r_e}{r_{\mathrm{ij}}}-1& 0\le r_{\mathrm{ij}}\le r_e\\ 0& r_{\mathrm{ij}}\ge r_e\end{array}& \left(1.1\mathrm{.2}\right)\\ \lambda =\frac{\sum _{j\ne i}^NW\left(r_{\mathrm{ij}}\right)r_{\mathrm{ij}}^2}{\sum _{j\ne i}^NW\left(r_{\mathrm{ij}}\right)}& \left(1.1\mathrm{.3}\right)\end{array}$

wherein, r_e represents an effective range of particles interaction, and the particle density n is defined in equation (1.1.4):

$\begin{array}{cc}n=\sum _{j\ne i}^NW\left(r_{\mathrm{ij}}\right)& \left(1.1\mathrm{.4}\right)\end{array}$

In order to keep the particle distribution more regular and improve the stability and accuracy of the calculation, the technology of particle displacement shifting and collision processing (mainly for free surface particles) are adopted respectively. In this technique, the displacement of the particles is adjusted directly, rather than by applying a force. After each time step, the particle position is slightly displaced to make its distribution regular. This technique can also be considered as a grid reconstruction process. In addition, since the pressure Poisson equation is derived from the incompressibility condition, the resulting pressure will keep the distance between adjacent particles approximately at the same value (i.e., the initial particle distance r0).

An adjustment amount of the particle spacing as shown in equation (1.1.5):

$\begin{array}{cc}\delta r_i=\sum _{j\ne i}\frac{r_{0-}r_{\mathrm{ij}}}{2}·\frac{r_{i-}r_j}{r_{\mathrm{ij}}},\mathrm{when}r_{\mathrm{ij}}\le r_0& \left(1.1\mathrm{.5}\right)\end{array}$

wherein, r₀ is an initial particle spacing.

For free surface particles distal from a main fluid domain, when a relative velocity between the particles is closer than a threshold before a prediction step S3, performing an operation as shown in equation (1.1.6) to set the relative velocity of the particle pair δu_i to be zero;

$\begin{array}{cc}\delta u_i=\sum _{j\ne i}-\square \left(r_{\mathrm{ij}}\right)u_{{\tau }_{\mathrm{ij}}},\mathrm{when}r_{\mathrm{ij}}-u_{{\tau }_{\mathrm{ij}}}\Delta t\le r_{\min },r_{\min }=0.3r_0& \left(1.1\mathrm{.6}\right)\end{array}$

wherein, u_rij is a relative tangential velocity of particle i and particle j, and the adjacent particles of solid and fluid are taken as 0.5 and 1.0 respectively; Δt is a size of the time step.

Specific steps of the step S2 are as follows: introducing a variable A describing the motion and deformation of an object, which is defined as follows:

A=[X_Rb,θ,q_b]  (2.1.1)

wherein, X_Rb represents a position of a center of mass of a rigid body motion, θ represents a deadrise angle of the structure, and q_b represents general coordinates corresponding to a modal function.

The corresponding position, velocity and acceleration of the particles at an interface of the structure and the fluid are Γ_fsi,0^k, {dot over (Γ)}_fsi,0^k and {umlaut over (Γ)}_fsi,0^k, wherein k represents a time step, a position where a subscript ‘0’ locates represents a number of iterations.

Estimating the position, velocity and acceleration of the current time step from the position, velocity and acceleration (that is, Γ_fsi,0^k−1, {dot over (Γ)}_fsi,0^k−1 and {umlaut over (Γ)}_fsi,0^k−1) of the structure boundary particles of the previous time step according to an assumption of uniform acceleration motion, and calculating an estimated value of the structural boundary position of the current time step based on the calculated position, velocity and acceleration of the current time step.

In the step 3, the Navier-Stokes equation is an incompressible and non-viscous Navier-Stokes equation based on a Lagrangian form, which is shown in equation (3.1.1):

$\begin{array}{cc}\frac{\mathrm{Du}}{\mathrm{Dt}}=g-\frac{\nabla p}{\rho }& \left(3.1\mathrm{.1}\right)\end{array}$

wherein, u, p and ρ represent the velocity, pressure and density of the fluid respectively, g is the vector pointing in a direction of gravity, and g represents the acceleration of gravity.

Specific steps of the step S3 are as follows:

Using a classic projection method to decouple the speed and pressure as follows:

Firstly, without considering the pressure, updating the liquid to an intermediate state only by inertia, as shown in equation (3.1.2):

u _* =u _k +Δtg  (1)

I _* =I _k +Δtu _*  (2)(3.1.2)

wherein, l is a position vector, variables with subscripts * and k represent variables at an intermediate time step and a k^th time step respectively, and u represents the velocity of the fluid.

Specific steps of the step S4 are as follows:

Deriving the pressure Poisson equation by taking a divergence of the Navier-Stokes equation on both sides according to a continuity equation, as shown in equation (4.1.2), wherein, the continuity equation is shown in equation (4.1.1),

$\begin{array}{cc}\nabla ·u=0& \left(4.1\mathrm{.1}\right)\\ {\nabla }^2p_{k+1}=\frac{\rho \nabla ·u_{*}}{\Delta t}+\mathrm{\alpha \rho }\frac{n_0-n_k}{n_0\Delta t^2}& \left(4.1\mathrm{.2}\right)\end{array}$

wherein, n₀ and n_k are the particle densities at the initial time and the k^th time step respectively, and the particle density is defined as shown in equation (1.1.4).

The coefficient α is defined as shown in (4.1.3):

$\begin{array}{cc}\alpha =\{\begin{array}{cc}\frac{n_0-n_k}{n_0}+\Delta t\nabla ·u_k& \left(n_0-n_k\right)\nabla ·u_k\ge 0\\ \frac{n_0-n_k}{n_0}& \left(n_0-n_k\right)\nabla ·u_k\le 0\end{array}& \left(4.1\mathrm{.3}\right)\end{array}$

The α model does not require special calibration for different problems, and has stronger stability than other similar models.

Boundary conditions for solving the pressure Poisson equation are as follows:

(1) Solid Boundary Condition

Applying Neumann boundary conditions to solid particles, as shown in equation (4.1.4),

n·∇p _k+1=ρ(n·g−n·{dot over (u)}_bk+1)≅β(n·g−n·{dot over (u)}_bk)  (4.1.4)

wherein, {dot over (u)}_bk and {dot over (u)}_bk+1 the accelerations of the solid boundary particles at the k^th and (k+1)^th time steps; since {dot over (u)}_bk+1 is unknown when solving the fluid motion at the (k+1)^th time step, the value {dot over (u)}_bk of the previous time step is used as an approximation.

Fluid particles proximal to the solid boundary need to be corrected for the Laplace operator, making it consistent with the Neumann boundary condition on the solid boundary, i.e. to make up for the deficiencies of adjacent particles. The pressure value for the virtual fluid particles outside the solid boundary is shown in equation (4.1.5); for near boundary fluid particles, the source term of the pressure Poisson equation also includes the role of solid boundary particles in the influence domain; therefore, due to the same consistency reasons, the intermediate velocity of the affected solid boundary particles is corrected as equation (4.1.6)-(1), which then are projected to tangential r and normal n directions respectively to obtain equations (4.1.6)-(2) and (4.1.6)-(3),

$\begin{array}{cc}{p}_v=p_s+{\rho }_0\left(n·g-n·u_b\right)r_0& \left(4.1\mathrm{.5}\right)\\ u_{b^{*}}=u_{\mathrm{bk}+1}+\Delta t\frac{\nabla p_k}{{\rho }_0}& \left(1\right)\\ n·u_{b^{*}}=n·u_{\mathrm{bk}+1}+\Delta t\left(n·g-n·{\stackrel{.}{u}}_{b_k}\right)& \left(2\right)\left(4.1\mathrm{.6}\right)\\ \tau ·u_{b^{*}}=\tau ·u_{\mathrm{bk}+1}+\frac{\Delta t}{{\rho }_0}\frac{\partial p_k}{\partial \tau }& \left(3\right)\end{array}$

wherein, p_v represents an intensity of pressure of virtual particles, p_s represents an intensity of pressure of responsive solid particles, u_b* and u_bk+1 are the velocity of the solid boundary particles at the intermediate time and the (k+1)^th time step, and n in equations (4.1.5) and (4.1.6) represents a normal direction.

(2) Free Surface Boundary Condition

Since solving the pressure Poisson equation requires applying a boundary condition with zero pressure on the free surface boundary, it is necessary to identify a position of the free surface particles:

For the case of two-dimensional cases, each particle is assigned a circle centered on itself with a radius of 1.05r₀, and the circle is discretized with 360 points evenly distributed on it; if all these points are covered by the circles of its neighboring particles, then the particle is an internal fluid particle, otherwise the particle is to be identified as a free surface particle; the inspection process is to identify the overlapping points on the arc (or patch in 3D case) covered by each adjacent particle by traversing each adjacent particle;

For the case of three-dimensional cases, a two-step method is used: firstly, checking a number of adjacent particles using equation (4.1.7) to detect potential free surface particles,

N _p<β_3dN_p0  (4.1.7)

wherein, β_3d is a parameter less than 1, taking values according to situations, and N_p0 represents a initial particle density within a delimited range.

Secondly, refining a search by establishing a spherical surface with a radius of 1.05r₀ centered on the particle, in detail: determining a vector pointing to the most sparse direction of a particle distribution around the particle by weighted average method, and discretizing a circular area on the spherical surface by uniformly distributed points, if all these points are covered by the spherical surface of its neighboring particles, the particle is an internal fluid particle, otherwise the particle will be identified as a free surface particle.

Calculation equation of step S5 is as shown in formula (5.1.1),

$\begin{array}{cc}{u}_{k+1}=u_{*}-\Delta t·\frac{\nabla p_{k+1}}{{\rho }_0}r_{k+1}=r_k+\Delta t·u_{k+1}& \left(5.1\mathrm{.1}\right)\end{array}$

The dynamic response of a typical ocean structure to wave impact can be expressed as a large rigid body motion and a small elastic deformation, which means that the modal superposition method can handle the elastic motion part. The present invention derives the time-domain model of the coupling rigid body and elastic modal based on the Lagrangian equation. The specific steps of step S6 are as follows:

S6.1. Establishment of Structural Motion Description Model

The structural motion description model uses a fixed global coordinate system XOY and an local coordinate system soη, wherein an origin of the attachment local coordinate system is selected at a center of gravity of a beam of the structural motion description model, thus a position X of any point on the beam of the structural motion description model can be described by equation (6.1.1):

X=X _Rb +Rξ ^T  (6.1.1)

wherein, X_Rb is a position of a center of mass of the beam of the structural motion description model, and equations (6.1.2) and (6.1.3) define local coordinates ξ of the beam of the structural motion description model and rotation matrix R of the beam of the structural motion description model,

$\begin{array}{cc}\xi =\left[s,{\eta }_l+{\eta }_0\right]& \left(6.1\mathrm{.2}\right)\\ R=\left[\begin{array}{cc}\cos \left(\theta \right)& -\sin \left(\theta \right)\\ \sin \left(\theta \right)& \cos \left(\theta \right)\end{array}\right]& \left(6.1\mathrm{.3}\right)\end{array}$

wherein, θ is an angle between O-X axis and o-s axis counterclockwise, η₁ is a coordinate of an undeformed structure, which is a constant for a specific point on the beam of the structural motion description model, and η₀ is a deflection of a mid-plane, the deflection of any plane parallel to the mid-plane taking the same value of the mid-plane; according to a modal superposition theory, the deflection η₀ of the mid-plane is expanded into equation (6.1.4):

η₀=φ^Tq_b  (6.1.4)

Equations (6.1.5) and (6.1.6) define a column vector φ of a modal function and corresponding general coordinates q_b,

φ=[φ₁,φ₂,φ₃, . . . ]^T  (6.1.5)

q _b=[q_b1,q_b2,q_b3, . . . ]^T  (6.1.6)

The modal function satisfies the following orthogonal relationship:

$\begin{array}{cc}\underset{s_1}{\overset{s_2}{\int }}{\mathrm{\varphi \rho }}_l{\varphi }^T\mathrm{ds}=I_u,I_u\mathrm{is}\mathrm{an}\mathrm{identity}\mathrm{matrix}& \left(6.1\mathrm{.7}\right)\\ \underset{s_1}{\overset{s_2}{\int }}\frac{d^2\varphi }{{\mathrm{ds}}^2}{\mathrm{EI}}_I\frac{d^2{\varphi }^T}{{\mathrm{ds}}^2}\mathrm{ds}=\Lambda ,\Lambda =\mathrm{diag}\left({\omega }_i^2\right)i=1,2,3,\dots & \left(6.1\mathrm{.8}\right)\end{array}$

wherein s₁ and s₂ are upper and lower limits integrated along o-s axis, ρ₁ is an line density of the beam of the structural motion description model, E represents an elastic modulus of a structure, and I_I is an inertia matrix of the structure;

S6.2. Establishment of the Governing Equation of Solid Response Based on Lagrange Equation

The main feature of the modal superposition method model is that the rigid body modal and the elastic modal are coupled in the time domain. The effects of interactions such as inertial effects of rigid body motion on elastic deformation are considered. This is an extension of the traditional hydroelastic theory. In the traditional hydroelastic theory, rigid body and elastic part are only connected by the corresponding cross-term hydrodynamic coefficient.

Establishing a model based on Lagrange equation:

$\begin{array}{cc}\frac{d}{\mathrm{dt}}\left(\frac{\partial T}{\partial {\stackrel{.}{q}}_j}\right)+\frac{\partial U}{\partial q_j}-\frac{\partial T}{\partial q_j}=Q_jj=1,2,3,\dots & \left(6.2\mathrm{.1}\right)\end{array}$

wherein T and U represent an kinetic energy and potential energy of the structural motion description model, q_j is a general coordinate corresponding to any rigid-flexible mode, and Q_j is a non-conservative force corresponding to the j^th coordinate.

The beam of the structural motion description model is an elastic non-uniform beam, and specific forms of T, U, q_j, Q_j of the elastic non-uniform beam are as follows:

$\begin{array}{cc}\begin{array}{cc}T=& \frac{1}{2}\underset{s_1}{\overset{s_2}{\int }}\underset{{\eta }_1}{\overset{{\eta }_2}{\int }}{\stackrel{.}{X}}^T{\rho }_s\stackrel{.}{X}d\eta \mathrm{ds}=\\ & \frac{1}{2}\underset{s_1}{\overset{s_2}{\int }}\underset{{\eta }_1}{\overset{{\eta }_2}{\int }}{\left({\stackrel{.}{X}}_R+\mathrm{RU}\stackrel{.}{\theta }\xi +R\stackrel{.}{\xi }\right)}^T{\rho }_s\left({\stackrel{.}{X}}_R+\mathrm{RU}\stackrel{.}{\theta }\xi +R\stackrel{.}{\xi }\right)d\eta \mathrm{ds}\\ =& \frac{1}{2}[M\left({\stackrel{.}{X}}_R^2+{\stackrel{.}{Y}}_R^2\right)+{\stackrel{.}{\theta }}^2I_b+{\stackrel{.}{\theta }}^2q_b^Tq_b+{\stackrel{.}{q}}_b^T{\stackrel{.}{q}}_b-\\ & 2\stackrel{.}{\theta }\left({\stackrel{.}{X}}_R\cos \left(\theta \right)+{\stackrel{.}{Y}}_R\sin \left(\theta \right)\right)q_b^T{\psi }_0+\\ & 2\left(-{\stackrel{.}{X}}_R\sin \left(\theta \right)+{\stackrel{.}{Y}}_R\cos \left(\theta \right)\right){\stackrel{.}{q}}_b^T{\psi }_0+2\stackrel{.}{\theta }{\stackrel{.}{q}}_b^T{\psi }_1]\end{array}& \left(6.2\mathrm{.2}\right)\\ \begin{array}{cc}U=& \frac{1}{2}\underset{s_1}{\overset{s_2}{\int }}\frac{d^2{\eta }_0}{{\mathrm{ds}}^2}\mathrm{EI}\frac{d^2{\eta }_0}{{\mathrm{ds}}^2}\mathrm{ds}+{\mathrm{MgY}}_R\\ =& \frac{1}{2}q_b^T\left(\underset{s_1}{\overset{s_2}{\int }}\frac{d^2\varphi }{{\mathrm{ds}}^2}\mathrm{EI}\frac{d^2{\varphi }^T}{{\mathrm{ds}}^2}\mathrm{ds}\right)q_b+\\ & {\mathrm{MgY}}_R=\frac{1}{2}q_b^T\Lambda q_b+{\mathrm{MgY}}_R\end{array}& \left(6.2\mathrm{.3}\right)\end{array}$

wherein ρ_s represents a linear density, M represents mass of a structure, X_R and Y_R correspond to X and Y coordinates of a center of mass of the structure respectively, g represents the acceleration of gravity, η₁ and η₂ are the upper and lower limits of integral along axis η, matrix integral U and column vector integrals Ψ₀ and ψ₁, the matrix integral U and column vector integrals Ψ₀, Ψ₁ defined in equations (6.2.4) and (6.2.6):

$\begin{array}{cc}U=\left[\begin{array}{cc}0& -1\\ 1& 0\end{array}\right]& \left(6.2\mathrm{.4}\right)\\ {\psi }_0={[{\psi }_{01},{\psi }_{02},{\psi }_{03},\dots ]}^T=\underset{s_1}{\overset{s_2}{\int }}{\mathrm{\varphi \rho }}_l\mathrm{ds}& \left(6.2\mathrm{.5}\right)\\ {\psi }_1={[{\psi }_{11},{\psi }_{12},{\psi }_{13},\dots ]}^T=\underset{s_1}{\overset{s_2}{\int }}s{\mathrm{\varphi \rho }}_l\mathrm{ds}& \left(6.2\mathrm{.6}\right)\end{array}$

The governing equations of the dynamic response for elastic non-uniform beam under a wave impact can be obtained by substituting equations (6.2.2) and (6.2.3) into equation (6.2.1), which is shown in equation (6.2.7):

M{umlaut over (X)} _Rb+{dot over (θ)}² sin(θ)ψ₀^Tq_b−2{dot over (θ)} cos(θ)Ψ₀^T{dot over (q)}_b−{umlaut over (θ)} cos(θ)Ψ₀^Tq_b−sin(θ)Ψ₀^T{umlaut over (q)}_b=Q_XRb

MŸ _Rb−{dot over (θ)}² cos(θ)ψ₀^Tq_b−2{dot over (θ)} sin(θ)ψ₀^T{dot over (q)}_b−{umlaut over (θ)} sin(θ)ψ₀^Tq_b+cos(θ)ψ₀^T{umlaut over (q)}_b+Mg=Q_YRb

−({umlaut over (X)}_Rb cos(θ)+Ÿ_Rb sin(θ))ψ₀^Tq_b+I_b{umlaut over (θ)}+{umlaut over (θ)}q_b^Tq_b+2{dot over (θ)}{dot over (q)}_b^Tq_b+ψ₁^T{umlaut over (q)}_b=Q_θ

(−{umlaut over (X)}_Rb sin(θ)+Ÿ_Rb cos(θ))ψ₀+{umlaut over (θ)}ψ₁+{umlaut over (q)}_b+Λq_b−{dot over (θ)}²q_b=Q_qb  (6.2.7)

Equations (6.2.8)-(6.2.11) provide the non-conservative forces corresponding to the general coordinates of the rigid body and the elastic modal:

Q _{X Rb} = pn _s dl  (6.2.8)

Q _{Y Rb} = pn _η dl  (6.2.9)

Q _θ = p(X_pn_p−Y_pn_s)dl  (6.2.10)

Q _{q b} = p(n·e_η)φdl  (6.2.11)

wherein, n_s and n_η are local normal components of s and η respectively, expressed by a fixed global coordinate system, X_p and Y_p are overall X and Y coordinates of a point of action of an intensity of pressure p, e_η is an unit vector of O-η axis, and the integration is carried out around a perimeter of the beam of the structural motion description model; and it can be seen from equation (6.2.7) that the interaction between the rigid body and the elastic mode, i.e. the cross term, is caused by angular motion.

S6.3. Solving the Governing Equation of Solid Response which Couples the Rigid Body and Elastic Mode.

For the interaction process of fluid and structure, standard iteration method is used until certain criteria are met.

Specific steps of the step S8 are as follows:

Assuming that all fluid and structural variables are known in time step t=t_k−1, and then the detailed interaction process at the next time t=t_k is as follows:

S8.1. assuming that an acceleration of a structure at the time step t=t_k is the same as the acceleration at the previous time step t=t_k−1, that is, {umlaut over (D)}^k={umlaut over (D)}^k−1, a position D^k and speed {dot over (D)}^k at the time step t=t_k can also be calculated by a finite difference method.

Calculating the corresponding position, velocity and acceleration of each point on the interface between fluid and structure through equations (8.1.1)˜(8.1.6), that is, Γ_fsi,0^(k), {dot over (Γ)}_fsi,0^(k) and {umlaut over (Γ)}_fsi,0^(k) can be calculated from the newly calculated structural motion and deformation parameters A, that is, equation (2.1.1),

$\begin{array}{cc}X_R=X_{\mathrm{CR}}+R_{{\mathrm{RX}}_R}& \left(8.1\mathrm{.1}\right)\\ X_{\mathrm{fi}}=X_{\mathrm{cfi}}+R_{\mathrm{fi}}{\xi }_i^T& \left(8.1\mathrm{.2}\right)\\ {\stackrel{.}{X}}_R={\stackrel{.}{X}}_{\mathrm{CR}}+{\stackrel{.}{R}}_{{\mathrm{RX}}_R}={\stackrel{.}{X}}_CR+{\stackrel{.}{R}}_RU{\stackrel{.}{\theta }}_R{\chi }_R& \left(8.1\mathrm{.3}\right)\\ {\stackrel{.}{X}}_{\mathrm{fi}}={\stackrel{.}{X}}_{\mathrm{cR}}+{\stackrel{.}{R}}_R{\chi }_{\mathrm{ofi}}+{\stackrel{.}{R}}_{\mathrm{fi}}{\xi }_i^T+{\stackrel{.}{R}}_{\mathrm{fi}}{\xi }_i^T={\stackrel{.}{X}}_{\mathrm{CR}}+R_RU{\stackrel{.}{\theta }}_R{\chi }_{\mathrm{ofi}}+R_{\mathrm{fi}}U{\stackrel{.}{\theta }}_{\mathrm{fi}}{\chi }_{\mathrm{ofi}}{\xi }_i+R_{\mathrm{fi}}{\stackrel{.}{\xi }}_i& \left(8.1\mathrm{.4}\right)\\ {\ddot{X}}_R={\ddot{X}}_{\mathrm{CR}}+{\ddot{R}}_R{\chi }_R={\ddot{X}}_R+\left(R_RU{\ddot{\theta }}_R-R_R{\ddot{\theta }}_R^2\right){\chi }_R& \left(8.1\mathrm{.5}\right)\\ \begin{array}{cc}{\ddot{X}}_{\mathrm{fi}}=& {\ddot{X}}_{\mathrm{CR}}+{\ddot{R}}_R{\chi }_{\mathrm{ofi}}+{\ddot{R}}_{\mathrm{fi}}{\xi }_i+2{\stackrel{.}{R}}_{\mathrm{fi}}{\stackrel{.}{\xi }}_i+R_{\mathrm{fi}}{\ddot{\xi }}_i\\ =& {\ddot{X}}_{\mathrm{CR}}+\left(R_RU{\ddot{\theta }}_R-R_R{\ddot{\theta }}_R^2\right){\chi }_{\mathrm{ofi}}+\\ & \left(R_{\mathrm{fi}}U{\ddot{\theta }}_{\mathrm{fi}}-R_{\mathrm{fi}}{\ddot{\theta }}_{\mathrm{fi}}^2\right){\xi }_i+2R_{\mathrm{fi}}U{\ddot{\theta }}_{\mathrm{fi}}{\stackrel{.}{\xi }}_i+R_{\mathrm{fi}}{\ddot{\xi }}_i\end{array}& \left(8.1\mathrm{.6}\right)\end{array}$

wherein, X_R is global coordinates of rigid body, X_CR is a center-of-mass coordinates, R_R is a rotation matrix that associates the local coordinate system soη with the fixed global coordinate system XOY; X_fi represents the coordinates on a center line; and ξ_i represents the coordinates of the beam center line of in the local coordinate system soη, as represented in equation (8.1.7), wherein η_i is deflection function of the beam;

ξ_i=[s_i,η_i]  (8.1.7)

S8.2. calculating the fluid motion at time t=t_k according to the modified MPS method by using updated information of the interface as a new boundary condition, by using the newly updated pressure p_fsi,i+1^k, the i^th iteration of the interface at time t=t_k can be conducted;

S8.3. comparing structural position differences Γ_fsi,i^k and {tilde over (Γ)}_fsi,i+1^k in step S2 and step S7, if a convergence condition is satisfied, as shown in equation (8.3.1), performing step S8.1 for the next time step t=t_k+1 operation,

$\begin{array}{cc}{\tilde{\Gamma }}_{\mathrm{fsi},i+1}^k-{\Gamma }_{\mathrm{fsi},i}^k\le ϵ& \left(8.3\mathrm{.1}\right)\end{array}$

Otherwise using equation (8.3.2) to correct the structural position D_i+1^k of the (i+1)^th iteration, and using Newmark method to update velocity {dot over (D)}_i+1^k and acceleration {umlaut over (D)}_i+1^k, calculating interface variables Γ_fsi,i+1^k, {dot over (Γ)}_fsi,i+1^k and {umlaut over (Γ)}_fsi,i+1^k, based on D_i+1^k, {dot over (D)}_i+1^k and {umlaut over (D)}_i+1^k; using these modified interface variables to perform the (i+1)^th iteration by returning to step S2,

D _i+1 ^k ={tilde over (D)} _i+1 ^k+(1−χ_i)D_i^k  (8.3.2)

wherein, χ_i is Aitken relaxation factor, a core idea of which is to improve the accuracy of prediction by using the first two iterations.

λ_i is calculated by equation (8.3.3):

$\begin{array}{cc}{\chi }_i=-{\chi }_{i-1}\frac{{\left({\mathrm{\Delta \Gamma }}_{\mathrm{fsi},i}^k\right)}^T\left({\mathrm{\Delta \Gamma }}_{\mathrm{fsi},i+1}^k-{\mathrm{\Delta \Gamma }}_{\mathrm{fsi},i}^k\right)}{{\left({\mathrm{\Delta \Gamma }}_{\mathrm{fsi},i+1}^k-{\mathrm{\Delta \Gamma }}_{\mathrm{fsi},i}^k\right)}^T\left({\mathrm{\Delta \Gamma }}_{\mathrm{fsi},i+1}^k-{\mathrm{\Delta \Gamma }}_{\mathrm{fsi},i}^k\right)}& \left(8.3\mathrm{.3}\right)\end{array}$

wherein, Γ_fsi,i^k={tilde over (Γ)}_fsi,i^k−Γ_fsi,i−1^k.

Disclosed is a method for analyzing hydroelastic effect of a marine structure base on a moving particle semi-implicit method and a modal superposition method. The method is based on a fluid dynamic solver of a moving particle semi-implicit method (MPS) and a structural response solver of a mode superposition method based on a coupling rigid body and an elastic modal, and the fluid dynamic solver and the structural response solver achieve strong coupling in a time domain in an iterative manner. According to the present disclosure, the deficiency that the traditional potential flow-based hydroelesticity calculation cannot directly simulate the severe deformation of the free surface and the slamming phenomenon of the ship body is improved, and the problem of mutual influence between the rigidity of the structure and the flexible motion modal is considered.

BRIEF DESCRIPTION OF THE DRAWINGS

In order to more clearly explain the embodiments of the present invention or the technical solutions in the prior art, the drawings required in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings in the following descriptions are some embodiments of the present invention. For those of ordinary skill in the art, without inventive effort, other drawings can be obtained based on these drawings.

FIG. 1 is two-dimensional and three-dimensional free surface particle recognition schemes.

FIG. 2 is a schematic diagram of an elastic non-uniform beam.

FIG. 3 is a schematic diagram of the initial configuration of the calculation (not to scale).

FIG. 4 is a cross-sectional and surface discretization diagram of a ship model.

FIG. 5 shows the distribution of the mass (a) and second moment (b) of the ship model.

FIG. 6 shows the first three modes of the ship model.

FIG. 7 is a typical time step free surface shape and ship deformation diagram ((a) t=0.1s, (b) t=0.2s, (c) t=0.418s, (d) t=0.684s, (e) t=0.91s, (f) t=1.05s).

FIG. 8 is a demonstration of the Laplace operator used to compensate for virtual particles near the solid boundary.

FIG. 9 is a flow chart of the iterative process of interaction between fluid and structure.

DESCRIPTION OF THE EMBODIMENTS

To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be described clearly and completely in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are parts of the embodiments of the present invention, but not all the embodiments. Based on the embodiments of the present invention, all the other embodiments obtained by those of ordinary skill in the art without inventive effort are within the scope of the present invention.

As shown in FIGS. 1-9, disclosed is a method for analyzing hydroelastic effect of a marine structure based on a moving particle semi-implicit method and a modal superposition method, and in the time domain, the solutions are decomposed into a number of time intervals with small intervals. Assuming that the information of the flow field and structural motion in the previous time step (or initial moment) is known, the following calculation is performed on the current time step, including the following steps:

S1. spatially discretizing the boundary of the flow field with mesh-less discrete points to avoid the difficulty of re-meshing caused by strong nonlinear free surface motion and large rigid body displacement of the structure, wherein if a previous time step is initial time, velocity and intensity of pressure of fluid particles need to be initialized;

S2. determining an estimated value of a structural boundary position of a current time step according to motion information of a solid boundary of the previous time step, wherein if the previous time step is the initial time, an initial value of a structural motion is adopted;

S3. updating a position and a velocity of fluid particles by Navier-Stokes equation according to flow field information of the previous time step, excluding fluid pressure;

S4. deriving a pressure Poisson equation and solving the pressure Poisson equation;

S5. updating the position and the velocity of the fluid particles using the pressure value obtained in S4;

S6. solving the governing equation of solid response using the pressure value obtained in step S4;

S7. updating the structure motion and deformation by using dynamic information obtained in step S6, which in turn provides the new fluid boundary position;

S8. comparing structural position differences in step S2 and step S7, wherein if a convergence condition is satisfied, a next time step calculation is performed, otherwise, repeating steps S2 to S8 at the current time step until the next time step calculation is performed after the convergence condition is satisfied.

Steps S1-S4 are part of the hydrodynamic solver (modified MPS method), and step S6 is a structural response solver.

Specific steps of the step S1 are as follows: discretizing based on uniformly distributed particles, as shown in equation (1.1.1).

$\begin{array}{cc}\nabla \phi =\frac{d_m}{n_0}\sum _{j\ne i}^N\frac{\phi \left(r_j\right)-\phi \left(r_i\right)}{r_{\mathrm{ij}}^2}\left(r_j-r_i\right)W\left(r_{\mathrm{ij}}\right)& \left(1\right)\\ {\nabla }^2\phi \left(r_i\right)=\frac{2d_m}{n_0\lambda }\sum _{j\ne i}^N\left[\phi \left(r_j\right)-\phi \left(r_i\right)\right]W\left(r_{\mathrm{ij}}\right)& \left(2\right)\left(1.1\mathrm{.1}\right)\\ \nabla ·\Phi \left(r_i\right)=\frac{d_m}{n_0}\sum _{j\ne i}^N\frac{\left(\Phi \left(r_j\right)-\Phi \left(r_i\right)\right)·\left(r_j-r_i\right)}{r_{\mathrm{ij}}^2}W\left(r_{\mathrm{ij}}\right)& \left(3\right)\end{array}$

wherein, ϕ and Φ represent arbitrary scalars and vectors respectively, d_m represents a number of dimensions of the problem under study: two dimensions or three dimensions, N is a number of particles in an affected area of the relevant particles; no is a particle density at initial time; gradient/divergence and the radius of Laplace operator r_e are 2.1r₀ and 4.0r₀, r₀ is an initial particle spacing; n₀ is the particle density at the initial moment.

A kernel function W (r_ij) and a parameter λ are defined as equation (1.1.2) and equation (1.1.3): an adjustment amount of the particle spacing as shown in equation (1.1.5):

$\begin{array}{cc}W\left(r_{\mathrm{ij}}\right)=\{\begin{array}{cc}\frac{r_e}{r_{\mathrm{ij}}}-1& 0\le r_{\mathrm{ij}}\le r_e\\ 0& r_{\mathrm{ij}}\ge r_e\end{array}& \left(1.1\mathrm{.2}\right)\\ \lambda =\frac{\sum _{j\ne i}^NW\left(r_{\mathrm{ij}}\right)r_{\mathrm{ij}}^2}{\sum _{j\ne i}^NW\left(r_{\mathrm{ij}}\right)}& \left(1.1\mathrm{.3}\right)\end{array}$

the particle density is defined in equation (1.1.4):

$\begin{array}{cc}n=\sum _{j\ne i}^NW\left(r_{\mathrm{ij}}\right)& \left(1.1\mathrm{.4}\right)\end{array}$

an adjustment amount of the particle spacing as shown in equation (1.1.5):

$\begin{array}{cc}\delta r_i=\sum _{j\ne i}\frac{r_{0-}r_{\mathrm{ij}}}{2}·\frac{r_{i-}r_j}{r_{\mathrm{ij}}},\mathrm{when}r_{\mathrm{ij}}\le r_0& \left(1.1\mathrm{.5}\right)\end{array}$

For the free surface particles distal from a main fluid domain, when a relative velocity between the particles is closer than a threshold before a prediction step, performing an operation as shown in equation (1.1.6) to set the relative velocity of the particle pair to be zero;

$\begin{array}{cc}\delta u_i=\sum _{j\ne i}-\square \left(r_{\mathrm{ij}}\right)u_{{\tau }_{\mathrm{ij}}},\mathrm{when}r_{\mathrm{ij}}-u_{{\tau }_{\mathrm{ij}}}\Delta t\le r_{\min },r_{\min }=0.3r_0& \left(1.1\mathrm{.6}\right)\end{array}$

wherein, u_rij is a relative tangential velocity of particle i and particle j, and the adjacent particles of solid and fluid are taken as 0.5 and 1.0 respectively.

Specific steps of the step S2 are as follows:

Introducing a variable A describing the motion and deformation of an object, which is defined as follows:

A=[X_Rb,θ,q_b]  (2.1.1)

The corresponding position, velocity and acceleration of the particles at an interface of the structure and the fluid are Γ_fsi,0^k, {dot over (Γ)}_fsi,0^k and {umlaut over (Γ)}_fsi,0^k, wherein k represents a time step, a position where a subscript ‘0’ locates represents a number of iterations.

Estimating the position, velocity and acceleration of the current time step from the position Γ_fsi,0^k−1, velocity {dot over (Γ)}_fsi,0^k−1 and acceleration {umlaut over (Γ)}_fsi,0^k−1 of the structure boundary particles of the previous time step according to an assumption of uniform acceleration motion, and calculating the estimated value of the structural boundary position of the current time step based on the calculated position, velocity and acceleration of the current time step.

The Lagrangian form non-viscous Navier-Stokes equation is shown in equation (3.1.1):

$\begin{array}{cc}\frac{\mathrm{Du}}{\mathrm{Dt}}=g-\frac{\nabla p}{\rho }& \left(3.1\mathrm{.1}\right)\end{array}$

wherein, u, p and ρ represent the velocity, pressure and density of the fluid respectively, which are all vectors pointing in a direction of gravity, and g is the acceleration of gravity.

Specific steps of the step S3 are as follows:

Firstly, using a classic projection method to decouple the speed and pressure as follows: without considering the pressure, updating the liquid to an intermediate state only by inertia, as shown in equation (3.1.2):

u _* =u _k +Δtg  (1)

l _* =l _k +Δtu _*  (2)(3.1.2)

wherein, l is a position vector, variables with subscripts * and k represent variables at an intermediate time step and a k^th time step respectively.

Specific steps of the step S4 are as follows:

Deriving the pressure Poisson equation by taking a divergence of the Navier-Stokes equation on both sides according to a continuity equation, as shown in (4.1.2), wherein, the continuity equation is shown in equation (4.1.1),

$\begin{array}{cc}\nabla ·u=0& \left(4.1\mathrm{.1}\right)\\ {\nabla }^2p_{k+1}=\frac{\rho \nabla ·u_{*}}{\Delta t}+\mathrm{\alpha \rho }\frac{n_0-n_k}{n_0\Delta t^2}& \left(4.1\mathrm{.2}\right)\end{array}$

wherein, n₀ and n_k are the particle densities at the initial time and the k^th time step respectively, and the particle density is defined as shown in equation (1.1.4).

The coefficient α is defined as shown in (4.1.3):

$\begin{array}{cc}\alpha =\{\begin{array}{cc}\frac{n_0-n_k}{n_0}+\Delta t\nabla ·u_k& \left(n_0-n_k\right)\nabla ·u_k\ge 0\\ \frac{n_0-n_k}{n_0}& \left(n_0-n_k\right)\nabla ·u_k\le 0\end{array}& \left(4.1\mathrm{.3}\right)\end{array}$

Boundary conditions for solving the pressure Poisson equation are as follows:

(1) Solid Boundary Conditions

Applying Neumann boundary conditions to solid particles, as shown in equation (4.1.4),

n·∇p _k+1=ρ(n·g−n·{dot over (u)}_bk+1)≈ρ(n·g−n·{dot over (u)}_bk)  (4.1.4)

wherein, {dot over (u)}_bk and {dot over (u)}_bk+1 the accelerations of the solid boundary particles at the k^th and (k+1)^th time steps; since {dot over (u)}_bk+1 is unknown when solving the fluid motion at the (k+1)^th time step, the value {dot over (u)}_bk of the previous time step is used as an approximation.

As shown in FIG. 8, fluid particles proximal to the solid boundary need to be corrected for the Laplace operator, wherein a pressure value is shown in equation (4.1.5), and a discrete term in equation (1.1.1) is corrected to make it consistent with equation (4.1.2); an intermediate velocity of the affected solid boundary particles is corrected as equation (4.1.6)-(1), which then are projected to tangential τ and normal n directions respectively to obtain equations (4.1.6)-(2) and (4.1.6)-(3),

$\begin{array}{cc}{p}_v=p_s+{\rho }_0\left(n·g-n·{\stackrel{.}{u}}_b\right)r_0& \left(4.1\mathrm{.5}\right)\\ u_{b^{*}}=u_{\mathrm{bk}+1}+\Delta t\frac{\nabla p_k}{{\rho }_0}& \left(1\right)\\ n·u_{b^{*}}=n·u_{\mathrm{bk}+1}+\Delta t\left(n·g-n·{\stackrel{.}{u}}_{b_k}\right)& \left(2\right)\left(4.1\mathrm{.6}\right)\\ \tau ·u_{b^{*}}=\tau ·u_{\mathrm{bk}+1}+\frac{\Delta t}{{\rho }_0}\frac{\partial p_k}{\partial \tau }& \left(3\right)\end{array}$

wherein, u_b* and u_bk+1 are the velocity of the solid boundary particles at the intermediate time and the (k+1)^th time step.

(2) Free Surface Boundary Condition

Since solving the pressure Poisson equation requires applying a boundary condition with zero pressure on the free surface boundary, it is necessary to identify a position of the free surface particles: for the case of two dimensions, each particle is assigned a circle centered on itself, with a radius of 1.05r₀, and discretized with 360 points evenly distributed on the circle; if all these points are covered by the circles of its neighboring particles, then the particle is an internal fluid particle (particle B in FIG. 1(a)), otherwise the particle is to be identified as a free surface particle (particle A in FIG. 1(a));

for the case of three dimensions, a two-step method is used: firstly, checking a number of adjacent particles using equation (4.1.7) to detect potential free surface particles,

N _p<δ_3dN_p0  (4.1.7)

in this embodiment, N_p0=32, β_3d=0.96.

Secondly, refining a search by establishing a spherical surface centered on the particle and having a radius of 1.05r₀. More specifically, finding the vector pointing to the most sparsely distributed direction of surrounding particles around the concerned particle by weighted average method, as shown in FIG. 1(b), and discretizing a circular area on the spherical surface by uniformly distributed points, if all these points are covered by the spherical surface of its neighboring particles, the particle is an internal fluid particle, otherwise the particle will be identified as a free surface particle.

Calculation equation of step S5 is as shown in equation (5.1.1),

$\begin{array}{cc}{u}_{k+1}=u_{*}-\Delta t·\frac{\nabla p_{k+1}}{{\rho }_0}r_{k+1}=r_k+\Delta t·u_{k+1}& \left(5.1\mathrm{.1}\right)\end{array}$

Specific steps of the step S6 are as follows:

S6.1. Establishment of Structural Motion Description Model

The structural motion description model uses a fixed global coordinate system XOY and an local coordinate system soη, wherein an origin of the local coordinate system is selected at a center of gravity of a beam of the structural motion description model, thus a position of any point on the beam of the structural motion description model can be described by equation (6.1.1), schematic diagram shown in FIG. 2.

X=X _Rb +Rξ ^T  (6.1.1)

wherein, X_Rb is a position of a center of mass, and equations (6.1.2) and (6.1.3) define local coordinates ξ of the beam of the structural motion description model and rotation matrix R of the beam of the structural motion description model,

$\begin{array}{cc}\xi =\left[s,{\eta }_l+{\eta }_0\right]& \left(6.1\mathrm{.2}\right)\\ R=\left[\begin{array}{cc}\cos \left(\theta \right)& -\sin \left(\theta \right)\\ \sin \left(\theta \right)& \cos \left(\theta \right)\end{array}\right]& \left(6.1\mathrm{.3}\right)\end{array}$

wherein, θ is an angle between O-X axis and o-s axis counterclockwise, η₁ is a coordinate of an undeformed structure, which is a constant for a specific point on the beam of the structural motion description model, and η₀ is a deflection of a mid-plane, the deflection of any plane parallel to the mid-plane taking the same value of the mid-plane; according to a modal superposition theory, the deflection η₀ of the mid-plane is expanded into equation (6.1.4).

Equations (6.1.5) and (6.1.6) define a column vector φ of a modal function and corresponding general coordinates q_b,

φ=[φ₁,φ₂,φ₃, . . . ]^T  (6.1.5)

q _b=[q_b1,q_b2,q_b3, . . . ]^T  (6.1.6)

The modal function satisfies the following orthogonal relationship:

$\begin{array}{cc}\underset{s_1}{\overset{s_2}{\int }}{\mathrm{\varphi \rho }}_l{\varphi }^T\mathrm{ds}=I_u,I_u\mathrm{is}\mathrm{an}\mathrm{identity}\mathrm{matrix}& \left(6.1\mathrm{.7}\right)\\ \underset{s_1}{\overset{s_2}{\int }}\frac{d^2\varphi }{{\mathrm{ds}}^2}\mathrm{EI}\frac{d^2{\varphi }^T}{{\mathrm{ds}}^2}\mathrm{ds}=\Lambda ,\Lambda =\mathrm{diag}\left({\omega }_i^2\right)i=1,2,3,\dots & \left(6.1\mathrm{.8}\right)\end{array}$

wherein s₁ and s₂ are upper and lower limits integrated along o-s axis, ρ₁ is an line density of the beam of the structural motion description model, E represents an elastic modulus of a structure, and I₁ is an inertia matrix of the structure;

S6.2. Establishment of the Governing Equation of Solid Response Based on Lagrange Equation

Establishing a model based on Lagrange equation:

$\begin{array}{cc}\frac{d}{\mathrm{dt}}\left(\frac{\partial T}{\partial {\stackrel{.}{q}}_j}\right)+\frac{\partial U}{\partial q_j}-\frac{\partial T}{\partial q_j}=Q_jj=1,2,3,\dots & \left(6.2\mathrm{.1}\right)\end{array}$

wherein T and U represent an kinetic energy and an potential energy of the structural motion description model, q_j is a general coordinate corresponding to any rigid-flexible mode, and Q_j is a non-conservative force corresponding to the j^th coordinate.

The beam of the structural motion description model is an elastic non-uniform beam, and specific forms of T, U, q_j, Q_j of the elastic non-uniform beam are as follows:

$\begin{array}{cc}\begin{array}{cc}T=& \frac{1}{2}\underset{s_1}{\overset{s_2}{\int }}\underset{{\eta }_1}{\overset{{\eta }_2}{\int }}{\stackrel{.}{X}}^T{\rho }_s\stackrel{.}{X}d\eta \mathrm{ds}=\\ & \frac{1}{2}\underset{s_1}{\overset{s_2}{\int }}\underset{{\eta }_1}{\overset{{\eta }_2}{\int }}{\left({\stackrel{.}{X}}_R+\mathrm{RU}\stackrel{.}{\theta }\xi +R\stackrel{.}{\xi }\right)}^T{\rho }_s\left({\stackrel{.}{X}}_R+\mathrm{RU}\stackrel{.}{\theta }\xi +R\stackrel{.}{\xi }\right)d\eta \mathrm{ds}\\ =& \frac{1}{2}[M\left({\stackrel{.}{X}}_R^2+{\stackrel{.}{Y}}_R^2\right)+{\stackrel{.}{\theta }}^2I_b+{\stackrel{.}{\theta }}^2q_b^Tq_b+{\stackrel{.}{q}}_b^T{\stackrel{.}{q}}_b-\\ & 2\stackrel{.}{\theta }\left({\stackrel{.}{X}}_R\cos \left(\theta \right)+{\stackrel{.}{Y}}_R\sin \left(\theta \right)\right)q_b^T{\psi }_0+\\ & 2\left(-{\stackrel{.}{X}}_R\sin \left(\theta \right)+{\stackrel{.}{Y}}_R\cos \left(\theta \right)\right){\stackrel{.}{q}}_b^T{\psi }_0+2\stackrel{.}{\theta }{\stackrel{.}{q}}_b^T{\psi }_1]\end{array}& \left(6.2\mathrm{.2}\right)\\ \begin{array}{cc}U=& \frac{1}{2}\underset{s_1}{\overset{s_2}{\int }}\frac{d^2{\eta }_0}{{\mathrm{ds}}^2}\mathrm{EI}\frac{d^2{\eta }_0}{{\mathrm{ds}}^2}\mathrm{ds}+{\mathrm{MgY}}_R\\ =& \frac{1}{2}q_b^T\left(\underset{s_1}{\overset{s_2}{\int }}\frac{d^2\varphi }{{\mathrm{ds}}^2}\mathrm{EI}\frac{d^2{\varphi }^T}{{\mathrm{ds}}^2}\mathrm{ds}\right)q_b+\\ & {\mathrm{MgY}}_R=\frac{1}{2}q_b^T\Lambda q_b+{\mathrm{MgY}}_R\end{array}& \left(6.2\mathrm{.3}\right)\end{array}$

wherein η₁ and η₂ are the upper and lower limits of integral along axis η, matrix integral U, column vector integrals Ψ₀ and Ψ₁, the matrix integral U and column vector integrals Ψ₀, Ψ₁ defined in equations (6.2.4) and (6.2.6):

$\begin{array}{cc}U=\left[\begin{array}{cc}0& -1\\ 1& 0\end{array}\right]& \left(6.2\mathrm{.4}\right)\\ {\psi }_0={[{\psi }_{01},{\psi }_{02},{\psi }_{03},\dots ]}^T=\underset{s_1}{\overset{s_2}{\int }}{\mathrm{\varphi \rho }}_l\mathrm{ds}& \left(6.2\mathrm{.5}\right)\\ {\psi }_1={[{\psi }_{11},{\psi }_{12},{\psi }_{13},\dots ]}^T=\underset{s_1}{\overset{s_2}{\int }}s{\mathrm{\varphi \rho }}_l\mathrm{ds}& \left(6.2\mathrm{.6}\right)\end{array}$

The governing equations of the dynamic response for elastic non-uniform beam under a wave impact can be obtained by substituting equations (6.2.2) and (6.2.3) into equation (6.2.1), which is shown in equation (6.2.7):

M{umlaut over (X)} _Rb+{dot over (θ)}² sin(θ)ψ₀^Tq_b−2{dot over (θ)} cos(θ)Ψ₀^T{dot over (q)}_b−{umlaut over (θ)} cos(θ)Ψ₀^Tq_b−sin(θ)Ψ₀^T{umlaut over (q)}_b=Q_XRb

MŸ _Rb−{dot over (θ)}² cos(θ)ψ₀^Tq_b−2{dot over (θ)} sin(θ)ψ₀^T{dot over (q)}_b−{umlaut over (θ)} sin(θ)ψ₀^Tq_b+cos(θ)ψ₀^T{umlaut over (q)}_b+Mg=Q_YRb

−({umlaut over (X)}_Rb cos(θ)+Ÿ_Rb sin(θ))ψ₀^Tq_b+I_b{umlaut over (θ)}+{umlaut over (θ)}q_b^Tq_b+2{dot over (θ)}{dot over (q)}_b^Tq_b+ψ₁^T{umlaut over (q)}_b=Q_θ

(−{umlaut over (X)}_Rb sin(θ)+Ÿ_Rb cos(θ))ψ₀+{umlaut over (θ)}ψ₁+{umlaut over (q)}_b+Λq_b−{dot over (θ)}²q_b=Q_qb  (6.2.7)

Equations (6.2.8)-(6.2.11) provide the non-conservative forces corresponding to the general coordinates of the rigid body and the elastic modal:

Q _{X Rb} = pn _s dl  (6.2.8)

Q _{Y Rb} = pn _η dl  (6.2.9)

Q _θ = p(X_pn_p−Y_pn_s)dl  (6.2.10)

Q _{q b} = p(n·e_η)φdl  (6.2.11)

wherein, n_s and n_η are local normal components of s and η respectively, expressed by a fixed global coordinate system, X_p and Y_p are overall X and Y coordinates of a point of action of an intensity of pressure p, e_η is an unit vector of o-η axis, and the integration is carried out around a perimeter of the beam of the structural motion description mode;

S6.3. Solving the Governing Equation of Solid Response which Couples the Rigid Body and the Elastic Mode.

Specific steps of the step S8 are as follows:

Assuming that all fluid and structural variables are known in time step t=t_k−1, and then the detailed interaction process at the next time t=t_k is as follows:

S8.1. assuming that an acceleration of a structure at the time step t=t_k is the same as the acceleration at the previous time step t=t_k−1, that is, {umlaut over (D)}^k={umlaut over (D)}^k−1, a position D^k and a speed {dot over (D)}^k of a position at the time step t=t_k can also be calculated by a finite difference method.

Calculating the corresponding position, velocity and acceleration of each point on the interface between fluid and structure through equations (8.1.1) (8.1.2) and (8.1.3)˜(8.1.6), that is, Γ_fsi,0^(k), {dot over (Γ)}_fsi,0^(k) and {umlaut over (Γ)}_fsi,0^(k) can be calculated from the newly calculated structural motion and deformation parameters D, that is equation (2.1.1),

$\begin{array}{cc}X_R=X_{\mathrm{CR}}+R_{{\mathrm{RX}}_R}& \left(8.1\mathrm{.1}\right)\\ X_{\mathrm{fi}}=X_{\mathrm{cfi}}+R_{\mathrm{fi}}{\xi }_i^T& \left(8.1\mathrm{.2}\right)\\ {\stackrel{.}{X}}_R={\stackrel{.}{X}}_{\mathrm{CR}}+{\stackrel{.}{R}}_{{\mathrm{RX}}_R}={\stackrel{.}{X}}_CR+{\stackrel{.}{R}}_RU{\stackrel{.}{\theta }}_R{\chi }_R& \left(8.1\mathrm{.3}\right)\\ {\stackrel{.}{X}}_{\mathrm{fi}}={\stackrel{.}{X}}_{\mathrm{cR}}+{\stackrel{.}{R}}_R{\chi }_{\mathrm{ofi}}+{\stackrel{.}{R}}_{\mathrm{fi}}{\xi }_i^T+{\stackrel{.}{R}}_{\mathrm{fi}}{\stackrel{.}{\xi }}_i^T={\stackrel{.}{X}}_{\mathrm{CR}}+R_RU{\stackrel{.}{\theta }}_R{\chi }_{\mathrm{ofi}}+R_{\mathrm{fi}}U{\stackrel{.}{\theta }}_{\mathrm{fi}}{\chi }_{\mathrm{ofi}}{\xi }_i+R_{\mathrm{fi}}{\stackrel{.}{\xi }}_i& \left(8.1\mathrm{.4}\right)\\ {\ddot{X}}_R={\ddot{X}}_{\mathrm{CR}}+{\ddot{R}}_R{\chi }_R={\ddot{X}}_R+\left(R_RU{\ddot{\theta }}_R-R_R{\ddot{\theta }}_R^2\right){\chi }_R& \left(8.1\mathrm{.5}\right)\\ \begin{array}{cc}{\ddot{X}}_{\mathrm{fi}}=& {\ddot{X}}_{\mathrm{CR}}+{\ddot{R}}_R{\chi }_{\mathrm{ofi}}+{\ddot{R}}_{\mathrm{fi}}{\xi }_i+2{\stackrel{.}{R}}_{\mathrm{fi}}{\stackrel{.}{\xi }}_i+R_{\mathrm{fi}}{\ddot{\xi }}_i\\ =& {\ddot{X}}_{\mathrm{CR}}+\left(R_RU{\ddot{\theta }}_R-R_R{\ddot{\theta }}_R^2\right){\chi }_{\mathrm{ofi}}+\left(R_{\mathrm{fi}}U{\ddot{\theta }}_{\mathrm{fi}}-R_{\mathrm{fi}}{\ddot{\theta }}_{\mathrm{fi}}^2\right){\xi }_i+\\ & 2R_{\mathrm{fi}}U{\ddot{\theta }}_{\mathrm{fi}}{\stackrel{.}{\xi }}_i+R_{\mathrm{fi}}{\ddot{\xi }}_i\end{array}& \left(8.1\mathrm{.6}\right)\end{array}$

wherein, X_R is a global coordinates of rigid body, X_CR is a center-of-mass coordinates, R_R is a rotation matrix that associates the local coordinate system soη with the fixed global coordinate system XOY; X_fi represents the coordinates on a center line; and represents the coordinates of the beam center line in the local coordinate system soη, as represented in equation (8.1.7), wherein 77; is deflection function of the beam of the structural motion description model;

ξ_i=[s_i,η_i]  (8.1.7)

S8.2. calculating the fluid motion at time t=t_k according to the modified MPS method by using the updated information of the interface as a new boundary condition, by using a newly updated pressure p_fsi,i+1^k; the i^th iteration of the interface at time t=t_k can be conducted;

S8.3. comparing structural position differences Γ_fsi,i^k and {tilde over (Γ)}_fsi,i+1^k in step S2 and step S7, if a convergence condition is satisfied, as shown in equation (8.3.1), performing step S8.1 for the next time step t=t_k+1 operation,

|{tilde over (Γ)}_fsi,i+1^k−Γ_fsi,i^k|≤ϵ  (8.3.1)

otherwise, using equation (8.3.2) to correct the structural position D_i+1^k of the (i+1)^th iteration, and using Newmark method to update velocity {dot over (D)}_i+1^k and acceleration {umlaut over (D)}_i+1^k, calculating interface variables Γ_fsi,i+1^k, {dot over (Γ)}_fsi,i+1^k and {umlaut over (Γ)}_fsi,i+1^k based on D_i+1^k, {dot over (D)}_i+1^k, and {umlaut over (D)}_i+1^k; using these modified interface variables, to perform the (i+1)^th iteration by returning to step S2,

D _i+1 ^k=χ_i{tilde over (D)}_i+1^k+(1−χ_i)D_i^k  (8.3.2)

wherein, χ_i is Aitken relaxation factor, a core idea of which is to improve the accuracy of prediction by using the first two iterations. χ_i is calculated by equation (8.3.3):

$\begin{array}{cc}{\chi }_i={\chi }_{i-1}\frac{{\left({\mathrm{\Delta \Gamma }}_{\mathrm{fsi},i}^k\right)}^T\left({\mathrm{\Delta \Gamma }}_{\mathrm{fsi},i+1}^k-{\mathrm{\Delta \Gamma }}_{\mathrm{fsi},i}^k\right)}{{\left({\mathrm{\Delta \Gamma }}_{\mathrm{fsi},i+1}^k-{\mathrm{\Delta \Gamma }}_{\mathrm{fsi},i}^k\right)}^T\left({\mathrm{\Delta \Gamma }}_{\mathrm{fsi},i+1}^k-{\mathrm{\Delta \Gamma }}_{\mathrm{fsi},i}^k\right)}& \left(8.3\mathrm{.3}\right)\end{array}$

wherein, Γ_fsi,i^k={tilde over (Γ)}_fsi,i^k−Γ_fsi,i−1^k. The initial value of each first iteration is 0.2.

The entire process is shown in FIG. 9.

Hereafter, the method based on the moving particle semi-implicit method and modal superposition method to solve the strong nonlinear time-domain hydroelastic problem is used to calculate the three-dimensional hull slamming problem, and the applicability of the method is explained.

Specific Calculation Process:

1. Model Selection

A 1:100 model of a 4,600-ton tanker was used to simulate an interaction between a dam-break wave and a ship. Assuming that the ship's dynamic response is two-dimensional (that is, only considering pitch, heave and vertical bending), this makes the elastic non-uniform beam suitable for this situation. The main parameters of the ship model are shown in Table 1. The ship model is placed in the pool facing the impact of the dam failure. The initial configuration is shown in FIG. 3. The width and length of the water tank are about three times the width and length of the ship.

TABLE 1
Main parameters of ship model
Length of ship      1.8 m
Width of ship       0.322 m
Ship body's draught 0.12 m
Mass of ship body   54.99 kg

2. Particle Discretization

The surface of the ship model is discretized using particles and the detailed procedure is as follows:

(1) the ship model is first discretized along the longitude direction with several cross-sections, i.e. stations (the norms of these cross-sections are parallel to the longitude direction) as shown in FIG. 4(a), the distance between each section is set as the value for initial fluid particles.

(2) for each section, the outline of the cross section is then discretized by evenly distributed particles.

The cross section and surface discretization of the ship model is shown in FIG. 4. The particle spacing on the surface of the ship model is approximately the same as the resolution of the flow field, which is 0.02 m (about 1/90 of the length of the ship model). The total number of particles is 289417, including 193284 fluid particles.

3. Parameter Selection

Assuming that the mass distribution of the ship model is the same as its water consumption, the thickness of the ship model is 0.05 m for the calculation of the second moment of section. FIG. 5 shows the distribution of the mass and second moment of the ship model.

Myklestad's method is used to calculate the mode function of the elastic non-uniform beam. The results of the first three order mode functions are shown in FIG. 6. Young's modulus E=2×10⁴ N/m², and the first three circular natural frequency modes are: ω₁=3.7857, ω₂=9.7777, ω₃=18.2063

4. Calculation Result

Some typical time-step free surface shapes and ship deformation diagrams are shown in FIG. 7. The free surface particles, the particles of the ship model part above the free surface and the particles of the ship model part below the free surface are marked as red, blue and black, respectively. It can be seen from FIG. 7 that no free surface particles are erroneously identified below the free surface. In this simulation, the deformation amplitude is significantly larger than the normal ship model structure. The lower stiffness is deliberately selected to test the coupling ability between the modified modal superposition and the MPS model. FIGS. 7 (a), (b), (c) show the shock wave motion processes at different times of dam failure. The fact that the draft of the ship model remains unchanged proves that the buoyancy calculated by MPS is accurate and stable for the support of the weight of the ship model. The impact of the waves on the ship's bow occurs at t=0.418 s in FIG. 7(c). When the wave propagates along the ship model, the maximum bending moment is generated at the peak of the wave (FIGS. 7 (d)-(f)), and the ship model accordingly generates the largest deformation. Qualitative analysis shows that the coupling between the modified MPS and the mentioned structural model can provide reasonable results for the simulation of wave-structure interaction.

At last, it should be stated that the above various embodiments are only used to illustrate the technical solutions of the present invention without limitation; and despite reference to the aforementioned embodiments to make a detailed description of the present invention, those of ordinary skilled in the art should understand: the described technical solutions in above various embodiments may be modified or the part of or all technical features may be equivalently substituted; while these modifications or substitutions do not make the essence of their corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for analyzing hydroelastic effect of a marine structure based on moving particle semi-implicit method and modal superposition method, comprising the following steps: S1. spatially discretizing a boundary of a flow field with mesh-less discrete points, wherein the boundary of the flow field is a surface of the marine structure in contact with a fluid, wherein if a previous time step is initial time, velocity and intensity of pressure of fluid particles need to be initialized;S2. determining an estimated value of a structural boundary position of a current time step according to motion information of a solid boundary of the previous time step, wherein if the previous time step is the initial time, an initial value of a structural motion is adopted;S3. updating a position and a velocity of fluid particles by Navier-Stokes equation according to flow field information of the previous time step, excluding fluid pressure;S4. deriving a pressure Poisson equation and solving the pressure Poisson equation;S5. updating the position and the velocity of the fluid particles using a pressure value in S4;S6. solving the governing equation of solid response, which couples the rigid body and elastic mode, by using the pressure value obtained in S4;S7. updating the structure motion and deformation by using dynamic information obtained in step S6, which in turn provides a new fluid boundary position;S8. comparing structural position differences in step S2 and step S7, wherein if a convergence condition is satisfied, a next time step calculation is performed, otherwise, repeating steps S2 to S8 at a current time step until the next time step calculation is performed after the convergence condition is satisfied.

2. The method according to claim 1, wherein, specific steps of the step S1 are as follows: discretizing based on uniformly distributed particles, as shown in equation (1.1.1)

3. The method according to claim 2, wherein, specific steps of the step S2 are as follows: introducing a variable A describing the motion and deformation of an object, which is defined as follows: A=[XRb,θ,qb]  (2.1.1)wherein, XRb represents a position of a center of mass of a rigid body motion, B represents a deadrise angle of the structure, and qb represents general coordinates corresponding to a modal function;the corresponding position, velocity and acceleration of the particles at an interface of the structure and the fluid are Γfsi,0k, {dot over (Γ)}fsi,0k, and {umlaut over (Γ)}fsi,0k, wherein k represents a time step, a position where a subscript ‘0’ locates represents a number of iterations;estimating a position, velocity and acceleration of the current time step from the position Γfsi,0k−1, velocity {dot over (Γ)}fsi,ik−1 and acceleration {umlaut over (Γ)}fsi,0k−1 of the structure boundary particles of the previous time step according to an assumption of uniform acceleration motion, and calculating an estimated value of the structural boundary position of the current time step based on the calculated position, velocity and acceleration of the current time step.

4. The method according to claim 3, wherein, in the step S3, the Lagrangian form non-viscous Navier-Stokes equation is shown in equation (3.1.1):

5. The method according to claim 4, wherein, specific steps of the step S4 are as follows: deriving the pressure Poisson equation by taking a divergence of the Navier-Stokes equation on both sides according to a continuity equation, as shown in equation (4.1.2), wherein, the continuity equation is shown in equation (4.1.1),

6. The method according to claim 5, wherein, calculation equation of step S5 is as shown in equation (5.1.1),

7. The method according to claim 6, wherein, specific steps of the step S6 are as follows: S6.1. establishment of structural motion description modelthe structural motion description model uses a fixed global coordinate system XOY and an attachment local coordinate system soη, wherein an origin of the attachment local coordinate system is selected at a center of gravity of a beam of the structural motion description model, a position X of any point on the beam of the structural motion description model can be described by equation (6.1.1), X=XRb+RξT  (6.1.1)wherein, XRb is a position of a center of mass of the beam of the structural motion description model, and equations (6.1.2) and (6.1.3) define local coordinates ξ of the beam of the structural motion description model and rotation matrix R of the beam of the structural motion description model,

8. The method according to claim 6, wherein, specific steps of the step S8 are as follows: assuming that all fluid and structural variables are known in a time step t=tk−1, and then the detailed interaction process at the next time step t=tk is as follows:S8.1. assuming that an acceleration of a structure at the time step t=tk is the same as the acceleration at the previous time step t=tk−1, that is, {umlaut over (D)}k={umlaut over (D)}k−1, a position Dk and a speed {dot over (D)}k at the time step t=tk can also be calculated by a finite difference method;calculating the corresponding position, velocity and acceleration of each point on the interface between fluid and structure through equations (8.1.1)˜(8.1.6), that is Γfsi,0(k), Γfsi,0(k) and Γfsi,0(k) can be calculated from the newly calculated structural motion and deformation parameters A, that is, equation (2.1.1),