Underwater object throwing track prediction method based on engineering algorithm

Abstract

The invention relates to the field of aerospace and ships, and discloses an underwater object throwing track prediction method based on an engineering algorithm, which comprises the following steps: acquiring parameters required by calculation; a pure elastic contact model is adopted to establish multi-body dynamics description of cylinder object discharge; establishing a six-degree-of-freedom motion equation of the cylindrical object, simplifying the six-degree-of-freedom motion equation into a centroid motion equation and a centroid-around motion equation, and establishing the centroid-around motion equation by adopting a quaternion method; establishing a water entry section object stress model based on an empirical formula; establishing an underwater stress model of the cylindrical object based on gravity, buoyancy, position force and inertia force; and selecting a corresponding model to calculate the stress condition of the cylindrical object and substituting the stress condition into a six-degree-of-freedom motion equation to solve by using an RK5 method based on the dynamic parameter judgment stage when the cylindrical object is taken out of the cylinder, judging whether the object reaches a target depth or not, and outputting a final result if the object reaches the target depth. The complexity of solving a Navier-Stokes (N-S) equation is avoided, and the calculation speed is remarkably improved.

Description

Technical Field

[0001] The invention relates to the fields of aviation and navigation and ships, in particular to an underwater cylindrical object trajectory prediction method based on an engineering algorithm in the ship field. Background Art

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

[0003] Currently, underwater trajectory simulation and prediction mainly rely on computational fluid dynamics (CFD) numerical simulation technology. CFD technology can provide accurate underwater trajectory prediction by simulating fluid flow and object interaction. However, this technology has extremely high requirements on time and computing resources during the calculation process, which limits 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 may be more stringent. In addition, the flexibility of CFD technology in dealing with sudden problems is also limited. Therefore, developing an underwater trajectory prediction method with high computational efficiency, fast response speed and stable operation in a changing environment is of great significance to meet the needs of modern maritime military technology. Summary of the invention

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

[0006] Based on the above purpose, the present invention provides a method for predicting underwater object trajectory based on an engineering algorithm, the method comprising the following steps:

[0007] S1. Collecting parameters required for calculation, including parameters of the object throwing tube, parameters of the cylindrical object, parameters of seawater and calculation parameters;

[0008] S2. Using a pure elastic contact model, a multi-body dynamics description of the cylindrical object coming out of the tube is established according to the Newton-Euler equation, and the dynamic parameters of the cylindrical object when it comes out of the tube are calculated;

[0009] S3, establish the six-degree-of-freedom motion equation of the cylindrical object, simplify it into the center-of-mass motion equation and the motion equation around the center of mass, wherein the quaternion method is used to transform the motion equation around the center of mass;

[0010] S4. Based on the empirical formula, simplify the process of cylindrical objects entering water and establish the force model of the objects in the water entry stage;

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

[0012] S6. Determine the stage of the cylindrical object based on the dynamic parameters when it comes out of the tube. If it is the water entry stage, select the force model of the object in the water entry stage. If it is the underwater stage, select the underwater force model of the cylindrical object, calculate the force condition of the cylindrical object and bring it into the six-degree-of-freedom motion equation.

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

[0014] S8. Determine whether the object has reached the target depth based on the dynamic parameters of the cylindrical object. If the object has reached the target depth, output the final result. Otherwise, repeat steps S6 and S7 until the target depth is reached.

[0015] Preferably, step S2, using a pure elastic contact model, establishing a multi-body dynamics description of the cylindrical object coming out of the tube according to the Newton-Euler equation, and calculating the dynamic parameters of the cylindrical object when it comes out of the tube, includes the following steps:

[0016] S2.1. Use the classic Newton-Euler equation to accurately describe the motion of the center of mass of a cylindrical object and its changing posture relative to the center of mass of the cylindrical object;

[0017] S2.2. For the contact force generated by the cylindrical object and the object throwing tube during the contact process, a pure elastic contact model is used, in which the contact force and the penetration depth are nonlinearly related:

[0018] F n =Kδ n

[0019] In the formula, F n is the contact force, δ represents the amount of pressure between the barrel and the cylindrical object; K represents the stiffness coefficient; n is the force index, which is determined by the material of the contact body and the geometric characteristics of the contact area; the calculation formula of the stiffness coefficient K is as follows:

[0020]

[0021] In the formula, R i , R j is the radius of curvature of the contact position between the throwing tube and the cylindrical object; σ i , σ j is the parameter related to the material of the throwing tube and the cylindrical object, and its calculation formula is as follows:

[0022]

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

[0024] S2.3. The dynamic model of the cylindrical object when it leaves the cylinder fully considers the gravity of the cylindrical object and the contact force provided by the cylinder wall constraint. The dynamic equation of the motion law of the center of mass of the cylindrical object is established based on the Newton equation as follows:

[0025]

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

[0027] In addition, in order to obtain the dynamic equation of the changing posture relative to the center of mass of the cylindrical object, the following formula is formulated based on the Euler equation:

[0028]

[0029] Wherein, J1, J2, J3 are the three principal moments of inertia of the cylindrical object in the central principal axis coordinate system; ω1, ω2, ω3 are the three components of the angular velocity vector of the cylindrical object in the object coordinate system; are the first-order derivatives of the three components respectively; M1, M2, M3 are the three components of the moment vector of the cylindrical object relative to the center of mass in the cylindrical object coordinate system; when r c When it is greater than the length L of the object throwing tube, the process of the cylindrical object coming out of the tube ends.

[0030] Preferably, in step S3, a six-degree-of-freedom motion equation of a cylindrical object is established, which is simplified into a center-of-mass motion equation and a motion equation around the center of mass, wherein the motion equation around the center of mass is established using a quaternion method, as follows:

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

[0032] For the conversion of coordinate axes, the ground coordinate system finally coincides with the object coordinate system through a finite number of rotations around the axis. The angle of rotation around the axis is set to be the attitude angle of the cylindrical object, which is marked as the roll angle γ, the yaw angle φ, and the pitch angle θ. The order of rotation transformation is set to zyx, and the transformation equation is as follows:

[0033]

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

[0035] S3.2, the six-degree-of-freedom dynamic equation of a cylindrical object is:

[0036] ma=F total

[0037]

[0038] Among them, m is the mass of the cylindrical object, a is the acceleration of the cylindrical object, and J is the moment of inertia of the cylindrical object. is the first-order derivative of the angular velocity of the cylindrical object, F total is the resultant force on the cylindrical object, M total is the resultant moment acting on the cylindrical object;

[0039] S3.3, the equation of motion of 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 is the velocity of the cylindrical object relative to the ground coordinate system in the x, y, and z directions. By integrating, we can get the position of the object's center of mass relative to the ground coordinate system (x, y, z);

[0043] S3.4. To confirm the motion attitude angle of the cylindrical object relative to the earth, the motion equation of the cylindrical object around the center of mass is established by the quaternion method:

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

[0045] From the equation of motion about the center of mass R, use the following formula o Calculate the quaternion q:

[0046]

[0047] Among them, r ij is the equation of motion about the center of mass R oThe elements in , 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 value of the rotation axis. x is the first imaginary component of the quaternion, which is related to the x component of the rotation axis. y is the second imaginary component of the quaternion, which is related to the y component of the rotation axis. z is the third imaginary component of the quaternion, which is related to the z component of the rotation axis.

[0048] The quaternion differential equation is as follows:

[0049]

[0050] Where ω x ,ω y ,ω z is the component of the rigid body's angular velocity in the inertial coordinate system.

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

[0052] S4.1. The calculation formula for the buoyancy and center of buoyancy of a cylindrical object in the ground coordinate system during its landing is:

[0053] B=ρgV

[0054]

[0055] Where B is the buoyancy, X B It's floating heart, L w is the length of the cylindrical object immersed in water, V is the volume of the cylindrical object immersed 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 semi-empirical formula is used as the basic model for calculating the fluid dynamics of a cylindrical object during the water landing process. The formula is as follows:

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

[0058]

[0059] Where x is the axial distance between the top of the cylindrical object and the water surface during the water discharge process, is the dimensionless quantity of the angular velocity of the cylindrical object, v is the velocity of the cylindrical object relative to the water flow, α is the angle of attack of the cylindrical object, Fz is the axial force on the cylinder, Fy is the normal force on the cylinder, and M is the moment on the cylinder;

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

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

[0062] S5.1. The gravity of a 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 the 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 immersed 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] Among them, lx, ly, and lz are the distances between the center of buoyancy and the center of mass in three directions;

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

[0072]

[0073] v=v i -V sea

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

[0075]

[0076] Among them, V1, V2, V3 are the set seawater flow rates, and x0 to x3 are the range of the set seawater flow rate for each layer;

[0077] After non-dimensionalization, the following expression is obtained:

[0078]

[0079] Where α, β, and δ are the angle of attack, sideslip angle, and roll angle of the cylindrical object, respectively:

[0080] Simplifying the above formula, we get:

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

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

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

[0084] The calculation of the position force coefficient for cylindrical objects at medium or large angles of attack uses special treatment, as follows:

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

[0086]

[0087] In the formula, S max is the maximum cross-sectional area of ​​the cylindrical object, A is the bottom area of ​​the cylindrical object, S p is the area of ​​the longitudinal projection plane of the cylindrical object, η is the cross-flow resistance proportional factor, F dn is the cross flow resistance;

[0088] At medium or high angles of attack:

[0089] F D =F NB sinα

[0090] F L =F NB cosα

[0091]

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

[0093] S5.4. Inertia force includes additional inertia force and additional inertia moment, which are expressed as follows:

[0094]

[0095] in,

[0096]

[0097] In the formula, a and is the first-order derivative of the object’s acceleration and angular velocity, F added is the additional inertia force, M added is the additional inertia moment, m added is the additional mass, I added is the additional moment of inertia, ρ object is 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, which are specifically as follows:

[0099] S7.1, the differential equation system is expressed as follows:

[0100]

[0101] The solution formula of 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 the present invention are as follows: (1) The present invention proposes a method for predicting underwater object trajectory based on an engineering algorithm, which simplifies the calculation process and avoids the complexity of solving the Navier-Stokes (NS) equation, thereby significantly improving the calculation speed. While maintaining reasonable prediction accuracy, the method of the present invention can quickly adapt to different underwater environments and conditions, effectively solving the limitations of the prior art.

[0110] (2) The present invention is based on an engineering algorithm and does not need to solve the NS equations, thereby greatly reducing the computing cost for 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 the object.

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

[0113] (5) The role of the present invention in underwater motion trajectory prediction is in line with the needs of various countries in this field and has great strategic value in the field of ocean exploration. BRIEF DESCRIPTION OF THE DRAWINGS

[0114] Figure 1 Analyze the process for software.

[0115] Figure 2 This is a graph showing the variation of the velocity in the z direction of the object in the vertical stationary cylinder in Example 1 with time (the velocities in the other two directions are 0).

[0116] Figure 3 For comparison between Example 1 and the prior art, a graph showing the variation of the y-direction displacement of an object falling in water with the falling depth is shown.

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

[0118] Figure 5 This is a graph showing how the velocity of the object in the z and x directions changes with time in Example 2 (the velocities in other directions are 0).

[0119] Figure 6 For comparison between Example 2 and the prior art, a graph showing the variation of the y-direction displacement of an object falling in water with the falling depth is shown.

[0120] Figure 7 This is a three-dimensional displacement diagram of the object in Example 2 when falling in water. DETAILED DESCRIPTION

[0121] In order to make the purpose, technical solution and advantages of the embodiments of the present application clearer, the technical solution in the embodiments of the present application will be clearly and completely described below in conjunction with the drawings in the embodiments of the present application. Obviously, the described embodiments are part of the embodiments of the present application, not all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of this application.

[0122] It should be noted that, unless otherwise defined, the technical terms or scientific terms used in the embodiments of the present invention should be understood by people with ordinary skills in the field to which the present invention belongs. The "first", "second" and similar words used in the present invention do not represent any order, quantity or importance, but are only used to distinguish different components. "Including" or "comprising" and similar words mean that the elements or objects appearing before the word cover the elements or objects listed after the word and their equivalents, without excluding other elements or objects. "Connecting" or "connected" and similar words are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. "Up", "down", "left", "right" and the like are only used 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 are further described in detail below with reference to the accompanying drawings.

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

[0125] In this embodiment, based on this algorithm, the "Cylindrical Object Underwater Motion Trajectory Calculation Software" was developed on the Visual Studio platform using the C++ language. The software analysis process is as follows: Figure 1 shown.

[0126] First, you need to input the parameters of the throwing tube and the cylindrical object to calculate the motion state of the cylindrical object when it leaves the throwing tube, such as Figure 2 and Figure 3 As shown in the figure, it serves as the initial state of the hydrodynamic analysis in the following steps; then it is necessary to input seawater parameters and calculation parameters, and draw the following Figure 4 The 3D displacement map shown, in this example the target depth is 60m.

[0127] S2. Using a pure elastic contact model, a multi-body dynamics description of the cylindrical object coming out of the tube is established according to the Newton-Euler equation, and the dynamic parameters of the cylindrical object when it comes out of the tube are calculated;

