A method for analyzing nonlinear dynamics and fatigue damage of a corrugated riser

Abstract

The application discloses a kind of non-linear dynamics and fatigue damage analysis method of slow wave shape riser, comprising: determining slow wave shape riser structure parameter and marine environment parameter, establishing three-dimensional model, the global position vector, velocity vector and acceleration vector of arbitrary point on unit are given after discretization processing is carried out to riser;Calculate the strain energy of unit mass matrix and riser unit, and deduce the generalized elastic force of unit;The environmental load suffered by riser is calculated, the system overall mass matrix is assembled, and the riser system dynamics equation set is established and converted into first-order ordinary differential equation set;The first-order ordinary differential equation set is solved, and the displacement, velocity and acceleration response time history of each node are obtained;The displacement, velocity and acceleration response time history under different working conditions are obtained, and the tension response of the key section of riser is obtained;The allowable cycle number corresponding to each tension cycle is determined.The application improves the accuracy and efficiency of the life cycle dynamic behavior and damage evolution prediction of SLWR.

Description

Technical Field

[0001] This invention relates to the technical field of dynamic analysis of marine engineering structures, and in particular to a nonlinear dynamics and fatigue damage analysis method for a gently sloping riser. Background Technology

[0002] Waveform risers, with their unique wavy configuration, effectively decouple platform motion and reduce dynamic response, making them a key piece of equipment in deepwater oil and gas development. However, their slender and flexible structural characteristics, under the coupled effects of waves, currents, and platform motion, exhibit strong geometric nonlinearity and complex dynamic behavior, posing significant challenges to structural design and safety assessment. Accurately predicting their dynamic response and the resulting fatigue damage remains a core challenge for the engineering community.

[0003] Currently, dynamic analysis of risers mainly relies on the finite element method (FEM) or the lumped mass method (LCM). The traditional FEM suffers from limitations in accuracy and efficiency when dealing with large rotations and deformations; while the LCM offers higher computational efficiency, it falls short in describing structural continuity and the details of local deformations. The absolute nodal coordinate method (ANCF), as an emerging flexible multibody dynamics approach, has inherent advantages in describing large deformation nonlinear problems and has been successfully applied in the analysis of structures such as steel catenary risers. However, existing research largely focuses on traditional riser forms. For structures like SLWRs with complex geometries (including suspended, buoyant, descending, and ground-lying sections), systematic overall dynamic and fatigue analysis based on ANCF remains a gap. In particular, existing methods often simplify platform motion to boundary conditions, failing to fully couple and consider the subtle influence of the amplitude, period, and waveform (regular and irregular waves) of the top surge excitation on the nonlinear dynamic transmission, stress distribution, and spatial evolution of fatigue damage throughout the SLWR. Furthermore, fatigue damage analysis is often treated as a subsequent processing step, failing to achieve deep integration and efficient collaboration with high-fidelity dynamic models. Summary of the Invention

[0004] Purpose of the invention: The purpose of this invention is to provide a nonlinear dynamics and fatigue damage analysis method for SLWR (Slow Wave Riser), which improves the accuracy and efficiency of predicting the dynamic behavior and damage evolution of SLWR throughout its entire life cycle.

[0005] Technical Solution: To achieve the above objectives, the present invention provides a nonlinear dynamics and fatigue damage analysis method for a gently sloping riser, comprising the following steps: (1) Determine the structural parameters of the wave-shaped riser system. The structural parameters of the wave-shaped riser include the coordinates of the upper end point, the coordinates of the lower end point, the riser length, the riser outer diameter, the riser wall thickness, the material density, the bending stiffness, the starting position of the buoyancy section, the arrangement length of the buoyancy section, the drag coefficient, and the inertia coefficient. (2) Set marine environmental parameters, including water depth, seawater density, ocean current velocity, peak surge excitation type, and peak surge excitation amplitude. and incentive cycle ; (3) In the coordinate system of the wave-shaped riser system, a three-dimensional model of the wave-shaped riser system is established based on the absolute node coordinate method ANCF. The riser is discretized by the reduced gradient ANCF beam element. The node coordinates and shape functions of the element are defined, and the global position vector, velocity vector and acceleration vector of any point on the element are given. (4) Construct the kinetic energy expression for the ANCF cable element and calculate the element mass matrix; (5) Based on the Euler-Bernoulli beam theory, calculate the strain energy of the riser unit, including axial strain energy and bending strain energy, and derive the generalized elastic force of the unit. (6) Calculate the environmental loads on the riser system. The environmental loads include gravity, buoyancy, seabed support force and ocean current force. The ocean current force is calculated using the Morison equation and includes tangential and normal components. (7) Assemble the overall mass matrix of the system. Based on the Lagrange equation and the principle of virtual work, establish the dynamic equations of the riser system. Use the Baumgarte method for numerical stabilization and transform it into a set of first-order ordinary differential equations. (8) The first-order ordinary differential equation system is solved by numerical integration to obtain the displacement, velocity and acceleration response time histories of each node of the riser system; (9) Modify the peak surge excitation amplitude set in step (2). With the excitation period T, steps (3) to (8) are repeated based on the modified sea state environmental parameters to obtain the displacement, velocity and acceleration response time histories under different working conditions, and to obtain the tension response of the key section of the riser. (10) Based on the tension data of the key section of the riser, the tension amplitude and tension mean are extracted by the rain flow counting method. Based on the tension-life curve of the riser material, fatigue damage is calculated to determine the allowable number of cycles corresponding to each tension cycle.

