Instability load analysis of a telescopic boom for an all-terrain crane

The instability load for the telescopic boom of an all-terrain crane is investigated in this paper. Combined with structural characteristics of the telescopic boom, each boom section is divided into several substructures, and the fixed-body coordinate system of each substructure is established based on the co-rotational method. A 3D Euler–Bernoulli eccentric beam element of the telescopic boom is derived. On the premise of considering the discretization of gravity and wind load, internal degrees of freedom of the substructure are condensed to the boundary nodes, forming a geometrical nonlinear super element. According to the nesting mode of the telescopic boom, a constraint way is established. The unstressed original length of the guy rope is calculated with a given preload so as to establish the equilibrium equations of the boom system with the external force of the guy rope and the corresponding tangent stiffness matrix. Regarding the above work, a new method for calculating the structural equilibrium path and instability load of telescopic boom structure is presented by solving the governing equations in a differential form. Finally, the method is validated by examples with different features.

1 Introduction

An all-terrain crane is a type of mobile crane. Due to the benefit of its lifting capacity and convenient mobility, it plays an important role in many construction fields (Ja et al., 2019; Yao et al., 2015). The telescopic boom of the all-terrain crane is a typical slender structure composed of multiple boom sections nested in one other. Each boom section is a box-type structure. When the telescopic boom carries a heavy load, its overall deformation exhibits a strong geometrical nonlinear effect. However, the integral stability of the slender structure is poor, which is prone to causing structural instability and result in accidents. The instability of the telescopic boom is one of the important reasons for all-terrain crane accidents (Neitzel et al., 2001). The instability load of the telescopic boom is one of the core indicators that determines the lifting capacity of all-terrain cranes. To improve the structural stability of the telescopic boom, a Y-shaped bracket is used to change its load conditions (Yao et al., 2020). The preload is applied to the telescopic boom through the guy rope to make the telescopic boom produce an initial deformation before lifting and reduces the deformation of telescopic boom during lifting. The Y-shaped bracket, the guy rope, and the back pull plate form the super lift system of the all-terrain crane.

In the structural instability analysis, external loads are commonly considered to be control parameters (Wang et al., 2015). Generally, the equilibrium path of the structures is followed by incremental methods as the control parameters changed. Response points, where tangent stiffness matrices become singular, are useful in the practical engineering applications, since they are related to the structural instability. There are generally two methods to determine the critical point of structural instability, i.e., the direct method and indirect method (Shi, 1996). The direct method is to add constraint equations according to the characteristics of the structure at the critical point (Fujii and Okazawa, 1997; Jari et al., 2012; Ding et al., 2014; Adnan and Mazen, 2000) and then directly solve the critical point based on the Newton–Raphson method. Most direct methods can converge from the region near the critical point to the critical point, and the point in the iterative process is not necessarily on the solution path. Most indirect methods determine whether there is a critical point according to the sign change in the determinant value of the structural tangent stiffness matrix on the structural equilibrium path (Bergan et al., 1978; Shi and Crisfield, 1994). A considerable number of constraint strategies are used to accurately determine the critical load and limit point, including load control, state control, and different kinds of arc length methods (Hellweg and Crisfield, 1998; Lu et al., 2005; Athisakul and Chucheepsakul, 2007). Although these methods have good adaptability in the structural stability calculation, they are also faced with the problem that the step size cannot be adjusted to achieve the required accuracy in the calculation process (Crisfield, 1983).

A correct and effective finite element model must be established for describing the geometrical nonlinear effect of the telescopic boom. At present, three approaches are often used for the finite element analysis of nonlinear solid and structural mechanics, namely total Lagrangian (TL; Pai et al., 2000; Nanakorn and Vu, 2006), updated Lagrangian (UL; Yang et al., 2007; Iu and Bradford, 2010), and co-rotational (CR) formulations (Crisfield and Moita, 1996; Felippa and Haugen, 2005; Li, 2007). Specifically, the CR formulation is suitable for describing the geometric nonlinearity of slender structures whose displacements and rotations may be arbitrarily large, while the local deformations are small. The CR formulations for a stability analysis of beams and shells have been studied by many researchers (Kisu, 1997; Hsiao and Lin, 2000; Battini and Pacoste, 2002; Verlinden et al., 2018).

Since the telescopic boom with a super lift system contains a large number of components, resulting in large-scale degrees of freedom to be solved, the CR formulation is more suitable for a finite element modeling of such slender structures composed of multiple boom sections. The whole structure is divided into several substructures, and the mechanical information of the internal nodes of a single sub structure is condensed to the boundary nodes, including the stiffness matrix (distributed forces), load matrix, etc. This substructure, formed by condensation, is called the super element, which is regarded as an individual element in the whole structural model during the modeling or analysis procedure (Li and Zhao, 2006). The substructure is defined by the characteristic points on a single boom section, and the nodal degrees of freedom in the substructure are statically condensed to form a super element, which reduces the computational burden for solving the global system variables (Mäkinen, 2007; Ghosh and Roy, 2009; He et al., 2010; Rantalainen et al., 2013). The precondition of adopting the static condensation method is that the nodal displacements and rotational angles are small in substructure, and the static condensation must be implemented in local coordinate systems.

The guy rope of the all-terrain crane connects the Y-shaped bracket and the telescopic boom head, which changes the initial deformation of the telescopic boom through a preload. Therefore, the unstressed original length of the guy rope must be calculated with the initial configuration of the telescopic boom without lifting the load, and then nonlinear equilibrium equations of the complete system can be established. Most existing work on cables mainly focuses on calculations with original length and without strain parameters (Jayaraman and Knudson, 1981; Gosling and Korban, 2001; Lee et al., 2003; Ju and Choo, 2005; Wang et al., 2015). However, calculating the unstressed original length of the guy rope under a known preload is necessary for the analysis of all-terrain cranes with super lift system. This calculation method has been studied in previous work (see Xu et al., 2022).

The purpose of this paper, therefore, is to establish the equilibrium equations of telescopic boom of an all-terrain crane, making use of co-rotational formulations, where a static condensation technique for the substructures and the nonlinear external forces from the guy rope with preload are integrated. It needs to be clarified that the material nonlinearity is not considered in this paper. After establishing the nonlinear equilibrium equations with the lifting load containing the control parameter, the accurate tangent stiffness matrix of the equations is derived based on the derivative of the equilibrium equations for the structural displacements. Finally, the load displacement curves of nodes in the telescopic boom are obtained by solving the differential form of the equilibrium equations. Combined with the singularity detection of the tangent stiffness matrix and the judgment criterion, the instability load is obtained.

The key points of this paper are mainly reflected in two aspects. First, the CR formulation is used to calculate the geometric nonlinear effect of the telescopic boom. Compared with the traditional method of establishing the CR formulations on each element (Wempner, 1969; Belytschko and Hsieh, 1973), the static condensation method is used to convert the corresponding substructure into a super element with two nodes, taking into account the influence of the self-weight and the external wind load, which has further improved the traditional static condensation method with zero external load (Przemieniecki, 1963; Bahar and Bahar, 2018). At the same time, the proposed super element can greatly reduce the dimension of the structural equilibrium equations. Second, the traditional load increment method needs to set a fixed load increment in the process of calculating the nonlinear equilibrium equations to obtain the structural equilibrium path (Wang et al., 2017). If the fixed increment is set too small, it will increase the calculation amount and reduce the calculation efficiency (Cheng et al., 1980). If the increment is too large, then the convergence may fail. Moreover, because it is impossible to accurately judge the range of the critical load in advance, the load increment is likely to directly cross the extreme point, which will lead to the failure of searching for the critical load. In this paper, the derivation of lifting the load parameter from the structural equilibrium equations of telescopic boom is converted into differential equations. With the advantage of automatic step size adjustment in the conventional differential equation solver, the function of automatically adjusting the load step size according to the corresponding nonlinearity of the current load state of the system is realized. Under the premise of ensuring the convergence of each step of solution, the structural equilibrium path can be quickly tracked and the critical load can be searched.

2 Substructure of telescopic boom and static condensation procedure

This paper takes a certain type of all-terrain crane as an example, as shown in Fig. 1. The telescopic boom is composed of eight boom sections, which can be combined into a boom with a length of 100 m. The combination of different lengths of the telescopic boom can be realized through the telescopic mechanism. The telescopic boom is connected to the slewing table with a shaft, and its luffing angle can be changed by a luffing cylinder. The Y-shaped bracket is installed on the first boom section of the telescopic boom. In accordance with the structural characteristics, the length of each boom is much larger than its section size, and the influence of shear deformation can be ignored. Therefore, a 3D Euler–Bernoulli beam element can be used for the finite element modeling of the telescopic boom. The boom section modulus in tension, bending, and torsion can be illustrated by the user-defined actual parameters of boom section (Dou et al., 2013). Unfortunately, this occurs under the assumption that cross sections of the beam element are rigid. In terms of the telescopic boom, some special conditions in engineering need to be calculated by establishing plate and shell element separately for local verification.

Figure 1 All-terrain crane structure.

2.1 Division of substructure and local coordinate system

Each boom section is divided into several substructures by the hinge point of the luffing cylinder, Y-shaped bracket, and overlap points between the boom sections, as shown in Fig. 2.

Figure 2 Substructures of a telescopic boom.

Unfortunately, the static condensation procedure cannot be used directly for the calculation of the geometrical nonlinear analysis of slender structures. However, if a co-rotational formulation for a substructure can be given, in which the elastic displacements and rotational angles are small, then the static condensation technique can be implemented for the telescopic boom structures. A substructure of the boom section after deflection is shown in Fig. 3. Euler angles and Cardan angles are widely used as generalized coordinates to describe the large rotation of the substructure (Wen, 1987; Cekus and Pawel, 2021). Due to the different rotation order of the Euler angle and the Cardan angle, in many practical applications in engineering, the latter is less singular than the former in a numerical calculation. Therefore, this paper uses Cardan angles to establish the transformation matrix between local coordinate system and global coordinate system.

Crisfield and Moita (1996) presented a unified formulation of the co-rotational approach for 3D elements with both translational and rotational degrees of freedom (Battini and Pacoste, 2002; Felippa and Haugen, 2005). The local coordinate system of the substructure is established at one side node of the substructure, and it is described by the global rotational angles (Cardan angles). A single substructure is a generalized beam element composed of multiple initially divided elements, and the origin of its section coordinate system is the node of the generalized beam element. The origin of the section coordinate system is not consistent with the section centroid, as shown in Fig. 3.

Figure 3 A local coordinate system of substructure.

$\left[{\mathit{g}}_{\mathrm{1}},{\mathit{g}}_{\mathrm{2}},{\mathit{g}}_{\mathrm{3}}\right]$ are base vectors of the global coordinate system of the telescopic boom system. The position vector of the origin of the substructure in the global coordinate system is r₁, and the right endpoint position vector is r₂. Global rotational angles are θ₁, θ₂, which are Cardan angles and can be used to describe the spatial rotation of the substructure coordinate system (Qi, 2008). A transformation matrix between local and global coordinate system is R₁. Its column vector can be determined by parameters ${\mathit{\theta }}_{\mathrm{1}}={\left[{\mathit{\alpha }}_{\mathrm{1}},{\mathit{\beta }}_{\mathrm{1}},{\mathit{\gamma }}_{\mathrm{1}}\right]}^{\mathrm{T}}$ as follows:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& {\mathit{e}}_{\mathrm{1}}=c_{\mathrm{2}}{c}_{\mathrm{3}}{\mathit{g}}_{\mathrm{1}}+(s_{\mathrm{1}}{s}_{\mathrm{2}}{c}_{\mathrm{3}}+c_{\mathrm{1}}{s}_{\mathrm{3}}){\mathit{g}}_{\mathrm{2}}+(s_{\mathrm{1}}{s}_{\mathrm{3}}-c_{\mathrm{1}}{s}_{\mathrm{2}}{c}_{\mathrm{3}}){\mathit{g}}_{\mathrm{3}}\\ & {\mathit{e}}_{\mathrm{2}}=-c_{\mathrm{2}}{s}_{\mathrm{3}}{\mathit{g}}_{\mathrm{1}}+(c_{\mathrm{1}}{c}_{\mathrm{3}}-s_{\mathrm{1}}{s}_{\mathrm{2}}{s}_{\mathrm{3}}){\mathit{g}}_{\mathrm{2}}+(s_{\mathrm{1}}{c}_{\mathrm{3}}+c_{\mathrm{1}}{s}_{\mathrm{2}}{s}_{\mathrm{3}}){\mathit{g}}_{\mathrm{3}}\\ & {\mathit{e}}_{\mathrm{3}}=s_{\mathrm{2}}{\mathit{g}}_{\mathrm{1}}-s_{\mathrm{1}}{c}_{\mathrm{2}}{\mathit{g}}_{\mathrm{2}}+c_{\mathrm{1}}{c}_{\mathrm{2}}{\mathit{g}}_{\mathrm{3}},\end{array}\end{array}$$

where s₁=sin α₁, c₁=cos α₁, s₂=sin β₁; c₂=cos β₁, s₃=sin γ₁, and c₃=cos γ₁.

$\mathit{\{}{\mathit{e}}_{\mathrm{1}},{\mathit{e}}_{\mathrm{2}},{\mathit{e}}_{\mathrm{3}}\mathit{\}}$ is a local coordinate system of a substructure. The position vector and derivative of right endpoint are given by the following:

$$\begin{array}{}\text{(2)}& {\mathit{r}}_{\mathrm{2}}={\mathit{r}}_{\mathrm{1}}+{\mathbf{R}}_{\mathrm{1}}({\widehat{\mathit{r}}}_{\mathrm{2}}+{\widehat{\mathit{u}}}_{\mathrm{2}}),{\dot{\mathit{r}}}_{\mathrm{2}}={\dot{\mathit{r}}}_{\mathrm{1}}+{\mathit{\omega }}_{\mathrm{1}}\times \left({\mathit{r}}_{\mathrm{2}}-{\mathit{r}}_{\mathrm{1}}\right)+{\mathbf{R}}_{\mathrm{1}}{\dot{\widehat{\mathit{u}}}}_{\mathrm{2}},\end{array}$$

