Sage Journals: Discover world-class research

Abstract

This article presents a numerical study of the stability of a hydraulic excavator during performing lifting operations. A planar dynamic model is developed with six degrees of freedom, which considers the base body elastic connection with the terrain, the front digging manipulator links, and the presence of the freely suspended payload. Differential equations describing the excavator dynamic behavior are obtained by using the Lagrange formalism. Numerical experiments are carried out to study the excavator dynamic stability under different operating conditions during the motion along a vertical straight-line trajectory. It is shown that the arising inertial loads during the movement of the links along the vertical trajectory, combined with the payload swinging and the motion of the base body, decreases the excavator stability. It was found that the excavator stability during following vertical straight-line trajectory decreases considerably in the lower part of the vertical trajectory. If the stability coefficient is close to 1, the payload swinging can cause the separation of a support from the terrain; nevertheless, the excavator stability can be restored. A method for tire stiffness and damping coefficients estimation is presented. The validation of the dynamical model is performed by the use of a small-scale elastically mounted manipulator.

Keywords

Hydraulic excavator lifting operations modeling overturning stability

Introduction

Due to the considerable efforts in the industry in the past decades, the hydraulic excavators have transformed from specialized heavy-duty earth-moving machines, whose primary function is the digging, into multifunctional devices. By the use of specialized attachments and working tools, they can perform a wide variety of technical tasks that are specific for other types of machines. This has allowed many operations at the construction and mining sites to be performed by a smaller number of machines and, thus, to increase the economic efficiency.

One of the widespread activities is the use of the hydraulic excavators as cranes.¹ The typically performed operations are lifting, moving, and positioning of freely suspended payloads, such as concrete pipes and beams, industrial equipment, and construction elements. The payload is attached by slings to the lifting points, hooks, lifting eyes, lashing points, or other special attachment that has been provided and approved by the excavator manufacturer.

Despite the undeniable convenience when using an excavator as a crane, the lifting operations performed have some peculiarities, leading to the reduction of the excavator exploitation reliability. Basically, excavators are not designed for lifting operations and the operator is trained to perform earth-moving operations but not to manipulate freely suspended payloads. Even in the case of a very skilled operator, this leads to an undesirable payload swinging and, as a consequence, to the increase of dynamical loads, which negatively influences the excavator exploitation reliability. The reduced ability of the excavators to perform lifting operations increases the risk of accidents with the staff near the excavator, and especially in the vicinity of the freely suspended payload. The statistical reports show a considerable risk of injuries caused by the strike by the bucket or swinging boom.² Edwards and Holt³ have studied the hazards associated with using construction excavators as cranes and have investigated the catastrophic failure of a link. The same authors have studied excavator overturn accidents.⁴ According to Lim et al.,⁵ 38% of the excavator accidents in Korea are due to performing lifting operations. Some requirements for the safety of the excavators are presented in the standards.^6,7

The literature survey reveals different aspects of the excavator exploitation as a crane. Some documents emphasize on the practical aspects of the problem,⁸ while other papers aim at the design and development of a research methodology and simulation models. Lim et al.⁵ presented the application of the Zero Moment Point theory to the computation of the excavator stability during performing lifting operations. Yu et al.⁹ introduced the Static Compensation Zero Moment Point algorithm that has a more accurate output compared to the general algorithm. Heep et al.¹⁰ proposed a method for payload estimation by the use of a multibody simulation.

Reduced overturning stability is the most critical indication of reduced exploitation reliability as it may even lead to accidental situations. The problem of using excavators as cranes and investigations of their overturning stability during performing lifting operations in the vertical plane of the front manipulator motion has not received enough attention over the years, mainly because such type of work operations are not typical for this type of machines, but are performed in practice. In this aspect, more attention has been given to some other machines that are similar to the excavators. Although ensuring the stability of the mobile machinery is an up-to-date subject and intensive research is carried out to improve it, much of the recent research does not take into account the dynamic nature of the loading due to the payload swinging, elastic supports, and moving links.^11,12 Rauch et al.¹³ developed a dynamic model of a mobile boom crane and studied the tip-over stability by a simplified semi-dynamic approach to account for the payload swing and inertia forces. Since the dynamic stability is influenced by the foundation–soil interaction, Messioud et al.¹⁴ used complex mathematical models, but Towarek¹⁵ obtained adequate results for the dynamic stability of a boom crane by the use of a two-parameter rheological model for the representation of the soil. Abo-Shanab and Sepehri¹⁶ presented the development of a model that can adequately simulate the dynamic stability of manipulators mounted on moveable platforms. Janosevic et al.¹⁷ investigated the influence of the hydraulic transmission oil temperature on the excavator dynamic stability. The presented mathematical model takes into account the hydraulic cylinders compliance due to the hydraulic fluid compressibility. Another research by Mitrev et al.¹⁸ resulted in a complex multibody model of a classical hydraulic excavator, based on the use of Lagrange multipliers, which allows the study of dynamic behavior during digging operations. The approaches for dynamic modeling and control of different types of cranes were presented in the very detailed and extensive study by Abdel-Rahman et al.¹⁹ Fodor et al.²⁰ presented a very interesting approach based on an input shaping control technique for oscillation reduction in the vertical plane of a forestry crane in order to improve the overturning stability of the machine.

Along with the widely used free or commercial computer codes for simulation of multibody systems, the dynamic models of a wide range of physical systems are formulated and studied by use of analytical^21–23 or numerical techniques.²⁴ Setting the mathematical model represents the very core of dynamic modeling of such type of systems composed of hydraulic and mechanical subsystems.²⁵ Recent investigations^26,27 show the feasibility of such an approach to study dynamic behavior of excavators, despite the fact that in certain cases, the resulting dynamic equations are quite demanding to resolve. However, the literature review also shows that the available dynamical models of excavators and similar mobile machines cannot be used reliably for all types of performed operations because the models do not include the peculiarities of the hydraulic excavator kinematic structure in performing lifting operations and especially the presence of a freely suspended payload attached to the bucket and elastic support. For this reason, the main aim of this article is to numerically investigate the dynamic behavior and overturning stability of the hydraulic excavator during performing lifting operations. For this purpose, kinematic and dynamic models of the excavator are developed, and numerical experiments are carried out to study dynamic behavior of the system in different operating conditions.

