A method for solving the nonlinear dynamic response of a towing system

Abstract

The present invention discloses a method for solving the nonlinear dynamic response of a towing system, including the following steps: initializing the time and the number of calculation steps, with the initialization value being 0; initializing the speed and displacement values, with the initialization value being a given initial value; using the speed and displacement value of time step i as the speed and displacement value of time step i+1 to initialize the next time step; constructing a towing system unit matrix based on the node coordinate finite element method; constructing a conversion matrix from local coordinates to global coordinates based on Euler angle transformation, and determining the towing system unit matrix in the global coordinate system; determining the external force matrix by considering external forces, and constructing the towing system dynamic control equation based on the towing system unit matrix in the global coordinate system; solving the dynamic control equation response by using a symplectic difference algorithm, and judging whether the solution meets the conditions. The present invention effectively reduces the cumulative error caused by long-term calculation and analysis, and improves the solution accuracy.

Description

Technical Field

[0001] The invention relates to a flexible cable towing system, and in particular to a method for solving the nonlinear dynamic response of the towing system. Background Art

[0002] In engineering structures, the space flexible cable structure is a typical component. From a physical point of view, its longitudinal dimension is much larger than its transverse dimension, showing great flexibility. From a mechanical point of view, the flexible cable bears tensile loads, and its ability to resist compression, bending, shearing and torsion is extremely low. The flexible cable structure has the characteristics of less material consumption, high strength and low flexibility, and is very easy to connect two systems at a long distance. As a medium for signal transmission and structural traction, it has been widely used in space tethered satellite systems, aerial refueling systems, drones, micro-aircraft systems, unmanned remote-controlled submersible systems, marine engineering, and large-span suspension bridge structures. Its application scope continues to expand and is closely related to modern science and technology. The dynamics of the flexible cable system is significantly affected by the working environment, showing strong nonlinear characteristics and complex dynamic response behavior. The application of theoretical analytical methods is limited to very simple working conditions and quasi-static situations. In most cases, there is no complete theoretical solution system that can integrate the medium resistance effect and the cable system movement, and the problem is mainly solved by numerical methods. The study of cable dynamic behavior has always been a hot issue. Among the methods for studying cable dynamics, the finite element method is based on the classical virtual work principle. The cable microelement is derived from the virtual work equilibrium equation to construct the dynamic control equation of the cable system, which is a set of second-order ordinary differential equations. For different types of structural units, the virtual work principle equilibrium equation of the corresponding form is constructed by selecting the interpolation shape function of the appropriate order, and the numerical algorithm is used to integrate and solve it. The finite element method has the advantages of dealing with complex boundary conditions and different cable element material properties, and is widely used in the analysis of flexible cable towing systems. However, the solution of this method cannot solve the energy dissipation problem in long-term integration. Traditional numerical integration methods (such as Runge-Kutta method, Newmark method, etc.) can effectively predict the dynamic response of cable structures. However, for flexible cable structures undergoing long-term and large displacement motions, the traditional integration algorithm gradually becomes numerically unstable due to the long-term coupling of low-frequency rigid body motion and high-frequency elastic deformation.

[0003] In 2019, Kan Ziyun et al. published "A method for analyzing the dynamic response of an aggregated tensegrity structure based on sliding rope units of a multi-body system". Based on the traditional rope units of the multi-body system, a sliding rope unit of a multi-body system was established to solve the multi-body dynamics differential algebraic equations to obtain the dynamic response of the aggregated tensegrity structure. This method cannot overcome the nonlinear problems in solving large displacement dynamic control equations. In the article "Characteristics of coupled orbital-attitude dynamics of flexible electric solar wind sail", G Li et al. used the finite element method to study the dynamic response of the cable system for the rope dragging problem of the solar sail, but the solution results could not solve the energy dissipation problem in the long-time integral. Summary of the invention

