An axial-lateral-torsional-bending fully coupled perforating string nonlinear dynamics modeling and solving method

Abstract

This invention discloses a nonlinear dynamic modeling and solution method for fully coupled axial-lateral-torsional-bending perforated tubing, relating to the field of tubing dynamics. The method includes: determining the analysis object of the perforated tubing and its structure, materials, wellbore constraints, and detonation load parameters; establishing a global coordinate system and a local coordinate system, discretizing them using extended three-dimensional Euler-Bernoulli beam elements; calculating the kinetic energy, potential energy, and resultant external force vector of the beam elements; establishing nonlinear dynamic equations based on the second type of Lagrange equations and setting boundary conditions; calculating the initial equilibrium configuration and initial stress state, and solving the nonlinear dynamic response under detonation impact loads. This invention addresses the problems of insufficient characterization of coupled vibrations, insufficient description of initial states, and insufficient convergence stability and computational efficiency under strongly nonlinear conditions in existing perforated tubing dynamic models. It achieves accurate prediction of the fully coupled axial-lateral-torsional-bending vibration response of the perforated tubing, thereby improving oil and gas well productivity and ensuring downhole operation safety.

Description

Technical Field

[0001] This invention relates to the field of tubing dynamics modeling and numerical solution technology, specifically to a nonlinear dynamics modeling and solution method for a fully coupled axial-lateral-torsional-bending perforated tubing. Background Technology

[0002] In oil and gas well completion operations, high-speed jets generated by high-energy explosive detonation penetrate the casing, cement sheath, and part of the formation to establish a connection between the oil and gas layer and the wellbore. This process is called perforation. The quality of perforation directly affects the production capacity and operational safety of oil and gas wells. The detonation wave generated during perforation propagates within the wellbore, creating transient impact loads that induce strong vibrations in the perforation string. If there are weaknesses in the structural design or connections, a series of structural failures can easily occur, such as string instability, buckling fracture, and packer release. These failures will compromise wellbore integrity, reduce production efficiency, and even endanger the lives of on-site personnel. Figure 2 As shown, the perforated string is located in a casing-constrained environment downhole and mainly consists of multiple components connected in series, including a packer, tubing, vibration damper, safety gun, perforating gun, and connecting tools. The perforated string is simultaneously subjected to detonation impact loads and casing constraints, causing its dynamic behavior to exhibit significant nonlinear characteristics.

[0003] In existing technologies, most perforated string dynamic models employ axial-lateral-torsional coupled dynamic models. These models typically introduce cross-sectional coordinates and bending angular displacement terms into the axial displacement field and establish nonlinear coupling relationships between axial, lateral, and torsional responses through axial nonlinear strain. However, existing models usually employ only four degrees of freedom per node, failing to adequately characterize lateral rotational kinetic energy, lateral bending strain energy, and the coupling effect between bending response and axial, lateral, and torsional responses. Furthermore, the string structure is often simplified, with the initial state typically assumed to coincide with the wellbore trajectory, neglecting the initial equilibrium configuration and initial stress state under wellbore hydrostatic pressure. Meanwhile, existing studies widely employ the Newmark-β isotime integration method, solving the corresponding governing equations through a full matrix linear equation system. For strongly nonlinear contact / impact problems, this method suffers from low computational efficiency and insufficient convergence stability. These shortcomings limit the accurate description of local bending, spatial buckling, bending-torsional coupling, mechanical response, and failure risk.

[0004] Therefore, there is an urgent need to establish a nonlinear dynamic modeling and solution method for perforated tubing that can consider the fully coupled relationships of axial, lateral, torsional, and bending responses, while taking into account the initial static pressure state of the wellbore, the stability of strongly nonlinear contact / impact solutions, and computational efficiency. Summary of the Invention

[0005] This invention provides a nonlinear dynamic modeling and solution method for perforated tubing with full axial-lateral-torsional-bending coupling. This method addresses the shortcomings of existing perforated tubing dynamic models in terms of insufficient characterization of coupled vibrations, insufficient description of initial states, and insufficient convergence stability and computational efficiency of strongly nonlinear contact / impact solutions. It enables accurate prediction of the axial-lateral-torsional-bending fully coupled vibration response of perforated tubing, thereby improving oil and gas well productivity and ensuring downhole operation safety.

[0006] This invention is achieved through the following technical solution:

[0007] A nonlinear dynamic modeling and solution method for a fully coupled axial-lateral-torsional-bending perforation string includes the following steps:

[0008] S1. Determine the dynamic analysis object of the perforation string, the structural parameters of the perforation string, the material parameters of the perforation string, the wellbore constraint parameters, and the perforation detonation load parameters;

[0009] S2. Establish the global and local coordinate systems of the perforation string, and use extended three-dimensional Euler-Bernoulli beam elements to discretize the perforation string using finite element methods.

[0010] S3. Calculate the kinetic energy, potential energy, and net external force vector of the beam element;

[0011] S4. Based on the second type of Lagrange equation, establish the nonlinear dynamic equation of the perforation string and set the boundary conditions;

[0012] S5. Solve for the initial equilibrium configuration and initial stress state of the perforated string under the action of wellbore static pressure, and use the initial equilibrium configuration and initial stress state as initial conditions to solve for the nonlinear dynamic response of the perforated string under the action of detonation impact load.

[0013] To address the shortcomings of existing perforation string dynamic models in terms of coupled vibration characterization, downhole initial state description, and strong nonlinear contact / impact solution, this invention proposes a fully coupled axial-lateral-torsional-bending nonlinear dynamic modeling and solution method for perforation strings. This method first determines the dynamic analysis object, structural parameters, material parameters, wellbore constraint parameters, and perforation detonation load parameters of the perforation string, providing a parameter basis for nonlinear dynamic modeling. Then, a global coordinate system is established to describe the spatial position and overall motion state of the perforation string, and a local coordinate system is established to describe the local deformation and nodal degrees of freedom of beam elements. The perforation string is then discretized using extended three-dimensional Euler-Bernoulli beam elements. The kinetic energy, potential energy, and resultant external force vector of each discretized beam element are then calculated. Finally, the kinetic energy, potential energy, and resultant external force vector of the beam elements are substituted into the second type of Lagrange equation, and the dynamic equations are assembled and transformed to obtain the nonlinear dynamic equations of the perforation string. Finally, top fixed boundary conditions and bottom free boundary conditions are set. Subsequently, the initial equilibrium configuration and initial stress state of the perforation string are solved under the action of wellbore static pressure. The initial equilibrium configuration and initial stress state are used as the initial dynamic conditions of the perforation string. Then, the nonlinear dynamic response of the perforation string under detonation impact load is solved based on the joint solution framework.

[0014] Furthermore, in step S1, the dynamic analysis object of the perforation string is the perforation string below the packer. The perforation string below the packer includes a perforating gun, tubing, damper, safety gun, and connecting tool. Each component of the perforation string is equivalently represented by a three-dimensional elastic beam, the casing is treated as a rigid body, and the contact collision between the perforation string and the casing is treated by an elastic contact collision method.

[0015] Furthermore, in step S2, the extended three-dimensional Euler-Bernoulli beam element retains the kinematic assumptions of the Euler-Bernoulli beam, introduces Rayleigh-type rotational inertia, and considers the nonlinear coupling between axial, lateral, torsional, and bending deformations. The extended three-dimensional Euler-Bernoulli beam element includes two nodes, each node having six degrees of freedom, including axial displacement, two lateral displacements, torsional angular displacement, and two bending angular displacements. The beam element node displacement vector is:

[0016] ;