[0006] Optionally, in step (1), the wave-shaped riser system includes a wave-shaped riser, a platform, and a seabed. The upper end of the wave-shaped riser is connected to the platform, and the lower end of the wave-shaped riser is connected to the seabed. The wave-shaped riser includes a suspended section, a buoyancy section, a descent section, and a ground-lying section. The contact point between the wave-shaped riser and the seabed is the contact point, and the riser between the contact point and the lower end is the ground-lying section.

[0007] Optionally, the peak surge excitation type in step (2) includes regular waves and irregular waves, and the peak surge excitation amplitude... The value range is from 0m to 15m, and the excitation period is... The value range is from 10s to 30s.

[0008] Optionally, step (3) specifically includes the following steps:

[0009] The coordinate system of the wave-shaped riser system is as follows: the center point where the platform intersects with the sea level is the origin of the XYZ coordinate system, the sea level is the XY plane, the ocean current direction is set along the positive X-axis, and the Z-axis is perpendicular to the XY plane; the top swell response of the platform under wave action is simulated by applying simple harmonic motion along the X-axis at the top of the platform.

[0010] The nodal coordinates of an ANCF reduced beam element consist of the position vectors of the nodes at both ends of the element and their gradient vectors with respect to the local coordinates of the element, and the global position vector of any point P on the element. Represented using shape function matrices and node coordinate vectors: ,

[0011] in, Represents the coordinates of the riser unit nodes. The shape function matrix representing the riser unit;

[0012] The shape function matrix of the unit From the local coordinates of the element The interpolation function is composed of, , The initial length of the cell is given by: ,

[0013] Where I represents the third-order identity matrix; Represents a dimensionless parameter. ;

[0014] Then interpolation function , , and The specific expression is: ,

[0015] Node coordinate vector of the element From the two end nodes of the unit and Position vector and its relationship to local coordinates The gradient vector is composed of, specifically, the following: , in, Representative node and The global position vector; These are the local coordinates of the position vector. The gradient vector; represent The transpose of the matrix;

[0016] Since the shape function is independent of time, the velocity vector at any point on the riser... and acceleration vector To express as: ,

[0017] in, , These are the first derivative (velocity) and second derivative (acceleration) of the element node coordinate vector with respect to time, respectively, which represent the global position vector of any point on the element. Velocity vector and acceleration vector .

[0018] Optionally, step (4) specifically includes the following steps:

[0019] Unit kinetic energy The expression is: ,

[0020] in, The cross-sectional area of ​​the representative unit; The density of the constituent elements represents the material density; the motion state of any point on the element is described by its vector in the global coordinate system: Let be the velocity vector at that point, and its transpose be denoted as . ; Let be the acceleration vector at that point, and let its transpose be denoted as . In the shape function representation based on the absolute nodal coordinate method, Let be the shape function matrix of the selected element, and its transpose be denoted as . ;

[0021] The element mass matrix is ​​derived based on the element characteristics. By integrating the shape function, we obtain: .

[0022] Optionally, step (5) specifically includes the following steps:

[0023] Unit strain energy It can be determined by axial strain and curvature express: ,

[0024] Where E represents the material's elastic modulus and J represents the moment of inertia of the reduced beam element;

[0025] in, and The calculation formulas are as follows: , ,

[0026] in, Represents the position vector of the cell local coordinates of the element The first-order partial derivative, This is the transpose of the gradient vector;

[0027] Generalized elastic force of a unit Through strain energy Find the coordinates of the nodes. The partial derivatives yield: ,

[0028] in, The axial strain ε is the absolute nodal coordinate vector. The partial derivatives yield the transpose matrix; For curvature absolute node coordinate vector The partial derivatives yield the transpose matrix. , , .

[0029] Optionally, step (6) specifically includes the following steps:

[0030] The Morison equation is used to calculate the ocean current force, and its tangential component is... With normal components The expression is: , ,

[0031] Where D represents the outer diameter of the riser; ρ s This indicates the density of the seawater in the environment where the riser is located; , , These represent the relative velocity vector between the seawater and the riser, its tangential component on the riser surface, and its normal component on the riser surface, respectively. , These are the fluid load loading functions corresponding to the tangential and normal directions, respectively, used to calculate the hydrodynamic loads acting on the riser.

[0032] Seabed support Simulation using a linear spring model: ,

[0033] in, This indicates the depth of the riser structure relative to the seabed reference level; It represents the equivalent contact stiffness of the contact area between the riser and the seabed, and is used to characterize the seabed's resistance to structural deformation; The depth of the seabed is the vertical distance from the sea level to the surface of the seabed. Represents the unit direction vector along the Z-axis in the global coordinate system; This represents the position vector of any point P on the riser structure, used to describe the spatial coordinates of that point in the global coordinate system;

[0034] For gravity and buoyancy The loads acting on the riser unit can be expressed as follows: , ,

[0035] Generalized force vector Including the system's generalized elastic force and system generalized external forces , ,

[0036] in, Including the tangential component of ocean current force Ocean current force normal component Seabed support ,gravity and buoyancy .

[0037] Optionally, step (7) specifically includes the following steps:

[0038] Taking the mass matrix of a unit as an example, after the unit assembly is completed, the mass matrix of the system composed of N units is... Represented as: ,

[0039] in, This represents the unit mass matrix of the Nth unit in the riser system;

[0040] The unit global position vector of N units Assembled into system generalized coordinates :

[0041] The dynamic equations of the riser system can be expressed in the following matrix form: ,

[0042] in, The Jacobian matrix (i.e., the first-order partial derivative matrix) of the complete constraint equations with respect to the system's generalized coordinates q describes the sensitivity of the constraints to changes in the system's configuration. This is the transpose of the Jacobian matrix, used to project the constraint reaction force onto the generalized force direction of the system; This represents the generalized acceleration vector of the system; For each constraint, represents the Lagrange multiplier vector, and its physical meaning is the magnitude of the constraint force; Q represents the generalized force vector acting on the system. The right-hand side of the system dynamics equations contains nonlinear coupling terms caused by velocity and contributions from constraints that are explicitly time-dependent. Together, they constitute the constraint compatibility conditions at the acceleration level in numerical solutions.