[0004] The purpose of the present invention is to provide a method for solving the nonlinear dynamic response of a towing system, which can provide a detailed and accurate solution to the dynamic performance of a flexible cable system in a fluid medium, effectively reduce the cumulative error caused by long-term calculation and analysis, and improve the solution accuracy.

[0005] The technical solution to achieve the purpose of the present invention is: a method for solving the nonlinear dynamic response of a towing system, comprising the steps of:

[0006] Time and calculation steps are initialized to 0; speed and displacement values ​​are initialized to the given initial values;

[0007] The velocity displacement value of time step i is used as the velocity displacement value of time step i+1 to initialize the next time step;

[0008] Construct the towing system unit array based on the nodal coordinate finite element method;

[0009] Based on the Euler angle transformation, a transformation matrix from local coordinates to global coordinates is constructed to determine the unit matrix of the towing system in the global coordinate system.

[0010] Considering the external force to determine the external force matrix, the dynamic control equation of the towing system is constructed based on the towing system unit matrix in the global coordinate system;

[0011] The response of the dynamic control equations is solved by the symplectic difference algorithm;

[0012] If the convergence of the solved displacement velocity does not meet the requirements, the displacement velocity matrix is ​​updated and the symplectic difference algorithm is used to continue to solve the dynamic control equation response; otherwise, it is determined whether the calculation time reaches the set total calculation time. If the total calculation time is reached, the result is output and the solution is completed; if the total calculation time is not reached, the next time step is entered for a new round of calculation.

[0013] Furthermore, the construction of the towing system unit array based on the node coordinate finite element method specifically includes:

[0014] Construct the element interpolation shape function N based on the nodal coordinate finite element method;

[0015] The node coordinates are used as variables to represent the energy of the unit and construct the unit matrix.

[0016] Furthermore, the unit matrix includes:

[0017] Unit mass matrix M e :

[0018]

[0019] Where L is the unit length, ρ is the material density of the unit, and A is the cross-sectional area of ​​the unit;

[0020] Element stiffness matrix K e :

[0021]

[0022] Where L0 is the initial length of the unit, θ X ,θ Y ,θ Z are the angles between the local coordinate system and the global coordinate system X, Y, and Z axes, namely the Euler angles, and E is the unit Young's modulus;

[0023] Unit displacement matrix R:

[0024] R=NX e

[0025] Where R is the position vector of any point on the unit, X e are the position vectors of the two nodes of the unit;

[0026] Unit damping array C e :

[0027] C e =αM e +βK e

[0028] Among them, α and β are damping coefficients.

[0029] Furthermore, the transformation matrix from local coordinates to global coordinates is constructed based on the Euler angle transformation:

[0030] Through the three Euler angles θ x -Roll, theta y - Pitch, θ z - Yaw describes the basis vectors of the local coordinate system (ex ,e y ,e z ); the order of transformation from the global coordinate system to the local coordinate system is defined as yaw, pitch and roll respectively; the transformation is divided into three steps:

[0031] (1) Rotate the OZ axis counterclockwise at an angle of θ z , 0≤θ z <π, so that the OY axis coincides with the ON axis of the local coordinate system;

[0032] (2) Rotate the new OY axis clockwise by an angle of θ y , 0≤θ y <π, so that the OX axis in the global coordinate system coincides with the ox axis in the local coordinate system;

[0033] (3) Rotate the new OX axis counterclockwise by an angle of θ z , 0≤θ z <π, so that the OY axis in the global coordinate system coincides with the oy axis in the local coordinate system, and at the same time, the OZ axis coincides with the oz axis in the local coordinate system;

[0034] Through the above rotation, the transformation matrix T from the global coordinate system to the local coordinate system is obtained.

[0035] Furthermore, the external force matrix F is:

[0036] F=F Ge +F de

[0037]

[0038]

[0039]

[0040] Among them, F Ge is the unit gravity vector in the global coordinate system, F de is the unit resistance in the global coordinate system, C d is the drag coefficient, ρ f is the medium density, L is the unit length, N G is the linear interpolation function, V f is the medium flow velocity of the cable structure at that point.