[0017] In the formula: Let the displacement vectors of the beam element nodes be the vectors of the beam elements. and They are nodes and nodes axial displacement; and They are nodes and nodes Lateral displacement in the y direction; and They are nodes and nodes Lateral displacement in the z-direction; and They are nodes and nodes axial torsional angular displacement; and They are nodes and nodes Angular displacement about the y-axis; and They are nodes and nodes Angular displacement about the z-axis; the superscript T indicates the transpose of the vector; This refers to the axial coordinate along the beam element axis in the local coordinate system of the beam element. and These are the two transverse coordinate directions perpendicular to the beam element axis; the beam element displacement interpolation relationship is:

[0018] ;

[0019] In the formula: Let be the displacement vector of any point on the axis of the beam element; The shape function matrix of the beam element includes a one-dimensional linear Lagrange shape function for axial displacement and axial torsional angular displacement, and a cubic Hermite shape function for lateral bending displacement. , , , , and These represent the axial displacement, lateral displacement in the y-direction, lateral displacement in the z-direction, axial torsional angular displacement, bending angular displacement about the y-axis, and bending angular displacement about the z-axis at any point on the beam element axis, respectively. These are the local dimensionless coordinates or natural coordinates of the beam element.

[0020] Furthermore, in step S3, the kinetic energy of the beam element is calculated using the following formula:

[0021] ;

[0022] In the formula, For the kinetic energy of the beam element; The translational kinetic energy of the beam element; The rotational kinetic energy of a beam element containing Rayleigh-type rotational inertia;

[0023] The translational kinetic energy of the beam element is calculated using the following formula:

[0024] ;

[0025] In the formula: The translational kinetic energy of the beam element; The translational velocity vector of the beam element; For beam element density; For the volume of a beam element; θ is the initial angle between the beam element and the positive x-axis; r is the eccentricity of the beam element; The length of the beam element; The cross-sectional area of ​​the beam element; The axial torsional angular velocity of the beam element; , and These are the axial velocities of the beam element, Lateral velocity in the direction and Lateral velocity in the direction of;

[0026] The rotational kinetic energy, including Rayleigh-type moment of inertia, is calculated using the following formula:

[0027] ;

[0028] In the formula: The rotational kinetic energy of a beam element containing Rayleigh-type rotational inertia; It is the angular velocity vector; , and These are the axial torsional angular velocity, the bending angular velocity about the y-axis, and the bending angular velocity about the z-axis, respectively. The inertia matrix of the beam element; , and They are respectively about , and Moment of inertia of the axis;

[0029] The potential energy of the beam element includes axial strain energy, transverse bending strain energy, torsional strain energy, and nonlinear coupled strain energy between axial, transverse, torsional, and bending responses. The potential energy of the beam element is:

[0030] ;

[0031] In the formula: , and These are the potential energy of the beam element, the potential energy of normal strain, and the potential energy of shear strain, respectively. The elastic modulus of the beam element; and The cross sections of the beam elements are respectively The second moment of cross-section and the polar moment of inertia of the shaft; Let be the shear modulus of the beam element.

[0032] Furthermore, in step S3, the resultant external force vector of the beam element includes the beam element gravity vector, the perforation fluid buoyancy vector, the pipe-casing contact collision force vector, and the detonation impact load vector. The construction of the resultant external force vector of the beam element includes the following steps:

[0033] S301, Based on the density, cross-sectional area, length, gravitational acceleration, and local coordinate system of the beam element. The angle between the axis and the direction of gravity is used to determine the gravity distribution load of the beam element, and the gravity distribution load of the beam element is equivalent to the gravity vector of the beam element based on the principle of virtual work.

[0034] S302, based on the density of the perforating fluid, gravitational acceleration, beam element length, beam element cross-sectional area, and the local coordinate system of the beam element. The angle between the axis and the direction of gravity is used to determine the buoyancy load distributed by the perforation fluid on the beam element, and the buoyancy load distributed by the perforation fluid is equivalent to the buoyancy vector of the perforation fluid based on the principle of virtual work.

[0035] S303. Based on the lateral displacement of the pipe column node, the outer radius of the pipe column and the inner radius of the casing, the contact state between the pipe column and the casing is determined, and the normal contact force is calculated using the penalty function method, and the tangential friction force is calculated using the Coulomb friction law to obtain the contact collision force vector between the pipe column and the casing.

[0036] S304. Based on the perforation projectile charge, the distance from the center of the perforation projectile charge to the reference point, the initial reference distance, and the shock wave propagation velocity, determine the detonation pressure-time relationship; based on the detonation pressure-time relationship, the perforation channel radius, the outer diameter of two adjacent tubular components at the variable cross-section connection, and the phase angle of the perforation projectile, calculate the axial detonation impact load and the transverse detonation impact load, and decompose the transverse detonation impact load to two transverse directions of the local coordinate system to obtain the detonation impact load vector.

[0037] S305. The gravity vector of the beam element, the buoyancy vector of the perforation fluid, the contact collision force vector of the pipe-casing and the detonation impact load vector are superimposed to obtain the resultant external force vector of the beam element.

[0038] Furthermore, the beam element gravity vector, perforation fluid buoyancy vector, tubing-casing contact collision force vector, detonation impact load vector, and beam element resultant external force vector are calculated as follows:

[0039] The gravity vector of the beam element is:

[0040] ;

[0041] In the formula: Let g be the gravity vector of the beam element; This represents the weight per unit length of the beam element. The length of the beam element; Beam element axis The angle between the axis and the direction of gravity;

[0042] The buoyancy vector of the perforation fluid is:

[0043] ;

[0044] In the formula: The vector of buoyancy of the perforating fluid; Density of the perforating fluid; For beam element density;

[0045] The contact collision force vector of the tubing-casing is:

[0046] ;

[0047] In the formula: This represents the collision force vector at the contact point between the tubing and the casing. and They are nodes Normal contact force and tangential friction force Component of direction; and They are nodes Normal contact force and tangential friction force Component of direction; and They are nodes Normal contact force and tangential friction force Component of direction; and They are nodes Normal contact force and tangential friction force Component of direction; and They are nodes and nodes The frictional torque around the axis of the tube column; and These are the two transverse coordinate directions perpendicular to the beam element axis;

[0048] The detonation impact load vector is:

[0049] ;

[0050] In the formula: This is the detonation impact load vector; and They are nodes and nodes Axial detonation impact load; and They are nodes and nodes exist Lateral detonation impact load in the direction; and They are nodes and nodes exist Lateral detonation impact load in the direction;

[0051] The resultant external force vector of the beam element is:

[0052] ;

[0053] In the formula: Let be the resultant external force vector of the beam element.

[0054] Furthermore, the detonation impact load vector and the tubing-casing contact collision force vector are determined as follows: the detonation pressure-time relationship is:

[0055] ;

[0056] In the formula: For detonation pressure; The duration of the detonation impact load; This refers to the propagation time of the detonation wave; This refers to the propellant charge of the perforation projectile; This is the distance from the center of the perforation projectile's charge to the reference point; This is the initial reference distance; Let be the propagation speed of the shock wave; e is the natural constant.

[0057] The transverse detonation impact load in the local coordinate system direction and The components in the directional direction and the axial detonation impact load are:

[0058] ;

[0059] In the formula: and These are the transverse detonation impact loads in the local coordinate system. direction and Component of direction; For axial detonation impact load; The radius of the perforation channel; and These are the outer diameters of two adjacent tubular components at the variable cross-section connection; The phase angle of the perforating projectile;

[0060] In the pipe-casing contact collision model, the normal contact force is:

[0061] ;

[0062] In the formula: Normal contact force; This is the radial distance between the center of the tubing string and the center of the casing. The outer radius of the tubular column; The inner radius of the sleeve; The penalty function is the contact stiffness coefficient;

[0063] The normal contact force in the local coordinate system direction and The components in the direction are:

[0064] ;

[0065] In the formula: Normal contact force and local coordinate system The included angle of the axis;