[0043] in, Expanding the right-hand side of the system dynamics equations corrected by Baumgarte's method, we get: ,

[0044] First-order partial derivatives of the constraint equations in time for: ,

[0045] Solving again yields the second-order partial derivatives. for: , in: , ,

[0046] in, This represents the explicit partial derivative of the constraint with respect to time t, reflecting the characteristics of the constraint itself as it changes over time. It represents the second-order explicit partial derivative of the constraint with respect to time, used to describe the acceleration effect of the constraint's time-varying rate; and These are the generalized velocity vector and the generalized acceleration vector of the system, respectively. Representing vectors The partial derivative with respect to the generalized coordinate q, this term reflects the nonlinear spatial variation of velocity in the acceleration equation; The mixed partial derivative is the constraint derivative taken with respect to coordinate q and then with respect to time t. It describes the rate of change of the constraint gradient with time and corresponds to the velocity-related coupling term in the dynamic equation.

[0047] Optionally, the method for solving the system dynamics equations in step (8) is as follows:

[0048] The Baumgarte solution method is used to solve the dynamic equations. In order to maintain the stability of the solution results, a feedback system can be introduced into the constraint equations in the Baumgarte solution method.

[0049] ,

[0050] Here, α and β are Baumgarte stabilization parameters.

[0051] The final dynamic equations of the riser model are expressed in the following form: ,

[0052] The dynamic equations of the final riser model are solved to obtain the displacement, velocity, and acceleration response time histories of each node in the riser system.

[0053] Optionally, in step (10), fatigue damage calculation uses the SN curve and Miner's linear cumulative damage criterion: ,

[0054] in, Represents the cumulative damage level; Represents the total number of load levels; Represents the current load level sequence number; This represents the actual number of cycles experienced at the i-th load level; Representative at the The critical number of cycles required to cause failure at a given load level.

[0055] Beneficial Effects: Compared with existing technologies, this invention has the following significant advantages: This invention integrates the advantages of the absolute nodal coordinate method in high-precision simulation of geometric nonlinearity, and systematically integrates environmental load calculation, constraint processing, and fatigue assessment, aiming to improve the accuracy and efficiency of predicting the dynamic behavior and damage evolution of SLWRs throughout their entire life cycle; This invention establishes a high-fidelity geometric nonlinear dynamic model of the gently sloping riser based on the absolute nodal coordinate method, which can accurately simulate its large displacement, large rotation, and large deformation behavior under the combined action of tip surge and ocean current, thereby significantly improving the simulation accuracy of riser dynamic response in complex marine environments; This invention combines nonlinear dynamic analysis and fatigue damage assessment. Seamless integration throughout the entire process enables integrated analysis from environmental load input to fatigue life prediction, avoiding data transmission errors in traditional multi-software coupled analysis and improving the completeness and reliability of the analysis. Specifically targeting the wavy configuration characteristics of gently sloping risers, this invention clearly reveals the differentiated influence of top excitation parameters on the dynamic response and stress distribution of different structural sections, providing a direct basis for local strengthening and overall optimization of SLWR. Through systematic research on the mechanical behavior of the buoyancy section under regular and irregular wave excitation, this invention quantitatively elucidates its dual mechanism of "local stress concentration and overall fatigue improvement," deepening the understanding of the design principles and performance advantages of SLWR. Attached Figure Description

[0056] Figure 1 A schematic diagram of the wave-shaped riser system provided by the present invention; Figure 2 This is a model diagram of the riser unit provided by the present invention; Figure 3 A schematic diagram of the overall displacement of the flexible pendulum at different times provided by the present invention; Figure 4 A schematic diagram showing the change of the Y-direction displacement of the end point of the flexible pendulum over time, provided by the present invention. Figure 5 This is a schematic diagram of the average tension of the wave-shaped riser provided by the present invention in the parameter space of excitation amplitude Af and riser length L; Figure 6 A schematic diagram of the average tension of the wave-shaped riser in the parameter space of excitation period T and riser length L provided by the present invention. Figure 7 A schematic diagram of fatigue damage of a wave-shaped riser in the parameter space of excitation amplitude Af and riser length L, provided by the present invention. Figure 8 This is a schematic diagram of fatigue damage of the wave-shaped riser provided by the present invention in the parameter space of excitation period T and riser length L. Detailed Implementation

[0057] The technical solution of the present invention will be further described below with reference to the accompanying drawings.

[0058] like Figure 1 As shown, the nonlinear dynamics and fatigue damage analysis method for a gently sloping riser of the present invention includes the following steps:

[0059] (1) Determine the structural parameters of the wave-shaped riser system. The wave-shaped riser system includes wave-shaped risers, a platform, and the seabed. The upper end of the wave-shaped riser is connected to the platform, and the lower end of the wave-shaped riser is connected to the seabed. The wave-shaped riser includes a suspended section, a buoyancy section, a descent section, and a ground-lying section. The contact point between the wave-shaped riser and the seabed is the contact point, and the riser between the contact point and the lower end is the ground-lying section. The structural parameters of the wave-shaped riser include the coordinates of the upper end and the lower end, the riser length, the riser outer diameter, the riser wall thickness, the material density, the bending stiffness, the starting position of the buoyancy section, the arrangement length of the buoyancy section, the drag coefficient, and the inertia coefficient. The polar moment of inertia of the section is J = 5.4094 x 10⁻⁶. -4 m 4 Elastic modulus E = 4.9001 x 10⁻⁶ 4 MPa, the specific parameters are shown in Table 1.

