Dynamic Modeling and Analysis of Loader Working Mechanism Considering Cooperative Motion with the Vehicle Body

Abstract

Achieving precise load detection for Intelligent Loaders is an important task, which directly affects the operation energy efficiency and the fatigue life analysis for the loader’s working mechanism. The operation of the mechanism is regarded as a 3-DOF (degree of freedom) planar motion process coordinated with the vehicle body. Affected by complex dynamic coupling in motion, the existing dynamic models of the mechanism have the problem of insufficient accuracy, which is not conducive to the precise calculation of load. Taking the reverse six-linkage loader as the research object, an accurate dynamic model of the mechanism is established considering its cooperative motion with the vehicle body. Firstly, the kinematic description of the mechanism is given by the Rodriguez method. Then, to overcome the coupling effect caused by the cooperative motion, the sufficient inertia forces of the mechanism are established in joint space using the Lagrange method. Furthermore, to overcome the coupling effect caused by the complex structure, the Newton–Euler method is used to establish the force mapping relations between the joint space and the drive space by multi-body modeling. Finally, the dynamic model of the mechanism in drive space is obtained, and the specific mapping relations between the bucket force, the vehicle driving force, and the drive parameters are given. Compared with existing dynamic models in simulation, the analysis shows that the average and maximum absolute errors of the vehicle driving force calculated by the established model do not exceed 20% of the existing model errors, and the corresponding errors of the bucket force do not exceed 10% of the existing model errors, which proves that the motions of vehicle body and front-end mechanism, as well as the force of the tilt hydraulic cylinder, play important roles in improving the model accuracy. The established model is superior to existing models and is more suitable for cooperative motion with the vehicle body.

1. Introduction

Loaders are important construction machines whose primary task is to load and transport bulk materials such as rock and soil. They are widely used in different applications, such as mines, ports, and construction sites. The operation cycle of the loader is mainly divided into four stages: shoveling, transportation, unloading, and flattening, which are mainly completed by the loader’s working mechanism [1,2,3]. The operating load directly affects the performance of the working mechanism, so achieving precise load detection has been an important task realted to loaders [1,4].

The front-end mechanism of the loader is mainly a 2-DOF (degree of freedom) mechanism driven by a lift hydraulic cylinder and a tilt hydraulic cylinder. In practical application, the front-end mechanism and the vehicle body work together, that is, the motion of the front-end mechanism is usually accompanied by the forward or backward motion of the vehicle body. So, the operation of the loader’s working mechanism can be considered as a 3-DOF planar motion process driven by the vehicle body and the hydraulic cylinders [5]. The vehicle’s motion will bring additional inertia forces to the front-end mechanism, and the inertia forces of the front-end mechanism will affect the precise measurement of the operating load [6], etc. Therefore, building an accurate dynamic model of the working mechanism in the drive space has important theoretical and practical value for the research of the precise online detection for the load, the driving force distribution for the vehicle’s motion, etc.

A series of achievements have been made in the dynamic modeling of the loader’s working mechanism. The existing dynamic models can be roughly divided into the following three categories:

(1) Without considering the influence of the vehicle motion, the working mechanism is simplified as a dynamic model of the two-bar mechanism. Some scholars [5,7,8,9,10] have established dynamic models of the working mechanism, which was simplified as a two-bar mechanism when the vehicle body was static, for the research purpose of shoveling trajectory planning and control, load weight measurement, structural parameter optimization, etc. In the representative research [8], Kang et al. established a dynamic model of the reverse six-linkage mechanism, which was simplified as a two-bar mechanism (only boom and bucket) on the assumption that the load was balanced and the vehicle body was static, for measuring the load weight. Such models only considered the dynamic characteristics of the two bars (boom and bucket) in the joint space without considering the dynamic effects caused by the vehicle motion and neglected the tilt cylinder action (affected by the offset load) when converting to the drive space, resulting in a loss of accuracy.

(2) Without considering the influence of the vehicle motion, a simplified dynamic or static model of the working mechanism with dynamic compensation can be created. Some scholars [6,11,12,13,14,15,16,17,18,19,20,21,22,23] have established simplified dynamic or static models of the reverse six-linkage, TP-linkage, or agricultural working mechanisms that include dynamic compensation for some components or the offset load in the drive space. Such models are used in the research of shoveling trajectory planning, resistance analysis, load measurement, and structure optimization, etc. In the representative research [22], Madau et al. established a static model of the reverse six-linkage working mechanism with the dynamic compensation of the boom and then proposed an online estimation algorithm to calculate the load force. Such models only considered the dynamic characteristics of some components or the offset load in the drive space without considering the dynamic effects caused by the vehicle motion, resulting in insufficient accuracy in the cooperative motion with the vehicle body.

(3) Considering the influence of the vehicle motion, a static model of the working mechanism can be created. Some scholars [2,24,25] have established static models of the reverse six-linkage or TP-linkage working mechanisms in the drive space to study vehicle control, torque distribution, inertia parameter estimation, etc. In the representative research [2], Gao et al. established a static model of the reverse six-linkage working mechanism considering the load force and proposed a torque distribution strategy based on this model for vehicle driving control. Such models mainly calculated the driving force of the vehicle body to achieve the driving control for the loader without considering the dynamic effects caused by the motion of the front-end mechanism, resulting in insufficient accuracy in the cooperative motion with the front-end mechanism.

In a word, the above models have the problem of insufficient accuracy in the cooperative motion of the vehicle body and the front-end mechanism (i.e., the working mechanism works in a 3-DOF motion), which is not conducive to the load calculation and directly affects the operation energy efficiency and the fatigue life analysis of the mechanism system [4,26]. Therefore, achieving precise load detection puts forward higher accuracy requirements for the dynamic model of the working mechanism during operation.

To improve the accuracy of the dynamic model of the loader’s working mechanism during operation, an accurate dynamic model of the working mechanism is established in the 3-DOF motion, which fully considers the inertia forces and the coupling of the mechanism system, including the influence of the vehicle’s motion. Specifically, taking the reverse six-linkage loader as the research object, first, the kinematic description of the mechanism system in the joint space is given by the Rodriguez method, and then the kinematic model in the drive space is given by the motion mapping relations, which can provide basic motion state information for the dynamic modeling of the mechanism system. Then, drawing on the idea of analytical mechanics, considering the strong dynamic coupling characteristics of the mechanism system caused by the cooperative motion of the vehicle body and the front-end mechanism comprehensively, the displacement of the vehicle body is introduced as an independent generalized coordinate; thus, the inertia and generalized forces of the mechanism system, including the inertia forces caused by the vehicle motion, are fully established by the Lagrange method, and the dynamic model of the mechanism system is finally obtained in the joint space. Furthermore, drawing on the idea of vector mechanics and considering the structural coupling of the mechanism and the action of the tilt hydraulic cylinder mainly affected by the offset load, the Newton–Euler method is used to model the whole body and the single body of the mechanism system to find the force mapping relations between the joint space and the drive space. Finally, the dynamic model of the mechanism system is obtained in the drive space, and the specific mapping relations between the vehicle driving force, the bucket force and the drive space parameters are given. In the validation and analysis section, the simulation platform is first built, and then the kinematic model and the dynamic model are verified separately. In the dynamic model’s validation, by comparing it with the existing dynamic models, the analysis shows that the dynamic model built in this paper has higher accuracy and applies to the cooperative motion of the vehicle body and the front-end mechanism.

2. Drive Space and Parameters

In Figure 1, the working mechanism system of the reverse six-linkage loader can be considered to consist of five rigid bodies and two driving devices interconnected by revolute joints. Among them, bodies 1∼5 represent the vehicle body, boom, bucket, rotating bar, and connecting bar, respectively. Driving device I represents the lift hydraulic cylinder used to drive the boom up and down. Driving device II represents the tilt hydraulic cylinder used to drive the bucket’s rotation. The working surface (ground) of the loader is generally considered as a plane (not a level surface) during the shoveling operation of the working mechanism, and in this working condition, the motion of the vehicle body (body 1) is considered as a straight-line translation. As the mechanism system is 3-DOF when working, in practical application, the operators can realize the loader’s operation by driving the vehicle body and the hydraulic cylinder devices. In addition, the vehicle body’s velocity and the hydraulic cylinders’ lengths are easy and economical to observe. Therefore, the kinematic and dynamic models of the loader’s working mechanism built in the drive space have practical significance and value. The drive space can be described as follows:

$${\mathbb{R}}_{drive}^{5\times 1}={\left[s_v\left(t\right),l_{lift}\left(t\right),l_{tilt}\left(t\right),F_{lift}\left(t\right),F_{tilt}\left(t\right)\right]}^T$$

(1)

where t is the time variable, and the displacement $s_v\left(t\right)$ of the vehicle body, the length $l_{lift}\left(t\right)$ of the lift cylinder, the length $l_{tilt}\left(t\right)$ of the tilt cylinder, the piston force $F_{lift}\left(t\right)$ of the lift cylinder, and the piston force $F_{tilt}\left(t\right)$ of the tilt cylinder are considered as the drive space parameters. Here, the three independent variables, $s_v\left(t\right)$, $l_{lift}\left(t\right)$, and $l_{tilt}\left(t\right)$, are taken as the generalized coordinates in the drive space.

Due to the complex structural coupling of the working mechanism system, it is difficult to directly establish the kinematic and dynamic models in the drive space. Therefore, the kinematic and dynamic models can be established in the joint space first, and then the models can be established in the drive space by seeking the mapping relations between the joint space and the drive space.

(For comprehensive nomenclature of all the mathematical symbols used in the paper, see Abbreviations at the end of the text.)

3. Kinematic Description of Working Mechanism System

