An underwater object trajectory prediction method based on engineering algorithm

Abstract

This invention relates to the fields of aerospace and shipbuilding, and discloses a method for predicting the trajectory of underwater objects launched based on engineering algorithms. The method includes: collecting and calculating the necessary parameters; establishing a multibody dynamics description of a cylindrical object exiting the cylinder using a purely elastic contact model; establishing a six-degree-of-freedom equation of motion for the cylindrical object, simplifying it into a center-of-mass motion equation and a motion equation about the center-of-mass, wherein the motion equation about the center-of-mass is established using the quaternion method; establishing a force model of the object in the water entry section based on empirical formulas; establishing an underwater force model of the cylindrical object based on gravity, buoyancy, position force, and inertial force; determining the stage based on the dynamic parameters of the cylindrical object exiting the cylinder, selecting the appropriate model to calculate the force situation of the cylindrical object, and substituting it into the six-degree-of-freedom equation of motion to solve using the RK5 method; determining whether the object has reached the target depth; and outputting the final result if the target depth has been reached. This method avoids the complexity of solving the Navier-Stokes (N-S) equations and significantly improves the calculation speed.

Description

Technical Field

[0001] This invention relates to the fields of aviation and navigation, and shipbuilding, and in particular to a method for predicting the trajectory of underwater cylindrical objects based on engineering algorithms in the field of shipbuilding. Background Technology

[0002] With the continuous development of marine technology, accurate prediction of the trajectory of underwater objects has become a crucial technical requirement. This requirement is particularly prominent in the field of underwater exploration, where accurate trajectory prediction has a significant impact on improving exploration efficiency and accuracy.

[0003] Currently, underwater trajectory simulation and prediction primarily rely on computational fluid dynamics (CFD) numerical simulation techniques. CFD technology can provide accurate underwater trajectory predictions by simulating fluid flow and object interactions. However, this technology has extremely high requirements for time and computational resources during the calculation process, limiting its application in fast-response and resource-constrained environments.

[0004] In some scenarios, the requirements for prediction accuracy may be relatively low, while the requirements for time consumption are more stringent. Furthermore, the flexibility of CFD technology in dealing with unexpected problems is also limited. Therefore, developing an underwater trajectory prediction method that is computationally efficient, has a fast response speed, and can operate stably in variable environments is of great significance for meeting the needs of modern maritime military technology. Summary of the Invention

[0005] In view of this, the purpose of this invention is to propose an underwater cylindrical object trajectory prediction software algorithm based on engineering algorithms to solve the problem of quickly calculating the underwater motion trajectory of simulated objects.

[0006] To achieve the above objectives, this invention provides a method for predicting the trajectory of underwater objects based on engineering algorithms, the method comprising the following steps:

[0007] S1. Collect the parameters required for calculation, including parameters of the discharge cylinder, cylindrical object, seawater, and calculation parameters;

[0008] S2. Using a purely elastic contact model, a multibody dynamics description of the cylindrical object exiting the cylinder is established based on Newton-Euler equations, and the dynamic parameters of the cylindrical object exiting the cylinder are calculated.

[0009] S3. Establish the six-degree-of-freedom motion equations of the cylindrical object and simplify them into the motion equations around the center of mass and the motion equations about the center of mass. The quaternion method is used to transform the motion equations about the center of mass.

[0010] S4. Based on empirical formulas, simplify the process of a cylindrical object entering water and establish a force model of the object in the water section.

[0011] S5. Based on gravity, buoyancy, position force, and inertial force, establish an underwater force model for a cylindrical object;

[0012] S6. Determine the stage based on the dynamic parameters of the cylindrical object when it exits the cylinder. If it is in the water entry stage, select the force model of the object in the water entry stage. If it is in the underwater stage, select the underwater force model of the cylindrical object. Calculate the force situation of the cylindrical object and substitute it into the six-degree-of-freedom motion equation.

[0013] S7. Solve the six-degree-of-freedom motion equations based on the RK5 method to obtain the dynamic parameters of the cylindrical object;

[0014] S8. Based on the dynamic parameters of the cylindrical object, determine whether the object has reached the target depth. If the target depth has been reached, output the final result; otherwise, repeat steps S6 and S7 until the target depth is reached.

[0015] Preferably, step S2 involves using a purely elastic contact model to establish a multibody dynamics description of the cylindrical object exiting the cylinder based on Newton's and Euler's equations, and calculating the dynamic parameters of the cylindrical object exiting the cylinder, including the following steps:

[0016] S2.1. The classic Newton-Euler equations are used to accurately describe the motion of the center of mass of a cylindrical object and its changing posture relative to the center of mass.

[0017] S2.2. For the contact force generated between the cylindrical object and the feeding cylinder during the contact process, a purely elastic contact model is adopted, in which the contact force and the indentation depth have a non-linear relationship:

[0018] F n =Kδ n

[0019] In the formula, F n Let δ represent the contact force, the amount of indentation between the discharge cylinder and the cylindrical object, K represent the stiffness coefficient, and n be the force exponent, determined by the material of the contacting bodies and the geometric properties of the contact area. The formula for calculating the stiffness coefficient K is as follows:

[0020]

[0021] In the formula, R i R j σ is the radius of curvature at the point of contact between the feeding tube and the cylindrical object; i σ j The parameters related to the material of the feeding cylinder and the cylindrical object are calculated using the following formulas:

[0022]

