A time-domain numerical method for calculating ship hydrodynamic performance

Abstract

This invention discloses a time-domain numerical method for calculating the hydrodynamic performance of ships in the field of marine hydrodynamics. The method includes: establishing a three-dimensional numerical model of the hull and free surface; decomposing the total velocity potential based on time-domain potential flow theory; constructing boundary integral equations; introducing a time-increasing monotonically increasing function to smooth the initial transient response; simultaneously iteratively solving the free surface evolution and the ship's six-degree-of-freedom motion equations; and achieving self-consistent coupling between flow field pressure and ship motion using a closed feedback loop. This method also supports nonlinear free surface boundary conditions, suppresses numerical dissipation through high-order extrapolation, and dynamically calculates the restoring force matrix based on surface element integrals, significantly improving the calculation accuracy and stability of ship motion response in confined waters. It is suitable for engineering prediction of complex conditions such as near-shore navigation, shallow water resonance, and shoreline reflection.

Description

Technical Field

[0001] This invention relates to the field of marine hydrodynamics, and in particular to a time-domain numerical method for calculating the hydrodynamic performance of ships. Background Technology

[0002] When ships navigate in such waters, their hydrodynamic performance differs significantly from that in deep water or unlimited water, exhibiting shallow water effects and shoreline effects. Therefore, the risk of accidents is greater for ships navigating in restricted waters than in unlimited waters. Although model testing provides a direct and accurate method for studying the hydrodynamic performance of ships in restricted waters, its high cost limits comprehensive studies under various operating conditions, such as different water depths, ship-to-shore distances, and speeds.

[0003] However, with the rapid development of computer technology and numerical simulation methods, the research scope of ship hydrodynamics and performance analysis is expanding, shifting from traditional linear frequency domain analysis to more complex nonlinear time domain simulation, and from single performance evaluation to comprehensive consideration of navigation performance. Computational fluid dynamics (CFD) is increasingly used in the study of ship-wave interactions, but its computational intensity in dealing with wave-related problems limits its application in rapid engineering design prediction.

[0004] Currently, in the field of ship and ocean engineering design, mainstream commercial software based on potential flow theory still dominates. Although it has certain advantages in hydrodynamic performance analysis, most of them are still based on linear frequency domain theory, which has limited adaptability to nonlinear problems. At the same time, most of them adopt unidirectional transmission or rigid coupling, that is, the flow field is calculated first and then the force is given, or a fixed number of iterations is used, which cannot guarantee the self-consistency of each step, resulting in numerical divergence or decreased accuracy. Therefore, we propose a time-domain numerical method for calculating the hydrodynamic performance of ships. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a time-domain numerical method for calculating the hydrodynamic performance of ships. This method employs a fluid-structure interaction solution strategy, ensures the self-consistency of each time step through synchronous iteration and convergence determination, avoids numerical divergence problems, and supports nonlinear free surface boundary conditions, enabling it to handle nonlinear phenomena such as large-amplitude motion.

[0006] The objective of this invention is achieved as follows: a time-domain numerical method for calculating the hydrodynamic performance of ships, comprising the following steps:

[0007] S1. Establish a three-dimensional numerical model of the hull and the free liquid surface. The model is constructed based on the hull geometry, fluid domain boundary conditions and environmental parameters, and characterizes the physical response boundary of the ship under wave action.

[0008] S2, based on the time-domain potential flow theory, decomposes the total velocity potential of the flow field into multiple physical components, solves the components related to ship motion and wave scattering, and pre-determines the remaining components through physical laws.

[0009] S3. Based on the dynamic variables determined in S2, construct the boundary integral equations and discretize them to form a simultaneous algebraic matrix system.

[0010] S4, during the time-domain propagation process, a buffer mechanism is introduced to smooth the initial transient response, and a high-order time integration algorithm is used to synchronously iteratively solve the free surface evolution and the six-degree-of-freedom motion equations of the ship.

[0011] S5, based on the flow field pressure distribution after the convergence of S4 iteration, performs integral calculations on the wetted surface of the hull to calculate the total hydrodynamic force and total torque on the ship, and feeds this load back into the ship's six-degree-of-freedom motion equations to update its motion state, driving the dynamic solution of the next moment until the full-time domain simulation is completed.

[0012] Optionally, in step S2, the decomposition of the total velocity potential of the flow field includes: background inflow, external wave incidence, wave generation from ship movement, and unsteady scattering potential.

[0013] The background flow and external wave incidence are known input conditions.

[0014] The wave-generating potential of a moving ship and the unsteady scattering potential are solved in real time as dynamic variables.

[0015] Optionally, the boundary integral equation includes the wetted surface boundary of the hull, the free liquid surface boundary, and the confined water area boundary;

[0016] Each boundary is meshed using high-order discretization.

[0017] Optionally, the restricted water boundary includes the seabed and the shoreline. The boundary is processed by introducing a mirror source effect in the integral kernel function to simulate the reflection and constraint effect of the boundary on the flow field.

[0018] Optionally, the step of performing mesh partitioning on each boundary through high-order discretization specifically involves:

[0019] The wetted surface of the hull, the free liquid surface, and the shore wall are divided into multiple small regions. Each small region is used as a calculation unit. The average velocity potential and normal rate of change of its surface jointly characterize the local fluid behavior, thereby transforming the continuous integral equation into a set of algebraic matrix systems that can be solved by a computer.

[0020] The boundary treatment of the free liquid surface adopts linear or nonlinear kinetic energy conservation kinematic boundary conditions, and combines high-order extrapolation techniques of normal velocity potential to suppress numerical dissipation.

[0021] Optionally, in step S4, the buffering mechanism specifically involves multiplying the scattering potential generated by the ship's motion by a monotonically increasing control function during the initial calculation period until the control function stabilizes, at which point the scattering potential participates in subsequent calculations with its full amplitude.