[0066] The tangential friction force in the local coordinate system direction and Components in the direction and frictional torque around the axis of the tube column for:

[0067] ;

[0068] In the formula: Tangential friction force and local coordinate system The included angle of the axis; and This is a parameter for correcting the friction coefficient. The coefficient of static friction; The coefficient of kinetic friction; ω is the angular velocity about the axis of the tube column.

[0069] Furthermore, in step S4, the nonlinear dynamic equations of the perforation string are established based on the second type of Lagrange equations as follows:

[0070] ;

[0071] In the formula: The linear stiffness matrix of the perforation string in the global coordinate system; , and These are the three nonlinear stiffness matrices of the perforated string in the global coordinate system; , and These are the displacement vector, velocity vector, and acceleration vector in the global coordinate system, respectively. and These are the translational mass matrix and rotational mass matrix of the perforation column in the global coordinate system, respectively. The damping matrix of the perforation string in the global coordinate system; The resultant external force vector of the perforation column in the global coordinate system;

[0072] The Rayleigh damping matrix is:

[0073] ;

[0074] In the formula: This is the mass damping coefficient; This is the stiffness damping coefficient.

[0075] Furthermore, the boundary conditions include a top fixed boundary condition and a bottom free boundary condition, wherein the top fixed boundary condition is:

[0076] ;

[0077] In the formula: Here is the global displacement vector for the top boundary nodes, which correspond to the fixed constraint of the packer; the bottom free boundary conditions are:

[0078] ;

[0079] In the formula: The shear force vector at the bottom boundary node; Here is the bending moment vector at the bottom boundary node, which corresponds to the free-hanging end unaffected by bending moment and shear force; the initial static equilibrium equation for the wellbore is:

[0080] ;

[0081] In the formula: and These are the resultant external force and the initial equilibrium displacement vector under wellbore static pressure, respectively; the initial equilibrium configuration and initial stress state of the perforated string under wellbore static pressure are determined based on the initial equilibrium displacement vector, and the initial equilibrium configuration and initial stress state serve as the initial conditions for the dynamic analysis of detonation impact load;

[0082] The solution to the nonlinear dynamic response of a perforated string under detonation impact load includes the following steps:

[0083] S501, set the time step, total computation time, NOCH-α time integration parameters, Newton-Raphson nonlinear iteration tolerance, contact constraint tolerance, and maximum number of iterations;

[0084] S502. Based on the displacement, velocity, and acceleration of the previous time step, the initial displacement, velocity, and acceleration of the current time step are predicted using the NOCH-α time integration method.

[0085] S503. Substitute the predicted displacement, velocity and acceleration of the current time step into the nonlinear dynamic equation of the perforation string to construct the nonlinear dynamic equilibrium residual vector of the current time step.

[0086] S504. Based on the nonlinear dynamic equilibrium residual vector, construct the Jacobian matrix and displacement correction equation for the Newton-Raphson nonlinear iteration;

[0087] S505. Update the pipe-casing contact state according to the current displacement state, determine whether the penetration amount of the contact node meets the contact constraint, and correct the penalty function contact stiffness when the contact constraint is not met.

[0088] S506. Based on the sparse band structure of the Jacobian matrix, the Newton-Raphson displacement correction equation is solved by the band Gaussian elimination method and the back substitution method to obtain the displacement correction vector of the current iteration step.

[0089] S507. Determine whether the current time step has converged based on the contact constraint conditions and the relative displacement increment convergence criterion.

[0090] S508. When the current time step meets the convergence condition, update the displacement, velocity, and acceleration of the current time step and proceed to the next time step; when the current time step does not meet the convergence condition, return to step S504 to continue the Newton-Raphson nonlinear iteration.

[0091] Furthermore, in steps S502 to S507, the NOCH-α time integration method, Newton-Raphson nonlinear iteration, contact constraint checking, relative displacement increment convergence judgment, and zonal matrix solution are performed as follows:

[0092] ;

[0093] In the formula: , , The first Displacement vector of the perforation string at each time step; velocity vector of the perforation string; acceleration vector of the perforation string. , , The first Displacement vector of the perforation string at each time step; velocity vector of the perforation string; acceleration vector of the perforation string. For time step; , , , , , , , , and All are NOCH— The method's integration parameter;

[0094] The Newton-Raphson displacement correction formula is:

[0095] ;

[0096] In the formula: and The first The first time step Displacement vector and displacement correction vector in the second Newton-Raphson nonlinear iteration; For the first The first time step The corrected displacement vector in the next Newton-Raphson nonlinear iteration; For the first The first time step Jacobian matrix in Newton-Raphson nonlinear iteration; For the first The first time step Nonlinear dynamic equilibrium residual vector in the second Newton-Raphson nonlinear iteration;

[0097] The nonlinear dynamic equilibrium residual vector is:

[0098] ;

[0099] In the formula: For the first The first time step The resultant external force vector of the next Newton-Raphson nonlinear iteration step; For the perforation string in the global coordinate system, the first... The first time step The linear stiffness matrix of the next Newton-Raphson nonlinear iteration step; , and For the perforation string in the global coordinate system, the first... The first time step The three nonlinear stiffness matrices of the next Newton-Raphson nonlinear iteration step; and These are the translational mass matrix and rotational mass matrix of the perforation string in the global coordinate system, respectively. For the perforation string in the global coordinate system, the first... The first time step Damping matrix for the next Newton-Raphson nonlinear iteration step;

[0100] The Jacobian matrix is:

[0101] ;

[0102] In steps S505 and S507, the contact constraint check formula is:

[0103] ;

[0104] In the formula: For the first The number of nodes penetrating each contact node; This represents the total number of nodes in the perforation string; To allow for contact penetration tolerance; when Less than or equal to When the contact constraint is satisfied, it is determined that the contact constraint is satisfied; when Greater than If the contact constraint is not satisfied, the penalty function contact stiffness correction is initiated; the adaptive penalty stiffness update formula is:

[0105] ;

[0106] In the formula: For the current number The penalty function contact stiffness coefficient used in the Newton-Raphson nonlinear iteration; The corrected penalty function contact stiffness coefficient; This represents the upper limit of the contact stiffness coefficient of the penalty function; For the first The set of contact nodes that do not satisfy the contact constraints during the second contact constraint check; This is the penalty stiffness amplification factor; Number the contact nodes;

[0107] The formula for calculating the displacement correction vector in the Newton-Raphson nonlinear iteration based on banded Gaussian elimination and back substitution is as follows:

[0108] ;

[0109] In the formula: The displacement correction vector is the first... One degree of freedom component; For the first The Jacobian matrix after second-order banded Gaussian elimination is the first... line, number Column elements; For the first The Jacobian matrix after second-order banded Gaussian elimination is the first... line, number The main diagonal element of the column; The first one obtained by back substitution One degree of freedom displacement correction component; For the first The time step, the first In the Newton-Raphson nonlinear iteration, after the th The right-hand residual vector after the second-order banded Gaussian elimination is the first... One degree of freedom component; This represents the total number of global degrees of freedom. This is the current degree of freedom number in the back-substitution solution process;

[0110] The current time step displacement correction vector is:

[0111] ;

[0112] The convergence criterion for relative displacement increment is:

[0113] ;

[0114] In the formula: The displacement convergence tolerance is defined; when the relative displacement increment is not greater than the displacement convergence tolerance, the displacement is determined to be converged; when the relative displacement increment is greater than the displacement convergence tolerance, the displacement is determined to be non-converged, and the process returns to step S504 to continue iteration.

[0115] Compared with the prior art, the present invention has at least the following advantages and beneficial effects:

[0116] (1) The present invention uses extended three-dimensional Euler-Bernoulli beam elements and introduces Rayleigh type rotational inertia, which can describe the nonlinear coupling relationship between axial, transverse, torsional and bending responses;