where ${\widehat{\mathit{r}}}_{\mathrm{2}}$ is initial position vector of the right node in local coordinate system, and ${\widehat{\mathit{u}}}_{\mathrm{2}}$ is a translational displacement vector of the right node in local coordinate system. The corresponding virtual velocity is as follows:

$$\begin{array}{}\text{(3)}& \mathit{\delta }{\dot{\mathit{r}}}_{\mathrm{2}}=\mathit{\delta }{\dot{\mathit{r}}}_{\mathrm{1}}+\mathit{\delta }{\mathit{\omega }}_{\mathrm{1}}\times \left({\mathit{r}}_{\mathrm{2}}-{\mathit{r}}_{\mathrm{1}}\right)+{\mathbf{R}}_{\mathrm{1}}\mathit{\delta }{\dot{\widehat{\mathit{u}}}}_{\mathrm{2}},\end{array}$$

where ω₁ is the global angular velocity of the left end section. Its conversion relationship with θ₁ is as follows:

$$\begin{array}{}\text{(4)}& {\mathit{\omega }}_{\mathrm{1}}=\left[\begin{array}{ccc}\mathrm{1}& \mathrm{0}& \sin {\mathit{\beta }}_{\mathrm{1}}\\ \mathrm{0}& \cos {\mathit{\alpha }}_{\mathrm{1}}& -\sin {\mathit{\alpha }}_{\mathrm{1}}\cos {\mathit{\beta }}_{\mathrm{1}}\\ \mathrm{0}& \sin {\mathit{\alpha }}_{\mathrm{1}}& \cos {\mathit{\alpha }}_{\mathrm{1}}\cos {\mathit{\beta }}_{\mathrm{1}}\end{array}\right]\left[\begin{array}{c}{\dot{\mathit{\alpha }}}_{\mathrm{1}}\\ {\dot{\mathit{\beta }}}_{\mathrm{1}}\\ {\dot{\mathit{\gamma }}}_{\mathrm{1}}\end{array}\right]\triangleq {\mathbf{T}}_{\mathrm{1}}{\dot{\mathit{\theta }}}_{\mathrm{1}}.\end{array}$$

The virtual velocity of the degree of freedom in the global coordinate system is decomposed into three parts, i.e., the translational virtual velocity, the rotational virtual velocity, and the virtual velocity of the local degree of freedom in a local coordinate system.

The global angular velocity of the right endpoint section is given according to Eq. (4), as follows:

$$\begin{array}{}\text{(5)}& {\mathit{\omega }}_{\mathrm{2}}=\left[\begin{array}{ccc}\mathrm{1}& \mathrm{0}& \sin {\mathit{\beta }}_{\mathrm{2}}\\ \mathrm{0}& \cos {\mathit{\alpha }}_{\mathrm{2}}& -\sin {\mathit{\alpha }}_{\mathrm{2}}\cos {\mathit{\beta }}_{\mathrm{2}}\\ \mathrm{0}& \sin {\mathit{\alpha }}_{\mathrm{2}}& \cos {\mathit{\alpha }}_{\mathrm{2}}\cos {\mathit{\beta }}_{\mathrm{2}}\end{array}\right]\left[\begin{array}{c}{\dot{\mathit{\alpha }}}_{\mathrm{2}}\\ {\dot{\mathit{\beta }}}_{\mathrm{2}}\\ {\dot{\mathit{\gamma }}}_{\mathrm{2}}\end{array}\right]\triangleq {\mathbf{T}}_{\mathrm{2}}{\dot{\mathit{\theta }}}_{\mathrm{2}}.\end{array}$$

According to the angular velocity superposition principle (Qi, 2008), ω₂ can be given with ${\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}$, which is the angular velocity in the local coordinate system, as follows:

$$\begin{array}{}\text{(6)}& {\mathit{\omega }}_{\mathrm{2}}={\mathit{\omega }}_{\mathrm{1}}+{\mathbf{R}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}.\end{array}$$

In a local coordinate system, the rotation of the section coordinate system is small (Betsch and Steinmann, 2003), combining Eq. (5), as follows:

$$\begin{array}{}\text{(7)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}={\mathbf{E}}_{\mathrm{3}\times \mathrm{3}}{\dot{\stackrel{\mathrm{‾}}{\mathit{\theta }}}}_{\mathrm{2}}={\dot{\stackrel{\mathrm{‾}}{\mathit{\theta }}}}_{\mathrm{2}},\end{array}$$

where ${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{2}}$ is the rotational angle of the right end section of the substructure in a local coordinate system.

The virtual velocity of Eq. (6) is as follows:

$$\begin{array}{}\text{(8)}& \mathit{\delta }{\mathit{\omega }}_{\mathrm{2}}=\mathit{\delta }{\mathit{\omega }}_{\mathrm{1}}+{\mathbf{R}}_{\mathrm{1}}\mathit{\delta }{\dot{\stackrel{\mathrm{‾}}{\mathit{\theta }}}}_{\mathrm{2}}.\end{array}$$

The transformation relationship between local and global degrees of freedom needs to be given for transforming the virtual power equations into the following algebraic equation:

$$\begin{array}{}\text{(9)}& {\widehat{\mathit{u}}}_{\mathrm{1}}=\mathbf{0},{\widehat{\mathit{u}}}_{\mathrm{2}}={\mathbf{R}}_{\mathrm{1}}^{\mathrm{T}}({\mathit{r}}_{\mathrm{2}}-{\mathit{r}}_{\mathrm{1}})-({\widehat{\mathit{r}}}_{\mathrm{2}}-{\widehat{\mathit{r}}}_{\mathrm{1}}),\end{array}$$

where ${\widehat{\mathit{u}}}_{\mathrm{1}}$, ${\widehat{\mathit{u}}}_{\mathrm{2}}$ are local translational displacement vectors in local coordinate system. Their derivatives are as follows:

$$\begin{array}{}\text{(10)}& {\dot{\widehat{\mathit{u}}}}_{\mathrm{1}}=\mathbf{0},{\dot{\widehat{\mathit{u}}}}_{\mathrm{2}}={\mathbf{R}}_{\mathrm{1}}^{\mathrm{T}}({\dot{\mathit{r}}}_{\mathrm{2}}-{\dot{\mathit{r}}}_{\mathrm{1}})+{\mathbf{R}}_{\mathrm{1}}^{\mathrm{T}}({\stackrel{\mathrm{̃}}{\mathit{r}}}_{\mathrm{2}}-{\stackrel{\mathrm{̃}}{\mathit{r}}}_{\mathrm{1}}){\mathit{\omega }}_{\mathrm{1}},\end{array}$$

where a vector with a symbol “∼” on the top is its skew symmetric matrix. For example, if a vector $\mathit{a}={\left[\begin{array}{ccc}{a}_{\mathrm{1}}& a_{\mathrm{2}}& a_{\mathrm{3}}\end{array}\right]}^{\mathrm{T}}$, then its skew symmetric matrix is expressed as follows:

$$\stackrel{\mathrm{̃}}{\mathit{a}}=\left[\begin{array}{ccc}\mathrm{0}& -a_{\mathrm{3}}& a_{\mathrm{2}}\\ a_{\mathrm{3}}& \mathrm{0}& -a_{\mathrm{1}}\\ -a_{\mathrm{2}}& a_{\mathrm{1}}& \mathrm{0}\end{array}\right].$$

${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{2}}$ can be obtained according to ${\stackrel{\mathrm{‾}}{\mathbf{R}}}_{\mathrm{2}}$, which is a rotational matrix of section coordinate system to the local coordinate system, as follows:

$$\begin{array}{}\text{(11)}& {\stackrel{\mathrm{‾}}{\mathbf{R}}}_{\mathrm{2}}={\mathbf{R}}_{\mathrm{1}}^{\mathrm{T}}{\mathbf{R}}_{\mathrm{2}}.\end{array}$$

According to Eq. (8), its derivative is as follows:

$$\begin{array}{}\text{(12)}& {\dot{\stackrel{\mathrm{‾}}{\mathit{\theta }}}}_{\mathrm{2}}={\mathbf{R}}_{\mathrm{1}}^{\mathrm{T}}({\mathit{\omega }}_{\mathrm{2}}-{\mathit{\omega }}_{\mathrm{1}}).\end{array}$$

The derivative of a substructure of nodal variables in a global coordinate system is as follows:

$$\begin{array}{}\text{(13)}& \dot{\mathit{q}}=\left[{\dot{\mathit{r}}}_{\mathrm{1}};{\dot{\mathit{\theta }}}_{\mathrm{1}};{\dot{\mathit{r}}}_{\mathrm{2}};{\dot{\mathit{\theta }}}_{\mathrm{2}}\right].\end{array}$$

Combining Eqs. (10) and (12) yields the following:

$$\begin{array}{}\text{(14)}& {\displaystyle }\begin{array}{rl}& {\dot{\mathit{r}}}_{\mathrm{1}}=\left[\begin{array}{cccc}{\mathbf{E}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}\end{array}\right]\dot{\mathit{q}}\triangleq {\mathbf{T}}_{r_{\mathrm{1}}}\dot{\mathit{q}},\\ & {\dot{\mathit{r}}}_{\mathrm{2}}=\left[\begin{array}{cccc}{\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{E}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}\end{array}\right]\dot{\mathit{q}}\triangleq {\mathbf{T}}_{r_{\mathrm{2}}}\dot{\mathit{q}}\end{array}\text{(15)}& {\displaystyle }\begin{array}{rl}& {\mathit{\omega }}_{\mathrm{1}}=\left[\begin{array}{cccc}{\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{T}}_{\mathrm{1}}& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}\end{array}\right]\dot{\mathit{q}}\triangleq {\mathbf{T}}_{{\mathit{\omega }}_{\mathrm{1}}}\dot{\mathit{q}},\\ & {\mathit{\omega }}_{\mathrm{2}}=\left[\begin{array}{cccc}{\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{T}}_{\mathrm{2}}\end{array}\right]\dot{\mathit{q}}\triangleq {\mathbf{T}}_{{\mathit{\omega }}_{\mathrm{2}}}\dot{\mathit{q}}\end{array}\text{(16)}& {\displaystyle }\begin{array}{rl}& {\dot{\widehat{\mathit{u}}}}_{\mathrm{1}}=\left[\begin{array}{cccc}{\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}\end{array}\right]\dot{\mathit{q}}\triangleq {\mathbf{T}}_{{\widehat{u}}_{\mathrm{1}}}\dot{\mathit{q}},\\ & {\dot{\stackrel{\mathrm{‾}}{\mathit{\theta }}}}_{\mathrm{1}}=\left[\begin{array}{cccc}{\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}\end{array}\right]\dot{\mathit{q}}\triangleq {\mathbf{T}}_{{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{1}}}\dot{\mathit{q}}\end{array}\text{(17)}& {\displaystyle }\begin{array}{rl}& {\dot{\widehat{\mathit{u}}}}_{\mathrm{2}}=\left({\mathbf{R}}_{\mathrm{1}}^{\mathrm{T}}({\mathbf{T}}_{r_{\mathrm{2}}}-{\mathbf{T}}_{r_{\mathrm{1}}})+{\mathbf{R}}_{\mathrm{1}}^{\mathrm{T}}({\stackrel{\mathrm{̃}}{\mathit{r}}}_{\mathrm{2}}-{\stackrel{\mathrm{̃}}{\mathit{r}}}_{\mathrm{1}}){\mathbf{T}}_{{\mathit{\omega }}_{\mathrm{1}}}\right)\dot{\mathit{q}}\triangleq {\mathbf{T}}_{{\widehat{u}}_{\mathrm{2}}}\dot{\mathit{q}},\\ & {\dot{\stackrel{\mathrm{‾}}{\mathit{\theta }}}}_{\mathrm{2}}={\mathbf{R}}_{\mathrm{1}}^{\mathrm{T}}({\mathbf{T}}_{{\mathit{\omega }}_{\mathrm{2}}}-{\mathbf{T}}_{{\mathit{\omega }}_{\mathrm{1}}})\dot{\mathit{q}}\triangleq {\mathbf{T}}_{{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{2}}}\dot{\mathit{q}}.\end{array}\end{array}$$

The transformation relationship between the variables in local coordinate system and global coordinate system is as follows:

$$\begin{array}{}\text{(18)}& \dot{\widehat{\mathbf{U}}}=\left[\begin{array}{c}{\mathbf{T}}_{{\widehat{u}}_{\mathrm{1}}}\\ {\mathbf{T}}_{{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{1}}}\\ {\mathbf{T}}_{{\widehat{u}}_{\mathrm{2}}}\\ {\mathbf{T}}_{{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{2}}}\end{array}\right]\left[\begin{array}{c}{\dot{\mathit{r}}}_{\mathrm{1}}\\ {\dot{\mathit{\theta }}}_{\mathrm{1}}\\ {\dot{\mathit{r}}}_{\mathrm{2}}\\ {\dot{\mathit{\theta }}}_{\mathrm{2}}\end{array}\right]\triangleq {\mathbf{T}}_{rg}\dot{\mathit{q}},\end{array}$$

where $\dot{\widehat{\mathbf{U}}}=\left[{\dot{\widehat{\mathit{u}}}}_{\mathrm{1}};{\dot{\stackrel{\mathrm{‾}}{\mathit{\theta }}}}_{\mathrm{1}};{\dot{\widehat{\mathit{u}}}}_{\mathrm{2}};{\dot{\stackrel{\mathrm{‾}}{\mathit{\theta }}}}_{\mathrm{2}}\right]$.

2.2 Eccentric beam element

Each boom section of the telescopic boom is a slender box structure, and its finite element model can be built by 3D Euler–Bernoulli beam element. The nodal parameters of Euler–Bernoulli beam elements are located at the centroid of the section, and its linear stiffness matrix can be written as follows:

$$\begin{array}{}\text{(19)}& {\stackrel{\mathrm{‾}}{\mathit{k}}}_e=\left[\begin{array}{cccccccccccc}{k}_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -k_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& k_{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& k_{\mathrm{7}}& \mathrm{0}& -k_{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& k_{\mathrm{7}}\\ \mathrm{0}& \mathrm{0}& k_{\mathrm{3}}& \mathrm{0}& -k_{\mathrm{8}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -k_{\mathrm{3}}& \mathrm{0}& -k_{\mathrm{8}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& k_{\mathrm{4}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -k_{\mathrm{4}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -k_{\mathrm{8}}& \mathrm{0}& k_{\mathrm{5}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& k_{\mathrm{8}}& \mathrm{0}& k_{\mathrm{5}}/\mathrm{2}& \mathrm{0}\\ \mathrm{0}& k_{\mathrm{3}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& k_{\mathrm{6}}& \mathrm{0}& -k_{\mathrm{7}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& k_{\mathrm{6}}/\mathrm{2}\\ -k_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& k_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -k_{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -k_{\mathrm{7}}& \mathrm{0}& k_{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -k_{\mathrm{7}}\\ \mathrm{0}& \mathrm{0}& -k_{\mathrm{3}}& \mathrm{0}& k_{\mathrm{8}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& k_{\mathrm{3}}& \mathrm{0}& k_{\mathrm{8}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& -k_{\mathrm{4}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& k_{\mathrm{4}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -k_{\mathrm{8}}& \mathrm{0}& k_{\mathrm{5}}/\mathrm{2}& \mathrm{0}& \mathrm{0}& \mathrm{0}& k_{\mathrm{8}}& \mathrm{0}& k_{\mathrm{5}}& \mathrm{0}\\ \mathrm{0}& k_{\mathrm{7}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& k_{\mathrm{6}}/\mathrm{2}& \mathrm{0}& -k_{\mathrm{7}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& k_{\mathrm{6}}\end{array}\right],\end{array}$$

where $k_{\mathrm{1}}=\frac{EA}{L_e}$; $k_{\mathrm{2}}=\frac{\mathrm{12}E{I}_z}{L_e^{\mathrm{3}}}$; $k_{\mathrm{3}}=\frac{\mathrm{12}E{I}_y}{L_e^{\mathrm{3}}}$; $k_{\mathrm{4}}=\frac{GJ}{L_e}$; $k_{\mathrm{5}}=\frac{\mathrm{4}E{I}_y}{L_e}$; $k_{\mathrm{6}}=\frac{\mathrm{4}E{I}_z}{L_e}$; $k_{\mathrm{7}}=\frac{\mathrm{6}E{I}_z}{L_e^{\mathrm{2}}}$; $k_{\mathrm{8}}=\frac{\mathrm{6}E{I}_y}{L_e^{\mathrm{3}}}$. The parameters E, A, I_y, I_z, G, J, and L_e are the corresponding Young's modulus, cross-sectional area, y-axis moment of inertia, z-axis moment of inertia, shear modulus, polar moment of inertia, and the length of a beam element.

The section of the telescopic boom section is symmetrical about the y axis and asymmetric about the z axis, so it belongs to a kind of eccentric beam element, as shown in Fig. 3. However, in engineering, the selection of the nodal parameters to reflect the section characteristics is more appropriate, which can avoid the non-coincidence of nodal parameters between nested boom sections. For a single beam element in a substructure, the nodal parameters of the nodes at both ends in the element coordinate system are $\left[{\widehat{\mathit{u}}}_{e_{\mathrm{1}}};{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{e_{\mathrm{1}}}\right]$, $\left[{\widehat{\mathit{u}}}_{e_{\mathrm{2}}};{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{e_{\mathrm{2}}}\right]$. Based on the assumption of a rigid section of the beam element, the rotational angles of its section are consistent, and the rotational angles in the element coordinate system are small. Without losing generality, the position vector of the section centroid in the section coordinate system is ${\left[\begin{array}{ccc}\mathrm{0}& y& z\end{array}\right]}^{\mathrm{T}}$. ${\stackrel{\mathrm{‾}}{\mathit{u}}}_{e_{\mathrm{1}}}$ and ${\stackrel{\mathrm{‾}}{\mathit{u}}}_{e_{\mathrm{2}}}$ are displacements at centroid in section coordinate system.

$$\begin{array}{}\text{(20)}& {\stackrel{\mathrm{‾}}{\mathit{u}}}_{e_{\mathrm{1}}}={\widehat{\mathit{u}}}_{e_{\mathrm{1}}}+{\mathbf{D}}_e{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{e_{\mathrm{1}}},{\stackrel{\mathrm{‾}}{\mathit{u}}}_{e_{\mathrm{2}}}={\widehat{\mathit{u}}}_{e_{\mathrm{2}}}+{\mathbf{D}}_e{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{e_{\mathrm{2}}},\end{array}$$

where

$${\mathbf{D}}_e=\left[\begin{array}{ccc}\mathrm{0}& z& -y\\ -z& \mathrm{0}& \mathrm{0}\\ y& \mathrm{0}& \mathrm{0}\end{array}\right].$$

The conversion between the nodal parameters at the centroid and the origin of the section coordinate system is as follows:

$$\begin{array}{}\text{(21)}& \left[\begin{array}{c}{\stackrel{\mathrm{‾}}{\mathbf{U}}}_{e_{\mathrm{1}}}\\ {\stackrel{\mathrm{‾}}{\mathbf{U}}}_{e_{\mathrm{2}}}\end{array}\right]=\left[\begin{array}{cccc}{\mathbf{E}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{D}}_e& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}\\ {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{E}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}\\ {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{E}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{D}}_e\\ {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{E}}_{\mathrm{3}\times \mathrm{3}}\end{array}\right]\triangleq \mathbf{D}\left[\begin{array}{c}{\widehat{\mathbf{U}}}_{e_{\mathrm{1}}}\\ {\widehat{\mathbf{U}}}_{e_{\mathrm{2}}}\end{array}\right],\end{array}$$

where

$$\begin{array}{rl}& {\stackrel{\mathrm{‾}}{\mathbf{U}}}_{e_{\mathrm{1}}}=\left[\begin{array}{c}{\stackrel{\mathrm{‾}}{\mathit{u}}}_{e_{\mathrm{1}}}\\ {\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{e_{\mathrm{1}}}\end{array}\right],{\stackrel{\mathrm{‾}}{\mathbf{U}}}_{e_{\mathrm{2}}}=\left[\begin{array}{c}{\stackrel{\mathrm{‾}}{\mathit{u}}}_{e_{\mathrm{2}}}\\ {\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{e_{\mathrm{2}}}\end{array}\right],\\ & {\widehat{\mathbf{U}}}_{e_{\mathrm{1}}}=\left[\begin{array}{c}{\widehat{\mathit{u}}}_{e_{\mathrm{1}}}\\ {\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{e_{\mathrm{1}}}\end{array}\right],{\widehat{\mathbf{U}}}_{e_{\mathrm{2}}}=\left[\begin{array}{c}{\widehat{\mathit{u}}}_{e_{\mathrm{2}}}\\ {\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{e_{\mathrm{2}}}\end{array}\right].\end{array}$$

For a beam element, the deformation virtual power of the centroid parameters of the beam element can be transformed into the parameters of the origin of the section coordinate system as follows:

$$\begin{array}{}\text{(22)}& \mathit{\delta }{\stackrel{\mathrm{‾}}{\mathit{w}}}_e=\mathit{\delta }{\stackrel{\mathrm{‾}}{\mathbf{U}}}_e^{\mathrm{T}}{\stackrel{\mathrm{‾}}{\mathit{k}}}_e{\stackrel{\mathrm{‾}}{\mathbf{U}}}_e=\mathit{\delta }{\dot{\widehat{\mathbf{U}}}}_e^{\mathrm{T}}{\widehat{\mathit{k}}}_e{\widehat{\mathbf{U}}}_e=\mathit{\delta }{\widehat{\mathit{w}}}_e,\end{array}$$

where ${\stackrel{\mathrm{‾}}{\mathbf{U}}}_e=\left[\begin{array}{c}{\stackrel{\mathrm{‾}}{\mathbf{U}}}_{e_{\mathrm{1}}}\\ {\stackrel{\mathrm{‾}}{\mathbf{U}}}_{e_{\mathrm{2}}}\end{array}\right]$, ${\widehat{\mathbf{U}}}_e=\left[\begin{array}{c}{\widehat{\mathbf{U}}}_{e_{\mathrm{1}}}\\ {\widehat{\mathbf{U}}}_{e_{\mathrm{2}}}\end{array}\right]$. The stiffness matrix of the eccentric beam element is ${\widehat{\mathit{k}}}_e={\mathbf{D}}^{\mathrm{T}}{\stackrel{\mathrm{‾}}{\mathit{k}}}_e\mathbf{D}$.

2.3 Influence of gravity and wind load of substructure

During the operation of an all-terrain crane, the luffing angle of the telescopic boom changes continuously, so that the influence direction of gravity on the telescopic boom also changes, which cannot meet the condition that the external force of the traditional static condensation method is zero in the internal degree of freedom. The gravity in each element of the substructure is a uniformly distributed force. When the static condensation method is used to condense the structural degrees of freedom, then the self-weight of the element needs to be dispersed to the nodes at both ends.

The position vector of any point in the element of a substructure and its derivative are as follows:

$$\begin{array}{}\text{(23)}& {\mathit{r}}_e={\mathit{r}}_{\mathrm{1}}+{\mathbf{R}}_{\mathrm{1}}({\stackrel{\mathrm{‾}}{\mathit{r}}}_e+{\stackrel{\mathrm{‾}}{\mathit{u}}}_e),{\dot{\mathit{r}}}_e={\dot{\mathit{r}}}_{\mathrm{1}}-({\stackrel{\mathrm{̃}}{\mathit{r}}}_e-{\stackrel{\mathrm{̃}}{\mathit{r}}}_{\mathrm{1}}){\mathit{\omega }}_{\mathrm{1}}+{\mathbf{R}}_{\mathrm{1}}{\dot{\stackrel{\mathrm{‾}}{\mathit{u}}}}_e,\end{array}$$

where ${\stackrel{\mathrm{‾}}{\mathit{r}}}_e$, ${\stackrel{\mathrm{‾}}{\mathit{u}}}_e$ are the position vector and displacements of any point in beam element in local coordinate system, R₁ is the transformation matrix of the local coordinate system, and ω₁ is angular velocity of the local coordinate system.

The rotational angles of any point in the element can be obtained by the interpolation of nodal parameters at both ends for Euler–Bernoulli beam as follows:

$$\begin{array}{}\text{(24)}& {\stackrel{\mathrm{‾}}{\mathit{u}}}_e={\mathbf{N}}_u{\stackrel{\mathrm{‾}}{\mathbf{U}}}_e,{\dot{\stackrel{\mathrm{‾}}{\mathit{u}}}}_e={\mathbf{N}}_u{\dot{\stackrel{\mathrm{‾}}{\mathbf{U}}}}_e,{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_e={\mathbf{N}}_{\mathit{\theta }}{\stackrel{\mathrm{‾}}{\mathbf{U}}}_e,{\dot{\stackrel{\mathrm{‾}}{\mathit{\theta }}}}_e={\mathbf{N}}_{\mathit{\theta }}{\dot{\stackrel{\mathrm{‾}}{\mathbf{U}}}}_e,\end{array}$$

where N_u, N_θ are displacements and the rotational angle interpolation shape function of a 3D Euler–Bernoulli beam element (Qi, 2008).

By integrating the length of the beam element, it is seen that the virtual power of gravity has the following substructure:

$$\begin{array}{}\text{(25)}& \mathit{\delta }{\mathit{w}}_g=\mathit{\delta }{\dot{\mathit{r}}}_{\mathrm{1}}^{\mathrm{T}}{\mathbf{F}}_g+\mathit{\delta }{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{T}}{\mathbf{T}}_g+\mathit{\delta }{\dot{\widehat{\mathbf{U}}}}_e^{\mathrm{T}}{\stackrel{\mathrm{‾}}{\mathit{f}}}_e.\end{array}$$

The resultant force and moment at the origin of the local coordinate system are as follows:

$$\begin{array}{}\text{(26)}& {\mathbf{F}}_g=m_e\mathit{g};{\mathbf{T}}_g={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{m}_e{\mathbf{R}}_{\mathrm{1}}\left({\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{\mathit{r}}}}_{e_{\mathrm{1}}}+{\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{\mathit{r}}}}_{e_{\mathrm{2}}}\right)\stackrel{\mathrm{‾}}{\mathit{g}},\end{array}$$

where $\stackrel{\mathrm{‾}}{\mathit{g}}$ is the component of the gravitational acceleration g in the global coordinate system, m_e is the mass of beam element with length L_e, and ${\stackrel{\mathrm{‾}}{\mathit{r}}}_{e_{\mathrm{1}}}$ and ${\stackrel{\mathrm{‾}}{\mathit{r}}}_{e_{\mathrm{2}}}$ are position vectors of nodes at both ends of beam element in local coordinate system.

Combing Eqs. (21) and (25) yields the following:

$$\begin{array}{}\text{(27)}& m_e\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\mathit{\delta }{\dot{\stackrel{\mathrm{‾}}{\mathit{u}}}}_e^{\mathrm{T}}\mathrm{d}\mathit{\eta }=m_e\mathit{\delta }{\dot{\widehat{\mathbf{U}}}}_e^{\mathrm{T}}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{D}^{\mathrm{T}}{\mathbf{N}}_u^{\mathrm{T}}\mathrm{d}\mathit{\eta }\triangleq \mathit{\delta }{\dot{\widehat{\mathbf{U}}}}_e^{\mathrm{T}}{\stackrel{\mathrm{‾}}{\mathit{f}}}_e.\end{array}$$

Generalized forces distributed at the nodes at both ends of the element become the following:

$$\begin{array}{}\text{(28)}& {\stackrel{\mathrm{‾}}{\mathit{f}}}_e={\mathbf{G}}_g\stackrel{\mathrm{‾}}{\mathit{g}}.\end{array}$$

The gravity influence coefficient matrix is as follows:

$$\begin{array}{}\text{(29)}& {\mathbf{G}}_g=m_e{\left[\begin{array}{cccccccccccc}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{2}}z& -\frac{\mathrm{1}}{\mathrm{2}}y& \frac{\mathrm{1}}{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{2}}z& -\frac{\mathrm{1}}{\mathrm{2}}y\\ \mathrm{0}& \frac{\mathrm{1}}{\mathrm{2}}& \mathrm{0}& -\frac{\mathrm{1}}{\mathrm{2}}z& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{12}}{L}_e& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{2}}& \mathrm{0}& -\frac{\mathrm{1}}{\mathrm{2}}z& \mathrm{0}& -\frac{\mathrm{1}}{\mathrm{12}}{L}_e\\ \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{2}}& \frac{\mathrm{1}}{\mathrm{2}}y& -\frac{\mathrm{1}}{\mathrm{12}}{L}_e& \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{2}}& \frac{\mathrm{1}}{\mathrm{2}}y& \frac{\mathrm{1}}{\mathrm{12}}{L}_e& \mathrm{0}\end{array}\right]}^{\mathrm{T}}.\end{array}$$

Base on the actual working conditions, the external force on the telescopic boom during operation needs to consider whether the influence of wind load, the factors of the boom section shape, and boom expansion nesting have been considered. A wind load diagram is shown in Fig. 3.

In each substructure, the wind load of each element is equivalently calculated to the origin of the local coordinate system, and the linear density of its resultant force and resultant moment is ${\stackrel{\mathrm{‾}}{\mathit{f}}}_{\mathrm{0}}$ and ${\stackrel{\mathrm{‾}}{\mathit{m}}}_{\mathrm{0}}$. The virtual power of the wind load in the element is as follows:

$$\begin{array}{}\text{(30)}& \mathit{\delta }{w}_w=\mathit{\delta }{\dot{\mathit{r}}}_{\mathrm{1}}^{\mathrm{T}}{\mathbf{F}}_w+\mathit{\delta }{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{T}}{\mathbf{T}}_w+\mathit{\delta }{\dot{\widehat{\mathbf{U}}}}_e^{\mathrm{T}}{\stackrel{\mathrm{‾}}{\mathit{f}}}_w,\end{array}$$

where the resultant force and moment at the origin of the local coordinate system are

$${\mathbf{F}}_w=L_e{\mathbf{R}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathit{f}}}_{\mathrm{0}};{\mathbf{T}}_w=L_e{\mathbf{R}}_{\mathrm{1}}\left(\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{\mathit{r}}}}_{e_{\mathrm{1}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{\mathit{r}}}}_{e_{\mathrm{2}}}\right){\stackrel{\mathrm{‾}}{\mathit{f}}}_{\mathrm{0}}+{\stackrel{\mathrm{‾}}{\mathit{m}}}_{\mathrm{0}}\right).$$

The physical meaning of the generalized force f_w distributed at the nodes is as follows:

$$\begin{array}{}\text{(31)}& \begin{array}{rl}& \mathit{\delta }{\dot{\widehat{\mathbf{U}}}}_e^{\mathrm{T}}{\stackrel{\mathrm{‾}}{\mathit{f}}}_w=L_e\left(\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\mathit{\delta }{\dot{\widehat{\mathit{u}}}}_e^{\mathrm{T}}\mathrm{d}\mathit{\eta }\right){\stackrel{\mathrm{‾}}{\mathit{f}}}_{\mathrm{0}}\\ & +L_e\left(\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\mathit{\delta }{\dot{\stackrel{\mathrm{‾}}{\mathit{\theta }}}}_e^{\mathrm{T}}\mathrm{d}s\right){\stackrel{\mathrm{‾}}{\mathit{m}}}_{\mathrm{0}}\triangleq \mathit{\delta }{\dot{\widehat{\mathbf{U}}}}_e^{\mathrm{T}}\left({\mathbf{G}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathit{f}}}_{\mathrm{0}}\right)\\ & +\mathit{\delta }{\dot{\widehat{\mathbf{U}}}}_e^{\mathrm{T}}\left({\mathbf{G}}_{\mathrm{2}}{\stackrel{\mathrm{‾}}{\mathit{m}}}_{\mathrm{0}}\right),\end{array}\end{array}$$

where ${\widehat{\mathit{u}}}_e=\left[{\widehat{\mathit{u}}}_{e_{\mathrm{1}}};{\widehat{\mathit{u}}}_{e_{\mathrm{2}}}\right]$, ${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_e=\left[{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{e_{\mathrm{1}}};{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{e_{\mathrm{2}}}\right]$.

The wind load influence coefficient matrix is obtained according to Eqs. (20), (21), (24), and (31) as follows:

$$\begin{array}{}\text{(32)}& {\displaystyle }{\mathbf{G}}_{\mathrm{1}}={\left[\begin{array}{cccccccccccc}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{L}_e& -y& -z& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{2}}z{L}_e& -\frac{\mathrm{1}}{\mathrm{2}}y{L}_e& \frac{\mathrm{1}}{\mathrm{2}}{L}_e& y& z& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{2}}z{L}_e& -\frac{\mathrm{1}}{\mathrm{2}}y{L}_e\\ \mathrm{0}& \frac{\mathrm{1}}{\mathrm{2}}{L}_e& \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{12}}{L}_e^{\mathrm{2}}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{2}}{L}_e& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\frac{\mathrm{1}}{\mathrm{12}}{L}_e^{\mathrm{2}}\\ \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{2}}{L}_e& \mathrm{0}& -\frac{\mathrm{1}}{\mathrm{12}}{L}_e^{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{2}}{L}_e& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{12}}{L}_e^{\mathrm{2}}& \mathrm{0}\end{array}\right]}^{\mathrm{T}}\text{(33)}& {\displaystyle }{\mathbf{G}}_{\mathrm{2}}={\left[\begin{array}{cccccccccccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{2}}{L}_e& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{2}}{L}_e& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& y& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& -y& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -\mathrm{1}& \mathrm{0}& z& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& -z& \mathrm{0}& \mathrm{0}\end{array}\right]}^{\mathrm{T}}.\end{array}$$

2.4 Static condensation of substructure and super element

The local coordinate system is established for a substructure in which the nodal displacements are small. The nodes in each substructure can be divided into two groups, namely internal nodes n_i and boundary nodes n_b. The gravity and wind load of all elements in the substructure have been converted to the element nodes; therefore, the internal node degrees of freedom of the substructure can be condensed to the boundary node degrees of freedom to form a super element. In the local coordinate system, each beam element stiffness matrix ${\widehat{\mathit{k}}}_e$ is assembled to form the global stiffness matrix of substructure and is divided into blocks according to the degrees of freedom of boundary nodes and internal nodes.

$$\begin{array}{}\text{(34)}& \stackrel{\mathrm{‾}}{\mathbf{K}}=\left[\begin{array}{cc}{\stackrel{\mathrm{‾}}{\mathbf{K}}}_{bb}& {\stackrel{\mathrm{‾}}{\mathbf{K}}}_{bi}\\ {\stackrel{\mathrm{‾}}{\mathbf{K}}}_{ib}& {\stackrel{\mathrm{‾}}{\mathbf{K}}}_{ii}\end{array}\right].\end{array}$$

The deformation virtual power of a substructure can be expressed in a local system as follows:

$$\begin{array}{}\text{(35)}& \mathit{\delta }{w}_e=\mathit{\delta }{\dot{\widehat{\mathbf{U}}}}_b^{\mathrm{T}}\left({\stackrel{\mathrm{‾}}{\mathbf{K}}}_{bb}{\widehat{\mathbf{U}}}_b+{\stackrel{\mathrm{‾}}{\mathbf{K}}}_{bi}{\widehat{\mathbf{U}}}_i\right)+\mathit{\delta }{\dot{\widehat{\mathbf{U}}}}_i^{\mathrm{T}}\left({\stackrel{\mathrm{‾}}{\mathbf{K}}}_{ib}{\widehat{\mathbf{U}}}_b+{\stackrel{\mathrm{‾}}{\mathbf{K}}}_{ii}{\widehat{\mathbf{U}}}_i\right),\end{array}$$

where ${\widehat{\mathbf{U}}}_b$ and ${\widehat{\mathbf{U}}}_i$ are the boundary and internal nodal degrees of freedom in a substructure.

The total virtual power of gravity and wind load in substructure can be written as follows:

$$\begin{array}{}\text{(36)}& \mathit{\delta }{w}_f=\mathit{\delta }{\dot{\mathit{r}}}_{\mathrm{1}}^{\mathrm{T}}{\mathbf{F}}_{\mathrm{0}}+\mathit{\delta }{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{T}}{\mathbf{T}}_{\mathrm{0}}+\mathit{\delta }{\dot{\widehat{\mathbf{U}}}}_b^{\mathrm{T}}{\mathbf{F}}_b+\mathit{\delta }{\dot{\widehat{\mathbf{U}}}}_i^{\mathrm{T}}\left({\mathbf{G}}_i\stackrel{\mathrm{‾}}{\mathit{g}}+{\stackrel{\mathrm{‾}}{\mathbf{F}}}_{wi}\right),\end{array}$$

where G_i is the gravity influence coefficient matrix corresponding to the internal node degrees of freedom in the substructure. ${\stackrel{\mathrm{‾}}{\mathbf{F}}}_{wi}$ is the wind load corresponding to the internal node degrees of freedom according to the wind load influence coefficient matrix. F₀ and T₀ are the resultant force and moment when the gravity and wind load are equivalent to the origin of the local coordinate system. F_b is the gravity and wind load at the boundary node of the substructure.

Since the substructure boundary conditions and external forces are independent of the internal degrees of freedom, $\mathit{\delta }{\dot{\widehat{\mathbf{U}}}}_i$ is independent, and the coefficients in the virtual power expression are not affected by substructure assembly and external forces. According to the principle of virtual power, combining Eqs. (35) and (36) can yield the following:

$$\begin{array}{}\text{(37)}& {\stackrel{\mathrm{‾}}{\mathbf{K}}}_{ib}{\widehat{\mathbf{U}}}_b+{\stackrel{\mathrm{‾}}{\mathbf{K}}}_{ii}{\widehat{\mathbf{U}}}_i={\mathbf{G}}_i\stackrel{\mathrm{‾}}{\mathit{g}}+{\stackrel{\mathrm{‾}}{\mathbf{F}}}_{wi}.\end{array}$$

Internal degrees of freedom can be written as follows:

$$\begin{array}{}\text{(38)}& {\widehat{\mathbf{U}}}_i={\stackrel{\mathrm{‾}}{\mathbf{K}}}_{ii}^{-\mathrm{1}}{\mathbf{G}}_i\stackrel{\mathrm{‾}}{\mathit{g}}+{\stackrel{\mathrm{‾}}{\mathbf{K}}}_{ii}^{-\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathbf{F}}}_{wi}-{\stackrel{\mathrm{‾}}{\mathbf{K}}}_{ii}^{-\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathbf{K}}}_{ib}{\widehat{\mathbf{U}}}_b.\end{array}$$

The internal nodal displacements in the local coordinate system are small, so the elements in the stiffness matrix are constant, and the virtual velocity of Eq. (38) is

$$\begin{array}{}\text{(39)}& \mathit{\delta }{\dot{\widehat{\mathbf{U}}}}_i={\mathbf{T}}_{ib}\mathit{\delta }{\dot{\widehat{\mathbf{U}}}}_b,\end{array}$$

where the incidence matrix is

$$\begin{array}{}\text{(40)}& {\mathbf{T}}_{ib}=-{\stackrel{\mathrm{‾}}{\mathbf{K}}}_{ii}^{-\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathbf{K}}}_{ib}.\end{array}$$

Substituting Eqs. (39) and (40) into Eqs. (35) and (36) yields the following:

$$\begin{array}{}\text{(41)}& {\displaystyle }\mathit{\delta }{w}_e=\mathit{\delta }{\dot{\widehat{\mathbf{U}}}}_b^{\mathrm{T}}\left({\stackrel{\mathrm{‾}}{\mathbf{K}}}_{bb}+{\mathbf{T}}_{ib}^{\mathrm{T}}{\stackrel{\mathrm{‾}}{\mathbf{K}}}_{ib}\right){\widehat{\mathbf{U}}}_b\text{(42)}& {\displaystyle }\mathit{\delta }{w}_f=\mathit{\delta }{\dot{\mathit{r}}}_{\mathrm{1}}^{\mathrm{T}}{\mathbf{F}}_{\mathrm{0}}+\mathit{\delta }{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{T}}{\mathbf{T}}_{\mathrm{0}}+\mathit{\delta }{\dot{\widehat{\mathbf{U}}}}_b^{\mathrm{T}}\left({\mathbf{F}}_b+{\mathbf{T}}_{ib}^{\mathrm{T}}\left({\mathbf{G}}_i\stackrel{\mathrm{‾}}{\mathit{g}}+{\stackrel{\mathrm{‾}}{\mathbf{F}}}_{wi}\right)\right).\end{array}$$

Considering Eqs. (14), (15), and (18), Eqs. (41) and (42) can be written as follows:

$$\begin{array}{}\text{(43)}& \mathit{\delta }{w}_e\triangleq \mathit{\delta }{\dot{\mathit{q}}}^{\mathrm{T}}{\mathbf{F}}_b^i,\mathit{\delta }{w}_f\triangleq \mathit{\delta }{\dot{\mathit{q}}}^{\mathrm{T}}{\mathbf{F}}_b^a.\end{array}$$

The generalized internal and generalized external force are, respectively, written as follows:

$$\begin{array}{}\text{(44)}& {\displaystyle }{\mathbf{F}}_b^i={\mathbf{T}}_{rg}^{\mathrm{T}}\left({\mathbf{K}}_{bb}+{\mathbf{T}}_{ib}^{\mathrm{T}}{\mathbf{K}}_{ib}\right){\widehat{\mathbf{U}}}_b{\mathbf{T}}_{rg}^{\mathrm{T}}{\mathbf{K}}_e{\widehat{\mathbf{U}}}_b\text{(45)}& {\displaystyle }{\mathbf{F}}_b^a={\mathbf{T}}_{r_{\mathrm{1}}}^{\mathrm{T}}{\mathbf{F}}_{\mathrm{0}}+{\mathbf{T}}_{{\mathit{\omega }}_{\mathrm{1}}}^{\mathrm{T}}{\mathbf{T}}_{\mathrm{0}}+{\mathbf{T}}_{rg}^{\mathrm{T}}\left({\mathbf{F}}_b+{\mathbf{T}}_{ib}^{\mathrm{T}}\left({\mathbf{G}}_i\stackrel{\mathrm{‾}}{\mathit{g}}+{\stackrel{\mathrm{‾}}{\mathbf{F}}}_{wi}\right)\right).\end{array}$$

The substructure composed of multiple beam elements is reduced to a super beam element expressed by degrees of freedom at both ends, and K_e is its equivalent stiffness matrix.

3 Constraints and boundary conditions of telescopic boom

The nested connection of the telescopic boom of an all-terrain crane is the form of a hydraulic cylinder and pin; that is, different combinations of multiple boom sections are realized successively through the telescopic hydraulic cylinder in the first boom section. Each boom section is connected through the boom pin at the tail, and a nylon slider is designed at the tail and head of the boom section for overlapping. Except for the first and last boom sections, four sections in each boom section generate relevant constraints between the outer and inner boom sections. The left overlap point, right overlap point, internal overlap point, and outer overlap point are shown in Fig. 4.

The divided substructure is a super element in each boom section. We take three boom sections as an example, as shown in Fig. 5.

Figure 4 Overlap point of boom sections.

Figure 5 Boom section division.

3.1 Constraints between boom sections

The first two boom sections are used to explain the constraints between them. The connection of the two boom sections corresponds to two types of constraints, where section 1 and section 3 correspond to the rotating joint in the multi-body theory and only have the rotational degrees of freedom about the axis of the boom pin, and section 2 and section 4 correspond to the prismatic joint, which constrains the displacements around the main axis of the section and the rotational degrees of freedom around the normal part of the section. The nodal position vectors and rotational angles in the global coordinate system for the left and right overlap points are ${\mathit{r}}_{b_i}$, α_i, β_i, and γ_i. The corresponding section coordinate system is $\left\{{\mathit{b}}_{\mathrm{1}}^i,{\mathit{b}}_{\mathrm{2}}^i,{\mathit{b}}_{\mathrm{3}}^i\right\}$, which is shown in Fig. 6.

Figure 6 Constraint relationships of the boom sections.

According to the structural form of the boom pin connection, the local part can be regarded as a rigid body. d_L is the horizontal distance between the origin of the two local coordinate systems in the initial configuration, as shown in Fig. 7. The constraint equations at the boom pin connection can be written as follows:

$$\begin{array}{}\text{(46)}& {\mathit{r}}_{b_{\mathrm{3}}}={\mathit{r}}_{b_{\mathrm{1}}}+d_L{\mathit{b}}_{\mathrm{3}}^{\mathrm{1}},{\mathit{b}}_{\mathrm{3}}^{\mathrm{3}}\cdot {\mathit{b}}_{\mathrm{1}}^{\mathrm{1}}=\mathrm{0},{\mathit{b}}_{\mathrm{3}}^{\mathrm{3}}\cdot {\mathit{b}}_{\mathrm{2}}^{\mathrm{1}}=\mathrm{0}.\end{array}$$

The constraint equations at the right end of the first boom section with sliders can be written as follows:

$$\begin{array}{}\text{(47)}& {\mathit{b}}_{\mathrm{3}}^{\mathrm{4}}\cdot {\mathit{b}}_{\mathrm{2}}^{\mathrm{2}}=\mathrm{0},({\mathit{r}}_{b_{\mathrm{4}}}-{\mathit{r}}_{b_{\mathrm{2}}})\cdot {\mathit{b}}_{\mathrm{2}}^{\mathrm{2}}=\mathrm{0},{\mathit{r}}_{b_{\mathrm{4}}}-{\mathit{r}}_{b_{\mathrm{2}}})\cdot {\mathit{b}}_{\mathrm{3}}^{\mathrm{2}}=\mathrm{0}.\end{array}$$

Figure 7 Overlap point of two boom sections.

3.2 Boundary constraints of telescopic boom

There is only an independent rotational degree of freedom around the boom foot pin axis, which can be considered to be the boom luffing angle. The local structure can be processed as a rigid part, according to the actual structure at the boom foot, as shown in Fig. 8a.

Figure 8 Foot pin and luffing hydraulic cylinder of a telescopic boom.

The rigid part has the same angular velocity as the local coordinate system of the first super element. Therefore, the luffing angle degree of freedom can be expressed by the Cardan angle β₁ in the global coordinate system. β₁ is used as the independent degree of freedom at the constraint of boom foot pin, and the corresponding rotational matrix can be written as follows:

$$\begin{array}{}\text{(48)}& {\mathbf{R}}_{p_{\mathrm{1}}}=\left[\begin{array}{ccc}\cos {\mathit{\beta }}_{\mathrm{1}}& \mathrm{0}& -\sin {\mathit{\beta }}_{\mathrm{1}}\\ & \mathrm{1}& \\ \sin {\mathit{\beta }}_{\mathrm{1}}& \mathrm{0}& \cos {\mathit{\beta }}_{\mathrm{1}}\end{array}\right].\end{array}$$

The position vector of the left end node in the global coordinate system is ${\mathit{r}}_{p_{\mathrm{1}}}$, and the rotational angles are ${\left[\begin{array}{ccc}{\mathit{\alpha }}_{\mathrm{1}}& {\mathit{\beta }}_{\mathrm{1}}& {\mathit{\gamma }}_{\mathrm{1}}\end{array}\right]}^{\mathrm{T}}$. The constraint equations of the boom foot pin are

$$\begin{array}{}\text{(49)}& {\mathit{r}}_{p_{\mathrm{1}}}={\mathbf{R}}_{p_{\mathrm{1}}}\mathrm{\Delta }{\stackrel{\mathrm{‾}}{\mathit{r}}}_{\mathrm{1}},{\mathit{\alpha }}_{\mathrm{1}}={\mathit{\gamma }}_{\mathrm{1}}=\mathrm{0}.\end{array}$$

The luffing of the telescopic boom is achieved by the luffing hydraulic cylinder, whose bottom and top hinge points are, respectively, connected with the slewing table and the first boom section, as shown in Fig. 8b. The luffing hydraulic cylinder is regarded as a constant length constraint in the finite element modeling and calculation. r_c and θ_c are the nodal variables of top hinge point, and the position vector of top hinge point is as follows:

$$\begin{array}{}\text{(50)}& {\mathit{r}}_t={\mathit{r}}_c+{\mathbf{R}}_c\mathrm{\Delta }{\stackrel{\mathrm{‾}}{\mathit{r}}}_{\mathrm{2}},\end{array}$$

where R_c is the transformation matrix at top hinge point section, which is formed by the base vectors of the local coordinate system.

The length of the luffing hydraulic cylinder under the corresponding luffing angle is a fixed value, r_h is the position vector of the luffing hydraulic cylinder in global coordinate system, and the constraint equation can be written as follows:

$$\begin{array}{}\text{(51)}& ∥{\mathit{r}}_t-{\mathit{r}}_b∥=L_h.\end{array}$$

4 External force of the telescopic boom

The external nodal forces in the telescopic boom system equations are applied to the nodes of the super element, including the self-weight, external wind load, guy rope force, back pull plate force, and lifting load, as shown in Fig. 9. The first two kind of loads have been considered in Sect. 3. The last two kind of loads can be referred to in a previous work (Xu et al., 2022). A method for calculating the unstressed original length of a guy rope with known preload is derived. After the unstressed original length of the guy rope is obtained by applying the preload under the unloaded state of the telescopic boom, the force of guy rope acting on boom system can be calculated during the lifting process.

Figure 9 Guy rope, back plate, and lifting load diagram.

4.1 Guy rope force

Base on the structural characteristics of all-terrain crane, the Y-shaped bracket remains perpendicular to the axial direction of the telescopic boom after installation and has good stiffness. It can be considered to be a rigid part in the calculation of the telescopic boom. The Y-shaped bracket and the guy rope are symmetrically connected on both sides.

As shown in Fig. 9a, ${\mathit{r}}_{s_{\mathrm{1}}}$ and ${\mathit{\theta }}_{s_{\mathrm{1}}}$ are the nodal position vectors and rotational angles of the Y-shaped bracket on the section of the telescopic boom installation point. ${\mathit{r}}_{s_{\mathrm{2}}}$ and ${\mathit{\theta }}_{s_{\mathrm{2}}}$ are the variables of the node, which is the boom head connecting point. e_b and ${\mathit{e}}_{wi(i=\mathrm{1},\mathrm{2})}$ are the unit vectors of the telescopic boom axis and Y-shaped bracket axis, respectively. ${\stackrel{\mathrm{‾}}{\mathit{r}}}_{bi(i=\mathrm{1},\mathrm{2})}$ are the position vectors of the guide pulley in local coordinate system. ${\mathit{\varphi }}_{i(i=\mathrm{1},\mathrm{2})}$ are the azimuth of the guy rope. A single side guy rope connection is selected for the description.

The derivative of the guide pulley's center in the global coordinate system is as follows:

$$\begin{array}{}\text{(52)}& \begin{array}{rl}{\dot{\mathit{r}}}_{b\mathrm{1}}& ={\dot{\mathit{r}}}_{s_{\mathrm{1}}}+{\mathit{\omega }}_{s_{\mathrm{1}}}\times {\mathbf{R}}_{s_{\mathrm{1}}}{\stackrel{\mathrm{‾}}{\mathit{r}}}_{b\mathrm{1}}\\ & ={\dot{\mathit{r}}}_{s_{\mathrm{1}}}+\left({\stackrel{\mathrm{̃}}{\mathit{r}}}_{s_{\mathrm{1}}}-{\stackrel{\mathrm{̃}}{\mathit{r}}}_{b\mathrm{1}}\right){\mathbf{T}}_{{\mathit{\omega }}_{s_{\mathrm{1}}}}{\dot{\mathit{\theta }}}_{s_{\mathrm{1}}}{\mathbf{T}}_{r_{b\mathrm{1}}}{\dot{\mathit{q}}}_s.\end{array}\end{array}$$

Here, ${\mathbf{R}}_{s_{\mathrm{1}}}$ is the transformation matrix at the Y-shaped bracket installation point section, which is formed by the base vector of the local coordinate system. ${\mathbf{T}}_{{\mathit{\omega }}_{s_{\mathrm{1}}}}$ is the coefficient matrix of the angular velocity of the local coordinate system.

$$\begin{array}{}\text{(53)}& {\displaystyle }{\dot{\mathit{q}}}_s=\left[{\dot{\mathit{r}}}_{s_{\mathrm{1}}};{\dot{\mathit{\theta }}}_{s_{\mathrm{1}}};{\dot{\mathit{r}}}_{s_{\mathrm{2}}};{\dot{\mathit{\theta }}}_{s_{\mathrm{2}}}\right]\text{(54)}& {\displaystyle }{\mathbf{T}}_{r_{b\mathrm{1}}}=\left[{\mathbf{E}}_{\mathrm{3}\times \mathrm{3}},\left({\stackrel{\mathrm{̃}}{\mathit{r}}}_{s_{\mathrm{1}}}-{\stackrel{\mathrm{̃}}{\mathit{r}}}_{b\mathrm{1}}\right){\mathbf{T}}_{{\mathit{\omega }}_{s_{\mathrm{1}}}},{\mathbf{0}}_{\mathrm{3}\times \mathrm{3}},{\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}\right].\end{array}$$

The position vector of the guy rope entry point at the guide pulley and its derivative are as follows:

$$\begin{array}{}\text{(55)}& \begin{array}{rl}& {\mathit{r}}_{p\mathrm{1}}={\mathit{r}}_{b\mathrm{1}}+\mathrm{\Delta }r{\mathit{n}}_{\mathrm{1}},\\ & {\dot{\mathit{r}}}_{p\mathrm{1}}=\left({\mathbf{T}}_{r_{b\mathrm{1}}}+\left[{\mathbf{0}}_{\mathrm{3}\times \mathrm{3}},\left({\stackrel{\mathrm{̃}}{\mathit{r}}}_{b\mathrm{1}}-{\stackrel{\mathrm{̃}}{\mathit{r}}}_{p\mathrm{1}}\right){\mathbf{T}}_{{\mathit{\omega }}_{s_{\mathrm{1}}}},{\mathbf{0}}_{\mathrm{3}\times \mathrm{3}},{\mathbf{0}}_{\mathrm{3}\times \mathrm{3}}\right]\right){\dot{\mathit{q}}}_s\triangleq {\mathbf{T}}_{r_{p\mathrm{1}}}{\dot{\mathit{q}}}_s\end{array}.\end{array}$$

Here, Δr is the radius of the guide pulley, and ${\mathit{n}}_{\mathrm{1}}=\mathrm{\Delta }r\cos {\mathit{\varphi }}_{\mathrm{1}}{\mathit{e}}_b+\mathrm{\Delta }r\sin {\mathit{\varphi }}_{\mathrm{1}}{\mathit{e}}_{\mathrm{w}\mathrm{1}}$.

The position vector of the guy rope connecting point at boom head and its derivative are as follows:

$$\begin{array}{}\text{(56)}& \begin{array}{rl}& {\mathit{r}}_{q\mathrm{1}}={\mathit{r}}_{s_{\mathrm{2}}}+{\mathbf{R}}_{s_{\mathrm{2}}}{\stackrel{\mathrm{‾}}{\mathit{r}}}_{q\mathrm{1}},\\ & {\dot{\mathit{r}}}_{q\mathrm{1}}=\left[{\mathbf{0}}_{\mathrm{3}\times \mathrm{3}},{\mathbf{0}}_{\mathrm{3}\times \mathrm{3}},{\mathbf{E}}_{\mathrm{3}\times \mathrm{3}},\left({\stackrel{\mathrm{̃}}{\mathit{r}}}_{s_{\mathrm{2}}}-{\stackrel{\mathrm{̃}}{\mathit{r}}}_{q\mathrm{1}}\right){\mathbf{T}}_{{\mathit{\omega }}_{s_{\mathrm{2}}}}\right]{\dot{\mathit{q}}}_s\triangleq {\mathbf{T}}_{r_{q\mathrm{1}}}{\dot{\mathit{q}}}_s.\end{array}\end{array}$$

The loads acting on the guide pulley and the connecting point of the guy rope are F_p1, F_p2, F_q1, and F_q2, respectively, which can be calculated on the premise that the original length of the guy rope is known under the corresponding working conditions. The virtual power of the guy rope force to the telescopic boom can be expressed as follows:

$$\begin{array}{}\text{(57)}& \mathit{\delta }{w}_q=\mathit{\delta }{\dot{\mathit{q}}}_s^{\mathrm{T}}\left({\mathbf{T}}_{r_{p\mathrm{1}}}^{\mathrm{T}}{\mathbf{F}}_{p\mathrm{1}}+{\mathbf{T}}_{r_{p\mathrm{2}}}^{\mathrm{T}}{\mathbf{F}}_{p\mathrm{2}}+{\mathbf{T}}_{r_{q\mathrm{1}}}^{\mathrm{T}}{\mathbf{F}}_{q\mathrm{1}}+{\mathbf{T}}_{r_{q\mathrm{2}}}^{\mathrm{T}}{\mathbf{F}}_{q\mathrm{2}}\right).\end{array}$$

4.2 Back pull plate force

The original length of the back pull plate is known, as shown in Fig. 9b. The forces of the plate on the connecting point of the Y-shaped bracket and the first boom section can be calculated with the method in the literature (Xu et al., 2022). Forces on both sides are F_j1, F_j2, F_h1, and F_h2, respectively. ${\mathit{r}}_{s_{\mathrm{3}}}$ and ${\mathit{\theta }}_{s_{\mathrm{3}}}={\left[\begin{array}{ccc}{a}_{\mathrm{3}}& {\mathit{\beta }}_{\mathrm{3}}& {\mathit{\gamma }}_{\mathrm{3}}\end{array}\right]}^{\mathrm{T}}$ are the nodal position vector and rotational angles of the back pull plate at the connection point. The hinge point at the foot of the boom is fixed axis rotation constraint, and the independent degree of freedom is β₃. ${\stackrel{\mathrm{‾}}{\mathit{r}}}_{ji}$ and ${\stackrel{\mathrm{‾}}{\mathit{r}}}_{hi}$ are the position vectors of the back pull plate connection point in a local coordinate system.

The derivative of the connection point between back pull plate in a global coordinate system is as follows:

$$\begin{array}{}\text{(58)}& \begin{array}{rl}{\dot{\mathit{r}}}_{j\mathrm{1}}& ={\dot{\mathit{r}}}_{s_{\mathrm{1}}}+{\mathit{\omega }}_{s_{\mathrm{1}}}\times {\mathbf{R}}_{s_{\mathrm{1}}}{\stackrel{\mathrm{‾}}{\mathit{r}}}_{j\mathrm{1}}\\ & ={\dot{\mathit{r}}}_{s_{\mathrm{1}}}+\left({\stackrel{\mathrm{̃}}{\mathit{r}}}_{s_{\mathrm{1}}}-{\stackrel{\mathrm{̃}}{\mathit{r}}}_{j\mathrm{1}}\right){\mathbf{T}}_{{\mathit{\omega }}_{s\mathrm{1}}}{\dot{\mathit{\theta }}}_{s_{\mathrm{1}}}\\ & =\left[{\mathbf{E}}_{\mathrm{3}\times \mathrm{3}},\left({\stackrel{\mathrm{̃}}{\mathit{r}}}_{s_{\mathrm{1}}}-{\stackrel{\mathrm{̃}}{\mathit{r}}}_{j\mathrm{1}}\right){\mathbf{T}}_{{\mathit{\omega }}_{s\mathrm{1}}}\right]\left[\begin{array}{c}{\dot{\mathit{r}}}_{s_{\mathrm{1}}}\\ {\dot{\mathit{\theta }}}_{s_{\mathrm{1}}}\end{array}\right]\triangleq {\mathbf{T}}_{r_{j\mathrm{1}}}\left[\begin{array}{c}{\dot{\mathit{r}}}_{s_{\mathrm{1}}}\\ {\dot{\mathit{\theta }}}_{s_{\mathrm{1}}}\end{array}\right].\end{array}\end{array}$$

The derivative of a connection point between the back pull plate and foot boom section in global coordinate system is as follows:

$$\begin{array}{}\text{(59)}& {\dot{\mathit{r}}}_{h\mathrm{1}}={\stackrel{\mathrm{̃}}{\mathit{g}}}_{\mathrm{2}}{\mathit{r}}_{h\mathrm{1}}{\dot{\mathit{\beta }}}_{\mathrm{3}}\triangleq {\mathbf{T}}_{r_{h\mathrm{1}}}{\dot{\mathit{\beta }}}_{\mathrm{3}}.\end{array}$$

The virtual power at the Y-shaped bracket connection point is as follows:

$$\begin{array}{}\text{(60)}& \begin{array}{rl}\mathit{\delta }{w}_j& ={\dot{\mathit{r}}}_{j\mathrm{1}}^{\mathrm{T}}{\mathbf{F}}_{j\mathrm{1}}+{\dot{\mathit{r}}}_{j\mathrm{2}}^{\mathrm{T}}{\mathbf{F}}_{j\mathrm{2}}\\ & =\left[\begin{array}{cc}\mathit{\delta }{\dot{\mathit{r}}}_{s_{\mathrm{1}}}^{\mathrm{T}}& \mathit{\delta }{\dot{\mathit{\theta }}}_{s_{\mathrm{1}}}^{\mathrm{T}}\end{array}\right]\left({\mathbf{T}}_{r_{j\mathrm{1}}}^{\mathrm{T}}{\mathbf{F}}_{j\mathrm{1}}+{\mathbf{T}}_{r_{j\mathrm{2}}}^{\mathrm{T}}{\mathbf{F}}_{j\mathrm{2}}\right).\end{array}\end{array}$$

The virtual power at the foot boom section connection point is as follows:

$$\begin{array}{}\text{(61)}& \mathit{\delta }{w}_h={\dot{\mathit{r}}}_{h\mathrm{1}}^{\mathrm{T}}{\mathbf{F}}_{h\mathrm{1}}+{\dot{\mathit{r}}}_{h\mathrm{2}}^{\mathrm{T}}{\mathbf{F}}_{h\mathrm{2}}=\mathit{\delta }{\dot{\mathit{\beta }}}_{\mathrm{3}}\left({\mathbf{T}}_{r_{h\mathrm{1}}}^{\mathrm{T}}{\mathbf{F}}_{h\mathrm{1}}+{\mathbf{T}}_{r_{h\mathrm{2}}}^{\mathrm{T}}{\mathbf{F}}_{h\mathrm{2}}\right).\end{array}$$

4.3 Lifting load and single lifting rope load

As shown in Fig. 9c, ${\stackrel{\mathrm{‾}}{\mathit{r}}}_{qg}$ and ${\stackrel{\mathrm{‾}}{\mathit{r}}}_{qs}$ are the position vectors of the equivalent action point of lifting load and lifting rope in the local coordinate system of a boom head connection. F_qg is the lifting load force under a corresponding working condition. e_s is the unit vector of the lifting rope between the boom guide pulley and the lifting winch. n_q is the number of lifting ropes corresponding to the lifting load, and the single lifting rope force is as follows:

$$\begin{array}{}\text{(62)}& {\mathbf{F}}_{qs}={\displaystyle \frac{∥{\mathbf{F}}_{qg}∥}{n_q}}{\mathit{e}}_s.\end{array}$$

The derivative of a lifting load action point in the global coordinate system is as follows:

$$\begin{array}{}\text{(63)}& \begin{array}{rl}{\dot{\mathit{r}}}_{qg}& ={\dot{\mathit{r}}}_{s_{\mathrm{4}}}+{\mathit{\omega }}_{s_{\mathrm{4}}}\times {\mathbf{R}}_{s_{\mathrm{4}}}{\stackrel{\mathrm{‾}}{\mathit{r}}}_{qg}\\ & ={\dot{\mathit{r}}}_{s_{\mathrm{4}}}+\left({\stackrel{\mathrm{̃}}{\mathit{r}}}_{s_{\mathrm{4}}}-{\stackrel{\mathrm{̃}}{\mathit{r}}}_{qg}\right){\mathbf{T}}_{{\mathit{\omega }}_{s\mathrm{4}}}{\dot{\mathit{\theta }}}_{s_{\mathrm{4}}}\triangleq {\mathbf{T}}_{r_{qg}}\left[\begin{array}{c}{\dot{\mathit{r}}}_{s_{\mathrm{4}}}\\ {\dot{\mathit{\theta }}}_{s_{\mathrm{4}}}\end{array}\right].\end{array}\end{array}$$

The derivative of a lifting rope action point in the global coordinate system is as follows:

$$\begin{array}{}\text{(64)}& \begin{array}{rl}{\dot{\mathit{r}}}_{qs}& ={\dot{\mathit{r}}}_{s_{\mathrm{4}}}+{\mathit{\omega }}_{s_{\mathrm{4}}}\times {\mathbf{R}}_{s_{\mathrm{4}}}{\stackrel{\mathrm{‾}}{\mathit{r}}}_{qs}\\ & ={\dot{\mathit{r}}}_{s_{\mathrm{4}}}+\left({\stackrel{\mathrm{̃}}{\mathit{r}}}_{s_{\mathrm{4}}}-{\stackrel{\mathrm{̃}}{\mathit{r}}}_{qs}\right){\mathbf{T}}_{{\mathit{\omega }}_{s_{\mathrm{4}}}}{\dot{\mathit{\theta }}}_{s_{\mathrm{4}}}\triangleq {\mathbf{T}}_{r_{qs}}\left[\begin{array}{c}{\dot{\mathit{r}}}_{s_{\mathrm{4}}}\\ {\dot{\mathit{\theta }}}_{s_{\mathrm{4}}}\end{array}\right]\end{array}.\end{array}$$

The virtual power of a lifting load and the single lifting rope acting on the telescopic boom is as follows:

$$\begin{array}{}\text{(65)}& \begin{array}{rl}\mathit{\delta }{w}_{qs}& =\mathit{\delta }{\dot{\mathit{r}}}_{qg}^{\mathrm{T}}{\mathbf{F}}_{qg}+\mathit{\delta }{\dot{\mathit{r}}}_{qs}^{\mathrm{T}}{\mathbf{F}}_{qs}\\ & =\left[\begin{array}{cc}\mathit{\delta }{\dot{\mathit{r}}}_{s_{\mathrm{4}}}^{\mathrm{T}}& \mathit{\delta }{\dot{\mathit{\theta }}}_{s_{\mathrm{4}}}^{\mathrm{T}}\end{array}\right]\left({\mathbf{T}}_{r_{qg}}^{\mathrm{T}}{\mathbf{F}}_{qg}+{\mathbf{T}}_{r_{qs}}^{\mathrm{T}}{\mathbf{F}}_{qs}\right).\end{array}\end{array}$$

5 Calculation of the instability load

5.1 System equilibrium equations and tangent stiffness matrix

Taking the degrees of freedom of the boundary nodes in each substructure as system variables and dividing them according to translation and rotation, the global virtual power equations of the telescopic boom can be obtained as follows:

$$\begin{array}{}\text{(66)}& \begin{array}{rl}& \sum _k(\mathit{\delta }{\dot{\widehat{\mathit{u}}}}_k^{\mathrm{T}}{\stackrel{\mathrm{‾}}{\mathit{f}}}_k^n+\mathit{\delta }{\dot{\stackrel{\mathrm{‾}}{\mathit{\theta }}}}_k^{\mathrm{T}}{\stackrel{\mathrm{‾}}{\mathit{m}}}_k^n-\mathit{\delta }{\dot{\mathit{r}}}_{k\mathrm{0}}^{\mathrm{T}}{\mathbf{F}}_{k\mathrm{0}}-\mathit{\delta }{\dot{\mathit{\theta }}}_{k\mathrm{0}}^{\mathrm{T}}{\mathbf{M}}_{k\mathrm{0}})\\ & =\sum _i(\mathit{\delta }{\dot{\mathit{r}}}_i^{\mathrm{T}}{\mathit{f}}_i^a+\mathit{\delta }{\dot{\mathit{\theta }}}_i^{\mathrm{T}}{\mathit{m}}_i^a),\end{array}\end{array}$$

where ${\mathit{f}}_i^a$ and ${\mathit{m}}_i^a$ are external forces and moments acting on nodes in the global coordinate system. r_k0 and θ_k0 are the position vector and rotational angles of the origin of the kth substructure coordinate system in global coordinate system. F_k0 and M_k0 are the resultant force and moment of gravity and wind load of the kth substructure on its coordinate system origin. ${\stackrel{\mathrm{‾}}{\mathit{f}}}_k^n$ and ${\stackrel{\mathrm{‾}}{\mathit{m}}}_k^n$ are the equivalent nodal force and moment in the kth substructure coordinate system, which can be expressed as follows:

$$\begin{array}{}\text{(67)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{f}}}_k^n=\left({\mathbf{K}}_{ek}^{uu}{\widehat{\mathit{u}}}_k+{\mathbf{K}}_{ek}^{u\mathit{\theta }}{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_k\right)-({\mathbf{G}}_{gk}^u\stackrel{\mathrm{‾}}{\mathit{g}}+{\stackrel{\mathrm{‾}}{\mathit{f}}}_{wk})\text{(68)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{m}}}_k^n=\left({\mathbf{K}}_{ek}^{\mathit{\theta }u}{\widehat{\mathit{u}}}_k+{\mathbf{K}}_{ek}^{\mathit{\theta }\mathit{\theta }}{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_k\right)-\left({\mathbf{G}}_{gk}^{\mathit{\theta }}\stackrel{\mathrm{‾}}{\mathit{g}}+{\stackrel{\mathrm{‾}}{\mathit{m}}}_{wk}\right),\end{array}$$

where the submatrixes ${\mathbf{K}}_{ek}^{uu}$, ${\mathbf{K}}_{ek}^{u\mathit{\theta }}$, ${\mathbf{K}}_{ek}^{\mathit{\theta }u}$, and ${\mathbf{K}}_{ek}^{\mathit{\theta }\mathit{\theta }}$ can be obtained from the condensed substructure stiffness matrix K_e. ${\mathbf{G}}_{gk}^u$ and ${\mathbf{G}}_{gk}^{\mathit{\theta }}$ can be obtained from the gravity influence coefficient matrix G_g. ${\stackrel{\mathrm{‾}}{\mathit{f}}}_{wk}$ and ${\stackrel{\mathrm{‾}}{\mathit{m}}}_{wk}$ are the force and moment generated by the wind load condensed to the substructure nodes. ${\widehat{\mathit{u}}}_k$ and ${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_k$ are the displacements and rotational angles of the nodes at both end section origins in the local coordinate system.

According to the transformation relationship between the local and global coordinate system variables in Eq. (18), the transformation matrix of the displacements and the rotational angles of the kth substructure can be transformed into the following:

$$\begin{array}{}\text{(69)}& {\mathbf{T}}_{krg}=\left[\begin{array}{cccc}{\mathbf{T}}_{k{u}_{\mathrm{1}}}& {\mathbf{T}}_{k{u}_{\mathrm{2}}}& {\mathbf{T}}_{k{u}_{\mathrm{3}}}& {\mathbf{T}}_{k{u}_{\mathrm{4}}}\\ {\mathbf{T}}_{k{\mathit{\theta }}_{\mathrm{1}}}& {\mathbf{T}}_{k{\mathit{\theta }}_{\mathrm{2}}}& {\mathbf{T}}_{k{\mathit{\theta }}_{\mathrm{3}}}& {\mathbf{T}}_{k{\mathit{\theta }}_{\mathrm{4}}}\end{array}\right].\end{array}$$

Equation (66) can be converted into algebraic equations. The virtual power equations of the equivalent nodal forces and moments of the kth substructure can be transformed into the following:

$$\begin{array}{}\text{(70)}& \begin{array}{rl}\mathit{\delta }{\dot{\widehat{\mathit{u}}}}_k^{\mathrm{T}}{\stackrel{\mathrm{‾}}{\mathit{f}}}_k^n+\mathit{\delta }{\dot{\stackrel{\mathrm{‾}}{\mathit{\theta }}}}_k^{\mathrm{T}}{\stackrel{\mathrm{‾}}{\mathit{m}}}_k^n& =\mathit{\delta }{\dot{\mathit{r}}}_{k\mathrm{0}}^{\mathrm{T}}\left({\mathbf{T}}_{k{u}_{\mathrm{1}}}^{\mathrm{T}}{\stackrel{\mathrm{‾}}{\mathit{f}}}_k^n+{\mathbf{T}}_{k{\mathit{\theta }}_{\mathrm{1}}}^{\mathrm{T}}{\stackrel{\mathrm{‾}}{\mathit{m}}}_k^n\right)\\ & +\mathit{\delta }{\dot{\mathit{\theta }}}_{k\mathrm{0}}^{\mathrm{T}}\left({\mathbf{T}}_{k{u}_{\mathrm{2}}}^{\mathrm{T}}{\stackrel{\mathrm{‾}}{\mathit{f}}}_k^n+{\mathbf{T}}_{k{\mathit{\theta }}_{\mathrm{2}}}^{\mathrm{T}}{\stackrel{\mathrm{‾}}{\mathit{m}}}_k^n\right)\\ & +\mathit{\delta }{\dot{\mathit{r}}}_{kn}^{\mathrm{T}}\left({\mathbf{T}}_{k{u}_{\mathrm{3}}}^{\mathrm{T}}{\stackrel{\mathrm{‾}}{\mathit{f}}}_k^n+{\mathbf{T}}_{k{\mathit{\theta }}_{\mathrm{3}}}^{\mathrm{T}}{\stackrel{\mathrm{‾}}{\mathit{m}}}_k^n\right)\\ & +\mathit{\delta }{\dot{\mathit{\theta }}}_{kn}^{\mathrm{T}}\left({\mathbf{T}}_{k{u}_{\mathrm{4}}}^{\mathrm{T}}{\stackrel{\mathrm{‾}}{\mathit{f}}}_k^n+{\mathbf{T}}_{k{\mathit{\theta }}_{\mathrm{4}}}^{\mathrm{T}}{\stackrel{\mathrm{‾}}{\mathit{m}}}_k^n\right),\end{array}\end{array}$$

where the degrees of freedom of the kth substructure in the global coordinate system are $\left[{\mathit{r}}_{k\mathrm{0}};{\mathit{\theta }}_{k\mathrm{0}};{\mathit{r}}_{kn};{\mathit{\theta }}_{kn}\right]$.

The set of system variables can be defined as q. Combining Eqs. (66) and (70), the global virtual power equations of the telescopic boom can be expressed as follows:

$$\begin{array}{}\text{(71)}& \mathit{\delta }{\dot{\mathit{q}}}^{\mathrm{T}}\mathbf{F}\left(\mathit{q}\right)=\mathrm{0},\end{array}$$

where F(q) is the generalized nodal force matrix of the system.

According to the constraints and boundary conditions between the boom sections of the telescopic boom, the non-independent and independent degrees of freedom in the system variables are divided as $\mathit{q}=\left[{\mathit{q}}_{\mathrm{nd}};{\mathit{q}}_{\mathrm{id}}\right]$.

The constraint equations of the telescopic boom system are established by combining Eqs. (46)–(47), (49), and (51), which include the constraints between the boom sections, the constraints between the telescopic boom and the slewing table, and the constraints between the luffing cylinder. Without losing generality, the constraint equations of the telescopic boom can be expressed as follows:

$$\begin{array}{}\text{(72)}& \mathbf{\Phi }\left(\mathit{q}\right)=\mathbf{0}.\end{array}$$

The variation constraint equations derived from Eq. (72) are obtained as follows:

$$\begin{array}{}\text{(73)}& {\mathbf{G}}_{\mathrm{nd}}\mathit{\delta }{\dot{\mathit{q}}}_{\mathrm{nd}}+{\mathbf{G}}_{\mathrm{id}}\mathit{\delta }{\dot{\mathit{q}}}_{\mathrm{id}}=\mathbf{0}.\end{array}$$

The virtual variation in the non-independent degrees of freedom can be expressed by independent degrees of freedom, as follows:

$$\begin{array}{}\text{(74)}& \mathit{\delta }\dot{\mathit{q}}=\left[\begin{array}{c}-{\mathbf{G}}_{\mathrm{nd}}^{-\mathrm{1}}{\mathbf{G}}_{\mathrm{id}}\\ {\mathbf{E}}_{\mathrm{id}\times \mathrm{id}}\end{array}\right]\mathit{\delta }{\dot{\mathit{q}}}_{\mathrm{id}}\triangleq {\mathbf{\Gamma }}_d\mathit{\delta }{\dot{\mathit{q}}}_{\mathrm{id}}.\end{array}$$

Substituting Eq. (74) into Eq. (71), the transformed algebraic equation can be obtained from the following:

$$\begin{array}{}\text{(75)}& {\mathbf{F}}_{\mathrm{id}}={\mathbf{\Gamma }}_{\mathrm{id}}^{\mathrm{T}}\mathbf{F}\left(\mathit{q}\right)=\mathbf{0}.\end{array}$$

After the equations are transformed, the number of equations is the same as the number of independent variables, and the number of non-independent variables is the same as the number of constraint equations added in the system. Combining Eqs. (72) and (75), the system equations for solving the system variables can be obtained as follows:

$$\begin{array}{}\text{(76)}& \left\{\begin{array}{l}{\mathbf{\Gamma }}_{\mathrm{id}}^{\mathrm{T}}\mathbf{F}\left(\mathit{q}\right)=\mathbf{0}\\ \mathbf{\Phi }\left(\mathit{q}\right)=\mathbf{0}.\end{array}\right.\end{array}$$

The above formula is a set of highly nonlinear equations, and the corresponding tangent stiffness matrix is given, which can greatly improve the computational efficiency of the system.

In order to obtain the tangent stiffness matrix of the system equations, the derivative of Eq. (76) can be obtained as follows:

$$\begin{array}{}\text{(77)}& \begin{array}{rl}{\dot{\mathit{F}}}_{\mathrm{id}}& =\left({\mathbf{\Gamma }}_{\mathrm{id}}^{\mathrm{T}}\right(\partial \mathbf{F}\left(\mathit{q}\right)/\partial \mathit{q})\\ & +(\partial {\mathbf{\Gamma }}_{\mathrm{id}}^{\mathrm{T}}/\partial \mathit{q})\mathbf{F}\left(\mathit{q}\right))\dot{\mathit{q}}\triangleq ({\mathbf{\Gamma }}_{\mathrm{id}}^{\mathrm{T}}{\mathbf{G}}_{fq}+{\mathbf{G}}_{\mathrm{id}q}\mathbf{F}\left(\mathit{q}\right))\dot{\mathit{q}},\end{array}\end{array}$$

where

$$\begin{array}{}\text{(78)}& {\mathbf{G}}_{\mathrm{fq}}=\partial \mathbf{F}\left(\mathit{q}\right)/\partial \mathit{q},{\mathbf{G}}_{\mathrm{id}q}=\partial {\mathbf{\Gamma }}_{\mathrm{id}}^{\mathrm{T}}/\partial \mathit{q}.\end{array}$$

Equation (78) is related to the time derivative of the equivalent nodal force and moment in the coordinate system of the kth substructure and the time derivative of the resultant force and moment generated by the gravity and wind load at the origin of the local coordinate system.

$$\begin{array}{}\text{(79)}& {\displaystyle }{\dot{\stackrel{\mathrm{‾}}{\mathit{f}}}}_k^n={\mathbf{K}}_{ek}^{uu}{\dot{\widehat{\mathit{u}}}}_k+{\mathbf{K}}_{ek}^{u\mathit{\theta }}{\dot{\stackrel{\mathrm{‾}}{\mathit{\theta }}}}_k-{\mathbf{G}}_{gk}^u\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{\mathit{g}}}{T}_{k{\mathit{\omega }}_{\mathrm{0}}}{\dot{\mathit{\theta }}}_{k\mathrm{0}}\text{(80)}& {\displaystyle }{\dot{\stackrel{\mathrm{‾}}{\mathit{m}}}}_k^n=\left({\mathbf{K}}_{ek}^{\mathit{\theta }u}{\dot{\widehat{\mathit{u}}}}_k+{\mathbf{K}}_{ek}^{\mathit{\theta }\mathit{\theta }}{\dot{\stackrel{\mathrm{‾}}{\mathit{\theta }}}}_k\right)-{\mathbf{G}}_{gk}^{\mathit{\theta }}\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{\mathit{g}}}{T}_{k{\mathit{\omega }}_{\mathrm{0}}}{\dot{\mathit{\theta }}}_{k\mathrm{0}}\text{(81)}& {\displaystyle }{\dot{\mathbf{F}}}_{k\mathrm{0}}=-{\stackrel{\mathrm{̃}}{\mathbf{F}}}_{k\mathrm{0}}{\mathbf{T}}_{k{\mathit{\omega }}_{\mathrm{0}}}{\dot{\mathit{\theta }}}_{k\mathrm{0}},{\dot{\mathbf{M}}}_{k\mathrm{0}}=-{\stackrel{\mathrm{̃}}{\mathbf{M}}}_{k\mathrm{0}}{\mathbf{T}}_{k{\mathit{\omega }}_{\mathrm{0}}}{\dot{\mathit{\theta }}}_{k\mathrm{0}}.\end{array}$$

Substituting the generalized force matrix F(q) into Eq. (78) yields the following:

$$\begin{array}{}\text{(82)}& \begin{array}{rl}(\partial {\mathbf{\Gamma }}_{\mathrm{id}}^{\mathrm{T}}/\partial \mathit{q})\mathbf{F}\left(\mathit{q}\right)& =-(\partial {\mathbf{G}}_{\mathrm{id}}^{\mathrm{T}}/\partial \mathit{q}){\mathbf{G}}_{\mathrm{nd}}^{-\mathrm{T}}{\mathit{f}}_{\mathrm{nd}}\\ & -{\mathbf{G}}_{\mathrm{id}}^{\mathrm{T}}(\partial {\mathbf{G}}_{\mathrm{nd}}^{-\mathrm{T}}/\partial \mathit{q}){\mathit{f}}_{\mathrm{nd}},\end{array}\end{array}$$

where f_nd is the corresponding submatrix of the non-independent degrees of freedom in the generalized force matrix F(q).

The time derivative of the constraint equations in Eq. (76) can be obtained as follows:

$$\begin{array}{}\text{(83)}& \dot{\mathbf{\Phi }}={\mathbf{G}}_{\mathrm{nd}}{\dot{\mathit{q}}}_{\mathrm{nd}}+{\mathbf{G}}_{\mathrm{id}}{\dot{\mathit{q}}}_{\mathrm{id}}=\left[\begin{array}{cc}{\mathbf{G}}_{\mathrm{nd}}& {\mathbf{G}}_{\mathrm{id}}\end{array}\right]\left[\begin{array}{c}{\dot{\mathit{q}}}_{\mathrm{nd}}\\ {\dot{\mathit{q}}}_{\mathrm{id}}\end{array}\right].\end{array}$$

Substituting Eqs. (78)–(82) into Eq. (77) and combining the results with Eq. (83) above, the tangent stiffness matrix of the system can be obtained.

5.2 Differential form of equilibrium equations and instability load

The instability load of the telescopic boom of all-terrain cranes, especially for the medium and long booms, is a key indicator in the lifting ability. R(q) is obtained by combining Eqs. (57), (60)–(61), and (76), which consists of the guy rope force, the back pull plate force, the generalized node internal force, and system constraint equations. The lifting load is the variable of external force in the system equations, and the related load is in Eq. (65). As shown in Fig. 15, since the lifting load and the single lifting rope load do not directly act on the node when the lifting load is a unit load, the generalized force related to the lifting load can be obtained, i.e., the generalized force matrix G₀ corresponding to the unit lifting load. The nonlinear equilibrium equations of the system can be established by R(q) and the generalized force matrix G₀ related to the lifting load through the virtual power equations of the assembly system. To sum up, the telescopic boom meets the structural equilibrium under any lifting load, and the nonlinear equilibrium equations of the telescopic boom can be expressed as follows:

$$\begin{array}{}\text{(84)}& \mathbf{R}\left(\mathit{q}\right)+\mathit{\lambda }{\mathbf{G}}_{\mathrm{0}}=\mathbf{0},\end{array}$$

where q is the system variable, including nodal displacements and rotational angles of all nodes in the telescopic boom, and λ is the lifting load parameter.

For Eq. (84), in order to obtain the instability load under the systematic consideration of structural nonlinear effects, the conventional method is to divide the load into multiple load steps, through the incremental method, and obtain the corresponding node displacement under each load through calculation, so as to establish the corresponding equilibrium path curve. During the solution process, we monitor the curve slope under all loads, i.e., $\frac{\mathrm{d}\mathit{q}}{\mathrm{d}\mathit{\lambda }}$. Theoretically, when the slope is infinite, then the load corresponding to this extreme point is the critical instability load of the structure. Therefore, it is necessary to set a threshold value as the judgment criterion of the instability load of the telescopic boom system, so as to obtain the instability load.

$$\begin{array}{}\text{(85)}& max\left(∥{\displaystyle \frac{\mathrm{d}\mathit{q}}{\mathrm{d}\mathit{\lambda }}}∥/∥{\displaystyle \frac{\mathrm{d}{\mathit{q}}_{\mathrm{0}}}{\mathrm{d}{\mathit{\lambda }}_{\mathrm{0}}}}∥\right)\ge K_{\mathrm{0}},\end{array}$$

where K₀ is the ratio of the current state parameters to the initial state parameters, i.e., the slope ratio of the equilibrium path curve.

Combined with the above contents, the key point is to obtain the slope of the equilibrium path curve in the process of the changing load, and then determine the instability load in combination with the judgment criterion. In the process of searching for the instability load, the lifting load parameter λ is an unknown variable, and then the derivative of λ in Eq. (84) can be obtained to obtain the corresponding differential equations, as follows:

$$\begin{array}{}\text{(86)}& \left[{\displaystyle \frac{\partial \mathbf{R}\left(\mathit{q}\right)}{\partial \mathit{q}}}\right]{\displaystyle \frac{\mathrm{d}\mathit{q}}{\mathrm{d}\mathit{\lambda }}}+{\mathbf{G}}_{\mathrm{0}}=\mathbf{0},\end{array}$$

where $\left[\frac{\partial \mathbf{R}\left(\mathit{q}\right)}{\partial \mathit{q}}\right]$ is the tangent stiffness matrix of the telescopic boom, which has been obtained in Sect. 5.1. $\frac{\mathrm{d}\mathit{q}}{\mathrm{d}\mathit{\lambda }}$ can reflect the nonlinearity of the load displacement curve under the current load.

Furthermore, the first derivative of system variables to the lifting load can be obtained by Eq. (86), as follows:

$$\begin{array}{}\text{(87)}& {\displaystyle \frac{\mathrm{d}\mathit{q}}{\mathrm{d}\mathit{\lambda }}}=-{\left[{\displaystyle \frac{\partial \mathbf{R}\left(\mathit{q}\right)}{\partial \mathit{q}}}\right]}^{-\mathrm{1}}{\mathbf{G}}_{\mathrm{0}}.\end{array}$$

With the advantage of an automatic step adjustment in ODE113, a conventional differential equation solver of MATLAB, Eq. (87) is solved to realize the function of automatically adjusting the step size according to the nonlinear degree corresponding to the current load state of the system with initial value, λ=0. When the structure is in the linear stage, the solution step can be automatically increased. When the structure is in the nonlinear stage, the step can be automatically reduced. Under the premise of ensuring the convergence of each step of solution, it can quickly track the equilibrium path and search for the instability load. Compared with the traditional incremental method and arc length method, this method can greatly reduce the number of solutions and improve the calculation efficiency.

Table 1 Length parameters of boom sections.

Table 2 Section parameters comparison of boom sections.

Figure 10 Cantilever beam model.

Figure 11 Telescopic boom combination modes.

Figure 12 Load–displacement relationship of right end node.

6 Numerical examples

The relevant numerical examples in this section are calculated based on the actual structural parameters of a certain type of all-terrain crane. The Young modulus of elasticity of boom sections material is $E=\mathrm{2.1}\times {\mathrm{10}}^{\mathrm{11}}$ Pa, material density is $\mathit{\rho }=\mathrm{7.85}\times {\mathrm{10}}^{\mathrm{3}}$ kg m⁻³, and Poisson's ratio is ν=0.3. The telescopic boom is composed of eight boom sections, each boom section contains one boom pin and four pin holes, and their positions along the length of the boom section are recorded as L₀, L₁, L₂, L₃, and L₄, as shown in Table 1.

In the commercial software Ansys19.1, the BEAM188 element is used to customize the section according to the section parameters of each boom section, and the section moment of inertia and area are calculated. The section parameters calculated in this paper, and their comparison, are shown in Table 2.

6.1 Example 1 of a cantilever beam about the boom section

In this example, a cantilever beam model is taken as the calculation object for relevant calculation comparison, as shown in Fig. 10. The combination of different lengths of the telescopic boom is realized by different connection modes of the boom pin and the pin hole between boom sections. The super element of the same boom section with different combinations is different. Four connection modes are defined here, corresponding to the lengths of different telescopic boom combinations shown in Fig. 11.

The boom section 1 is taken as calculation model, which is constrained as a cantilever beam. The end load is applied to the origin of the section coordinate system. The z-direction displacement and rotational angle around the y axis of the right end node of four combination modes are calculated under different super elements, as shown in Fig. 11.

Through the comparison of calculation results in Fig. 12, it can be seen that, for the same boom section, when different super elements are used, the error in calculation results is very small, and the calculation results are consistent.

6.2 Example 2 of cantilever beam about boom section

Taking boom section 7 as the calculation model, the division of a super element is established by the combination of mode 3, which is divided into three super elements with a total of 24 degrees of freedom. As shown in Fig. 11, the load application mode and constraint form are established. A model consistent with its section and length parameters is established in the Ansys software, which is divided into 15 elements with a total of 96 degrees of freedom, according to the basic length of the subelement in a super element of this paper. The two model diagrams are shown in Fig. 13.

Figure 13 Ansys and present finite element model of boom section 7.

Different discrete loads are applied at the section centroid of the right end section of the model. Considering the geometrical nonlinear effect of the structure, Ansys software sets the load step as 20 in a large displacement calculation. The vertical z-direction displacements and rotational angles around the y axis of the end nodes of the two models are calculated. The load displacements and rotational angles comparison curves are shown in Fig. 14.

Figure 14 Force–displacement relationship of the right end section.

Figure 15 Telescopic boom node number.

Figure 16 Nodal displacements of telescopic boom under different preloads.

The calculation model is divided into several substructures, the calculation degrees of freedom are reduced by static condensation, and the co-rotational formulation method is used to consider the geometrical nonlinear effect of the whole structure, which improves the efficiency of numerical calculation. At the same time, through a curve comparison, it can be seen that the calculation results obtained by the two different models have small errors, and the errors in the displacements and rotational angles can meet the requirements of the calculation errors in engineering.

6.3 Instability load of the all-terrain crane with a super lift system

In this example, eight boom sections are assembled into a telescopic boom, with the length being 55.5 m in accordance with connection mode 2. The node number of the telescopic boom is shown in Fig. 15. The preload of the guy rope is applied on the telescopic boom to resist the deformation of self-weight before lifting the load. From the engineering point of view, the larger the initial luffing angle of the telescopic boom, the smaller the preload of the guy rope to overcome the deformation will be. In this example, the luffing angle of the telescopic boom is 83^∘.

The preload is discretized from 10 to 40 kN into four preloads for application, and the nodal displacements of the telescopic boom relative to the horizontal and vertical directions of the initial position under different preload conditions are obtained, as shown in Fig. 16.

Figure 17 Load displacements of markers on a telescopic boom.

Figure 18 Instability loads of a telescopic boom under different preloads.

Figure 19 Telescopic boom deformation in the luffing plane and slewing plane.

Combining with the engineering practice, taking 20 kN as the preload of the guy rope, the original length of the guy rope is calculated, and the subsequent calculation is carried out on this basis. In Fig. 15, the displacements of marker 1 and marker 2 points change with the increasing load, and the final load displacement curves are shown in Fig. 17. It can be seen from the curve that the lifting load is 6.9×10⁵ kg, and the slope of the load displacement curve increases sharply at that point, which can be considered to be the critical instability load of the telescopic boom.

In practical engineering, the critical instability load cannot be used as the rated load of the all-terrain crane. During the calculation, we monitor the slope change in the load displacement curve of marker 1 and marker 2 points and take the ratio of the real-time slope to the initial slope as the judgment criterion, which is 3 in this calculation example. The most reasonable value of this criterion needs to be selected through a large number of calculations and tests. The instability load curve of the telescopic boom with different preloads, and the range of luffing angle at 70 to 85^∘, is shown in Fig. 18.

It can be seen from Fig. 18 that, with the increase in preload in the initial state, the deformation of telescopic boom in the initial state is reduced to a certain extent, and even the reverse bending deformation appears. At the same time, the axial load of the telescopic boom increases, and the instability load shows a downward trend in Fig. 18. For the lifting capacity determined by the whole telescopic boom, it also includes the load determined by the structural strength of the boom. The final lifting capacity is data determined by comprehensive factors, and the structural strength load will be studied in follow-up research.

The preload of guy rope is 20 kN, and three lifting loads are selected for calculation, with the lateral load ratio being 3 % of the lifting load, which acts on the boom head. The deformations in the luffing plane and slewing plane are shown in Fig. 19.

7 Conclusion

The instability load of the telescopic boom of an all-terrain crane can be obtained based on the differential form of system governing equations. First, for each boom section, the corresponding substructure is established by selecting reasonable boundary nodes, and the internal degrees of freedom of the substructure are reduced by static condensation method to form a super element. The approximately specific beam relationship between the corresponding boom sections and the boundary conditions of the telescopic boom are given. Second, based on the proposed geometrical nonlinear calculation method of the cable element, the nonlinear external force at the guy rope connection node of the telescopic boom with initial preload is presented. The equilibrium equations of telescopic boom with the control parameters of load are derived based on co-rotational procedures, and the tangent stiffness of the equilibrium equations are formulated. Finally, the equilibrium equations are transformed into differential form, and the load displacement curves are illustrated by solving the differential equations with existing numerical methods. A method to calculate the structural equilibrium path and instability load of the telescopic boom of the all-terrain crane with a guy rope preload is given, which provides a certain theoretical support for the design of the all-terrain crane.

Appendix A: Nomenclature

${\mathit{g}}_{\mathrm{1}},{\mathit{g}}_{\mathrm{2}},{\mathit{g}}_{\mathrm{3}}$	Base vectors of the global coordinate system for telescopic boom
${\mathit{e}}_{\mathrm{1}},{\mathit{e}}_{\mathrm{2}},{\mathit{e}}_{\mathrm{3}}$	Base vectors of the local coordinate system for a substructure
g	Gravity acceleration in global coordinate system
r1,r2	Current global position vectors of the two end-section origins of a substructure
θ1,θ2	Global rotational angles of the two end-sections of a substructure (Cardan angles)
${\widehat{\mathit{r}}}_{\mathrm{1}},{\widehat{\mathit{r}}}_{\mathrm{2}}$	Initial position vectors of the two end-section origins for a substructure in local coordinate system
${\widehat{\mathit{u}}}_{\mathrm{1}},{\widehat{\mathit{u}}}_{\mathrm{2}}$	Translational displacement vectors of the two end-section origins for a substructure in local coordinate system
ω1,ω2	Global angular velocities of the two end-sections of a substructure
${\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}$	Angular velocities in local coordinate system of a right section of a substructure
${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{1}},{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{2}}$	Rotational angle vectors of two end-sections of a substructure in a local coordinate system
R1,R2	Global rotational matrices of the two end-sections of a substructure
${\stackrel{\mathrm{‾}}{\mathbf{R}}}_{\mathrm{2}}$	Transformation matrix between a right section and a local coordinate system of a substructure
${\widehat{\mathit{u}}}_{e\mathrm{1}},{\widehat{\mathit{u}}}_{e\mathrm{2}}$	Local displacement vectors of the two section origins of a beam element in beam coordinate system
${\stackrel{\mathrm{‾}}{\mathit{u}}}_{e\mathrm{1}},{\stackrel{\mathrm{‾}}{\mathit{u}}}_{e\mathrm{2}}$	Local translational displacement vectors of the two-section centroid of a beam element in a beam coordinate system
${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{e\mathrm{1}},{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{e\mathrm{2}}$	Local rotational angle vectors of the two end-sections of a beam element in a beam coordinate system
${\stackrel{\mathrm{‾}}{\mathit{k}}}_e,{\widehat{\mathit{k}}}_e$	Element stiffness matrix based on a section centroid and section origin parameters
Ke	Stiffness matrix of a super element
${\mathit{r}}_e,{\stackrel{\mathrm{‾}}{\mathit{r}}}_e$	Current global and initial local position vectors of any point in a beam element
me,Le	Mass per unit length and length of a beam element
Gg	Gravity influence coefficient matrix
G1,G2	Wind load influence coefficient matrix
${\mathit{b}}_{\mathrm{1}}^i,{\mathit{b}}_{\mathrm{2}}^i,{\mathit{b}}_{\mathrm{3}}^i$	Base vectors of the local coordinate system of a boom section

${\mathbf{T}}_{{\mathit{\omega }}_{si}}$	Overlap point angular velocities coefficient matrix of boom section overlap point
${\mathit{e}}_b,{\mathit{e}}_{wi},{\mathit{e}}_s$	Telescopic boom, Y-shaped bracket, and lifting rope axial base vectors
q	Telescopic boom system variables with independent and dependent variables
Γid	Transformation matrix between independent variables and system variables
Φ(q)	Constraint equations of the telescopic boom system
G0	Generalized force matrix corresponding to the unit lifting load
λ	Lifting load parameter