WO2020192126A1 - 一种基于改进的移动粒子半隐式法和模态叠加方法求解强非线性时域水弹性问题设计方法 - Google Patents

一种基于改进的移动粒子半隐式法和模态叠加方法求解强非线性时域水弹性问题设计方法 Download PDF

Info

Publication number
WO2020192126A1
WO2020192126A1 PCT/CN2019/115597 CN2019115597W WO2020192126A1 WO 2020192126 A1 WO2020192126 A1 WO 2020192126A1 CN 2019115597 W CN2019115597 W CN 2019115597W WO 2020192126 A1 WO2020192126 A1 WO 2020192126A1
Authority
WO
WIPO (PCT)
Prior art keywords
equation
particle
particles
time step
formula
Prior art date
Application number
PCT/CN2019/115597
Other languages
English (en)
French (fr)
Inventor
孙哲
邓彦增
张桂勇
宗智
王浩
Original Assignee
大连理工大学
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 大连理工大学 filed Critical 大连理工大学
Publication of WO2020192126A1 publication Critical patent/WO2020192126A1/zh
Priority to US16/950,886 priority Critical patent/US20210086877A1/en

Links

Images

Classifications

    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B63—SHIPS OR OTHER WATERBORNE VESSELS; RELATED EQUIPMENT
    • B63B—SHIPS OR OTHER WATERBORNE VESSELS; EQUIPMENT FOR SHIPPING 
    • B63B71/00—Designing vessels; Predicting their performance
    • B63B71/10—Designing vessels; Predicting their performance using computer simulation, e.g. finite element method [FEM] or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/25—Design optimisation, verification or simulation using particle-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Definitions

  • the invention relates to a design method based on an improved moving particle semi-implicit method (MPS) and a modal superposition method for solving strong nonlinear time-domain hydroelastic problems, and belongs to the field of ship engineering.
  • MPS moving particle semi-implicit method
  • Structural integrity analysis in this case requires time domain simulation.
  • the traditional potential flow-based solver cannot directly simulate the severe deformation of the free surface during the slamming process. This requires the CFD-based solver to perform flow calculations. Partial combination or full application of the CFD model based on viscous flow is hydroelastic calculation The ultimate goal.
  • particle methods such as SPH and MPS have no grid distortion problems and the flow field is updated by 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 violent fluid and structure interaction problems.
  • the typical motion of a ship is a large rigid body motion plus relatively small elastic deformation. This means that the modal superposition theory is sufficient to accurately calculate the elastic deformation.
  • the interaction between the rigid body and the elastic mode 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.
  • step S7 Use the dynamic information obtained in step S6 and the various elastic modes of the solid structure to calculate the overall elastic deformation of the structure through a linear superposition method to update the structure motion and deformation, and determine the boundary position of the new fluid;
  • step S8 Compare the difference between the structure position in step S2 and step S7, if the convergence condition is met, proceed to the next time step calculation, otherwise repeat step S2 to step S8 at the current time step, until the convergence condition is met, proceed to the next time step Calculation.
  • step S1 The specific steps of step S1 are as follows:
  • the gradient, divergence, and Laplacian of physical quantities use the weighted average method, which is based on uniformly distributed particles, as shown in equation (1.1.1).
  • ⁇ and ⁇ represent arbitrary scalar and vector, respectively
  • d m represents the dimension of the research problem: two-dimensional or three-dimensional
  • N is the number of particles in the affected area of the relevant particle
  • n 0 is the particle density at the initial moment
  • r is the particle
  • Free surface particle technology.
  • the displacement of the particle is directly adjusted, not by applying a force.
  • the particle position is slightly shifted to make its distribution regular. This technique can also be seen as a process of mesh reconstruction.
  • the pressure Poisson equation is derived based on the incompressibility condition, the resulting pressure will keep the distance between adjacent particles approximately the same value (ie, the initial particle distance r 0 ).
  • step S2 introduce the variable A describing the movement and deformation of the object, which is defined as follows:
  • X Rb represents the position of the center of mass of the rigid body motion
  • ⁇ represents the bottom lift angle of the structure
  • q b represents the general coordinates corresponding to the modal function.
  • the position, velocity and acceleration of the boundary particles are with Estimate the position, speed and acceleration of the current time step according to the assumption of uniformly accelerating motion, and calculate the estimated value of the structure boundary position at the current time step with the calculated position, speed and acceleration at this moment.
  • step S3 the Navier-Stokes equation is an incompressible and inviscid Navier-Stokes equation based on the Lagrangian form, which is shown in equation (3.1.1):
  • u, p and ⁇ represent the velocity, pressure and density of the fluid respectively, which are all vectors pointing in the direction of gravity, and g represents the acceleration of gravity;
  • step S3 The specific steps of step S3 are as follows:
  • the classical projection method is used to decouple the velocity and pressure as follows:
  • l is the position vector
  • the variables with the subscript * and k respectively represent the variables of the intermediate time step and the kth time step
  • u represents the velocity of the fluid
  • step S4 The specific steps of step S4 are as follows:
  • n 0 and n k are the particle density at the initial time and k-th time step, respectively, and the particle density is defined as shown in formula (1.1.4);
  • the ⁇ model does not require special calibration according to different problems, and has stronger stability compared with other similar models.
  • Fluid particles close to the solid boundary need to modify the Laplacian to make it consistent with the Neumann boundary condition on the solid boundary to make up for the lack of adjacent particles.
  • the pressure value is shown in equation (4.1.5) , Considering the pressure value of the virtual particle outside the solid boundary, modify the discrete term in equation (1.1.1) to make it consistent with equation (4.1.2).
  • the source term also includes the effect of solid boundary particles in the influence domain. Therefore, due to the same consistency reason, the intermediate velocity of the affected solid boundary particles is corrected to formula (4.1.6)-(1), which respectively project To the tangential ⁇ and the normal direction n, the equations (4.1.6)-(2) and (4.1.6)-(3) are obtained,
  • n the normal direction
  • ⁇ 3d is a parameter less than 1, which is selected according to the situation.
  • N p0 represents the initial particle density within the circled range;
  • the search is refined by establishing a spherical surface with a particle as the center and a radius of 1.05r 0 , specifically: a weighted average method is used to determine a vector pointing to the sparsest direction of particle distribution around the particle, and the spherical surface is discretized by uniformly distributed points In the upper circular area, if all these points are covered by the spherical surfaces of their neighboring particles, the particle is an internal fluid particle, otherwise the particle will be recognized as a free surface particle.
  • step S5 The calculation formula of step S5 is shown in formula (5.1.1),
  • step S6 The dynamic response of a typical ocean structure to wave impact can be expressed as rigid body motion with relatively 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 coupling rigid body and elastic mode based on Lagrangian equation, and the specific steps of step S6 are as follows:
  • the structure motion description model adopts the fixed global coordinate system XOY and the attachment local coordinate system so ⁇ .
  • the origin of the attachment local coordinate system is selected at the center of gravity of the beam of the structure motion description model. Therefore, any point on the beam of the structure motion description model
  • Equations (6.1.2) and (6.1.3) define the local coordinate ⁇ of the beam of the structure motion description model and the rotation of the beam of the structure motion description model.
  • ⁇ is the counterclockwise angle between the OX axis and the os axis
  • ⁇ l is the coordinate of the undeformed structure
  • a specific point on the beam of the structural motion description model is a constant
  • ⁇ 0 is the deflection of the midplane, arbitrarily parallel to The deflection of the mid-plane is the same value as the mid-plane.
  • the deflection ⁇ 0 of the mid-plane is expanded into equation (6.1.4):
  • Equations (6.1.5) and (6.1.6) define the column vector of the modal function And the corresponding general coordinates q b ,
  • the modal function satisfies the following orthogonal relationship:
  • I u is the identity matrix (6.1.7)
  • ⁇ l is the linear density of the beam of the structure motion description model
  • E is the elastic modulus of the structure
  • I I is the inertia matrix of the structure
  • the main feature of the improved modal superposition method model is that the rigid body mode and the elastic mode are coupled in the time domain.
  • the influence of interactions such as the inertial effect of rigid body motion on elastic deformation is considered.
  • This is an extension of the traditional hydroelastic theory.
  • the rigid body and the elastic part are only connected by the hydrodynamic coefficient of the corresponding cross term.
  • T and U represent the kinetic energy and potential energy of the structural motion description model
  • q j is the general coordinate corresponding to any rigid-flexible mode
  • Q j is the non-conservative force corresponding to the j-th coordinate
  • the beam of the structural motion description model is an elastic non-uniform beam.
  • the specific forms of T, U, q j and Q j of the elastic non-uniform beam are as follows:
  • ⁇ s represents the linear density
  • M represents the mass of the structure
  • X R and Y R respectively correspond to the X and Y coordinates of the center of mass of the structure
  • g represents the acceleration of gravity
  • ⁇ 1 , ⁇ 2 are the integral ⁇ along the axis
  • the matrix integral U The upper and lower limits of the column vector integrals ⁇ 0 and ⁇ 1.
  • the matrix integral U and column vector integrals ⁇ 0 and ⁇ 1 are defined in equations (6.2.4) and (6.2.6):
  • Equations (6.2.8)-(6.2.11) give the non-conservative force corresponding to the general coordinates of rigid body and elastic mode:
  • n s and n ⁇ are the local normal components of s and ⁇ respectively, expressed in a fixed global coordinate system
  • X p and Y p are the overall X and Y coordinates of the point of pressure p
  • e ⁇ is the o- ⁇ axis
  • the unit vector of, the integral is carried out around the circumference of the beam of the structural motion description model. It can be seen from equation (6.2.7) that the interaction between the rigid body and the elastic mode, that is, the cross term is caused by angular motion.
  • step S8 The specific steps of step S8 are as follows:
  • X R is the global coordinate of the rigid body
  • X CR is the center of mass coordinate
  • R R is the rotation matrix linking the attachment local coordinate system so ⁇ with the fixed global coordinate system XOY
  • X fi represents the coordinates on the center line
  • ⁇ i represents The coordinates of the beam of the structural motion description model on the centerline of the attachment local coordinate system so ⁇ are expressed in formula (8.1.7), where ⁇ i is the deflection function of the beam of the structural motion description model;
  • ⁇ i is the Aitken relaxation factor, and its core idea is to improve the accuracy of prediction by using the first two iterations. Its value is calculated by formula (8.3.3):
  • the invention discloses a design method based on an improved moving particle semi-implicit method and a modal superposition method for solving strong nonlinear time-domain hydroelastic problems, which is based on an improved moving particle semi-implicit method (MPS) ) Hydrodynamic solver, structural response solver based on the modal superposition method of coupled rigid body and elastic mode, the two adopt an iterative method to achieve strong coupling in the time domain.
  • MPS moving particle semi-implicit method
  • Hydrodynamic solver structural response solver based on the modal superposition method of coupled rigid body and elastic mode, the two adopt an iterative method to achieve strong coupling in the time domain.
  • the invention improves the traditional hydroelasticity calculation based on potential flow that cannot directly simulate the severe deformation of the free surface and the slamming phenomenon of the hull, and considers the mutual influence between the rigidity of the structure and the flexible motion mode.
  • Figure 1 is a two-dimensional and three-dimensional free surface particle identification scheme diagram.
  • Figure 2 is a schematic diagram of an elastic non-uniform beam.
  • Figure 3 is a schematic diagram of the calculated initial configuration (not to scale).
  • Figure 4 shows the cross-section and surface discretization of the ship model.
  • Figure 5 shows the mass (a) and second moment (b) distribution diagram of the ship model.
  • Figure 6 shows the first three modes of the ship model.
  • Figure 8 is a demonstration of the Laplacian used to compensate for virtual particles near the solid boundary.
  • Figure 9 is a flow chart of the iterative process of fluid-structure interaction.
  • step S7 Use the dynamic information obtained in step S6 and the various elastic modes of the solid structure to calculate the overall elastic deformation of the structure through a linear superposition method to update the structure motion and deformation, and determine the boundary position of the new fluid;
  • step S8 Compare the difference between the structure position in step S2 and step S7, if the convergence condition is met, proceed to the next time step calculation, otherwise repeat step S2 to step S8 at the current time step, until the convergence condition is met, proceed to the next time step Calculation.
  • Steps S1-S4 are the hydrodynamic solver part (modified MPS method), and step S6 is the structural response solver.
  • step S1 Discrete based on uniformly distributed particles, as shown in formula (1.1.1).
  • ⁇ and ⁇ represent arbitrary scalar and vector, respectively
  • d m represents the dimension of the research problem: two-dimensional or three-dimensional
  • N is the number of particles in the affected area of the relevant particle
  • the domain radius r e takes the values 2.1r 0 and 4.0r 0
  • r 0 is the initial particle distance
  • n 0 is the particle density at the initial moment
  • the particle density is defined in formula (1.1.4):
  • Is the tangential relative velocity of particle i and particle j, and the adjacent particles of solid and fluid are respectively 0.5 and 1.0.
  • step S2 The specific steps of step S2 are as follows:
  • the position, velocity and acceleration of the boundary particles are with Estimate the position, speed and acceleration of the current time step according to the assumption of uniformly accelerating motion, and calculate the estimated value of the structure boundary position at the current time step with the calculated position, speed and acceleration at this moment.
  • Navier-Stokes equation is an incompressible and inviscid Navier-Stokes equation based on the Lagrangian form, which is shown in equation (3.1.1):
  • u, p and ⁇ represent fluid velocity, pressure and density respectively, which are all vectors pointing in the direction of gravity, and g represents the acceleration of gravity;
  • step S3 The specific steps of step S3 are as follows:
  • l is the position vector
  • the variables with subscripts * and k respectively represent the variables of the intermediate time step and the kth time step.
  • step S4 The specific steps of step S4 are as follows:
  • n 0 and n k are the particle density at the initial time and the k-th time step, respectively, and the particle density is defined as shown in formula (1.1.4);
  • the particle is an internal fluid particle (particle B in Figure 1(a)), otherwise it will be recognized as a free surface particle (particle A in Figure 1(a));
  • the search is refined by establishing a spherical surface with a particle as the center and a radius of 1.05r 0 , specifically: a weighted average method is used to determine a vector pointing to the most sparse direction of particle distribution around the particle, as shown in Figure 1(b) , Discrete circular areas on the sphere by uniformly distributed points. If all these points are covered by the spheres of their neighboring particles, the particle is an internal fluid particle, otherwise the particle will be recognized as a free surface particle.
  • step S5 The calculation formula of step S5 is shown in formula (5.1.1),
  • step S6 The specific steps of step S6 are as follows:
  • the structure motion description model adopts the fixed global coordinate system XOY and the attachment local coordinate system so ⁇ .
  • the origin of the attachment local coordinate system is selected at the center of gravity of the beam of the structure motion description model. Therefore, any point on the beam of the structure motion description model
  • the position of can be described by equation (6.1.1), the schematic diagram is shown in Figure 2.
  • X X Rb + R ⁇ T (6.1.1)
  • Equations (6.1.2) and (6.1.3) define the local coordinate ⁇ of the beam of the structure motion description model and the rotation matrix R of the beam of the structure motion description model.
  • ⁇ is the counterclockwise angle between the OX axis and the os axis
  • ⁇ l is the coordinate of the undeformed structure
  • a specific point on the beam of the structural motion description model is a constant
  • ⁇ 0 is the deflection of the midplane, arbitrarily parallel to The deflection of the mid-plane is the same value as the mid-plane.
  • the deflection ⁇ 0 of the mid-plane is expanded into equation (6.1.4),
  • Equations (6.1.5) and (6.1.6) define the column vector of the modal function And the corresponding general coordinates q b ,
  • the modal function satisfies the following orthogonal relationship:
  • I u is the identity matrix (6.1.7)
  • ⁇ l is the linear density of the beam of the structure motion description model
  • E is the elastic modulus of the structure
  • I I is the moment of inertia of the structure
  • T and U represent the kinetic energy and potential energy of the structural motion description model
  • q j is the general coordinate corresponding to any rigid-flexible mode
  • Q j is the non-conservative force corresponding to the j-th coordinate
  • the beam of the structural motion description model is an elastic non-uniform beam.
  • the specific forms of T, U, q j and Q j of the elastic non-uniform beam are as follows:
  • ⁇ 1 and ⁇ 2 are the upper and lower limits of the axis integral ⁇ , the matrix integral U, and the column vector integral ⁇ 0 and ⁇ 1.
  • the matrix integral U and the column vector integral ⁇ 0 and ⁇ 1 are defined in equation (6.2.4) and (6.2.6):
  • Equations (6.2.8)-(6.2.11) give the non-conservative force corresponding to the general coordinates of rigid body and elastic mode:
  • n s and n ⁇ are the local normal components of s and ⁇ respectively, expressed in a fixed global coordinate system
  • X p and Y p are the overall X and Y coordinates of the point of pressure p
  • e ⁇ is the o- ⁇ axis
  • the unit vector of, the integral is performed around the circumference of the beam of the structural motion description model.
  • step S8 The specific steps of step S8 are as follows:
  • X R is the global coordinate of the rigid body
  • X CR is the center of mass coordinate
  • R R is the rotation matrix linking the attachment local coordinate system so ⁇ with the fixed global coordinate system XOY
  • X fi represents the coordinates on the center line
  • ⁇ i represents The coordinates of the beam of the structural motion description model on the centerline of the attachment local coordinate system so ⁇ are expressed in formula (8.1.7), where ⁇ i is the deflection function of the beam of the structural motion description model;
  • ⁇ i is the Aitken relaxation factor, and its core idea is to improve the accuracy of prediction by using the first two iterations. Its value is calculated by formula (8.3.3):
  • the initial value for each first iteration is 0.2.
  • the three-dimensional hull slamming problem is calculated by a design method based on the improved moving particle semi-implicit method and modal superposition method to solve the strong nonlinear time-domain hydroelasticity problem to illustrate the applicability of the method.
  • a 1:100 model of a 4600-ton tanker was used to simulate the interaction between the dam break wave and the 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 break.
  • the initial configuration is shown in Figure 3.
  • the width and length of the water tank are respectively about 3 times the width and length of the ship.
  • the profile and surface discretization of the ship model are shown in Figure 4.
  • the particle spacing on the surface of the ship model is roughly the same as the resolution of the flow field, that is, 0.02m (about 1/90 of the length of the ship model).
  • the total number of particles is 289417, of which 193284 are fluid particles.
  • the Myklestad's method is used to calculate the mode function of an elastic non-uniform beam.
  • the results of the first three-order mode function are shown in Figure 6.
  • Young's modulus E 2 ⁇ 10 4 N/m 2
  • Figure 7 Some typical time steps of free surface shape and ship deformation are shown in Figure 7.
  • the particles of the free surface, the particles of the part of the ship model above the free surface, and the particles of the part of the ship model below the free surface are marked as red, blue, and black, respectively. It can be seen from Figure 7 that there are no wrong free surface particles identified below the free surface. In this simulation, the deformation amplitude is significantly larger than the normal ship model structure. The purpose of choosing a lower stiffness is to test the coupling ability between the modified modal superposition and the MPS model.
  • Figure 7 (a), (b), (c) shows the shock wave motion process at different moments of the dam break.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Mathematical Analysis (AREA)
  • Fluid Mechanics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Chemical & Material Sciences (AREA)
  • Combustion & Propulsion (AREA)
  • Mechanical Engineering (AREA)
  • Ocean & Marine Engineering (AREA)
  • Automation & Control Theory (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Databases & Information Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于改进的移动粒子半隐式法和模态叠加方法求解强非线性时域水弹性问题设计方法，其基于改进的移动粒子半隐式法(Moving Particle Semi-implicit,MPS)的流体动力求解器、基于耦合刚体与弹性模态的模态叠加方法的结构响应求解器，二者采用迭代的方式在时域实现强耦合。本发明改进了传统的基于势流的水弹性计算无法直接模拟自由表面的剧烈变形和船体砰击现象的不足，并且考虑了结构的刚性与柔性运动模态之间的相互影响问题。

Description

一种基于改进的移动粒子半隐式法和模态叠加方法求解强非线性时域水弹性问题设计方法 技术领域
本发明涉及一种基于改进的移动粒子半隐式法(MPS)和模态叠加方法求解强非线性时域水弹性问题设计方法，属于船舶工程领域。
背景技术
船舶和海洋结构物的体积越来越大、速度越来越快的趋势使得结构的固有频率和遇到的波浪频率趋于相同的范围。这意味着水弹性效应对于海洋结构的整体疲劳分析越来越重要。
这种情况下的结构完整性分析需要时域仿真。传统的基于势流的求解器无法直接模拟砰击过程中自由表面的剧烈变形，这就要求基于CFD的求解器进行流动计算，部分的结合或全面的应用基于粘流的CFD模型是水弹性计算的终极目标。与传统的基于网格的CFD方法相比，SPH、MPS等粒子法由于不存在网格畸变问题，且流场采用拉格朗日方法更新，非常适合模拟大的自由表面变形。移动粒子半隐式法(Moving Particle Semi-implicit,MPS)作为主流粒子方法之一，已经被改进并应用于各种剧烈流体和结构相互作用问题。此外，船舶的典型运动是大刚体运动加上相对较小的弹性变形。这意味着模态叠加理论足以准确计算弹性变形。然而，在传统的水弹性分析中，不考虑刚体与弹性模态之间的相互作用，刚体和弹性部分仅在交叉项的水动力系数意义上相互联系。
发明内容
根据上述提出的技术问题，而提供一种基于改进的移动粒子半隐式法和模态叠加方法求解强非线性时域水弹性问题设计方法。本发明采用的技术手段如下：
一种基于改进的移动粒子半隐式法和模态叠加方法求解强非线性时域水弹性问题设计方法，在时域上将求解分解为若干间隔很小的时间段，假定前一时间步(或初始时刻)流场和结构运动的信息已知，对当前时间步实施如下的计算，具体步骤如下：
S1、用无网格离散点对流场和边界进行空间离散，避免了强非线性自由面运动和结构大刚体位移导致的网格处理困难，若前一时间步为初始时刻，则还需对流体粒子的速度和压强进行初始化；
S2、根据前一时间步固体边界的运动信息，确定当前时间步的结构边界位置的估计值，若前一时间步为初始时刻，则采用结构运动的初始值；
S3、不计流体压力，基于前一时间步的流场信息，按照Navier-Stokes方程更新流体粒子的位置和速度；
S4、推导出压力泊松方程并求解压力泊松方程；
S5、通过求解压力泊松方程得到的压力值更新流体粒子的位置和速度；
S6、通过求解压力泊松方程得到的压力值，求解耦合刚体和弹性模态的固体响应控制方程；
S7、利用步骤S6得到的动力学信息以及固体结构的各阶弹性模态并通过线性叠加的方法计算出结构总体的弹性变形更新结构运动和变形，确定新的流体的边界位置；
S8、比较步骤S2和步骤S7中的结构位置差值，若满足收敛条件则进行下一时间步运算，否则在当前时间步重复步骤S2到步骤S8，直至收敛条件满足后进行下一时间步的计算。
所述步骤S1的具体步骤如下：
在MPS方法中，物理量的梯度、散度和拉普拉斯算子采用加权平均方法，基于均匀分布的粒子进行离散，如式(1.1.1)所示。
Figure PCTCN2019115597-appb-000001
其中φ和Φ分别表示任意的标量和矢量，d m表示所研究问题的维数：二维或者三维，N为相关粒子影响区域内的粒子数；n 0为初始时刻的粒子密度，r表示粒子间距，r i和r j分别表示粒子i和粒子j的坐标r ij＝|r i-r j|；
核函数W(r ij)和参数λ定义如式(1.1.2)和式(1.1.3)：
Figure PCTCN2019115597-appb-000002
Figure PCTCN2019115597-appb-000003
自由表面粒子)技术，在这种技术中，是直接调整粒子的位移，而不是通过施加一个力来调整的。在每一个时间步后，对粒子的位置进行微小的位移，使其分布规律。这种技术也可以看作是一种网格重构过程。此外，由于压力泊松方程是根据不可压缩性条件推导出来的，由此产生的压力将使相邻粒子之间的距离大致保持在相同的值附近(即初始粒子距离r 0)。
粒子间距的调整量如式(1.1.5)所示：
Figure PCTCN2019115597-appb-000004
当|r ij|≤r 0时(1.1.5)，其中，r 0为初始粒子间距，
对于远离主流体域的自由表面粒子，当粒子之间的相对速度比预测步骤S3之前的阈值更接近时，进行如式(1.1.6)的操作，将粒子对相对速度δu i设置为零。
Figure PCTCN2019115597-appb-000005
当
Figure PCTCN2019115597-appb-000006
r min＝0.3r 0时   (1.1.6)
式中，
Figure PCTCN2019115597-appb-000007
为粒子i与粒子j的切向相对速度，固体与流体相邻粒子分别取0.5和1.0，Δt为时间步的大小。
所述步骤S2的具体步骤如下：引入描述物体运动和变形的变量A，其定义如下：
A＝[X Rb,θ,q b]      (2.1.1)
式中，X Rb代表刚体运动的质心位置，θ表示结构的底升角，q b表示模态函数对应的一般坐标。
结构与流体交界面的粒子的相应的位置、速度和加速度为
Figure PCTCN2019115597-appb-000008
和
Figure PCTCN2019115597-appb-000009
其中k表示时间步，下标‘0’所在的位置表示迭代的次数；
由前一时间步结构边界粒子的位置、速度、加速度即
Figure PCTCN2019115597-appb-000010
和
Figure PCTCN2019115597-appb-000011
按照匀加速运动的假设估算出当前时间步的位置、速度和加速度，并以此刻计算出的位置、速度、加速度计算出当前时间步的结构边界位置的估计值。
步骤S3中，Navier-Stokes方程为基于拉格朗日形式的不可压缩无粘Navier-Stokes方程，其如式(3.1.1)所示：
Figure PCTCN2019115597-appb-000012
其中，u，p和ρ分别表示流体的速度、压力和密度，均为指向重力方向的向量，g表示重力加速度；
所述步骤S3的具体步骤如下：
采用经典投影法对速度和压力解耦如下：
在不考虑压力的情况，仅通过惯性将流体推动到中间状态，如式(3.1.2)所示：
Figure PCTCN2019115597-appb-000013
其中，l为位置向量，下标带*和k的变量分别表示中间时间步和第k个时间步的变量，u表示流体的速度。
所述步骤S4的具体步骤如下：
根据连续性方程对Navier-Stokes方程两侧取散度可推导出压力泊松方程，如式(4.1.2)所示，其中，连续性方程如式(4.1.1)所示，
Figure PCTCN2019115597-appb-000014
Figure PCTCN2019115597-appb-000015
其中n 0和n k分别为初始时刻和第k个时间步的粒子密度，粒子密度定义如式(1.1.4)所示；
系数α定义如式(4.1.3)所示：
Figure PCTCN2019115597-appb-000016
该α模型不需要根据不同的问题进行特殊的标定，与其他相似的模型相比，具有更强的稳定性。
求解压力泊松方程的边界条件如下：
(1)固体边界条件
将Neumann边界条件应用于固体粒子，如式(4.1.4)所示，
Figure PCTCN2019115597-appb-000017
式中
Figure PCTCN2019115597-appb-000018
和
Figure PCTCN2019115597-appb-000019
为固体边界粒子在第k和(k+1)个时间步处的加速度；由于求解第(k+1)个时间步的流体运动时
Figure PCTCN2019115597-appb-000020
是未知的，所以用上一个时间步的值
Figure PCTCN2019115597-appb-000021
作为近似；
靠近固体边界的流体粒子需要对拉普拉斯算子进行修正，使其与固体边界上的Neumann边界条件相一致，以弥补相邻粒子的不足，其压力值如式(4.1.5)所示，考虑固体边界外虚拟粒子的压力值，对式(1.1.1)中的离散项进行修正，使其与式(4.1.2)一致，对近边界流体粒子来说，在压力泊松方程的源项中，也包括在影响域的固体边界粒子的作用，因此，由于相同的一致性原因，受影响固体边界粒子的中间速度被修正为式(4.1.6)-(1)，其分别投影到切向τ和法向n方向得到式(4.1.6)-(2)和式(4.1.6)-(3)，
Figure PCTCN2019115597-appb-000022
Figure PCTCN2019115597-appb-000023
其中p v表示虚粒子的压强，p s表示响应固体粒子的压强，u b*和
Figure PCTCN2019115597-appb-000024
为固体边界粒子在中间时刻和第(k+1)个时间步长时的速度，公式(4.1.5)和(4.1.6)中n均表示法向；
(2)自由表面边界条件
由于求解压力泊松方程需要在自由面边界施加压力为零的边界条件，因此需要识别自由液面粒子位置：
对于二维情况，为每个粒子指定一个以自身为中心的圆，半径为1.05r 0，并用360个均匀分布在圆上的点离散，若所有这些点全部被其相邻粒子的圆所覆盖，则该粒子为内流体粒子，否则该粒子为将被识别为自由表面粒子；检查过程是通过遍历每个相邻粒子，识别每个相邻粒子所覆盖的弧(或3D情况下的patch)上的重叠点。
对于三维的情况，采用两步方法，首先，利用公式(4.1.7)检查相邻粒子数，检测出潜在的自由表面粒子：N p＜β 3dN p0    (4.1.7)
式中，β 3d是一个小于1的参数，根据情况取值，N p0表示圈定的范围内初始粒子密度；
其次，通过建立以粒子为中心、半径为1.05r 0的球面来细化搜索，具体为：通过加权平均法确定一个指向该粒子周围粒子分布最稀疏方向的向量，通过均匀分布的点来离散球面上圆形区域，若所有这些点全部被其相邻粒子的球面所覆盖，则该粒子为内流体粒子，否则该粒子为将被识别为自由表面粒子。
所述步骤S5的计算公式如式(5.1.1)所示，
Figure PCTCN2019115597-appb-000025
典型海洋结构对波浪冲击的动力响应可以表现为刚体运动较大，弹性变形较小，这意味着模态叠加法能够处理弹性运动部分。本发明基于拉格朗日方程，推导了耦合刚体和弹性模态的时域模型，所述步骤S6的具体步骤如下：
S6.1、结构运动描述模型建立
结构运动描述模型采用固定整体坐标系XOY和附着体局部坐标系soη，其中，附着体局部坐标系的原点选择在结构运动描述模型的梁的重心处，因此，结构运动描述模型的梁上任意点的位置X可由式(6.1.1)描述：X＝X Rb+Rξ T    (6.1.1)
其中，X Rb为结构运动描述模型的梁的质心位置，式(6.1.2)和式(6.1.3)中定义了结构运动描述模型的梁的局部坐标ξ和结构运动描述模型的梁的旋转矩阵R，
ξ＝[s,η l+η 0]    (6.1.2)
Figure PCTCN2019115597-appb-000026
其中，θ为O-X轴与o-s轴逆时针的夹角，η l为未变形结构的坐标，对于结构运动描述模型的梁上的一个特定点是常数，η 0为中平面的挠度，任意平行于中平面的面的挠度均取中平面的相同值，根据模态叠加理论，将中平面的挠度η 0展开为式(6.1.4)：
Figure PCTCN2019115597-appb-000027
式(6.1.5)和(6.1.6)中定义了模态函数的列向量
Figure PCTCN2019115597-appb-000028
及对应的一般坐标q b，
Figure PCTCN2019115597-appb-000029
q b＝[q b1,q b2,q b3,...] T    (6.1.6)
模态函数满足以下正交关系：
Figure PCTCN2019115597-appb-000030
I u是单位矩阵    (6.1.7)
Figure PCTCN2019115597-appb-000031
Λ＝diag(ω i 2)i＝1,2,3,...    (6.1.8)
其中s 1和s 2为沿o-s轴积分的上下限，ρ l为结构运动描述模型的梁的线密度，E表示结构的弹性模量，I I为结构的惯性矩阵；
S6.2、基于拉格朗日方程的固体响应控制方程建立
该改进的模态叠加法模型的主要特点是刚体模态和弹性模态在时域内耦合。考虑了刚体运动的惯性效应等相互作用对弹性变形的影响。这是对传统水弹性理论的一种扩展，在传统水弹性理论中，刚体和弹性部件只通过对应的交叉项的水动力系数进行连接。
基于拉格朗日方程建立模型：
Figure PCTCN2019115597-appb-000032
其中T和U分别表示结构运动描述模型的动能和势能，q j为任意刚柔模态对应的一般坐标，Q j为第j个坐标对应的非保守力；
结构运动描述模型的梁为弹性的非均匀梁，弹性的非均匀梁的T、U、q j、Q j的具体形式如下：
Figure PCTCN2019115597-appb-000033
Figure PCTCN2019115597-appb-000034
式中，ρ s表示线密度，M表示结构的质量，X R和Y R分别对应结构质心的X和Y坐标，g表示重力加速度，η 1、η 2为沿轴积分η、矩阵积分U、列向量积分ψ 0和ψ 1的上下限，矩阵积分U和列向量积分ψ 0、ψ 1定义在式(6.2.4)和(6.2.6)中：
Figure PCTCN2019115597-appb-000035
Figure PCTCN2019115597-appb-000036
Figure PCTCN2019115597-appb-000037
将式(6.2.2)和(6.2.3)代入式(6.2.1)，可得在波浪冲击作用下弹性的非均匀梁的动力响应控制方程，即得到耦合刚体和弹性模态的固体响应控制方程，如式(6.2.7)所示：
Figure PCTCN2019115597-appb-000038
式(6.2.8)-(6.2.11)给出刚体和弹性模态的一般坐标对应的非保守力：
Figure PCTCN2019115597-appb-000039
Figure PCTCN2019115597-appb-000040
Figure PCTCN2019115597-appb-000041
Figure PCTCN2019115597-appb-000042
式中，n s、n η分别为s和η的局部法向分量，用固定整体坐标系表示，X p和Y p是压强p作用点的整体X和Y坐标，e η是o-η轴的单位向量，积分是围绕结构运动描述模型的梁的周长进行的，由式(6.2.7)可以看出，刚体与弹性模态之间的相互作用即交叉项是由角运动引起的。
S6.3、求解耦合刚体和弹性模态的固体响应控制方程。
对于流体与结构的相互作用过程，采用标准迭代法，直到满足一定的准则。
所述步骤S8的具体步骤如下：
假设所有的流体和结构变量在时间步t＝t k-1，而后，在下一时间步t＝t k的交互详细过程如下：
S8.1、假设结构在时间步t＝t k的加速度与上一时间步t＝t k-1的加速度是相同的，即
Figure PCTCN2019115597-appb-000043
则在时间步t＝t k的位置D k和速度
Figure PCTCN2019115597-appb-000044
也可以通过有限差分法计算出来；
通过式(8.1.1)(8.1.2)和式(8.1.3)～(8.1.6)计算出流体和结构交界面上每个点的相应的位置、速度和加速度，即
Figure PCTCN2019115597-appb-000045
和
Figure PCTCN2019115597-appb-000046
可以从新计算出来的结构运动和变形参数A计算得到，即式(2.1.1)，
Figure PCTCN2019115597-appb-000047
Figure PCTCN2019115597-appb-000048
Figure PCTCN2019115597-appb-000049
Figure PCTCN2019115597-appb-000050
Figure PCTCN2019115597-appb-000051
Figure PCTCN2019115597-appb-000052
其中，X R为刚体的全局坐标，X CR为质心坐标，R R是将附着体局部坐标系soη与固定整体坐标系XOY联系起来的旋转矩阵；X fi表示中心线上的坐标；ξ i表示结构运动描述模型的梁在附着体局部坐标系soη下的中心线上的坐标，其表示方法为式(8.1.7)，式中的η i为结构运动描述模型的梁的挠度函数；
ξ i＝[s i,η i]    (8.1.7)
S8.2、利用交界面新的动力学信息作为新的边界条件，根据修正的MPS方法计算t＝t k时刻的流体运动，然后由此得到新的压力
Figure PCTCN2019115597-appb-000053
并将之应用于交界面在t＝t k时刻的第i次迭代；
S8.3、比较步骤S2和步骤S7中的结构位置差值
Figure PCTCN2019115597-appb-000054
和
Figure PCTCN2019115597-appb-000055
若满足收敛条件，如式8.3.1)所示，则转到步骤S8.1进行下一时间步t＝t k+1运算，
Figure PCTCN2019115597-appb-000056
否则，利用式(8.3.2)修正第i+1次迭代的结构位置
Figure PCTCN2019115597-appb-000057
并且利用Newmark方法，更新速度
Figure PCTCN2019115597-appb-000058
和加速度
Figure PCTCN2019115597-appb-000059
根据
Figure PCTCN2019115597-appb-000060
和
Figure PCTCN2019115597-appb-000061
计算交界面的变量
Figure PCTCN2019115597-appb-000062
和
Figure PCTCN2019115597-appb-000063
使用这些修正的交界面变量，通过返回步骤S2进行第i+1次迭代。
Figure PCTCN2019115597-appb-000064
其中，χ i为Aitken弛豫因子，其核心思想是通过使用前两次迭代来提高预测的准确性。其值由式(8.3.3)计算：
Figure PCTCN2019115597-appb-000065
其中，
Figure PCTCN2019115597-appb-000066
本发明公开了一种基于改进的移动粒子半隐式法和模态叠加方法求解强非线性时域水弹性问题设计方法，其基于改进的移动粒子半隐式法(Moving Particle Semi-implicit,MPS)的流体动力求解器、基于耦合刚体与弹性模态的模态叠加方法的结构响应求解器，二者采用迭代的方式在时域实现强耦合。本发明改进了传统的基于势流的水弹性计算无法直接模拟自由表面的剧烈变形和船体砰击现象的不足，并且考虑了结构的刚性与柔性运动模态之间的相互影响问题。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图做以简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为二维和三维自由表面粒子识别方案图。
图2为弹性的非均匀梁示意图。
图3为计算的初始配置示意图(不按比例)。
图4为船舶模型的剖面和表面离散化图。
图5为船舶模型的质量(a)和二阶矩(b)分布图。
图6为船舶模型前三阶振型图。
图7为典型时间步自由表面形状和船舶变形图((a)t＝0.1s,(b)t＝0.2s,(c)t＝0.418s,(d)t＝0.684s,(e)t＝0.91s,(f)t＝1.05s)。
图8为用于补偿固体边界附近的虚拟粒子的拉普拉斯算子的演示。
图9为流体与结构的相互作用迭代过程流程图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图1-图9所示，一种基于改进的移动粒子半隐式法和模态叠加方法求解强非线性时域水弹性问题设计方法，在时域上将求解分解为若干间隔很小的时间段，假定前一时间步(或初始时刻)流场和结构运动的信息已知，对当前时间步实施如下的计算，具体步骤如下：
S1、用无网格离散点对流场和边界进行空间离散，若前一时间步为初始时刻，则还需对流体粒子的速度和压强进行初始化；
S2、根据前一时间步固体边界的运动信息，确定当前时间步的结构边界位置的估计值，若前一时间步为初始时刻，则采用结构运动的初始值；
S3、不计流体压力，基于前一时间步的流场信息，按照Navier-Stokes方程更新流体粒子的位置和速度；
S4、推导出压力泊松方程并求解压力泊松方程；
S5、通过求解压力泊松方程得到的压力值更新流体粒子的位置和速度；
S6、通过求解压力泊松方程得到的压力值，求解耦合刚体和弹性模态的固体响应控制方程；
S7、利用步骤S6得到的动力学信息以及固体结构的各阶弹性模态并通过线性叠加的方法计算出结构总体的弹性变形更新结构运动和变形，确定新的流体的边界位置；
S8、比较步骤S2和步骤S7中的结构位置差值，若满足收敛条件则进行下一时间步运算，否则在当前时间步重复步骤S2到步骤S8，直至收敛条件满足后进行下一时间步的计算。
步骤S1-S4为流体动力求解器部分(修正的MPS方法)，步骤S6为结构响应求解器。
所述步骤S1的具体步骤如下：基于均匀分布的粒子进行离散，如式(1.1.1)所示。
Figure PCTCN2019115597-appb-000067
其中φ和Φ分别表示任意的标量和矢量，d m表示所研究问题的维数：二维或者三维，N为相关粒子影响区域内的粒子数，梯度/散度和拉普拉斯算子圆域半径r e取值为2.1r 0和4.0r 0，r 0为初始粒子间距；n 0为初始时刻的粒子密度；
核函数W(r ij)和参数λ定义如式(1.1.2)和式(1.1.3)：
Figure PCTCN2019115597-appb-000068
Figure PCTCN2019115597-appb-000069
粒子密度定义在式(1.1.4)中:
Figure PCTCN2019115597-appb-000070
粒子间距的调整量如式(1.1.5)所示：
Figure PCTCN2019115597-appb-000071
当|r ij|≤r 0时    (1.1.5)
对于远离主流体域的自由表面粒子，当粒子之间的相对速度比预测步骤之前的阈值更接近时，进行如式(1.1.6)的操作，将粒子对相对速度设置为零。
Figure PCTCN2019115597-appb-000072
当
Figure PCTCN2019115597-appb-000073
r min＝0.3r 0时，(1.1.6)
式中
Figure PCTCN2019115597-appb-000074
为粒子i与粒子j的切向相对速度，固体与流体相邻粒子分别取0.5和1.0。
所述步骤S2的具体步骤如下：
引入描述物体运动和变形的变量A，其定义如下：A＝[X Rb,θ,q b]   (2.1.1)
结构与流体交界面的粒子的相应的位置、速度和加速度为
Figure PCTCN2019115597-appb-000075
和
Figure PCTCN2019115597-appb-000076
其中k表示时间步，下标‘0’所在的位置表示迭代的次数；
由前一步结构边界粒子的位置、速度、加速度即
Figure PCTCN2019115597-appb-000077
和
Figure PCTCN2019115597-appb-000078
按照匀加速运动的假设估算出当 前时间步的位置、速度和加速度，并以此刻计算出的位置、速度、加速度计算出当前时间步的结构边界位置的估计值。
Navier-Stokes方程为基于拉格朗日形式的不可压缩无粘Navier-Stokes方程，其如式(3.1.1)所示：
Figure PCTCN2019115597-appb-000079
其中，u，p和ρ分别表示流体速度、压力和密度，均为指向重力方向的向量，g表示重力加速度；
所述步骤S3的具体步骤如下：
采用经典投影法对速度和压力解耦如下：在不考虑压力的情况，仅通过惯性将流体推动到中间状态，如式(3.1.2)所示：
Figure PCTCN2019115597-appb-000080
其中，l为位置向量，下标带*和k的变量分别表示中间时间步和第k个时间步的变量。
所述步骤S4的具体步骤如下：
根据连续性方程对Navier-Stokes方程两侧取散度可推导出压力泊松方程，如式(4.1.2)所示，其中，连续性方程如式(4.1.1)所示，
Figure PCTCN2019115597-appb-000081
Figure PCTCN2019115597-appb-000082
其中n 0和n k分别为初始时刻和第k个时间步的粒子密度，粒子密度定义如式(1.1.4)所示；
系数α定义如式(4.1.3)所示：
Figure PCTCN2019115597-appb-000083
求解压力泊松方程的边界条件如下：
(1)固体边界条件
将Neumann边界条件应用于固体粒子，如式(4.1.4)所示，
Figure PCTCN2019115597-appb-000084
式中
Figure PCTCN2019115597-appb-000085
和
Figure PCTCN2019115597-appb-000086
为固体边界粒子在第k和(k+1)个时间步处的加速度；由于求解第(k+1)个时间步的流体运动时
Figure PCTCN2019115597-appb-000087
是未知的，所以用上一个时间步的值
Figure PCTCN2019115597-appb-000088
作为近似；
图8所示，靠近固体边界的流体粒子需要对拉普拉斯算子进行修正，其压力值如式(4.1.5)所示，对式(1.1.1)中的离散项进行修正，使其与式(4.1.2)一致，受影响固体边界粒子的中间速度被修正为式(4.1.6)-(1)，其分别投影到切向τ和法向n方向得到式(4.1.6)-(2)和式(4.1.6)-(3)，
Figure PCTCN2019115597-appb-000089
Figure PCTCN2019115597-appb-000090
其中u b*和
Figure PCTCN2019115597-appb-000091
为固体边界粒子在中间时刻和第(k+1)个时间步长时的速度；
(2)自由表面边界条件
由于求解压力泊松方程需要在自由面边界施加压力为零的边界条件，因此需要识别自由液面粒子位置：
对于二维情况，为每个粒子指定一个以自身为中心的圆，半径为1.05r 0，并用360个均匀分布在圆上的点离散，若所有这些点全部被其相邻粒子的圆所覆盖，则该粒子为内流体粒子(图1(a)中的粒子B)，否则该粒子为将被识别为自由表面粒子(图1(a)中的粒子A)；
对于三维的情况，采用两步方法，首先，利用公式(4.1.7)检查相邻粒子数，检测出潜在的自由表面粒子，
N p＜β 3dN p0    (4.1.7)
在本实施例中N p0＝32,β 3d＝0.96。
其次，通过建立以粒子为中心、半径为1.05r 0的球面来细化搜索，具体为：通过加权平均法确定一个指向该粒子周围粒子分布最稀疏方向的向量，如图1(b)所示，通过均匀分布的点来离散球面上圆形区域，若所有这些点全部被其相邻粒子的球面所覆盖，则该粒子为内流体粒子，否则该粒子为将被识别为自由表面粒子。
所述步骤S5的计算公式如式(5.1.1)所示，
Figure PCTCN2019115597-appb-000092
所述步骤S6的具体步骤如下：
S6.1、结构运动描述模型建立
结构运动描述模型采用固定整体坐标系XOY和附着体局部坐标系soη，其中，附着体局部坐标系的原点选择在结构运动描述模型的梁的重心处，因此，结构运动描述模型的梁上任意点的位置可由式(6.1.1)描述，示意图如图2所示。X＝X Rb+Rξ T    (6.1.1)
其中，X Rb为的质心位置，式(6.1.2)和式(6.1.3)中定义了结构运动描述模型的梁的局部坐标ξ和结构运动描述模型的梁的旋转矩阵R，
ξ＝[s,η l+η 0]    (6.1.2)
Figure PCTCN2019115597-appb-000093
其中，θ为O-X轴与o-s轴逆时针的夹角，η l为未变形结构的坐标，对于结构运动描述模型的梁上的一个特定点是常数，η 0为中平面的挠度，任意平行于中平面的面的挠度均取中平面的相同值，根据模态叠加理论，将中平面的挠度η 0展开为式(6.1.4)，
Figure PCTCN2019115597-appb-000094
式(6.1.5)和(6.1.6)中定义了模态函数的列向量
Figure PCTCN2019115597-appb-000095
及对应的一般坐标q b，
Figure PCTCN2019115597-appb-000096
q b＝[q b1,q b2,q b3,...] T    (6.1.6)
模态函数满足以下正交关系：
Figure PCTCN2019115597-appb-000097
I u是单位矩阵    (6.1.7)
Figure PCTCN2019115597-appb-000098
Λ＝diag(ω i 2)i＝1,2,3,...    (6.1.8)
其中s 1和s 2为沿o-s轴积分的上下限，ρ l为结构运动描述模型的梁的线密度，E表示结构的弹性模量，I I为结构的惯性矩；
S6.2、基于拉格朗日方程的固体响应控制方程建立
基于拉格朗日方程建立模型：
Figure PCTCN2019115597-appb-000099
其中T和U分别表示结构运动描述模型的动能和势能，q j为任意刚柔模态对应的一般坐标，Q j为第j个坐标对应的非保守力；
结构运动描述模型的梁为弹性的非均匀梁，弹性的非均匀梁的T、U、q j、Q j的具体形式如下：
Figure PCTCN2019115597-appb-000100
Figure PCTCN2019115597-appb-000101
式中η 1、η 2为沿轴积分η、矩阵积分U、列向量积分ψ 0和ψ 1的上下限，矩阵积分U和列向量积分ψ 0、ψ 1定义在式(6.2.4)和(6.2.6)中：
Figure PCTCN2019115597-appb-000102
Figure PCTCN2019115597-appb-000103
Figure PCTCN2019115597-appb-000104
将式(6.2.2)和(6.2.3)代入式(6.2.1)，可得在波浪冲击作用下弹性的非均匀梁的动力响应控制方程，即得到耦合刚体和弹性模态的固体响应控制方程，如式(6.2.7)所示：
Figure PCTCN2019115597-appb-000105
式(6.2.8)-(6.2.11)给出刚体和弹性模态的一般坐标对应的非保守力：
Figure PCTCN2019115597-appb-000106
Figure PCTCN2019115597-appb-000107
Figure PCTCN2019115597-appb-000108
Figure PCTCN2019115597-appb-000109
式中，n s、n η分别为s和η的局部法向分量，用固定整体坐标系表示，X p和Y p是压强p作用点的整体X和Y坐标，e η是o-η轴的单位向量，积分是围绕结构运动描述模型的梁的周长进行的。
S6.3、求解耦合刚体和弹性模态的固体响应控制方程。
所述步骤S8的具体步骤如下：
假设所有的流体和结构变量在时间步t＝t k-1，而后，在下一时间步t＝t k的交互详细过程如下：
S8.1、假设结构在时间步t＝t k的加速度与上一时间步t＝t k-1的加速度是相同的，即
Figure PCTCN2019115597-appb-000110
则在时间步t＝t k的位置D k和速度
Figure PCTCN2019115597-appb-000111
也可以通过有限差分法计算出来；
通过式(8.1.1)(8.1.2)和式(8.1.3)～(8.1.6)计算出流体和结构交界面上每个点的相应的位置、速度和加速度，即
Figure PCTCN2019115597-appb-000112
和
Figure PCTCN2019115597-appb-000113
可以从新计算出来的结构运动和变形参数A计算得到，即式(2.1.1)，
Figure PCTCN2019115597-appb-000114
Figure PCTCN2019115597-appb-000115
Figure PCTCN2019115597-appb-000116
Figure PCTCN2019115597-appb-000117
Figure PCTCN2019115597-appb-000118
Figure PCTCN2019115597-appb-000119
其中，X R为刚体的全局坐标，X CR为质心坐标，R R是将附着体局部坐标系soη与固定整体坐标系XOY联系起来的旋转矩阵；X fi表示中心线上的坐标；ξ i表示结构运动描述模型的梁在附着体局部坐标系soη下的中心线上的坐标，其表示方法为式(8.1.7)，式中的η i为结构运动描述模型的梁的挠度函数；
ξ i＝[s i,η i]    (8.1.7)
S8.2、利用交界面新的动力学信息作为新的边界条件，根据修正的MPS方法计算t＝t k时刻的流体运动，然后由此得到新的压力
Figure PCTCN2019115597-appb-000120
并将之应用于交界面在t＝t k时刻的第i次迭代；
S8.3、比较步骤S2和步骤S7中的结构位置差值
Figure PCTCN2019115597-appb-000121
和
Figure PCTCN2019115597-appb-000122
若满足收敛条件，如式8.3.1)所示，则转到步骤S8.1进行下一时间步t＝t k+1运算，
Figure PCTCN2019115597-appb-000123
否则，利用式(8.3.2)修正第i+1次迭代的结构位置
Figure PCTCN2019115597-appb-000124
并且利用Newmark方法，更新速度
Figure PCTCN2019115597-appb-000125
和加速度
Figure PCTCN2019115597-appb-000126
根据
Figure PCTCN2019115597-appb-000127
和
Figure PCTCN2019115597-appb-000128
计算交界面的变量
Figure PCTCN2019115597-appb-000129
和
Figure PCTCN2019115597-appb-000130
使用这些修正的交界面变量，通过返回步骤S2进行第i+1次迭代。
Figure PCTCN2019115597-appb-000131
其中，χ i为Aitken弛豫因子，其核心思想是通过使用前两次迭代来提高预测的准确性。其值由式(8.3.3)计算：
Figure PCTCN2019115597-appb-000132
其中，
Figure PCTCN2019115597-appb-000133
每次第一次迭代的初始值为0.2。
整个过程如图9所示。
下面通过基于改进的移动粒子半隐式法和模态叠加方法求解强非线性时域水弹性问题设计方法对三维船体砰击的问题进行计算，说明该方法的适用性。
具体计算过程：
1.模型选择
采用某4600吨级油轮1:100模型，模拟了溃坝波与船舶的相互作用。假设船舶动力响应为二维(即只考虑俯仰、升沉和垂直弯曲)，这使得弹性的非均匀梁适用于这种情况。船舶模型的主要参数如表1所示。船舶模型放置在面临溃坝冲击的水池中，初始配置如图3所示。水箱的宽度和长度分别为船舶宽度和长度的3倍左右。
表1 船舶模型主要参数
Figure PCTCN2019115597-appb-000134
2.粒子离散化
船舶模型的剖面和表面离散化如图4所示。船舶模型表面的粒子间距与流场分辨率大致相同，即0.02m(约为船舶模型长度的1/90)。总粒子数为289417，其中流体粒子193284个。
3.参数选择
假设船舶模型的质量分布与其占水量相同，截面二阶矩的计算选用船舶模型厚度为0.05m。图5为船舶模型的质量和二阶矩的分布。
采用Myklestad’s方法计算弹性的非均匀梁的振型函数，前三阶振型函数的结果如图6所示。杨氏模量E＝2×10 4N/m 2，和前三的圆形固有频率模式是：ω 1＝3.7857,ω 2＝9.7777,ω 3＝18.2063。
4.计算结果
某些典型时间步自由表面形状和船舶变形图如图7所示。自由表面粒子、自由表面上方的船舶模型部分的粒子和自由表面下方的船舶模型部分的粒子分别标记为红色、蓝色和黑色。从图7中可以看出，在自由表面下方没有识别错误的自由表面粒子。在此仿真中，变形幅度明显大于正常船舶模型结构。特意选择较低的刚度是为了测试修正模态叠加与MPS模型之间的耦合能力。图7(a)、(b)、(c)为溃坝不同时刻的冲击波运动过程。船舶模型吃水保持不变的事实证明，MPS计算的浮力对船舶模型重量的支撑是准确稳定的。波浪对船头的冲击发生在图7(c)中t＝0.418s处。当波浪沿船舶模型传播时，波浪的波峰处产生最大的弯矩(图7(d)-(f))，船舶模型相应地产生最大的变形。定性分析表明，修正的MPS与所提结构模型的耦合能够为波浪与结构相互作用的模拟提供合理的结果。
最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (8)

  1. 一种基于改进的移动粒子半隐式法和模态叠加方法求解强非线性时域水弹性问题设计方法，其特征在于具体步骤如下：
    S1、用无网格离散点对流场和边界进行空间离散，若前一时间步为初始时刻，则还需对流体粒子的速度和压强进行初始化；
    S2、根据前一时间步固体边界的运动信息，确定当前时间步的结构边界位置的估计值，若前一时间步为初始时刻，则采用结构运动的初始值；
    S3、不计流体压力，基于前一时间步的流场信息，按照Navier-Stokes方程更新流体粒子的位置和速度；
    S4、推导出压力泊松方程并求解压力泊松方程；
    S5、通过求解压力泊松方程得到的压力值更新流体粒子的位置和速度；
    S6、通过求解压力泊松方程得到的压力值，求解耦合刚体和弹性模态的固体响应控制方程；
    S7、利用步骤S6得到的动力学信息以及固体结构的各阶弹性模态并通过线性叠加的方法计算出的结构总体的弹性变形来更新结构运动和变形，确定新的流体的边界位置；
    S8、比较步骤S2和步骤S7中的结构位置差值，若满足收敛条件则进行下一时间步运算，否则在当前时间步重复步骤S2到步骤S8，直至收敛条件满足后进行下一时间步的计算。
  2. 根据权利要求1所述的方法，其特征在于：所述步骤S1的具体步骤如下：
    基于均匀分布的粒子进行离散，如式(1.1.1)所示：
    Figure PCTCN2019115597-appb-100001
    其中φ和Φ分别表示任意的标量和矢量，d m表示所研究问题的维数：二维或者三维，N为相关粒子影响区域内的粒子数；n 0为初始时刻的粒子密度，r表示粒子间距，r i和r j分别表示粒子i和粒子j的坐标r ij＝|r i-r j|；
    核函数W(ri j)和参数λ定义如式(1.1.2)和式(1.1.3)：
    Figure PCTCN2019115597-appb-100002
    Figure PCTCN2019115597-appb-100003
    其中，r e表示粒子相互作用的有效范围，粒子密度n定义在式(1.1.4)中:
    Figure PCTCN2019115597-appb-100004
    粒子间距r的调整量如式(1.1.5)所示：
    Figure PCTCN2019115597-appb-100005
    当|r ij|≤r 0时  (1.1.5)
    其中，r 0为初始粒子间距；
    对于远离主流体域的自由表面粒子，当粒子之间的相对速度比预测步骤S3之前的阈值更接近时，进行如式(1.1.6)的操作，将粒子对相对速度δu i设置为零；
    Figure PCTCN2019115597-appb-100006
    当
    Figure PCTCN2019115597-appb-100007
    时,r min＝0.3r 0  (1.1.6)
    式中
    Figure PCTCN2019115597-appb-100008
    为粒子i与粒子j的切向相对速度，固体与流体相邻粒子分别取0.5和1.0；Δt为时间步的大小。
  3. 根据权利要求2所述的方法，其特征在于：所述步骤S2的具体步骤如下：
    引入描述物体运动和变形的变量A，其定义如下：
    A＝[X Rb,θ,q b]  (2.1.1)
    式中，X Rb代表刚体运动的质心位置，θ表示结构的底升角，q b表示模态函数对应的一般坐标；
    结构与流体交界面的粒子的相应的位置、速度和加速度为
    Figure PCTCN2019115597-appb-100009
    和
    Figure PCTCN2019115597-appb-100010
    其中k表示时间步，下标‘0’所在的位置表示迭代的次数；
    由前一时间步结构边界粒子的位置、速度、加速度即
    Figure PCTCN2019115597-appb-100011
    和
    Figure PCTCN2019115597-appb-100012
    按照匀加速运动的假设估算出当前时间步的位置、速度和加速度，并以计算出的当前时间步的位置、速度、加速度计算出当前时间步的结构边界位置的估计值。
  4. 根据权利要求3所述的方法，其特征在于：步骤S3中，Navier-Stokes方程为基于拉格朗日形式的不可压缩无粘Navier-Stokes方程，其如式(3.1.1)所示：
    Figure PCTCN2019115597-appb-100013
    其中，u，p和ρ分别表示流体的速度、压力和密度，均为指向重力方向的向量，g表示重力加速度；
    所述步骤S3的具体步骤如下：
    采用经典投影法对速度和压力解耦如下：
    在不考虑压力的情况，仅通过惯性将流体推动到中间状态，如式(3.1.2)所示：
    Figure PCTCN2019115597-appb-100014
    其中，l为位置向量，下标带*和k的变量分别表示中间时间步和第k个时间步的变量，u表示流体的速度。
  5. 根据权利要求4所述的方法，其特征在于：所述步骤S4的具体步骤如下：
    根据连续性方程对Navier-Stokes方程两侧取散度可推导出压力泊松方程，如式(4.1.2)所示，其中，连续性方程如式(4.1.1)所示，
    Figure PCTCN2019115597-appb-100015
    Figure PCTCN2019115597-appb-100016
    其中n 0和n k分别为初始时刻和第k个时间步的粒子密度，粒子密度定义如式(1.1.4)所示；
    系数α定义如式(4.1.3)所示：
    Figure PCTCN2019115597-appb-100017
    求解压力泊松方程的边界条件如下：
    (1)固体边界条件
    将Neumann边界条件应用于固体粒子，如式(4.1.4)所示，
    Figure PCTCN2019115597-appb-100018
    式中
    Figure PCTCN2019115597-appb-100019
    和
    Figure PCTCN2019115597-appb-100020
    为固体边界粒子在第k和(k+1)个时间步处的加速度；由于求解第(k+1)个时间步的流体运动时
    Figure PCTCN2019115597-appb-100021
    是未知的，所以用上一个时间步的值
    Figure PCTCN2019115597-appb-100022
    作为近似；
    靠近固体边界的流体粒子需要对拉普拉斯算子进行修正，其压力值如式(4.1.5)所示，对式(1.1.1)中的离散项进行修正，使其与式(4.1.2)一致，受影响固体边界粒子的中间速度被修正为式(4.1.6)-(1)，其分别投影到切向τ和法向n方向得到式(4.1.6)-(2)和式(4.1.6)-(3)，
    Figure PCTCN2019115597-appb-100023
    Figure PCTCN2019115597-appb-100024
    其中p v表示虚粒子的压强，p s表示响应固体粒子的压强，u b*和
    Figure PCTCN2019115597-appb-100025
    为固体边界粒子在中间时刻和第(k+1)个时间步长时的速度，公式(4.1.5)和(4.1.6)中n均表示法向；
    (2)自由表面边界条件
    由于求解压力泊松方程需要在自由面边界施加压力为零的边界条件，因此需要识别自由液面粒子位置：
    对于二维情况，为每个粒子指定一个以自身为中心的圆，半径为1.05r 0，并用360个均匀分布在圆上的点离散，若所有这些点全部被其相邻粒子的圆所覆盖，则该粒子为内流体粒子，否则该粒子为将被识别为自由表面粒子；
    对于三维的情况，采用两步方法，首先，利用公式(4.1.7)检查相邻粒子数，检测出潜在的自由表面粒子，
    N p＜β 3dN p0  (4.1.7)
    β 3d是一个小于1的参数，根据情况取值，N p0表示圈定的范围内初始粒子密度；
    其次，通过建立以粒子为中心、半径为1.05r 0的球面来细化搜索，具体为：通过加权平均法确定一个指向该粒子周围粒子分布最稀疏方向的向量，通过均匀分布的点来离散球面上圆形区域，若所有这些点全部被其相邻粒子的球面所覆盖，则该粒子为内流体粒子，否则该粒子为将被识别为自由表面粒子。
  6. 根据权利要求5所述的方法，其特征在于：所述步骤S5的计算公式如式(5.1.1)所示，
    Figure PCTCN2019115597-appb-100026
  7. 根据权利要求6所述的方法，其特征在于：所述步骤S6的具体步骤如下：
    S6.1、结构运动描述模型建立
    结构运动描述模型采用固定整体坐标系XOY和附着体局部坐标系soη，其中，附着体局部坐标系的原点选择在结构运动描述模型的梁的重心处，因此，结构运动描述模型的梁上任意点的位置X可由式(6.1.1)描述，
    X＝X Rb+Rξ T  (6.1.1)
    其中，X Rb为结构运动描述模型的梁的质心位置，式(6.1.2)和式(6.1.3)中定义了结构运动描述模型的梁的局部坐标ξ和结构运动描述模型的梁的旋转矩阵R，
    ξ＝[s,η l+η 0]  (6.1.2)
    Figure PCTCN2019115597-appb-100027
    其中，θ为O-X轴与o-s轴逆时针的夹角，η l为未变形结构的坐标，对于结构运动描述模型的梁上的一个特定点是常数，η 0为中平面的挠度，任意平行于中平面的面的挠度均取中平面的相同值，根据模态叠加理论，将中平面的挠度η 0展开为式(6.1.4)，
    Figure PCTCN2019115597-appb-100028
    式(6.1.5)和(6.1.6)中定义了模态函数的列向量
    Figure PCTCN2019115597-appb-100029
    及对应的一般坐标q b，
    Figure PCTCN2019115597-appb-100030
    q b＝[q b1,q b2,q b3,...] T  (6.1.6)
    模态函数满足以下正交关系：
    Figure PCTCN2019115597-appb-100031
    I u是单位矩阵  (6.1.7)
    Figure PCTCN2019115597-appb-100032
    其中s 1和s 2为沿o-s轴积分的上下限，ρ l为结构运动描述模型的梁的线密度，E表示结构的弹性模量，I I为结构的惯性矩阵；
    S6.2、基于拉格朗日方程的固体响应控制方程建立
    基于拉格朗日方程建立模型：
    Figure PCTCN2019115597-appb-100033
    其中T和U分别表示结构运动描述模型的动能和势能，q j为任意刚柔模态对应的一般坐标，Q j为第j个坐标对应的非保守力；
    结构运动描述模型的梁为弹性的非均匀梁，弹性的非均匀梁的T、U、q j、Q j的具体形式如下：
    Figure PCTCN2019115597-appb-100034
    Figure PCTCN2019115597-appb-100035
    式中ρ s表示线密度，M表示结构的质量，X R和Y R分别对应结构质心的X和Y坐标，g表示重力加速度，η 1、η 2为沿轴积分η、矩阵积分U、列向量积分ψ 0和ψ 1的上下限，矩阵积分U和列向量积分ψ 0、ψ 1定义在式(6.2.4)和(6.2.6)中：
    Figure PCTCN2019115597-appb-100036
    Figure PCTCN2019115597-appb-100037
    Figure PCTCN2019115597-appb-100038
    将式(6.2.2)和(6.2.3)代入式(6.2.1)，可得在波浪冲击作用下弹性的非均匀梁的动力响应控制方程，即得到耦合刚体和弹性模态的固体响应控制方程，如式(6.2.7)所示：
    Figure PCTCN2019115597-appb-100039
    式(6.2.8)-(6.2.11)给出刚体和弹性模态的一般坐标对应的非保守力：
    Figure PCTCN2019115597-appb-100040
    Figure PCTCN2019115597-appb-100041
    Figure PCTCN2019115597-appb-100042
    Figure PCTCN2019115597-appb-100043
    式中，n s、n η分别为s和η的局部法向分量，用固定整体坐标系表示，X p和Y p是压强p作用点的整体X和Y坐标，e η是o-η轴的单位向量，积分是围绕结构运动描述模型的梁的周长进行的；
    S6.3、求解耦合刚体和弹性模态的固体响应控制方程。
  8. 根据权利要求6所述的方法，其特征在于：所述步骤S8的具体步骤如下：
    假设所有的流体和结构变量在时间步t＝t k-1，而后，在下一时间步t＝t k的交互详细过程如下：
    S8.1、假设结构在时间步t＝t k的加速度与上一时间步t＝t k-1的加速度是相同的，即
    Figure PCTCN2019115597-appb-100044
    则在时间步t＝t k的位置的D k和速度
    Figure PCTCN2019115597-appb-100045
    也可以通过有限差分法计算出来；
    通过式(8.1.1)(8.1.2)和式(8.1.3)～(8.1.6)计算出流体和结构交界面上每个点的相应的位置、速度和加速度，即
    Figure PCTCN2019115597-appb-100046
    和
    Figure PCTCN2019115597-appb-100047
    可以从新计算出来的结构运动和变形参数A计算得到，即式(2.1.1)，
    Figure PCTCN2019115597-appb-100048
    Figure PCTCN2019115597-appb-100049
    Figure PCTCN2019115597-appb-100050
    Figure PCTCN2019115597-appb-100051
    Figure PCTCN2019115597-appb-100052
    Figure PCTCN2019115597-appb-100053
    其中，X R为刚体的全局坐标，X CR为质心坐标，R R是将附着体局部坐标系soη与固定整体坐标系XOY联系起来的旋转矩阵；X fi表示中心线上的坐标；ξ i表示结构运动描述模型的梁在附着体局部坐标系soη下的中心线上的坐标，其表示方法为式(8.1.7)，式中的η i为结构运动描述模型的梁的挠度函数；
    ξ i＝[s i,η i]  (8.1.7)
    S8.2、利用交界面新的动力学信息作为新的边界条件，根据修正的MPS方法计算t＝t k时刻的流体运动，然后由此得到新的压力
    Figure PCTCN2019115597-appb-100054
    并将之应用于交界面在t＝t k时刻的第i次迭代；
    S8.3、比较步骤S2和步骤S7中的结构位置差值
    Figure PCTCN2019115597-appb-100055
    和
    Figure PCTCN2019115597-appb-100056
    若满足收敛条件，如式8.3.1)所示，则转到步骤S8.1进行下一时间步t＝t k+1运算，
    Figure PCTCN2019115597-appb-100057
    否则，利用式(8.3.2)修正第i+1次迭代的结构位置
    Figure PCTCN2019115597-appb-100058
    并且利用Newmark方法，更新速度
    Figure PCTCN2019115597-appb-100059
    和加速度
    Figure PCTCN2019115597-appb-100060
    根据
    Figure PCTCN2019115597-appb-100061
    和
    Figure PCTCN2019115597-appb-100062
    计算交界面的变量
    Figure PCTCN2019115597-appb-100063
    和
    Figure PCTCN2019115597-appb-100064
    使用这些修正的交界面变量，通过返回步骤S2进行第i+1次迭代，
    Figure PCTCN2019115597-appb-100065
    其中，χ i为Aitken弛豫因子，其值由式(8.3.3)计算：
    Figure PCTCN2019115597-appb-100066
    其中，
    Figure PCTCN2019115597-appb-100067