[0022] Optionally, in step S4, the synchronous iterative solution specifically involves:

[0023] Within each time step, the free surface height, velocity potential and its normal derivative, as well as the ship's displacement, velocity and acceleration are solved sequentially through a multi-stage calculation process. The intermediate calculation results of each stage are then linearly combined according to a preset convergence weighting coefficient to generate the updated values ​​for the next time step.

[0024] Optionally, the six-degree-of-freedom equations of motion for the ship are constructed based on the Newton-Euler equations;

[0025] The six-degree-of-freedom equations of motion for the ship include a mass matrix and a restoring force matrix;

[0026] The mass matrix and restoring force matrix are calculated based on the total mass of the hull, the coordinates of the center of gravity, the inertia tensor, and the still water buoyancy distribution.

[0027] Optionally, the synchronous iterative solution in step S4 and the feedback update in step S5 constitute a closed feedback loop;

[0028] Specifically, within each time step, the six-degree-of-freedom motion state of the ship directly affects the setting of the flow field boundary conditions, and the hydrodynamic loads calculated from the flow field pressure distribution directly drive the update of the ship's motion equations.

[0029] The two achieve state self-consistency through repeated coupled calculations within the time step until the entire simulation time history is completed.

[0030] Optionally, within each time step, the following iterative process is performed:

[0031] 1) Based on the ship's current six-degree-of-freedom position and attitude, update the hull boundary conditions and solve the hydrodynamic equations to obtain the flow field pressure distribution;

[0032] 2) Based on the pressure distribution, calculate the total hydrodynamic force and total torque acting on the ship by integrating over the hull surface;

[0033] 3) Using the total hydrodynamic force and total torque as external force inputs, solve the six-degree-of-freedom motion equations of the ship to obtain the ship's motion state at the next iteration time;

[0034] 4) Repeat steps 1) to 3) until the change in the ship's motion state is less than the preset convergence threshold, and determine that the coupled system has reached self-consistency within this time step;

[0035] 5) Using the self-consistent state as the initial condition for the next time step, a higher-order time integration algorithm is used to solve the problem.

[0036] Compared with the prior art, the beneficial effects of the present invention are as follows: First, by decomposing the flow field velocity potential into and solving only the dynamic variables in real time, computational redundancy is significantly reduced and solution efficiency is improved; Second, by introducing a buffer mechanism based on monotonically increasing cosine, the scattering potential smoothly transitions from zero to the full amplitude, effectively suppressing numerical oscillations caused by sudden changes in boundary conditions at the initial moment and reducing the risk of divergence.

[0037] By employing joint processing of the free surface boundary, this invention suppresses numerical dissipation in large waves using traditional linear methods and improves waveform preservation accuracy. Within the time domain framework, this invention achieves a two-way, self-consistent, iteratively coupled closed loop between the ship's six-degree-of-freedom motion and the flow and pressure fields. By pre-setting dynamic weighting coefficients to coordinate the evolution of the free surface, the solution of the potential function, and the iteration of the motion equations, the system achieves stable convergence within each time step. Furthermore, it ensures that the restoring force matrix is ​​dynamically updated based on the Gaussian integral of the surface element method, truly reflecting the buoyancy center shift caused by changes in ship attitude, thus overcoming the fatal flaw of "static fitting" of the still water restoring force in traditional methods.

[0038] By using the mirror source Green's function and the singularity analytical stripping method to jointly process the seabed and shoreline boundary, we can achieve the simulation of flow field reflection under complex boundary conditions without additional domain, thereby reducing the prediction error of near-shore navigation response. Attached Figure Description

[0039] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort.

[0040] Figure 1 This is a schematic diagram of the water surface grid provided by the present invention.

[0041] Figure 2 This is a schematic diagram of the ship grid provided by the present invention.

[0042] Figure 3 This is a schematic diagram of the time-domain numerical method for calculating the hydrodynamic performance of ships provided by the present invention.

[0043] Figure 4 This is a schematic diagram comparing the heave motion results of the ship provided by this invention with those obtained from experiments and other numerical calculations.

[0044] Figure 5This is a schematic diagram comparing the results of the pitching motion of a sailing vessel provided by the present invention with those obtained from experiments and other numerical calculations.

[0045] Figure 6 This is a schematic diagram comparing the heave motion results of near-shore navigation vessels provided by the present invention with those obtained from other numerical calculations.

[0046] Figure 7 This is a schematic diagram comparing the pitching motion results of near-shore navigation vessels provided by the present invention with those obtained from other numerical calculations. Detailed Implementation

[0047] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0048] like Figures 1 to 7 The time-domain numerical method shown includes the following steps for calculating the hydrodynamic performance of a ship:

[0049] S1. Establish a three-dimensional numerical model of the hull and the free liquid surface. The model is constructed based on the hull geometry, fluid domain boundary conditions and environmental parameters, and characterizes the physical response boundary of the ship under wave action.

[0050] In step S1, the process of establishing the three-dimensional numerical model involves three main steps: defining the coordinate system, determining the computational domain, and generating the mesh.

[0051] It should be noted that this invention provides a time-domain coupled solution method based on physical modeling and numerical stability control for the field of ship engineering design, which solves the technical problem that existing commercial software cannot calculate the large motion response of ships in restricted waters. The method provided by this invention can be implemented on a general-purpose computer platform in any programming language, and its technical solution does not depend on a specific operating system, programming language or hardware structure.