This article is organized as follows: “Introduction” section lists some peculiarities during the performance of lifting operations by hydraulic excavators. Some studies concerning different aspects of the dynamic modeling of the excavators and similar machines are reported. “Kinematic model and trajectory planning” section presents the solution to the inverse kinematics task and describes the trajectory planning. “Dynamic model” section is focused on dynamic modeling of the excavator and the derivation of the system of differential equations. “Overturning stability of the excavator” section points out some peculiarities concerning the excavator overturning stability. The results of the numerical experiments are shown in “Numerical experiments and discussions” section. In “Parameter estimation and validation of the mathematical model” section, a method for determination of the tire characteristics and experimental model validation are presented. The article concludes with a summary in “Conclusion” section.

Kinematic model and trajectory planning

Figure 1 shows a hydraulic excavator performing typical lifting operation—lifting/lowering of the payload R along the vertical straight-line trajectory. The cutting edge, denoted by H, serves as a reference point for the start $(x_{H,} y_{H}^{s})$ and end $(x_{H,} y_{H}^{f})$ points of the vertical lifting trajectory. The payload is attached to the lifting point P by slings, allowing its free swinging in the vertical plane. There are two options for controlling the movement of the front digging manipulator.

Figure 1.

Schematic view of an excavator performing lifting operations.

In the first case, the cutting edge H moves along a straight line and the angle of the bucket tilt is kept constant relative to the stick. In this case, the manipulator is considered as a system with two controllable degrees of freedom—the boom and the stick together with the bucket—and one uncontrollable degree of freedom—the free payload swinging.

The motion of the cutting edge along the trajectory leads to the change of the bucket orientation φ with respect to the horizontal line (see Figure 1), resulting in a change of direction of the payload weight force according to the lifting point causing its unfavorable loading. The other drawback due to the bucket rotation is the occurrence of the horizontal movement of the lifting point P which induces a payload swinging and, as a consequence, poor positioning of the payload and increased risk of accidents, especially in confined spaces.

In the second case, during the motion along the vertical trajectory, the bucket hydraulic cylinder is additionally controlled to maintain a constant bucket angle φ. The system has four degrees of freedom, three of which are controlled and one is uncontrolled. The above described motion along a predefined trajectory can be considerably facilitated by the use of the widely used²⁸ and well-studied excavator automatic digging system.^29–31

The current article investigates the second case. The values of the angles $θ_{2}$ , $θ_{3}$ , and $θ_{4}$ (see Figure 1) determine the horizontal and vertical position of the cutting edge in the plane. In order to find the angle values corresponding to the motion of the cutting edge along the straight line with start point with coordinates $(x_{H,} y_{H}^{s})$ and end point with coordinates $(x_{H,} y_{H}^{f})$ , one must solve the inverse kinematics task for a three link planar manipulator.³² The angular velocities and accelerations of the links are determined by using the mechanical system Jacobian. Due to the constant orientation of the bucket $\overset{\cdot}{φ} = \overset{\cdot\cdot}{φ} = 0$ , the velocity of the lifting point P is the same as the velocity of the cutting edge H. The determination of the values of the angles $θ_{2}$ , $θ_{3}$ , and $θ_{4}$ and their derivatives as a function of the cutting edge position and bucket orientation completely ends the inverse kinematics task solution.

The goal of the trajectory planning is to transport the payload from the initial point $y_{H}^{s}$ to the final point $y_{H}^{f}$ (see Figure 1) keeping the motion of the bucket reference point H along a vertical straight line with a predefined constant position $x_{H}$ inside the excavator workspace. The desired trajectory for the vertical motion of point H is proposed in the form of a fifth-order polynomial.³² Imposing zero initial and final conditions for the velocity and acceleration of the point H, the polynomial takes the following form

y_{H}^{d} (t) = y_{H}^{s} + \frac{t^{3} (y_{H}^{f} - y_{H}^{s}) (6 t^{2} - 15 t t_{f} + 10 t_{f}^{2})}{t_{f}^{5}}

(1)

where t_f denotes the duration of the motion along the trajectory.

Dynamic model

Figure 2 shows a geometrical layout of the excavator with the consideration of the payload swinging. The excavator mechanical system is modeled as an open kinematic chain undergoing planar motion and consisting of a base body (composed of traveling and swing bodies, represented as a single body), boom, stick, bucket, and payload. The links are interconnected by rotational joints and are subjected to forces generated by the hydraulic actuators and gravity. Each link is characterized by its geometrical and inertia parameters. The base body is connected to the ground by spring-damper elements which represent the elastic and damping properties of the excavator compliant tires, outriggers, or undercarriage. The presented system has six degrees of freedom, and thus, six generalized coordinates are used to define the system geometrical configuration. The base body can move in the vertical direction—generalized coordinate s (heave)—and rotate according to its mass center C₁—generalized coordinate θ₁ (pitch). The boom, stick, bucket, and payload can rotate according to the connecting joints (points E, F, G, and P)—the corresponding generalized coordinates are θ₂, θ₃, θ₄, and θ₅. All angles are measured counterclockwise.

Figure 2.

Dynamic model of the hydraulic excavator.

An up-to-date analysis of methods suitable for derivation of the mathematical models of cranes is presented in Abdel-Rahman et al.¹⁹ In order to derive a dynamic model suitable for the study of the lifting operations, in this article, we make the following assumptions: (1) the payload is considered as a point mass; (2) the stiffness of the excavator elements is neglected, and they are considered as rigid bodies; (3) the inertia properties of the hydraulic cylinders and linkages are included into the inertial parameters of the boom, stick, and bucket; (4) the backlashes and the friction forces in the joints are neglected; and (5) the compliance of the hydraulic cylinders is neglected.

The dynamic equations of motion of the excavator during the lifting operations are derived in a systematic way using the Lagrange formalism³³

\frac{d}{dt} \frac{\partial L_{a}}{\partial {\overset{\cdot}{q}}_{i}} - \frac{\partial L_{a}}{\partial q_{i}} + \frac{\partial D}{\partial {\overset{\cdot}{q}}_{i}} = Q_{i} (i = 1 / 6)

(2)

where the Lagrangian L_a represents the difference between the kinetic T and potential P energies of the system studied, D denotes the dissipation energy of the system, and Q_i are the generalized forces associated with the generalized coordinates.

