- Email: [email protected]
Post-buckling analysis of compressed rods in cylinders by using dynamic relaxation method
International Journal of Mechanical Sciences 159 (2019) 103–115
Contents lists available at ScienceDirect
International Journal of Mechanical Sciences journal homepage: www.elsevier.com/locate/ijmecsci
Post-buckling analysis of compressed rods in cylinders by using dynamic relaxation method Qiang Zhang a,b, Bao Jiang a, Zhongmin Xiao b,∗, Wei Cui a,∗, Jubao Liu a a b
School of Mechanical Science and Engineering, Northeast Petroleum University, Daqing 163318, PR China School of Mechanical and Aerospace Engineering, Nanyang Technological University, 50 Nanyang Avenue, 639798, Singapore
a r t i c l e
i n f o
Keywords: Post-buckling Transition section Finite element Gap element Dynamic relaxation method
a b s t r a c t Buckling of rods with cylindrical constraints is an essential problem in engineering and medical fields. Completely different to the buckling of classical Euler rods, the post-buckling modes for rods with cylindrical constraints involve highly complicated deformations and geometric configurations. In this paper, the rods are discretized into beam elements by finite element method, and the constraint relation between the rods and cylinders is described by gap element. A dynamic relaxation method for static buckling of compressed rods in cylinders is introduced. Numerical simulations have been carried out to investigate the effects of damping, element length, initial contact stiffness and penetration on post-buckling, etc. The characteristic buckling modes include deflection curves with single point, two points, three points and point-line-point contact. For the helical buckling, it is found that the transition section consists of two noncontact sections only, while the perturbed-helix section does not exist. Except the critical load and shear force of helical buckling, numerical simulation results are in good agreement with the analytical solution. By using the dynamic relaxation method, stable helical buckling modes can be obtained directly without undergoing sinusoidal buckling deformation under given compressive loads.
1. Introduction Slender rods in finite space exist in various petroleum engineering and medical applications. In petroleum engineering, the drill pipe in the borehole of several kilometers is subjected to great axial compressive force during drilling. When the axial force exceeds the threshold value, the drill pipe may bend into a three-dimensional shape, resulting in rapid increase in contact and friction force, and ultimately hinder its further drilling work, known as “lock-up” [1,2]. While in medical applications, stent therapy is often used for patients with coronary artery disease. It involves the insertion of a long flexible guide wire [3,4]. A great deal of effort is needed to drive the front end of the guide wire through a narrow lesion area with severe vascular occlusion. However, the guide wire may buckle due to large axial force. In addition, microtubules confined by microfilaments and cytoplasm around the protrusions of living cells may be helical, which may be caused by buckling [5,6]. Due to the importance of slender rods’ stability in cylinders, their buckling behavior has attracted extensive attention in the engineering field for decades. Lubinski [7] studied the two-dimensional buckling of drill pipe in wellbore for the first time and derived the calculation formula of the critical buckling load. ∗
Later, researchers extensively studied the buckling behavior of slender rods in cylinders from different aspects. For example, cylindrical constraint types include vertical [8,9], inclined and horizontal wellbores [2,10] in petroleum engineering. They also analyzed the effects of various factors on buckling modes and critical loads, including compression load [11–15], torque [16–19] gravity [9,20,21], boundary conditions [11,15,22], friction [23,24] and joints [25,26], etc. With their efforts, the following research methods have been developed: (1) theoretical methods such as beam-column model [7,26,27], buckling differential equation [17,28] and energy method [29–31]; (2) finite element methods such as implicit [9] and explicit [32–34] schemes; (3) experimental methods [31,35,36]. This paper focuses on the buckling of weightless compression rods in cylinders. Researchers have found that when the applied axial force reaches a certain critical value, the rods (restrained in cylinders) first buckle from the initial linear state into a sinusoidal shape. This phenomenon is called sinusoidal buckling. In the subsequent post-buckling process, with the increasing compressive load, more complex buckling configuration known as helical buckling occurs. Lubinski et al. [29] assumed that the tubing buckled into a space helix. Without considering the effect of boundary conditions, the mathematical relationship between the pitch and the axial load
Corresponding authors. E-mail addresses: [email protected] (Z. Xiao), [email protected] (W. Cui).
https://doi.org/10.1016/j.ijmecsci.2019.05.040 Received 9 November 2018; Received in revised form 10 May 2019; Accepted 26 May 2019 Available online 28 May 2019 0020-7403/© 2019 Published by Elsevier Ltd.
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
was derived by the energy method. The relationship is given as follows 𝐹hel =
8π2 𝐸𝐼 . 𝑝2
(1)
where Fhel is the axially critical load of helical buckling, EI is flexural rigidity of the tubing, p is the pitch of the buckled tubing. Mitchell [11] considered the boundary conditions effect and obtained another formula for the critical load. Under the clamped boundary condition, he calculated the length of the rod end to the helix section. He divided the helical buckling mode into two parts: the end transition section and the helical section. He pointed out that the end boundary conditions had an important effect on the buckling behavior. The critical load he obtained was 2.8 times that of the conventional analysis. Liu et al. [22] divided the transition section into a non-contact section and a perturbed-helix section. Sorenson et al. [12], Mitchell [36] and Huang et al. [15] divided the transition section into two no-contact sections and a perturbed-helix section. The two no-contact sections are called the first no-contact section and the second no-contact section respectively. The two no-contact sections were described by the beam-column model [37]. In the perturbed-helix section, the deformation curve was described by the buckling differential equation [22]. However, the above work considered only the final configuration of the buckled rod. They used different deflection assumptions to divide the ultimate deflection curve of the helical buckling into different sections and calculated the length of each section, but some of these length values are different, the corresponding critical loads are also different. In fact, there are some intermediate deformation stages between sinusoidal and helical shape. Based on the deflection curve hypothesis, the buckling modes and critical loads of sinusoidal to helical buckling were studied by using the small deformation hypothesis [13], the elastic model [14], the beam-column model and the differential equation [8]. In this paper, the assumption of deflection curve used in conventional theoretical analysis is abandoned. With the aid of finite element simulation, a dynamic relaxation method for post-buckling analysis of compressed rods in cylinders is proposed by combining the gap element with the beam element. We investigate the effects of damping, element size, contact stiffness and penetration on buckling. The critical conditions of various buckling modes from sinusoidal to helical buckling are analyzed. The relations between the buckling modes and shear force, bending moment, contact force and critical load are investigated. In the configuration of helical buckling, an exception that is inconsistent with the previous conclusions about the relation between the boundary condition and transition section is discovered, and a new critical load of helical buckling is obtained.
Fig. 1. Physical model of compressed rods in cylinders: (a) initial configuration of compressed rods in cylinders before buckling, (b) the deformation after buckling under the axial compression load Fc , (c) cross section before buckling, and (d) cross section after buckling. L is the rod length, 𝛽 G is the direction angle of the buckled rod with respect to global coordinates (X, Y, Z).
rigid body rotation, the rotation in the x-axis direction is restricted at the lower end. Under a given axial compressive load at the end of the rods, they will buckle and contact the inner wall of the rigid cylinders. The shape of the buckling deformation is shown in Fig. 1(b). The cross section before buckling is shown in Fig. 1(c), and the cross section after buckling deformation is shown in Fig. 1(d). The direction angle of deformation is 𝛽 d . 3. Beam element The finite element method is used to divide the compressive rod along the axis into several 3D straight beam elements. Bernoulli–Euler hypothesis is adopted. Fig. 2 shows the nodal displacement and nodal force of any 3D beam element. The node displacement vector is [ ]T 𝑞𝑒 = 𝑢𝑖 , 𝑣𝑖 , 𝑤𝑖 , 𝜃𝑖𝑥 , 𝜃𝑖𝑦 , 𝜃𝑖𝑧 , 𝑢𝑗 , 𝑣𝑗 , 𝑤𝑗 , 𝜃𝑗𝑥 , 𝜃𝑗𝑦 , 𝜃𝑗𝑧 . (2)
2. Physical model
where u is the axial displacement; v and w are the lateral displacement; 𝜃 is the angular displacement; i and j is the nodes of the beam element. The element displacement represented by nodal displacement is
The established physical model on buckling of compressed rods in cylinders is shown in Fig. 1, and the following assumptions are made: (i) The cross sections of the rods and cylinders are both circular. (ii) The initial annular gap between the rods and the cylinders exists. The axis of the rods before deformation coincides with the axis of the cylinders, and the axes are all straight lines. (iii) The rods are flexible/slender and are in a linear elastic state. The small deformation assumption is adopted. The annular gap between the rod and the cylinder is far lower than the axial length of the rod, so the lateral deformation of the rod is small. (iv) The cylinders are rigid. (v) After buckling, the rods will contact with the inner wall of the cylinders. The effect of friction is ignored. (vi) The ultimate buckling state of the rods subjected to a given axial compression load is studied. The transient effects of inertia force and damping force are neglected.
𝑢 = [𝑢, 𝑣, 𝑤, 𝜃]T = 𝑁 𝑞𝑒 ,
(3)
where N is the Hermitian shape function matrix given by Zienkiewicz et al. [38] ⎧ ⎫ ⎡ ⎪ 𝑁 𝑢 ⎪ ⎢𝑁 1 ⎪ 𝑁𝑣 ⎪ ⎢ 0 𝑁 =⎨ ⎬= ⎪ 𝑁𝑤 ⎪ ⎢⎢ 0 ⎪𝑁𝜃𝑥 ⎪ ⎣ 0 ⎩ ⎭
0 𝑁3 0 0
0 0 𝑁3 0
0 0 0 0 0 −𝑁4 𝑁1 0
0 𝑁4 0 0
where 𝑁1 = 1 − 𝑥𝑙 , 𝑁2 = 𝑥𝑙 , 𝑁3 = 1 − 3𝑥2
2𝑥3
𝑥2
𝑥3
𝑁2 0 0 0
3𝑥2 𝑙2
0 𝑁5 0 0
+
0 0 𝑁5 0
2𝑥3 , 𝑙3
0 0 0 0 0 −𝑁6 𝑁2 0
𝑁4 = 𝑥 −
0 𝑁6 0 0
⎤ ⎥ ⎥, ⎥ ⎥ ⎦
2𝑥2 𝑙
+
(4)
𝑥3 , 𝑙2
𝑁5 =
− 𝑙3 , 𝑁6 = − 𝑙 + 𝑙2 , l is the length of one beam element. The post-buckling deformation of the rods are subjected to both axial and bending actions. Based on Green–Lagrange strain and neglecting higher order terms, we get the generalized strain 𝜀 due to these actions 𝑙2
The two ends of the rods are hinged or clamped, and the mechanical model of the hinged boundary is shown in Fig. 1(a). To prevent 104
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
Fig. 2. Schematic diagram of beam element: (a) nodal displacements, and (b) nodal forces with respect to local coordinates (x, y, z). F and M is the force and moment, respectively.
[39]: ⎧ d𝑢 ⎧ ⎫ ⎪ d𝑥 + ⎪ 𝜀𝑥 ⎪ ⎪ ⎪𝜀 ⎪ ⎪ 𝜀 = ⎨ 𝑏𝑦 ⎬ = ⎨ ⎪𝜀𝑏𝑧 ⎪ ⎪ ⎪ 𝛾𝑥 ⎪ ⎪ ⎩ ⎭ ⎪ ⎩
1 2
(
d𝑣 d𝑥
)2 +
−𝑦 dd𝑥𝑣2 2
−𝑧 dd𝑥𝑤2 2
𝜌 dd𝑥𝜃
1 2
(
)2 ⎫ d𝑤 ⎪ d𝑥 ⎪ ⎪ ⎬ = 𝜀𝐿 + 𝜀𝑁𝐿 , ⎪ ⎪ ⎪ ⎭
of the cylinders, respectively. The shape is a ring with a certain thickness. The gap element has the following physical characteristics: (5)
(i) When the rods are not in contact with the inner wall of the cylinders, the compressive stiffness of the gap element approaches zero, which does not affect the free motion of the rods. (ii) When the rods are in contact with the inner wall of the cylinders, the gap element is completely compressed in a certain radial direction of the circumference, and the outer boundary of the gap element is tangent at the contact point. Its normal contact stiffness becomes quite large, which prevents the rods from intruding into the inner wall of the cylinders, but allows the rods to slide along the inner wall of the cylinders and produce contact forces, as shown in Fig. 3.
where 𝜀x , 𝜀by , and 𝜀bz are the axial strain along x axis, bending strain around y and z axis, respectively; 𝛾 x is the shear strain; 𝜌 is radial distance from the x axis of the rod. Generalized strain 𝜀 can be decomposed into linear strain 𝜀L and non-linear strain 𝜀NL . Matrix BL of kinematric relation corresponds to the linear part of −𝑁 ′′ 𝑣 −𝑁 ′′ 𝑤 𝑁 ′ 𝜃𝑥 ]T , the non-linear parts of Eq. (5): 𝐵𝐿 = [ 𝑁 ′ 𝑢 Eq. (5) are expressed by the vector of displacement gradients [40]: { } 𝑔̂ = 𝑣′ (6) 𝑤′ 𝜃𝑥 𝑣′′ 𝑤′′ 𝜃′ 𝑥 ,
The local coordinate system oyz of the gap element coincides with the local coordinate system of the beam element. The gap element is a twodimensional element whose displacements are the lateral displacements of the beam element. The displacement vector of the gap element at any position x of the ith node of the beam element is
and the shape function derivative matrix Γ is associated with the gradient of Eq. (6) { } 𝑁 ′𝑤 𝑁𝜃𝑥 𝑁𝑣 ′′ 𝑁𝑤 ′′ 𝑁 ′ 𝜃𝑥 . Γ = 𝑁 ′𝑣 (7) The linear stiffness matrix 𝑘𝑒0 and the geometric matrix 𝑘𝑒𝜎 have the following components [40]: 𝑘𝑒0
𝑙
=
∫0
𝑘𝑒𝜎 = 𝐹𝑐
𝐵𝐿 𝐷𝐵𝐿 d𝑥, T
𝑙
∫0
[ ]T 𝑢𝐺 = 𝑣𝑥 , 𝑤𝑥 = 𝑁𝐺 𝑞𝑒 ,
(11)
(8)
ΓT Γd𝑥,
(9)
where Fc the axial compression load, D is an elastic matrix and the expression is ⎡𝐸𝐴 ⎢ 𝐷=⎢ ⎢ ⎢ ⎣
𝐸 𝐼𝑧
𝐸𝐼𝑦 0
⎤ ⎥ ⎥, ⎥ 𝐺𝐼𝑝 ⎥⎦
(10)
in which E and G are, respectively, the elastic modulus and shear modulus of the beam; A is the cross-sectional area; Iz and Iy is second moments of area about the neutral axes z, y, respectively; Ip is the polar moment of inertia of area. 4. Gap element Fig. 3. Deformation of gap element: (a) initial position before buckling, and (b) contact between rods and the inner wall of the cylinders after buckling. FGN is the normal contact force, vx and wx are lateral displacement, 𝛽 G is the direction angle, D and d are the outer diameter of the rod and the inner diameter of the cylinder, respectively.
4.1. Physical properties and strains of gap elements A gap element [41,42] is a virtual element, its inner and outer boundaries are connected with the outer surface of the rods and the inner wall 105
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
𝑁𝑣 } is a 2 × 12 ma𝑁𝑤 trix, and the explicit expressions are given in Eq. (4). The lateral displacement UG and the direction angle 𝛽 G of the gap element at any position are described below: √ { 𝑈𝐺 = 𝑈𝑑 = 𝑣2𝑥 + 𝑤2𝑥 (12) ( ), 𝛽𝐺 = 𝛽𝑑 = 180𝑛 + arctan 𝑣𝑥 ∕𝑤𝑥
(i) The free state is in accordance with the following relations
where shape function matrix of gap element 𝑁𝐺 = {
⎧𝜀 − 1.0 < −𝜀𝑁,min ⎪| 𝐺𝑁 ⎨|𝐹𝐺𝑁 − 𝐹𝑁0 || ≤ 0.0 , ⎪𝐺𝑘 = 𝐺𝑘𝑗 ⎩
(ii) The elastic contact state is in accordance with the following relations ⎧−𝜀 ≤ 𝜀𝐺𝑁 − 1.0 ≤ 𝜀𝑁,max ⎪| 𝑁,min , ⎨|𝐹𝐺𝑁 − 𝐹𝑁0 || ≥ 0.0 ⎪𝐺𝑘 = 𝐺𝑘𝑐 ⎩
where Ud and 𝛽 d are the lateral displacement and the direction angle of the beam element, n is the half circle number, the initial gap rG of the gap element is expressed as 𝑟𝐺 = (𝐷 − 𝑑 )∕2,
(17)
where 𝜀N,min , 𝜀N,max and FN0 are small quantities larger than zero to control the convergence accuracy of gap element iteration. 𝜀N,min and 𝜀N,max are the minimum and maximum iteration errors for gap element strain, respectively. FN0 is the iterative error of contact force of gap element. In order to make the contact between the rods and the inner wall of the cylinders satisfy the definite solution of the gap element, the stiffness Gk of the gap element should be revised continuously in the iterative calculation. The following principles should be adopted in the correction:
(13)
where D and d are the outer diameter of the rod and the inner diameter of the cylinder, respectively. Thus, the normal strain 𝜀GN of the gap element is 𝜀𝐺𝑁 = 𝑈𝐺 ∕𝑟𝐺 .
(16)
(14)
4.2. Criterion and definite formula for the contact state of gap element One of the two states after the deformation of the rods is
(i) The initial contact stiffness Gk0 and the initial penetration UPe,0 of the gap element are given for correction in the iteration: { 𝐺𝑘0 = 𝑓Kn ⋅ 12𝑙𝐸𝐼 3 , (18) 𝑈Pe,0 = 𝑓Pe ⋅ 𝑑∕2
(i) Free state: the rods are separated from the inner surface of the cylinders and do not contact with the wall of the cylinder. (ii) Elastic contact: the rods are in contact with the inner wall of the cylinders, a small amount of penetration between them exists. The contact force is linear with penetration.
where fKn is the coefficient of contact stiffness, and fPe is the penetration coefficient.
As shown in Fig. 4, the criterion for describing the above state by gap element is (1) free state: 𝜀GN < 1.0, GK = Gkj ,SG1 ∈ S; (2) elastic contact: 𝜀GN ≥ 1.0, GK = Gkc ,SG2 ∈ S, here SG1 and SG2 represent the region where the gap element is located, S represents the contact region. The relation between the stiffness Gk of the gap element and the strain 𝜀GN is shown in Fig. 4, in which Gk0 is the initial value of the contact stiffness of the gap element, Gkj is the compressive stiffness of the gap element, and Gkc is the normal contact stiffness of the gap element after the contact between the rods and the inner wall of the cylinders. Obviously, the gap element is actually a nonlinear element with variable stiffness. The contact force between the rode and the inner surface of the cylinder is { , 𝑈𝐺 − 𝑟𝐺 < 0 0 ( ) 𝐹𝐺𝑁 = . (15) 𝐺𝑘 𝑈𝐺 − 𝑟𝐺 , 𝑈𝐺 − 𝑟𝐺 ≥ 0
(i) In the free state, the compressive stiffness Gkj of the gap element should be minimized, so that the free movement of the rod will not be affected. (ii) For elastic contact, the stiffness of gap element Gk is related to the bending stiffness of the rods. (iii) The contact between the rods and the cylinder wall is multidirectional. The stiffness of the gap elements between the adjacent beam elements is mutually affected among the gap elements. Therefore, this effect should be taken into account when correcting the stiffness of the gap element. Each gap element should be corrected repeatedly until all the gap elements satisfy the definite solution. The introduction of gap element enables the connection of beam elements and cylinders into a continuum.
where FGN is the contact force, Gk is the stiffness of the gap element. When contact occurs, the small penetration is UPe = UG − rG . The contact force FGN can be further transformed to the additional nodal force of the beam element. Based on the strain 𝜀GN in Eq. (14) and contact force FGN in Eq. (15), the definite formula of iteration calculation for gap element can be written as:
4.3. Static buckling equilibrium equation Combing of beam element and gap element, the equilibrium equation is (𝑘𝑒0 + 𝑘𝑒𝜎 + 𝑘𝑒𝐺 )𝑞𝑒 = 𝐹e + 𝐹𝐺𝑒 ,
(19)
Fig. 4. Relation between strain and stiffness of gap element. (a) The normal contact stiffness less than the initial contact stiffness: Gkc < Gk0 , and (b) the normal contact stiffness greater than the initial contact stiffness: Gkc > Gk0 . Gkj is the compressive stiffness of the gap element.
106
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
where Fe and 𝐹𝐺𝑒 is the equivalent nodal force of beam and gap element, respectively; 𝑘𝑒𝐺 is the gap stiffness matrix given by 𝑘𝑒𝐺 = 𝑁𝐺T 𝐺𝑘 𝑁𝐺 .
The dynamic relaxation method is a numerical method which transforms the static problem into a quasi-dynamic problem according to the D’Alembert principle and finds the equilibrium form of the structure through parameter iteration. The structure is discretized into elements with certain virtual mass. Under the action of unbalanced force, these discrete elements will inevitably move along the direction of unbalanced force and reduce the overall unbalanced force of the structure. It also discretizes the vibration process of a spatially discretized structural system in time and traces the vibration process of the system step by step. By introducing the dynamic damping to consume the oscillation energy, the kinetic energy of the structure motion is gradually reduced and tends to zero, that is, to reach the static equilibrium state. Dynamic relaxation method is also widely used in static buckling analysis of structures [44,45]. However, the post-buckling analysis of cylindrically constrained rods has rarely been reported by using this dynamic relaxation method.
(20)
After coordinates transformation and assembling of all beam elements and gap elements, we obtain the global equilibrium equation for the geometric and contact nonlinear static analysis of the rods in cylinders: (𝐾0 + 𝐾𝜎 + 𝐾𝐺 )𝑑 = 𝐹b + 𝐹𝐺 ,
(21)
where K0 , K𝜎 and KG are, respectively, the global linear stiffness matrix, geometric stiffness matrix and gap stiffness matrix of the rods; d is the node displacement vector; Fb and FG are, respectively, the equivalent nodal force vector of beam element and gap element with respect to global coordinates. The nonlinear Eq. (21) shows that with the increase of axial compressive load Fc , the geometric stiffness matrix KG decreases, resulting in the rod buckling. The lateral displacement of the rods suddenly increases, leading to the contact between the rods and the inner cylinder wall. In solving the nonlinear Eq. (21), it is necessary to distinguish and modify the definite solution of the gap element in the same iteration step. Only when all these conditions are satisfied, the generalized displacement and internal force can converge to the exact solution.
5.3. Dynamic relaxation method for static buckling of compressed rods in cylinders In this paper, with the effect of cylindrical constraint, the dynamic relaxation method is applied to solve the static buckling problem of compressed rods in cylinders. The static buckling model of Fig. 1(a) is transformed into the dynamic one shown in Fig. 5. In order to reduce the impact of axial compression load on the rods, the sinusoidal wave loading is adopted, and the time history of axial compression load is as follows { 𝐹c sin (π𝑡∕𝑇 ), 0 ≤ 𝑡 ≤ 𝑡1 𝐹 c (𝑡 ) = , (22) 𝐹c , 𝑡 > 𝑡1
5. Dynamic relaxation method 5.1. Problems in static buckling calculation Eq. (21) is a complex equation involving the dual nonlinearity of geometry and contact. For static buckling calculation, some problems such as convergence difficulty and algorithm instability still exist. They are manifested in the following aspects:
where the axial compression load 𝐹c = 𝑘 𝜋 𝐿𝐸𝐼 2 , k is a dimensionless compressive load; T is the characteristic time for loading; t is the response time; t1 = T/4. The initial geometric imperfections will affect the critical buckling load. The small perturbation load is initially applied and revoked in subsequent analysis. Thus, the perturbation load will not affect the critical buckling load. 2
(i) In the process of buckling, the jumping change between different buckling modes always exists. The strong nonlinear mechanical behaviors such as contact and separation between the rods and the cylindrical wall cause a sudden change in the contact state in the calculation process. These physical phenomena are easy to cause abnormal termination of calculation and have difficulty in convergence. Especially in the critical buckling state, the deflection suddenly increases and the contact state cannot be recognized. The rods penetrate through the cylindrical wall, causing the calculation to fail. (ii) The algorithm instability shows that the convergence solution is not unique, and the configuration of the rod buckling is arbitrary. For example, when the plane sinusoidal buckling mode is formed, this mode rotates in the cylinders with the increase of axial compression load, but the helical buckling mode cannot be formed. Because of the sensitivity to iterative load step, buckling modes can sometimes be formed, but sometimes they cannot. After the helical buckling mode is formed, the buckling modes of the right and the left helix may exist simultaneously at different axial positions. To solve the problems of convergence difficulty and algorithm instability in the buckling of rods in cylinders, a dynamic relaxation method is proposed. 5.2. Basic idea of dynamic relaxation method The basic idea of dynamic relaxation [43] originates from the steadystate solution of a single degree of freedom damping system. Large damping is used to attenuate the dynamic response of the system, and the static part of the transient dynamic solution (i.e., the steady state) is found.
Fig. 5. Mechanical model based on dynamic relaxation method: (a) the position where these loads [axial compression force Fc (t) and perturbation loads Fp,i (t)] are applied, (b) the plane projection of the perturbation loads, and (c) the curve of the applied loads vs. time. 107
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
In order to easily form the spatial deformation, lateral perturbation loads Fp,1 (t), Fp,2 (t) and Fp,3 (t) are applied at the positions of 0.25L intervals along the rod length: { 𝐹p,𝑖 ⋅ sin (2π𝑡∕𝑇 ), 0 ≤ 𝑡 ≤ 𝑡2 𝐹p,𝑖 (𝑡) = , (23) 0, 𝑡 > 𝑡2 where i = 1, 2, 3, t2 = T/2. The directions of perturbation loads are shown in Fig. 5(a) and (b). Considering the inertia force and damping force in the dynamic model, the perturbation load and axial compression load are applied in Eq. (21). The basic dynamic buckling equation of the compression rods in cylinders is obtained: 𝑀𝑑 ′′ + 𝐶𝑑 ′ + 𝐾𝑑 = 𝐹 (𝑡),
(24) Fig. 6. Effect of compressive load on the lateral displacement Ud of the compression rod without cylindrical constraints (𝛼 = 𝛼 0 , 𝛼 0 = 2𝜔).
where M, C and K are, respectively, the mass, damping and stiffness matrix of the compression rods, K = K0 + K𝜎 + KG ; d′′ and d′ are node acceleration and velocity vectors of the rods, respectively; F(t) is the load vector of the rods, F(t) = Fb (t) + FG (t). In this paper, Rayleigh damping is used to set the Rayleigh damping as a large value to attenuate the dynamic response of the system. Considering the effect of 𝛼 damping in the Rayleigh damping, then the damping matrix C = 𝛼M. The damping force 𝐶 𝑑 ′ = α𝑀𝑑 ′ . Thus, the faster the movement speed of the rods, the greater the damping force. The bigger damping force can effectively control the jumping change of buckling mode and the sudden change of contact state, which is helpful to improve the calculation efficiency and stability. Newmark numerical integration method is used to calculate Eq. (24). When t > t2 , the perturbation load is equal to 0, the axial compression load remains unchanged. Let the left item of the static buckling Eq. (21) be ( ) (25) 𝐾0 + 𝐾𝜎 + 𝐾𝐺 𝑑 = 𝐹k ,
Fig. 7. Effect of 𝛼 damping on the lateral displacement Ud of the compression rod without cylindrical constraints (k = 0.998, 𝛼 0 = 2𝜔).
and the right item of Eq. (21) be 𝐹b + 𝐹𝐺 = 𝐹J .
6.1. The effect of 𝛼 damping without cylindrical constraints
(26)
The kinetic Eq. (24) is rewritten into 𝑀𝑑 + 𝐶𝑑 = 𝐹J (𝑡) − 𝐹k (𝑡) = 𝐹R (𝑡), ′′
′
𝜋2 𝐿2
(27)
The natural frequency of the hinged rod at both ends is 𝜔 = √ 𝐸𝐼 = 80.40 rad. The reference 𝛼 damping is 𝛼 0 = 2𝜔. Figs. 6 and 𝜌𝐴
7 show the effect of the compression load 𝑘 = 𝐹c ∕( 𝜋 𝐿𝐸𝐼 2 ) and 𝛼 damping on the lateral displacement of the middle of the compression rod, respectively. As can be seen from Fig. 6, when k<0.998, the lateral displacement Ud decreases with the withdrawal of the perturbation load, and tends to zero gradually, so the rod is in a stable state. When k = 0.998, Ud remains unchanged as perturbation load is removed. The compressed rod is in a critical state of equilibrium, and the corresponding load 𝐹c = 2
where FR (t) is the residual force vector. With the extension of calculation time, the axial compression load Fc remains unchanged (when t > t1 ), and the perturbation load is 0 (when t > t2 ). Due to the existence of damping, the transient response of the rods tends to be static, and the rod deformations d will also tend to be constant after buckling. Residual force vector FR (t) containing inertia force Md′′ and damping force Cd′ decreases to 0, then Fk = FJ , the contact force FG is no longer changing. As a result, the solution of nonlinear static buckling of rods in cylinders is achieved.
0.998 𝜋 𝐿𝐸𝐼 is the critical value. When k>0.998, Ud increases with the 2 withdrawal of perturbation load, so the rod is in an unstable state. The dimensionless critical load (k = 0.998) calculated by us is 0.2% less 2
6. Analysis of effect factors on buckling mode of compressed rods
than the classical Euler formula 𝐹c = 1.000 𝜋 𝐿𝐸𝐼 2 , which indicates that the dynamic relaxation method has a high calculation accuracy. From Fig. 7, with different 𝛼 damping, the lateral displacement is still in the equilibrium state with the withdrawal of the lateral force. The larger the 𝛼 damping is, the smaller the ultimate displacement will be. Therefore, the 𝛼 damping has no effect on the critical load of the compression rod. It only affects the displacement of the ultimate equilibrium. 2
Firstly, the effect of 𝛼 damping on buckling of a compression rod is analyzed without the constraint of the cylinder. Then, with the cylindrical constraints, the effects of 𝛼 damping, element length, contact stiffness and penetration on post-buckling of the compression rod are analyzed under a given load. The calculation parameters of the rod are as follows: the length L = 0.864 m, diameter d = 4.76E−3 m, density is 7850 kg/m3 , elastic modulus E = 204.78 GPa. The inner diameter of the cylinder is D = 38.1E−3 m. The initial gap between the cylinder and the rod is rG = 16.67E−3 m. The beam element length is l = 8E−3 m. The response time t = 20 s, and the characteristic time for loading T = 5 s. Lateral perturbation load Fp,1 = Fp,2 = Fp,3 = 1 N. Compressive stiffness of gap element Gkj = 0. The minimum iterative error of gap element strain is 𝜀N,min = 1E−5. The maximum value of permissible penetration UPe,max ≤ 1E−3 m. The maximum iteration error of gap element strain resulted is 𝜀N,max = 0.06.
6.2. The effect of 𝛼 damping with cylindrical constraints Given the dimensionless compression load (k = 30), the dynamic relaxation method is used to calculate the post-buckling of the compression rod under different 𝛼 damping conditions. Fig. 8 shows the effect of the 𝛼 damping on lateral displacement, velocity and acceleration at the 0.25L position. For 𝛼 = 300𝛼 0 , Fig. 9 shows the buckling process of the post-buckling mode. 108
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
Fig. 8. The effect of 𝛼 damping on (a) lateral displacement Ud , (b) velocity Ud ′ and (c) acceleration U′′d at the 0.25L position (𝛼 0 = 2𝜔). Fig. 9. Stabilization process of post-buckling mode of the rod with cylindrical constraints (k = 30, 𝛼 = 300𝛼 0 ): (a) 3D deformation curve, and (b) top view of the deformation at certain times.
In Fig. 8, during the loading stage of axial compression and lateral perturbation load, the lateral displacement Ud increases with time. When Ud reaches the annular gap (rG = 0.01667 m), the rod contacts with the inner wall of the cylinder, and the lateral displacement remains unchanged. The larger the 𝛼 damping, the longer the time needed to reach the contact between the rod and the cylindrical wall. If the 𝛼 damping is too small, it will lead to excessive lateral velocity and acceleration of the rod, and produce a greater impact vibration. In the stable stage of axial compression load, the lateral displacement Ud tends to be constant, and the lateral velocity and acceleration tend to be zero. The lateral velocity and acceleration can be restrained by choosing a larger 𝛼 damping. Under the condition of long response time, the effect of lateral velocity and acceleration is eliminated, and the residual force FR in Eq. (27) is reduced to 0. As a result, a stable buckling deformation can be obtained. In Fig. 9, during the loading stage of axial compression load and lateral perturbation load (t<5 s), the compression rod deforms into a three-dimensional deflection curve, and the contact state of two points (t = 3.25 s), three points (t = 3.5 s) and point-line-point (t = 4.25 s) with the cylinder wall gradually occurs. In the stable stage of compressive load, the buckling mode does not change with time when t ≥ 10 s. Therefore, the deformation of the rod does not undergo sinusoidal buckling mode, the stable and helical post-buckling modes can be directly calculated by using the dynamic relaxation method with its loading style of axial compression and perturbation load.
flection curve. When the element length is less than or equal to 8E−3 m, i.e., when the element number is greater than or equal to 108, both the lateral displacement curves and the angular displacement curves coincide with each other. 6.4. The effect of contact stiffness and penetration With the given 𝛼 damping (𝛼 = 300𝛼 0 ) and compressive load (k = 30), the effect of the initial contact stiffness factor fKn and penetration factor fPe of the gap element on the calculation results is given in Fig. 11. From Fig. 11, after nonlinear iterative calculation, the normal contact stiffness Gkc and penetration UPe are both related to fKn and fPe . Gkc is directly proportional to the initial contact stiffness factor fKn , and is inversely proportional to the penetration factor fPe . UPe is inversely proportional to fKn , and is proportional to fPe . For fKn = 0.01 and fPe = 1.0, the value of UPe is close 1E−3 m, reaching the maximum value. When fKn is very small, the contact force FGN is inversely proportional to fPe . When fKn ≥ 0.3, FGN is no longer changing with fKn , fPe , and UPe after the iteration is less than 0.1E−3 m. 7. The analysis of characteristic buckling mode 7.1. Characteristic buckling modes and critical loads The dynamic relaxation method is used to carry out the post-buckling finite element calculation of the compressive rods. Various characteristic buckling modes and their contact positions obtained by varying the compressive load are shown in Fig. 12, and the corresponding deflection curves are shown in Fig. 13. The point contact in the pictures is represented by the first contact points c1 and pentagram. The end of the
6.3. The effect of element length Under the given 𝛼 damping (𝛼 = 300𝛼 0 ) and compressive load (k = 30), Fig. 10 shows the effect of beam element length on the de109
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
Fig. 10. The effect of element length on (a) lateral displacement Ud , and (b) angular displacement 𝛽 d . The element length from l1 to l8 is 2E−3 m, 4E−3 m, 6E-3 m, 8E−3 m, 12E−3 m, 16E−3 m, 24E−3 m and 32E−3 m, respectively, i.e., the rods are divided into 432, 216, 144, 108, 72, 54, 36 and 27 beam elements, respectively.
Fig. 11. The effect of contact stiffness and penetration (𝛼 = 300𝛼 0 , k = 30): (a) initial contact stiffness factor fKn = Gk0 / Gkc , (c) penetration UPe , and (d) contact force FGN .
continuous contact section is represented by the second contact points c2 and solid dots. According to the contact state between the compression rods and the cylinders, the characteristic buckling modes with single point contact (k = 1.00), two points contact (k = 4.00), three points contact (k = 7.83) and point-line-point contact are given in Figs. 12 and 13, respectively. In the buckling configuration of point-line-point contact, they are divided into just point-line-point deformation (k = 8.83), helix pitch measured by contact point c1 (k = 20.43) and continuous contact section (k = 30.16), respectively.
•
(12EI/l3 ), (b) normal contact stiffness
The characteristics of various buckling modes are as follows: (i) For buckling modes with single point contact (1.00≤ k ≤ 3.99), the contact position is located in the middle of the compressive rod. The value 𝜉 c1 = 0.5. (ii) For buckling modes with two points contact (4.00≤ k ≤ 7.82), the middle part of the rod is detached from the wall of the cylinder, and the contact points on both sides of the middle part appear. With the increase of compressive load, the two contact points gradually move along to the rod end. The position of the contact 110
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
Fig. 12. Characteristic buckling modes and their contact positions (𝜉 = x/L) for hinged boundary: (a) buckling modes, and (b) location of the contact point. Mode A: single point contact; Mode B: two points contact; Mode C: three points contact; Mode D: just point-line-point contact; Mode E: helix pitch measured by contact point c1; and Mode F: helix pitch measured by continuous contact section. The last three buckling modes are point-line-point contact deformation.
Fig. 13. Deflection curves of buckling modes for hinged boundary: (a) lateral displacement, and (b) direction angle vs. axial distance, respectively.
111
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
Table 1 Dimensionless shear force, bending moment and contact force in six post-buckling modes. Boundary
Post-buckling mode
k
qmax
q|𝜉 = 0.5
mmax
m|𝜉 = 0.5
n1
nh,max
nh |𝜉 = 0.5
Hinged
Mode A Mode B Mode C Mode D Mode E Mode F
1.00 4.00 7.83 8.83 20.43 30.16
0.936 1.033 1.732 1.741 1.756 1.715
0.006 0.450 1.056 1.072 1.043 1.017
0.981 0.722 1.214 1.218 1.230 1.213
0.981 0.048 0.935 0.945 1.037 1.030
0.007 0.903, 0.903 1.671, 0.003, 1.671 1.708 1.758 1.773
– – – 0.149 1.611 1.375
– – – 0.014 1.021 1.063
Clamped
Mode A Mode B Mode C Mode D Mode E Mode F
4.00 8.19 12.00 12.49 28.16 36.67
1.412 1.149 1.425 1.427 1.440 1.443
0.013 1.025 1.045 1.049 1.015 1.019
0.998 0.812 1.007 1.009 1.018 1.023
0.998 0.411 0.902 0.909 1.016 1.023
0.004 0.410, 0.410 1.020, 0.033, 1.020 1.039 1.102 1.108
– – – 0.223 1.029 1.055
– – – 0.223 1.029 1.055
point c1 decreases to 𝜉 c1 = 0.278, and the angle between the two contact points c1–c1 increases to 𝛽 d,c1 − c1 = 135.8°. (iii) For buckling modes with three points contact (7.83≤ k ≤ 8.82), the middle part of the rod contacts with the cylinder wall again. The contact point c1 decreases to 𝜉 c1 = 0.267, and the angle increases to 𝛽 d,c1 − c1 = 152.2°. (iv) For buckling modes with point-line-point contact, the contact point in the middle of the rod is changed into continuous contact section. For helix pitch measured by contact point c1, contact points c1 and c2 are located at 𝜉 c1 = 0.171 and 𝜉 c2 = 0.302, respectively. The values 𝜉 c1,c1 = 0.657, 𝛽 d,c1 − c1 = 360°, 𝛽 d,c2 − c2 = 228°. For helix pitch measured by continuous contact section, the values 𝜉 c1 = 0.142, 𝜉 c2 = 0.246, 𝜉 c2,c2 = 0.508, 𝛽 d,c1 − c1 = 360° +133.6°, 𝛽 d,c2 − c2 = 360°.
middle of the rods, except for the first two modes of hinged boundary and the first mode of clamped restraint, the shear force q|𝜉 = 0.5 for the other modes is about 1.0. The maximum bending moment mmax in mode A with hinged boundary is in the middle of the rod, and the mmax in the other modes are in the suspended section, as shown in Fig. 14(c), (d) and Table 1. mmax is 1.0 and 0.7, respectively, in the first two post-buckling modes, and the values in the others is about 1.2. For clamped boundary, the maximum bending moment exists at the rod ends, and the maximum values also exist in the suspended section and in the middle of the rod. mmax is 0.8 for the mode B, and is about 1.0 for other post-buckling modes. Regardless of the boundaries, m|𝜉 = 0.5 is about 1.0 for the last four modes. The contact force n1 at the first contact point c1 increases gradually in the first three modes as shown in Table 1. The value n1 increases from a small value in mode A to 0.9 and 0.4 in Mode B, and finally to 1.7 and 1.0 in Mode C for hinged and clamped boundary, respectively. The latter three modes are the point-line-point contact deformations, the contact force n1 increase slightly, and become 1.8 and 1.1. The contact force nh in the continuous contact section is not uniform along the rod, as shown in Fig. 14(e), (f) and Table 1. The maximum values exist at the ends of the continuous contact section for hinged boundary; however, the maximum values are in the middle of the continuous contact section for clamped boundary. The contact force distribution is more uniform in the middle of the continuous contact section. The contact force nh |𝜉 = 0.5 in the middle is about 1.1 for the helical buckling mode (mode F). It is emphasized that the mechanical parameters in Fig. 14 and Table 1 are dimensionless values, and the actual values (in SI unit) can be obtained by introducing the values in Table 1 into Formula (28). For example, the contact force F1 in mode E are respectively 237.9 N and 241.3 N for hinged and clamped boundary; the contact force Fh in the middle of the continuous contact section are 1605.4 N/m and 3074.0 N/m, respectively. In mode F, F1 are 430.3 N and 360.5 N, Fh are 3642.6 N/m and 5344.3 N/m, respectively.
It should be emphasized that in Fig. 12 we call the section from the first contact point to the end the first suspended section, Lc1 = 𝜉 c1 L; the section from the first contact point to the second contact point is referred to as the second suspended section, Lc2 = (𝜉 c2 − 𝜉 c1 )L; the length of pitch measured by contact point c1 is Lh1 = 𝜉 c1,c1 L; the length of pitch measured by contact point c2 is Lh2 = 𝜉 c2,c2 L. 7.2. Mechanical parameters of various post-buckling modes Introduce dimensionless shear force q, bending moment m and contact forces n1 , nh as √ 𝐹c ⎧𝐹 = 𝑞 𝐹 c 𝑟 𝐺 𝑠 2 2𝐸𝐼 ⎪ 𝐹c 𝑟𝐺 ⎪ ⎪𝑀 = 𝑚 2 √ , (28) ⎨ 𝐹c ⎪𝐹 1 = 𝑛 1 𝐹 c 𝑟 𝐺 2 2 𝐸𝐼 ⎪ ⎪𝐹 = 𝑛 𝐹𝑐2 𝑟𝐺 ℎ 4𝐸𝐼 ⎩ h where Fs is the shear force (in N); q is the dimensionless shear force; M is the bending moment (in N•m); m is the dimensionless bending moment; F1 is the contact force of the contact point c1 (in N); n1 is the dimensionless contact force of the contact point c1; Fh is the contact force of the continuous contact section, N/m; nh is the dimensionless contact force of the continuous contact section. In post-buckling modes for hinged and clamped boundaries, the distribution curves of dimensionless shear force, bending moment and contact force in continuous contact section are shown in Fig. 14, and the dimensionless shear force, bending moment and contact force of point and line contact are shown in Table 1. From Fig. 14(a) and (b) and Table 1, the maximum shear force qmax of the first two post-buckling modes is about 1.0 and that of the last four post-buckling modes is about 1.7 for hinged boundary. For clamped boundary, the maximum shear force is in the suspension section, qmax is about 1.1 for the mode B, and qmax is about 1.4 for the others. In the
7.3. Verification of critical load for various post-buckling modes In Table 2, the dimensionless critical loads of the characteristic buckling modes of the hinged and clamped boundary conditions are given and compared with the results in the literatures. For the hinged boundary, the critical loads of the first two modes in this paper are similar to those given in open literature [8], except the critical value of helical buckling (mode F). For the clamped boundary, the critical load obtained in this paper is slightly different from that obtained by the elastic model [14], and is in good agreement with that obtained by the small deformation hypothesis [13]. It should be noted that more critical loads of buckling modes are calculated in this paper than in the literatures. 112
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
Fig. 14. Non-dimensional load curves for post-buckling modes: (a) and (b) dimensionless shear force q vs. x/L, (c) and (d) dimensionless moment m vs. x/L, (e) and (f) dimensionless contact force nh vs. x/L for continuous contact sections.
Table 2 Dimensionless critical loads of the characteristic buckling modes. Post-buckling mode
Hinged boundary
Clamped boundary
Present
Huang et al. [8]
% Difference
Present
Chen and Fang [13]
% Difference
Fang et al. [14]
% Difference
Mode A Mode B Mode C Mode D Mode E Mode F
1.00 4.00 7.83 8.83 20.43 30.16
1.00 4.00 – 9.09 – 34.15
0 0 – −2.86 – −11.68
4.00 8.19 12.00 12.49 28.16 36.67
4.00 8.18 11.92 12.52 – –
0 0.12 0.67 −0.24 – –
4.00 8.08 10.80 12.03 – –
0 1.36 11.11 3.82 – –
113
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
Fig. 15. Deflection curves of helical buckling modes (mode F in Fig. 12) for (a) hinged boundary (k = 30.16, 𝜂 1 = 1.735, 𝜂 2 = 1.269, 𝜂 h = 6.192, 𝜂 L = 12.200), and (b) clamped boundary (k = 36.67, 𝜂 1 = 2.709, 𝜂 2 = 0.888, 𝜂 h = 6.258, 𝜂 L = 13.452) conditions. 𝜂 1 , 𝜂 2 , 𝜂 h and 𝜂 L is the dimensionless length of the first non-contact √ 2 𝐹c section, the second non-contact section, the continuous contact section and the compressed rod. 𝑘 = 𝐹c ∕( 𝜋 𝐿𝐸𝐼 ), 𝜂 = ⋅ 𝑥. 2 2𝐸𝐼 Table 3 Length of transition section in helical buckling modes (mode F in Fig. 12). Boundary
Author
𝜂1
𝜂2
𝜂1 + 𝜂2
𝜂3
𝜂1 + 𝜂2 + 𝜂3
Hinged
Sorenson and Cheathamjr [12], Mitchell [36], Huang et al. [15] This paper
1.747 1.735
1.347 1.269
3.094 3.004
1.677 0
4.771 3.004
Clamped
Mitchell [11] Liu et al. [22] Sorenson and Cheathamjr [12], Mitchell [36], Huang et al. [15] This paper
3.178 2.718 2.700 2.709
0 0 0.903 0.888
3.178 2.718 3.603 3.597
0 1.211 1.846 0
3.178 3.929 5.449 3.597
Table 4 Relevant mechanical parameters of helical buckling modes (mode F in Fig. 12). Boundary
Author
𝑚t _max
mh
q1
n1
nh
Hinged
Sorenson and Cheathamjr [12], Huang et al. [15] This paper
1.198 1.213
1.000 1.030
2.870 1.715
1.734 1.773
1.000 1.063
Clamped
Liu et al. [22] Mitchell [36] Sorenson and Cheathamjr [12], Huang et al. [15] This paper
1.000 0.766 1.010 1.020
1.000 – 1.000 1.023
0.543 0.457 0.342 1.443
1.253 – 1.069 1.108
– – 1.000 1.055
eratures [12,15,36]. The dimensionless total length 𝜂 L in this paper is 𝜂 L = 2(𝜂 1 + 𝜂 2 ) + 𝜂 h . According to Formula (22) and Formula (29), the relation between dimensionless total length 𝜂 L and dimensionless load k is
7.4. Verification of the lengths and relevant mechanical parameters of helical buckling modes Introduce dimensionless length 𝜂 √ 𝐹c ⋅ 𝑥, 𝜂= 2𝐸𝐼
𝑘 = 2𝜂𝐿2 ∕𝜋 2 ,
(29)
(30)
therefore, the critical load of helical buckling is less than that obtained by Sorenson and Cheathamjr [12], Huang et al. [15], and Mitchell [36]. Table 4 gives the relevant mechanical parameters of helical buckling modes and the results from literatures. In Table 4, 𝑚t _max is the maximum bending moment in transition section, mh is the bending moment in helical section, q1 is the maximum shear force in the transition section, n1 is the contact force at the contact point c1, nh is the contact force at the helical section, respectively. The values 𝑚t _max , mh , n1 and nh are in good agreement with the results of Sorenson et al. [12] and Huang et al. [15]. However, the maximum shear force in the transition section q1 is quite different from that in the literature. The value q1 of the clamped boundary is higher than those in the literatures [12,15,22,36], while this value of the hinged boundary is lower than those in literatures [12,15].
Under hinged and clamped boundary conditions, the deflection curves of helical buckling (mode F) are shown in Fig. 15. The dimensionless length of transition section is given and compared with the literatures results in Table 3. Both the first and second non-contact sections 𝜂 1 and 𝜂 2 exist in the helical buckling configuration, which are basically consistent with the results in the literatures [12,15,36]. In Fig. 15, the slope of angle 𝛽 d increases gradually in the first noncontact section, and tends to be constant in the second non-contact section, and finally remains unchanged in the continuous contact section. On the other hand, the length 𝜂 h of the continuous contact section calculated in this paper is about 2𝜋, which is basically consistent with the analytical results of the continuous contact hypothesis in the literature [29]. Therefore, the calculated transition section includes only the two non-contact sections 𝜂 1 and 𝜂 2 , does not have perturbed-helix section 𝜂 3 . Since there is no perturbed-helix section in this paper, the total length of the transition section (𝜂 1 + 𝜂 2 ) in this paper is significantly less than that of the three-section model (𝜂 1 + 𝜂 2 + 𝜂 3 ) used in the lit-
8. Conclusion Considering the geometric nonlinearity of the compressed rods and the contact nonlinearity with the external cylinder, the rods are discretized into beam elements, and the gap element is constructed by the 114
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
finite element method. To solve the problems of convergence difficulty and algorithm instability in static buckling analysis, a dynamic relaxation method for static buckling analysis of rods in cylinders is introduced. Based on the classical Euler formula, the effectiveness of dynamic relaxation method in calculating the static buckling of the compressed rod is investigated by choosing a larger damping coefficient. Then, the effect factors of dynamic relaxation method in post-buckling analysis of the compressed rods in cylinders are analyzed, with the compressive loads being far greater than that in Euler’s formula. To reduce the dynamic impact effect, a larger 𝛼 damping should be chosen. Under the condition of applying the initial lateral perturbation load, the post-buckling mode with helical shape and stable configuration can be obtained directly by avoiding the sinusoidal buckling process. In order to obtain the invariant deflection curve in the helical configuration, the number of beam elements should be at least 100. The effects of initial contact stiffness and initial penetration of the gap element on the final contact stiffness, penetration and contact force are also investigated. When the initial contact stiffness coefficient is greater than 0.3, the constant contact force can be obtained, and the final penetration is less than 0.1E−3 m. We studied the characteristic post-buckling modes with single point, two points, three points and point-line-point contact, respectively. In the configuration of point-line-point contact, it is further divided into just point-line-point contact and helical buckling deformation measured by the first contact point and continuous contact section. Except the critical load of helical buckling, the critical values of other configurations are in good agreement with the results of literatures, which indicates that the critical values obtained by the dynamic relaxation method have high accuracy. We also analyzed the distribution of shear force, bending moment and contact force in various post-buckling modes. In the helical buckling mode, we have found exceptions which are significantly different from those in the literatures. The transition section consists of the first non-contact section and the second non-contact section, but no perturbed-helix section. The more interesting conclusion is that, although there is no perturbed-helix, the length of two noncontact sections and the helical section obtained are in good agreement with the results in the literatures. It is found that the critical load of helical buckling obtained in this paper is smaller than that of the literature. The bending moment and contact force obtained are consistent with the results in the literatures, except the shear force.
[7] Lubinski A. A study of the buckling of rotary drilling strings. Drilling and production practice. American Petroleum Institute; 1950. [8] Huang W, Gao D, Liu Y. A study of tubular string buckling in vertical wells. Int J Mech Sci 2016;118:231–53. [9] Zhang Q, Jiang B, Huang W, Cui W, Liu J. Effect of wellhead tension on buckling load of tubular strings in vertical wells. J Petrol Sci Eng 2018;164:351–61. [10] Yue Q, Liu J, Zhang L, Zhang Q. The posting-buckling analysis and evaluations of limit drilling length for coiled tubing in the sidetrack horizontal well. J Petrol Sci Eng 2018;164:559–70. [11] Mitchell RF. Buckling behavior of well tubing: the packer effect. SPE J 1982;22(5):616–24. [12] Sorenson KG, Cheathamjr JB. Post-buckling behavior of a circular rod constrained within a circular cylinder. J Appl Mech 1986;53(4):929–34. [13] Chen JS, Fang J. Deformation sequence of a constrained spatial buckled beam under edge thrust. Int J Non Linear Mech 2013;55(10):98–101. [14] Fang J, Li SY, Chen JS. On a compressed spatial elastica constrained inside a tube. Acta Mech 2013;224(11):2635–47. [15] Huang W, Gao D, Wei S, Chen P. Boundary condition: a key factor in tubular string buckling. SPE J 2015;20(6):1409–20. [16] van der Heijden GHM, Champneys AR, Thompson JMT. Spatially complex localisation in twisted elastic rods constrained to a cylinder. Int J Solids Struct 2002;39(7):1863–83. [17] Gulyayev VI, Gaidaichuk VV, Solovjov IL, Gorbunovich IV. The buckling of elongated rotating drill strings. J Petrol Sci Eng 2009;67(3):140–8. [18] Zdvizhkov A, Miska S, Mitchell R. Measurement and analysis of induced torsion in helically buckled tubing. SPE Drill Complet 2005;24(2):266–75. [19] Chen JS, Li HC. On an elastic rod inside a slender tube under end twisting moment. J Appl Mech 2011;78(4):041009. [20] Kwon YW. Analysis of helical buckling. SPE Drill Eng 1988;3(2):211–16. [21] Zhang YL. A study of helical buckling of the compressed portion of a drillstring during straight-hole drilling using a down-hole motor. Soc Pet Eng 1989:1–23 Document ID: SPE-19369-MS. Available from SPE, Richardson, Texas. [22] Liu F, Xu B, Gao D. Packer effect analysis of helical buckling of well tubing. J Tsinghua Univ 1999;8(39):105–8 (in Chinese). [23] Gao G, Miska S. Effects of friction on post-buckling behavior and axial load transfer in a horizontal well. SPE J 2010;15(4):1110–24. [24] Wang X, Yuan Z. Investigation of frictional effects on the nonlinear buckling behavior of a circular rod laterally constrained in a horizontal rigid cylinder. J Petrol Sci Eng 2012;90–91:70–8. [25] Gao G, Di Q, Miska S, Wang W. Stability analysis of pipe with connectors in horizontal wells. SPE J 2012;17:931–41. [26] Huang W, Gao D. Helical buckling of a thin rod with connectors constrained in a cylinder. Int J Mech Sci 2014;84:189–98. [27] Mitchell RF, Miska S. Helical buckling of pipe with connectors and torque. SPE Drill Complet 2006;21(2):108–15. [28] Gulyayev VI, Gaidaichuk VV, Andrusenko EN, Shlyun NV. Critical buckling of drill strings in curvilinear channels of directed bore-holes. J Pet Sci Eng 2015;129:168–77. [29] Lubinski A, Althouse WS, Logan JL. Helical buckling of tubing sealed in packers. J Pet Technol 1962;14(6):655–70. [30] Yue Q, Liu J, Dong R. The research of post-buckling about slender rod string in wellbore based on energy method and experiment. J Pet Sci Eng 2017;156:732– 739. [31] Hajianmaleki M, Daily JS. Critical-buckling-load assessment of drillstrings in different wellbores by use of the explicit finite-element method. SPE Drill Complet 2014;29(2):256–64. [32] Liang Z, Zhu ZL. Critical helical buckling load assessment of coiled tubing under axial force by use of the explicit finite-element method. J Pet Sci Eng 2018;169:51–7. [33] Liu JP, Zhong XY, Cheng ZB, Feng XQ, Ren GX. Buckling of a slender rod confined in a circular tube: theory, simulation, and experiment. Int J Mech Sci 2018;140:288–305. [34] Xiao J, Chen Y, Lu X, Xu B, Chen X, Xu J. Three dimensional buckling beam under cylindrical constraint. Int J Mech Sci 2019;150:348–55. [35] Wang C, Li Z, Li Y, Qi M, Wang L. Experimental study on the rotation of compression-buckling rod column in a liquid-filled cylinder. J Petrol Sci Eng 2018;164:459–66. [36] Mitchell RF. The pitch of helically buckled pipe. Soc Pet Eng 2005 SPE-146959-MS. [37] Timoshenko SP, Gere JM. Theory of elastic stability. McGraw-Hill; 1961. [38] Zienkiewicz OC, Taylor RL, Nithiarasu P, Zhu JZ. The finite element method. McGraw-hill; 1977. [39] Ekhande SG, Selvappalam M, Madugula MK. Stability functions for three-dimensional beam-columns. J Struct Eng 1989;115(2):467–79. [40] Waszczyszyn Z, Cichon C, Radwanska M. Stability of structures by finite element methods. Elsevier; 1994. [41] Zhang X, Liu J. Gap-element method of three-dimension contact problem for drillstring in horizontal well. Comput Struct Mech Appl 1992;9(4):405–10 (in Chinese). [42] Liu J, Ding H, Zhang X. Application of gap element to nonlinear mechanics analysis of drillstring. Zhejiang Univ-Sci 2002;3(4):440–4. [43] Otter JRH. Computations for prestressed concrete reactor pressure vessel using dynamic relaxation. Nucl Struct Eng 1965;1(1):61–75. [44] Alamatian J. Displacement-based methods for calculating the buckling load and tracing the post-buckling regions with dynamic relaxation method. Comput Struct 2013;114–115:84–97. [45] Rezaiee-Pajand M, Estiri H. Computing the structural buckling limit load by using dynamic relaxation method. Int J Non Linear Mech 2016;81:245–60.
Acknowledgments The authors gratefully acknowledge the financial support under Nanyang Technological University, Singapore’s Academic Research Fund (AcRF) Tier 1 Grant no. RG 168/16. This research is also supported by National Natural Science Foundation of China (Grant numbers: 11502051, 51607035, 51674088), China Postdoctoral Science Foundation (2018M641804), Heilongjiang Postdoctoral Research Foundation (LBH-Q18029), and Heilongjiang Youth Innovation Talents of Ordinary Undergraduate Colleges and Universities (UNPYSCT-2018046). References [1] He X, Kyllingstad A. Helical buckling and lock-up conditions for coiled tubing in curved wells. SPE Drill Complet 1955;10(1):10–15. [2] Miller JT, Su T, Dussan VEB, Pabon J, Wicks N, Bertoldi K, Reis PM. Buckling-induced lock-up of a slender rod injected into a horizontal cylinder. Int J Solids Struct 2015;72:153–64. [3] Hirose H, Yamane K, Marhefka G, Cavarocchi N. Right ventricular rupture and tamponade caused by malposition of the avalon cannula for venovenous extracorporeal membrane oxygenation. J Cardiothorac Surg 2012;7(1):36. [4] Bao X, Guo S, Xiao N, Li Y, Yang C, Shen R, Cui J, Jiang Y, Liu X, Liu K. Operation evaluation in-human of a novel remote-controlled vascular interventional robot. Biomed Microdevices 2018;20(2):34. [5] Li T. A mechanics model of microtubule buckling in living cells. J Biomech 2008;41(8):1722–9. [6] Afrin T, Kabir AM, Sada K, Kakugo A, Nitta T. Buckling of microtubules on elastic media via breakable bonds. Biochem Bioph Res Co 2016;480(1):132–8.
115
Contents lists available at ScienceDirect
International Journal of Mechanical Sciences journal homepage: www.elsevier.com/locate/ijmecsci
Post-buckling analysis of compressed rods in cylinders by using dynamic relaxation method Qiang Zhang a,b, Bao Jiang a, Zhongmin Xiao b,∗, Wei Cui a,∗, Jubao Liu a a b
School of Mechanical Science and Engineering, Northeast Petroleum University, Daqing 163318, PR China School of Mechanical and Aerospace Engineering, Nanyang Technological University, 50 Nanyang Avenue, 639798, Singapore
a r t i c l e
i n f o
Keywords: Post-buckling Transition section Finite element Gap element Dynamic relaxation method
a b s t r a c t Buckling of rods with cylindrical constraints is an essential problem in engineering and medical fields. Completely different to the buckling of classical Euler rods, the post-buckling modes for rods with cylindrical constraints involve highly complicated deformations and geometric configurations. In this paper, the rods are discretized into beam elements by finite element method, and the constraint relation between the rods and cylinders is described by gap element. A dynamic relaxation method for static buckling of compressed rods in cylinders is introduced. Numerical simulations have been carried out to investigate the effects of damping, element length, initial contact stiffness and penetration on post-buckling, etc. The characteristic buckling modes include deflection curves with single point, two points, three points and point-line-point contact. For the helical buckling, it is found that the transition section consists of two noncontact sections only, while the perturbed-helix section does not exist. Except the critical load and shear force of helical buckling, numerical simulation results are in good agreement with the analytical solution. By using the dynamic relaxation method, stable helical buckling modes can be obtained directly without undergoing sinusoidal buckling deformation under given compressive loads.
1. Introduction Slender rods in finite space exist in various petroleum engineering and medical applications. In petroleum engineering, the drill pipe in the borehole of several kilometers is subjected to great axial compressive force during drilling. When the axial force exceeds the threshold value, the drill pipe may bend into a three-dimensional shape, resulting in rapid increase in contact and friction force, and ultimately hinder its further drilling work, known as “lock-up” [1,2]. While in medical applications, stent therapy is often used for patients with coronary artery disease. It involves the insertion of a long flexible guide wire [3,4]. A great deal of effort is needed to drive the front end of the guide wire through a narrow lesion area with severe vascular occlusion. However, the guide wire may buckle due to large axial force. In addition, microtubules confined by microfilaments and cytoplasm around the protrusions of living cells may be helical, which may be caused by buckling [5,6]. Due to the importance of slender rods’ stability in cylinders, their buckling behavior has attracted extensive attention in the engineering field for decades. Lubinski [7] studied the two-dimensional buckling of drill pipe in wellbore for the first time and derived the calculation formula of the critical buckling load. ∗
Later, researchers extensively studied the buckling behavior of slender rods in cylinders from different aspects. For example, cylindrical constraint types include vertical [8,9], inclined and horizontal wellbores [2,10] in petroleum engineering. They also analyzed the effects of various factors on buckling modes and critical loads, including compression load [11–15], torque [16–19] gravity [9,20,21], boundary conditions [11,15,22], friction [23,24] and joints [25,26], etc. With their efforts, the following research methods have been developed: (1) theoretical methods such as beam-column model [7,26,27], buckling differential equation [17,28] and energy method [29–31]; (2) finite element methods such as implicit [9] and explicit [32–34] schemes; (3) experimental methods [31,35,36]. This paper focuses on the buckling of weightless compression rods in cylinders. Researchers have found that when the applied axial force reaches a certain critical value, the rods (restrained in cylinders) first buckle from the initial linear state into a sinusoidal shape. This phenomenon is called sinusoidal buckling. In the subsequent post-buckling process, with the increasing compressive load, more complex buckling configuration known as helical buckling occurs. Lubinski et al. [29] assumed that the tubing buckled into a space helix. Without considering the effect of boundary conditions, the mathematical relationship between the pitch and the axial load
Corresponding authors. E-mail addresses: [email protected] (Z. Xiao), [email protected] (W. Cui).
https://doi.org/10.1016/j.ijmecsci.2019.05.040 Received 9 November 2018; Received in revised form 10 May 2019; Accepted 26 May 2019 Available online 28 May 2019 0020-7403/© 2019 Published by Elsevier Ltd.
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
was derived by the energy method. The relationship is given as follows 𝐹hel =
8π2 𝐸𝐼 . 𝑝2
(1)
where Fhel is the axially critical load of helical buckling, EI is flexural rigidity of the tubing, p is the pitch of the buckled tubing. Mitchell [11] considered the boundary conditions effect and obtained another formula for the critical load. Under the clamped boundary condition, he calculated the length of the rod end to the helix section. He divided the helical buckling mode into two parts: the end transition section and the helical section. He pointed out that the end boundary conditions had an important effect on the buckling behavior. The critical load he obtained was 2.8 times that of the conventional analysis. Liu et al. [22] divided the transition section into a non-contact section and a perturbed-helix section. Sorenson et al. [12], Mitchell [36] and Huang et al. [15] divided the transition section into two no-contact sections and a perturbed-helix section. The two no-contact sections are called the first no-contact section and the second no-contact section respectively. The two no-contact sections were described by the beam-column model [37]. In the perturbed-helix section, the deformation curve was described by the buckling differential equation [22]. However, the above work considered only the final configuration of the buckled rod. They used different deflection assumptions to divide the ultimate deflection curve of the helical buckling into different sections and calculated the length of each section, but some of these length values are different, the corresponding critical loads are also different. In fact, there are some intermediate deformation stages between sinusoidal and helical shape. Based on the deflection curve hypothesis, the buckling modes and critical loads of sinusoidal to helical buckling were studied by using the small deformation hypothesis [13], the elastic model [14], the beam-column model and the differential equation [8]. In this paper, the assumption of deflection curve used in conventional theoretical analysis is abandoned. With the aid of finite element simulation, a dynamic relaxation method for post-buckling analysis of compressed rods in cylinders is proposed by combining the gap element with the beam element. We investigate the effects of damping, element size, contact stiffness and penetration on buckling. The critical conditions of various buckling modes from sinusoidal to helical buckling are analyzed. The relations between the buckling modes and shear force, bending moment, contact force and critical load are investigated. In the configuration of helical buckling, an exception that is inconsistent with the previous conclusions about the relation between the boundary condition and transition section is discovered, and a new critical load of helical buckling is obtained.
Fig. 1. Physical model of compressed rods in cylinders: (a) initial configuration of compressed rods in cylinders before buckling, (b) the deformation after buckling under the axial compression load Fc , (c) cross section before buckling, and (d) cross section after buckling. L is the rod length, 𝛽 G is the direction angle of the buckled rod with respect to global coordinates (X, Y, Z).
rigid body rotation, the rotation in the x-axis direction is restricted at the lower end. Under a given axial compressive load at the end of the rods, they will buckle and contact the inner wall of the rigid cylinders. The shape of the buckling deformation is shown in Fig. 1(b). The cross section before buckling is shown in Fig. 1(c), and the cross section after buckling deformation is shown in Fig. 1(d). The direction angle of deformation is 𝛽 d . 3. Beam element The finite element method is used to divide the compressive rod along the axis into several 3D straight beam elements. Bernoulli–Euler hypothesis is adopted. Fig. 2 shows the nodal displacement and nodal force of any 3D beam element. The node displacement vector is [ ]T 𝑞𝑒 = 𝑢𝑖 , 𝑣𝑖 , 𝑤𝑖 , 𝜃𝑖𝑥 , 𝜃𝑖𝑦 , 𝜃𝑖𝑧 , 𝑢𝑗 , 𝑣𝑗 , 𝑤𝑗 , 𝜃𝑗𝑥 , 𝜃𝑗𝑦 , 𝜃𝑗𝑧 . (2)
2. Physical model
where u is the axial displacement; v and w are the lateral displacement; 𝜃 is the angular displacement; i and j is the nodes of the beam element. The element displacement represented by nodal displacement is
The established physical model on buckling of compressed rods in cylinders is shown in Fig. 1, and the following assumptions are made: (i) The cross sections of the rods and cylinders are both circular. (ii) The initial annular gap between the rods and the cylinders exists. The axis of the rods before deformation coincides with the axis of the cylinders, and the axes are all straight lines. (iii) The rods are flexible/slender and are in a linear elastic state. The small deformation assumption is adopted. The annular gap between the rod and the cylinder is far lower than the axial length of the rod, so the lateral deformation of the rod is small. (iv) The cylinders are rigid. (v) After buckling, the rods will contact with the inner wall of the cylinders. The effect of friction is ignored. (vi) The ultimate buckling state of the rods subjected to a given axial compression load is studied. The transient effects of inertia force and damping force are neglected.
𝑢 = [𝑢, 𝑣, 𝑤, 𝜃]T = 𝑁 𝑞𝑒 ,
(3)
where N is the Hermitian shape function matrix given by Zienkiewicz et al. [38] ⎧ ⎫ ⎡ ⎪ 𝑁 𝑢 ⎪ ⎢𝑁 1 ⎪ 𝑁𝑣 ⎪ ⎢ 0 𝑁 =⎨ ⎬= ⎪ 𝑁𝑤 ⎪ ⎢⎢ 0 ⎪𝑁𝜃𝑥 ⎪ ⎣ 0 ⎩ ⎭
0 𝑁3 0 0
0 0 𝑁3 0
0 0 0 0 0 −𝑁4 𝑁1 0
0 𝑁4 0 0
where 𝑁1 = 1 − 𝑥𝑙 , 𝑁2 = 𝑥𝑙 , 𝑁3 = 1 − 3𝑥2
2𝑥3
𝑥2
𝑥3
𝑁2 0 0 0
3𝑥2 𝑙2
0 𝑁5 0 0
+
0 0 𝑁5 0
2𝑥3 , 𝑙3
0 0 0 0 0 −𝑁6 𝑁2 0
𝑁4 = 𝑥 −
0 𝑁6 0 0
⎤ ⎥ ⎥, ⎥ ⎥ ⎦
2𝑥2 𝑙
+
(4)
𝑥3 , 𝑙2
𝑁5 =
− 𝑙3 , 𝑁6 = − 𝑙 + 𝑙2 , l is the length of one beam element. The post-buckling deformation of the rods are subjected to both axial and bending actions. Based on Green–Lagrange strain and neglecting higher order terms, we get the generalized strain 𝜀 due to these actions 𝑙2
The two ends of the rods are hinged or clamped, and the mechanical model of the hinged boundary is shown in Fig. 1(a). To prevent 104
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
Fig. 2. Schematic diagram of beam element: (a) nodal displacements, and (b) nodal forces with respect to local coordinates (x, y, z). F and M is the force and moment, respectively.
[39]: ⎧ d𝑢 ⎧ ⎫ ⎪ d𝑥 + ⎪ 𝜀𝑥 ⎪ ⎪ ⎪𝜀 ⎪ ⎪ 𝜀 = ⎨ 𝑏𝑦 ⎬ = ⎨ ⎪𝜀𝑏𝑧 ⎪ ⎪ ⎪ 𝛾𝑥 ⎪ ⎪ ⎩ ⎭ ⎪ ⎩
1 2
(
d𝑣 d𝑥
)2 +
−𝑦 dd𝑥𝑣2 2
−𝑧 dd𝑥𝑤2 2
𝜌 dd𝑥𝜃
1 2
(
)2 ⎫ d𝑤 ⎪ d𝑥 ⎪ ⎪ ⎬ = 𝜀𝐿 + 𝜀𝑁𝐿 , ⎪ ⎪ ⎪ ⎭
of the cylinders, respectively. The shape is a ring with a certain thickness. The gap element has the following physical characteristics: (5)
(i) When the rods are not in contact with the inner wall of the cylinders, the compressive stiffness of the gap element approaches zero, which does not affect the free motion of the rods. (ii) When the rods are in contact with the inner wall of the cylinders, the gap element is completely compressed in a certain radial direction of the circumference, and the outer boundary of the gap element is tangent at the contact point. Its normal contact stiffness becomes quite large, which prevents the rods from intruding into the inner wall of the cylinders, but allows the rods to slide along the inner wall of the cylinders and produce contact forces, as shown in Fig. 3.
where 𝜀x , 𝜀by , and 𝜀bz are the axial strain along x axis, bending strain around y and z axis, respectively; 𝛾 x is the shear strain; 𝜌 is radial distance from the x axis of the rod. Generalized strain 𝜀 can be decomposed into linear strain 𝜀L and non-linear strain 𝜀NL . Matrix BL of kinematric relation corresponds to the linear part of −𝑁 ′′ 𝑣 −𝑁 ′′ 𝑤 𝑁 ′ 𝜃𝑥 ]T , the non-linear parts of Eq. (5): 𝐵𝐿 = [ 𝑁 ′ 𝑢 Eq. (5) are expressed by the vector of displacement gradients [40]: { } 𝑔̂ = 𝑣′ (6) 𝑤′ 𝜃𝑥 𝑣′′ 𝑤′′ 𝜃′ 𝑥 ,
The local coordinate system oyz of the gap element coincides with the local coordinate system of the beam element. The gap element is a twodimensional element whose displacements are the lateral displacements of the beam element. The displacement vector of the gap element at any position x of the ith node of the beam element is
and the shape function derivative matrix Γ is associated with the gradient of Eq. (6) { } 𝑁 ′𝑤 𝑁𝜃𝑥 𝑁𝑣 ′′ 𝑁𝑤 ′′ 𝑁 ′ 𝜃𝑥 . Γ = 𝑁 ′𝑣 (7) The linear stiffness matrix 𝑘𝑒0 and the geometric matrix 𝑘𝑒𝜎 have the following components [40]: 𝑘𝑒0
𝑙
=
∫0
𝑘𝑒𝜎 = 𝐹𝑐
𝐵𝐿 𝐷𝐵𝐿 d𝑥, T
𝑙
∫0
[ ]T 𝑢𝐺 = 𝑣𝑥 , 𝑤𝑥 = 𝑁𝐺 𝑞𝑒 ,
(11)
(8)
ΓT Γd𝑥,
(9)
where Fc the axial compression load, D is an elastic matrix and the expression is ⎡𝐸𝐴 ⎢ 𝐷=⎢ ⎢ ⎢ ⎣
𝐸 𝐼𝑧
𝐸𝐼𝑦 0
⎤ ⎥ ⎥, ⎥ 𝐺𝐼𝑝 ⎥⎦
(10)
in which E and G are, respectively, the elastic modulus and shear modulus of the beam; A is the cross-sectional area; Iz and Iy is second moments of area about the neutral axes z, y, respectively; Ip is the polar moment of inertia of area. 4. Gap element Fig. 3. Deformation of gap element: (a) initial position before buckling, and (b) contact between rods and the inner wall of the cylinders after buckling. FGN is the normal contact force, vx and wx are lateral displacement, 𝛽 G is the direction angle, D and d are the outer diameter of the rod and the inner diameter of the cylinder, respectively.
4.1. Physical properties and strains of gap elements A gap element [41,42] is a virtual element, its inner and outer boundaries are connected with the outer surface of the rods and the inner wall 105
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
𝑁𝑣 } is a 2 × 12 ma𝑁𝑤 trix, and the explicit expressions are given in Eq. (4). The lateral displacement UG and the direction angle 𝛽 G of the gap element at any position are described below: √ { 𝑈𝐺 = 𝑈𝑑 = 𝑣2𝑥 + 𝑤2𝑥 (12) ( ), 𝛽𝐺 = 𝛽𝑑 = 180𝑛 + arctan 𝑣𝑥 ∕𝑤𝑥
(i) The free state is in accordance with the following relations
where shape function matrix of gap element 𝑁𝐺 = {
⎧𝜀 − 1.0 < −𝜀𝑁,min ⎪| 𝐺𝑁 ⎨|𝐹𝐺𝑁 − 𝐹𝑁0 || ≤ 0.0 , ⎪𝐺𝑘 = 𝐺𝑘𝑗 ⎩
(ii) The elastic contact state is in accordance with the following relations ⎧−𝜀 ≤ 𝜀𝐺𝑁 − 1.0 ≤ 𝜀𝑁,max ⎪| 𝑁,min , ⎨|𝐹𝐺𝑁 − 𝐹𝑁0 || ≥ 0.0 ⎪𝐺𝑘 = 𝐺𝑘𝑐 ⎩
where Ud and 𝛽 d are the lateral displacement and the direction angle of the beam element, n is the half circle number, the initial gap rG of the gap element is expressed as 𝑟𝐺 = (𝐷 − 𝑑 )∕2,
(17)
where 𝜀N,min , 𝜀N,max and FN0 are small quantities larger than zero to control the convergence accuracy of gap element iteration. 𝜀N,min and 𝜀N,max are the minimum and maximum iteration errors for gap element strain, respectively. FN0 is the iterative error of contact force of gap element. In order to make the contact between the rods and the inner wall of the cylinders satisfy the definite solution of the gap element, the stiffness Gk of the gap element should be revised continuously in the iterative calculation. The following principles should be adopted in the correction:
(13)
where D and d are the outer diameter of the rod and the inner diameter of the cylinder, respectively. Thus, the normal strain 𝜀GN of the gap element is 𝜀𝐺𝑁 = 𝑈𝐺 ∕𝑟𝐺 .
(16)
(14)
4.2. Criterion and definite formula for the contact state of gap element One of the two states after the deformation of the rods is
(i) The initial contact stiffness Gk0 and the initial penetration UPe,0 of the gap element are given for correction in the iteration: { 𝐺𝑘0 = 𝑓Kn ⋅ 12𝑙𝐸𝐼 3 , (18) 𝑈Pe,0 = 𝑓Pe ⋅ 𝑑∕2
(i) Free state: the rods are separated from the inner surface of the cylinders and do not contact with the wall of the cylinder. (ii) Elastic contact: the rods are in contact with the inner wall of the cylinders, a small amount of penetration between them exists. The contact force is linear with penetration.
where fKn is the coefficient of contact stiffness, and fPe is the penetration coefficient.
As shown in Fig. 4, the criterion for describing the above state by gap element is (1) free state: 𝜀GN < 1.0, GK = Gkj ,SG1 ∈ S; (2) elastic contact: 𝜀GN ≥ 1.0, GK = Gkc ,SG2 ∈ S, here SG1 and SG2 represent the region where the gap element is located, S represents the contact region. The relation between the stiffness Gk of the gap element and the strain 𝜀GN is shown in Fig. 4, in which Gk0 is the initial value of the contact stiffness of the gap element, Gkj is the compressive stiffness of the gap element, and Gkc is the normal contact stiffness of the gap element after the contact between the rods and the inner wall of the cylinders. Obviously, the gap element is actually a nonlinear element with variable stiffness. The contact force between the rode and the inner surface of the cylinder is { , 𝑈𝐺 − 𝑟𝐺 < 0 0 ( ) 𝐹𝐺𝑁 = . (15) 𝐺𝑘 𝑈𝐺 − 𝑟𝐺 , 𝑈𝐺 − 𝑟𝐺 ≥ 0
(i) In the free state, the compressive stiffness Gkj of the gap element should be minimized, so that the free movement of the rod will not be affected. (ii) For elastic contact, the stiffness of gap element Gk is related to the bending stiffness of the rods. (iii) The contact between the rods and the cylinder wall is multidirectional. The stiffness of the gap elements between the adjacent beam elements is mutually affected among the gap elements. Therefore, this effect should be taken into account when correcting the stiffness of the gap element. Each gap element should be corrected repeatedly until all the gap elements satisfy the definite solution. The introduction of gap element enables the connection of beam elements and cylinders into a continuum.
where FGN is the contact force, Gk is the stiffness of the gap element. When contact occurs, the small penetration is UPe = UG − rG . The contact force FGN can be further transformed to the additional nodal force of the beam element. Based on the strain 𝜀GN in Eq. (14) and contact force FGN in Eq. (15), the definite formula of iteration calculation for gap element can be written as:
4.3. Static buckling equilibrium equation Combing of beam element and gap element, the equilibrium equation is (𝑘𝑒0 + 𝑘𝑒𝜎 + 𝑘𝑒𝐺 )𝑞𝑒 = 𝐹e + 𝐹𝐺𝑒 ,
(19)
Fig. 4. Relation between strain and stiffness of gap element. (a) The normal contact stiffness less than the initial contact stiffness: Gkc < Gk0 , and (b) the normal contact stiffness greater than the initial contact stiffness: Gkc > Gk0 . Gkj is the compressive stiffness of the gap element.
106
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
where Fe and 𝐹𝐺𝑒 is the equivalent nodal force of beam and gap element, respectively; 𝑘𝑒𝐺 is the gap stiffness matrix given by 𝑘𝑒𝐺 = 𝑁𝐺T 𝐺𝑘 𝑁𝐺 .
The dynamic relaxation method is a numerical method which transforms the static problem into a quasi-dynamic problem according to the D’Alembert principle and finds the equilibrium form of the structure through parameter iteration. The structure is discretized into elements with certain virtual mass. Under the action of unbalanced force, these discrete elements will inevitably move along the direction of unbalanced force and reduce the overall unbalanced force of the structure. It also discretizes the vibration process of a spatially discretized structural system in time and traces the vibration process of the system step by step. By introducing the dynamic damping to consume the oscillation energy, the kinetic energy of the structure motion is gradually reduced and tends to zero, that is, to reach the static equilibrium state. Dynamic relaxation method is also widely used in static buckling analysis of structures [44,45]. However, the post-buckling analysis of cylindrically constrained rods has rarely been reported by using this dynamic relaxation method.
(20)
After coordinates transformation and assembling of all beam elements and gap elements, we obtain the global equilibrium equation for the geometric and contact nonlinear static analysis of the rods in cylinders: (𝐾0 + 𝐾𝜎 + 𝐾𝐺 )𝑑 = 𝐹b + 𝐹𝐺 ,
(21)
where K0 , K𝜎 and KG are, respectively, the global linear stiffness matrix, geometric stiffness matrix and gap stiffness matrix of the rods; d is the node displacement vector; Fb and FG are, respectively, the equivalent nodal force vector of beam element and gap element with respect to global coordinates. The nonlinear Eq. (21) shows that with the increase of axial compressive load Fc , the geometric stiffness matrix KG decreases, resulting in the rod buckling. The lateral displacement of the rods suddenly increases, leading to the contact between the rods and the inner cylinder wall. In solving the nonlinear Eq. (21), it is necessary to distinguish and modify the definite solution of the gap element in the same iteration step. Only when all these conditions are satisfied, the generalized displacement and internal force can converge to the exact solution.
5.3. Dynamic relaxation method for static buckling of compressed rods in cylinders In this paper, with the effect of cylindrical constraint, the dynamic relaxation method is applied to solve the static buckling problem of compressed rods in cylinders. The static buckling model of Fig. 1(a) is transformed into the dynamic one shown in Fig. 5. In order to reduce the impact of axial compression load on the rods, the sinusoidal wave loading is adopted, and the time history of axial compression load is as follows { 𝐹c sin (π𝑡∕𝑇 ), 0 ≤ 𝑡 ≤ 𝑡1 𝐹 c (𝑡 ) = , (22) 𝐹c , 𝑡 > 𝑡1
5. Dynamic relaxation method 5.1. Problems in static buckling calculation Eq. (21) is a complex equation involving the dual nonlinearity of geometry and contact. For static buckling calculation, some problems such as convergence difficulty and algorithm instability still exist. They are manifested in the following aspects:
where the axial compression load 𝐹c = 𝑘 𝜋 𝐿𝐸𝐼 2 , k is a dimensionless compressive load; T is the characteristic time for loading; t is the response time; t1 = T/4. The initial geometric imperfections will affect the critical buckling load. The small perturbation load is initially applied and revoked in subsequent analysis. Thus, the perturbation load will not affect the critical buckling load. 2
(i) In the process of buckling, the jumping change between different buckling modes always exists. The strong nonlinear mechanical behaviors such as contact and separation between the rods and the cylindrical wall cause a sudden change in the contact state in the calculation process. These physical phenomena are easy to cause abnormal termination of calculation and have difficulty in convergence. Especially in the critical buckling state, the deflection suddenly increases and the contact state cannot be recognized. The rods penetrate through the cylindrical wall, causing the calculation to fail. (ii) The algorithm instability shows that the convergence solution is not unique, and the configuration of the rod buckling is arbitrary. For example, when the plane sinusoidal buckling mode is formed, this mode rotates in the cylinders with the increase of axial compression load, but the helical buckling mode cannot be formed. Because of the sensitivity to iterative load step, buckling modes can sometimes be formed, but sometimes they cannot. After the helical buckling mode is formed, the buckling modes of the right and the left helix may exist simultaneously at different axial positions. To solve the problems of convergence difficulty and algorithm instability in the buckling of rods in cylinders, a dynamic relaxation method is proposed. 5.2. Basic idea of dynamic relaxation method The basic idea of dynamic relaxation [43] originates from the steadystate solution of a single degree of freedom damping system. Large damping is used to attenuate the dynamic response of the system, and the static part of the transient dynamic solution (i.e., the steady state) is found.
Fig. 5. Mechanical model based on dynamic relaxation method: (a) the position where these loads [axial compression force Fc (t) and perturbation loads Fp,i (t)] are applied, (b) the plane projection of the perturbation loads, and (c) the curve of the applied loads vs. time. 107
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
In order to easily form the spatial deformation, lateral perturbation loads Fp,1 (t), Fp,2 (t) and Fp,3 (t) are applied at the positions of 0.25L intervals along the rod length: { 𝐹p,𝑖 ⋅ sin (2π𝑡∕𝑇 ), 0 ≤ 𝑡 ≤ 𝑡2 𝐹p,𝑖 (𝑡) = , (23) 0, 𝑡 > 𝑡2 where i = 1, 2, 3, t2 = T/2. The directions of perturbation loads are shown in Fig. 5(a) and (b). Considering the inertia force and damping force in the dynamic model, the perturbation load and axial compression load are applied in Eq. (21). The basic dynamic buckling equation of the compression rods in cylinders is obtained: 𝑀𝑑 ′′ + 𝐶𝑑 ′ + 𝐾𝑑 = 𝐹 (𝑡),
(24) Fig. 6. Effect of compressive load on the lateral displacement Ud of the compression rod without cylindrical constraints (𝛼 = 𝛼 0 , 𝛼 0 = 2𝜔).
where M, C and K are, respectively, the mass, damping and stiffness matrix of the compression rods, K = K0 + K𝜎 + KG ; d′′ and d′ are node acceleration and velocity vectors of the rods, respectively; F(t) is the load vector of the rods, F(t) = Fb (t) + FG (t). In this paper, Rayleigh damping is used to set the Rayleigh damping as a large value to attenuate the dynamic response of the system. Considering the effect of 𝛼 damping in the Rayleigh damping, then the damping matrix C = 𝛼M. The damping force 𝐶 𝑑 ′ = α𝑀𝑑 ′ . Thus, the faster the movement speed of the rods, the greater the damping force. The bigger damping force can effectively control the jumping change of buckling mode and the sudden change of contact state, which is helpful to improve the calculation efficiency and stability. Newmark numerical integration method is used to calculate Eq. (24). When t > t2 , the perturbation load is equal to 0, the axial compression load remains unchanged. Let the left item of the static buckling Eq. (21) be ( ) (25) 𝐾0 + 𝐾𝜎 + 𝐾𝐺 𝑑 = 𝐹k ,
Fig. 7. Effect of 𝛼 damping on the lateral displacement Ud of the compression rod without cylindrical constraints (k = 0.998, 𝛼 0 = 2𝜔).
and the right item of Eq. (21) be 𝐹b + 𝐹𝐺 = 𝐹J .
6.1. The effect of 𝛼 damping without cylindrical constraints
(26)
The kinetic Eq. (24) is rewritten into 𝑀𝑑 + 𝐶𝑑 = 𝐹J (𝑡) − 𝐹k (𝑡) = 𝐹R (𝑡), ′′
′
𝜋2 𝐿2
(27)
The natural frequency of the hinged rod at both ends is 𝜔 = √ 𝐸𝐼 = 80.40 rad. The reference 𝛼 damping is 𝛼 0 = 2𝜔. Figs. 6 and 𝜌𝐴
7 show the effect of the compression load 𝑘 = 𝐹c ∕( 𝜋 𝐿𝐸𝐼 2 ) and 𝛼 damping on the lateral displacement of the middle of the compression rod, respectively. As can be seen from Fig. 6, when k<0.998, the lateral displacement Ud decreases with the withdrawal of the perturbation load, and tends to zero gradually, so the rod is in a stable state. When k = 0.998, Ud remains unchanged as perturbation load is removed. The compressed rod is in a critical state of equilibrium, and the corresponding load 𝐹c = 2
where FR (t) is the residual force vector. With the extension of calculation time, the axial compression load Fc remains unchanged (when t > t1 ), and the perturbation load is 0 (when t > t2 ). Due to the existence of damping, the transient response of the rods tends to be static, and the rod deformations d will also tend to be constant after buckling. Residual force vector FR (t) containing inertia force Md′′ and damping force Cd′ decreases to 0, then Fk = FJ , the contact force FG is no longer changing. As a result, the solution of nonlinear static buckling of rods in cylinders is achieved.
0.998 𝜋 𝐿𝐸𝐼 is the critical value. When k>0.998, Ud increases with the 2 withdrawal of perturbation load, so the rod is in an unstable state. The dimensionless critical load (k = 0.998) calculated by us is 0.2% less 2
6. Analysis of effect factors on buckling mode of compressed rods
than the classical Euler formula 𝐹c = 1.000 𝜋 𝐿𝐸𝐼 2 , which indicates that the dynamic relaxation method has a high calculation accuracy. From Fig. 7, with different 𝛼 damping, the lateral displacement is still in the equilibrium state with the withdrawal of the lateral force. The larger the 𝛼 damping is, the smaller the ultimate displacement will be. Therefore, the 𝛼 damping has no effect on the critical load of the compression rod. It only affects the displacement of the ultimate equilibrium. 2
Firstly, the effect of 𝛼 damping on buckling of a compression rod is analyzed without the constraint of the cylinder. Then, with the cylindrical constraints, the effects of 𝛼 damping, element length, contact stiffness and penetration on post-buckling of the compression rod are analyzed under a given load. The calculation parameters of the rod are as follows: the length L = 0.864 m, diameter d = 4.76E−3 m, density is 7850 kg/m3 , elastic modulus E = 204.78 GPa. The inner diameter of the cylinder is D = 38.1E−3 m. The initial gap between the cylinder and the rod is rG = 16.67E−3 m. The beam element length is l = 8E−3 m. The response time t = 20 s, and the characteristic time for loading T = 5 s. Lateral perturbation load Fp,1 = Fp,2 = Fp,3 = 1 N. Compressive stiffness of gap element Gkj = 0. The minimum iterative error of gap element strain is 𝜀N,min = 1E−5. The maximum value of permissible penetration UPe,max ≤ 1E−3 m. The maximum iteration error of gap element strain resulted is 𝜀N,max = 0.06.
6.2. The effect of 𝛼 damping with cylindrical constraints Given the dimensionless compression load (k = 30), the dynamic relaxation method is used to calculate the post-buckling of the compression rod under different 𝛼 damping conditions. Fig. 8 shows the effect of the 𝛼 damping on lateral displacement, velocity and acceleration at the 0.25L position. For 𝛼 = 300𝛼 0 , Fig. 9 shows the buckling process of the post-buckling mode. 108
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
Fig. 8. The effect of 𝛼 damping on (a) lateral displacement Ud , (b) velocity Ud ′ and (c) acceleration U′′d at the 0.25L position (𝛼 0 = 2𝜔). Fig. 9. Stabilization process of post-buckling mode of the rod with cylindrical constraints (k = 30, 𝛼 = 300𝛼 0 ): (a) 3D deformation curve, and (b) top view of the deformation at certain times.
In Fig. 8, during the loading stage of axial compression and lateral perturbation load, the lateral displacement Ud increases with time. When Ud reaches the annular gap (rG = 0.01667 m), the rod contacts with the inner wall of the cylinder, and the lateral displacement remains unchanged. The larger the 𝛼 damping, the longer the time needed to reach the contact between the rod and the cylindrical wall. If the 𝛼 damping is too small, it will lead to excessive lateral velocity and acceleration of the rod, and produce a greater impact vibration. In the stable stage of axial compression load, the lateral displacement Ud tends to be constant, and the lateral velocity and acceleration tend to be zero. The lateral velocity and acceleration can be restrained by choosing a larger 𝛼 damping. Under the condition of long response time, the effect of lateral velocity and acceleration is eliminated, and the residual force FR in Eq. (27) is reduced to 0. As a result, a stable buckling deformation can be obtained. In Fig. 9, during the loading stage of axial compression load and lateral perturbation load (t<5 s), the compression rod deforms into a three-dimensional deflection curve, and the contact state of two points (t = 3.25 s), three points (t = 3.5 s) and point-line-point (t = 4.25 s) with the cylinder wall gradually occurs. In the stable stage of compressive load, the buckling mode does not change with time when t ≥ 10 s. Therefore, the deformation of the rod does not undergo sinusoidal buckling mode, the stable and helical post-buckling modes can be directly calculated by using the dynamic relaxation method with its loading style of axial compression and perturbation load.
flection curve. When the element length is less than or equal to 8E−3 m, i.e., when the element number is greater than or equal to 108, both the lateral displacement curves and the angular displacement curves coincide with each other. 6.4. The effect of contact stiffness and penetration With the given 𝛼 damping (𝛼 = 300𝛼 0 ) and compressive load (k = 30), the effect of the initial contact stiffness factor fKn and penetration factor fPe of the gap element on the calculation results is given in Fig. 11. From Fig. 11, after nonlinear iterative calculation, the normal contact stiffness Gkc and penetration UPe are both related to fKn and fPe . Gkc is directly proportional to the initial contact stiffness factor fKn , and is inversely proportional to the penetration factor fPe . UPe is inversely proportional to fKn , and is proportional to fPe . For fKn = 0.01 and fPe = 1.0, the value of UPe is close 1E−3 m, reaching the maximum value. When fKn is very small, the contact force FGN is inversely proportional to fPe . When fKn ≥ 0.3, FGN is no longer changing with fKn , fPe , and UPe after the iteration is less than 0.1E−3 m. 7. The analysis of characteristic buckling mode 7.1. Characteristic buckling modes and critical loads The dynamic relaxation method is used to carry out the post-buckling finite element calculation of the compressive rods. Various characteristic buckling modes and their contact positions obtained by varying the compressive load are shown in Fig. 12, and the corresponding deflection curves are shown in Fig. 13. The point contact in the pictures is represented by the first contact points c1 and pentagram. The end of the
6.3. The effect of element length Under the given 𝛼 damping (𝛼 = 300𝛼 0 ) and compressive load (k = 30), Fig. 10 shows the effect of beam element length on the de109
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
Fig. 10. The effect of element length on (a) lateral displacement Ud , and (b) angular displacement 𝛽 d . The element length from l1 to l8 is 2E−3 m, 4E−3 m, 6E-3 m, 8E−3 m, 12E−3 m, 16E−3 m, 24E−3 m and 32E−3 m, respectively, i.e., the rods are divided into 432, 216, 144, 108, 72, 54, 36 and 27 beam elements, respectively.
Fig. 11. The effect of contact stiffness and penetration (𝛼 = 300𝛼 0 , k = 30): (a) initial contact stiffness factor fKn = Gk0 / Gkc , (c) penetration UPe , and (d) contact force FGN .
continuous contact section is represented by the second contact points c2 and solid dots. According to the contact state between the compression rods and the cylinders, the characteristic buckling modes with single point contact (k = 1.00), two points contact (k = 4.00), three points contact (k = 7.83) and point-line-point contact are given in Figs. 12 and 13, respectively. In the buckling configuration of point-line-point contact, they are divided into just point-line-point deformation (k = 8.83), helix pitch measured by contact point c1 (k = 20.43) and continuous contact section (k = 30.16), respectively.
•
(12EI/l3 ), (b) normal contact stiffness
The characteristics of various buckling modes are as follows: (i) For buckling modes with single point contact (1.00≤ k ≤ 3.99), the contact position is located in the middle of the compressive rod. The value 𝜉 c1 = 0.5. (ii) For buckling modes with two points contact (4.00≤ k ≤ 7.82), the middle part of the rod is detached from the wall of the cylinder, and the contact points on both sides of the middle part appear. With the increase of compressive load, the two contact points gradually move along to the rod end. The position of the contact 110
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
Fig. 12. Characteristic buckling modes and their contact positions (𝜉 = x/L) for hinged boundary: (a) buckling modes, and (b) location of the contact point. Mode A: single point contact; Mode B: two points contact; Mode C: three points contact; Mode D: just point-line-point contact; Mode E: helix pitch measured by contact point c1; and Mode F: helix pitch measured by continuous contact section. The last three buckling modes are point-line-point contact deformation.
Fig. 13. Deflection curves of buckling modes for hinged boundary: (a) lateral displacement, and (b) direction angle vs. axial distance, respectively.
111
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
Table 1 Dimensionless shear force, bending moment and contact force in six post-buckling modes. Boundary
Post-buckling mode
k
qmax
q|𝜉 = 0.5
mmax
m|𝜉 = 0.5
n1
nh,max
nh |𝜉 = 0.5
Hinged
Mode A Mode B Mode C Mode D Mode E Mode F
1.00 4.00 7.83 8.83 20.43 30.16
0.936 1.033 1.732 1.741 1.756 1.715
0.006 0.450 1.056 1.072 1.043 1.017
0.981 0.722 1.214 1.218 1.230 1.213
0.981 0.048 0.935 0.945 1.037 1.030
0.007 0.903, 0.903 1.671, 0.003, 1.671 1.708 1.758 1.773
– – – 0.149 1.611 1.375
– – – 0.014 1.021 1.063
Clamped
Mode A Mode B Mode C Mode D Mode E Mode F
4.00 8.19 12.00 12.49 28.16 36.67
1.412 1.149 1.425 1.427 1.440 1.443
0.013 1.025 1.045 1.049 1.015 1.019
0.998 0.812 1.007 1.009 1.018 1.023
0.998 0.411 0.902 0.909 1.016 1.023
0.004 0.410, 0.410 1.020, 0.033, 1.020 1.039 1.102 1.108
– – – 0.223 1.029 1.055
– – – 0.223 1.029 1.055
point c1 decreases to 𝜉 c1 = 0.278, and the angle between the two contact points c1–c1 increases to 𝛽 d,c1 − c1 = 135.8°. (iii) For buckling modes with three points contact (7.83≤ k ≤ 8.82), the middle part of the rod contacts with the cylinder wall again. The contact point c1 decreases to 𝜉 c1 = 0.267, and the angle increases to 𝛽 d,c1 − c1 = 152.2°. (iv) For buckling modes with point-line-point contact, the contact point in the middle of the rod is changed into continuous contact section. For helix pitch measured by contact point c1, contact points c1 and c2 are located at 𝜉 c1 = 0.171 and 𝜉 c2 = 0.302, respectively. The values 𝜉 c1,c1 = 0.657, 𝛽 d,c1 − c1 = 360°, 𝛽 d,c2 − c2 = 228°. For helix pitch measured by continuous contact section, the values 𝜉 c1 = 0.142, 𝜉 c2 = 0.246, 𝜉 c2,c2 = 0.508, 𝛽 d,c1 − c1 = 360° +133.6°, 𝛽 d,c2 − c2 = 360°.
middle of the rods, except for the first two modes of hinged boundary and the first mode of clamped restraint, the shear force q|𝜉 = 0.5 for the other modes is about 1.0. The maximum bending moment mmax in mode A with hinged boundary is in the middle of the rod, and the mmax in the other modes are in the suspended section, as shown in Fig. 14(c), (d) and Table 1. mmax is 1.0 and 0.7, respectively, in the first two post-buckling modes, and the values in the others is about 1.2. For clamped boundary, the maximum bending moment exists at the rod ends, and the maximum values also exist in the suspended section and in the middle of the rod. mmax is 0.8 for the mode B, and is about 1.0 for other post-buckling modes. Regardless of the boundaries, m|𝜉 = 0.5 is about 1.0 for the last four modes. The contact force n1 at the first contact point c1 increases gradually in the first three modes as shown in Table 1. The value n1 increases from a small value in mode A to 0.9 and 0.4 in Mode B, and finally to 1.7 and 1.0 in Mode C for hinged and clamped boundary, respectively. The latter three modes are the point-line-point contact deformations, the contact force n1 increase slightly, and become 1.8 and 1.1. The contact force nh in the continuous contact section is not uniform along the rod, as shown in Fig. 14(e), (f) and Table 1. The maximum values exist at the ends of the continuous contact section for hinged boundary; however, the maximum values are in the middle of the continuous contact section for clamped boundary. The contact force distribution is more uniform in the middle of the continuous contact section. The contact force nh |𝜉 = 0.5 in the middle is about 1.1 for the helical buckling mode (mode F). It is emphasized that the mechanical parameters in Fig. 14 and Table 1 are dimensionless values, and the actual values (in SI unit) can be obtained by introducing the values in Table 1 into Formula (28). For example, the contact force F1 in mode E are respectively 237.9 N and 241.3 N for hinged and clamped boundary; the contact force Fh in the middle of the continuous contact section are 1605.4 N/m and 3074.0 N/m, respectively. In mode F, F1 are 430.3 N and 360.5 N, Fh are 3642.6 N/m and 5344.3 N/m, respectively.
It should be emphasized that in Fig. 12 we call the section from the first contact point to the end the first suspended section, Lc1 = 𝜉 c1 L; the section from the first contact point to the second contact point is referred to as the second suspended section, Lc2 = (𝜉 c2 − 𝜉 c1 )L; the length of pitch measured by contact point c1 is Lh1 = 𝜉 c1,c1 L; the length of pitch measured by contact point c2 is Lh2 = 𝜉 c2,c2 L. 7.2. Mechanical parameters of various post-buckling modes Introduce dimensionless shear force q, bending moment m and contact forces n1 , nh as √ 𝐹c ⎧𝐹 = 𝑞 𝐹 c 𝑟 𝐺 𝑠 2 2𝐸𝐼 ⎪ 𝐹c 𝑟𝐺 ⎪ ⎪𝑀 = 𝑚 2 √ , (28) ⎨ 𝐹c ⎪𝐹 1 = 𝑛 1 𝐹 c 𝑟 𝐺 2 2 𝐸𝐼 ⎪ ⎪𝐹 = 𝑛 𝐹𝑐2 𝑟𝐺 ℎ 4𝐸𝐼 ⎩ h where Fs is the shear force (in N); q is the dimensionless shear force; M is the bending moment (in N•m); m is the dimensionless bending moment; F1 is the contact force of the contact point c1 (in N); n1 is the dimensionless contact force of the contact point c1; Fh is the contact force of the continuous contact section, N/m; nh is the dimensionless contact force of the continuous contact section. In post-buckling modes for hinged and clamped boundaries, the distribution curves of dimensionless shear force, bending moment and contact force in continuous contact section are shown in Fig. 14, and the dimensionless shear force, bending moment and contact force of point and line contact are shown in Table 1. From Fig. 14(a) and (b) and Table 1, the maximum shear force qmax of the first two post-buckling modes is about 1.0 and that of the last four post-buckling modes is about 1.7 for hinged boundary. For clamped boundary, the maximum shear force is in the suspension section, qmax is about 1.1 for the mode B, and qmax is about 1.4 for the others. In the
7.3. Verification of critical load for various post-buckling modes In Table 2, the dimensionless critical loads of the characteristic buckling modes of the hinged and clamped boundary conditions are given and compared with the results in the literatures. For the hinged boundary, the critical loads of the first two modes in this paper are similar to those given in open literature [8], except the critical value of helical buckling (mode F). For the clamped boundary, the critical load obtained in this paper is slightly different from that obtained by the elastic model [14], and is in good agreement with that obtained by the small deformation hypothesis [13]. It should be noted that more critical loads of buckling modes are calculated in this paper than in the literatures. 112
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
Fig. 14. Non-dimensional load curves for post-buckling modes: (a) and (b) dimensionless shear force q vs. x/L, (c) and (d) dimensionless moment m vs. x/L, (e) and (f) dimensionless contact force nh vs. x/L for continuous contact sections.
Table 2 Dimensionless critical loads of the characteristic buckling modes. Post-buckling mode
Hinged boundary
Clamped boundary
Present
Huang et al. [8]
% Difference
Present
Chen and Fang [13]
% Difference
Fang et al. [14]
% Difference
Mode A Mode B Mode C Mode D Mode E Mode F
1.00 4.00 7.83 8.83 20.43 30.16
1.00 4.00 – 9.09 – 34.15
0 0 – −2.86 – −11.68
4.00 8.19 12.00 12.49 28.16 36.67
4.00 8.18 11.92 12.52 – –
0 0.12 0.67 −0.24 – –
4.00 8.08 10.80 12.03 – –
0 1.36 11.11 3.82 – –
113
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
Fig. 15. Deflection curves of helical buckling modes (mode F in Fig. 12) for (a) hinged boundary (k = 30.16, 𝜂 1 = 1.735, 𝜂 2 = 1.269, 𝜂 h = 6.192, 𝜂 L = 12.200), and (b) clamped boundary (k = 36.67, 𝜂 1 = 2.709, 𝜂 2 = 0.888, 𝜂 h = 6.258, 𝜂 L = 13.452) conditions. 𝜂 1 , 𝜂 2 , 𝜂 h and 𝜂 L is the dimensionless length of the first non-contact √ 2 𝐹c section, the second non-contact section, the continuous contact section and the compressed rod. 𝑘 = 𝐹c ∕( 𝜋 𝐿𝐸𝐼 ), 𝜂 = ⋅ 𝑥. 2 2𝐸𝐼 Table 3 Length of transition section in helical buckling modes (mode F in Fig. 12). Boundary
Author
𝜂1
𝜂2
𝜂1 + 𝜂2
𝜂3
𝜂1 + 𝜂2 + 𝜂3
Hinged
Sorenson and Cheathamjr [12], Mitchell [36], Huang et al. [15] This paper
1.747 1.735
1.347 1.269
3.094 3.004
1.677 0
4.771 3.004
Clamped
Mitchell [11] Liu et al. [22] Sorenson and Cheathamjr [12], Mitchell [36], Huang et al. [15] This paper
3.178 2.718 2.700 2.709
0 0 0.903 0.888
3.178 2.718 3.603 3.597
0 1.211 1.846 0
3.178 3.929 5.449 3.597
Table 4 Relevant mechanical parameters of helical buckling modes (mode F in Fig. 12). Boundary
Author
𝑚t _max
mh
q1
n1
nh
Hinged
Sorenson and Cheathamjr [12], Huang et al. [15] This paper
1.198 1.213
1.000 1.030
2.870 1.715
1.734 1.773
1.000 1.063
Clamped
Liu et al. [22] Mitchell [36] Sorenson and Cheathamjr [12], Huang et al. [15] This paper
1.000 0.766 1.010 1.020
1.000 – 1.000 1.023
0.543 0.457 0.342 1.443
1.253 – 1.069 1.108
– – 1.000 1.055
eratures [12,15,36]. The dimensionless total length 𝜂 L in this paper is 𝜂 L = 2(𝜂 1 + 𝜂 2 ) + 𝜂 h . According to Formula (22) and Formula (29), the relation between dimensionless total length 𝜂 L and dimensionless load k is
7.4. Verification of the lengths and relevant mechanical parameters of helical buckling modes Introduce dimensionless length 𝜂 √ 𝐹c ⋅ 𝑥, 𝜂= 2𝐸𝐼
𝑘 = 2𝜂𝐿2 ∕𝜋 2 ,
(29)
(30)
therefore, the critical load of helical buckling is less than that obtained by Sorenson and Cheathamjr [12], Huang et al. [15], and Mitchell [36]. Table 4 gives the relevant mechanical parameters of helical buckling modes and the results from literatures. In Table 4, 𝑚t _max is the maximum bending moment in transition section, mh is the bending moment in helical section, q1 is the maximum shear force in the transition section, n1 is the contact force at the contact point c1, nh is the contact force at the helical section, respectively. The values 𝑚t _max , mh , n1 and nh are in good agreement with the results of Sorenson et al. [12] and Huang et al. [15]. However, the maximum shear force in the transition section q1 is quite different from that in the literature. The value q1 of the clamped boundary is higher than those in the literatures [12,15,22,36], while this value of the hinged boundary is lower than those in literatures [12,15].
Under hinged and clamped boundary conditions, the deflection curves of helical buckling (mode F) are shown in Fig. 15. The dimensionless length of transition section is given and compared with the literatures results in Table 3. Both the first and second non-contact sections 𝜂 1 and 𝜂 2 exist in the helical buckling configuration, which are basically consistent with the results in the literatures [12,15,36]. In Fig. 15, the slope of angle 𝛽 d increases gradually in the first noncontact section, and tends to be constant in the second non-contact section, and finally remains unchanged in the continuous contact section. On the other hand, the length 𝜂 h of the continuous contact section calculated in this paper is about 2𝜋, which is basically consistent with the analytical results of the continuous contact hypothesis in the literature [29]. Therefore, the calculated transition section includes only the two non-contact sections 𝜂 1 and 𝜂 2 , does not have perturbed-helix section 𝜂 3 . Since there is no perturbed-helix section in this paper, the total length of the transition section (𝜂 1 + 𝜂 2 ) in this paper is significantly less than that of the three-section model (𝜂 1 + 𝜂 2 + 𝜂 3 ) used in the lit-
8. Conclusion Considering the geometric nonlinearity of the compressed rods and the contact nonlinearity with the external cylinder, the rods are discretized into beam elements, and the gap element is constructed by the 114
Q. Zhang, B. Jiang and Z. Xiao et al.
International Journal of Mechanical Sciences 159 (2019) 103–115
finite element method. To solve the problems of convergence difficulty and algorithm instability in static buckling analysis, a dynamic relaxation method for static buckling analysis of rods in cylinders is introduced. Based on the classical Euler formula, the effectiveness of dynamic relaxation method in calculating the static buckling of the compressed rod is investigated by choosing a larger damping coefficient. Then, the effect factors of dynamic relaxation method in post-buckling analysis of the compressed rods in cylinders are analyzed, with the compressive loads being far greater than that in Euler’s formula. To reduce the dynamic impact effect, a larger 𝛼 damping should be chosen. Under the condition of applying the initial lateral perturbation load, the post-buckling mode with helical shape and stable configuration can be obtained directly by avoiding the sinusoidal buckling process. In order to obtain the invariant deflection curve in the helical configuration, the number of beam elements should be at least 100. The effects of initial contact stiffness and initial penetration of the gap element on the final contact stiffness, penetration and contact force are also investigated. When the initial contact stiffness coefficient is greater than 0.3, the constant contact force can be obtained, and the final penetration is less than 0.1E−3 m. We studied the characteristic post-buckling modes with single point, two points, three points and point-line-point contact, respectively. In the configuration of point-line-point contact, it is further divided into just point-line-point contact and helical buckling deformation measured by the first contact point and continuous contact section. Except the critical load of helical buckling, the critical values of other configurations are in good agreement with the results of literatures, which indicates that the critical values obtained by the dynamic relaxation method have high accuracy. We also analyzed the distribution of shear force, bending moment and contact force in various post-buckling modes. In the helical buckling mode, we have found exceptions which are significantly different from those in the literatures. The transition section consists of the first non-contact section and the second non-contact section, but no perturbed-helix section. The more interesting conclusion is that, although there is no perturbed-helix, the length of two noncontact sections and the helical section obtained are in good agreement with the results in the literatures. It is found that the critical load of helical buckling obtained in this paper is smaller than that of the literature. The bending moment and contact force obtained are consistent with the results in the literatures, except the shear force.
[7] Lubinski A. A study of the buckling of rotary drilling strings. Drilling and production practice. American Petroleum Institute; 1950. [8] Huang W, Gao D, Liu Y. A study of tubular string buckling in vertical wells. Int J Mech Sci 2016;118:231–53. [9] Zhang Q, Jiang B, Huang W, Cui W, Liu J. Effect of wellhead tension on buckling load of tubular strings in vertical wells. J Petrol Sci Eng 2018;164:351–61. [10] Yue Q, Liu J, Zhang L, Zhang Q. The posting-buckling analysis and evaluations of limit drilling length for coiled tubing in the sidetrack horizontal well. J Petrol Sci Eng 2018;164:559–70. [11] Mitchell RF. Buckling behavior of well tubing: the packer effect. SPE J 1982;22(5):616–24. [12] Sorenson KG, Cheathamjr JB. Post-buckling behavior of a circular rod constrained within a circular cylinder. J Appl Mech 1986;53(4):929–34. [13] Chen JS, Fang J. Deformation sequence of a constrained spatial buckled beam under edge thrust. Int J Non Linear Mech 2013;55(10):98–101. [14] Fang J, Li SY, Chen JS. On a compressed spatial elastica constrained inside a tube. Acta Mech 2013;224(11):2635–47. [15] Huang W, Gao D, Wei S, Chen P. Boundary condition: a key factor in tubular string buckling. SPE J 2015;20(6):1409–20. [16] van der Heijden GHM, Champneys AR, Thompson JMT. Spatially complex localisation in twisted elastic rods constrained to a cylinder. Int J Solids Struct 2002;39(7):1863–83. [17] Gulyayev VI, Gaidaichuk VV, Solovjov IL, Gorbunovich IV. The buckling of elongated rotating drill strings. J Petrol Sci Eng 2009;67(3):140–8. [18] Zdvizhkov A, Miska S, Mitchell R. Measurement and analysis of induced torsion in helically buckled tubing. SPE Drill Complet 2005;24(2):266–75. [19] Chen JS, Li HC. On an elastic rod inside a slender tube under end twisting moment. J Appl Mech 2011;78(4):041009. [20] Kwon YW. Analysis of helical buckling. SPE Drill Eng 1988;3(2):211–16. [21] Zhang YL. A study of helical buckling of the compressed portion of a drillstring during straight-hole drilling using a down-hole motor. Soc Pet Eng 1989:1–23 Document ID: SPE-19369-MS. Available from SPE, Richardson, Texas. [22] Liu F, Xu B, Gao D. Packer effect analysis of helical buckling of well tubing. J Tsinghua Univ 1999;8(39):105–8 (in Chinese). [23] Gao G, Miska S. Effects of friction on post-buckling behavior and axial load transfer in a horizontal well. SPE J 2010;15(4):1110–24. [24] Wang X, Yuan Z. Investigation of frictional effects on the nonlinear buckling behavior of a circular rod laterally constrained in a horizontal rigid cylinder. J Petrol Sci Eng 2012;90–91:70–8. [25] Gao G, Di Q, Miska S, Wang W. Stability analysis of pipe with connectors in horizontal wells. SPE J 2012;17:931–41. [26] Huang W, Gao D. Helical buckling of a thin rod with connectors constrained in a cylinder. Int J Mech Sci 2014;84:189–98. [27] Mitchell RF, Miska S. Helical buckling of pipe with connectors and torque. SPE Drill Complet 2006;21(2):108–15. [28] Gulyayev VI, Gaidaichuk VV, Andrusenko EN, Shlyun NV. Critical buckling of drill strings in curvilinear channels of directed bore-holes. J Pet Sci Eng 2015;129:168–77. [29] Lubinski A, Althouse WS, Logan JL. Helical buckling of tubing sealed in packers. J Pet Technol 1962;14(6):655–70. [30] Yue Q, Liu J, Dong R. The research of post-buckling about slender rod string in wellbore based on energy method and experiment. J Pet Sci Eng 2017;156:732– 739. [31] Hajianmaleki M, Daily JS. Critical-buckling-load assessment of drillstrings in different wellbores by use of the explicit finite-element method. SPE Drill Complet 2014;29(2):256–64. [32] Liang Z, Zhu ZL. Critical helical buckling load assessment of coiled tubing under axial force by use of the explicit finite-element method. J Pet Sci Eng 2018;169:51–7. [33] Liu JP, Zhong XY, Cheng ZB, Feng XQ, Ren GX. Buckling of a slender rod confined in a circular tube: theory, simulation, and experiment. Int J Mech Sci 2018;140:288–305. [34] Xiao J, Chen Y, Lu X, Xu B, Chen X, Xu J. Three dimensional buckling beam under cylindrical constraint. Int J Mech Sci 2019;150:348–55. [35] Wang C, Li Z, Li Y, Qi M, Wang L. Experimental study on the rotation of compression-buckling rod column in a liquid-filled cylinder. J Petrol Sci Eng 2018;164:459–66. [36] Mitchell RF. The pitch of helically buckled pipe. Soc Pet Eng 2005 SPE-146959-MS. [37] Timoshenko SP, Gere JM. Theory of elastic stability. McGraw-Hill; 1961. [38] Zienkiewicz OC, Taylor RL, Nithiarasu P, Zhu JZ. The finite element method. McGraw-hill; 1977. [39] Ekhande SG, Selvappalam M, Madugula MK. Stability functions for three-dimensional beam-columns. J Struct Eng 1989;115(2):467–79. [40] Waszczyszyn Z, Cichon C, Radwanska M. Stability of structures by finite element methods. Elsevier; 1994. [41] Zhang X, Liu J. Gap-element method of three-dimension contact problem for drillstring in horizontal well. Comput Struct Mech Appl 1992;9(4):405–10 (in Chinese). [42] Liu J, Ding H, Zhang X. Application of gap element to nonlinear mechanics analysis of drillstring. Zhejiang Univ-Sci 2002;3(4):440–4. [43] Otter JRH. Computations for prestressed concrete reactor pressure vessel using dynamic relaxation. Nucl Struct Eng 1965;1(1):61–75. [44] Alamatian J. Displacement-based methods for calculating the buckling load and tracing the post-buckling regions with dynamic relaxation method. Comput Struct 2013;114–115:84–97. [45] Rezaiee-Pajand M, Estiri H. Computing the structural buckling limit load by using dynamic relaxation method. Int J Non Linear Mech 2016;81:245–60.
Acknowledgments The authors gratefully acknowledge the financial support under Nanyang Technological University, Singapore’s Academic Research Fund (AcRF) Tier 1 Grant no. RG 168/16. This research is also supported by National Natural Science Foundation of China (Grant numbers: 11502051, 51607035, 51674088), China Postdoctoral Science Foundation (2018M641804), Heilongjiang Postdoctoral Research Foundation (LBH-Q18029), and Heilongjiang Youth Innovation Talents of Ordinary Undergraduate Colleges and Universities (UNPYSCT-2018046). References [1] He X, Kyllingstad A. Helical buckling and lock-up conditions for coiled tubing in curved wells. SPE Drill Complet 1955;10(1):10–15. [2] Miller JT, Su T, Dussan VEB, Pabon J, Wicks N, Bertoldi K, Reis PM. Buckling-induced lock-up of a slender rod injected into a horizontal cylinder. Int J Solids Struct 2015;72:153–64. [3] Hirose H, Yamane K, Marhefka G, Cavarocchi N. Right ventricular rupture and tamponade caused by malposition of the avalon cannula for venovenous extracorporeal membrane oxygenation. J Cardiothorac Surg 2012;7(1):36. [4] Bao X, Guo S, Xiao N, Li Y, Yang C, Shen R, Cui J, Jiang Y, Liu X, Liu K. Operation evaluation in-human of a novel remote-controlled vascular interventional robot. Biomed Microdevices 2018;20(2):34. [5] Li T. A mechanics model of microtubule buckling in living cells. J Biomech 2008;41(8):1722–9. [6] Afrin T, Kabir AM, Sada K, Kakugo A, Nitta T. Buckling of microtubules on elastic media via breakable bonds. Biochem Bioph Res Co 2016;480(1):132–8.
115