[0023] In the formula, E I υ IFor Poisson's ratio and Young's modulus of the material;

[0024] S2.3. For the dynamic model of the cylindrical object exiting the cylinder, the effects of gravity and the contact force provided by the cylinder wall constraint on the cylindrical object are fully considered. Based on Newton's equations, the dynamic equations for the motion of the center of mass of the cylindrical object are established as follows:

[0025]

[0026] Where m represents the mass of the cylindrical object; ω0 is the second derivative of the axial displacement of the center of mass of the cylindrical object in the global coordinate system; ω0 is the rotational velocity of the throwing cylinder; t is the current physical time; F is the resultant force of the contact force and gravity acting on the object in the axial direction.

[0027] Furthermore, to obtain the dynamic equations for the changing posture relative to the center of mass of the cylindrical object, the following equations are formulated based on the Euler equations:

[0028]

[0029] In the formula, J1, J2, and J3 are the three principal moments of inertia of the cylindrical object in the central principal axis coordinate system; ω1, ω2, and ω3 are the three components of the angular velocity vector of the cylindrical object in the object coordinate system. These are the first derivatives of the three components; M1, M2, and M3 are the three components of the torque vector acting on the cylindrical object relative to the center of mass in the cylindrical object's coordinate system; when r c When the length L of the feeding tube is greater than the length of the feeding tube, the process of the cylindrical object exiting the tube ends.

[0030] Preferably, in step S3, the six-degree-of-freedom motion equations of the cylindrical object are established and simplified into the motion equations around the center of mass and the motion equations about the center of mass. The motion equations about the center of mass are established using the quaternion method, as detailed below:

[0031] S3.1 First, define a coordinate system to determine the direction of relative position, velocity, and external force vector parameters. The coordinate system includes the ground coordinate system oxyz and the object coordinate system o1x1y1z1.

[0032] For the coordinate axis transformation, the ground coordinate system eventually coincides with the object coordinate system through a finite number of rotations around the axis. The angles of rotation around the axis are defined as the attitude angles of the cylindrical object, denoted as roll angle γ, yaw angle φ, and pitch angle θ, respectively. The order of rotation transformation is set as zyx, resulting in the following transformation equations:

[0033]

[0034] Among them, R z (φ), R y (γ), Rx (θ) are the transformation matrices for rotations about the z, y, and x axes, respectively;

[0035] S3.2 The six-degree-of-freedom dynamic equations of a cylindrical object are:

[0036] ma = F total

[0037]

[0038] Where m is the mass of the cylinder, a is the acceleration of the cylinder, and J is the moment of inertia of the cylinder. F is the first derivative of the rotational angular velocity of the cylindrical object. total M is the net force acting on the cylindrical object. total The net torque acting on the cylindrical object;

[0039] S3.3 The equation of motion for the center of mass of a cylindrical object is:

[0040] In the ground coordinate system, the motion components of the cylindrical object are:

[0041]

[0042] Among them, v x v y v z Let x, y, z be the velocity of the cylindrical object relative to the ground coordinate system. By integrating, we can obtain the position (x, y, z) of the object's center of mass relative to the ground coordinate system.

[0043] S3.4 To determine the attitude angle of the cylindrical object relative to the Earth, the equation of motion of the cylindrical object about its center of mass is established using the quaternion method:

[0044] R o =R z (φ)R y (γ)R x (θ)

[0045] Use the following formula from the equation of motion about the center of mass R o Calculate the quaternion q:

[0046]

[0047] Where, r ij The equation of motion R about the center of mass oIn the quaternion, w is the real part of the quaternion, corresponding to the scaling factor of the rotation. When describing the rotation, w is related to the angle of rotation and the direction cosine of the rotation axis. x is the first imaginary part of the quaternion, related to the x-component of the rotation axis. y is the second imaginary part of the quaternion, related to the y-component of the rotation axis. z is the third imaginary part of the quaternion, related to the z-component of the rotation axis.

[0048] The quaternion differential equation is obtained as follows:

[0049]

[0050] In the formula ω x ω y ω z This represents the component of the rigid body's angular velocity in the inertial coordinate system.

[0051] Preferably, in step S4, based on empirical formulas, the process of the cylindrical object entering the water is simplified, and a force model of the object in the water section is established, as follows:

[0052] S4.1 Formulas for calculating buoyancy and center of buoyancy of a cylindrical object during its immersion in water, in a ground coordinate system:

[0053] B=ρgV

[0054]

[0055] In the formula, B is buoyancy, and X... B It is the center of gravity, L w ρ is the length of the cylindrical object submerged in water, V is the volume of the cylindrical object submerged in water, ρ is the density of seawater, g is the acceleration due to gravity, and R is the radius of the cylindrical object.

[0056] S4.2 The hydrodynamics of the water-immersion process adopts a semi-empirical formula as the basic model for calculating the hydrodynamics of a cylindrical object during the water-immersion process. The formula is as follows:

[0057] Fz=-0.0076ρπR 2 v 2 (10-x)

[0058]

[0059] In the formula, x is the axial distance from the top of the cylindrical object to the water surface during the water discharge process. Let be the dimensionless quantity of the angular velocity of the cylinder, v be the velocity of the cylinder relative to the water flow, α be the angle of attack of the cylinder, Fz be the axial force on the cylinder, Fy be the normal force on the cylinder, and M be the torque on the cylinder.

[0060] The resultant force of the force obtained from the empirical formula and the buoyancy force is used as the force on the cylindrical object during the water entry phase.