In this section, firstly, a kinematic model of the mechanism system considering vehicle motion is established in joint space, and then motion mapping relations between the joint space and drive space are given. Finally, the kinematic model of the mechanism system is obtained in drive space, which provides a basis for accurate dynamic modeling of the mechanism system.

3.1. Kinematic Description in Joint Space

Based on the previous assumptions, in Figure 2, since the boom (body 2) and the bucket (body 3) of the loader are the most important working parts in the mechanism’s operation and the bucket tip $P_{tip}$ is the key point, to facilitate the following dynamic modeling, the loader’s reverse six-linkage working mechanism system can be modeled as two parallel kinematic chains 1 and 2. Kinematic chain 1 is $O_0$-$O_1$-$O_2$-$O_3$-$P_{tip}$, and kinematic chain 2 is $O_0$-$O_1$-$O_2$-$O_4$-$O_5$-$O_6$. Without loss of generality, for each kinematic chain, the following kinematic descriptions can be given.

The motion of the loader with respect to the fixed reference frame (ground) $O_0{x}_0{y}_0{z}_0$ is described by the joint coordinates $s_v$ and ${\theta }_i$$\left(i=1,\cdots ,n\right)$. The coordinate $s_v$ (seen as a moving joint) represents the displacement of body 1. The coordinate ${\theta }_i$ represents the relative rotation of the body i with respect to the body i-1 carried out about the joint axis determined by the unit vector ${\mathbf{e}}_i$. This vector is fixed to the body i-1. The point $O_i$ is placed on the direction ${\mathbf{e}}_i$ and represents the center point of the i-th joint connecting the bodies i-1 and i. The local coordinate frames $O_i{x}_i{y}_i{z}_i$$\left(i=1,\cdots ,n\right)$ that are fixed to the bodies $\left(i=1,\cdots ,n\right)$, respectively, are introduced (see Figure 2). As the local coordinate frame $O_1{x}_1{y}_1{z}_1$ is considered a moving joint with respect to $O_0{x}_0{y}_0{z}_0$, here ${\theta }_1\equiv 0$, ${}^i{\mathbf{s}}_v={\left[s_v,0,0\right]}^T$$\left(i=0\right)$. Note that the left superscript i indicates that the components of the corresponding vectors and matrices are given in the $O_i{x}_i{y}_i{z}_i$ local frame. The non-bold symbol corresponding to a vector in the following text represents the amplitude of the vector. Among them, ${\mathbf{l}}_{\overrightarrow{ij}}$ represents the length vector from point i to point j, and ${\gamma }_{ijk}$ represents the angle $\angle ijk$.

Without loss of generality, it is assumed that the configuration ${\theta }_i=0$$\left(i=1,\cdots ,n\right)$ is a reference configuration of the loader and that, in this configuration, the axes of all local reference frames are parallel to the corresponding axes of the fixed reference frame, that is, $x_i\Vert x_0$, $y_i\Vert y_0$ and $z_i\Vert z_0$. According to this, the transformation matrix ${\mathbf{R}}_{i,j}$$\left(i=0,\cdots ,n;j=1,\cdots ,n\right)$ from $O_j{x}_j{y}_j{z}_j$ to $O_i{x}_i{y}_i{z}_i$ reference frames have the form [27,28]:

$${\mathbf{R}}_{i,j}=\prod _{k=i+1}^j{\mathbf{R}}_k^r=\prod _{k=i+1}^j\left[\mathbf{I}+\left(1-cos{\theta }_k\right){\left({}^k{\tilde{\mathbf{e}}}_k\right)}^2+{}^k{\tilde{\mathbf{e}}}_{k}sin{\theta }_k\right],i<j$$

(2)

where ${\mathbf{R}}_k^r\in {\mathbb{R}}^{3\times 3}$ is the Rodriguez matrix [27], $\mathbf{I}\in {\mathbb{R}}^{3\times 3}$ is the identity matrix, and ${}^k{\tilde{\mathbf{e}}}_k\in {\mathbb{R}}^{3\times 3}$ is the skew-symmetric matrix [27,28] associated with the vector ${}^k{\mathbf{e}}_k$. In addition, since the Rodriguez matrix is orthogonal [27], the following holds:

$${\mathbf{R}}_{i,j}^{-1}={\mathbf{R}}_{i,j}^T\Rightarrow {\mathbf{R}}_{j,i}={\mathbf{R}}_{i,j}^T$$

(3)

where the right superscript T denotes matrix transposition.

Now, in regard to the literature [28,29], the angular velocity ${\omega }_i$ and the angular acceleration ${\alpha }_i$ of the body i can be calculated recursively, respectively, as follows:

$${}^i{\omega }_i={\mathbf{R}}_{i-1,i}^T{}^{i-1}{\omega }_{i-1}+{\dot{\theta }}_i{}^i{\mathbf{e}}_i,\left(i=1,\cdots ,n\right)$$

(4)

$${}^i{\alpha }_i={\mathbf{R}}_{i-1,i}^T{}^{i-1}{\alpha }_{i-1}+{\ddot{\theta }}_i{}^i{\mathbf{e}}_i+{\dot{\theta }}_i{\mathbf{R}}_{i-1,i}^T{}^{i-1}{\tilde{\omega }}_{i-1}{}^{i-1}{\mathbf{e}}_i,\left(i=1,\cdots ,n\right)$$

(5)

Based on the previous relations, the velocity ${\mathbf{v}}_{P_i}$ and the acceleration ${\mathbf{a}}_{P_i}$ of the any point $P_i$ of the body i, respectively, are given by the following recurrence relations (for proof, see ref. [29]):