[0052] This invention employs two coordinate systems to describe the motion of a ship in waves: a fixed coordinate system and a body-following coordinate system. The origin of the fixed coordinate system is set on the still water surface, with the X-axis pointing in the direction of wave propagation, the Y-axis perpendicular to the X-axis, and the Z-axis pointing upwards as positive. The origin of the body-following coordinate system is set at the ship's center of gravity, moving with the ship, and is used to describe the ship's six degrees of freedom motion. During numerical calculations, a coordinate transformation matrix is ​​needed to convert between the two coordinate systems. Physical quantities in the fixed coordinate system are mapped to the body-following coordinate system through this transformation matrix, and vice versa.

[0053] The determination of the computational domain needs to be based on the actual engineering requirements and the required computational accuracy to determine the size of the virtual wave field. For calculations involving infinite water areas, the computational domain is typically set with a truncated boundary far from the hull, and appropriate wave-damping measures are used at the truncated boundary to simulate the outward propagation of waves. For calculations involving confined water areas, the computational domain also needs to include the seabed and shoreline boundaries. The seabed boundary is processed using mirror source technology, and the shoreline boundary is simulated by setting a surface condition with zero normal velocity to simulate the constraint effect of a solid wall. The horizontal range of the computational domain is typically 5 to 10 times the ship's length, and the vertical range extends from the seabed to 2 to 3 times the wave height above the still water surface.

[0054] Mesh generation is a core step in numerical model building, directly impacting computational accuracy and efficiency. The wetted surface of the hull is meshed using quadrilateral high-order elements, with the velocity potential on each element approximated using bilinear or biquadratic interpolation functions to achieve computational accuracy of second order or higher. Mesh generation for the free surface needs to consider the influence of wavelength; sufficient mesh nodes should be ensured within the wavelength range to capture wave phase information, typically no less than 20 nodes per wavelength.

[0055] Furthermore, the mesh near the hull is locally refined, while areas far from the hull can use a sparser mesh to save computational resources. Figure 1 and Figure 2 The figures represent typical forms of free surface grids and hull grids, respectively. As can be seen from the figures, the grid distribution gradually becomes sparser from near the hull to the distance.

[0056] S2, based on the time-domain potential flow theory, decomposes the total velocity potential of the flow field into multiple physical components, solves the components related to ship motion and wave scattering, and pre-determines the remaining components through physical laws.

[0057] In step S2, the decomposition of the total velocity potential of the flow field includes: background inflow, external wave incidence, wave generation from ship movement, and unsteady scattering potential.

[0058] Among them, the background inflow and the external wave incidence are known input conditions;

[0059] The wave-generating potential of a moving ship and the unsteady scattering potential are solved in real time as dynamic variables.

[0060] Background flow refers to the uniform flow field generated when a ship sails at a constant speed in still water. This component is determined by the ship's speed and is used as a known input condition in the calculation.

[0061] External wave incidence refers to the wave state of an incoming wave before it is disturbed by a ship; its velocity potential and wave height both have analytical expressions.

[0062] Furthermore, in step S2, based on the time-domain potential flow theory, the total velocity potential of the flow field is decomposed into multiple physical components. This decomposition method can simplify complex flow field problems into a superposition of several relatively simple subproblems, making them easier to solve and process separately.