[0061] Preferably, in step S5, based on the theoretical algorithms for gravity, buoyancy, position force, and inertial force, the specific steps for establishing the underwater force model of the cylindrical object are as follows:

[0062] S5.1 The weight of the cylindrical object is:

[0063] G = mg

[0064] The direction of gravity is perpendicular to the ground, and gravity does not produce torque;

[0065] S5.2 When a cylindrical object is completely submerged in water, the buoyancy will no longer change. The calculation formula is as follows:

[0066] B=ρgV

[0067] Where V is the volume of the cylindrical object submerged in water, ρ is the density of seawater, and g is the acceleration due to gravity.

[0068] The direction of buoyancy is always vertically upward, and the torque generated by buoyancy is as follows:

[0069]

[0070] Where lx, ly, and lz are the distances between the center of buoyancy and the center of mass in the three directions;

[0071] S5.3 The expression for position force under small angle of attack is as follows:

[0072]

[0073] v = v i -V sea

[0074] Among them, V sea Let v be the velocity of the seawater current, and v be the velocity of the cylindrical object relative to the water current. i Let be the velocity of the cylindrical object relative to the ground, and its expression is as follows:

[0075]

[0076] Where V1, V2, and V3 are the set seawater flow velocities, and x0 to x3 are the set range of seawater flow velocities for each layer;

[0077] After dimensionless transformation, the following expression is obtained:

[0078]

[0079] In the formula, α, β, and δ represent the angle of attack, sideslip angle, and roll angle of the cylindrical object, respectively.

[0080] Simplifying the above equation, we get:

[0081] C l =C0+C α α+C δ δ

[0082] M l =M0+M β β+M δ δ

[0083] In the formula, C0 is the normal force coefficient when α and δ are equal to zero, C α C is the position derivative of the normal force coefficient with respect to the angle of attack; δ M0 is the positional derivative of the normal force coefficient with respect to the roll angle δ, and M0 is the normal moment coefficient when β and δ are equal to zero. β M is the position derivative of the normal moment coefficient with respect to the sideslip angle β; δ This is the positional derivative of the normal moment coefficient with respect to the roll angle δ.

[0084] The calculation of the position force coefficient for cylindrical objects under medium or large angles of attack is handled specially, as follows:

[0085] Normal force F of a cylindrical object NB The calculation formula is:

[0086]

[0087] In the formula, S max Let S be the maximum cross-sectional area of ​​the cylindrical object, A be the area of ​​the base of the cylindrical object, and S be the area of ​​the cylinder. p Let F be the area of ​​the longitudinal projection plane of the cylindrical object, η be the crossflow drag scaling factor, and F be the crossflow drag scaling factor. dn For crossflow resistance;

[0088] Under medium or high angle of attack conditions:

[0089] F D =F NB sinα

[0090] F L =F NB cosα

[0091]

[0092] In the formula, L is the length of the cylindrical object, and x m The position of the center of mass of the cylindrical object, x c The position of the center of the longitudinal projection plane of the cylindrical object;

[0093] S5.4 Inertial forces include additional inertial forces and additional inertial moments, and their expressions are as follows:

[0094]

[0095] in,

[0096]

[0097] In the formula, a and F is the first derivative of the object's acceleration and angular velocity. added To add inertial force, M added To add an inertial torque, m added For added mass, I added To add the moment of inertia, ρ object The density of the cylindrical object.

[0098] Preferably, in step S7, the quaternion differential equation is solved based on the RK5 method to obtain the rotation parameters of the cylindrical object, as follows:

[0099] S7.1 The expressions for the system of differential equations are as follows:

[0100]

[0101] The solution formula for the RK5 method is:

[0102]

[0103] in,

[0104] k1=h·f(t n ,y n )

[0105]

[0106] k5=h·f(t n +h,y n -5k1+3k4)

[0107]

[0108] The error is: Where h is the iteration step size.

[0109] The beneficial effects of this invention are as follows: (1) This invention proposes an underwater project trajectory prediction method based on engineering algorithms. This method simplifies the calculation process and avoids the complexity of solving the Navier-Stokes (NS) equations, thereby significantly improving the calculation speed. While maintaining reasonable prediction accuracy, the method of this invention can quickly adapt to different underwater environments and conditions, effectively solving the limitations of existing technologies.

[0110] (2) Based on engineering algorithms, this invention eliminates the need to solve the Navier-Stokes equations, thereby greatly reducing the computational cost of calculating the underwater motion trajectory of an object.

[0111] (3) The present invention simplifies the water entry problem, and all calculation processes adopt engineering algorithms, which can greatly reduce the time cost of calculating the underwater motion trajectory of objects.

[0112] (4) The present invention has great universality in dealing with problems, and thus can efficiently and quickly calculate the underwater trajectory of a cylindrical object in various sea conditions.

[0113] (5) The present invention is currently very much in line with the needs of various countries in the field of underwater motion trajectory prediction, and has great strategic value in the field of marine exploration. Attached Figure Description

[0114] Figure 1 This refers to the software analysis process.

[0115] Figure 2 The graph shows the change of the z-direction velocity of the object in the vertical stationary cylinder in Example 1 over time (the velocities in the other two directions are 0).

[0116] Figure 3 This is a graph showing the change in the y-direction displacement of an object falling in water with the depth of fall, in comparison with the prior art in Example 1.

[0117] Figure 4 This is a three-dimensional displacement diagram of the object falling in water in Example 1.