PCT/CN2019/115597 2019-03-22 2019-11-05 一种基于改进的移动粒子半隐式法和模态叠加方法求解强非线性时域水弹性问题设计方法 WO2020192126A1 (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US16/950,886 US20210086877A1 (en) 2019-03-22 2020-11-17 Method for analyzing hydroelastic effect of a marine structure

Applications Claiming Priority (4)

Application Number Priority Date Filing Date Title
CN201910223524.X 2019-03-22
CN201910223524 2019-03-22
CN201910413625.3 2019-05-17
CN201910413625.3A CN110750833A (zh) 2019-03-22 2019-05-17 一种基于改进的移动粒子半隐式法和模态叠加方法求解强非线性时域水弹性问题设计方法

Related Child Applications (1)

Application Number Title Priority Date Filing Date
US16/950,886 Continuation-In-Part US20210086877A1 (en) 2019-03-22 2020-11-17 Method for analyzing hydroelastic effect of a marine structure

Publications (1)

Publication Number Publication Date
WO2020192126A1 true WO2020192126A1 (zh) 2020-10-01

Family

ID=69275736

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2019/115597 WO2020192126A1 (zh) 2019-03-22 2019-11-05 一种基于改进的移动粒子半隐式法和模态叠加方法求解强非线性时域水弹性问题设计方法

Country Status (3)

Country Link
US (1) US20210086877A1 (zh)
CN (1) CN110750833A (zh)
WO (1) WO2020192126A1 (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113361173A (zh) * 2021-06-17 2021-09-07 中国空气动力研究与发展中心超高速空气动力研究所 一种完全基于当地流场参数对转捩模型可压缩修正的方法
CN113807034A (zh) * 2021-08-30 2021-12-17 西安交通大学 基于移动粒子半隐式法的轴对称流场二维模拟方法
CN116205152A (zh) * 2022-12-12 2023-06-02 中广核风电有限公司 一种海上漂浮式风机的数值模拟方法及装置
CN116911135A (zh) * 2023-07-24 2023-10-20 武汉理工大学 计及波浪下砰击载荷的非线性水弹性时域计算评估方法

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR102660215B1 (ko) * 2020-06-01 2024-04-23 스미토모 긴조쿠 고잔 가부시키가이샤 시뮬레이션 장치, 시뮬레이션 방법, 프로그램
US20220083702A1 (en) * 2020-09-11 2022-03-17 Autodesk, Inc. Techniques for designing structures using torsion-deformable spatial beam elements
CN112507600B (zh) * 2020-11-24 2024-03-29 西安交通大学 一种移动粒子半隐式法的对称边界条件的构建方法
CN113191066B (zh) * 2021-04-30 2022-12-09 西安交通大学 基于无网格法的核反应堆燃料元件失效分析方法
CN114638143B (zh) * 2022-03-21 2023-03-21 中国空气动力研究与发展中心计算空气动力研究所 一种适用于模拟水弹性问题的耦合数值计算方法
CN114912333B (zh) * 2022-07-19 2022-11-04 上海索辰信息科技股份有限公司 一种轻流体下和结构的声振耦合的模拟方法
CN115374730B (zh) * 2022-09-29 2023-04-14 中国海洋大学 一种风暴作用下液化海床土中管线稳定性的评估方法
CN115659522B (zh) * 2022-12-27 2023-03-28 中国空气动力研究与发展中心计算空气动力研究所 一种飞行器转捩位置预测方法、装置、设备及介质

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100161298A1 (en) * 2008-12-22 2010-06-24 Electronics And Telecommunications Research Institute Method for calculating force acting on interface between immiscible fluids in fluid simulation
CN102867094A (zh) * 2012-09-19 2013-01-09 西安交通大学 一种移动粒子半隐式算法中自由表面流动模型的构建方法
CN109374254A (zh) * 2018-11-21 2019-02-22 北京理工大学 一种航行体入水空泡特性的分析方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102012953B (zh) * 2010-11-04 2013-05-08 西北工业大学 Cfd/csd耦合求解非线性气动弹性仿真方法
CN103853884A (zh) * 2014-02-24 2014-06-11 昆明理工大学 一种水轮机活动导叶振动特性预测方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100161298A1 (en) * 2008-12-22 2010-06-24 Electronics And Telecommunications Research Institute Method for calculating force acting on interface between immiscible fluids in fluid simulation
CN102867094A (zh) * 2012-09-19 2013-01-09 西安交通大学 一种移动粒子半隐式算法中自由表面流动模型的构建方法
CN109374254A (zh) * 2018-11-21 2019-02-22 北京理工大学 一种航行体入水空泡特性的分析方法

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113361173A (zh) * 2021-06-17 2021-09-07 中国空气动力研究与发展中心超高速空气动力研究所 一种完全基于当地流场参数对转捩模型可压缩修正的方法
CN113807034A (zh) * 2021-08-30 2021-12-17 西安交通大学 基于移动粒子半隐式法的轴对称流场二维模拟方法
CN116205152A (zh) * 2022-12-12 2023-06-02 中广核风电有限公司 一种海上漂浮式风机的数值模拟方法及装置
CN116205152B (zh) * 2022-12-12 2024-06-07 中广核风电有限公司 一种海上漂浮式风机的数值模拟方法及装置
CN116911135A (zh) * 2023-07-24 2023-10-20 武汉理工大学 计及波浪下砰击载荷的非线性水弹性时域计算评估方法
CN116911135B (zh) * 2023-07-24 2024-02-13 武汉理工大学 计及波浪下砰击载荷的非线性水弹性时域计算评估方法

Also Published As

Publication number Publication date
CN110750833A (zh) 2020-02-04
US20210086877A1 (en) 2021-03-25

Similar Documents

Publication Publication Date Title
WO2020192126A1 (zh) 一种基于改进的移动粒子半隐式法和模态叠加方法求解强非线性时域水弹性问题设计方法
Sun et al. Numerical analysis of violent hydroelastic problems based on a mixed MPS-mode superposition method
Cummings et al. An integrated computational/experimental approach to UCAV stability & control estimation: overview of NATO RTO AVT-161
Huang et al. The simulation of deformation and vibration characteristics of a flexible hydrofoil based on static and transient FSI
Rao et al. Numerical study of the wave-induced slamming force on the elastic plate based on MPS-FEM coupled method
Chen et al. Numerical study of 3-D liquid sloshing in an elastic tank by MPS-FEM coupled method
Shobeiri The topology optimization design for cracked structures
Chen et al. A multi-resolution SPH-FEM method for fluid–structure interactions
Malik et al. Transient numerical simulations for hydrodynamic derivatives predictions of an axisymmetric submersible vehicle
Volpi et al. Composite bottom panel slamming of a fast planing hull via tightly coupled fluid-structure interaction simulations and sea trials
Cavallaro et al. Phenomenology of nonlinear aeroelastic responses of highly deformable joined wings
Zhang et al. Moving particle semi-implicit method coupled with finite element method for hydroelastic responses of floating structures in waves
Zhang et al. A two-phase flow model coupling with volume of fluid and immersed boundary methods for free surface and moving structure problems
Huang et al. Hydroelastic responses of an elastic cylinder impacting on the free surface by MPS-FEM coupled method
Wackers et al. Adaptive grid refinement for ship resistance computations
Wang et al. Extended variable-time-step Adams–Bashforth–Moulton method for strongly coupled fluid–structure interaction simulation
Lakshmynarayanana et al. Hydroelastic analysis of a flexible barge in regular waves using coupled CFD-FEM modelling
Gurugubelli et al. A variational projection scheme for nonmatching surface-to-line coupling between 3D flexible multibody system and incompressible turbulent flow
Shin Simulation of two-dimensional internal waves generated by a translating and pitching foil
Makarov et al. Parametric sensitivity of crane ship numerical model with experimental verification in a wave basin
Wan et al. Overset-rans computations of two surface ships moving in viscous fluids
Palacios et al. Shape sensitivity of free-surface interfaces using a level set methodology
Qu et al. Numerical simulation of sphere impacting water by SPH with hydrodynamics
Hung et al. Numerical investigation of dynamics of the flexible riser by applying absolute nodal coordinate formulation
Zhang et al. Prediction of roll hydrodynamic characteristics for a ship hull section with a high-order finite volume solver

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 19921743

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 19921743

Country of ref document: EP

Kind code of ref document: A1