[0060] Table 1 Parameter Table

[0061]

[0062] (2) Set marine environmental parameters, including water depth, seawater density, ocean current velocity, peak surge excitation type, and peak surge excitation amplitude. and incentive cycle The peak surge excitation types include regular waves (Airy waves) and irregular waves (Jonswap waves); the peak surge excitation amplitude... The value range is from 0m to 15m, and the excitation period is... The value ranges from 10s to 30s, and the incident angle is 180°.

[0063] The starting position and arrangement length of the buoyancy section are determined based on the buoyancy ratio. The buoyancy ratio is defined as the ratio of the net buoyancy of the buoyancy section to the net underwater weight of the riser, with a value range of 10%–50%. The equivalent buoyancy distribution of the buoyancy section is calculated based on the buoyancy ratio, and the arrangement range of the buoyancy section in the riser axis is determined.

[0064] In the design of the gently sloping riser, the length of the buoyancy section should be selected within the range of 150 to 250 meters, and the starting position is recommended to be within the range of 700 to 900 meters from the upper end.

[0065] Specifically, increasing the arrangement length significantly improves net buoyancy and increases the riser waveform curvature, which helps reduce top tension. However, it also increases local tension and bending curvature at the end of the buoyancy section and adjacent descending sections. Lowering the starting position of the buoyancy section expands the horizontal span of the riser, making the configuration more gradual, thus increasing tension at the upper end, while correspondingly reducing tension and maximum curvature at the end of the buoyancy section. Therefore, the arrangement length and starting position need to be optimized synergistically to balance the overall tension distribution, control local bending stress, and effectively block the transmission of upper dynamic response to the contact area, ultimately improving the safety and fatigue performance of the gentle-wave riser system under complex sea conditions.

[0066] (3) such as Figure 3 As shown, in the coordinate system of the wave-shaped riser system, a three-dimensional model of the wave-shaped riser system is established based on the absolute node coordinate method ANCF. The riser is discretized using reduced gradient ANCF beam elements. The node coordinates and shape functions of the elements are defined, and the global position vector, velocity vector and acceleration vector of any point on the element are given.

[0067] Step (3) specifically includes the following steps:

[0068] The coordinate system of the wave-shaped riser system is as follows: the center point where the platform intersects with the sea level is the origin of the XYZ coordinate system, the sea level is the XY plane, the ocean current direction is set along the positive X-axis, and the Z-axis is perpendicular to the XY plane; the top swell response of the platform under wave action is simulated by applying simple harmonic motion along the X-axis at the top of the platform.

[0069] The nodal coordinates of an ANCF reduced beam element consist of the position vectors of the nodes at both ends of the element and their gradient vectors with respect to the local coordinates of the element, and the global position vector of any point P on the element. Represented using shape function matrices and node coordinate vectors: ,

[0070] in, Represents the coordinates of the riser unit nodes. The shape function matrix representing the riser unit;

[0071] The shape function matrix of the unit From the local coordinates of the element The interpolation function is composed of, , The initial length of the cell is given by: ,

[0072] Where I represents the third-order identity matrix; Represents a dimensionless parameter. ; Then interpolation function , , and The specific expression is: ,

[0073] Node coordinate vector of the element From the two end nodes of the unit and Position vector and its relationship to local coordinates The gradient vector is composed of, specifically, the following: , in, Representative node and The global position vector; These are the local coordinates of the position vector. The gradient vector; represent The transpose of the matrix;

[0074] Since the shape function is independent of time, the velocity vector at any point on the riser... and acceleration vector To express as: ,

[0075] in, , These are the first derivative (velocity) and second derivative (acceleration) of the element node coordinate vector with respect to time, respectively, which represent the global position vector of any point on the element. Velocity vector and acceleration vector ;

[0076] (4) Construct the kinetic energy expression for the ANCF cable element and calculate the element mass matrix;

[0077] Step (4) specifically includes the following steps:

[0078] Unit kinetic energy The expression is: ,

[0079] in, The cross-sectional area of ​​the representative unit; The density of the constituent elements represents the material density; the motion state of any point on the element is described by its vector in the global coordinate system: Let be the velocity vector at that point, and its transpose be denoted as . ; Let be the acceleration vector at that point, and let its transpose be denoted as . In the shape function representation based on the absolute nodal coordinate method, Let be the shape function matrix of the selected element, and its transpose be denoted as . ;

[0080] The element mass matrix is ​​derived based on the element characteristics. By integrating the shape function, we obtain: ,

[0081] (5) Based on the Euler-Bernoulli beam theory, calculate the strain energy of the riser unit, including axial strain energy and bending strain energy, and derive the generalized elastic force of the unit.

[0082] Step (5) specifically includes the following steps:

[0083] Unit strain energy It can be determined by axial strain and curvature express: , Where E represents the material's elastic modulus and J represents the moment of inertia of the reduced beam element; in, and The calculation formulas are as follows: , , in, Represents the position vector of the cell local coordinates of the element The first-order partial derivative, This is the transpose of the gradient vector; Generalized elastic force of a unit Through strain energy Find the coordinates of the nodes. The partial derivatives yield: , in, The axial strain ε is the absolute nodal coordinate vector. The partial derivatives yield the transpose matrix; For curvature absolute node coordinate vector The partial derivatives yield the transpose matrix. , , , (6) Calculate the environmental loads on the riser system. The environmental loads include gravity, buoyancy, seabed support force and ocean current force. The ocean current force is calculated using the Morison equation and includes tangential and normal components.