[0128] In S2, a dynamic model of a cylindrical object in a throwing tube is constructed to solve the dynamic parameters of the cylindrical object when it leaves the tube, including the following steps:

[0129] S2.1. Use the classic Newton-Euler equation to accurately describe the motion of the center of mass of a cylindrical object and its changing posture relative to the center of mass of the cylindrical object;

[0130] S2.2. For the contact force generated by the cylindrical object and the object throwing tube during the contact process, a pure elastic contact model is used, in which the contact force and the penetration depth have a nonlinear relationship, and the nonlinear relationship is as follows:

[0131] F n =Kδ n

[0132] In the formula, F n is the contact force, δ represents the amount of pressure between the barrel and the cylindrical object; K represents the stiffness coefficient; n is the force index, which defaults to 1.5 and is determined by the material of the contact body and the geometric characteristics of the contact area; the calculation formula of the stiffness coefficient K is as follows:

[0133]

[0134] In the formula, R i , R j is the radius of curvature of the contact position between the throwing tube and the cylindrical object; σ i , σ j is the parameter related to the material of the throwing tube and the cylindrical object, and its calculation formula is as follows:

[0135]

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

[0137] S2.3. The dynamic model of the cylindrical object when it leaves the cylinder fully considers the gravity of the cylindrical object and the contact force provided by the cylinder wall constraint. The dynamic equation of the motion law of the center of mass of the cylindrical object is established based on the Newton equation as follows:

[0138]

[0139] Where m represents the mass of the cylindrical object; is the second-order derivative of the axial displacement of the center of mass of the cylindrical object in the global coordinate system; ω0 is the rotation speed of the cylinder; t is the current physical time; F is the resultant force of the contact force and gravity on the object in the axial direction;

[0140] In addition, in order to obtain the dynamic equation of the changing posture relative to the center of mass of the cylindrical object, the following formula is formulated based on the Euler equation:

[0141]

[0142] Wherein, J1, J2, J3 are the three principal moments of inertia of the cylindrical object in the central principal axis coordinate system; ω1, ω2, ω3 are the three components of the angular velocity vector of the cylindrical object in the object coordinate system; are the first-order derivatives of the three components; M1, M2, and M3 are the three components of the moment vector of the cylindrical object relative to the center of mass in the cylindrical object coordinate system; when r c When it is greater than the length L of the object throwing tube, the object out of the tube process ends.

[0143] S3. Establish the six-degree-of-freedom motion equation of the cylindrical object, including the center-of-mass motion equation and the motion equation around the center of mass, wherein the quaternion method is used to establish the motion equation around the center of mass, as follows:

[0144] S3.1. First, define the coordinate system to determine the relative position, velocity, and direction of the 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 takes a point on the ground (usually the launch position of the cylindrical object) as its origin, the oz axis is in the plumb plane, and the positive direction is defined as upward; the oy axis is in the forward direction of the ship; and the ox axis is determined according to the right-hand rule.

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

[0147] For the conversion of coordinate axes, the ground coordinate system finally coincides with the object coordinate system through a finite number of rotations around the axis. The angle of rotation around the axis is set to be the attitude angle of the cylindrical object, which is marked as the roll angle γ, the yaw angle φ, and the pitch angle θ. The order of rotation transformation is set to ZYX, and the transformation equation is as follows:

[0148]

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

[0150] S3.2, the six-degree-of-freedom dynamic equation of a cylindrical object is:

[0151] ma=F total

[0152]

[0153] Among them, m is the mass of the cylindrical object, a is the acceleration of the cylindrical object, and J is the moment of inertia of the cylindrical object. is the first-order derivative of the angular velocity of the cylindrical object, F total is the resultant force on the cylindrical object, M total is the resultant moment acting on the cylindrical object;

[0154] S3.3, the equation of motion of 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 is the velocity of the cylindrical object relative to the ground coordinate system in the x, y, and z directions. By integrating, we can get the position of the object's center of mass relative to the ground coordinate system (x, y, z);

[0158] S3.4. To confirm the motion attitude angle of the cylindrical object relative to the earth, the motion equation of the cylindrical object around the center of mass is established by the quaternion method:

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

[0160] From the equation of motion about the center of mass R, use the following formula o Calculate the quaternion q:

[0161]

[0162] Among them, r ij is the equation of motion about the center of mass R o The elements in , 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 value of the rotation axis. x is the first imaginary component of the quaternion, which is related to the x component of the rotation axis. y is the second imaginary component of the quaternion, which is related to the y component of the rotation axis. z is the third imaginary component of the quaternion, which is related to the z component of the rotation axis.