[0117] (2) The present invention uses the initial equilibrium configuration and initial stress state under the action of wellbore static pressure as the initial conditions for the dynamic analysis of detonation impact load, so that the nonlinear dynamic model of the perforation string is more consistent with the actual downhole loading state.

[0118] (3) The present invention comprehensively considers detonation impact load, casing constraint and tubing-casing contact collision, and can obtain the nonlinear dynamic response of the perforated tubing under complex boundary constraints and transient loads;

[0119] (4) This invention combines the NOCH-α time integration method, Newton-Raphson nonlinear iteration, penalty function contact processing method, and banded Gaussian elimination and back substitution method, which improves the solution stability and computational efficiency of strongly nonlinear contact / impact problems. Attached Figure Description

[0120] The accompanying drawings, which are included to provide a further understanding of embodiments of the invention and form part of this application, do not constitute a limitation thereof. In the drawings:

[0121] Figure 1 This is a flowchart illustrating a specific embodiment of the present invention;

[0122] Figure 2 This is a schematic diagram of the perforation string structure in a specific embodiment of the present invention;

[0123] Figure 3 This is a comparison chart of the detonation pressure at the upper pressure gauge support position and the measured data in a specific embodiment of the present invention;

[0124] Figure 4 This is a comparison diagram of the axial acceleration of the upper pressure gauge support position and the measured data in a specific embodiment of the present invention;

[0125] Figure 5 This is a von Mises stress cloud diagram of the perforated string in a specific embodiment of the present invention;

[0126] Figure 6 This is a graph showing the maximum von Mises stress of the perforated string in a specific embodiment of the present invention. Detailed Implementation

[0127] To make the objectives, technical solutions, and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the embodiments and accompanying drawings. The illustrative embodiments and descriptions of the present invention are only used to explain the present invention and are not intended to limit the present invention.

[0128] Example 1: A nonlinear dynamic modeling and solution method for a fully coupled axial-lateral-torsional-bending perforation string, based on the following assumptions: each component of the perforation string is represented by a three-dimensional elastic beam, the casing is treated as a rigid body, and the contact collision between the perforation string and the casing is treated as an elastic contact collision.

[0129] like Figure 1 As shown, the method includes the following steps:

[0130] S1. Determine the dynamic analysis object of the perforation string, the structural parameters of the perforation string, the material parameters of the perforation string, the wellbore constraint parameters, and the perforation detonation load parameters;

[0131] S2. Establish the global and local coordinate systems of the perforation string, and use extended three-dimensional Euler-Bernoulli beam elements to discretize the perforation string using finite element methods.

[0132] S3. Calculate the kinetic energy, potential energy, and net external force vector of the beam element;

[0133] S4. Based on the second type of Lagrange equation, establish the nonlinear dynamic equation of the perforation string and set the boundary conditions;

[0134] S5. Solve for the initial equilibrium configuration and initial stress state of the perforated string under the action of wellbore static pressure, and use the initial equilibrium configuration and initial stress state as initial conditions to solve for the nonlinear dynamic response of the perforated string under the action of detonation impact load.

[0135] In step S1, the dynamic analysis object of the perforation string is the perforation string below the packer. The perforation string below the packer includes a perforating gun, tubing, damper, safety gun, and connecting tool. Each component of the perforation string is equivalently represented by a three-dimensional elastic beam, and the casing is treated as a rigid body. The contact collision between the perforation string and the casing is treated by an elastic contact collision method.

[0136] In step S2, the extended three-dimensional Euler-Bernoulli beam element retains the kinematic assumptions of the Euler-Bernoulli beam, introduces Rayleigh-type rotational inertia, and considers the nonlinear coupling between axial, lateral, torsional, and bending deformations. The extended three-dimensional Euler-Bernoulli beam element includes two nodes, each node having six degrees of freedom, including axial displacement, two lateral displacements, torsional angular displacement, and two bending angular displacements. The beam element nodal displacement vector is:

[0137] ;

[0138] In the formula: Let the displacement vectors of the beam element nodes be the vectors of the beam elements. and They are nodes and nodes axial displacement; and They are nodes and nodes Lateral displacement in the y direction; and They are nodes and nodes Lateral displacement in the z-direction; and They are nodes and nodes axial torsional angular displacement; and They are nodes and nodes Angular displacement about the y-axis; and They are nodes and nodes Angular displacement about the z-axis; This refers to the axial coordinate along the beam element axis in the local coordinate system of the beam element. and These are the two transverse coordinate directions perpendicular to the beam element axis; the superscript T indicates the transpose of the vector; the beam element displacement interpolation relationship is:

[0139] ;

[0140] In the formula: Let be the displacement vector of any point on the axis of the beam element; The shape function matrix of the beam element includes a one-dimensional linear Lagrange shape function for axial displacement and axial torsional angular displacement, and a cubic Hermite shape function for lateral bending displacement. , , , , and These represent the axial displacement, lateral displacement in the y-direction, lateral displacement in the z-direction, axial torsional angular displacement, bending angular displacement about the y-axis, and bending angular displacement about the z-axis at any point on the beam element axis, respectively. These are the local dimensionless coordinates or natural coordinates of the beam element.

[0141] In step S3, the kinetic energy of the beam element is calculated using the following formula:

[0142] ;

[0143] In the formula, For the kinetic energy of the beam element; The translational kinetic energy of the beam element; The rotational kinetic energy of a beam element containing Rayleigh-type rotational inertia;

[0144] Among them, the translational kinetic energy of the beam element Calculated using the following formula:

[0145] ;

[0146] In the formula: The translational kinetic energy of the beam element; The translational velocity vector of the beam element; For beam element density; For the volume of a beam element; θ is the initial angle between the beam element and the positive x-axis; r is the eccentricity of the beam element; The length of the beam element; The cross-sectional area of ​​the beam element; The axial torsional angular velocity of the beam element; , and These are the axial velocities of the beam element, Lateral velocity in the direction and Lateral velocity in the direction of;

[0147] The rotational kinetic energy, including Rayleigh-type moment of inertia, is calculated using the following formula:

[0148] ;

[0149] In the formula: The rotational kinetic energy of a beam element containing Rayleigh-type rotational inertia; It is the angular velocity vector; , and These are the axial torsional angular velocity, the bending angular velocity about the y-axis, and the bending angular velocity about the z-axis, respectively. The inertia matrix of the beam element; , and They are respectively about , and Moment of inertia of the axis;

[0150] The potential energy of the beam element includes axial strain energy, transverse bending strain energy, torsional strain energy, and nonlinear coupled strain energy between axial, transverse, torsional, and bending responses. The potential energy of the beam element is:

[0151] ;

[0152] In the formula: , and These are the potential energy of the beam element, the potential energy of normal strain, and the potential energy of shear strain, respectively. The elastic modulus of the beam element; and The cross sections of the beam elements are respectively The second moment of cross-section and the polar moment of inertia of the shaft; Let be the shear modulus of the beam element.

[0153] In step S3, the resultant external force vector of the beam element includes the beam element gravity vector, the perforation fluid buoyancy vector, the pipe-casing contact collision force vector, and the detonation impact load vector. The construction of the resultant external force vector of the beam element includes the following steps:

[0154] S301, Based on the density, cross-sectional area, length, gravitational acceleration, and local coordinate system of the beam element. The angle between the axis and the direction of gravity is used to determine the gravity distribution load of the beam element, and the gravity distribution load of the beam element is equivalent to the gravity vector of the beam element based on the principle of virtual work.

[0155] S302, based on the density of the perforating fluid, gravitational acceleration, beam element length, beam element cross-sectional area, and the local coordinate system of the beam element. The angle between the axis and the direction of gravity is used to determine the buoyancy load distributed by the perforation fluid on the beam element, and the buoyancy load distributed by the perforation fluid is equivalent to the buoyancy vector of the perforation fluid based on the principle of virtual work.