[0084] Step (6) specifically includes the following steps: The Morison equation is used to calculate the ocean current force, and its tangential component is... With normal components The expression is: , , Where D represents the outer diameter of the riser; ρ s This indicates the density of the seawater in the environment where the riser is located; , , These represent the relative velocity vector between the seawater and the riser, its tangential component on the riser surface, and its normal component on the riser surface, respectively. , These are the fluid load loading functions corresponding to the tangential and normal directions, respectively, used to calculate the hydrodynamic loads acting on the riser. Seabed support Simulation using a linear spring model: , in, This indicates the depth of the riser structure relative to the seabed reference level; It represents the equivalent contact stiffness of the contact area between the riser and the seabed, and is used to characterize the seabed's resistance to structural deformation; The depth of the seabed is the vertical distance from the sea level to the surface of the seabed. Represents the unit direction vector along the Z-axis in the global coordinate system; This represents the position vector of any point P on the riser structure, used to describe the spatial coordinates of that point in the global coordinate system; For gravity and buoyancy The loads acting on the riser unit can be expressed as follows: , , Generalized force vector Including the system's generalized elastic force and system generalized external forces , , in, Including the tangential component of ocean current force Ocean current force normal component Seabed support ,gravity and buoyancy ;

[0085] (7) Assemble the overall mass matrix of the system. Based on the Lagrange equation and the principle of virtual work, establish the dynamic equations of the riser system. Use the Baumgarte method for numerical stabilization and transform it into a set of first-order ordinary differential equations.

[0086] Step (7) specifically includes the following steps: Taking the mass matrix of a unit as an example, after the unit assembly is completed, the mass matrix of the system composed of N units is... Represented as: , in, This represents the unit mass matrix of the Nth unit in the riser system; The unit global position vector of N units Assembled into system generalized coordinates : The dynamic equations of the riser system can be expressed in the following matrix form: ,

[0087] in, The Jacobian matrix (i.e., the first-order partial derivative matrix) of the complete constraint equations with respect to the system's generalized coordinates q describes the sensitivity of the constraints to changes in the system's configuration. This is the transpose of the Jacobian matrix, used to project the constraint reaction force onto the generalized force direction of the system; This represents the generalized acceleration vector of the system; For each constraint, represents the Lagrange multiplier vector, and its physical meaning is the magnitude of the constraint force; Q represents the generalized force vector acting on the system. The right-hand side of the system dynamics equations contains nonlinear coupling terms caused by velocity and contributions from constraints that are explicitly time-dependent. Together, these terms constitute the constraint compatibility conditions at the acceleration level in the numerical solution.

[0088] in, Expanding the right-hand side of the system dynamics equations corrected by Baumgarte's method, we get: , First-order partial derivatives of the constraint equations in time for: , Solving again yields the second-order partial derivatives. for: , in: , ,

[0089] in, This represents the explicit partial derivative of the constraint with respect to time t, reflecting the characteristics of the constraint itself as it changes over time. It represents the second-order explicit partial derivative of the constraint with respect to time, used to describe the acceleration effect of the constraint's time-varying rate; and These are the generalized velocity vector and the generalized acceleration vector of the system, respectively. Representing vectors The partial derivative with respect to the generalized coordinate q, this term reflects the nonlinear spatial variation of velocity in the acceleration equation; The mixed partial derivative is the constraint derivative taken with respect to coordinate q and then with respect to time t. It describes the rate of change of the constraint gradient with time and corresponds to the velocity-related coupling term in the dynamic equation.

[0090] (9) The first-order ordinary differential equation system is solved by numerical integration to obtain the displacement, velocity and acceleration response time histories of each node of the riser system;

[0091] In step (9), the system dynamic equations are solved as follows:

[0092] The Baumgarte solution method is used to solve the dynamic equations. In order to maintain the stability of the solution results, a feedback system can be introduced into the constraint equations in the Baumgarte solution method.

[0093] ,

[0094] Here, α and β are Baumgarte stabilization parameters.

[0095] The final dynamic equations of the riser model are expressed in the following form: ,

[0096] The dynamic equations of the final riser model are solved to obtain the displacement, velocity, and acceleration response time histories of each node in the riser system;

[0097] (10) Modify the peak surge excitation amplitude set in step (2). With the excitation period T, steps (3) to (9) are repeated based on the modified sea state environmental parameters to obtain the displacement, velocity and acceleration response time histories under different working conditions, and to obtain the tension response of the key section of the riser.

[0098] (10) Based on the tension data of the key section of the riser, the tension amplitude and tension mean are extracted by the rain flow counting method. Based on the tension-life curve of the riser material, fatigue damage is calculated to determine the allowable number of cycles corresponding to each tension cycle.

[0099] In step (10), fatigue damage calculation uses the SN curve and the Miner linear cumulative damage criterion: ,

[0100] in, Represents the cumulative damage level; Represents the total number of load levels; Represents the current load level sequence number; This represents the actual number of cycles experienced at the i-th load level; Representative at the The critical number of cycles required to cause failure at a given load level.

[0101] like Figure 3 The diagram shows the overall displacement of the flexible pendulum at different times provided by the present invention. Figure 4 The diagram shows the Y-direction displacement of the end point of the flexible pendulum provided by this invention over time. The program developed in this invention is compared with an existing ANCF-based flexible pendulum model, verifying the correctness of the ANCF model. Furthermore, even with a smaller number of elements, the ANCF flexible pendulum exhibits high accuracy, further demonstrating the applicability of the ANCF method in describing the large displacement and deformation behavior of flexible bodies.