[0063] The decomposition formula for the total velocity potential of the flow field is as follows: The total velocity potential function can be expressed as the velocity potential of the uniform inflow (- The sum of the absolute velocity potential of the flow field relative to the reference coordinate system and ( Its mathematical expression is:

[0064]

[0065] In the formula, It represents the total velocity potential of the entire flow field relative to a fixed coordinate system when a ship is moving in waves;

[0066] U is the ship's forward speed;

[0067] x is the position coordinate along the length of the ship;

[0068] It is the absolute velocity potential;

[0069] Furthermore, the total velocity potential of the flow field It can be divided into three parts: the wave-making potential of ship movement. Wave incidence potential and unsteady scattering potential ,

[0070] The formula is as follows:

[0071]

[0072] The incident wave velocity potential and incident wave height are usually expressed analytically, as shown in the formulas below. Therefore, in actual calculations, only the solution needs to be found. and These correspond to the wave-making problem of ship movement and the time-domain motion problem, respectively;

[0073] The expression for the incident wave velocity potential is:

[0074]

[0075] Incident wave height The expression is:

[0076]

[0077]

[0078] Where A is the wave amplitude, g is the gravitational acceleration, ω is the incident wave frequency, k is the wave number, and β is the incident wave angle. Let d be the encounter frequency and d be the water depth.

[0079] Ship movement creates wave potential This refers to the steady wave-making problem generated by a ship traveling at a constant speed in still water. The potential function satisfies the Laplace equation, as well as the linearized boundary conditions of the free surface and the impenetrable boundary conditions of the object surface. The moving wave-making potential, as a steady-state component, is solved once at the beginning of the calculation and remains unchanged in subsequent calculations.

[0080] Unsteady scattering potential It is the time-varying wave component generated when a ship moves in waves. This component is affected by both the incident wave and the ship's motion, and is a dynamic variable that needs to be solved in real time in time-domain calculations.

[0081] Corresponding to the velocity potential, the free liquid surface Elevation can also be decomposed into moving waves. Incident wave unsteady scattering The three parts are decomposed into the following formula:

[0082]

[0083] Furthermore, based on the velocity potential decomposition form determined in step S2, corresponding boundary conditions need to be established for each velocity potential component.

[0084] The boundary conditions for the moving wave-making velocity potential include the following categories: Within the fluid domain, the velocity potential satisfies the Laplace equation:

[0085]

[0086] In the formula It is the unit normal vector pointing inward from the wetted surface.

[0087] The boundary conditions of the unsteady scattering potential are solved by a joint boundary value problem around radiation, and the boundary conditions are stated as follows:

[0088]

[0089] In the formula n j Let the j-th component of the normal vector be defined as:

[0090]

[0091] m j The j-th component of the so-called m terms can be represented as:

[0092]

[0093] Assuming the steady wave is very small, Neumann-Kelvin linearization can be used to simplify the m-terms.

[0094]

[0095] For the initial boundary value problem when a ship is navigating near a shore wall, it is only necessary to add a surface condition on the shore wall, namely:

[0096]

[0097] S3. Based on the dynamic variables determined in S2, construct the boundary integral equations and discretize them to form a simultaneous algebraic matrix system.

[0098] The boundary integral equations include the wetted surface boundary of the hull, the free surface boundary, and the confined water boundary;

[0099] The boundaries are meshed using high-order discretization.

[0100] Furthermore, the Rankine source and its image about the seabed are represented by the Green's function, as follows:

[0101]

[0102] In the formula: x is the field point, and its coordinates are (x0, y0, z0);

[0103] x0 is the source point, and its coordinates are (x, y, z);

[0104] R and R1 represent the distances from the source point to the field point and its mirror image, respectively:

[0105] in:

[0106]

[0107] The boundary integral equation is established based on Green's second identity, transforming the Laplace equation into a boundary integral equation concerning the velocity potential on the boundary of the computational domain.

[0108] For any point P in the flow field, we have:

[0109]

[0110] In the formula: α is the fixed angle coefficient, the value of which depends on the position of the field point P relative to the boundary S, and can be obtained by direct calculation;

[0111] The discrete form of the boundary integral equation can be written as:

[0112]

[0113] S is the boundary of the entire computational domain, including the wetted surface of the hull.B Free liquid surface S F 、Seabed S D and the shoreline boundary S W .

[0114] Furthermore, in the numerical implementation of the boundary integral equation, when the distance between the computation field point and the source point is extremely close, the Green's function will exhibit strong numerical singular behavior due to the geometric distance approaching zero. If the standard numerical integration method is used directly, it will lead to severe distortion or even divergence of the computation matrix.

[0115] To address this issue, this invention employs a partitioning strategy, subdividing the computational domain into near-field and far-field regions. In the near-field region, where the field point and source point are adjacent or extremely close, the system avoids conventional numerical integration methods. Instead, it utilizes an analytical correction technique based on the separation of geometric characteristics and physical properties. This method first extracts the principal singularity component that can be precisely integrated from the integrator kernel and, based on the fixed angle coefficients corresponding to this component, directly assigns it an analytical closed-form solution, thereby completely avoiding numerical errors.

[0116] The remaining non-singular components are relatively smooth functions with gradual changes, allowing for numerical processing using conventional high-order interpolation integration methods without affecting overall accuracy. This separation mechanism ensures accurate handling of singular terms while preserving the numerical efficiency of non-singular terms, enabling the entire boundary integral system to converge stably under arbitrary grid densities and ship geometries, free from matrix ill-conditioning or computational oscillations.

[0117] Furthermore, this processing method does not depend on a specific mesh shape or size, and is applicable to any complex hull shape and free surface boundary, exhibiting good versatility and engineering robustness.

[0118] Specifically, the confined water boundary consists of two parts: the seabed and the shoreline. The seabed boundary is set according to the actual water depth, while the shoreline boundary is determined according to the geometry of the actual shoreline structure. In numerical calculations, these boundaries are simulated by introducing a mirror source effect into the integral kernel function to model the reflection and constraint effects of the boundaries on the flow field.

[0119] The mirror source method used in this invention is only used to process near-shore reflection simulation under conditions of low-amplitude incident waves (wave amplitude less than 5% of ship width) and small disturbance motion (heave / roll less than 5°). When the ship motion exceeds the above threshold, the system automatically switches to the PML (Perfect Matching Layer) boundary processing method based on physical boundary layer absorption to ensure energy conservation and numerical stability.

[0120] Specifically, for seabed boundaries, a mirror image of the source point is placed symmetrically below the actual seabed; for shoreline boundaries, a mirror image of the shoreline is placed symmetrically on the other side of the shoreline. This approach can accurately simulate the impact of finite water areas on wave propagation without increasing the computational domain.

[0121] The surface condition of the shore wall boundary is set to zero normal velocity, indicating that the shore wall is an impenetrable solid boundary and waves will be completely reflected when they encounter it. The combination of the mirror source method and the surface boundary condition can accurately simulate the reflection and constraint effect of the shore wall on the flow field.

[0122] The seabed boundary is also handled using the mirror source method. This method eliminates the need to explicitly set seabed boundary conditions in the computational domain, as the Green function itself automatically satisfies the solid wall conditions of the seabed.

[0123] Specifically, the boundaries are meshed using high-order discretization, as follows:

[0124] The wetted surface of the hull, the free liquid surface, and the shore wall are divided into multiple small regions. Each small region is used as a calculation unit. The average velocity potential and normal rate of change of its surface jointly characterize the local fluid behavior, thereby transforming the continuous integral equation into a set of algebraic matrix systems that can be solved by a computer.

[0125] The boundary treatment of the free liquid surface adopts linear or nonlinear kinetic energy conservation kinematic boundary conditions, combined with normal velocity potential to suppress numerical dissipation.

[0126] The high-order discretization of the boundaries is one of the key features of this invention. During the mesh generation process, the wetted surface of the hull, the free liquid surface, and the shoreline boundary are divided into multiple small regions, each serving as an independent computational unit. The velocity potential on each unit is approximated using a high-order interpolation function, specifically biquadratic interpolation or a higher-order interpolation method.

[0127] Furthermore, quadrilateral nine-node elements are used for the wetted surface of the hull, while quadrilateral or triangular higher-order elements are used for the free surface. The velocity potential of each element is characterized by the average velocity potential and normal rate of change on its surface. This approach transforms the continuous integral equations into a set of algebraic matrix systems that can be solved by a computer.

[0128] The handling of free surface boundary conditions is a crucial step in time-domain computation. This invention supports two forms of kinematic boundary conditions: linear kinetic energy conservation kinematic boundary conditions and nonlinear kinetic energy conservation kinematic boundary conditions.

[0129] Nonlinear boundary conditions retain more nonlinear terms, enabling more accurate simulation of large-amplitude wave motion. To suppress numerical dissipation during numerical calculation, this invention also employs a high-order extrapolation technique for the normal velocity potential. Using the normal velocity potential data from the first three time steps, a quadratic interpolation polynomial for the time series is constructed to extrapolate the normal derivative at the next time step.

[0130] This method reduces numerical dissipation error by two orders of magnitude without incurring additional computational overhead, and is applicable to complex free surface topologies. Extrapolation is applied only to free surface mesh nodes, and linear interpolation is used at boundary edge nodes to avoid extrapolation divergence.

[0131] Specifically, the six-degree-of-freedom equations of motion for ships are constructed based on the Newton-Euler equations;

[0132] The six-degree-of-freedom equations of motion for a ship include the mass matrix and the restoring force matrix;

[0133] The mass matrix and restoring force matrix are calculated based on the total mass of the hull, the coordinates of the center of gravity, the inertia tensor, and the still water buoyancy distribution.

[0134] Furthermore, as the ship moves forward in the waves, it is subjected to wave excitation forces and torques, resulting in displacements of six degrees of freedom. According to Newton's second law, the equations of motion are as follows:

[0135]

[0136] In the formula: M is the rigid body mass matrix, C is the restoring force matrix, and F is the total pressure and torque components. The general form of the M matrix is ​​given below:

[0137] In the formula: M is the total mass of the ship, including the total mass of the hull structure, equipment, load and fuel;

[0138] (x c y c , z c ) represents the three-dimensional spatial coordinates of the ship's center of gravity relative to the origin of the fixed coordinate system. When the ship undergoes six-degree-of-freedom motion, this coordinate will change with the attitude. However, in the mass matrix, the position of the center of gravity in the initial still water state is usually used as a reference and a fixed input parameter.

[0139] Calculated by the following formula:

[0140]

[0141] The restoring force matrix C describes the hydrostatic restoring forces and moments of an object during heaving, pitching, and rolling motions. The general form of the restoring force matrix C is as follows:

[0142] In the formula: Indicates the density of seawater;

[0143] A. The waterline area of ​​a ship in still water, that is, the horizontal cross-sectional area where the hull intersects the still water surface.

[0144] V represents the displacement volume of the ship in still water, which is determined by the total volume enclosed by the underwater parts of the hull.

[0145] , The centroid of the waterline area is represented by its horizontal and vertical coordinates in the ship's coordinate system, which is the offset of the center of symmetry of the waterline area relative to the ship's center of gravity.

[0146] , It represents the second moment of the waterline area about its own centroid, and is used to calculate the restoring moment of a ship due to changes in buoyancy during roll and pitch.

[0147] The definitions of each physical quantity are as follows:

[0148]

[0149] To accurately calculate the still water restoring force and restoring moment of a ship under wave action, this invention uses numerical discretization to solve the key physical quantities in the restoring force matrix based on the underwater shape changes of the ship under different six-degree-of-freedom attitudes. Each element of this matrix originates from the geometric distribution characteristics of the interface between the submerged part of the hull and the free surface, including parameters such as displacement volume, center of gravity shift, and inertial distribution.

[0150] In practice, the hull surface is divided into multiple discrete small surface elements, each with a definite position and normal in space. By traversing all these surface elements, the system evaluates their buoyancy contribution under the hull's current attitude. For physical quantities involving area integration, such as the moment of inertia and the projection of the center of gravity at the water surface, the system performs element-by-element accumulation calculations based on the projected area of ​​each surface element on the liquid surface and its spatial position. For terms involving volume integration, such as the inertia tensor and the displacement volume, the system performs internal accumulation using the surface element method based on the underwater volume enclosed by the hull's closed surface, ensuring that the buoyancy center and inertial distribution dynamically update in real time to follow the attitude when the hull rolls, pitches, or heaves.

[0151] The calculation process is entirely automated, relying on the relationship between the hull mesh data and attitude transformation, without requiring any analytical assumptions or simplified models. Even under complex conditions such as large movements or shallow water near the shore, the system can maintain the continuity and physical consistency of the restoring force matrix, thereby ensuring the realism and stability of the subsequent six-degree-of-freedom motion equations.

[0152] Specifically, in S4, a buffer mechanism is introduced to smooth the initial transient response during the time-domain propagation process, and a high-order time integration algorithm is used to synchronously and iteratively solve the free surface evolution and the six-degree-of-freedom motion equations of the ship.

[0153] S5, based on the flow field pressure distribution after the convergence of S4 iteration, performs integral calculations on the wetted surface of the hull to calculate the total hydrodynamic force and total torque on the ship, and feeds this load back into the ship's six-degree-of-freedom motion equations to update its motion state, driving the dynamic solution of the next moment until the full-time domain simulation is completed.

[0154] Specifically, in step S4, the buffering mechanism is as follows: during the initial time period of the calculation, the scattering potential generated by the ship's motion is weighted by multiplying it by a control function that increases monotonically with time until the control function tends to stabilize, and the scattering potential participates in the subsequent calculation with its full amplitude.

[0155] Furthermore, the buffering mechanism is introduced to avoid numerical oscillations or divergence caused by the sudden appearance of the scattering potential at the initial moment in numerical calculations.

[0156] The buffer function takes the form of a cosine function that monotonically increases with time, and its mathematical expression is:

[0157]

[0158] Where T m The buffer time is an integer multiple of the wave period, typically 3 to 5 wave periods. When the time is less than the buffer time, the buffer function monotonically decreases from 1 to 0.

[0159] Once the buffer time is reached or exceeded, the value of the buffer function remains constant at 1, and the scattering potential participates in subsequent calculations at its full amplitude.

[0160] The buffering mechanism handles the scattering potential as follows: During the initial calculation period, the scattering potential generated by the ship's motion needs to be weighted, and the weighted scattering potential participates in subsequent hydrodynamic calculations. Specifically, the scattering potential participating in the calculation at any time t is:

[0161]

[0162] in Let W(t) be the theoretical value of the scattering potential, and W(t) be the buffer function. This represents the effective scattering potential after buffering.

[0163] When the buffer function tends to stabilize (i.e. When the scattering potential is at its full amplitude, it is used in subsequent calculations.

[0164] The advantages of using a buffer mechanism are: it allows the scattering potential to gradually increase from zero, avoiding numerical instability caused by abrupt changes at the initial moment; at the same time, the cosine function form ensures the continuity of the buffer function and its first derivative, further improving the stability of numerical calculation.

[0165] Specifically, in step S4, the synchronous iterative solution is as follows:

[0166] Within each time step, the free surface height, velocity potential and its normal derivative, as well as the ship's displacement, velocity and acceleration are solved sequentially through a multi-stage calculation process. The intermediate calculation results of each stage are then linearly combined according to a preset convergence weighting coefficient to generate the updated values ​​for the next time step.

[0167] Furthermore, synchronous iterative solving enables the simultaneous solution of the free surface evolution and the ship's six-degree-of-freedom motion equations, avoiding the insufficient coupling problem caused by the independent calculation of flow field and ship motion in traditional methods.

[0168] The multi-stage computation process is the basic framework for synchronous iterative solutions. Within each time step Δt, the computation process is divided into the following three stages performed sequentially:

[0169] The first stage is the free surface height solution stage. Based on the flow field velocity potential distribution at the current time t, the free surface boundary condition equations are solved to obtain the distribution of the free surface height ζ. The free surface boundary conditions are solved using kinetic energy conservation kinematic boundary conditions, and their expression is:

[0170] The time integration method for the boundary conditions of the object surface and the iterative solution steps for the equations of motion both adopt the fourth-order Runge-Kutta numerical integration method.

[0171] For the time integration process, the formula is as follows:

[0172]

[0173]

[0174] The second stage involves solving for the velocity potential and its normal derivative. Based on the updated free surface boundary conditions, the boundary integral equations are solved to obtain the velocity potential and its normal derivative across the entire flow field boundary.

[0175] The third stage is the ship motion solution stage. Based on the obtained velocity potential distribution, the pressure distribution on the wetted surface of the hull is calculated using Bernoulli's equation. Then, the surface integral of the pressure is used to obtain the hydrodynamic forces and torques acting on the ship. Finally, these are substituted into the ship's six-degree-of-freedom motion equations to obtain the ship's displacement, velocity, and acceleration.

[0176] Determining the convergence weighting coefficients is a key parameter in the multi-stage computation process. This invention uses preset convergence weighting coefficients to linearly combine the intermediate computation results of each stage to generate the updated value for the next time step. The values ​​of the convergence weighting coefficients are determined based on the convergence characteristics and numerical stability requirements of the computation results at each stage.

[0177] Furthermore, the preset convergence weighting coefficients are based on the multi-objective optimization results of balancing system energy transfer efficiency and numerical dissipation. In this invention, after verification by a large number of numerical experiments, a stable convergence parameter range is finally obtained.

[0178] In the implementation of this invention, the first stage (solving the height of the free liquid surface) has a fast convergence speed but low accuracy because it is less affected by the velocity potential.

[0179] The second stage (solving the velocity potential) is the core of the algorithm; it converges slowly but determines the physical reality.

[0180] The third stage (solution of ship motion) relies on the results of the first two stages and has a lag effect; based on numerous numerical tests (such as the Wigley III ship), , (0.5~1.5, water depth d / L=0.1~0.5), it was found that when the weighting coefficients were (0.2, 0.6, 0.2), the overall convergence speed of the system increased by 37%, and there was no oscillation divergence phenomenon.

[0181] This combination ensures the controllability of free surface renewal, the dominance of potential function solution, and the stability of motion response. In a preferred embodiment of the invention, the weighting coefficients can be set to the range (0.1~0.3, 0.5~0.7, 0.1~0.2), more preferably (0.2, 0.6, 0.2).

[0182] The time integration method employs the fourth-order Runge-Kutta 4-RK method, which boasts high computational accuracy and numerical stability. For the time progression of the free surface and velocity potential, the recursive formula is as follows:

[0183]

[0184]

[0185] The formulas for calculating the derivatives of each order are as follows:

[0186]

[0187]

[0188]

[0189]

[0190]

[0191]

[0192]

[0193]

[0194] The normal derivative of the scattering potential on the object's surface is obtained by calculating the incident potential. The scattering potential on the free water surface is then calculated using a time-step process. This allows us to derive the velocity potential integral equation for the new moment, thus obtaining the fluid motion at that moment. Subsequently, the next new moment can be calculated, and this process is repeated until the calculation period ends.

[0195] Determining the velocity of an object requires utilizing the boundary conditions of its surface, while determining its displacement requires calculating wave forces. These physical quantities can be obtained from the object's equations of motion. In numerical simulations, the equations of motion are also solved iteratively using the standard 4-RK method, specifically:

[0196]

[0197] The displacement and velocity of an object can be written in the following form:

[0198]

[0199]

[0200] In the formula, M1, M2, M3, and M4 are respectively:

[0201]

[0202]

[0203]

[0204]

[0205] For fluid-structure interaction problems, iterative solutions are typically required. In the traditional 4-RK algorithm, each iteration requires four calculations. First, the flow field equations are solved using the displacement and velocity at the current time t, yielding the result at time t, which is M1. Then, M1 is used to calculate the object's displacement and velocity for the next iteration, and the flow field equations are solved again to obtain M2. This process is repeated to obtain M3 and M4.

[0206] Specifically, the synchronous iterative solution in step S4 and the feedback update in step S5 constitute a closed feedback loop. Specifically, in each time step, the six-degree-of-freedom motion state of the ship directly affects the setting of the flow field boundary conditions, and the hydrodynamic load calculated by the flow field pressure distribution directly drives the update of the ship's motion equations.

[0207] The two achieve state self-consistency through repeated coupled calculations within the time step until the entire simulation time history is completed;

[0208] In this invention, the maximum number of coupling iterations per time step is set to 8. If convergence is not achieved within 8 iterations, the Δt of that time step is automatically reduced to 0.5 times the original value and recalculated. If convergence exceeds the threshold for 3 consecutive iterations, the Δt is gradually extended. Under the premise of ensuring accuracy, the average time consumption is controlled within 2 to 3 times that of the traditional one-way transmission method.

[0209] Specifically, within each time step, the following iterative process is executed:

[0210] 1) Based on the ship's current six-degree-of-freedom position and attitude, update the hull boundary conditions and solve the hydrodynamic equations to obtain the flow field pressure distribution;

[0211] 2) Based on the pressure distribution, calculate the total hydrodynamic force and total torque acting on the ship by integrating over the hull surface;

[0212] 3) Using the total hydrodynamic force and total torque as external force inputs, solve the six-degree-of-freedom motion equations of the ship to obtain the ship's motion state at the next iteration time;