[0118] Figure 5 The graph shows the change of the velocity of the object in the z and x directions over time in Example 2 (the velocity in other directions is 0).

[0119] Figure 6 This is a graph showing the change in the y-direction displacement of an object falling in water with the depth of fall, in comparison with the prior art in Example 2.

[0120] Figure 7 This is a three-dimensional displacement diagram of an object falling in water as shown in Example 2. Detailed Implementation

[0121] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0122] It should be noted that, unless otherwise defined, the technical or scientific terms used in the embodiments of this invention should have the ordinary meaning understood by those skilled in the art to which this invention pertains. The terms "first," "second," and similar terms used in this invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed following the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0123] The embodiments of the present invention will be further described in detail below with reference to the accompanying drawings.

[0124] S1. Input the parameters required for calculation, including the parameters of the sample tube, the parameters of the cylindrical object (such as an underwater sampler or an underwater sonar scanner), the seawater parameters, and the calculation parameters (such as the iteration time).

[0125] In this embodiment, based on this algorithm, a software program for calculating the underwater motion trajectory of a cylindrical object was developed using C++ on the Visual Studio platform. The software analysis process is as follows: Figure 1 As shown.

[0126] First, you need to input the parameters of the feeding cylinder and the cylindrical object to calculate the motion state of the cylindrical object when it leaves the feeding cylinder. Figure 2 and Figure 3 As shown, this serves as the initial state for subsequent hydrodynamic analysis; then, seawater parameters and calculation parameters need to be input, and plotted as shown... Figure 4 The three-dimensional displacement diagram shown here has a target depth of 60m in this example.

[0127] S2. Using a purely elastic contact model, a multibody dynamics description of the cylindrical object exiting the cylinder is established based on Newton-Euler equations, and the dynamic parameters of the cylindrical object exiting the cylinder are calculated.

[0128] In S2, construct a dynamic model of the cylindrical object in the feeding cylinder and solve for the dynamic parameters of the cylindrical object when it exits the cylinder, including the following steps:

[0129] S2.1. The classic Newton-Euler equations are used to accurately describe the motion of the center of mass of a cylindrical object and its changing posture relative to the center of mass.

[0130] S2.2. For the contact force generated between the cylindrical object and the delivery cylinder during the contact process, a purely elastic contact model is adopted, in which the contact force and the indentation depth have a non-linear relationship, as follows:

[0131] F n =Kδ n

[0132] In the formula, F n The contact force is δ, which represents the amount of indentation between the discharge cylinder and the cylindrical object; K represents the stiffness coefficient; n is the force exponent, which defaults to 1.5 and is determined by the material of the contacting bodies and the geometric properties of the contact area; the formula for calculating the stiffness coefficient K is as follows:

[0133]

[0134] In the formula, R i R j σ is the radius of curvature at the point of contact between the feeding tube and the cylindrical object; i σ j The parameters related to the material of the feeding cylinder and the cylindrical object are calculated using the following formulas:

[0135]

[0136] In the formula, E I υ I For Poisson's ratio and Young's modulus of the material;

[0137] S2.3. For the dynamic model of the cylindrical object exiting the cylinder, the effects of gravity and the contact force provided by the cylinder wall constraint on the cylindrical object are fully considered. Based on Newton's equations, the dynamic equations for the motion of the center of mass of the cylindrical object are established as follows:

[0138]

[0139] Where m represents the mass of the cylindrical object; ω0 is the second derivative of the axial displacement of the center of mass of the cylindrical object in the global coordinate system; ω0 is the rotational velocity of the cartridge; t is the current physical time; F is the resultant force of the contact force and gravity acting on the object in the axial direction.

[0140] Furthermore, to obtain the dynamic equations for the changing posture relative to the center of mass of the cylindrical object, the following equations are formulated based on the Euler equations:

[0141]

[0142] In the formula, J1, J2, and J3 are the three principal moments of inertia of the cylindrical object in the central principal axis coordinate system; ω1, ω2, and ω3 are the three components of the angular velocity vector of the cylindrical object in the object coordinate system. Let M1, M2, and M3 be the first derivatives of the three components; M1, M2, and M3 are the three components of the torque vector acting on the cylindrical object relative to the center of mass in the cylindrical object's coordinate system; when r c When the length L of the feeding tube is greater than the length of the feeding tube, the process of the object exiting the tube ends.

[0143] S3. Establish the six-degree-of-freedom equations of motion for the cylindrical object, including the equations of motion around the center of mass and the equations of motion about the center of mass. The equations of motion about the center of mass are established using the quaternion method, as detailed below:

[0144] S3.1 First, define a coordinate system to determine the direction of relative position, velocity, and external force vector parameters. The coordinate system includes the ground coordinate system oxyz and the object coordinate system o1x1y1z1.

[0145] The ground coordinate system oxyz has a point on the ground (usually the launch position of the cylindrical object) as its origin, the oz axis is in the vertical plane and upward is defined as the positive direction; the oy axis is along the direction of the ship's movement; and the ox axis is determined according to the right-hand rule.

[0146] The object's coordinate system is o1x1y1z1, with the center of mass of the cylindrical object as its origin o1. The o1x1 axis coincides with the longitudinal axis of the cylindrical object, and the positive direction of the axis is defined as the positive direction of the cylindrical object. The o1y1 axis is perpendicular to the axis of the object and lies in its longitudinal plane, with upward being positive. The direction of the o1z1 axis can be determined using the right-hand rule. The object's coordinate system moves with the movement of the cylindrical object.