[0156] S303. Based on the lateral displacement of the pipe column node, the outer radius of the pipe column and the inner radius of the casing, the contact state between the pipe column and the casing is determined, and the normal contact force is calculated using the penalty function method, and the tangential friction force is calculated using the Coulomb friction law to obtain the contact collision force vector between the pipe column and the casing.

[0157] S304. Based on the perforation projectile charge, the distance from the center of the perforation projectile charge to the reference point, the initial reference distance, and the shock wave propagation velocity, determine the detonation pressure-time relationship; based on the detonation pressure-time relationship, the perforation channel radius, the outer diameter of two adjacent tubular components at the variable cross-section connection, and the phase angle of the perforation projectile, calculate the axial detonation impact load and the transverse detonation impact load, and decompose the transverse detonation impact load to two transverse directions of the local coordinate system to obtain the detonation impact load vector.

[0158] S305. The gravity vector of the beam element, the buoyancy vector of the perforation fluid, the contact collision force vector of the pipe-casing and the detonation impact load vector are superimposed to obtain the resultant external force vector of the beam element.

[0159] The beam element gravity vector, perforation fluid buoyancy vector, tubing-casing contact collision force vector, detonation impact load vector, and beam element resultant external force vector are calculated as follows:

[0160] The gravity vector of the beam element is:

[0161] ;

[0162] In the formula: Let g be the gravity vector of the beam element; This represents the weight per unit length of the beam element. The length of the beam element; Beam element axis The angle between the axis and the direction of gravity;

[0163] The buoyancy vector of the perforation fluid is:

[0164] ;

[0165] In the formula: The vector of buoyancy of the perforating fluid; Density of the perforating fluid; For beam element density;

[0166] The contact collision force vector of the tubing-casing is:

[0167] ;

[0168] In the formula: This represents the collision force vector at the contact point between the tubing and the casing. and They are nodes Normal contact force and tangential friction force Component of direction; and They are nodes Normal contact force and tangential friction force Component of direction; and They are nodes Normal contact force and tangential friction force Component of direction; and They are nodes Normal contact force and tangential friction force Component of direction; and They are nodes and nodes The frictional torque around the axis of the tube column; and These are the two transverse coordinate directions perpendicular to the beam element axis;

[0169] The detonation impact load vector is:

[0170] ;

[0171] In the formula: This is the detonation impact load vector; and They are nodes and nodes Axial detonation impact load; and They are nodes and nodes exist Lateral detonation impact load in the direction; and They are nodes and nodes exist Lateral detonation impact load in the direction;

[0172] The resultant external force vector of the beam element is:

[0173] ;

[0174] In the formula: Let be the resultant external force vector of the beam element.

[0175] The detonation impact load vector and the tubing-casing contact collision force vector are determined as follows: the detonation pressure-time relationship is:

[0176] ;

[0177] In the formula: For detonation pressure; The duration of the detonation impact load; This refers to the propagation time of the detonation wave; This refers to the propellant charge of the perforation projectile; This is the distance from the center of the perforation projectile's charge to the reference point; This is the initial reference distance; Let be the propagation speed of the shock wave; e is the natural constant.

[0178] The transverse detonation impact load in the local coordinate system direction and The components in the directional direction and the axial detonation impact load are:

[0179] ;

[0180] In the formula: and These are the transverse detonation impact loads in the local coordinate system. direction and Component of direction; For axial detonation impact load; The radius of the perforation channel; and These are the outer diameters of two adjacent tubular components at the variable cross-section connection; The phase angle of the perforating projectile;

[0181] In the pipe-casing contact collision model, the normal contact force is:

[0182] ;

[0183] In the formula: Normal contact force; This is the radial distance between the center of the tubing string and the center of the casing. The outer radius of the tubular column; The inner radius of the sleeve; The penalty function is the contact stiffness coefficient;

[0184] The normal contact force in the local coordinate system direction and The components in the direction are:

[0185] ;

[0186] In the formula: Normal contact force and local coordinate system The included angle of the axis;

[0187] The tangential friction force in the local coordinate system direction and Components in the direction and frictional torque around the axis of the tube column for:

[0188] ;

[0189] In the formula: Tangential friction force and local coordinate system The included angle of the axis; and This is a parameter for correcting the friction coefficient. The coefficient of static friction; The coefficient of kinetic friction; ω is the angular velocity about the axis of the tube column.

[0190] In step S4, the nonlinear dynamic equations of the perforation string are established based on the second type of Lagrange equations as follows:

[0191] ;

[0192] In the formula: The linear stiffness matrix of the perforation string in the global coordinate system; , and These are the three nonlinear stiffness matrices of the perforated string in the global coordinate system; , and These are the displacement vector, velocity vector, and acceleration vector in the global coordinate system, respectively. and These are the translational mass matrix and rotational mass matrix of the perforation column in the global coordinate system, respectively. The damping matrix of the perforation string in the global coordinate system; The resultant external force vector of the perforation column in the global coordinate system;

[0193] The Rayleigh damping matrix is:

[0194] ;

[0195] In the formula: This is the mass damping coefficient; This is the stiffness damping coefficient.

[0196] The boundary conditions include a fixed top boundary condition and a free bottom boundary condition. The fixed top boundary condition is as follows:

[0197] ;

[0198] In the formula: Here is the global displacement vector for the top boundary nodes, which correspond to the fixed constraint of the packer; the bottom free boundary conditions are:

[0199] ;

[0200] In the formula: The shear force vector at the bottom boundary node; Here is the bending moment vector at the bottom boundary node, which corresponds to the free-hanging end unaffected by bending moment and shear force; the initial static equilibrium equation for the wellbore is:

[0201] ;

[0202] In the formula: and These are the resultant external force and the initial equilibrium displacement vector under wellbore static pressure, respectively; the initial equilibrium configuration and initial stress state of the perforated string under wellbore static pressure are determined based on the initial equilibrium displacement vector, and the initial equilibrium configuration and initial stress state serve as the initial conditions for the dynamic analysis of detonation impact load;

[0203] The solution to the nonlinear dynamic response of a perforated string under detonation impact load includes the following steps:

[0204] S501, set the time step, total computation time, NOCH-α time integration parameters, Newton-Raphson nonlinear iteration tolerance, contact constraint tolerance, and maximum number of iterations;

[0205] S502. Based on the displacement, velocity, and acceleration of the previous time step, the initial displacement, velocity, and acceleration of the current time step are predicted using the NOCH-α time integration method.

[0206] S503. Substitute the predicted displacement, velocity and acceleration of the current time step into the nonlinear dynamic equation of the perforation string to construct the nonlinear dynamic equilibrium residual vector of the current time step.

[0207] S504. Based on the nonlinear dynamic equilibrium residual vector, construct the Jacobian matrix and displacement correction equation for the Newton-Raphson nonlinear iteration;

[0208] S505. Update the pipe-casing contact state according to the current displacement state, determine whether the penetration amount of the contact node meets the contact constraint, and correct the penalty function contact stiffness when the contact constraint is not met.

[0209] S506. Based on the sparse band structure of the Jacobian matrix, the Newton-Raphson displacement correction equation is solved by the band Gaussian elimination method and the back substitution method to obtain the displacement correction vector of the current iteration step.

[0210] S507. Determine whether the current time step has converged based on the contact constraint conditions and the relative displacement increment convergence criterion.

[0211] S508. When the current time step meets the convergence condition, update the displacement, velocity, and acceleration of the current time step and proceed to the next time step; when the current time step does not meet the convergence condition, return to step S504 to continue the Newton-Raphson nonlinear iteration.

[0212] In steps S502 to S507, the NOCH-α time integration method, Newton-Raphson nonlinear iteration, contact constraint checking, relative displacement increment convergence judgment, and strip matrix solution are performed as follows:

[0213] ;

[0214] In the formula: , , The first Displacement vector of the perforation string at each time step; velocity vector of the perforation string; acceleration vector of the perforation string. , , The first Displacement vector of the perforation string at each time step; velocity vector of the perforation string; acceleration vector of the perforation string. For time step; , , , , , , , , and All are NOCH— The method's integration parameter;

[0215] The Newton-Raphson displacement correction formula is:

[0216] ;

[0217] In the formula: and The first The first time step Displacement vector and displacement correction vector in the second Newton-Raphson nonlinear iteration; For the first The first time step The corrected displacement vector in the next Newton-Raphson nonlinear iteration; For the first The first time step Jacobian matrix in Newton-Raphson nonlinear iteration; For the first The first time step Nonlinear dynamic equilibrium residual vector in the second Newton-Raphson nonlinear iteration;

[0218] The nonlinear dynamic equilibrium residual vector is:

[0219] ;

[0220] In the formula: For the first The first time step The resultant external force vector of the next Newton-Raphson nonlinear iteration step; For the perforation string in the global coordinate system, the first... The first time step The linear stiffness matrix of the next Newton-Raphson nonlinear iteration step; , and For the perforation string in the global coordinate system, the first... The first time step The three nonlinear stiffness matrices of the next Newton-Raphson nonlinear iteration step; and These are the translational mass matrix and rotational mass matrix of the perforation string in the global coordinate system, respectively. For the perforation string in the global coordinate system, the first... The first time step Damping matrix for the next Newton-Raphson nonlinear iteration step;

[0221] The Jacobian matrix is:

[0222] ;

[0223] In steps S505 and S507, the contact constraint check formula is:

[0224] ;

[0225] In the formula: For the first The number of nodes penetrating each contact node; This represents the total number of nodes in the perforation string; To allow for contact penetration tolerance; when Less than or equal to When the contact constraint is satisfied, it is determined that the contact constraint is satisfied; when Greater than If the contact constraint is not satisfied, the penalty function contact stiffness correction is initiated; the adaptive penalty stiffness update formula is:

[0226] ;

[0227] In the formula: For the current number The penalty function contact stiffness coefficient used in the Newton-Raphson nonlinear iteration; The corrected penalty function contact stiffness coefficient; This represents the upper limit of the contact stiffness coefficient of the penalty function; For the first The set of contact nodes that do not satisfy the contact constraints during the second contact constraint check; This is the penalty stiffness amplification factor; Number the contact nodes;

[0228] The formula for calculating the displacement correction vector in the Newton-Raphson nonlinear iteration based on banded Gaussian elimination and back substitution is as follows:

[0229] ;

[0230] In the formula: The displacement correction vector is the first... One degree of freedom component; For the first The Jacobian matrix after second-order banded Gaussian elimination is the first... line, number Column elements; For the first The Jacobian matrix after second-order banded Gaussian elimination is the first... line, number The main diagonal element of the column; The first one obtained by back substitution One degree of freedom displacement correction component; For the first The time step, the first In the Newton-Raphson nonlinear iteration, after the th The right-hand residual vector after the second-order banded Gaussian elimination is the first... One degree of freedom component; This represents the total number of global degrees of freedom. This is the current degree of freedom number in the back-substitution solution process;

[0231] The current time step displacement correction vector is:

[0232] ;

[0233] The convergence criterion for relative displacement increment is:

[0234] ;

[0235] In the formula: The displacement convergence tolerance is defined; when the relative displacement increment is not greater than the displacement convergence tolerance, the displacement is determined to be converged; when the relative displacement increment is greater than the displacement convergence tolerance, the displacement is determined to be non-converged, and the process returns to step S504 to continue iteration.

[0236] Example 2:

[0237] A nonlinear dynamic modeling and solution method for a fully coupled axial-lateral-torsional-bending perforation string is proposed. Based on Example 1, this example is verified using an oil well in an oilfield.

[0238] A schematic diagram of the perforation string structure is shown below. Figure 2 As shown, the packer setting position is 3000m. The packer and subsequent perforation string structure includes, in sequence, the packer, tubing, upper screen tube, upper pressure gauge support, bidirectional vibration damper, lower pressure gauge support, lower screen tube, detonator, safety gun, perforating gun, and pressure transmission gun tail. The required calculation parameters are shown in Table 1.

[0239] Table 1 Model Calculation Parameters

[0240]

[0241] Using the nonlinear dynamic modeling and solution method for the fully coupled axial-lateral-torsional-bending perforation string proposed in this application, the detonation pressure and axial acceleration at the upper pressure gauge support position were calculated and compared with field measured data. The comparison results are as follows: Figure 3 and Figure 4 As shown in the figure. For detonation pressure, the calculated peak value is 63.713 MPa, corresponding to the field measured data of 67.945 MPa, with an error of 6.22%. For axial acceleration, the errors for positive and negative peak values ​​are 9.28% and 4.42%, respectively. From the overall time history variation trend and main peak characteristics, the simulated detonation pressure response and axial acceleration response are generally consistent with the measured data. The errors in detonation pressure and axial acceleration are both within 10%, which is within an acceptable range, verifying the accuracy and engineering applicability of the method in this application.

[0242] Figure 5 This is a von Mises stress contour plot of a perforated string under detonation load. Figure 6 This is the maximum von Mises stress curve of a perforated string under detonation loading. Figure 5 and Figure 6 It can be seen that the von Mises stress exhibits a clear oscillating decay trend, with a maximum value of 310.86 MPa, occurring at the tubing-packer connection (node ​​5). The evolution law and dominant frequency of the von Mises stress are generally consistent with the axial stress response, indicating that axial stress plays a dominant role in the von Mises stress, while transverse shear stress makes a secondary contribution to the local equivalent stress. Under the operating conditions of this embodiment, the maximum von Mises stress is lower than the yield strength of the tubing material (756 MPa), indicating that yield failure was not predicted in this embodiment; the tubing-packer connection is a major high-risk area and should be given special attention in structural design and strength verification.

[0243] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

[0244] It should be noted that, in this document, terms such as “comprising,” “including,” or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article, or apparatus.

Claims

1. A nonlinear dynamic modeling and solution method for a fully coupled axial-lateral-torsional-bending perforation string, characterized in that, Includes the following steps: S1. Determine the dynamic analysis object of the perforation string, the structural parameters of the perforation string, the material parameters of the perforation string, the wellbore constraint parameters, and the perforation detonation load parameters; S2. Establish the global and local coordinate systems of the perforation string, and use extended three-dimensional Euler-Bernoulli beam elements to discretize the perforation string using finite element methods. S3. Calculate the kinetic energy, potential energy, and net external force vector of the beam element; S4. Based on the second type of Lagrange equation, establish the nonlinear dynamic equation of the perforation string and set the boundary conditions; S5. Solve for the initial equilibrium configuration and initial stress state of the perforated string under the action of wellbore static pressure, and use the initial equilibrium configuration and initial stress state as initial conditions to solve for the nonlinear dynamic response of the perforated string under the action of detonation impact load.

2. The nonlinear dynamic modeling and solution method for a fully coupled axial-lateral-torsional-bending perforation string according to claim 1, characterized in that, In step S1, the dynamic analysis object of the perforation string is the perforation string below the packer. The perforation string below the packer includes a perforating gun, tubing, damper, safety gun, and connecting tool. Each component of the perforation string is equivalently represented by a three-dimensional elastic beam, and the casing is treated as a rigid body. The contact collision between the perforation string and the casing is treated by an elastic contact collision method.