[0041] Furthermore, the linear interpolation function is:

[0042]

[0043] where ξ is the local coordinate of the element.

[0044] Furthermore, the construction of the towing system dynamics control equation includes:

[0045] Establish the Hamiltonian canonical equations of the towing system unit;

[0046] Based on the continuity conditions of node displacement and velocity, the unit array in the overall coordinate system is assembled into the dynamic control equation of the overall towing system.

[0047] Furthermore, the Hamiltonian canonical equation of the towing system unit is:

[0048]

[0049] Among them, the conjugate momentum of the unit The velocity vector of the cell is The kinetic energy of the unit is He is the Hamiltonian function of the unit, F ke and F Ge They are the internal force matrix vectors and gravity vector of the two nodes of the unit respectively.

[0050] Furthermore, the dynamic control equation of the towing system is:

[0051]

[0052] Among them, M is the system mass matrix, K is the system stiffness matrix, C is the system damping matrix, and F k is the system internal force matrix, F G is the system gravity matrix.

[0053] Furthermore, the symplectic difference algorithm decomposes the Hamiltonian canonical equation at the k+1 time step to obtain the matrix:

[0054]

[0055]

[0056] Here, h is the computation time step.

[0057] Compared with the prior art, the present invention has the following remarkable effects:

[0058] (1) Compared with the traditional finite element method, the nodal coordinate finite element method adopts the full-quantity format, which is simple to express and has high solution accuracy, effectively improving the cumbersome solution process of the traditional finite element method;

[0059] (2) The symplectic difference algorithm is used to solve the dynamic control equations of the cable system. Compared with the traditional integral algorithm, it avoids the problem that the cumulative calculation error gradually increases over time, which will significantly reduce the prediction accuracy of the dynamic response, and improves the solution accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0060] Figure 1 The present invention is a flow chart for the derivation of the towing system matrix and dynamic control equations based on the node coordinate finite element method and the symplectic difference algorithm.

[0061] Figure 2 The figure is a flow chart of the method for solving the nonlinear dynamic response of the towing system of the present invention. DETAILED DESCRIPTION

[0062] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0063] In order to effectively solve the dynamic response of the tethered towing system in the fluid medium, the present invention combines the node coordinate finite element method with simple expression and high solution accuracy and the symplectic difference algorithm with high accuracy in long-term simulation. The method includes the following steps:

[0064] Combination Figure 1 The method for constructing the matrix and dynamic control equation of the towing system based on the node coordinate finite element method and the symplectic difference algorithm provided by the present invention comprises the following steps:

[0065] A1. Constructing unit interpolation shape functions based on nodal coordinate finite element method

[0066] Specifically, consider a straight cable element with one or two nodes. The position of the element is described by its node coordinates. The position of any point in the cable element can be expressed by the shape function of the element in the overall frame and the coordinates of the two nodes of the element.

[0067] R=NX e

[0068] Where R is the position vector of any point on the flexible cable unit, X e is the position vector of the two nodes of the unit, and N is the interpolation shape function. Similarly, the velocity and acceleration of any point on the unit will also be interpolated by the same unit shape function and.

[0069]

[0070]

[0071] Where v and a are the velocity and acceleration vectors of any point on the cable unit, respectively. and are the velocity and acceleration vectors of the two nodes of the element, respectively.

[0072] A2. Construction unit matrix

[0073] Specifically, the energy of the unit is represented by the node coordinates as variables, including kinetic energy, elastic potential energy, and gravitational potential energy.

[0074]

[0075]

[0076]