[0147] For the coordinate axis transformation, the ground coordinate system eventually coincides with the object coordinate system through a finite number of rotations around the axis. The angles of rotation around the axis are defined as the attitude angles of the cylindrical object, denoted as roll angle γ, yaw angle φ, and pitch angle θ, respectively. The order of the rotation transformation is set as ZYX, resulting in the following transformation equations:

[0148]

[0149] Among them, R z (φ), R y (γ), R x (θ) are the transformation matrices for rotations about the z, y, and x axes, respectively;

[0150] S3.2 The six-degree-of-freedom dynamic equations of a cylindrical object are:

[0151] ma = F total

[0152]

[0153] Where m is the mass of the cylinder, a is the acceleration of the cylinder, and J is the moment of inertia of the cylinder. F is the first derivative of the rotational angular velocity of the cylindrical object. total M is the net force acting on the cylindrical object. total The net torque acting on the cylindrical object;

[0154] S3.3 The equation of motion for the center of mass of a cylindrical object is:

[0155] In the ground coordinate system, the motion components of the cylindrical object are:

[0156]

[0157] Among them, v x v y v z Let x, y, z be the velocity of the cylindrical object relative to the ground coordinate system. By integrating, we can obtain the position (x, y, z) of the object's center of mass relative to the ground coordinate system.

[0158] S3.4 To determine the attitude angle of the cylindrical object relative to the Earth, the equation of motion of the cylindrical object about its center of mass is established using the quaternion method:

[0159] R o =R z (φ)R y (γ)R x (θ)

[0160] Use the following formula from the equation of motion about the center of mass R o Calculate the quaternion q:

[0161]

[0162] Where, r ij The equation of motion R about the center of mass o In the quaternion, w is the real part of the quaternion, corresponding to the scaling factor of the rotation. When describing the rotation, w is related to the angle of rotation and the direction cosine of the rotation axis. x is the first imaginary part of the quaternion, related to the x-component of the rotation axis. y is the second imaginary part of the quaternion, related to the y-component of the rotation axis. z is the third imaginary part of the quaternion, related to the z-component of the rotation axis.

[0163] The quaternion differential equation is obtained as follows:

[0164]

[0165] In the formula ω x ω y ω z The components of the rigid body angular velocity in the inertial coordinate system

[0166] S4. Based on empirical formulas, simplify the process of a cylindrical object entering water and establish a force model of the object in the water entry section, specifically as follows:

[0167] S4.1 Formulas for calculating buoyancy and center of buoyancy of a cylindrical object during its immersion in water, in a ground coordinate system:

[0168] B=ρgV

[0169]

[0170] In the formula, B is buoyancy, and X... B It is the center of gravity, L w ρ is the length of the cylindrical object submerged in water, V is the volume of the cylindrical object submerged in water, ρ is the density of seawater, g is the acceleration due to gravity, and R is the radius of the cylindrical object.

[0171] S4.2 The hydrodynamics of the water-immersion process adopts a semi-empirical formula as the basic model for calculating the hydrodynamics of a cylindrical object during the water-immersion process. The formula is as follows:

[0172] Fz=-0.0076ρπR 2 v 2 (10-x)

[0173]

[0174] In the formula, x is the axial distance from the top of the cylindrical object to the water surface during the water discharge process, and R is the radius of the cylindrical object. Let be the dimensionless quantity of the angular velocity of the cylinder, v be the velocity of the cylinder relative to the water flow, Fz be the axial force on the cylinder, Fy be the normal force on the cylinder, and M be the torque on the cylinder.

[0175] The resultant force of the force obtained from the empirical formula and the buoyancy force is used as the force on the cylindrical object during the water entry phase.

[0176] S5. The underwater force situation of the cylindrical object is divided into gravity, buoyancy, position force, and inertial force. A theoretical algorithm is established to create an underwater force model of the cylindrical object, and each component is solved to obtain the underwater force distribution. The specific steps are as follows:

[0177] S5.1 The weight of the cylindrical object is:

[0178] G = mg

[0179] The direction of gravity is perpendicular to the ground, and gravity does not produce torque;

[0180] S5.2 When a cylindrical object is completely submerged in water, the buoyancy will no longer change. The calculation formula is as follows:

[0181] B=ρgV

[0182] The direction of buoyancy is always vertically upward, and the torque generated by buoyancy is as follows:

[0183]

[0184] Where lx, ly, and lz are the distances between the center of buoyancy and the center of mass in the three directions;

[0185] S5.3 The expression for position force under small angle of attack is as follows:

[0186]

[0187] v = v i -V sea

[0188] Among them, V sea v is the velocity of the seawater current. i The expression for the velocity of the cylindrical object relative to the ground is as follows:

[0189]

[0190] Where V1, V2, and V3 are the set seawater flow velocities, and x0 to x3 are the set range of seawater flow velocities for each layer;

[0191] After dimensionless transformation, the following expression is obtained:

[0192]

[0193] In the formula, α, β, and δ represent the angle of attack, sideslip angle, and roll angle of the cylindrical object, respectively.

[0194] Simplifying the above equation, we get:

[0195] C l =C0+C α α+C δ δ

[0196] M l =M0+M β β+M δ δ