[0102] like Figure 5 The diagram shows the average tension of the wave-shaped riser provided by this invention in the parameter space of excitation amplitude Af and riser length L. The tension initially decreases linearly from the upper end; a significant increase in tension is observed in the buoyancy section, followed by an approximately linear decay and eventually stabilization. As the excitation amplitude Af increases, the tension level of SLWR increases overall at all positions along its length, but the rate of increase shows obvious spatial distribution differences.

[0103] like Figure 6 The diagram shows the average tension of the wave-shaped riser provided by this invention in the parameter space of excitation period T and riser length L. According to the tension distribution, it can be observed that the tension decreases approximately linearly from the top end along the riser axis, but a significant tension surge occurs in the buoyancy section, forming a local peak, after which the decreasing trend resumes and gradually stabilizes. Analysis shows that the change in period T has no significant impact on the overall tension distribution of the SLWR, and no systematic increasing trend with changing period was observed.

[0104] like Figure 7The diagram shows the fatigue damage of the gently sloping riser provided by this invention in the parameter space of excitation amplitude Af and riser length L. The fatigue damage is mainly concentrated near the upper end, with secondary local fatigue damage concentration areas also existing at L = 700 m and 2098 m. Furthermore, under conditions of larger Af (especially Af > 8 m), a significant fatigue damage area appears near the upper end of the riser. Analysis shows that although adding a buoyancy section creates local stress concentration areas, this design effectively reduces the overall fatigue damage of the riser, achieving optimization of the relationship between local and global factors.

[0105] like Figure 8 The diagram shows the fatigue damage of the wave-shaped riser provided by this invention in the parameter space of excitation period T and riser length L. The maximum fatigue damage still occurs near the upper end, while two other minor local fatigue damage concentration areas exist at L = 700 m and 2098 m. The fatigue damage is alleviated as T increases, and its growth rate also gradually decreases. Meanwhile, under smaller T conditions (especially T < 15 s), a significant fatigue damage area appears near the suspension point at the top of the riser. This indicates that although the addition of the buoyancy section introduces local fatigue damage, it brings a global improvement in the overall fatigue resistance of the riser.

[0106] This invention is the first to systematically introduce the Absolute Nodal Coordinates (ANCF) reduced gradient beam element system into the nonlinear dynamics and fatigue damage analysis of sloping wave risers (SLWR). It uses the nodal position vector and its first-order gradient vector with respect to the element's local coordinates as the nodal coordinate vector, and achieves high-precision interpolation of the global position vector at any point within the element through a cubic Hermite shape function. This modeling method eliminates the need for rotational coordinate linearization approximation, fundamentally avoiding the rotation singularity and stiffness matrix distortion problems inherent in traditional finite element methods for large displacement and rotation analyses. It can accurately describe the large displacement, rotation, and deformation behavior of SLWR under the combined effects of tip surge and ocean currents. Compared to the commercial software OrcaFlex lumped mass method, this invention requires only 90 ANCF reduced beam elements to achieve convergence accuracy, while the lumped mass method typically requires approximately 360 elements. This invention achieves equal or even better computational accuracy with less than one-third the number of elements, significantly improving the computational efficiency of long-term dynamic simulation and fatigue analysis of deep-water risers.

[0107] In terms of environmental load and excitation modeling, this invention systematically integrates two peak surge excitation modes: regular waves (Airy waves) and irregular waves (JONSWAP waves), fully preserving the complete parameter space of excitation amplitude Af (0m~15m) and excitation period T (10s~30s). Based on the Morison equation, it accurately calculates the tangential and normal components of the ocean current force and uses a linear spring model to simulate the normal support force of the seabed, achieving high-fidelity dynamic response simulation under strong fluid-structure-seabed coupling. Based on this platform, this invention systematically reveals for the first time that an increase in excitation amplitude Af significantly exacerbates tension fluctuations at the suspension point, contact point, and anchoring point, and has a dominant effect on the dynamic response of the suspension segment; changes in wave period T have no significant impact on the distribution of tension and bending moment, but have a clear inhibitory effect on fatigue damage. The above refined parameter research fills the gap in the systematic study of peak excitation parameter space in existing SLWR dynamic analysis.

[0108] This invention seamlessly integrates nonlinear dynamic analysis and fatigue damage assessment throughout the entire process. Within a single ANCF numerical framework, it achieves a closed-loop analysis encompassing everything from top surge excitation input, nonlinear dynamic response solving, and key section stress time history extraction, to rainflow counting, SN curve matching, and Miner's linear cumulative damage criterion. This integrated strategy completely eliminates data transfer errors and model mismatch risks inherent in traditional multi-software coupling models, significantly improving the integrity and traceability of long-term fatigue assessment. Based on this, this invention provides the first quantitative explanation of the dual mechanism of "local stress concentration and overall fatigue improvement" in the buoyancy section: although significant compressive stress concentration areas form locally in the buoyancy section, the transmission of large-scale top movement to the ground contact section is effectively isolated by reconstructing the overall riser configuration and power transmission path. This results in a significant reduction in fatigue damage in key areas such as the suspension point and contact point of the SLWR compared to traditional steel catenary risers (SCR) under the same operating conditions. This quantitative disclosure of the mechanism provides a clear mechanical basis for the trade-off optimization of SLWR buoyancy section parameters.

[0109] To address the inherent index-3 differential-algebraic equations of the ANCF dynamic equations, this invention employs an augmented method and introduces a Baumgarte violation stabilization term. By introducing feedback corrections (stability coefficients α and β) for position-level and velocity-level constraint violations into the right-hand side of the acceleration constraint equations, the original index-3 DAEs are transformed into a directly integrable set of first-order ordinary differential equations. This effectively suppresses constraint drift during long-term numerical integration and significantly improves the numerical stability of long-term dynamic simulations under wave-platform coupled excitation. This method does not rely on specific commercial software; the entire solution process is autonomously implemented on general-purpose numerical platforms such as MATLAB, facilitating customized extensions by engineers based on actual project needs. It provides a reliable technical tool for the autonomous design, integrity management, and life extension assessment of deep-sea SLWR systems.