3. The nonlinear dynamic modeling and solution method for a fully coupled axial-lateral-torsional-bending perforation string according to claim 1, characterized in that, In step S2, the extended three-dimensional Euler-Bernoulli beam element retains the kinematic assumptions of the Euler-Bernoulli beam, introduces Rayleigh-type rotational inertia, and considers the nonlinear coupling between axial, lateral, torsional, and bending deformations. The extended three-dimensional Euler-Bernoulli beam element includes two nodes, each node having six degrees of freedom, including axial displacement, two lateral displacements, torsional angular displacement, and two bending angular displacements. The beam element nodal displacement vector is: ； In the formula: Let the displacement vectors of the beam element nodes be the vectors of the beam elements. and They are nodes and nodes axial displacement; and They are nodes and nodes Lateral displacement in the y direction; and They are nodes and nodes Lateral displacement in the z-direction; and They are nodes and nodes axial torsional angular displacement; and They are nodes and nodes Angular displacement about the y-axis; and They are nodes and nodes Angular displacement about the z-axis; The superscript T denotes the transpose of a vector; This refers to the axial coordinate along the beam element axis in the local coordinate system of the beam element. and These are the two transverse coordinate directions perpendicular to the beam element axis; the beam element displacement interpolation relationship is: ； In the formula: Let be the displacement vector of any point on the axis of the beam element; The shape function matrix of the beam element includes a one-dimensional linear Lagrange shape function for axial displacement and axial torsional angular displacement, and a cubic Hermite shape function for lateral bending displacement. , , , , and These represent the axial displacement, lateral displacement in the y-direction, lateral displacement in the z-direction, axial torsional angular displacement, bending angular displacement about the y-axis, and bending angular displacement about the z-axis at any point on the beam element axis, respectively. These are the local dimensionless coordinates or natural coordinates of the beam element.

4. The nonlinear dynamic modeling and solution method for axial-lateral-torsional-bending fully coupled perforated tubing as described in claim 3, characterized in that, In step S3, the kinetic energy of the beam element is calculated using the following formula: ； In the formula, For the kinetic energy of the beam element; The translational kinetic energy of the beam element; The rotational kinetic energy of a beam element containing Rayleigh-type rotational inertia; The translational kinetic energy of the beam element is calculated using the following formula: ； In the formula: The translational kinetic energy of the beam element; The translational velocity vector of the beam element; For beam element density; For the volume of a beam element; θ is the initial angle between the beam element and the positive x-axis; r is the eccentricity of the beam element; The length of the beam element; The cross-sectional area of ​​the beam element; The axial torsional angular velocity of the beam element; , and These are the axial velocities of the beam element, Lateral velocity in the direction and Lateral velocity in the direction of; The rotational kinetic energy, including Rayleigh-type moment of inertia, is calculated using the following formula: ； In the formula: The rotational kinetic energy of a beam element containing Rayleigh-type rotational inertia; It is the angular velocity vector; , and These are the axial torsional angular velocity, the bending angular velocity about the y-axis, and the bending angular velocity about the z-axis, respectively. The inertia matrix of the beam element; , and They are respectively about , and Moment of inertia of the axis; The potential energy of the beam element includes axial strain energy, transverse bending strain energy, torsional strain energy, and nonlinear coupled strain energy between axial, transverse, torsional, and bending responses. The potential energy of the beam element is: ； In the formula: , and These are the potential energy of the beam element, the potential energy of normal strain, and the potential energy of shear strain, respectively. The elastic modulus of the beam element; and The cross-sections of the beam elements are respectively The second moment of cross-section and the polar moment of inertia of the shaft; This is the shear modulus of the beam element.

5. The nonlinear dynamic modeling and solution method for a fully coupled axial-lateral-torsional-bending perforation string according to claim 4, characterized in that, In step S3, the resultant external force vector of the beam element includes the beam element gravity vector, the perforation fluid buoyancy vector, the pipe-casing contact collision force vector, and the detonation impact load vector. The construction of the resultant external force vector of the beam element includes the following steps: S301, Based on the density, cross-sectional area, length, gravitational acceleration, and local coordinate system of the beam element. The angle between the axis and the direction of gravity is used to determine the gravity distribution load of the beam element, and the gravity distribution load of the beam element is equivalent to the gravity vector of the beam element based on the principle of virtual work. S302, based on the density of the perforating fluid, gravitational acceleration, beam element length, beam element cross-sectional area, and the local coordinate system of the beam element. The angle between the axis and the direction of gravity is used to determine the buoyancy load distributed by the perforation fluid on the beam element, and the buoyancy load distributed by the perforation fluid is equivalent to the buoyancy vector of the perforation fluid based on the principle of virtual work. S303. Based on the lateral displacement of the pipe column node, the outer radius of the pipe column and the inner radius of the casing, the contact state between the pipe column and the casing is determined, and the normal contact force is calculated using the penalty function method, and the tangential friction force is calculated using the Coulomb friction law to obtain the contact collision force vector between the pipe column and the casing. S304. Based on the perforation projectile charge, the distance from the center of the perforation projectile charge to the reference point, the initial reference distance, and the shock wave propagation velocity, determine the detonation pressure-time relationship; based on the detonation pressure-time relationship, the perforation channel radius, the outer diameter of two adjacent tubular components at the variable cross-section connection, and the phase angle of the perforation projectile, calculate the axial detonation impact load and the transverse detonation impact load, and decompose the transverse detonation impact load to two transverse directions of the local coordinate system to obtain the detonation impact load vector. S305. The gravity vector of the beam element, the buoyancy vector of the perforation fluid, the contact collision force vector of the pipe-casing and the detonation impact load vector are superimposed to obtain the resultant external force vector of the beam element.

6. The nonlinear dynamic modeling and solution method for a fully coupled axial-lateral-torsional-bending perforation string according to claim 5, characterized in that, The beam element gravity vector, perforation fluid buoyancy vector, tubing-casing contact collision force vector, detonation impact load vector, and beam element resultant external force vector are calculated as follows: The gravity vector of the beam element is: ； In the formula: Let g be the gravity vector of the beam element; This represents the weight per unit length of the beam element. The length of the beam element; Beam element axis The angle between the axis and the direction of gravity; The buoyancy vector of the perforation fluid is: ； In the formula: The vector of buoyancy of the perforating fluid; Density of the perforating fluid; For beam element density; The contact collision force vector of the tubing-sleeve is: ； In the formula: The tube-casing contact collision force vector; and They are nodes Normal contact force and tangential friction force Component of direction; and They are nodes Normal contact force and tangential friction force Component of direction; and They are nodes Normal contact force and tangential friction force Component of direction; and They are nodes Normal contact force and tangential friction force Component of direction; and They are nodes and nodes The frictional torque around the axis of the tube column; and These are the two transverse coordinate directions perpendicular to the beam element axis; The detonation impact load vector is: ； In the formula: This is the detonation impact load vector; and They are nodes and nodes Axial detonation impact load; and They are nodes and nodes exist Lateral detonation impact load in the direction; and They are nodes and nodes exist Lateral detonation impact load in the direction; The resultant external force vector of the beam element is: ； In the formula: Let be the resultant external force vector of the beam element.