[0213] 4) Repeat steps 1) to 3) until the change in the ship's motion state is less than a preset convergence threshold, which is set to a displacement of less than or equal to 10. -4 m, speed less than or equal to 10 -4 m / s, indicating that the coupled system is self-consistent within this time step;

[0214] 5) Using the self-consistent state as the initial condition for the next time step, a higher-order time integration algorithm is used to solve the problem.

[0215] The synchronous iterative solution in step S4 and the feedback update in step S5 constitute a closed feedback loop, which is the key mechanism for achieving true fluid-structure interaction in this invention.

[0216] Furthermore, the implementation process of the closed feedback loop is as follows: In each time step, the six-degree-of-freedom motion state of the ship directly affects the setting of the flow field boundary conditions;

[0217] Specifically, when a ship undergoes displacement and attitude changes, the position of its wetted hull surface is updated accordingly, and the corresponding boundary conditions need to be adjusted based on the new ship position. Simultaneously, the hydrodynamic loads calculated from the flow field pressure distribution directly drive the update of the ship's motion equations, and the new ship motion state in turn affects the flow field boundary conditions at the next moment. This two-way data exchange process repeats at each time step until a self-consistent state is achieved.

[0218] The convergence criterion is crucial for ensuring the self-consistency of each time step. This invention uses the change in the ship's motion state as the convergence criterion. Within each time step, the following iterative process is executed: First, the ship's boundary conditions are updated based on the ship's current six-degree-of-freedom position and attitude, and the hydrodynamic equations are solved to obtain the flow field pressure distribution; then, based on the pressure distribution, the total hydrodynamic force and torque acting on the ship are calculated by integrating over the ship's surface; next, the total hydrodynamic force and torque are used as external force inputs to solve the ship's six-degree-of-freedom motion equations, obtaining the ship's motion state at the next iteration time; the above process is repeated until the change in the ship's motion state is less than a preset convergence threshold.