[0077] Among them, T e and U e , U ge are the unit kinetic energy, elastic potential energy, and gravitational potential energy, respectively. L and ε are the unit length and strain, respectively. ρ and E are the unit material density and Young's modulus, respectively. A is the cross-sectional area of ​​the unit, and g is the gravitational acceleration. In the above expressions, Me and Ke are the unit mass matrix and stiffness matrix, respectively. ke and F Ge They are the internal force matrix vectors of the two nodes of the unit and the gravity vector. The mass matrix and stiffness matrix of the unit can be obtained from the mechanical expression of the unit energy.

[0078]

[0079]

[0080] Where L0 is the initial length of the unit, θ X ,θ Y ,θ Z They are the angles between the local coordinate system and the global coordinate system's X, Y, and Z axes.

[0081] A3. Establish Hamilton's canonical equation

[0082] Specifically, the conjugate momentum of the unit is defined as The velocity vector of the cell is The kinetic energy of the unit is

[0083] The Hamiltonian function of the unit is

[0084]

[0085] The Hamiltonian canonical equation of the cable element can be obtained

[0086]

[0087] Among them, C e =αM e +βK e, α, β are damping coefficients.

[0088] A4. Using symplectic difference equations based on Hamiltonian system

[0089] Specifically, the n second-order differential equations of the cable system dynamics in the configuration space are simplified to 2n first-order canonical equations in the phase space by the Hamiltonian form. The obtained canonical equations can be naturally solved by integration using the symplectic integration format. The Hamiltonian canonical equations can be decomposed at k+1 time steps, and the following second-order symplectic scheme is obtained.

[0090]

[0091]

[0092] Where h is the computation time step. X and H P Respectively represent the partial differentials of the system Hamiltonian function with respect to the canonical variables X and P. The unit matrix is ​​in the form of:

[0093]

[0094]

[0095] A5. Conversion from local coordinates to global coordinates

[0096] Specifically, O-XYZ defines a global Cartesian inertial coordinate system with unit basis vectors (i, j, k) at the origin O. For any element e, the origin of the local coordinate system o-xyz is located at the starting point of the element. The basis vectors (e x ,e y ,e z ) consists of three Euler angles (θ x -Roll, theta y - Pitch, θ z -yaw). The order of transformation from the global coordinate system to the local coordinate system is defined as yaw, pitch, and roll respectively. The transformation is divided into three steps:

[0097] (1) Rotate the OZ axis counterclockwise at an angle of θ z (0≤θ z <π), so that the OY axis coincides with the ON axis.

[0098] (2) Rotate the new OY axis clockwise by an angle of θ y (0≤θ y <π), so that the OX axis in the global coordinate system coincides with the ox axis in the local coordinate system.

[0099] (3) Rotate the new OX axis counterclockwise by an angle of θ z (0≤θz <π), so that the OY axis in the global coordinate system coincides with the oy axis in the local coordinate system, and at the same time, the OZ axis coincides with the oz axis.

[0100] Through the above rotation, we can get the transformation matrix T from the global coordinate system to the local coordinate system

[0101] Once the transformation matrix from the global coordinate system to the local coordinate system is determined, the kinetic energy T of the unit e and elastic potential energy U e It can be expressed in the global coordinate system, and the displacement and velocity continuity boundary constraints are imposed in the unit equation to assemble the overall mass matrix and stiffness matrix of the cable system. According to the conversion relationship of the formula, the velocity vector of the two nodes of the unit in the local coordinate system can be converted into the expression in the global coordinate system, as shown below

[0102]

[0103] Among them, V e is the expression of the velocity vector of the two nodes of the unit in the global coordinate system, T e is a 6×6 unit transformation matrix. The unit kinetic energy can be expressed in the global coordinate system as

[0104]

[0105] And the unit mass matrix in the global coordinate system is obtained by this formula

[0106]

[0107] A6. Consider the work done by external forces

[0108] Specifically, for cable structures, the weight of the cable structure itself cannot be ignored. The OZ axis of the global coordinate system is in the vertical direction and opposite to the gravity direction. Obviously, the gravity vector of the unit node can be easily derived in the global coordinate system. Assuming that the undeformed cable structure is free to sag, the horizontal plane where the end node is located is the zero potential energy plane. The gravitational potential energy of the unit can be expressed in the global coordinate system as