[0110] In summary, this invention significantly outperforms existing technologies in terms of modeling accuracy, geometric nonlinearity adaptability, environmental load completeness, parameter research depth, integrated fatigue assessment, and numerical stability. It has important theoretical value and broad engineering application prospects for promoting the transformation of deep-sea wave-shaped riser design methods from empirical analogy to refined numerical simulation.

Claims

1. A nonlinear dynamics and fatigue damage analysis method for a gently sloping riser, characterized in that, Includes the following steps: (1) Determine the structural parameters of the wave-shaped riser system. The structural parameters of the wave-shaped riser include the coordinates of the upper end point, the coordinates of the lower end point, the riser length, the riser outer diameter, the riser wall thickness, the material density, the bending stiffness, the starting position of the buoyancy section, the arrangement length of the buoyancy section, the drag coefficient, and the inertia coefficient. (2) Set marine environmental parameters, including water depth, seawater density, ocean current velocity, peak surge excitation type, and peak surge excitation amplitude. and incentive cycle ; (3) In the coordinate system of the wave-shaped riser system, a three-dimensional model of the wave-shaped riser system is established based on the absolute node coordinate method ANCF. The riser is discretized by the reduced gradient ANCF beam element. The node coordinates and shape functions of the element are defined, and the global position vector, velocity vector and acceleration vector of any point on the element are given. (4) Construct the kinetic energy expression for the ANCF cable element and calculate the element mass matrix; (5) Based on the Euler-Bernoulli beam theory, calculate the strain energy of the riser unit, including axial strain energy and bending strain energy, and derive the generalized elastic force of the unit. (6) Calculate the environmental loads on the riser system. The environmental loads include gravity, buoyancy, seabed support force and ocean current force. The ocean current force is calculated using the Morison equation and includes tangential and normal components. (7) Assemble the overall mass matrix of the system. Based on the Lagrange equation and the principle of virtual work, establish the dynamic equations of the riser system. Use the Baumgarte method for numerical stabilization and transform it into a set of first-order ordinary differential equations. (8) The first-order ordinary differential equation system is solved by numerical integration to obtain the displacement, velocity and acceleration response time histories of each node of the riser system; (9) Modify the peak surge excitation amplitude set in step (2). With the excitation period T, steps (3) to (8) are repeated based on the modified sea state environmental parameters to obtain the displacement, velocity and acceleration response time histories under different working conditions, and to obtain the tension response of the key section of the riser. (10) Based on the tension data of the key section of the riser, the tension amplitude and tension mean are extracted by the rain flow counting method. Based on the tension-life curve of the riser material, fatigue damage is calculated to determine the allowable number of cycles corresponding to each tension cycle.

2. The nonlinear dynamics and fatigue damage analysis method for a gently sloping riser according to claim 1, characterized in that, The wave-shaped riser system in step (1) includes a wave-shaped riser, a platform and a seabed. The upper end of the wave-shaped riser is connected to the platform and the lower end of the wave-shaped riser is connected to the seabed. The wave-shaped riser includes a suspended section, a buoyancy section, a descent section and a ground-lying section. The contact point between the wave-shaped riser and the seabed is the contact point, and the riser between the contact point and the lower end is the ground-lying section.

3. The nonlinear dynamics and fatigue damage analysis method for a gently sloping riser according to claim 2, characterized in that, The peak surge excitation type in step (2) includes regular waves and irregular waves, and the peak surge excitation amplitude... The value range is from 0m to 15m, and the excitation period is... The value range is from 10s to 30s.

4. The nonlinear dynamics and fatigue damage analysis method for a gently sloping riser according to claim 3, characterized in that, Step (3) specifically includes the following steps: The coordinate system of the wave-shaped riser system is as follows: the center point where the platform intersects with the sea level is the origin of the XYZ coordinate system, the sea level is the XY plane, the ocean current direction is set along the positive X-axis, and the Z-axis is perpendicular to the XY plane; The tip surge response of the platform under wave action is simulated by applying a simple harmonic motion along the X-axis at the top of the platform. The nodal coordinates of an ANCF reduced beam element consist of the position vectors of the nodes at both ends of the element and their gradient vectors with respect to the local coordinates of the element, and the global position vector of any point P on the element. Represented using shape function matrices and node coordinate vectors: ， in, Represents the coordinates of the riser unit nodes. The shape function matrix representing the riser unit; The shape function matrix of the unit From the local coordinates of the element The interpolation function is composed of, , The initial length of the cell is given by: ， Where I represents the third-order identity matrix; Represents a dimensionless parameter. ; Then interpolation function , , and The specific expression is: ， Node coordinate vector of the element From the two end nodes of the unit and Position vector and its relationship to local coordinates The gradient vector is composed of, specifically, the following: ， in, Representative node and The global position vector; These are the local coordinates of the position vector. The gradient vector; represent The transpose of the matrix; Since the shape function is independent of time, the velocity vector at any point on the riser... and acceleration vector To express as: ， in, , These are the first derivative (velocity) and second derivative (acceleration) of the element node coordinate vector with respect to time, respectively, which represent the global position vector of any point on the element. Velocity vector and acceleration vector .