[0163] The quaternion differential equation is as follows:

[0164]

[0165] Where ω x ,ω y ,ω z is the component of the rigid body's angular velocity in the inertial coordinate system

[0166] S4. Based on the empirical formula, the process of cylindrical objects entering water is simplified and the force model of the object entering water is established, specifically:

[0167] S4.1. The calculation formula for the buoyancy and center of buoyancy of a cylindrical object in the ground coordinate system during its landing is:

[0168] B=ρgV

[0169]

[0170] Where B is the buoyancy, X B It's floating heart, L w is the length of the cylindrical object immersed in water, V is the volume of the cylindrical object immersed 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 semi-empirical formula is used as the basic model for calculating the fluid dynamics of a cylindrical object during the water landing process. The formula is as follows:

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

[0173]

[0174] Where x is the axial distance between the top of the cylindrical object and the water surface during the process of exiting the water, R is the radius of the cylindrical object, is the dimensionless quantity of the angular velocity of the cylindrical object, v is the velocity of the cylindrical object relative to the water flow, Fz is the axial force on the cylinder, Fy is the normal force on the cylinder, and M is the moment on the cylinder;

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

[0176] S5. Divide the underwater force of the cylindrical object into theoretical algorithms of gravity, buoyancy, position force, and inertia force, establish an underwater force model of the cylindrical object, and solve each part to obtain the underwater force distribution of the cylindrical object. The specific steps are as follows:

[0177] S5.1. The gravity of a 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 the 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] Among them, lx, ly, and lz are the distances between the center of buoyancy and the center of mass in three directions;

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

[0186]

[0187] v=v i -V sea

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

[0189]

[0190] Among them, V1, V2, V3 are the set seawater flow rates, and x0 to x3 are the range of the set seawater flow rate for each layer;

[0191] After non-dimensionalization, the following expression is obtained:

[0192]

[0193] Where α, β, and δ are the angle of attack, sideslip angle, and roll angle of the cylindrical object, respectively:

[0194] Simplifying the above formula, we get:

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

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

[0197] Where C0 is the normal force coefficient when α and δ are equal to zero, C α is the position derivative of the normal force coefficient with respect to the angle of attack; C δ is the position 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 β M is the position derivative of the normal moment coefficient with respect to the sideslip angle β; δ is the position derivative of the normal moment coefficient with respect to the rolling angle δ;

[0198] The calculation of the position force coefficient for cylindrical objects at medium or large angles of attack uses special treatment, as follows:

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

[0200]

[0201] In the formula, S max is the maximum cross-sectional area of ​​the cylindrical object, A is the bottom area of ​​the cylindrical object, S p is the area of ​​the longitudinal projection plane of the cylindrical object, η is the cross-flow resistance proportional factor, F dn is the cross flow resistance;

[0202] At medium or high angles of attack:

[0203] F D =F NB sinα

[0204] F L =F NB cosα

[0205]

[0206] Where L is the length of the cylindrical object, x m The location of the center of mass of the cylindrical object, x c The position of the center point of the longitudinal projection plane of a cylindrical object.

[0207] S5.4. Inertia force includes additional inertia force and additional inertia moment, which are expressed as follows:

[0208]

[0209] in,

[0210]

[0211] In the formula, a and is the first-order derivative of the object’s acceleration and angular velocity, F added is the additional inertia force, M added is the additional inertia moment, m added is the additional mass, I added is the additional moment of inertia, ρ object is the density of the cylindrical object.

[0212] S7, based on the RK5 method, solve the quaternion differential equation to obtain the rotation parameters of the cylindrical object, as follows:

[0213] S7.1, the differential equation system is expressed as follows:

[0214]

[0215] The solution formula of 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 it has reached the target depth, output the final result. Otherwise, repeat processes 6 and 7.

[0224] Example 1

[0225] Select the case where the throwing tube is stationary and the object is thrown vertically downward. The length of the throwing tube is 2m and the mass of the object is 133kg. The distribution of seawater flow velocity is: between 0-10m water depth, the seawater velocity is 0; between 10-20m water depth, the seawater velocity is 10m / s, and the direction is along the positive and negative direction of the y-axis; between 20m water depth and the bottom of the water, the seawater velocity is 30m / s, and the direction is along the positive direction of the y-axis. The falling velocity change of the object out of the tube is obtained as follows Figure 2 (The speeds in the other two directions are 0). The specific input parameters are shown in Table 1. 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 parameters of Example 1