[0197] In the formula, C0 is the normal force coefficient when α and δ are equal to zero, C α C is the position derivative of the normal force coefficient with respect to the angle of attack; δ M0 is the positional derivative of the normal force coefficient with respect to the roll angle δ; M0 is the normal moment coefficient when β and δ are equal to zero. β M is the position derivative of the normal moment coefficient with respect to the sideslip angle β; δ This is the positional derivative of the normal moment coefficient with respect to the roll angle δ;

[0198] The calculation of the position force coefficient for cylindrical objects under medium or large angles of attack is handled specially, as follows:

[0199] Normal force F of a cylindrical object NB The calculation formula is:

[0200]

[0201] In the formula, S max Let S be the maximum cross-sectional area of ​​the cylindrical object, A be the area of ​​the base of the cylindrical object, and S be the area of ​​the cylinder. p Let F be the area of ​​the longitudinal projection plane of the cylindrical object, η be the crossflow drag scaling factor, and F be the crossflow drag scaling factor. dn For crossflow resistance;

[0202] Under medium or high angle of attack conditions:

[0203] F D =F NB sinα

[0204] F L =F NB cosα

[0205]

[0206] In the formula, L is the length of the cylindrical object, and x m The position of the center of mass of the cylindrical object, x c The position of the center of the longitudinal projection plane of the cylindrical object.

[0207] S5.4 Inertial forces include additional inertial forces and additional inertial moments, and their expressions are as follows:

[0208]

[0209] in,

[0210]

[0211] In the formula, a and F is the first derivative of the object's acceleration and angular velocity. added To add inertial force, M added To add an inertial torque, m added For added mass, I added To add the moment of inertia, ρ object The density of the cylindrical object.

[0212] S7. Based on the RK5 method, the quaternion differential equation is solved to obtain the rotation parameters of the cylindrical object, as follows:

[0213] S7.1 The expressions for the system of differential equations are as follows:

[0214]

[0215] The solution formula for the RK5 method is:

[0216]

[0217] in,

[0218] k1=h·f(t n ,y n )

[0219]

[0220] k5=h·f(t n +h,y n -5k1+3k4)

[0221]

[0222] The error is: Where h is the iteration step size.

[0223] S8. Based on the dynamic parameters of the cylindrical object, determine whether the object has reached the target depth. If the target depth has been reached, output the final result; otherwise, repeat steps 6 and 7.

[0224] Example 1

[0225] The scenario involves a vertically downward drop of an object from a stationary drop tube. The tube is 2 meters long, and the object weighs 133 kg. The seawater velocity distribution is as follows: 0 m at depths of 0-10 m; 10 m / s at depths of 10-20 m, in both directions along the y-axis; and 30 m / s from depth 20 m to the bottom, in the positive y-axis direction. The change in the object's falling velocity as it exits the tube is shown below. Figure 2 (The velocity in the other two directions is 0). The specific input parameters are shown in Table 1 below. The displacement of the object in the y-axis direction is compared with the calculation results of the commercial software STARCCM. Figure 3 The calculation error and time consumption are shown in Table 2.

[0226] Table 1 Input Parameter Table for Example 1

[0227] Material parameters of cylindrical objects <![CDATA[σ j ]]> 200Gpa Feeding cylinder rotation speed <![CDATA[ω0]]> 0 rad / s Length of feeding tube <![CDATA[L i ]]> 2m Mass of a cylindrical object m 133kg First layer of seawater depth x1 -10m The speed of the first layer of seawater <![CDATA[V1]]> 0m / s Second layer of seawater depth x2 -20m Second layer of seawater speed <![CDATA[V2]]> 10m / s Third layer of seawater depth x3 -60m The speed of the third layer of seawater <![CDATA[V3]]> 30m / s Iteration time step h 0.001s

[0228] Table 2 Comparison of Results

[0229] Time consumption 3.82h 0.25s 3.82h Final displacement calculation results 66.72m 64.12m 2.60m

[0230] Example 2

[0231] A projectile cylinder is selected to undergo a regular oscillating motion with an oscillation angle of 30°, along the positive x-axis, and an oscillation frequency of 0.1 rad / s. The cylinder length is 2 m, and the mass of the object is 133 kg. The seawater velocity distribution is as follows: 0 m / s between depths of 0-10 m; 10 m / s between depths of 10-20 m, along both the positive and negative y-axis; and 30 m / s between depths of 20 m and the bottom, along the negative y-axis. The input parameters are shown in Table 3 below. The changes in the falling velocity and the y-direction velocity of the object during its exit from the cylinder are obtained as follows: Figure 5 (Velocity in other directions is 0), the displacement of the object in the y-axis direction is as follows: Figure 6 The calculation error and time consumption are shown in Table 4.

[0232] Table 3 Input Parameter Table for Example 2

[0233] Material parameters of cylindrical objects <![CDATA[σ j ]]> 200Gpa Feeding cylinder rotation speed <![CDATA[ω0]]> 0.1 rad / s Length of feeding tube <![CDATA[L i ]]> 2m Mass of a cylindrical object m 133kg First layer of seawater depth x1 -10m The speed of the first layer of seawater <![CDATA[V1]]> 0m / s Second layer of seawater depth x2 -20m Second layer of seawater speed <![CDATA[V2]]> 10m / s Third layer of seawater depth x3 -60m The speed of the third layer of seawater <![CDATA[V3]]> -30m / s Iteration time step h 0.001s

[0234] Table 4 Comparison of Results