5. The nonlinear dynamics and fatigue damage analysis method for a gently sloping riser according to claim 4, characterized in that, Step (4) specifically includes the following steps: Unit kinetic energy The expression is: ， in, The cross-sectional area of ​​the representative unit; The density of the constituent elements represents the material density; the motion state of any point on the element is described by its vector in the global coordinate system: Let be the velocity vector at that point, and its transpose be denoted as . ; Let be the acceleration vector at that point, and let its transpose be denoted as . In the shape function representation based on the absolute nodal coordinate method, Let be the shape function matrix of the selected element, and its transpose be denoted as . ; The element mass matrix is ​​derived based on the element characteristics. By integrating the shape function, we obtain: 。 6. The nonlinear dynamics and fatigue damage analysis method for a gently sloping riser according to claim 5, characterized in that, Step (5) specifically includes the following steps: Unit strain energy It can be determined by axial strain and curvature express: ， Where E represents the material's elastic modulus and J represents the moment of inertia of the reduced beam element; in, and The calculation formulas are as follows: ， ， in, Represents the position vector of the cell local coordinates of the element The first-order partial derivative, This is the transpose of the gradient vector; Generalized elastic force of a unit Through strain energy Find the coordinates of the nodes. The partial derivatives yield: ， in, The axial strain ε is the absolute nodal coordinate vector. The partial derivatives yield the transpose matrix; For curvature absolute node coordinate vector The partial derivatives yield the transpose matrix. ， ， 。 7. The nonlinear dynamics and fatigue damage analysis method for a gently sloping riser according to claim 6, characterized in that, Step (6) specifically includes the following steps: The Morison equation is used to calculate the ocean current force, and its tangential component is... With normal components The expression is: ， ， Where D represents the outer diameter of the riser; ρ s This indicates the density of the seawater in the environment where the riser is located; , , These represent the relative velocity vector between the seawater and the riser, its tangential component on the riser surface, and its normal component on the riser surface, respectively. , These are the fluid load loading functions corresponding to the tangential and normal directions, respectively, used to calculate the hydrodynamic loads acting on the riser. Seabed support Simulation using a linear spring model: ， in, This indicates the depth of the riser structure relative to the seabed reference level; It represents the equivalent contact stiffness of the contact area between the riser and the seabed, and is used to characterize the seabed's resistance to structural deformation; The depth of the seabed is the vertical distance from the sea level to the surface of the seabed. Represents the unit direction vector along the Z-axis in the global coordinate system; This represents the position vector of any point P on the riser structure, used to describe the spatial coordinates of that point in the global coordinate system; For gravity and buoyancy The loads acting on the riser unit can be expressed as follows: ， ， Generalized force vector Including the system's generalized elastic force and system generalized external forces , ， in, Including the tangential component of ocean current force Ocean current force normal component Seabed support ,gravity and buoyancy .

8. The nonlinear dynamics and fatigue damage analysis method for a gently sloping riser according to claim 7, characterized in that, Step (7) specifically includes the following steps: Taking the mass matrix of a unit as an example, after the unit assembly is completed, the mass matrix of the system composed of N units is... Represented as: ， in, This represents the unit mass matrix of the Nth unit in the riser system; The unit global position vector of N units Assembled into system generalized coordinates : The dynamic equations of the riser system can be expressed in the following matrix form: ， in, The Jacobian matrix (i.e., the first-order partial derivative matrix) of the complete constraint equations with respect to the system's generalized coordinates q describes the sensitivity of the constraints to changes in the system's configuration. This is the transpose of the Jacobian matrix, used to project the constraint reaction force onto the generalized force direction of the system; This represents the generalized acceleration vector of the system; For each constraint, represents the Lagrange multiplier vector, and its physical meaning is the magnitude of the constraint force; Q represents the generalized force vector acting on the system. The right-hand side of the system dynamics equations contains nonlinear coupling terms caused by velocity and contributions from constraints that are explicitly time-dependent. Together, they constitute the constraint compatibility conditions at the acceleration level in numerical solutions. in, Expanding the right-hand side of the system dynamics equations corrected by Baumgarte's method, we get: ， First-order partial derivatives of the constraint equations in time for: ， Solving again yields the second-order partial derivatives. for: ， in: ， ， in, This represents the explicit partial derivative of the constraint with respect to time t, reflecting the characteristics of the constraint itself as it changes over time. It represents the second-order explicit partial derivative of the constraint with respect to time, used to describe the acceleration effect of the constraint's time-varying rate; and These are the generalized velocity vector and the generalized acceleration vector of the system, respectively. Representing vectors The partial derivative with respect to the generalized coordinate q, this term reflects the nonlinear spatial variation of velocity in the acceleration equation; The mixed partial derivative is the constraint derivative taken with respect to coordinate q and then with respect to time t. It describes the rate of change of the constraint gradient with time and corresponds to the velocity-related coupling term in the dynamic equation.

9. The nonlinear dynamics and fatigue damage analysis method for a gently sloping riser according to claim 8, characterized in that, The method for solving the system dynamics equations in step (8) is as follows: The Baumgarte solution method is used to solve the dynamic equations. In order to maintain the stability of the solution results, a feedback system can be introduced into the constraint equations in the Baumgarte solution method. ， Where α and β are Baumgarte stabilization parameters; The final dynamic equations of the riser model are expressed in the following form: ， The dynamic equations of the final riser model are solved to obtain the displacement, velocity, and acceleration response time histories of each node in the riser system.

10. The nonlinear dynamics and fatigue damage analysis method for a gently sloping riser according to claim 9, characterized in that, In step (10), fatigue damage calculation uses the SN curve and the Miner linear cumulative damage criterion: ， in, Represents the cumulative damage level; Represents the total number of load levels; Represents the current load level sequence number; This represents the actual number of cycles experienced at the i-th load level; Representative at the The critical number of cycles required to cause failure at a given load level.

Images