[0109]

[0110]

[0111] Among them, N G is the linear interpolation function

[0112]

[0113] g,F Ge and L totalare the gravitational acceleration, the unit gravity vector and the total length of the initial cable structure. ξ is the unit local coordinate. X e is the unit node position vector in the global coordinate system. The Morrison equation is used to calculate the medium resistance. When the diameter is much smaller than the cable unit length, the corresponding unit resistance vector can be directly obtained in the global coordinate system. Applying the same unit linear interpolation shape function N, the velocity vector of the cable structure and the flow velocity of the medium are interpolated along the unit arc length, and the following equation can be obtained:

[0114] V f =NV fe

[0115] in, and V f are the movement speed of the cable structure and the flow velocity of the medium at that point respectively. fe is the medium velocity vector at the two nodes of the unit. The virtual work done by the drag force in the global coordinates can be expressed as

[0116]

[0117]

[0118] Among them, C d is the drag coefficient, ρ f is the medium density and L is the unit length.

[0119] A7. Unit assembly

[0120] Specifically, based on the continuity conditions of node displacement and velocity, the system unit array is assembled into the dynamic control equation of the overall cable system in the global coordinate system.

[0121]

[0122] M is the system mass matrix, K is the system stiffness matrix, C is the system damping matrix, F k is the system internal force matrix, F G is the system gravity matrix

[0123] exist Figure 1 Based on the nonlinear dynamic control equation of the towing system, combined with Figure 2 , the method for solving the nonlinear dynamic response of the towing system of the present invention is described, and the specific steps are as follows:

[0124] S1: Initialize the time and calculation steps with the initial value of 0; initialize the speed and displacement values ​​with the assigned initial value.

[0125] S2: Use the velocity displacement value of time step i as the velocity displacement value of time step i+1 to initialize the next time step.

[0126] S3: Construct a transformation matrix from local coordinates to global coordinates based on Euler angle transformation.

[0127] S4: Use node coordinates as variables to describe unit energy, and then further derive the unit's mass matrix, stiffness matrix, displacement matrix, damping matrix, and external force matrix.

[0128] Quality Array

[0129]

[0130] Stiffness Matrix

[0131]

[0132] Displacement Array

[0133] R=NX e

[0134] Damping Array

[0135] C e =αM e +βK e

[0136] External Force Array

[0137] F=F Ge +F de

[0138] S5: According to the continuity condition that the displacement velocities of different units at the connection nodes should remain consistent, the units are assembled into a cable-drag system and solved using the second-order symplectic difference algorithm.

[0139] S6: Determine the convergence of the solved displacement velocity result. If the convergence meets the requirement, proceed to step S7; if not, update the displacement velocity matrix and return to S5 to continue a new round of solution.

[0140] S7: Determine whether the calculation time reaches the total calculation time. If so, output the result and the solution ends. If not, enter the next time step for a new round of calculation.

Claims