$$\left\{\begin{array}{c}{}^i{\mathbf{v}}_{P_i}={\mathbf{R}}_{i-1,i}^T\left({}^{i-1}{\dot{\mathbf{s}}}_v+{}^{i-1}{\mathbf{v}}_{P_{i-1}}\right),\left(i=1\right)\hfill \\ {}^i{\mathbf{v}}_{P_i}={\mathbf{R}}_{i-1,i}^T\left({}^{i-1}{\mathbf{v}}_{P_{i-1}}+{}^{i-1}{\tilde{\omega }}_{i-1}\left({}^{i-1}{\mathbf{l}}_{i-1}-{}^{i-1}{\mathbf{l}}_{P_{i-1}}\right)\right)+{}^i{\tilde{\omega }}_i{}^i{\mathbf{l}}_{P_i},\left(i=2,\cdots ,n\right)\hfill \end{array}\right.$$

(6)

$$\left\{\begin{array}{ll}{}^i{\mathbf{a}}_{P_i}& ={\mathbf{R}}_{i-1,i}^T\left({}^{i-1}{\ddot{\mathbf{s}}}_v+{}^{i-1}{\mathbf{a}}_{P_{i-1}}\right),\left(i=1\right)\hfill \\ \hfill {}^i{\mathbf{a}}_{P_i}& ={\mathbf{R}}_{i-1,i}^T\left({}^{i-1}{\mathbf{a}}_{P_{i-1}}+{}^{i-1}{\tilde{\alpha }}_{i-1}\left({}^{i-1}{\mathbf{l}}_{i-1}-{}^{i-1}{\mathbf{l}}_{P_{i-1}}\right)\right)+{}^i{\tilde{\alpha }}_i{}^i{\mathbf{l}}_{P_i}+{\left({}^i{\tilde{\omega }}_i\right)}^2{}^i{\mathbf{l}}_{P_i}\hfill \\ & +{\left({}^{i-1}{\tilde{\omega }}_{i-1}\right)}^2\left({}^{i-1}{\mathbf{l}}_{i-1}-{}^{i-1}{\mathbf{l}}_{P_{i-1}}\right),\left(i=2,\cdots ,n\right)\hfill \end{array}\right.$$

(7)

where ${\mathbf{l}}_i={\mathbf{l}}_{\overrightarrow{O_i{O}_{i+1}}}$, ${\mathbf{l}}_{P_i}={\mathbf{l}}_{\overrightarrow{O_i{P}_i}}$.

In Figure 2, it is obvious that the vectors ${\mathbf{e}}_i$$\left(i=1,\cdots ,5\right)$ have the following components:

$${}^i{\mathbf{e}}_i={}^{i-1}{\mathbf{e}}_i=\left\{\begin{array}{c}{\left[0,0,1\right]}^T,i=1,3,5\hfill \\ {\left[0,0,-1\right]}^T,i=2,4\hfill \end{array}\right.$$

(8)

Note that for body 0, ${}^0{\omega }_0={\left[0,0,0\right]}^T$, ${}^0{\alpha }_0={\left[0,0,0\right]}^T$, ${}^0{\mathbf{v}}_{P_0}={\left[0,0,0\right]}^T$, and ${}^0{\mathbf{a}}_{P_0}={\left[0,0,0\right]}^T$ hold.

Therefore, for kinematic chain 1 $\left(i=1,2,3\right)$, the motion data of any point $P_i$ of the body i can be expressed by joint coordinates $\left(s_v,{\theta }_2,{\theta }_3\right)$. Similarly, for kinematic chain 2 $\left(i=1,2^{\prime },4,5\right)$, the motion data of any point $P_i$ of the body i can be expressed by joint coordinates $\left(s_v,{\theta }_{2^{\prime }},{\theta }_4,{\theta }_5\right)$. In fact, the front-end mechanism of the loader has 2-DOF. Considering the importance of kinematic chain 1 and facilitating the dynamic modeling, any joint angle can be expressed by selecting $\left({\theta }_2,{\theta }_3\right)$ as the generalized coordinates in the joint space. Specifically, ${\theta }_{2^{\prime }}={\theta }_2+{\gamma }_{O_3{O}_2{O}_4}$, ${\theta }_4$ and ${\theta }_5$ can be expressed by ${\theta }_3$, either by the methods in the next section or in the literature [30]. For this, there is the motion mapping relation in the following form:

$$\left\{\begin{array}{c}{\mathsf{\Theta }}_{joint}={\mathsf{\Phi }}_{k1}\left(s_v,{\theta }_2,{\theta }_3\right)\hfill \\ {\mathsf{\Theta }}_{joint}={\left[s_v,{\theta }_2,{\theta }_{2^{\prime }},{\theta }_3,{\theta }_4,{\theta }_5\right]}^T\hfill \end{array}\right.$$

(9)

In this section, the kinematic model of the mechanism is established in joint space, and the mapping relations between the motion data of any point $P_i$ and the joint coordinates $\left(s_v,{\theta }_2,{\theta }_3\right)$ are given.

3.2. Kinematic Description in Drive Space

Based on the explicit derivation method of the analytic geometry, this section establishes the motion mapping relations between joint space parameters $\left({\theta }_2,{\theta }_3\right)$ and drive space parameters $\left(l_{lift},l_{tilt}\right)$, so as to establish the kinematic model of the mechanism in the drive space. In Figure 2, when the design dimensions of each body are known, the motion mapping relations between the joint space and the drive space are derived as follows.

Mapping functions of the joint space parameter ${\theta }_2$:

$$\left\{\begin{array}{c}{\gamma }_{A_0{O}_2{A}_1}=co{s}^{-1}\left(\left({l_{O_2{A}_0}}^2+{l_{O_2{A}_1}}^2-{l_{lift}}^2\right)/\left(2l_{O_2{A}_0}{l}_{O_2{A}_1}\right)\right)\hfill \\ {\dot{\gamma }}_{A_0{O}_2{A}_1}=\frac{d{\gamma }_{A_0{O}_2{A}_1}}{d{l}_{lift}}{\dot{l}}_{lift}\hfill \\ {\ddot{\gamma }}_{A_0{O}_2{A}_1}=\frac{d^2{\gamma }_{A_0{O}_2{A}_1}}{d^2{l}_{lfit}}{{\dot{l}}_{lift}}^2+\frac{d{\gamma }_{A_0{O}_2{A}_1}}{d{l}_{lift}}{\ddot{l}}_{lift}\hfill \end{array}\right.$$

(10)

$$\left\{\begin{array}{c}{\theta }_2=\pi /2-{\gamma }_{A_0{O}_2{y}_1}+{\gamma }_{A_0{O}_2{A}_1}-{\gamma }_{O_4{O}_2{A}_1}-{\gamma }_{O_3{O}_2{O}_4}\hfill \\ {\dot{\theta }}_2=\frac{d{\theta }_2}{d{\gamma }_{A_0{O}_2{A}_1}}{\dot{\gamma }}_{A_0{O}_2{A}_1}\hfill \\ {\ddot{\theta }}_2=\frac{d^2{\theta }_2}{d^2{\gamma }_{A_0{O}_2{A}_1}}{{\dot{\gamma }}_{A_0{O}_2{A}_1}}^2+\frac{d{\theta }_2}{d{\gamma }_{A_0{O}_2{A}_1}}{\ddot{\gamma }}_{A_0{O}_2{A}_1}\hfill \end{array}\right.$$

(11)

Mapping functions of the joint space parameter ${\theta }_3$:

Since the velocity and acceleration of each transfer variable are derived in the same way, only the position mapping functions are given below:

$$\left\{\begin{array}{c}{\gamma }_{O_4{O}_2{B}_0}={\theta }_2+{\gamma }_{O_3{O}_2{O}_4}+{\gamma }_{B_0{O}_2{y}_1}-\pi /2\hfill \\ l_{O_4{B}_0}={\left({l_{O_2{O}_4}}^2+{l_{O_2{B}_0}}^2-2l_{O_2{O}_4}{l}_{O_2{B}_0}cos{\gamma }_{O_4{O}_2{B}_0}\right)}^{1/2}\hfill \\ {\gamma }_{O_2{O}_4{B}_0}=co{s}^{-1}\left(\left({l_{O_2{O}_4}}^2+{l_{O_4{B}_0}}^2-{l_{O_2{B}_0}}^2\right)/\left(2l_{O_2{O}_4}{l}_{O_4{B}_0}\right)\right)\hfill \\ {\gamma }_{B_0{O}_4{B}_1}=co{s}^{-1}\left(\left({l_{O_4{B}_0}}^2+{l_{O_4{B}_1}}^2-{l_{tilt}}^2\right)/\left(2l_{O_4{B}_0}{l_{O_{4}B}}_{{}_1}\right)\right)\hfill \\ {\gamma }_{O_2{O}_4{B}_1}={\gamma }_{B_0{O}_4{B}_1}-{\gamma }_{O_2{O}_4{B}_0}\hfill \\ {\gamma }_{O_3{O}_4{O}_5}={\gamma }_{O_2{O}_4{O}_3}+{\gamma }_{O_2{O}_4{B}_1}+{\gamma }_{O_5{O}_4{B}_1}-2\pi \hfill \\ l_{O_3{O}_5}={\left({l_{O_3{O}_4}}^2+{l_{O_4{O}_5}}^2-2l_{O_3{O}_4}{l}_{O_4{O}_5}cos{\gamma }_{O_3{O}_4{O}_5}\right)}^{1/2}\hfill \\ {\gamma }_{O_5{O}_3{O}_6}=co{s}^{-1}\left(\left({l_{O_3{O}_5}}^2+{l_{O_3{O}_6}}^2-{l_{O_5{O}_6}}^2\right)/\left(2l_{O_3{O}_5}{l}_{O_3{O}_6}\right)\right)\hfill \\ {\gamma }_{O_4{O}_3{O}_5}=co{s}^{-1}\left(\left({l_{O_3{O}_4}}^2+{l_{O_3{O}_5}}^2-{l_{O_4{O}_5}}^2\right)/\left(2l_{O_3{O}_4}{l}_{O_3{O}_5}\right)\right)\hfill \\ {\gamma }_{O_4{O}_3{O}_6}={\gamma }_{O_5{O}_3{O}_6}-{\gamma }_{O_4{O}_3{O}_5}\hfill \end{array}\right.$$

(12)

Then:

$$\left\{\begin{array}{c}{\theta }_3=\pi -{\gamma }_{O_4{O}_3{O}_6}-{\gamma }_{O_6{O}_3{P}_{tip}}-{\gamma }_{O_2{O}_3{O}_4}\hfill \\ {\dot{\theta }}_3=\frac{d{\theta }_3}{d{\gamma }_{O_4{O}_3{O}_6}}{\dot{\gamma }}_{O_4{O}_3{O}_6}\hfill \\ {\ddot{\theta }}_3=\frac{d^2{\theta }_3}{d^2{\gamma }_{O_4{O}_3{O}_6}}{{\dot{\gamma }}_{O_4{O}_3{O}_6}}^2+\frac{d{\theta }_3}{d{\gamma }_{O_4{O}_3{O}_6}}{\ddot{\gamma }}_{O_4{O}_3{O}_6}\hfill \end{array}\right.$$

(13)

The accelerations of the transfer variables in the above derivation are complicated, so methods from the literature [30] can be used. Finally, the following form of the motion mapping relation is obtained by Equations (9)–(13):

$$\left\{\begin{array}{c}{\mathsf{\Theta }}_{joint}={\mathsf{\Phi }}_{k2}\left(s_v,l_{lift},l_{tilt}\right)\hfill \\ {\mathsf{\Theta }}_{joint}={\left[s_v,{\theta }_2,{\theta }_{2^{\prime }},{\theta }_3,{\theta }_4,{\theta }_5\right]}^T\hfill \end{array}\right.$$

(14)

In this section, the kinematic model of the working mechanism in the drive space is established (Equations (6), (7) and (14)). In practical application, when the vehicle displacement and the cylinder lengths are observed, the motion data of any point $P_i$ in the mechanism system can be calculated based on the kinematic model, which provides the basic conditions for the following dynamic modeling of the mechanism.

4. Dynamic Modeling of Working Mechanism System

The complex structure of the loader’s working mechanism and the vehicle’s motion lead to strong dynamic coupling. Firstly, the Lagrange method is used in analytical mechanics to fully establish a mechanism dynamic model in joint space considering the vehicle’s motion. Then, the Newton–Euler method is used in vector mechanics to establish the force mapping relations between the joint space and the drive space. Finally, the dynamic model of the mechanism in drive space is obtained.

4.1. Dynamic Model Based on Lagrange Method

The mechanism system is considered as a complete and ideal restraint rigid system with the generalized coordinates $\left(s_v,{\theta }_2,{\theta }_3\right)$. As mentioned at the beginning, the working surface (ground) on which the vehicle is traveling is not a level surface, i.e., there may be uphill or downhill roads in practical application. This condition results in the pitch motion of the vehicle, which has a great influence on the effect of the system’s gravity in the dynamic model. Therefore, it is necessary to consider the pitch angle of the vehicle in the Earth’s inertial system. Figure 3 shows the pitch motion of the vehicle relative to the Earth’s inertial system $\tilde{O}\tilde{x}\tilde{y}\tilde{z}$. Angle ${\phi }_0$ is the inclination of the working surface to the level surface and can be regarded as the pitch angle of the vehicle. In practical application, the angle ${\phi }_0$ can be obtained by an inclination sensor installed on the vehicle.

In Figure 4, in order to make the mathematical symbols neat and uniform, $m_i$ represents the mass of the body i$\left(i=1,\cdots ,5\right)$; $M_i$ represents the mass center of the body i; ${\mathbf{d}}_i$ represents the vector ${\mathbf{l}}_{\overrightarrow{O_i{M}_i}}$; and ${\phi }_i$ represents the rotation angle of ${\mathbf{d}}_i$ relative to $x_i$-axis. The bucket force can be expressed as follows:

$$\left\{\begin{array}{c}{\mathbf{F}}_{buck}={\left[F_{buck,x_0},F_{buck,y_0}\right]}^T\hfill \\ F_{buck,x_0}=F_{buck}cos{\phi }_R\hfill \\ F_{buck,y_0}=F_{buck}sin{\phi }_R\hfill \end{array}\right.$$

(15)

where ${\mathbf{F}}_{buck}$ is considered as the equivalent force exerted by the operating load on the bucket; $F_{buck,x_0}$ and $F_{buck,y_0}$ are the component forces of ${\mathbf{F}}_{buck}$ on $x_0$-axis and $y_0$-axis in the coordinate frame $O_0{x}_0{y}_0{z}_0$, respectively; Point R is the equivalent action point of ${\mathbf{F}}_{buck}$ on the bucket; and Angle ${\phi }_R$ is the direction angle of ${\mathbf{F}}_{buck}$ in the coordinate frame $O_0{x}_0{y}_0{z}_0$.

For any mechanical systems, the Lagrange function L is defined as follows [31]:

$$L=E_K-E_P$$

(16)

where $E_K$ represents the total kinetic energy of the system and $E_P$ represents the total potential energy of the system. Then, the Lagrange dynamic equation of the system’s dynamic state, described by the Lagrange function L, is given by [31]:

$$Q_i=\frac{d}{dt}\left(\frac{\partial L}{\partial {\dot{\mathsf{\Theta }}}_i}\right)-\frac{\partial L}{\partial {\mathsf{\Theta }}_i}\left(i=1,\cdots ,n\right)$$

(17)

where ${\mathsf{\Theta }}_i$ is the generalized coordinate in the system; ${\dot{\mathsf{\Theta }}}_i$ is the generalized velocity; and $Q_i$ is the non-conservative generalized force on the i-th generalized coordinate.

Considering that the Equation (17) does not explicitly contain ${\dot{\mathsf{\Theta }}}_i$, Equation (17) can be expressed as follows:

$$Q_i=\frac{d}{dt}\left(\frac{\partial E_k}{\partial {\dot{\mathsf{\Theta }}}_i}\right)-\frac{\partial E_k}{\partial {\mathsf{\Theta }}_i}+\frac{\partial E_p}{\partial {\mathsf{\Theta }}_i}$$

(18)

Note that the kinetic energy and the potential energy of each body can be obtained separately. So, Equation (18) can be further expressed as follows:

$$Q_i=\sum _{j=1}^n\frac{d}{dt}\left(\frac{\partial {E_k}_j}{\partial {\dot{\mathsf{\Theta }}}_i}\right)-\sum _{j=1}^n\frac{\partial {E_k}_j}{\partial {\mathsf{\Theta }}_i}+\sum _{j=1}^n\frac{\partial {E_p}_j}{\partial {\mathsf{\Theta }}_i}$$

(19)

The kinetic energy ${E_k}_j$ and the potential energy ${E_p}_j$ of each body can be given by [27,31]:

$$\left\{\begin{array}{c}{E_k}_j=\frac{1}{2}{}^j{\mathbf{v}}_{M_j}^T{\mathbf{m}}_j{}^j{\mathbf{v}}_{M_j}+\frac{1}{2}{}^j{\omega }_j^T{}^M{\mathbf{I}}_j{}^j{\omega }_j\hfill \\ {E_p}_j=m_j{\tilde{\mathbf{g}}}^T{\tilde{\mathbf{l}}}_{\overrightarrow{O_0{M}_j}}\hfill \end{array}\right.$$

(20)

where ${}^j{\mathbf{v}}_{M_j}$ is the velocity of the mass center $M_j$ of the body j in the coordinate frame $O_j{x}_j{y}_j{z}_j$; ${}^j{\omega }_j$ is the angular velocity of the body j in the coordinate frame $O_j{x}_j{y}_j{z}_j$; ${\mathbf{m}}_j\in {\mathbb{R}}^{3\times 3}$ is the diagonal matrix associated with the mass $m_j$; ${}^M{\mathbf{I}}_j\in {\mathbb{R}}^{3\times 3}$ is the inertia tensor of the body j in the coordinate frame ${}^M{O}_j{}^M{x}_j{}^M{y}_j{}^M{z}_j$$\left({}^M{x}_j\Vert x_j,{}^M{y}_j\Vert y_j,{}^M{y}_j\Vert y_j\right)$ at its mass center $M_j$, which is a symmetric matrix; ${\tilde{\mathbf{l}}}_{\overrightarrow{O_0{M}_j}}$ is the vector ${\mathbf{l}}_{\overrightarrow{O_0{M}_j}}$ in the Earth’s inertial frame $\tilde{O}\tilde{x}\tilde{y}\tilde{z}$; and $\tilde{\mathbf{g}}$ is the acceleration vector of gravity in the Earth’s inertial frame $\tilde{O}\tilde{x}\tilde{y}\tilde{z}$. Thus, there are the following expressions:

$$\left\{\begin{array}{c}{\mathbf{m}}_j=\left[\begin{array}{ccc}{m}_j& 0& 0\\ 0& m_j& 0\\ 0& 0& m_j\end{array}\right]\hfill \\ {}^M{\mathbf{I}}_j=\left[\begin{array}{ccc}{}^M{I_j}_x& -{}^M{I_j}_{xy}& -{}^M{I_j}_{xz}\\ -{}^M{I_j}_{yx}& {}^M{I_j}_y& -{}^M{I_j}_{yz}\\ -{}^M{I_j}_{zx}& -{}^M{I_j}_{zy}& {}^M{I_j}_z\end{array}\right]\hfill \\ \tilde{\mathbf{g}}={\left[0,g,0\right]}^T\hfill \\ {\tilde{\mathbf{l}}}_{\overrightarrow{O_0{M}_j}}={\mathbf{R}}_{-1,0}\left({\displaystyle \sum _{i=1}^j}{\mathbf{R}}_{0,i}{}^{i-1}{\mathbf{l}}_{i-1}+{\mathbf{R}}_{0,j}{}^j{\mathbf{d}}_j\right)\hfill \end{array}\right.$$

(21)

where ${}^M{I_j}_x$, ${}^M{I_j}_y$, and ${}^M{I_j}_z$ are the moments of inertia of body j relative to the ${}^M{x}_j$-axis, ${}^M{y}_j$-axis, and ${}^M{z}_j$-axis, respectively, and ${}^M{I_j}_{xy}\left({}^M{I_j}_{yx}\right)$, ${}^M{I_j}_{xz}\left({}^M{I_j}_{zx}\right)$, and ${}^M{I_j}_{yz}\left({}^M{I_j}_{zy}\right)$ are the products of inertia; g is the magnitude of the gravity acceleration vector; and ${\mathbf{R}}_{-1,0}$ is the transformation matrix from the ground coordinate frame $O_0{x}_0{y}_0{z}_0$ to the Earth’s inertial reference frame $\tilde{O}\tilde{x}\tilde{y}\tilde{z}$. Combined with Equation (2) and Figure 3, ${\mathbf{R}}_{-1,0}$ can be obtained by the following equation:

$${\mathbf{R}}_{-1,0}=\mathbf{I}+\left(1-cos{\phi }_0\right){\left({}^0{\tilde{\mathbf{e}}}_0\right)}^2+{}^0{\tilde{\mathbf{e}}}_{0}sin{\phi }_0$$

(22)

In Equation (22), note that ${}^0{\mathbf{e}}_0={\left[0,0,1\right]}^T$ and ${}^0{\tilde{\mathbf{e}}}_0\in {\mathbb{R}}^{3\times 3}$ is the skew-symmetric matrix [27,28] associated with the vector ${}^0{\mathbf{e}}_0$.

Combined with Equations (4), (6), and (20)–(22), in Figure 4, for the kinematic chains 1 $\left(i=1,2,3\right)$ and 2 $\left(i=1,2^{\prime },4,5\right)$, the kinetic energy and the potential energy of the mechanism system can be expressed in the following forms:

$${\mathbf{E}}_{\mathbf{k}}=\left[\begin{array}{c}{E_k}_1\\ ⋮\\ {E_k}_5\end{array}\right]={}_1^{{\mathbf{E}}_{\mathbf{k}}}{\mathbf{D}}^{5\times 5}\left[\begin{array}{c}{\dot{s}}_v^2\\ {\dot{\theta }}_2^2\\ ⋮\\ {\dot{\theta }}_5^2\end{array}\right]+{}_2^{{\mathbf{E}}_{\mathbf{k}}}{\mathbf{D}}^{5\times 10}\left[\begin{array}{c}{\dot{s}}_v{\dot{\theta }}_2\\ {\dot{s}}_v{\dot{\theta }}_3\\ ⋮\\ {\dot{\theta }}_4{\dot{\theta }}_5\end{array}\right]$$

(23)

$${\mathbf{E}}_{\mathbf{p}}=\left[\begin{array}{c}{E_p}_1\\ ⋮\\ {E_p}_5\end{array}\right]={}^{{\mathbf{E}}_{\mathbf{p}}}{\mathbf{D}}^{5\times 5}\left[\begin{array}{c}{m}_{1}g\\ ⋮\\ m_{5}g\end{array}\right]$$

(24)

where ${}_1^{{\mathbf{E}}_{\mathbf{k}}}{\mathbf{D}}^{5\times 5}$ and ${}^{{\mathbf{E}}_{\mathbf{p}}}{\mathbf{D}}^{5\times 5}$ are the ${\mathbb{R}}^{5\times 5}$ matrices, respectively, and ${}_2^{{\mathbf{E}}_{\mathbf{k}}}{\mathbf{D}}^{5\times 10}$ is the ${\mathbb{R}}^{5\times 10}$ matrix.

To facilitate the characterization of the dynamic model, in the following equations, $s\left(i+j+\cdots n\right)$ and $c\left(i+j+\cdots n\right)$ represent $sin\left(i+j+\cdots n\right)$ and $cos\left(i+j+\cdots n\right)$, respectively. $s_{ij\cdots n}$ and $c_{ij\cdots n}$ represent $sin({\theta }_i+{\theta }_j+\cdots +{\theta }_n)$ and $cos({\theta }_i+{\theta }_j+\cdots +{\theta }_n)$, respectively. $s_{ij\cdots n+ij\cdots n}$ and $c_{ij\cdots n+ij\cdots n}$ represent $sin\left({\theta }_i+{\theta }_j+\cdots +{\theta }_n+{\phi }_i+{\phi }_j+\cdots +{\phi }_n\right)$ and $cos\left({\theta }_i+{\theta }_j+\cdots +{\theta }_n+{\phi }_i+{\phi }_j+\cdots +{\phi }_n\right)$, respectively. Note the directions of the vectors and angles. In further considerations, the notations ${\theta }_i$ and ${\phi }_i$ will be represented as signed angles, that is, it is specified that they are positive in the positive direction of the $z_i$-axis.

Unlike open kinematic chains, the mechanism system in this paper belongs to the category of closed six-link chains. Although it is divided into two kinematic chains to facilitate dynamic modeling and characterization, there are still many couplings, e.g., ${\theta }_4$ and ${\theta }_5$ have nonlinear relations with ${\theta }_3$. As a result, the analytical equations of Lagrange function L and its partial derivatives are still complicated. In practical application, the overall size and mass of body 5 are relatively tiny, so its inertia force has little effect on the system compared with the driving force, load and other bodies. In addition, the booms rise and the bucket flips inward during the shoveling operation of the mechanism, so the attitude of body 4 changes slightly (tiny swing) with respect to the working surface. Therefore, the motion of body 4 is considered as a plane translation and the kinetic energy of body 5 is neglected to facilitate the solution, which is validated in the final section that the neglected parts have little effect on the whole system. Thus, the non-zero items of the matrices ${}_1^{{\mathbf{E}}_{\mathbf{k}}}{\mathbf{D}}^{5\times 5}$ and ${}_2^{{\mathbf{E}}_{\mathbf{k}}}{\mathbf{D}}^{5\times 10}$ in Equation (23) are, respectively, as follows:

$$\left\{\begin{array}{c}{}_1^{E_k}{D}_{i1}=m_i/2,\left(i=1,\cdots ,4\right)\hfill \\ {}_1^{E_k}{D}_{ii}=m_i\left({d_i}^2+{}^M{I_i}_z/m_i\right)/2,\left(i=2,3\right)\hfill \\ {}_1^{E_k}{D}_{32}=m_3\left({l_2}^2+{d_3}^2+2l_2{d}_3{c}_{3+3}+{}^M{I_3}_z/m_3\right)/2\hfill \\ {}_1^{E_k}{D}_{42}=m_4{l_{2^{\prime }}}^2/2\hfill \end{array}\right.$$

(25)

$$\left\{\begin{array}{c}{}_2^{E_k}{D}_{21}=-m_2{d}_2{s}_{2+2}\hfill \\ {}_2^{E_k}{D}_{31}=-m_3\left(l_2{s}_2+d_3{s}_{23+3}\right)\hfill \\ {}_2^{E_k}{D}_{32}=-m_3{d}_3{s}_{23+3}\hfill \\ {}_2^{E_k}{D}_{35}=m_3\left({d_3}^2+l_2{d}_3{c}_{3+3}+{}^M{I_3}_z/m_3\right)\hfill \\ {}_2^{E_k}{D}_{41}=-m_4{l}_{2^{\prime }}{s}_{2^{\prime }}\hfill \end{array}\right.$$

(26)

The non-zero items of the matrix ${}^{{\mathbf{E}}_{\mathbf{p}}}{\mathbf{D}}^{5\times 5}$ in Equation (24) are as follows:

$$\left\{\begin{array}{c}{}^{E_p}{D}_{11}=s_{v}s{\phi }_0+d_{1}s\left({\phi }_0+{\phi }_1\right)\hfill \\ {}^{E_p}{D}_{22}=s_{v}s{\phi }_0+l_{1}c{\phi }_0+d_2{s}_{2+02}\hfill \\ {}^{E_p}{D}_{33}=s_{v}s{\phi }_0+l_{1}c{\phi }_0+l_2{s}_{2+0}+d_3{s}_{23+03}\hfill \\ {}^{E_p}{D}_{44}=s_{v}s{\phi }_0+l_{1}c{\phi }_0+l_{2^{\prime }}{s}_{2^{\prime }+0}+d_4{s}_{2^{\prime }4+04}\hfill \\ {}^{E_p}{D}_{55}=s_{v}s{\phi }_0+l_{1}c{\phi }_0+l_{2^{\prime }}{s}_{2^{\prime }+0}+l_4{s}_{2^{\prime }4+0}+d_5{s}_{2^{\prime }45+05}\hfill \end{array}\right.$$

(27)

Combined with Equations (19) and (23)–(27), after strict derivation, the dynamic equations of the non-conservative generalized forces $\mathbf{Q}\left(Q_1,Q_2,Q_3\right)$ in the generalized coordinates $\mathsf{\Theta }\left(s_v,{\theta }_2,{\theta }_3\right)$ are obtained as follows:

$$\begin{array}{c}\mathbf{Q}=\left[\begin{array}{c}{Q}_1\\ Q_2\\ Q_3\end{array}\right]={}_1{\mathbf{D}}^{3\times 3}\left[\begin{array}{c}{\ddot{s}}_v\\ {\ddot{\theta }}_2\\ {\ddot{\theta }}_3\end{array}\right]+{}_2{\mathbf{D}}^{3\times 3}\left[\begin{array}{c}{\dot{s}}_v^2\\ {\dot{\theta }}_2^2\\ {\dot{\theta }}_3^2\end{array}\right]+{}_3{\mathbf{D}}^{3\times 3}\left[\begin{array}{c}{\dot{s}}_v{\dot{\theta }}_2\hfill \\ {\dot{s}}_v{\dot{\theta }}_3\hfill \\ {\dot{\theta }}_2{\dot{\theta }}_3\hfill \end{array}\right]+{}_4{\mathbf{D}}^{3\times 5}\left[\begin{array}{c}{m}_{1}g\\ ⋮\\ m_{5}g\end{array}\right]\hfill \end{array}$$

(28)

where ${}_1{\mathbf{D}}^{3\times 3}$, ${}_2{\mathbf{D}}^{3\times 3}$ and ${}_3{\mathbf{D}}^{3\times 3}$ are the ${\mathbb{R}}^{3\times 3}$ matrices, respectively; ${}_4{\mathbf{D}}^{3\times 5}$ is the ${\mathbb{R}}^{3\times 5}$ matrix.

The items of the matrices ${}_1{\mathbf{D}}^{3\times 3}$, ${}_2{\mathbf{D}}^{3\times 3}$, and ${}_3{\mathbf{D}}^{3\times 3}$ are, respectively, as follows:

$$\left\{\begin{array}{c}{}_1{D}_{11}=\sum m_i\hfill \\ {}_1{D}_{22}=m_2\left({d_2}^2+{}^M{I_2}_z/m_2\right)+m_3\left({l_2}^2+{d_3}^2+2l_2{d}_3{c}_{3+3}+{}^M{I_3}_z/m_3\right)+m_4{l_{2^{\prime }}}^2\hfill \\ {}_1{D}_{33}=m_3\left({d_3}^2+{}^M{I_3}_z/m_3\right)\hfill \\ {}_1{D}_{12}={}_1{D}_{21}=-m_2{d}_2{s}_{2+2}-m_3\left(l_2{s}_2+d_3{s}_{23+3}\right)-m_4{l}_{2^{\prime }}{s}_{2^{\prime }}\hfill \\ {}_1{D}_{13}={}_1{D}_{31}=-m_3{d}_3{s}_{23+3}\hfill \\ {}_1{D}_{23}={}_1{D}_{32}=m_3\left({d_3}^2+l_2{d}_3{c}_{3+3}+{}^M{I_3}_z/m_3\right)\hfill \end{array}\right.$$

(29)

$$\left\{\begin{array}{c}{}_2{D}_{12}=-m_2{d}_2{c}_{2+2}-m_3\left(l_2{c}_2+d_3{c}_{23+3}\right)-m_4{l}_{2^{\prime }}{c}_{2^{\prime }}\hfill \\ {}_2{D}_{13}=-m_3{d}_3{c}_{23+3}\hfill \\ {}_2{D}_{23}=-{}_2{D}_{32}=-m_3{l}_2{d}_3{s}_{3+3}\hfill \\ {}_2{D}_{ii}={}_2{D}_{i1}=0,\left(i=1,\cdots ,3\right)\hfill \end{array}\right.$$

(30)

$$\left\{\begin{array}{c}{}_3{D}_{13}=-2m_3{d}_3{c}_{23+3}\hfill \\ {}_3{D}_{23}=-2m_3{l}_2{d}_3{s}_{3+3}\hfill \\ {}_3{D}_{ii}={}_3{D}_{i1}={}_3{D}_{i2}=0,\left(i=1,\cdots ,3\right)\hfill \end{array}\right.$$

(31)

The items of the matrix ${}_4{\mathbf{D}}^{3\times 5}$ are as follows:

$$\left\{\begin{array}{c}{}_4{D}_{1i}=s{\phi }_0,\left(i=1,\cdots ,5\right)\hfill \\ {}_4{D}_{22}=d_2{c}_{2+02}\hfill \\ {}_4{D}_{23}=l_2{c}_{2+0}+d_3{c}_{23+03}\hfill \\ {}_4{D}_{24}=l_{2^{\prime }}{c}_{2^{\prime }+0}+d_4{c}_{2^{\prime }4+04}\hfill \\ {}_4{D}_{25}=l_{2^{\prime }}{c}_{2^{\prime }+0}+l_4{c}_{2^{\prime }4+0}+d_5{c}_{2^{\prime }45+05}\hfill \\ {}_4{D}_{33}=d_3{c}_{23+03}\hfill \\ {}_4{D}_{34}=d_4{c}_{2^{\prime }4+04}\frac{d{\theta }_4}{d{\theta }_3}\hfill \\ {}_4{D}_{35}=\left(l_4{c}_{2^{\prime }4+0}+d_5{c}_{2^{\prime }45+05}\right)\frac{d{\theta }_4}{d{\theta }_3}+d_5{c}_{2^{\prime }45+05}\frac{d{\theta }_5}{d{\theta }_3}\hfill \\ {}_4{D}_{21}={}_4{D}_{31}={}_4{D}_{32}=0\hfill \end{array}\right.$$

(32)

Furthermore, combined with Equations (14) and (28), the non-conservative generalized force $\mathbf{Q}$ can be expressed as follows:

$$\mathbf{Q}=\left[\begin{array}{c}{Q}_1\\ Q_2\\ Q_3\end{array}\right]={\mathsf{\Phi }}_{d1}(s_v,{\dot{s}}_v,{\ddot{s}}_v,l_{lift},{\dot{l}}_{lift},{\ddot{l}}_{lift},l_{tilt},{\dot{l}}_{tilt},{\ddot{l}}_{tilt})$$

(33)

Next, the Newton–Euler method in vector mechanics will be used to derive the model to build the force mapping relations between the joint space and drive space.

4.2. Dynamic Model Derivation Based on Newton–Euler Method

In Figure 4, the generalized force $\mathbf{Q}$ in the joint space has the following form:

$$\mathbf{Q}=\left[\begin{array}{c}{Q}_1\\ Q_2\\ Q_3\end{array}\right]=\left[\begin{array}{c}{F}_1\\ {\tau }_2\\ {\tau }_3\end{array}\right]$$

(34)

where $F_1$ is the total force of all non-conservative forces driving the vehicle’s motions and ${\tau }_2$ and ${\tau }_3$ are the total torques of all non-conservative forces driving the rotation of joints $O_2$ and $O_3$, respectively.

Then, the Newton–Euler method in vector mechanics is used to derive the dynamic model. Based on the overall analysis of the mechanism system (see Figure 5a), the following relations can be obtained:

$$\left\{\begin{array}{c}{\mathbf{F}}_1={\mathbf{F}}_{vehicle}+{\mathbf{F}}_{buck,x_0}\hfill \\ {\tau }_2={}^2{\mathbf{M}}_{lift}+{}^2{\mathbf{M}}_{tilt}+{}^2{\mathbf{M}}_{buck}\hfill \end{array}\right.$$

(35)

Based on the analysis for the body 3 (see Figure 5b), the following relation can be obtained:

$${\tau }_3={}^3{\mathbf{M}}_{O_6{O}_5}+{}^3{\mathbf{M}}_{buck}$$

(36)

in which:

$$\left\{\begin{array}{c}{}^2{\mathbf{M}}_{lift}={\mathbf{l}}_{\overrightarrow{O_2{A}_0}}\times {\mathbf{F}}_{lift}\hfill \\ {}^2{\mathbf{M}}_{tilt}={\mathbf{l}}_{\overrightarrow{O_2{B}_0}}\times {\mathbf{F}}_{tilt}\hfill \\ {}^2{\mathbf{M}}_{buck}={}^3{\mathbf{M}}_{buck}+{\mathbf{l}}_{\overrightarrow{O_2{O}_3}}\times {\mathbf{F}}_{buck}\hfill \\ {}^3{\mathbf{M}}_{O_6{O}_5}={\mathbf{l}}_{\overrightarrow{O_3{O}_6}}\times {\mathbf{F}}_{O_6{O}_5}\hfill \end{array}\right.$$

(37)

In Equations (35)–(37), ${\mathbf{F}}_{vehicle}$ is the driving force for the vehicle’s motion; ${\mathbf{F}}_{O_6{O}_5}$ is the generalized force exerted by body 5 on body 3; ${}^2{\mathbf{M}}_{lift}$, ${}^2{\mathbf{M}}_{tilt}$, and ${}^2{\mathbf{M}}_{buck}$ are the torques of the lift cylinder force ${\mathbf{F}}_{lift}$, the tilt cylinder force ${\mathbf{F}}_{tilt}$, and the bucket force ${\mathbf{F}}_{buck}$ to the $z_2$-axis, respectively; and ${}^3{\mathbf{M}}_{O_6{O}_5}$ and ${}^3{\mathbf{M}}_{buck}$ are the torques of the forces ${\mathbf{F}}_{O_6{O}_5}$ and ${\mathbf{F}}_{buck}$ to the $z_3$-axis, respectively.

Furthermore, based on the analysis of bodies 4 and 5 (see Figure 5c,d), the following relations can be obtained, respectively:

$$\left\{\begin{array}{c}{\mathbf{l}}_{\overrightarrow{O_4{B}_1}}\times {\mathbf{F}}_{tilt}+{\mathbf{l}}_{\overrightarrow{O_4{O}_5}}\times {\mathbf{F}}_{O_5{O}_6}=0\hfill \\ {\mathbf{F}}_{O_5{O}_6}+{\mathbf{F}}_{O_6{O}_5}=0\hfill \end{array}\right.$$

(38)

where ${\mathbf{F}}_{O_5{O}_6}$ is the generalized force exerted by body 5 on body 4. In addition, in Equations (37) and (38), there are the following expressions:

$$\left\{\begin{array}{c}{\mathbf{l}}_{\overrightarrow{O_2{A}_0}}\times {\mathbf{F}}_{lift}=-F_{lift}{l}_{O_2{A}_0}sin{\gamma }_{O_2{A}_0{A}_1}{}^2{\mathbf{e}}_2\hfill \\ {\mathbf{l}}_{\overrightarrow{O_2{B}_0}}\times {\mathbf{F}}_{tilt}=-F_{tilt}{l}_{O_2{B}_0}sin{\gamma }_{O_2{B}_0{B}_1}{}^2{\mathbf{e}}_2\hfill \\ {\mathbf{l}}_{\overrightarrow{O_2{O}_3}}\times {\mathbf{F}}_{buck}=F_{buck}{l}_{O_2{O}_3}sin\left({\phi }_R-{\theta }_2\right){}^2{\mathbf{e}}_2\hfill \\ {\mathbf{l}}_{\overrightarrow{O_3{O}_6}}\times {\mathbf{F}}_{O_6{O}_5}=F_{O_6{O}_5}{l}_{O_3{O}_6}sin{\gamma }_{O_3{O}_6{O}_5}{}^3{\mathbf{e}}_3\hfill \\ {\mathbf{l}}_{\overrightarrow{O_4{B}_1}}\times {\mathbf{F}}_{tilt}=F_{tilt}{l}_{O_4{B}_1}sin{\gamma }_{O_4{B}_1{B}_0}{}^4{\mathbf{e}}_4\hfill \\ {\mathbf{l}}_{\overrightarrow{O_4{O}_5}}\times {\mathbf{F}}_{O_5{O}_6}=-F_{O_5{O}_6}{l}_{O_4{O}_5}sin{\gamma }_{O_4{O}_5{O}_6}{}^4{\mathbf{e}}_4\hfill \end{array}\right.$$

(39)

where ${\gamma }_{O_2{A}_0{A}_1}$, ${\gamma }_{O_2{B}_0{B}_1}$, ${\gamma }_{O_3{O}_6{O}_5}$, ${\gamma }_{O_4{B}_1{B}_0}$, and ${\gamma }_{O_4{O}_5{O}_6}$ can be expressed by the cylinder lengths $l_{lift}$ and $l_{tilt}$ (refer to the method in Kinematic description Section, for proof see the literature [6]).

In practical application, the direction of the bucket force ${\mathbf{F}}_{buck}$ can be evaluated according to the materials in the bucket at different operation stages [22], that is, ${\phi }_R$ can be estimated. In addition, as the lifting height of the boom is limited in shovel operation, ${\mathbf{F}}_{buck}$ and ${\mathbf{l}}_{\overrightarrow{O_2{O}_3}}$ are not co-linear, that is, $sin({\phi }_R-{\theta }_2)\ne 0$. Based on this, combined with Equations (33)–(39), ${\mathbf{F}}_{buck}$ can be calculated as follows:

$$F_{buck}=\frac{{}^2{M}_{buck}-{}^3{M}_{buck}}{l_{O_2{O}_3}sin\left({\phi }_R-{\theta }_2\right)}$$

(40)

Then, based on Equations (15), (34), (35), and (40), there are the following expressions:

$$\left\{\begin{array}{c}{F}_{vehicle}=Q_1-F_{buck,x_0}\hfill \\ F_{buck,x_0}=F_{buck}cos{\phi }_R\hfill \\ F_{buck,y_0}=F_{buck}sin{\phi }_R\hfill \end{array}\right.$$

(41)

Finally, combined with Equations (33) and (41), the following form of dynamics mapping relation can be obtained:

$$\begin{array}{c}\mathbf{F}=\left[\begin{array}{c}{F}_{vehicle}\hfill \\ F_{buck,x_0}\hfill \\ F_{buck,y_0}\hfill \end{array}\right]={\mathsf{\Phi }}_{d2}(s_v,{\dot{s}}_v,{\ddot{s}}_v,l_{lift},{\dot{l}}_{lift},{\ddot{l}}_{lift},l_{tilt},{\dot{l}}_{tilt},{\ddot{l}}_{tilt},F_{lift},F_{tilt})\hfill \end{array}$$

(42)

In this section, the dynamic model of the working mechanism in drive space is established (Equation (42)). In practical application, when the vehicle displacement, the cylinder lengths and the hydraulic pressures of cylinders are observed, information such as the vehicle’s driving force and the bucket force of the loader can be accurately calculated based on the dynamic model.

5. Validation and Analysis

In practical application, it is difficult to directly and accurately observe the bucket force. So, the simulation method is adopted to validate and analyze the dynamic model by building a reliable simulation platform (Figure 6).

5.1. Validation of the Kinematic Model

Firstly, the kinematic model is validated separately to prove that the kinematic model has no disturbance effect on the dynamic modeling. Following the outline of the counterclockwise outer chain in Figure 6, a virtual prototype based on a 4-${\mathrm{m}}^3$ underground loader is created in the dynamic simulation software ADAMS (see Figure 7). The setting parameters, such as the dimensions and constraints of the virtual prototype, can be seen in the literature [6]. Considering the cooperative motion of the vehicle and mechanism, the Cartesian curve is taken as an example of the planned trajectory of the bucket tip $P_{tip}$, and the motion data can be calculated by the inverse solution of the kinematic model (Equations (6), (7), and (14)). Then, the motion data are used to move the vehicle body and the cylinder pistons of the virtual prototype to achieve the kinematic simulation. The post-processing module of simulation can directly give the trajectory data (called ’calculated trajectory’) of the bucket tip $P_{tip}$. Finally, the calculated trajectory is compared with the planned trajectory to validate the kinematic model. The Cartesian curve function is as follows:

$$\rho =\zeta \left(1-sin\theta \right)$$

(43)

In Figure 8, the bucket tip $P_{tip}$ is planned to start from the lowest point $P_{tip}^{\left(0\right)}$ (initial position) of the Cartesian curve, move clockwise for one cycle and return to the initial position. At the same time, the attitude angle ${\theta }_{P_{tip}}$ of the bucket rotates counterclockwise from the initial angle ${\theta }_{P_{tip}}^{\left(0\right)}$ (bucket bottom is flat), reaches the maximum angle ${\theta }_{P_{tip}}^{\left(1\right)}$ at the center point $P_{tip}^{\left(1\right)}$ of the Cartesian curve, and then rotates clockwise to return to the initial angle. For this, the planned trajectory function of the bucket tip $P_{tip}$ is as follows:

$$\left\{\begin{array}{c}{x}_{P_{tip}}=\zeta \left(1-sin\theta \right)cos\theta +x_{P_{tip}}^{\left(0\right)}\hfill \\ y_{P_{tip}}=\zeta \left(1-sin\theta \right)sin\theta +2\zeta +y_{P_{tip}}^{\left(0\right)}\hfill \\ \theta =3\pi /2-2\pi t/t_{max}\hfill \\ {\theta }_{P_{tip}}=-\left({\theta }_{P_{tip}}^{\left(1\right)}-{\theta }_{P_{tip}}^{\left(0\right)}\right){\left(2t/t_{max}-1\right)}^2+{\theta }_{P_{tip}}^{\left(1\right)}\hfill \end{array}\right.$$

(44)

where $x_{P_{tip}}$, $y_{P_{tip}}$, and ${\theta }_{P_{tip}}$ are the position and attitude parameters of the bucket tip $P_{tip}$ at any t time in the coordinate frame $O_0{x}_0{y}_0{z}_0$, and $x_{P_{tip}}^{\left(0\right)}$, $y_{P_{tip}}^{\left(0\right)}$, and ${\theta }_{P_{tip}}^{\left(0\right)}$ are the initial values. Then, the characteristic coefficient $\zeta$ of the Cartesian curve is set to 1, the total time $t_{max}$ is set to 10 s, and the maximum attitude angle ${\theta }_{P_{tip}}^{\left(1\right)}$ of the bucket is set to 0.6 radians.

Based on the kinematic model (Equations (6), (7), and (14)), the motion data (e.g., see the acceleration data in Figure 9) of the vehicle body and the lift and tilt cylinders of the loader are calculated by the inverse solution of the planned trajectory (Equation (44)) of the bucket tip $P_{tip}$. Then, the motion data are added to the drivers of the vehicle body and the cylinder pistons of the virtual prototype for kinematic simulation. The post-processing module of simulation directly gives the motion trajectory (calculated trajectory) of the bucket tip $P_{tip}$. The calculated trajectory is compared with the planned trajectory (see Figure 10), which shows that the whole process is highly consistent. Figure 11 shows the absolute errors of the x and y coordinates of the calculated trajectory of bucket tip $P_{tip}$. It can be seen that the average absolute error of the x coordinate of calculated trajectory is $1.38\times {10}^{-4}$ m with a maximum absolute error of $2.48\times {10}^{-4}$ m, not exceeding 0.04‰ of the x coordinate of the planned trajectory. The average absolute error of the y coordinate of the calculated trajectory is $5.22\times {10}^{-5}$ m with a maximum absolute error of $1.14\times {10}^{-4}$ m, not exceeding 0.05‰ of the y coordinate of the planned trajectory. These errors can be neglected; thus, the accuracy of the kinematic model is validated.

5.2. Validation and Analysis of the Dynamic Model

Furthermore, the dynamic model is validated and analyzed. Following the outline of the clockwise inner chain in Figure 6, in the motion process of the above virtual prototype, dynamic simulation is carried out by adding a variable force (called ’added force’) on the bucket at any position in the bucket. The post-processing module of simulation can directly obtain the driving force data of the vehicle body and the cylinder pistons. Then, the driving forces data and the motion data of the cylinder pistons are written into the dynamic model (Equation (42)) to calculate the vehicle driving force and the bucket force (called ’calculated force’). Finally, the calculated forces are compared and analyzed with the added force to validate the dynamic model. In addition, the existing dynamic models are also added for comparative analysis to validate the higher accuracy of the dynamic model built in this paper.

Corresponding to the introduction, the existing dynamic models of the loader’s working mechanism can be roughly divided into the following three categories:

(1) Without considering the influence of the vehicle motion, the working mechanism is simplified as the dynamic model of a two-bar mechanism (only the boom and the bucket) [5,7,8,9,10] called ’Model 1’. This model is usually built in the joint space and often neglects the tilt cylinder’s action (affected by the offset load) when switching to the drive space. For this, the Equations (29), (30), (35), and (36) corresponding to such a model can be described as follows:

$$\left\{\begin{array}{c}{}_1{D}_{12}=-m_2{d}_2{s}_{2+2}-m_3\left(l_2{s}_2+d_3{s}_{23+3}\right)\hfill \\ {}_1{D}_{13}=-m_3{d}_3{s}_{23+3}\hfill \\ {}_1{D}_{22}=m_2\left({d_2}^2+{}^M{I_2}_z/m_2\right)+m_3\left({l_2}^2+{d_3}^2+2l_2{d}_3{c}_{3+3}+{}^M{I_3}_z/m_3\right)\hfill \\ {}_1{D}_{33}=m_3\left({d_3}^2+{}^M{I_3}_z/m_3\right)\hfill \\ {}_1{D}_{23}={}_1{D}_{32}=m_3\left({d_3}^2+l_2{d}_3{c}_{3+3}+{}^M{I_3}_z/m_3\right)\hfill \\ {}_1{D}_{i1}=0,\left(i=1,\cdots ,3\right)\hfill \end{array}\right.$$

(45)

$$\left\{\begin{array}{c}{}_2{D}_{12}=-m_2{d}_2{c}_{2+2}-m_3\left(l_2{c}_2+d_3{c}_{23+3}\right)\hfill \\ {}_2{D}_{13}=-m_3{d}_3{c}_{23+3}\hfill \\ {}_2{D}_{23}=-{}_2{D}_{32}=-m_3{l}_2{d}_3{s}_{3+3}\hfill \\ {}_2{D}_{ii}={}_2{D}_{i1}=0,\left(i=1,\cdots ,3\right)\hfill \end{array}\right.$$

(46)

$$\left\{\begin{array}{c}{\mathbf{F}}_1={\mathbf{F}}_{vehicle}+{\mathbf{F}}_{buck,x_0}\hfill \\ {\tau }_2={}^2{\mathbf{M}}_{lift}+{}^2{\mathbf{M}}_{buck}\hfill \\ {\tau }_3={}^3{\mathbf{M}}_{buck}\hfill \end{array}\right.$$

(47)

(2) Without considering the influence of the vehicle motion, a simplified dynamic or static model of the working mechanism with some dynamic compensation [6,11,12,13,14,15,16,17,18,19,20,21,22,23] called ’Model 2’ is established. This model is built in the drive space and incorporates the dynamic effects of the boom and the offset load. For this, the Equations (29)–(31) corresponding to such a model can be described as follows:

$$\left\{\begin{array}{c}{}_1{D}_{12}=-m_2{d}_2{s}_{2+2}-m_3\left(l_2{s}_2+d_3{s}_{23+3}\right)-m_4{l}_{2^{\prime }}{s}_{2^{\prime }}\hfill \\ {}_1{D}_{22}=m_2\left({d_2}^2+{}^M{I_2}_z/m_2\right)+m_3\left({l_2}^2+{d_3}^2+2l_2{d}_3{c}_{3+3}+{}^M{I_3}_z/m_3\right)+m_4{l_{2^{\prime }}}^2\hfill \\ {}_1{D}_{32}=m_3\left({d_3}^2+l_2{d}_3{c}_{3+3}+{}^M{I_3}_z/m_3\right)\hfill \\ {}_1{D}_{i1}={}_1{D}_{i3}=0,\left(i=1,\cdots ,3\right)\hfill \end{array}\right.$$

(48)

$$\left\{\begin{array}{c}{}_2{D}_{12}=-m_2{d}_2{c}_{2+2}-m_3\left(l_2{c}_2+d_3{c}_{23+3}\right)-m_4{l}_{2^{\prime }}{c}_{2^{\prime }}\hfill \\ {}_2{D}_{32}=m_3{l}_2{d}_3{s}_{3+3}\hfill \\ {}_2{D}_{ii}={}_2{D}_{i1}={}_2{D}_{i3}=0,\left(i=1,\cdots ,3\right)\hfill \end{array}\right.$$

(49)

$${}_3{D}_{ij}=0,\left(i,j=1,\cdots ,3\right)$$

(50)

(3) Considering the influence of the vehicle motion, a static model of the working mechanism is established [2,24,25]. The working mechanism of the loader in references [2,24,25] is not exactly the same as that in this study, but their dynamic mechanism models can be classified as a kind of static model that only considers the influence of the vehicle motion but not the dynamic influence of the front-end mechanism. So, a similar model (called ’Model 3’) is built for the reverse six-linkage mechanism in this study drawing on the kind of modeling method in the references [2,24,25] to describe such models. As a comparison group, it is also convenient to analyze the effect caused by various factors more comprehensively and sufficiently in the following text. For this, the Equations (29)–(31) corresponding to such a model can be described as follows:

$$\left\{\begin{array}{c}{}_1{D}_{11}=\sum m_i\hfill \\ {}_1{D}_{21}=-m_2{d}_2{s}_{2+2}-m_3\left(l_2{s}_2+d_3{s}_{23+3}\right)-m_4{l}_{2^{\prime }}{s}_{2^{\prime }}\hfill \\ {}_1{D}_{31}=-m_3{d}_3{s}_{23+3}\hfill \\ {}_1{D}_{i2}={}_1{D}_{i3}=0,\left(i=1,\cdots ,3\right)\hfill \end{array}\right.$$

(51)

$${}_2{D}_{ij}=0,\left(i,j=1,\cdots ,3\right)$$

(52)

$${}_3{D}_{ij}=0,\left(i,j=1,\cdots ,3\right)$$

(53)

The dynamic model (Equation (42)) built in this study is called ’Model Novel’ (Model N).

The variable force added on the bucket of the prototype can be expressed as follows:

$$\left\{\begin{array}{c}{F}_{buck}=\xi sin\left(\pi t\right)\hfill \\ {\phi }_R=-3\pi /2\hfill \end{array}\right.$$

(54)

where the characteristic coefficient $\xi$ is set to 50 KN.

The dynamic simulation is carried out (see Figure 12). The driving force data of the vehicle body and the lift/tilt cylinders are obtained directly by the post-processing module (see Figure 13). Then, the driving force data and the motion data of the vehicle body and lift and tilt cylinders are written into Models 1∼3 and Model N, respectively; thus, the vehicle driving force and the bucket force of each model are calculated. Finally, the calculated forces are compared with the added force and the absolute errors are given (see Figure 14, Figure 15 and Figure 16).

Table 1. Absolute error comparison for each model.

		Model 1	Model 2	Model 3	Model N
	$F_{vehicle}$	11,018.17	11,314.74	328.40	55.72
$Erro{r}_{avg}$/N	$F_{buck,x_0}$	1271.63	322.52	467.70	23.26
	$F_{buck,y_0}$	2202.52	558.62	810.08	40.30
	$F_{vehicle}$	29,347.76	27,032.53	748.03	117.95
$Erro{r}_{max}$/N	$F_{buck,x_0}$	4372.93	946.42	1074.77	63.96
	$F_{buck,y_0}$	7574.14	1639.25	1861.56	110.79

shows the comparison of the average absolute errors ($erro{r}_{avg}$) and the maximum absolute errors ($erro{r}_{max}$) for all models.

In the whole process, the vehicle body and the front-end mechanism move together, and the added force value changes sinusoidally from 0 KN to 50 KN. Combined with Figure 14 and Table 1, in Models 1 and 2, the average absolute errors of the calculated driving force of the vehicle body are 11,018.17 N and 11,314.74 N, respectively, and the maximum absolute errors are 29,347.76 N and 27,032.53 N, respectively. Combined with Figure 9, it can be seen that the excessive errors are sensitive to the acceleration of the vehicle body, which is mainly caused by ignoring the inertia forces generated by the vehicle body motion. In Model 3, the average absolute error of the calculated driving force of the vehicle body is 328.40 N, and the maximum absolute error is 748.03 N. Due to the consideration of the influence of the vehicle motion, the errors perform well, but combined with Figure 9, it can be seen that the errors are sensitive to the acceleration of the front-end mechanism, which is mainly caused by ignoring the inertia forces generated by the front-end mechanism motion. In Model N, the calculated driving force of the vehicle body is highly consistent with the reference driving force in the whole process, with an average absolute error of 55.72 N and a maximum absolute error of 117.95 N. Compared with Models 1∼3, the average and maximum absolute errors of the calculated driving force of the vehicle body of Model N do not exceed 20% of the corresponding errors.

Combined with Figure 15 and Figure 16 and Table 1, in Model 1, the average absolute error of the calculated bucket force in the x direction is 1271.63 N with a maximum absolute error of 4372.93 N, and the average absolute error of the calculated bucket force in the y direction is 2202.52 N with a maximum absolute error of 7574.14 N. The reason for the excessive absolute errors is mainly caused by ignoring the tilt cylinder force (affected by the offset load) and the inertia force generated by the motions of the vehicle body and some components. In Model 2, the average absolute error in the x direction is 322.52 N with a maximum absolute error of 946.42 N, and the average absolute error in the y direction is 558.62 N with a maximum absolute error of 1639.25 N. Due to the dynamic compensation of the front-end mechanism, the errors perform well, but combined with Figure 9, it can be seen that the errors are sensitive to the acceleration of the vehicle body and front-end mechanism, which is mainly caused by ignoring the inertial forces generated by the motions of the vehicle body and some of the components. In Model 3, the average absolute error in the x direction is 467.70 N with a maximum absolute error of 1074.77 N, and the average absolute error in the y direction is 810.08 N with a maximum absolute error of 1861.56 N. Combined with Figure 9, it can be seen that the errors are sensitive to the acceleration of the front-end mechanism, which is mainly caused by ignoring the inertia forces generated by the front-end mechanism motion. In Model N, the calculated bucket force is highly consistent with the added force in the whole process. The average absolute error in the x direction is 23.26 N with a maximum absolute error of only 63.96 N, and the average absolute error in the y direction is 40.30 N with a maximum absolute error of only 110.79 N. Compared with Models 1∼3, the average and maximum absolute errors of the calculated bucket force of Model N do not exceed 10% of the corresponding errors.

In conclusion, the above comparison test fully proves that the dynamics of the vehicle body and front-end mechanism and the force of the tilt cylinder play important roles in improving the model’s accuracy. The dynamic model of the working mechanism system built in this study is superior to the existing dynamic models since it is of higher accuracy and is more suitable for the cooperative motion of the vehicle body and the front-end mechanism.

6. Conclusions

Taking a reverse six-linkage loader as the research object, a kinematic model of the loader’s working mechanism is firstly given. Then, on this basis, an accurate dynamic model of the loader’s working mechanism considering the influence of vehicle motion is established in the 3-DOF motion. Finally, the mapping relations between the vehicle driving force, the bucket force, and the drive space parameters are given.

Since it is difficult to directly and accurately observe the bucket force in practical application, the model in this paper is validated by adding a variable force on the bucket in a simulation platform. In the simulation, the vehicle body and the front-end mechanism move together, and the added force on the bucket changes sinusoidally from 0 KN to 50 KN. The calculated force values based on the established model are highly consistent with the added force values in the whole process. Compared with the existing models, the average and maximum absolute errors of the vehicle driving force calculated by the established model do not exceed 20% of the existing model errors, and the average and maximum absolute errors of the bucket force calculated by the established dynamic model do not exceed 10% of the existing model errors. So, it also proves that the factors, e.g., the dynamics of the vehicle body and front-end mechanism and the force of the tilt cylinder, play important roles in improving the model’s accuracy. The established dynamic model is superior to the existing dynamic models, since it is of higher accuracy and more suitable for the cooperative motion with the vehicle body.

This paper provides important theoretical and practical value for research on the precise online detection of load, the precise identification of mechanical parameters (e.g., system inertia), the optimization of structural design, the identification of vehicle driving force distribution, etc. It provides a reliable basic model for the load detection system of an Intelligent Loader, which can directly and quickly realize the engineering application. The modeling method is also applicable to other mobile multi-bar vehicles. In future work, to better match the conditions of actual loaders, the disturbance effects caused by some non-rigid factors (e.g., the elastic deformation and joint gap of the mechanism, the tire deformation of the vehicle body, changes in the ground, etc.) need to be considered in the model in combination with the mechanism characteristics and operation data of the actual loader, so as to optimize the model to a greater extent. Moreover, in practical application, the model will certainly be applied to the control system of the mechanism for load detection, so the robustness and stability of the model in the control system also need to be further analyzed. The references [32,33] have given good methods and ideas to deal with the above problems.