[0227] Parameter name expression value Cylindrical object material parameters <![CDATA[σ j ]]> 200Gpa Throwing tube speed <![CDATA[ω0]]> 0rad / s Length of the tube <![CDATA[L i ]]> 2m Mass of cylindrical object m 133kg The first layer of seawater depth x1 -10m The first layer of seawater speed <![CDATA[V1]]> 0m / s The second layer of seawater depth x2 -20m Second layer seawater speed <![CDATA[V2]]> 10m / s The third layer of seawater depth x3 -60m The third layer of seawater speed <![CDATA[V3]]> 30m / s Iteration time step h 0.001s

[0228] Table 2 Results comparison table

[0229] Commercial CFD Methods Prediction method of the present invention Difference Time consumption 3.82h 0.25s 3.82h Final displacement calculation results 66.72m 64.12m 2.60m

[0230] Example 2

[0231] The throwing tube is selected to perform regular swinging motion, with a swing angle of 30°, a direction along the positive direction of the x-axis, a swing frequency of 0.1 rad / s, a length of 2m, and an object mass of 133kg. The seawater velocity distribution is as follows: between 0-10m in water depth, the seawater velocity is 0; between 10-20m in water depth, the seawater velocity is 10m / s, and the direction is along the positive and negative directions of the y-axis; between 20m in water depth and the bottom of the water, the seawater velocity is 30m / s, and the direction is in the opposite direction of the y-axis. The input parameters are shown in Table 3 below, and the falling velocity change and the y-direction velocity change of the object in the process of leaving the tube are obtained as shown in Table 3 below. Figure 5 (the speed 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 parameters of Example 2

[0233] Parameter name expression value Cylindrical object material parameters <![CDATA[σ j ]]> 200Gpa Throwing tube speed <![CDATA[ω0]]> 0.1rad / s Length of the tube <![CDATA[L i ]]> 2m Mass of cylindrical object m 133kg The first layer of seawater depth x1 -10m The first layer of seawater speed <![CDATA[V1]]> 0m / s The second layer of seawater depth x2 -20m Second layer seawater speed <![CDATA[V2]]> 10m / s The third layer of seawater depth x3 -60m The third layer of seawater speed <![CDATA[V3]]> -30m / s Iteration time step h 0.001s

[0234] Table 4 Results comparison table

[0235] Commercial CFD Methods Prediction method of the present invention Difference Time consumption 4.62h 0.28s 4.62h Final displacement calculation results -66.01m -59.77m 6.24m

Claims

1. A method for predicting underwater object trajectory based on engineering algorithm, characterized in that: The method comprises the following steps: S1. Collecting parameters required for calculation, including parameters of the object throwing tube, parameters of the cylindrical object, parameters of seawater and calculation parameters; S2. Using a pure elastic contact model, a multi-body dynamics description of the cylindrical object coming out of the tube is established according to the Newton-Euler equation, and the dynamic parameters of the cylindrical object when it comes out of the tube are calculated; S3, establish the six-degree-of-freedom motion equation of the cylindrical object, simplify it into the center-of-mass motion equation and the motion equation around the center of mass, wherein the quaternion method is used to transform the motion equation around the center of mass; S4. Based on the empirical formula, simplify the process of cylindrical objects entering water and establish the force model of the objects in the water entry stage; S5. Based on gravity, buoyancy, position force and inertia force, establish an underwater force model for cylindrical objects; S6. Determine the stage of the cylindrical object based on the dynamic parameters when it comes out of the tube. If it is the water entry stage, select the force model of the object in the water entry stage. If it is the underwater stage, select the underwater force model of the cylindrical object, calculate the force condition of the cylindrical object and bring it into the six-degree-of-freedom motion equation. S7. Solve the six-degree-of-freedom motion equation based on the RK5 method to obtain the dynamic parameters of the cylindrical object; S8. Determine whether the object has reached the target depth based on the dynamic parameters of the cylindrical object. If the object has reached the target depth, output the final result. Otherwise, repeat steps S6 and S7 until the target depth is reached.

2. The underwater object trajectory prediction method based on engineering algorithm according to claim 1 is characterized in that: Step S2, using a pure elastic contact model, establishing a multi-body dynamic description of the cylindrical object coming out of the tube according to the Newton-Euler equation, and calculating the dynamic parameters of the cylindrical object when it comes out of the tube, including the following steps: S2.

1. Use the classic Newton-Euler equation to accurately describe the motion of the center of mass of a cylindrical object and its changing posture relative to the center of mass of the cylindrical object; S2.

2. For the contact force generated by the cylindrical object and the object throwing tube during the contact process, a pure elastic contact model is used, in which the contact force and the penetration depth are nonlinearly related: F n =Kδ n In the formula, F n is the contact force, δ represents the amount of pressure between the barrel and the cylindrical object; K represents the stiffness coefficient; n is the force index, which is determined by the material of the contact body and the geometric characteristics of the contact area; the calculation formula of the stiffness coefficient K is as follows: In the formula, R i , R j is the radius of curvature of the contact position between the throwing tube and the cylindrical object; σ i , σ j is the parameter related to the material of the throwing tube and the cylindrical object, and its calculation formula is as follows: In the formula, E I ,υ I are the Poisson’s ratio and Young’s modulus of the material; S2.

3. The dynamic model of the cylindrical object when it leaves the cylinder fully considers the gravity of the cylindrical object and the contact force provided by the cylinder wall constraint. The dynamic equation of the motion law of the center of mass of the cylindrical object is established based on the Newton equation as follows: Where m represents the mass of the cylindrical object; is the second-order derivative of the axial displacement of the center of mass of the cylindrical object in the global coordinate system; ω0 is the rotation speed of the cylinder; t is the current physical time; F is the resultant force of the contact force and gravity on the object in the axial direction; In addition, in order to obtain the dynamic equation of the changing posture relative to the center of mass of the cylindrical object, the following formula is formulated based on the Euler equation: Wherein, J1, J2, J3 are the three principal moments of inertia of the cylindrical object in the central principal axis coordinate system; ω1, ω2, ω3 are the three components of the angular velocity vector of the cylindrical object in the object coordinate system; are the first-order derivatives of the three components respectively; M1, M2, M3 are the three components of the moment vector of the cylindrical object relative to the center of mass in the cylindrical object coordinate system; when r c Greater than the length of the tube L i When , the process of cylindrical object coming out of the tube is completed.

3. The underwater object trajectory prediction method based on engineering algorithm according to claim 1 is characterized in that: Step S3, establish the six-degree-of-freedom motion equation of the cylindrical object, simplify it into the center-of-mass motion equation and the motion equation around the center of mass, wherein the quaternion method is used to establish the motion equation around the center of mass, as follows: S3.

1. First, define the coordinate system to determine the relative position, velocity, and direction of the external force vector parameters. The coordinate system includes the ground coordinate system oxyz and the object coordinate system o1x1y1z1. For the conversion of coordinate axes, the ground coordinate system finally coincides with the object coordinate system through a finite number of rotations around the axis. The angle of rotation around the axis is set to be the attitude angle of the cylindrical object, which is marked as the roll angle γ, the yaw angle φ, and the pitch angle θ. The order of rotation transformation is set to zyx, and the transformation equation is as follows: Among them, R z (φ), R y (γ), R x (θ) are the transformation matrices for rotation around the z, y, and x axes respectively; S3.2, the six-degree-of-freedom dynamic equation of a cylindrical object is: in=F total Among them, m is the mass of the cylindrical object, a is the acceleration of the cylindrical object, and J is the moment of inertia of the cylindrical object. is the first-order derivative of the angular velocity of the cylindrical object, F total is the resultant force on the cylindrical object, M total is the resultant moment acting on the cylindrical object; S3.3, the equation of motion of 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 is the velocity of the cylindrical object relative to the ground coordinate system in the x, y, and z directions. By integrating, we can get the position of the object's center of mass relative to the ground coordinate system (x, y, z); S3.

4. To confirm the motion attitude angle of the cylindrical object relative to the earth, the motion equation of the cylindrical object around the center of mass is established by the quaternion method: R o =R z (φ)R y (c)R x (i) From the equation of motion about the center of mass R, use the following formula o Calculate the quaternion q: Among them, r ij is the equation of motion about the center of mass R o The elements in , 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 value of the rotation axis. x is the first imaginary component of the quaternion, which is related to the x component of the rotation axis. y is the second imaginary component of the quaternion, which is related to the y component of the rotation axis. z is the third imaginary component of the quaternion, which is related to the z component of the rotation axis. The quaternion differential equation is as follows: Where ω x ,ω y ,ω z is the component of the rigid body's angular velocity in the inertial coordinate system.

4. The underwater object trajectory prediction method based on engineering algorithm according to claim 1 is characterized in that: Step S4, based on the empirical formula, simplify the process of the cylindrical object entering the water and establish the force model of the object in the water entry section, which is as follows: S4.

1. The calculation formula for the buoyancy and center of buoyancy of a cylindrical object in the ground coordinate system during its landing is: B=ρgV Where B is the buoyancy, X B It's floating heart, L w is the length of the cylindrical object immersed in water, V is the volume of the cylindrical object immersed in water, ρ is the density of seawater, g is the acceleration due to gravity, and R is the radius of the cylindrical object; S4.

2. The semi-empirical formula is used as the basic model for calculating the fluid dynamics of a cylindrical object during the water landing process. The formula is as follows: Fz=-0.0076ρπR 2 v 2 (10-x) Where x is the axial distance between the top of the cylindrical object and the water surface during the water discharge process, is the dimensionless quantity of the angular velocity of the cylindrical object, v is the velocity of the cylindrical object relative to the water flow, α is the angle of attack of the cylindrical object, Fz is the axial force on the cylinder, Fy is the normal force on the cylinder, and M is the moment on the cylinder; The resultant force of the force obtained from the empirical formula and the buoyancy is used as the force acting on the cylindrical object during the entry into water stage.

5. The underwater object trajectory prediction method based on engineering algorithm according to claim 1 is characterized in that: Step S5, based on the theoretical algorithms of gravity, buoyancy, position force and inertia force, establish the underwater force model of the cylindrical object. The specific steps are as follows: S5.

1. The gravity of a cylindrical object is: G=mg The direction of gravity is perpendicular to the ground, and gravity does not produce torque; S5.

2. When the cylindrical object is completely submerged in water, the buoyancy will no longer change. The calculation formula is as follows: B=ρgV Where V is the volume of the cylindrical object immersed in water, ρ is the density of seawater, and g is the acceleration due to gravity; The direction of buoyancy is always vertically upward, and the torque generated by buoyancy is as follows: Among them, lx, ly, and lz are the distances between the center of buoyancy and the center of mass in three directions; S5.

3. The expression of position force at small angle of attack is as follows: v=v i -V sea Among them, V sea is the velocity of seawater, v is the velocity of the cylindrical object relative to the water flow, v i is the speed of the cylindrical object relative to the ground, and its expression is as follows: Among them, V1, V2, V3 are the set seawater flow rates, and x0 to x3 are the range of the set seawater flow rate for each layer; After non-dimensionalization, the following expression is obtained: Where α, β, and δ are the angle of attack, sideslip angle, and roll angle of the cylindrical object, respectively: Simplifying the above formula, we get: C l =C0+C α α+C δ δ M l =M0+M β β+M δ δ Where C0 is the normal force coefficient when α and δ are equal to zero, C α is the position derivative of the normal force coefficient with respect to the angle of attack; C δ is the position derivative of the normal force coefficient with respect to the roll angle δ, M0 is the normal moment coefficient when β and δ are equal to zero, and M β M is the position derivative of the normal moment coefficient with respect to the sideslip angle β; δ is the position derivative of the normal moment coefficient with respect to the roll angle δ, The calculation of the position force coefficient for cylindrical objects at medium or large angles of attack uses special treatment, as follows: Normal force F of a cylindrical object NB The calculation formula is: In the formula, S max is the maximum cross-sectional area of ​​the cylindrical object, A is the bottom area of ​​the cylindrical object, S p is the area of ​​the longitudinal projection plane of the cylindrical object, η is the cross-flow resistance proportional factor, F dn is the cross flow resistance; At medium or high angles of attack: F D =F NB sinα F L =F NB cosα Where L is the length of the cylindrical object, x m The location of the center of mass of the cylindrical object, x c The position of the center of the plane of longitudinal projection of the cylindrical object; S5.

4. Inertia force includes additional inertia force and additional inertia moment, which are expressed as follows: in, In the formula, a and is the first-order derivative of the object’s acceleration and angular velocity, F added is the additional inertia force, M added is the additional inertia moment, m added is the additional mass, I added is the additional moment of inertia, ρ object is the density of the cylindrical object.

6. The underwater object trajectory prediction method based on engineering algorithm according to claim 1 is characterized in that: Step S7, solving the six-degree-of-freedom motion equation based on the RK5 method to obtain the dynamic parameters of the cylindrical object, which are as follows: S7.1, the differential equation system is expressed as follows: The solution formula of the RK5 method is: in, k1=h·f(t n ,y n ) k5=h·f(t n +h,y n -5k1+3k4) The error is: Where h is the iteration step size.

Images