The kinetic and potential energies of the system are positive definite quadratic forms of the, respectively, generalized velocities $\overset{\cdot}{q} = [\begin{matrix} {\overset{\cdot}{θ}}_{1} & {\overset{\cdot}{θ}}_{2} & {\overset{\cdot}{θ}}_{3} & {\overset{\cdot}{θ}}_{4} & {\overset{\cdot}{θ}}_{5} & \overset{\cdot}{s} \end{matrix}]^{T}$ and generalized coordinates $q = [\begin{matrix} θ_{1} & θ_{2} & θ_{3} & θ_{4} & θ_{5} & s \end{matrix}]^{T}$ .

Kinetic energy of the system

The total kinetic energy of the system comprises the kinetic energies of the base body T₁, boom T₂, stick T₃, bucket T₄, and swinging payload T₅. Using the notations in Figure 2 and according to König’s theorem, the total kinetic energy of the mechanical system is obtained as follows

T = \sum_{i = 1}^{5} T_{i}

(3)

where $2 T_{1} = m_{1} {\overset{\cdot}{s}}^{2} + J_{1} {\overset{\cdot}{θ}}_{1}^{2}$ , $2 T_{2} = m_{2} ({\overset{\cdot}{x}}_{C_{2}}^{2} + {\overset{\cdot}{y}}_{C_{2}}^{2}) + J_{2} {\overset{\cdot}{θ}}_{12}^{2}$ , $2 T_{3} = m_{3} ({\overset{\cdot}{x}}_{C_{3}}^{2} + {\overset{\cdot}{y}}_{C_{3}}^{2}) + J_{3} {\overset{\cdot}{θ}}_{123}^{2}$ , $2 T_{4} = m_{4} ({\overset{\cdot}{x}}_{C_{4}}^{2} + {\overset{\cdot}{y}}_{C_{4}}^{2}) + J_{4} {\overset{\cdot}{θ}}_{1234}^{2}$ , and $2 T_{5} = m_{5} ({\overset{\cdot}{x}}_{R}^{2} + {\overset{\cdot}{y}}_{R}^{2})$ . Equation (12) uses the following short notations: ${\overset{\cdot}{θ}}_{12}^{2} = ({\overset{\cdot}{θ}}_{1} + {\overset{\cdot}{θ}}_{2})^{2}$ , ${\overset{\cdot}{θ}}_{123}^{2} = ({\overset{\cdot}{θ}}_{1} + {\overset{\cdot}{θ}}_{2} + {\overset{\cdot}{θ}}_{3})^{2}$ , and so on; ${\overset{\cdot}{x}}_{Ci}$ and ${\overset{\cdot}{y}}_{Ci}$ denote the absolute velocities of the mass centers of the links. Similarly, ${\overset{\cdot}{x}}_{R}$ and ${\overset{\cdot}{y}}_{R}$ denote the absolute velocities of the payload. All velocities are obtained by a time differentiation of the corresponding position vectors according to the global coordinate system (CS).

As is shown in Figure 2, a local CS {x_iy_i} is connected to every link of the kinematic chain. The CS {x₁y₁} is connected to the gravity center C₁ of the base body, the CS {x₂y₂}, {x₃y₃}, {x₄y₄}, and {x₆y₆} are connected to the joints E, F, G, and P, respectively. The transition between the local CS is performed by the use of transformation matrices in the following form

{}_{i}^{j}T (α, a_{x}, a_{y}) = [\begin{matrix} \cos α & - \sin α & a_{x} \\ \sin α & \cos α & a_{y} \\ 0 & 0 & 1 \end{matrix}]

(4)

where i and j denote the numbers of the adjacent CSs, α is the angle between the CSs, $[\begin{matrix} a_{x} & a_{y} \end{matrix}]^{T}$ represents the vector from the origin of CS i to the origin of CS j expressed in the CS i. The kinetic energy of the system depends on the system geometrical configuration and the velocities of the links.

Potential and dissipative energies of the system

The potential energy of the system is computed as a sum of the potential energies of the links P_G and the deformations of the springs P_k

P = P_{G} + P_{k}

(5)

The potential energy of the spring deformation is easily determined by use of Figure 3, where the dashed line depicts the position of the base body when the springs are undeformed. In this position, the gravity center C₁ of the base body coincides with the beginning O of the global CS. The normal line shows the position of the base body when the springs are deformed due to the vertical translation s and rotation θ₁.

Figure 3.

Scheme for determining spring deformation.

If Δy_B denotes the vertical spring deformation in point B, then