7. The nonlinear dynamic modeling and solution method for a fully coupled axial-lateral-torsional-bending perforation string according to claim 6, characterized in that, The detonation impact load vector and the tubing-casing contact collision force vector are determined as follows: the detonation pressure-time relationship is: ； In the formula: For detonation pressure; The duration of the detonation impact load; This refers to the propagation time of the detonation wave; This refers to the propellant charge of the perforation projectile; This is the distance from the center of the perforating projectile's charge to the reference point; This is the initial reference distance; Let be the propagation speed of the shock wave; e is the natural constant. The transverse detonation impact load in the local coordinate system direction and The components in the directional direction and the axial detonation impact load are: ； In the formula: and These are the transverse detonation impact loads in the local coordinate system. direction and Component of direction; For axial detonation impact load; The radius of the perforation channel; and These are the outer diameters of two adjacent tubular components at the variable cross-section connection; The phase angle of the perforating projectile; In the pipe-casing contact collision model, the normal contact force is: ； In the formula: Normal contact force; This is the radial distance between the center of the tubing string and the center of the casing. The outer radius of the tubular column; The inner radius of the sleeve; The penalty function is the contact stiffness coefficient; The normal contact force in the local coordinate system direction and The components in the direction are: ； In the formula: Normal contact force and local coordinate system The included angle of the axis; The tangential friction force in the local coordinate system direction and Components in the direction and frictional torque around the axis of the tube column for: ； In the formula: Tangential friction force and local coordinate system The included angle of the axis; and This is a parameter for correcting the friction coefficient. The coefficient of static friction; The coefficient of kinetic friction; ω is the angular velocity about the axis of the tube column.

8. The nonlinear dynamic modeling and solution method for a fully coupled axial-lateral-torsional-bending perforation string according to claim 7, characterized in that, In step S4, the nonlinear dynamic equations of the perforation string are established based on the second type of Lagrange equations as follows: ； In the formula: The linear stiffness matrix of the perforation string in the global coordinate system; , and These are the three nonlinear stiffness matrices of the perforated string in the global coordinate system; , and These are the displacement vector, velocity vector, and acceleration vector in the global coordinate system, respectively. and These are the translational mass matrix and rotational mass matrix of the perforation column in the global coordinate system, respectively. The damping matrix of the perforation string in the global coordinate system; The resultant external force vector of the perforation column in the global coordinate system; The Rayleigh damping matrix is: ； In the formula: This is the mass damping coefficient; This represents the stiffness damping coefficient.

9. The nonlinear dynamic modeling and solution method for a fully coupled axial-lateral-torsional-bending perforation string according to claim 8, characterized in that, The boundary conditions include a fixed top boundary condition and a free bottom boundary condition. The fixed top boundary condition is as follows: ； In the formula: Here is the global displacement vector for the top boundary nodes, which correspond to the fixed constraint of the packer; the bottom free boundary conditions are: ； In the formula: The shear force vector at the bottom boundary node; Here is the bending moment vector at the bottom boundary node, which corresponds to the free-hanging end unaffected by bending moment and shear force; the initial static equilibrium equation for the wellbore is: ； In the formula: and These are the resultant external force and the initial equilibrium displacement vector under wellbore static pressure, respectively. Based on the initial equilibrium displacement vector, the initial equilibrium configuration and initial stress state of the perforated string under wellbore static pressure are determined, and the initial equilibrium configuration and initial stress state serve as the initial conditions for the dynamic analysis of detonation impact load. The solution to the nonlinear dynamic response of a perforated string under detonation impact load includes the following steps: S501, set the time step, total computation time, NOCH-α time integration parameters, Newton-Raphson nonlinear iteration tolerance, contact constraint tolerance, and maximum number of iterations; S502. Based on the displacement, velocity, and acceleration of the previous time step, the initial displacement, velocity, and acceleration of the current time step are predicted using the NOCH-α time integration method. S503. Substitute the predicted displacement, velocity and acceleration of the current time step into the nonlinear dynamic equation of the perforation string to construct the nonlinear dynamic equilibrium residual vector of the current time step. S504. Based on the nonlinear dynamic equilibrium residual vector, construct the Jacobian matrix and displacement correction equation for the Newton-Raphson nonlinear iteration; S505. Update the pipe-casing contact state according to the current displacement state, determine whether the penetration amount of the contact node meets the contact constraint, and correct the penalty function contact stiffness when the contact constraint is not met. S506. Based on the sparse band structure of the Jacobian matrix, the Newton-Raphson displacement correction equation is solved by the band Gaussian elimination method and the back substitution method to obtain the displacement correction vector of the current iteration step. S507. Determine whether the current time step has converged based on the contact constraint conditions and the relative displacement increment convergence criterion. S508. When the current time step meets the convergence condition, update the displacement, velocity, and acceleration of the current time step and proceed to the next time step; when the current time step does not meet the convergence condition, return to step S504 to continue the Newton-Raphson nonlinear iteration.

10. The nonlinear dynamic modeling and solution method for a fully coupled axial-lateral-torsional-bending perforation string according to claim 9, characterized in that, In steps S502 to S507, the NOCH-α time integration method, Newton-Raphson nonlinear iteration, contact constraint checking, relative displacement increment convergence judgment, and strip matrix solution are performed as follows: ； In the formula: , , The first Displacement vector of the perforation string at each time step; velocity vector of the perforation string; acceleration vector of the perforation string. , , The first Displacement vector of the perforation string at each time step; velocity vector of the perforation string; acceleration vector of the perforation string. For time step; , , , , , , , , and All are NOCH— The method's integration parameter; The Newton-Raphson displacement correction formula is: ； In the formula: and The first The first time step Displacement vector and displacement correction vector in the second Newton-Raphson nonlinear iteration; For the first The first time step The corrected displacement vector in the next Newton-Raphson nonlinear iteration; For the first The first time step Jacobian matrix in Newton-Raphson nonlinear iteration; For the first The first time step Nonlinear dynamic equilibrium residual vector in the second Newton-Raphson nonlinear iteration; The nonlinear dynamic equilibrium residual vector is: ； In the formula: For the first The first time step The resultant external force vector of the next Newton-Raphson nonlinear iteration step; For the perforation string in the global coordinate system, the first... The first time step The linear stiffness matrix of the next Newton-Raphson nonlinear iteration step; , and For the perforation string in the global coordinate system, the first... The first time step The three nonlinear stiffness matrices of the next Newton-Raphson nonlinear iteration step; and These are the translational mass matrix and rotational mass matrix of the perforation string in the global coordinate system, respectively. For the perforation string in the global coordinate system, the first... The first time step Damping matrix for the next Newton-Raphson nonlinear iteration step; The Jacobian matrix is: ； In steps S505 and S507, the contact constraint check formula is: ； In the formula: For the first The number of nodes penetrating each contact node; This represents the total number of nodes in the perforation string. To allow for contact penetration tolerance; when Less than or equal to When the contact constraint is satisfied, it is determined that the contact constraint is satisfied; when Greater than If the contact constraint is not satisfied, the penalty function contact stiffness correction is initiated; the adaptive penalty stiffness update formula is: ； In the formula: For the current number The penalty function contact stiffness coefficient used in the Newton-Raphson nonlinear iteration; The corrected penalty function contact stiffness coefficient; This represents the upper limit of the contact stiffness coefficient of the penalty function; For the first The set of contact nodes that do not satisfy the contact constraints during the second contact constraint check; This is the penalty stiffness amplification factor; Number the contact nodes; The formula for calculating the displacement correction vector in the Newton-Raphson nonlinear iteration based on banded Gaussian elimination and back substitution is as follows: ； In the formula: The displacement correction vector is the first... One degree of freedom component; For the first The Jacobian matrix after second-order banded Gaussian elimination is the first... line, number Column elements; For the first The Jacobian matrix after second-order banded Gaussian elimination is the first... line, number The main diagonal element of the column; The first one obtained by back substitution One degree of freedom displacement correction component; For the first The time step, the first In the Newton-Raphson nonlinear iteration, after the th The right-hand residual vector after the second-order banded Gaussian elimination is the first... One degree of freedom component; This represents the total number of global degrees of freedom. This is the current degree of freedom number in the back-substitution solution process; The current time step displacement correction vector is: ； The convergence criterion for relative displacement increment is: ； In the formula: The displacement convergence tolerance is defined; when the relative displacement increment is not greater than the displacement convergence tolerance, the displacement is determined to be converged; when the relative displacement increment is greater than the displacement convergence tolerance, the displacement is determined to be non-converged, and the process returns to step S504 to continue iteration.

Images