[0235] Time consumption 4.62h 0.28s 4.62h Final displacement calculation results -66.01m -59.77m 6.24m

Claims

1. A method for predicting the trajectory of underwater objects based on engineering algorithms, characterized in that, The method includes the following steps: S1. Collect the parameters required for calculation, including parameters of the discharge cylinder, cylindrical object, seawater, and calculation parameters; S2. Using a purely elastic contact model, a multibody dynamics description of the cylindrical object exiting the cylinder is established based on Newton's and Euler's equations, and the dynamic parameters of the cylindrical object exiting the cylinder are calculated; specifically, the following steps are included: S2.

1. The classic Newton-Euler equations are used to accurately describe the motion of the center of mass of a cylindrical object and its changing posture relative to the center of mass. S2.

2. For the contact force generated between the cylindrical object and the feeding cylinder during the contact process, a purely elastic contact model is adopted, in which the contact force and the indentation depth have a non-linear relationship: ； In the formula, F n Let δ represent the contact force, the amount of indentation between the discharge cylinder and the cylindrical object, K represent the stiffness coefficient, and n be the force exponent, determined by the material of the contacting bodies and the geometric properties of the contact area. The formula for calculating the stiffness coefficient K is as follows: ； In the formula, , The radius of curvature at the point of contact between the feeding tube and the cylindrical object; , The parameters related to the material of the feeding cylinder and the cylindrical object are calculated using the following formulas: ； In the formula, , For Poisson's ratio and Young's modulus of the material; S2.

3. For the dynamic model of the cylindrical object exiting the cylinder, the effects of gravity and the contact force provided by the cylinder wall constraint on the cylindrical object are fully considered. Based on Newton's equations, the dynamic equations for the motion of the center of mass of the cylindrical object are established as follows: ； Where m represents the mass of the cylindrical object; Let be the second derivative of the axial displacement of the center of mass of the cylindrical object in the global coordinate system; t is the rotational speed of the discharge tube; t is the current physical time; F is the resultant force of the contact force and gravity acting on the object along the axis. Furthermore, to obtain the dynamic equations for the changing posture relative to the center of mass of the cylindrical object, the following equations are formulated based on the Euler equations: ； In the formula, , , These are the three principal moments of inertia of the cylindrical object in the central principal axis coordinate system; , , These are the three components of the angular velocity vector of a cylindrical object in the object's coordinate system. , , These are the first derivatives of the three components; , , These are the three components of the torque vector acting on a cylindrical object relative to its center of mass in the cylindrical object's coordinate system. Let be the axial displacement of the center of mass of the cylindrical object in the global coordinate system; when Greater than the length L of the feeding tube i At that time, the process of the cylindrical object exiting the cylinder ends; S3. Establish the six-degree-of-freedom motion equations of the cylindrical object and simplify them into the motion equations around the center of mass and the motion equations about the center of mass. The quaternion method is used to transform the motion equations about the center of mass. S4. Based on empirical formulas, simplify the process of a cylindrical object entering water and establish a force model of the object in the water section. S5. Based on gravity, buoyancy, position force, and inertial force, establish an underwater force model for a cylindrical object; S6. Determine the stage based on the dynamic parameters of the cylindrical object when it exits the cylinder. If it is in the water entry stage, select the force model of the object in the water entry stage. If it is in the underwater stage, select the underwater force model of the cylindrical object. Calculate the force situation of the cylindrical object and substitute it into the six-degree-of-freedom motion equation. S7. Solve the six-degree-of-freedom motion equations based on the RK5 method to obtain the dynamic parameters of the cylindrical object; S8. Based on the dynamic parameters of the cylindrical object, determine whether the object has reached the target depth. If the target depth has been reached, output the final result; otherwise, repeat steps S6 and S7 until the target depth is reached.

2. The underwater project trajectory prediction method based on engineering algorithms according to claim 1, characterized in that, Step S3: Establish the six-degree-of-freedom motion equations for the cylindrical object, simplifying them into equations of motion around the center of mass and equations of motion about the center of mass. The quaternion method is used to establish the equations of motion about the center of mass, as detailed below: S3.1 First, define a coordinate system to determine the direction of relative position, velocity, and external force vector parameters. The coordinate system includes the ground coordinate system oxyz and the object coordinate system o1x1y1z1. For coordinate axis transformation, the ground coordinate system eventually coincides with the object coordinate system through a finite number of rotations around the axis. The angle of rotation around the axis is defined as the attitude angle of the cylindrical object, denoted as roll angle γ and yaw angle, respectively. and pitch angle If we set the order of the rotation transformation as zyx, then we obtain the following transformation equation: ； ； ； ； in, , , These are the transformation matrices for rotations about the z, y, and x axes, respectively; S3.2 The six-degree-of-freedom dynamic equations of a cylindrical object are: ； Where m is the mass of the cylindrical object. Let J be the acceleration of the cylindrical object, and J be the moment of inertia of the cylindrical object. Let be the first derivative of the rotational angular velocity of the cylindrical object. The net force acting on the cylindrical object. The net torque acting on the cylindrical object; S3.3 The equation of motion for the center of mass of a cylindrical object is: In the ground coordinate system, the motion components of the cylindrical object are: ； Among them, v x v y v z Let be the velocity of the cylindrical object relative to the ground coordinate system in the x, y, and z directions. By integrating, the position (x, y, z) of the object's center of mass relative to the ground coordinate system can be obtained. S3.4 To determine the attitude angle of the cylindrical object relative to the Earth, the equation of motion of the cylindrical object about its center of mass is established using the quaternion method: ； Use the following formula from the equation of motion about the center of mass R o Calculate the quaternion q: ； ； ； ； ； in, The equation of motion R about the center of mass o In the quaternion, w is the real part of the quaternion, corresponding to the scaling factor of the rotation. When describing the rotation, w is related to the angle of rotation and the direction cosine of the rotation axis. x is the first imaginary part of the quaternion, related to the x-component of the rotation axis. y is the second imaginary part of the quaternion, related to the y-component of the rotation axis. z is the third imaginary part of the quaternion, related to the z-component of the rotation axis. The quaternion differential equation is obtained as follows: ； In the formula , , This represents the component of the rigid body's angular velocity in the inertial coordinate system.