Δ y_{B} = {{}^{0}y}_{B}^{″} - {{}^{0}y}_{B}^{'}

(6)

where ${{}^{0}y}_{B}^{″}$ and ${{}^{0}y}_{B}^{'}$ are y coordinates of point B in the global CS for a deformed and undeformed spring, respectively.

Denoting the local coordinates of point B in the base body local CS by ${{}^{1}x}_{B}$ and ${{}^{1}y}_{B}$ and taking into account that

{{}^{0}y}_{B}^{″} = s + {{}^{1}y}_{B} \cos θ_{1} + {{}^{1}x}_{B} \sin θ_{1}

(7)

and

{{}^{0}y}_{B}^{'} = {{}^{1}y}_{B}

(8)

then for the spring deformation, one can write

Δ y_{B} = s + {{}^{1}y}_{B} (\cos θ_{1} - 1) + {{}^{1}x}_{B} \sin θ_{1}

(9)

In a similar manner, for the deformation of the spring in point A, we have

Δ y_{A} = s + {{}^{1}y}_{A} (\cos θ_{1} - 1) + {{}^{1}x}_{A} \sin θ_{1}

(10)

Then, the total potential energy of the spring deformation is

2 P_{k} = k (Δ y_{B}^{2} + Δ y_{A}^{2})

(11)

The potential energy of a single excavator link is computed as the work required to raise its gravity center according to the x axis of the global CS, that is, it is a function of the system geometrical configuration. The potential energy of the links is the sum of the potential energies of the base body $P_{G}^{1}$ , boom $P_{G}^{2}$ , stick $P_{G}^{3}$ , bucket $P_{G}^{4}$ , and payload $P_{G}^{5}$ and is given as

P_{G} = \sum_{i = 1}^{5} P_{G}^{i}

(12)

where $P_{G}^{1} = m_{1} gs$ , $P_{G}^{2} = m_{2} g (L_{1 x} s_{1} + L_{1 y} c_{1} + L_{C_{2} x} s_{12} + L_{C_{2} y} c_{12} + s)$ , $P_{G}^{3} = m_{3} g (L_{1 x} s_{1} + L_{1 y} c_{1} + L_{2 x} s_{12} + L_{C_{3} x} s_{123} + L_{C_{3} y} c_{123} + s)$ , $P_{G}^{4} = m_{4} g (s + L_{1 x} s_{1} + L_{1 y} c_{1} + L_{2 x} s_{12} + L_{3 x} s_{123} + L_{C_{4} x} s_{1234} + L_{C_{4} y} c_{1234})$ , $P_{G}^{5} = m_{5} g (L_{1 x} s_{1} + L_{1 y} c_{1} + L_{2 x} s_{12} + L_{3 x} s_{123} + L_{5 x} s_{δ 1234} + L_{6 x} s_{δ 12345} + s)$ , $s_{δ 1234} = \sin (δ + θ_{1} + θ_{2} + θ_{3} + θ_{4})$ , and so on.

The dissipative energy of the system has the same structure as the potential energy (equation (11)) of the spring deformation, thus

2 D = b (Δ {\overset{\cdot}{y}}_{B}^{2} + \overset{\cdot}{Δ} y_{A}^{2})

(13)

Using the already determined deformations of the spring-damper elements and their first derivatives, it is possible to determine the dynamic forces in the spring-damper elements, attached to the points A and B

F_{A (B)} = k Δ y_{A (B)} + b Δ {\overset{\cdot}{y}}_{A (B)}

(14)

One can take into account that during the normal operating conditions, the computed forces must have only negative or zero values, which corresponds to the negative or zero values of the spring-damper element deformations. In case of positive values of the spring-damper element deformations, a separation of the tires from the terrain is observed and the spring-damper elements are not active anymore.

Differential equations of the system

Using the derived equations for the kinetic, potential, and dissipative energies and performing the mathematical operations in equation (2), the nonlinear dynamic equations of motion of the excavator are obtained in the form of

M (q) \overset{\cdot\cdot}{q} + V (q, \overset{\cdot}{q}) + G (q) = Q

(15)

The notations used are as follows: inertia matrix $M (q)_{6 \times 6}$ , consisting of inertia terms, it should be noted that the inertia matrix is symmetric and positive definite; vector $V (q, \overset{\cdot}{q})_{6 \times 1}$ , consisting of centrifugal and Coriolis terms; vector $G (q)_{6 \times 1}$ , consisting of terms proportional to the weight of the links; vector $Q_{6 \times 1}$ , consisting of generalized forces and moments, associated with the respective generalized coordinates: $Q_{1} = 0$ , $Q_{2} = τ_{2}$ , $Q_{3} = τ_{3}$ , $Q_{4} = τ_{4}$ , $Q_{5} = 0$ , and $Q_{6} = 0$ where $τ_{i}$ denotes the driving torques applied to the corresponding joints. The elements of the matrices are not presented here due to their rather large size.

Due to the already available inverse kinematics solution for the generalized coordinates θ₂, θ₃, and θ₄ and their corresponding velocities and accelerations (see “Kinematic model and trajectory planning” section), the system of differential equations (15) are divided into two sets.

The first set of equations consists of the three equations which describe the motion of the system for coordinates θ₁, θ₅, and s:

{\overset{\cdot\cdot}{q}}_{I} = - M {(q)}_{I}^{- 1} (V {(q, \overset{\cdot}{q})}_{I} + G {(q)}_{I})

(16)

where ${\overset{\cdot\cdot}{q}}_{I} = [\begin{matrix} {\overset{\cdot\cdot}{θ}}_{1} & {\overset{\cdot\cdot}{θ}}_{5} & \overset{\cdot\cdot}{s} \end{matrix}]^{T}$ denotes the vector of the generalized accelerations and $M (q)_{I}$ , $V (q, \overset{\cdot}{q})_{I}$ , and $G (q)_{I}$ denote the matrices, corresponding to the first set of equations. The system of equations (16) is considered as an initial-value problem and is solved by a standard numerical routine.

The second set consists of the rest of the equations used to compute the necessary driving torques τ₂, τ₃, and τ₄

Q_{II} = M {(q)}_{II} {\overset{\cdot\cdot}{q}}_{II} + V {(q, \overset{\cdot}{q})}_{II} + G {(q)}_{II}

(17)

where $Q_{II} = [\begin{matrix} τ_{2} & τ_{3} & τ_{4} \end{matrix}]^{T}$ denotes the vector of driving torques and ${\overset{\cdot\cdot}{q}}_{II} = [\begin{matrix} {\overset{\cdot\cdot}{θ}}_{2} & {\overset{\cdot\cdot}{θ}}_{3} & {\overset{\cdot\cdot}{θ}}_{4} \end{matrix}]^{T}$ denotes the vector of the known accelerations. $M (q)_{II}$ , $V (q, \overset{\cdot}{q})_{II}$ , and $G (q)_{II}$ denote the matrices, corresponding to the second set of equations.

Overturning stability of the excavator

As Alexandrov³⁴ states, the overturning stability of the mobile machine is defined as its ability to resist the overturn. The excavator performing lifting operations will overturn around the tipping axis (point B, see Figure 2) if the overturning moment is greater than the restoring moment for a sufficient period of time. The excavator overturning stability is estimated by the overturning stability coefficient

u_{s} = \frac{M_{res}}{M_{ov}}

(18)

whose value must be equal or greater than 1.15. In equation (18), M_res and M_ov denote the sums of the restoring and the overturning moments around the excavator tipping axis, respectively. When u_s = 1, the balanced condition is reached and the reaction in the left support (point A) becomes equal to 0.

The arising inertial loads during the movement of the links, combined with the payload swinging and motion of the base body due to the elastic mounting and soft ground, can contribute to the decreasing of the excavator overturning stability. That is why during the performance of lifting operations, the dynamic character of M_res and M_ov must be considered and the dynamic stability coefficient $u_{s}^{dyn}$ must be used to determine the stability of the excavator.

Numerical experiments and discussions

The developed system of differential equations (15) is used in order to perform numerical experiments for different operating conditions and study the excavator static and dynamic stability. The inertia and geometrical data are obtained by CAD modeling of a real-life excavator design. The numerical values used for the inertia, geometrical, stiffness, and damping data of the excavator elements are shown in Table 1, where the measurement units are as follows: for mass (kg), for mass moment of inertia (kg m²), for length (m), for angles (degrees), for spring constant (N/m), and for damping coefficient (Ns/m).

Table 1.

Values for the inertia and geometrical data.

Link	Value
Base body	$m_{1} = 18 \times 10^{3}$ , $J_{1} = 7.53 \times 10^{3}$ , $L_{1 x} = 0.51$ , $L_{1 y} = 0.98$ , $L_{C_{1} x} = 0$ , $L_{C_{1} y} = 0$ , ${{}^{1}x}_{B} = 1.5$ , ${{}^{1}y}_{B} = - 1.04$ , ${{}^{1}x}_{A} = - 1.1$ , ${{}^{1}y}_{A} = - 1.04$
Boom	$m_{2} = 2.52 \times 10^{3}$ , $J_{2} = 4.208 \times 10^{3}$ , $L_{2 x} = 4.5$ , $L_{C_{2} x} = 2.2$ , $L_{C_{2} y} = 1.19$
Stick	$m_{3} = 1.83 \times 10^{3}$ , $J_{3} = 1.26 \times 10^{3}$ , $L_{3 x} = 3$ , $L_{C_{3} x} = 0.93$ , $L_{C_{3} y} = 0.19$
Bucket	$m_{4} = 540$ , $J_{4} = 148$ , $L_{4 x} = 1$ , $L_{5 x} = 0.6$ , $L_{C_{4} x} = 0.6$ , $L_{C_{4} y} = 0.52$ , $δ = 70$
Payload	$m_{5} = 2.5 \times 10^{3}$ , $L_{6 x} = 1.5$
Spring and damping constants	$k = 4.61 \times 10^{6}$ , $b = 45.01 \times 10^{3}$

To study the dynamic behavior of the system during the different operating conditions, a few simulations are carried out for the following parameters of the trajectory: $y_{H}^{s} = 4 m$ , $y_{H}^{f} = - 4.5 m$ , x_H = 4 m, t_f = 12 s, and the bucket angle is kept constant at φ = 180°. The simulated motion corresponds to the payload lowering. The length of the simulations is divided into two time intervals. The first interval is [0, t_f] and it corresponds to the motion along the trajectory. The second interval [t_f, t_f + 3] is intended to capture the dynamic effects in the system after the end of the motion along the trajectory. The initial values of the coordinates θ₁(0) and s(0) correspond to the stable equilibrium position of the system, determined by the gradient of the potential energy (equation (5)) with respect to the generalized coordinates.

In the first simulation, the initial angle of the payload is set to θ₅(0) = 20° (0° from the vertical position), which corresponds to the motion along the vertical trajectory without a payload swinging. In the second simulation, the initial angle of the payload is set to θ₅(0) = 10° (10° from the vertical position), which corresponds to the presence of a payload swinging.

To clarify the dynamical behavior of the excavator, frames from the motion animation of the two considered cases are shown in Figure 4. The following notations are used: 1—mechanical structure; 2—trajectory of the gravity center CG of the whole structure; 3, 4, 5, and 6—the trajectories of the gravity centers C₂, C₃, C₄, and R of the boom, stick, bucket, and payload, respectively. It is remarkable that despite the constant distance between the payload R and the right support B during the vertical motion along the trajectory 7, the gravity center CG of the whole system moves toward the right support, but still remains between the two supports. The reasons for that are the variable boom and stick x-positions (curves 3 and 4), also the payload swinging and base body motion, thus reducing the stability reserve. This fact must be taken into account during the excavator exploitation and used to create the mobile machine load charts for the considered case.

Figure 4.

Frames of the motion animation and trajectories of the specific points: (a) in the absence of payload swinging and (b) in the presence of payload swinging.

For the case without a payload swinging, in Figure 5(a), the static coefficient of stability $u_{s}^{st}$ is presented jointly with the dynamic one $u_{s}^{dyn}$ as a function of the position $y_{H}^{d}$ along the trajectory. The static stability coefficient is calculated without considering the support deformations and by taking into account only the weights of the links, the payload, and the base body. As can be seen from the presented results, the dynamic coefficient is slightly smaller than the static one. The main reason for that is the deformation of the supports in the equilibrium position leading to the machine tilting to the right (see Figure 2) and, as a consequence, an increase of the overturning moments and decrease of the restoring moment is observed. Figure 5(b) shows the graph of the dynamic stability coefficient in the case of the presence of a payload swinging. Its values have oscillatory behavior, induced by the interaction of the payload swinging (see Figure 4(b), curve 6) and the oscillatory structure. For example, for $y_{H}^{d} = 2.5 m$ , the value of the dynamic coefficient of the overturning stability is smaller than the static one by approximately 8%. The reduction of the excavator stability due to the presence of payload swinging should be taken into account during the performance of lifting operations, especially when the coefficient of stability is close to its limit value.

Figure 5.

Static $u_{s}^{st}$ and dynamic $u_{s}^{dyn}$ stability coefficients: (a) in the absence of payload swinging and (b) in the presence of payload swinging.

To study the system dynamic behavior when the stability coefficient is close to 1, additional numerical experiments are conducted. For this purpose, the mass of the payload is chosen to be m₅ = 4300 kg and the simulation model is modified to allow the separation of the left support (see Figure 2, point A) from the ground. The total length of the simulation is set to 17 s. The first numerical experiment is performed without a payload swinging. Figure 6(a) shows the time evolution of the deformation of the spring in the left support Δy_A, computed according to equation (10). The shape of the deformation graph is influenced by the vertical motion of the payload combined with the vibrations due to the elastic mounting. The values of the deformation are negative during the simulation interval. Around 7.5 s, the deformation is very close to 0 which means that balance point is reached, but the left support is still in contact with the ground. The time evolution of the dynamic stability coefficient for the considered case is shown in Figure 6(b).

Figure 6.

Time evolution of the system characteristics (in the absence of payload swinging): (a) spring deformation and (b) dynamic stability coefficient.

The second experiment is performed in the presence of payload swinging, achieved by setting the initial value of the payload angle equal to 7°, measured from the vertical position. Figure 7(a) clearly shows the separation of the left support from the ground—the maximum value of Δy_A is 0.08 m. After the first separation, one can see three consecutive impacts between the base body and the terrain. Although the left support is separated for a certain period of time, the excavator initial position has been restored. Figure 7(b) shows the graph of the dynamic overturning stability coefficient corresponding to the case studied. Compared to Figure 6(b), the overturning stability coefficient is smaller and at t = 8.5 s, its value is less than 1.

Figure 7.

Time evolution of the system characteristics (in the presence of payload swinging): (a) spring deformation and (b) dynamic stability coefficient.

A typical real-world case is the lifting of the payload laying on the ground with loose slings, that is, the load is lifted with a non-zero initial bucket speed. The simulation results for motion with loose slings for 2 s accompanied by the payload swinging are shown in Figure 8. In this case, during the sling tightening due to the vertical motion, a sudden displacement of the system gravity center to the right support (see Figure 2, point B) is occurred, accompanied by the left support separation (see Figure 8(a)) and a sudden overturning stability decrease which is restored later (Figure 8(b)).

Figure 8.

Time evolution of the system characteristics (in the presence of payload swinging): (a) spring deformation and (b) dynamic stability coefficient.

Another operating condition causing the change in the excavator dynamic behavior is the sudden change of the payload weight due to the payload spillage or slings and hanging elements failure. In Figure 9, the results of the simulation during the payload raising and sudden change of the payload weight are shown, and in Figure 10, the same during the payload lowering is shown. As one can observe, both cases lead to oscillatory motions of the mechanical system, but without extremal values and support separation.

Figure 9.

Time evolution of the system characteristics (in the presence of payload swinging): (a) spring deformation and (b) dynamic stability coefficient.

Figure 10.

Time evolution of the system characteristics (in the presence of payload swinging): (a) spring deformation and (b) dynamic stability coefficient.

Parameter estimation and validation of the mathematical model

Tire stiffness and damping characteristics estimation

The accuracy of the simulation results strongly depends on the accuracy of the dynamic model parameter values. While the inertia and geometrical parameters are determined relatively easy by proper measurements, calculations, and CAD modeling, the stiffness and damping properties have to be determined by an experimental research. Although the excavator tires are complex structures with nonlinear characteristics that depend on a number of factors,³⁵ under certain load conditions, their stiffness and damping properties in the vertical direction can be adequately modeled as a linear spring connected in parallel combined with a viscous damper. A number of methods are used to estimate the numerical values of the stiffness and damping coefficients, among which is the well-known Free-Vibration Logarithm Decay Method. Typical usage of this method includes the application of an impulse loading to a single tire and the use of the measured free vibrations signal for estimation of the stiffness and damping coefficients values. The main advantage of this method is that it allows estimating the stiffness and damping coefficients of the tire together with the contacting surface. Due to difficulties to test the single tire in this study, a new method is presented here to carry out the study of the tire elastic and damping characteristics.

The excavator base body jointly with the fixed in a certain configuration digging manipulator without payload (see Figure 2) is represented as a single rigid body with two degrees of freedom, mounted on linear springs and viscous dampers in parallel. Two generalized coordinates represent the oscillatory motion of the excavator—vertical motion s of the gravity center C and in-plane rotation θ according to the same point (see Figure 11(a)).

Figure 11.

(a) Dynamic model of the excavator, represented as a single rigid body, and (b) experimental and fitted curves of the gravity center vertical motion acceleration.

An impulse loading is applied to the oscillatory structure by the excavator boom lowering at maximal velocity followed by a sudden stopping caused by the corresponding hydraulic valve closing. The appeared free decaying oscillations of the structure are measured by a two-axis accelerometer mounted on the chassis in the point M, shown in Figure 11(a). Under the assumption of small oscillations, the following equation is used to determine the generalized acceleration ${\overset{\cdot\cdot}{s}}_{e}$ of the gravity center (see Figure 11(b))

{\overset{\cdot\cdot}{s}}_{e} = \frac{a_{x} x_{M} + a_{y} y_{M}}{y_{M}}

(19)

whereby a_x and a_y are denoted the x_M and y_M accelerations, captured by the accelerometer.

Based on the theory of the free decaying oscillations, the motion of the body along the vertical translation coordinate is presented in the following form

\overset{\cdot\cdot}{s} (t) = \sum_{k = 1}^{2} A_{k}^{s} e^{- B_{k} t} \cos (ω_{dk} t - ψ_{k})

(20)

where ω_dk is the kth damped natural angular frequency, $A_{k}^{s}$ is the initial amplitude, B_k is the decay rate, and ψ_k is the initial phase. The theoretical solution (equation (20)) is fitted to the experimental data (equation (19)) by the use of the least-squares method and the following numerical values of the parameters are found: A₁ = 26.29, A₂ = –6.04, B₁ = 2.37, B₂ = 0.55, ω_d₁ = 21.9, ω_d₂ = 10.55, ψ₁ = –0.11, ψ₂ = –0.05. The coefficient of multiple determination is R² ≈ 0.97, which can be explained by the noise of the measured signal. Figure 11(b) depicts the experimental and the fitted curves. The obtained numerical values are used to compose the roots of the characteristic polynomial

λ_{1, 2}^{k} = - B_{k} \pm i ω_{dk}

(21)

and to calculate the characteristic polynomial f_e(λ) coefficients $a_{j}^{e}$

f_{e} (λ) = \sum_{j = 0}^{4} a_{j}^{e} λ^{j} = Π_{k = 1}^{2} (λ - λ_{1}^{k}) (λ - λ_{2}^{k})

(22)

The system of differential equations, describing the free decaying oscillations of the presented two degrees of freedom system, is

M \overset{\cdot\cdot}{q} + B \overset{\cdot}{q} + Cq = 0

(23)

where $M = [\begin{matrix} m & 0 \\ 0 & J \end{matrix}]$ , $C = [\begin{matrix} 2 k & k (L_{B} - L_{A}) \\ k (L_{B} - L_{A}) & k (L_{A}^{2} + L_{B}^{2}) \end{matrix}]$ , $B = [\begin{matrix} 2 b & b (L_{B} - L_{A}) \\ b (L_{B} - L_{A}) & b (L_{A}^{2} + L_{B}^{2}) \end{matrix}]$ , m is the mass of the body, J is the mass moment of inertia of the body according to the point C, k is the tire stiffness coefficient, and b is the damping coefficient.

Comparing the experimentally determined coefficients $a_{j}^{e}$ to the derived from the theoretical model (equation (23)) characteristic polynomial coefficients $a_{j}^{t}$ , it is possible to estimate the numerical values of the stiffness and damping coefficients.

Inserting notations $L_{1} = L_{B} - L_{A}$ and $L_{2} = L_{A}^{2} + L_{B}^{2}$ , the theoretically derived coefficients $a_{j}^{t}$ for the system (equation (23)) are as follows: $a_{0}^{t} = k^{2} (2 L_{2} - L_{1}^{2}) / Jm$ , $a_{1}^{t} = 2 bk (2 L_{2} - L_{1}^{2}) / Jm$ , $a_{2}^{t} = [k (m L_{2} + 2 J) - b^{2} (L_{1}^{2} - 2 L_{2})] / J m$ , $a_{3}^{t} = b (2 J + L_{2} m) / Jm$ , and $a_{4}^{t} = 1$ . Thus, comparing $a_{0}^{t}$ to $a_{0}^{e}$ , as well as $a_{3}^{t}$ to $a_{3}^{e}$ , the following equations for calculation of k and b are obtained

k = \sqrt{\frac{a_{0}^{e} Jm}{2 L_{2} - L_{1}^{2}}}

(24)

b = \frac{a_{3}^{e} Jm}{2 J + L_{2} m}

(25)

For m = 22,890 kg, J = 154,320 kg m², L_A = 2.5 m, and L_B=0.5 m, the calculated values are k=4610 kN/m and b = 45,098 Ns/m.

Validation of the mathematical model

The mathematical model is validated by the comparison between the simulation and the experimentally obtained data.³⁶ For this purpose, a small-scale elastically mounted manipulator, depicted in Figure 12, with known geometrical and inertial parameters and having the same kinematic structure as the considered excavator is used. The manipulator consists of an arm and forearm and allows free suspension of the load to the forearm distal end. It consists of the following main elements: a rotating column 1, an arm 2, a forearm 3, and a hydraulic power unit 4. The structure is equipped with sensors for measuring six different quantities as follows: (1) the pressures at the cap end and the rod end of the arm 5 and forearm 6 driving cylinders 5 and (2) the displacement of the arm and the forearm driving cylinders 7 and 8. For the displacement measurement, two Linear Variable Inductance Transducers are used, attached to the hydraulic cylinders, while for the pressure measurement, four pressure sensors to measure absolute pressure are utilized. Analog outputs from the sensors are measured, converted, and transmitted to the computer measurement system 9 by an analog-to-digital converter.

Figure 12.

View of the small-scale manipulator and sensor placement.

The lifting process starts with a fixed in a certain forearm, payload laying on the ground and loose slings. The process corresponds to the following sequence of motions: (1) arm raising—from initial angle θ₁ = –12° to angle θ₁ = 20.5°; (2) dwell phase in the upper position; and (3) arm lowering to the initial position.

The experimental data for the measured quantities are shown in Figure 13(a): (1) the pressure at the cap end of the arm, (2) the pressure at the rod end of the arm, (3) the pressure at the cap end of the forearm, (4) the pressure at the rod end of the forearm, (5) linear displacement of the arm cylinder, and (6) the linear displacement of the forearm cylinder. As can be expected, the payload swinging inserts an additional low-frequency variation in the pressures—see Figure 13(a), curves 1–4. The following phases with time length t_i can clearly be distinguished in the displacement of the arm cylinder (see Figure 13(a), line 5): 1—arm raising without payload (loose slings), 2—arm raising with payload (tight slings), 3—fixed arm (dwell phase), and 4—arm lowering.

Figure 13.

(a) The experimental and (b) the comparison between the experimental and the simulation data.

The measured and filtered linear displacement signal of the arm hydraulic cylinder (Figure 13(a), line 5) is used to perform a simulation by the use of the theoretical mathematical model (equation (15)) with the corresponding known numerical values of the small-scale manipulator parameters. In Figure 13(b), the obtained by the simulation torques of the arm and forearm are used to determine the driving hydraulic cylinders forces (arm—line 1, forearm—line 3) shown together with the experimental forces (arm—line 2, forearm—line 4), determined by the use of the measured hydraulic cylinders pressures, presented in Figure 13(a), lines 1–4. Although some differences are observed during the duty cycle, the good correlation between the experimental data and simulation results allows us to conclude that the adopted mathematical model may confidently be used to simulate the dynamic problems presented.

Conclusion

In this article, we have studied the dynamic stability of a hydraulic excavator during performing lifting operations. The developed dynamic model with six degrees of freedom considers the base body elastic connection with the terrain, the front digging manipulator links, and the presence of the freely suspended payload swinging. A system of non-linear differential equations describing the dynamic behavior is obtained by using the Lagrange formalism. Numerical experiments are carried out to study the excavator dynamic overturning stability during the motion along a vertical straight-line vertical trajectory under various operating conditions. Finally, the mechanical system model has been validated by a small-scale experimental model.

The revealed basic insights allow a deeper understanding of the dynamic behavior of the excavator during performing lifting operations along the vertical trajectory. It is shown that the arising inertial loads due to the movement of the links, combined with the payload swinging, motion of the system gravity center, and the motion of the base body, decrease the excavator overturning stability. It was found that the excavator overturning stability while following a vertical straight-line trajectory decreases during the motion from the higher to the lower part of the trajectory. If the stability coefficient is close to 1, the payload swinging can cause the separation of a support from the terrain; nevertheless, it was found that the excavator overturning stability can be restored. The low-frequency vibration, induced in the hydraulic cylinders by the payload swinging and the base body elastic mounting, should be taken into account during the hydraulic system design and its improvement. The proposed calculation method for estimation of the elastic and damping characteristics of the supports is based on the experimental data fitting and can be used reliably.

The decrease of the overturning stability coefficient due to the payload swinging and moving gravity center of the system is very important from the exploitation point of view and this knowledge should be used for the following:

Correction of the existing and creation of new mobile machine load diagrams in order to keep the overturning stability reserve demanded by the standards;

Education of the excavator operators to account for the specific aspects of the excavator dynamical behavior during performing lifting operations.

Footnotes

Handling Editor: James Baldwin

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: We acknowledge support by the German Research Foundation and the Open Access Publication Fund of TU Berlin.

ORCID iD

Dragan Marinković

References

ISO 10567:2007. Earth—moving machinery—hydraulic excavators—lift capacity.

National Institute for Occupational Safety and Health (NIOSH). Preventing injuries when working with hydraulic excavators and backhoe loaders. Washington, DC: NIOSH, https://www.cdc.gov/niosh/docs/wp-solutions/2004-107/ (2004, accessed 13 September 2017).

Edwards

Holt

. Case study analysis of risk from using excavators as “cranes.” Automat Constr 2010; 19: 127–133.

Edwards

Holt

DG.

Case study analysis of construction excavator H&S overturn incidents. Eng Construct Arch Manage 2010; 17: 493–511.

Lim

Kim

Hwan

et al . Development of tipping-over rate computation system for hydraulic excavator having crane function. In: Proceedings of the 8th Russian-Korean International Symposium on Science and Technology (KORUS), Tomsk, 26 June–03 July 2004. IEEE.

BS EN 474—1:2006. Earth—moving machinery—safety—part 1: general requirements.

BS EN 474—5:2006. Earth—moving machinery—safety—part 5: requirements for hydraulic excavators.

Off-Highway Plant and Equipment Research Centre (OPERC). Safety alert (ALT—017)—excavators used as cranes. Dudley: OPERC, https://www.operc.com/safetyalert/alt_017_a4.pdf (2017, accessed 9 September 2017).

Park

Lee

et al . Static compensation ZMP algorithm preventing tips-over of a tele-operation excavator. In: Proceedings of the 29th international symposium on automation and robotics in construction, Eindhoven, 26–29 June 2012. Technische Universiteit Eindhoven.

10.

Heep

Schneider

Walawalkar

et al . Method for payload estimation in excavators. In: Proceedings of the 4th commercial vehicle technology symposium CVT 2016, Kaiserslautern, 8–10 March 2016. Aachen: Shaker.

11.

Gang

Hang

Yang

et al . The stability analysis of luffing jib tower crane with different working conditions. In: Proceedings of the 2nd international conference on civil engineering and rock engineering (ICCERE 2017), Guangzhou, China, 2–3 December 2017.

12.

Shaikh

Dineesh

KD.

Lifting capacity enhancement of a crawler crane by improving stability. J Theor Appl Mech 2016; 54: 219–227.

13.

Rauch

Singhose

Fujioka

et al . Tip-over stability analysis of mobile boom cranes with swinging payloads. J Dyn Syst-T ASME 2013; 35: 0310081.

14.

Messioud

Sbartai

Dias

Soil by taking into account the soil-foundation interaction. Struct Eng Mech 2016; 58: 887–903.

15.

Towarek

Dynamic stability of a crane standing on soil during the rotation of the boom. Int J Mech Sci 1998; 40: 557–574.

16.

Abo-Shanab

Sepehri

. On dynamic stability of manipulators on moveable platforms. In: Proceedings of the 2001 IEEE International symposium on computational intelligence in robotics and automation, Banff, Alberta, Canada, 29 July–01 August 2001. IEEE.

17.

Janosevic

Pavlovic

Jovanovic

et al . A numerical and experimental analysis of the dynamic stability of hydraulic excavators. Facta Univ Ser: Mech Eng 2018; 16: 157–170.

18.

Mitrev

Janošević

Marinković

Dynamical modelling of hydraulic excavator considered as a multibody system. Teh Vjesn 2017; 24(Suppl. 2): 327–338.

19.

Abdel-Rahman

Nayfeh

Masoud

ZN.

Dynamics and control of cranes: a review. Modal Anal 2003; 9: 863–908.

20.

Fodor

Vázquez

Freidovich

et al . Towards oscillation reduction in forestry cranes. In: Proceedings of the symposium on fluid power and motion control BATH/ASME 2016, Bath, 7–9 September 2016.

21.

Alkamachi

Ercelebi

Modelling and control of H-shaped racing quadcopter with tilting propellers. Facta Univ Ser: Mech Eng 2017; 15: 201–215.

22.

Mitrev

Savov

A theoretical-experimental approach for elasto-damping parameters estimation of cone inertial crusher mounting. Facta Univ Ser: Mech Eng 2017; 15: 73–83.

23.

Edalati

Bayat

Pakar

et al . A novel approximate solution for nonlinear problems of vibratory systems. Struct Eng Mech 2016; 57: 1039–1049.

24.

Gabbert

Duvigneau

Ringwelski

Noise control of vehicle drive system. Facta Univ Ser: Mech Eng 2017; 15: 183–200.

25.

Yoon

HS.

A review on mechanical and hydraulic system modeling of excavator manipulator system. J Construct Eng 2016; 2016: 1–11.

26.

Vujic

Lazarevic

Batinic

Development of dynamic-mathematical model of hydraulic excavator. J Central South Univ 2017; 24: 2010–2018.

27.

Mitrev

Tudjarov

Todorov

Cloud based expert system for synthesis and evolutionary optimization of planar linkages. Facta Univ Ser: Mech Eng 2018; 16: 139–155.

28.

Caterpillar, http://www.caterpillar.com (2017, accessed 13 September 2017).

29.

Zhang

Wang

Liu

et al . Research on trajectory planning and autodig of hydraulic excavator. Math Probl Eng 2017; 2017: 1–10.

30.

Taylor

Seward

Proportional-integral-plus control of an intelligent excavator. Comput-Aided Civ Inf 2004; 19: 16–27.

31.

Le Hanh

Ahn

Kha

et al . Trajectory control of electro-hydraulic excavator using fuzzy self tuning algorithm with neural network. J Mech Sci Technol 2009; 23: 149–160.

32.

Siciliano

Sciavicco

Villanni

et al . Robotics: modelling, planning and control. London: Springer, 2009, pp.91–94.

33.

Merkin

DR.

Introduction to the theory of stability. New York: Springer, 1996, pp.75–76.

34.

Alexandrov

MP.

Lifting machines. Moscow: Mashinostroenie, 1986, pp.234–236.

35.

Taylor

Bashford

Schrock

Methods for measuring vertical tire stiffness. T ASAE 2000; 43: 1415–1419.

36.

Petrović

Maneski

Trišović

et al . Identification of crack initiation cause in pylons construction of the excavator SchRs630. Teh Vjesn 2018; 25: 486–491.

Numerical study of the hydraulic excavator overturning stability during performing lifting operations

Abstract

Keywords

Introduction

Kinematic model and trajectory planning

Dynamic model

Kinetic energy of the system

Potential and dissipative energies of the system

Differential equations of the system

Overturning stability of the excavator

Numerical experiments and discussions

Parameter estimation and validation of the mathematical model

Tire stiffness and damping characteristics estimation

Validation of the mathematical model

Conclusion

Footnotes

Declaration of conflicting interests

Funding

ORCID iD

References