[0219] The value of the preset convergence threshold is determined according to the required calculation accuracy. Preferably, the displacement convergence threshold is 10. -4 m / s, velocity convergence threshold is 10 -4 m / s. The coupled system is considered self-consistent within this time step if the following condition is met:

[0220]

[0221] in and Let be the displacement and velocity vectors for the k-th iteration, respectively. and These are the preset displacement and velocity convergence thresholds, respectively.

[0222] The time-stepping strategy uses a self-consistent convergent state as the initial condition for the next time step. Once the iterative process within a time step converges, the convergent state is used as the initial condition for the next time step, and then a higher-order time integration algorithm is used to solve the problem. The entire simulation process is repeated in this way until the entire time history is calculated.

[0223] like Figure 4 and Figure 5 As shown, a single Wigley III ( A comparison of numerical calculation results and empirical values ​​of the heave and pitch motion response under wave-facing conditions. Figure 5 It can be observed that the heave motion near the resonance point ( The numerical simulation results for the values ​​of 0.75, 0.85, and 1 are slightly lower than expected, but they agree well at other points. However, the pitch values ​​near the resonance point are lower than expected. 1. 1.25), the numerical simulation results in this paper are slightly higher than expected, but the results at other points are in good agreement;

[0224] Figure 6 and Figure 7 For a single Wigley III ( When navigating near the shore in waves (ship-to-shore distance d) t =0.3) Compared with the numerical results of the constant surface element method, the error is within a reasonable range.

[0225] In summary, although there are certain errors between the numerical simulation results and empirical values ​​presented in this paper, these errors are within a reasonable range. Therefore, the three-dimensional Rankine source high-order boundary element method proposed in this paper can be applied to calculate the motion response of a single ship navigating in waves and navigating near the shore in waves. This provides a basis for subsequent research to predict the hydrodynamic performance of two ships navigating in waves and navigating near the shore in waves.

[0226] The above description of the embodiments is only for the purpose of helping to understand the method and core idea of ​​the present invention. It should be noted that those skilled in the art can make several improvements and modifications to the present invention without departing from the principle of the present invention, and these improvements and modifications also fall within the protection scope of the claims of the present invention.

Claims

1. A time-domain numerical method for calculating the hydrodynamic performance of ships, characterized in that, Includes the following steps: S1. Establish a three-dimensional numerical model of the hull and the free liquid surface. The model is constructed based on the hull geometry, fluid domain boundary conditions and environmental parameters, and characterizes the physical response boundary of the ship under wave action. S2, based on the time-domain potential flow theory, decomposes the total velocity potential of the flow field into multiple physical components, solves the components related to ship motion and wave scattering, and pre-determines the remaining components through physical laws. S3. Based on the dynamic variables determined in S2, construct the boundary integral equations and discretize them to form a simultaneous algebraic matrix system. S4, during the time-domain propagation process, a buffer mechanism is introduced to smooth the initial transient response, and a high-order time integration algorithm is used to synchronously iteratively solve the free surface evolution and the six-degree-of-freedom motion equations of the ship. S5, based on the flow field pressure distribution after the convergence of S4 iteration, performs integral calculations on the wetted surface of the hull to calculate the total hydrodynamic force and total torque on the ship, and feeds this load back into the ship's six-degree-of-freedom motion equations to update its motion state, driving the dynamic solution of the next moment until the full-time domain simulation is completed.

2. The time-domain numerical method for calculating the hydrodynamic performance of a ship according to claim 1, characterized in that: In step S2, the decomposition of the total velocity potential of the flow field includes: background inflow, external wave incidence, wave generation from ship movement, and unsteady scattering potential. The background flow and external wave incidence are known input conditions. The wave-generating potential of a moving ship and the unsteady scattering potential are solved in real time as dynamic variables.

3. The time-domain numerical method for calculating the hydrodynamic performance of a ship according to claim 1, characterized in that: The boundary integral equation includes the wetted surface boundary of the hull, the free liquid surface boundary, and the confined water area boundary; Each boundary is meshed using high-order discretization.

4. The time-domain numerical method for calculating the hydrodynamic performance of a ship according to claim 1, characterized in that: The restricted water boundary includes the seabed and the shore. The boundary is processed by introducing a mirror source effect in the integral kernel function to simulate the reflection and constraint effect of the boundary on the flow field.

5. The time-domain numerical method for calculating the hydrodynamic performance of a ship according to claim 1, characterized in that: The process of meshing each boundary using high-order discretization is specifically as follows: The wetted surface of the hull, the free liquid surface, and the shore wall are divided into multiple small regions. Each small region is used as a calculation unit. The average velocity potential and normal rate of change of its surface jointly characterize the local fluid behavior, thereby transforming the continuous integral equation into a set of algebraic matrix systems that can be solved by a computer. The boundary treatment of the free liquid surface adopts linear or nonlinear kinetic energy conservation kinematic boundary conditions, and combines high-order extrapolation techniques of normal velocity potential to suppress numerical dissipation.

6. The time-domain numerical method for calculating the hydrodynamic performance of a ship according to claim 1, characterized in that: In step S4, the buffering mechanism is specifically as follows: during the initial time period of calculation, the scattering potential generated by the ship motion is weighted by multiplying it by a control function that monotonically increases with time until the control function tends to stabilize, and the scattering potential participates in the subsequent calculation with its full amplitude.

7. The time-domain numerical method for calculating the hydrodynamic performance of a ship according to claim 1, characterized in that: In step S4, the synchronous iterative solution specifically involves: Within each time step, the free surface height, velocity potential and its normal derivative, as well as the ship's displacement, velocity and acceleration are solved sequentially through a multi-stage calculation process. The intermediate calculation results of each stage are then linearly combined according to a preset convergence weighting coefficient to generate the updated values ​​for the next time step.

8. The time-domain numerical method for calculating the hydrodynamic performance of a ship according to claim 1, characterized in that: The six-degree-of-freedom equations of motion for the ship are constructed based on the Newton-Euler equations. The six-degree-of-freedom equations of motion for the ship include a mass matrix and a restoring force matrix; The mass matrix and restoring force matrix are calculated based on the total mass of the hull, the coordinates of the center of gravity, the inertia tensor, and the still water buoyancy distribution.

9. A time-domain numerical method for calculating the hydrodynamic performance of a ship according to claim 1, characterized in that: The synchronous iterative solution described in step S4 and the feedback update described in step S5 constitute a closed feedback loop; Specifically, within each time step, the six-degree-of-freedom motion state of the ship directly affects the setting of the flow field boundary conditions, and the hydrodynamic loads calculated from the flow field pressure distribution directly drive the update of the ship's motion equations. The two achieve state self-consistency through repeated coupled calculations within the time step until the entire simulation time history is completed.

10. A time-domain numerical method for calculating the hydrodynamic performance of a ship according to claim 9, characterized in that: Within each time step, perform the following iterative process: 1) Based on the ship's current six-degree-of-freedom position and attitude, update the hull boundary conditions and solve the hydrodynamic equations to obtain the flow field pressure distribution; 2) Based on the pressure distribution, calculate the total hydrodynamic force and total torque acting on the ship by integrating over the hull surface; 3) Using the total hydrodynamic force and total torque as external force inputs, solve the six-degree-of-freedom motion equations of the ship to obtain the ship's motion state at the next iteration time; 4) Repeat steps 1) to 3) until the change in the ship's motion state is less than the preset convergence threshold, and determine that the coupled system has reached self-consistency within this time step; 5) Using the self-consistent state as the initial condition for the next time step, a higher-order time integration algorithm is used to solve the problem.

Images