3. The underwater project trajectory prediction method based on engineering algorithms according to claim 1, characterized in that, Step S4: Based on empirical formulas, simplify the process of the cylindrical object entering the water and establish a force model of the object in the water section, as follows: S4.1 Formulas for calculating buoyancy and center of buoyancy of a cylindrical object during its immersion in water, in a ground coordinate system: ； ； In the formula, B is buoyancy, and X... B It is the center of gravity, L w V is the length of the cylindrical object submerged in water, and V is the volume of the cylindrical object submerged in water. Let g be the density of seawater, g be the acceleration due to gravity, and R be the radius of the cylindrical object. S4.2 The hydrodynamics of the water-immersion process adopts a semi-empirical formula as the basic model for calculating the hydrodynamics of a cylindrical object during the water-immersion process. The formula is as follows: ； ； In the formula, x is the axial distance from the top of the cylindrical object to the water surface during the water discharge process. Let be the dimensionless angular velocity of the cylindrical object, and v be the velocity of the cylindrical object relative to the water flow. Let Fz be the angle of attack of the cylindrical object, Fy be the axial force on the cylinder, M be the normal force on the cylinder, and M be the torque on the cylinder. The resultant force of the force obtained from the empirical formula and the buoyancy force is used as the force on the cylindrical object during the water entry phase.

4. The underwater project trajectory prediction method based on engineering algorithms according to claim 1, characterized in that, Step S5, based on the theoretical algorithms for gravity, buoyancy, position force, and inertial force, establishes the underwater force model of the cylindrical object. The specific steps are as follows: S5.1 The weight of the cylindrical object is: ； The direction of gravity is perpendicular to the ground, and gravity does not produce torque; S5.2 When a cylindrical object is completely submerged in water, the buoyancy will no longer change. The calculation formula is as follows: ； Where V is the volume of the cylindrical object submerged in water. Let g be the density of seawater, and g be the acceleration due to gravity. The direction of buoyancy is always vertically upward, and the torque generated by buoyancy is as follows: ； Where lx, ly, and lz are the distances between the center of buoyancy and the center of mass in the three directions; S5.3 The expression for position force under small angle of attack is as follows: ； ； in, Let v be the velocity of the seawater current, and v be the velocity of the cylindrical object relative to the water current. Let be the velocity of the cylindrical object relative to the ground, and its expression is as follows: ； Where V1, V2, and V3 are the set seawater flow velocities, and x0 to x3 are the set range of seawater flow velocities for each layer; After dimensionless transformation, the following expression is obtained: ； In the formula, These are the angle of attack, sideslip angle, and roll angle of a cylindrical object: Simplifying the above equation, we get: ； ； In the formula, for , The normal force coefficient when it is equal to zero This is the position derivative of the normal force coefficient with respect to the angle of attack; Normal force coefficient versus roll angle The positional derivative, for , The normal moment coefficient when it is equal to zero The normal moment coefficient is relative to the slip angle. The positional derivative; Normal moment coefficient with respect to roll angle The positional derivative, The calculation of the position force coefficient for cylindrical objects under medium or large angles of attack is handled specially, as follows: Normal force F of a cylindrical object NB The calculation formula is: ； In the formula, Let S be the maximum cross-sectional area of ​​the cylindrical object, A be the area of ​​the base of the cylindrical object, and S be the area of ​​the cylinder. p Let be the area of ​​the longitudinal projection plane of the cylindrical object. F is the crossflow drag proportionality factor. dn For crossflow resistance; Under medium or high angle of attack conditions: ； ； ； In the formula, L is the length of the cylindrical object, and x m The position of the center of mass of the cylindrical object, x c The position of the center of the longitudinal projection plane of the cylindrical object; F D F is the lateral force acting on the cylindrical object. L The axial force acting on it is S; the cross-sectional area of ​​the cylinder is S; and the radius of the cylinder is R. S5.4 Inertial forces include additional inertial forces and additional inertial moments, and their expressions are as follows: ； in, ； In the formula, and Let be the first derivative of the object's acceleration and angular velocity. To add inertial force, To add an inertial torque, For added mass, To add a moment of inertia, The density of the cylindrical object.

5. The underwater project trajectory prediction method based on engineering algorithms according to claim 1, characterized in that, Step S7: Solve the six-degree-of-freedom motion equations using the RK5 method to obtain the dynamic parameters of the cylindrical object, as follows: S7.1 The expressions for the system of differential equations are as follows: ； The solution formula for the RK5 method is: ； in, ； ； ； ； ； ； The error is: , where h is the iteration step size, w is the real part of the quaternion, and x, y, z are the imaginary parts of the quaternion.

Images