1. A method for solving the nonlinear dynamic response of a towing system, characterized in that: Includes steps: Time and calculation steps are initialized to 0; speed and displacement values ​​are initialized to the given initial values; The velocity displacement value of time step i is used as the velocity displacement value of time step i+1 to initialize the next time step; Construct the towing system unit array based on the nodal coordinate finite element method; Based on the Euler angle transformation, a transformation matrix from local coordinates to global coordinates is constructed to determine the unit matrix of the towing system in the global coordinate system. Considering the external force to determine the external force matrix, the dynamic control equation of the towing system is constructed based on the towing system unit matrix in the global coordinate system; The response of the dynamic control equations is solved by the symplectic difference algorithm; If the convergence of the solved displacement velocity does not meet the requirements, the displacement velocity matrix is ​​updated and the symplectic difference algorithm is used to continue solving the dynamic control equation response; Otherwise, determine whether the calculation time reaches the set total calculation time. If it reaches the total calculation time, output the result and the solution ends. If the total calculation time has not been reached, a new round of calculation will be performed at the next time step; The construction of the towing system unit array based on the node coordinate finite element method specifically includes: Construct the element interpolation shape function N based on the nodal coordinate finite element method; The node coordinates are used as variables to represent the energy of the unit and construct the unit matrix; The unit matrix includes: Unit mass matrix M e : Where L is the unit length, ρ is the material density of the unit, and A is the cross-sectional area of ​​the unit; Element stiffness matrix K e : Where L0 is the initial length of the unit, θ X ,θ Y ,θ Z are the angles between the local coordinate system and the global coordinate system X, Y, and Z axes, namely the Euler angles, and E is the unit Young's modulus; Unit displacement matrix R: R=NX e Where R is the position vector of any point on the unit, X e are the position vectors of the two nodes of the unit; Unit damping array C e : C e =αM e +βK e Among them, α, β are damping coefficients; The external force matrix F is: F=F Ge +F de Among them, F Ge is the unit gravity vector in the global coordinate system, F de is the unit resistance in the global coordinate system, C d is the drag coefficient, ρ f is the medium density, L is the unit length, N G is the linear interpolation function, V f is the medium flow velocity of the cable structure at the node.

2. A method for solving the nonlinear dynamic response of a towing system according to claim 1, characterized in that: The transformation matrix from local coordinates to global coordinates based on Euler angle transformation is: Through the three Euler angles θ x -Roll, theta y - Pitch, θ z - Yaw describes the basis vectors of the local coordinate system (e x ,e y ,e z ); the order of transformation from the global coordinate system to the local coordinate system is defined as yaw, pitch and roll respectively; the transformation is divided into three steps: (1) Rotate the OZ axis counterclockwise at an angle of θ z , 0≤θ z <π, so that the OY axis coincides with the ON axis of the local coordinate system; (2) Rotate the new OY axis clockwise by an angle of θ y , 0≤θ y <π, so that the OX axis in the global coordinate system coincides with the ox axis in the local coordinate system; (3) Rotate the new OX axis counterclockwise by an angle of θ z , 0≤θ z <π, so that the OY axis in the global coordinate system coincides with the oy axis in the local coordinate system, and at the same time, the OZ axis coincides with the oz axis in the local coordinate system; Through the above rotation, the transformation matrix T from the global coordinate system to the local coordinate system is obtained.

3. A method for solving the nonlinear dynamic response of a towing system according to claim 1, characterized in that: The linear interpolation function is: where ξ is the local coordinate of the element.

4. A method for solving the nonlinear dynamic response of a towing system according to claim 1, characterized in that: The construction of the towing system dynamics control equation includes: Establish the Hamiltonian canonical equations of the towing system unit; Based on the continuity conditions of node displacement and velocity, the unit array in the overall coordinate system is assembled into the dynamic control equation of the overall towing system.

5. A method for solving the nonlinear dynamic response of a towing system according to claim 4, characterized in that: The Hamiltonian canonical equation of the towing system unit is: Among them, the conjugate momentum of the unit The velocity vector of the cell is The kinetic energy of the unit is He is the Hamiltonian function of the unit, F ke and F Ge They are the internal force matrix vectors and gravity vector of the two nodes of the unit respectively.

6. A method for solving the nonlinear dynamic response of a towing system according to claim 5, characterized in that: The dynamic control equation of the towing system is: Among them, M is the system mass matrix, K is the system stiffness matrix, C is the system damping matrix, and F k is the system internal force matrix, F G is the system gravity matrix.

7. A method for solving the nonlinear dynamic response of a towing system according to claim 6, characterized in that: The symplectic difference algorithm decomposes the Hamiltonian canonical equation at the k+1 time step to obtain the matrix: Here, h is the computation time step.

Images