Multi Body Dynamic Equations of Belt Conveyor and the Reasonable Starting Mode

Based on the rigid finite element method and multibody dynamics, a discrete model of a flexible conveyor belt considering the material viscoelasticity is established. RFE (rigid finite element) and SDE (spring damping element) are used to describe the rigidity and flexibility of a conveyor belt. The dynamic differential equations of the RFE are derived by using Lagrange’s equation of the second kind of the non-conservative system. The generalized elastic potential capacity and generalized dissipation force of the SDE are considered. The forward recursive formula is used to construct the conveyor belt model. The validity of dynamic equations of conveyor belt is verified by field test. The starting mode of the conveyor is simulated by the model.


Introduction
Belt conveyors are widely used in the field of transportation engineering due to their advantages of a strong conveying capacity, a long conveying distance, high speed, stable operation and low noise [1]. A typical belt conveyor is composed of conveyor belts, drive rollers, redirecting rollers, idlers, tensioning devices and braking devices. It is a complex long-distance rigid-flexible coupling mechanical system. Therefore, the establishment of an easy-to-solve and accurate calculation system model has significant engineering significance for the analysis of the dynamic performance of a conveyor and the prediction of fatigue damage [2,3]. The core component of a conveyor-the material of a conveyor belt-has relatively good mechanical properties. It cannot be treated simply as a rigid or elastic body [4]. The model that closely resembles its essential characteristics is the viscoelastic model [5,6], which makes the dynamic analysis model of the conveyor belt much more complex than that of general mechanical parts. The dynamic analysis model of conveyor belt mainly includes continuous model and discrete model. The continuation model refers to the various dynamic differential equations commonly used by early researchers studying belt conveyors [7]. The establishment of these models promotes the improvement of the design of belt conveyors. However, the differential equations for the boundary conditions and initial conditions have become more complex. Discrete models are derived from continuous models. The modeling accuracy and engineering application value are continuously improved [8,9].
The earliest finite element discrete model is a one-dimensional longitudinal vibration analysis model [10,11], which is essentially based on the concentrated mass finite element method. It regards the continuous conveyor belt as a system composed of a series of particles connected by springs and damping elements. The quality and force of the conveyor belt, including the mass of the bulk material, the running resistance of the rollers, etc. [12,13], are concentrated at a node. The dynamic equations of the entire system are linear differential equations. This type of model has achieved good results in the longitudinal vibration analysis of a conveyor belt. Li [14] used the Kelvin-Voigt viscoelastic model to establish continuous dynamic equations for tail hammer tension belt conveyors. Various factors influencing the longitudinal vibration of the belt conveyor were simulated. The two-dimensional finite Symmetry 2020, 12, 1489 2 of 21 element model mainly refers to the mixed finite element model of rod and beam elements [15]. This method can comprehensively consider the coupling problem associated with the longitudinal and lateral vibrations of a conveyor belt. To investigate belt deflection for variable dynamic conditions, Shen [16] developed a coupled finite element method (FEM) and discrete element method (DEM) model to predict the dynamic belt deflection. Ilic [17] established a conveyor belt model with the DEM to study the interaction between bulk solids and the belt during transportation. There was a good correlation between the experimental data and the simulation results. The calculation method of the conveyor roller load has improved.
Although the finite element method has achieved great success in the analysis of the longitudinal and transverse vibration characteristics of a conveyor belt, in essence, the various belt conveyor dynamic analysis methods based on the finite element method still belong to the structural dynamic analysis. Most existing models do not fully consider the nonlinear problem of conveyor belts. The flexible conveyor belt exhibits viscoelasticity in terms of mechanical properties. During the working process, there is a relatively large deformation relative to other parts. This coupling of large displacement, large deformation and viscoelasticity on the conveyor belt makes the nonlinear problem in its dynamic analysis very prominent. Many researchers use a variety of relatively simple linear models to study the dynamic characteristics of belt conveyors, which will inevitably have an adverse effect on the accuracy of its analysis. The dynamic analysis method for a belt conveyor system must still be established based on the multibody system dynamics method.
The rigid finite body element method is a flexible multibody dynamic analysis method based on the theory of multibody dynamics. It divides the flexible or rigid flexible multibody system into a series of rigid finite elements (RFE) connected by spring damping elements (SDE). The main difference between the RFE method and the FE method is that the RFE method uses a rigid finite body element to describe the inertia forces of a flexible body. The element flexibility is described relative to the rigid body element's conjoined coordinate system, instead of the flexible body's body coordinate system. The small deformation between rigid body elements is used to describe the large deformation of the entire flexible body. This naturally takes into account the geometric nonlinearity effects of the dynamic stiffness, without the need for a geometric stiffness matrix. Flexible multibody systems and rigid-flexible hybrid multibody systems can be handled in a unified way. Nonlinear spring damping element parameters can be used to conveniently solve material nonlinearity problems. For example, the relationship between the rolling resistance of a roller and the parameters of the equipment system can be described using spring damping elements [18]. Therefore, the rigid finite body element method is particularly suitable for the computer simulation modelling of a rigid-flexible complex mechanical system composed of a flexible body (conveyor belt) and a rigid body (roller, roller system, etc.), such as a belt conveyor.
In this paper, a flexible multibody dynamic model of a conveyor belt is established based on the RFE method. The model is used to study the dynamic performance of the conveyor under different starting modes.

Discrete Model of the Conveyor Belt
As shown in Figure 1, a conveyor system consists of p-1 parts and conveyor belt p. The conveyor belt p can be discretized into m p rigid body elements (RFEs) reflecting the rigid characteristics of the conveyor belt, and the rigid body elements are connected by spring damping elements (SDE) reflecting the flexible characteristics of the conveyor belt. Then, the whole conveyor system is simplified into a rigid flexible multi-body system composed of p-1 rigid bodies, m p RFEs and m p SDEs. The simplified dynamic model of conveyor belt is symmetrical in structure.

Coordinate Description
(1) Rigid finite element  The constant transformation matrix from the local coordinate system l i e of the RFE (p, i) to the reference coordinate system r e is as follows:  Figure 2 shows the coordinate system of the RFE of the conveyor belt. It includes the reference coordinate system e r , local coordinate system e l i and floating coordinate system e b i of the RFE (p, i) [19].

Coordinate Description
(1) Rigid finite element    The reference coordinate system e r is selected at the first RFE (p, 1), and each RFE has a local coordinate system e l i . The elastic displacement of the RFE (p, 1) is described by the position and altitude of e l i relative to e r before belt deformation. The origin of e l i is at the initial position C 0i of the centroid of the RFE. The three coordinate axes coincide with the main inertia axis. Before deformation, the floating coordinate system e b i coincides with the local coordinate system e l i and moves with the RFE. The constant transformation matrix from the local coordinate system e l i of the RFE (p, i) to the reference coordinate system e r is as follows: where A rl i is the direction cosine matrix of e l i with respect to e r , and r i is the coordinate vector of the origin of e l i in e r . Unlike rigid bodies, the conveyor belt is elastic, and the relative positions of various points in the body are constantly changing. Therefore, the position and altitude of the floating coordinate system e b i relative to the local coordinate system e l i can be described by the following generalized coordinates: where x i is a coordinate vector of the origin of e b i in e l i , and ϕ i is an altitude coordinate vector of e b i in e l i . According to the above equation, the transformation matrix of e b i to e l i for RFE (p, i) is: where, where, C i1 = cos ϕ i1 , S i1 = sin ϕ i1 . Therefore, the uniform transformation matrix of e b i to e r for RFE (p, i) is:

Spring Damping Element (SDE)
The deformation of the conveyor belt caused by external and inertial forces is described by the SDE. The coordinate system of SDE (p, j), connecting RFE (p, j) and RFE (p, k), is shown in Figure 3. To describe the flexible characteristics of the conveyor belt, a local coordinate system e j is established at SDE (p, j) before its deformation. Coordinate systems e s j and e d k are established at the connection points s (between RFE (p, j) and SDE (p, j)) and d (between RFE (p, k) and SDE (p, j)), respectively. The three coordinate systems are now coincident.   The positions and altitudes of the SDE (p, j) coordinate system e j at the local coordinate systems e l j and e l k of RFE (p, j) and RFE (p, k) are represented by the following constant transformation matrices: It should be noted that the constant transformation matrices of the coordinate systems e s j and e d k at points s and d with respect to their respective floating coordinate systems are equal to the above equations. Due to the load on the conveyor belt, RFE (p,j) and RFE (p,k) will produce elastic deformation. Therefore, the position and altitude of the floating coordinate system e b i relative to the respective local coordinate systems e l j and e l k can be described by the generalized coordinate given in Equation (2): From Equations (6) and (7), and considering Equation (3), the position coordinate vectors of points s and d relative to their respective local coordinate systems e l j and e l k can be calculated by the following equation: Considering Equation (6) and its inverse matrix, from Equation (8), the position coordinate vector of points s and d relative to the SDE (p,j) coordinate system e j can be obtained ( Figure 4):  From the above equation, the relative displacement matrix of points s and d at SDE (p,j) caused by the deformation of the conveyor belt can be obtained: Similarly, the relative rotation angular displacement matrix is:

Spring Damping Element Parameters
(1) Potential energy and generalized stiffness matrix The potential energy of the SDE (p,j) is divided into two terms. One is caused by the From the above equation, the relative displacement matrix of points s and d at SDE (p,j) caused by the deformation of the conveyor belt can be obtained: Similarly, the relative rotation angular displacement matrix is: Symmetry 2020, 12, 1489 6 of 21

Potential Energy and Generalized Stiffness Matrix
The potential energy of the SDE (p,j) is divided into two terms. One is caused by the displacement changes between elements j and k: The other term is caused by the rotation change: where ∆x j and ∆ϕ j are defined by Equations (10) and (11). C j x and C j ϕ are generalized stiffness matrices: where c j xi and c j ϕi are the translation stiffness coefficient and rotation stiffness coefficient, respectively. Then, the total elastic potential energy of SDE (p,j) is: Similarly, the total energy loss function of SDE (p, j) is: where D Each RFE can be simplified into an equal length beam element with a rectangular section shape; ii.
The deformation mode and velocity of the actual conveyor belt are the same as that of its equivalent SDE; iii. The stress of each section in the conveyor belt section is equal, and its mechanical properties can be described by the Voigt model [20]: where σ(t) is the belt stress, ε(t) is the belt strain, . ε(t) is the strain rate, E is the Elastic modulus, and η is the viscosity coefficient.
Based on the above assumptions, the following relationships can be obtained: where F 1 (t) is the tension of the belt element, ∆w 1 (t) is the longitudinal deformation of the belt element, ∆ . w 1 (t) is the longitudinal deformation rate of the belt element, A is the cross-sectional area of the belt element, and ∆l is the belt element length.

Kinetic Energy and Lagrangian Operator
There are two kinds of simplified rigid bodies that can be used when modelling rigid flexible hybrid multibody systems or flexible multibody systems with the rigid finite body element method. One is the natural rigid body, that is, the actual object whose elastic deformation is negligible compared with that of other objects in the system. The other is the RFE formed after the discretization of the flexible body, which must be considered after deformation. Although these two kinds of rigid bodies have different choices for generalized coordinates, the method of establishing the kinematic equation is the same. In order to simplify the derivation process, the two kinds of rigid bodies will not be distinguished. The generalized coordinates of a rigid body in the system are uniformly expressed as follows: where q k is the generalized coordinate component.
The dynamic differential equations of a rigid body given by Lagrange's equations of the second kind for a nonconservative system are [21]: where T is the kinetic energy, V is the potential energy, including the gravitational energy and elastic energy, . q k is the generalized velocity component, and Q k is the generalized force components of nonconservative systems. The Lagrange operator is: The matrix form of Equation (18) is: As shown in Figure 5, the mass of a particle on a rigid body is dm, and the position vector diameter in the rigid body connected coordinate system e b is r . In the global coordinate system, the vector of e is r. e r is a moving reference frame. For a natural rigid body, e r is the coordinate system between adjacent bodies. For rigid elements of a flexible body, e r is the coordinate system of the starting rigid element.
where T is the kinetic energy, V is the potential energy, including the gravitational energy and elastic energy, k q  is the generalized velocity component, and k Q is the generalized force components of nonconservative systems. The Lagrange operator is: The matrix form of Equation (18) is: As shown in Figure 5, the mass of a particle on a rigid body is m d , and the position vector diameter in the rigid body connected coordinate system b e is r . In the global coordinate system, the vector of e is r .  Then, the relationship between the position and diameter of a particle in different coordinate systems can be expressed as: is the coordinate matrix for a particle in the rigid body coordinate Then, the relationship between the position and diameter of a particle in different coordinate systems can be expressed as: where, T is the coordinate matrix for a particle in the rigid body coordinate system e b ; T is the coordinate matrix for a particle in the global coordinate system e.
B b is the uniform transformation matrix of the rigid body coordinate system e b relative to global coordinate system e, which is abbreviated as B.
Then, the kinetic energy of a rigid body is: where H is the so-called pseudo inertia matrix at the origin C of the rigid body relative coordinate system e b . Equation (23) can be obtained by substituting Equation (22) into Equation (19) and omitting the derivation process: where, Equation (23) can be written in matrix form as: where A = (a k,i ) k,i=1,...,n and e = (e k ) k=1,...,n

Gravitational Potential Energy and Generalized Force
If the mass of the rigid body element is m, gravity is g, and the vector diameter of the rigid body centre in coordinate system e b is expressed by r C , then the gravitational potential energy can be calculated as: where θ 3 = 0 0 1 0 is the transformation matrix. Matrix B is the same as that the above section. Deriving the generalized coordinate q k for the rigid body from the above equation gives: Equation (26) written in matrix form is: Symmetry 2020, 12, 1489 9 of 21 where g = (g k ) k=1,...,n and g k = mgθ 3 B k r C If the coordinate matrix of the external force P acting on a point on the rigid body in coordinate system e b is P and the coordinate matrix of the position vector is r , then the generalized force corresponding to the generalized coordinate q k of the rigid body is: If the coordinate matrix of the external force couple M acting on the rigid body in coordinate system e b is: Then, the generalized force corresponding to the generalized coordinate q k of the rigid body is: The generalized force matrix of the rigid body is: Substituting Equations (24), (27), and (31) into Equation (18), the dynamic equation of the rigid body can be obtained as follows: where F = Q − e − g.

Elastic Potential Energy
For a flexible body such as a conveyor belt, its deformation is described by a spring damping element. Therefore, the dynamic equation of the rigid body element should also include the generalized force term corresponding to the potential energy caused by elastic deformation. From Section 2.3, the elastic deformation potential energy of SDE (p, j) connecting RFE (p, j) and RFE (p, k) of flexible body p is: Then, the total elastic potential energy of the flexible body p can be obtained from the sum of the elastic deformation potential energy of all m p spring damping elements. Then, the total elastic potential energy of p can be calculated from the sum of the potential energy of all m p SDEs: The generalized force corresponding to the generalized coordinateq The generalized elastic potential capacity matrix corresponding to all the generalized coordinates q (p) of the flexible body p can be written as: Similarly, the generalized dissipative force matrix corresponding to all generalized coordinates of p is:

Forward Recursive Formulation
Theoretically, the system's equations of motion can be obtained by combining the dynamic equations of all the rigid bodies or elements in the system. However, the classical rigid finite element method uses two types of generalized coordinates for a rigid-flexible hybrid system (hinged Lagrange coordinates for the adjacent objects in the system and generalized coordinates of Cartesian and Euler angles for the rigid elements of the belt). Thus, the right term of dynamic Equation (42) for each object in the system's motion chain and matrix A will contain the topological information of all the low-order bodies in its path. The form is quite complicated. Therefore, the dynamic equation of the entire system needs to be established using recursive groups, starting from the low-order body of the system [22].
Consider a system consisting of p − 1 rigid bodies and flexible body p. To simplify the description, flexible body p is not initially discrete, and is treated in the same way as the rigid bodies. The generalized coordinate vector of the flexible body p in the global coordinate system and the transformation matrix of its local coordinate system to the global coordinate system can be expressed as: where q (p) is the generalized coordinate vector of the flexible body p relative to the rigid body p − 1, which can be expressed by:q is the transformation matrix of the local coordinate system of the flexible body p relative to the local coordinate system of the rigid body p − 1; B (p−1) is the transformation matrix of the rigid body p − 1 in the local coordinate system relative to the global coordinate system.
Note that the number of generalized coordinates n p of the flexible body p is dependent on the number of generalized coordinates of all its low-order bodies and the number of generalized coordinates of itself relative to the adjacent low-order rigid body p − 1, i.e., The dynamic equation of the flexible body p described by Equation (32) can be rewritten as follows: It can also be expanded into: q (1) . . .
Equation (41) reveals the relationship between the motion of the flexible body p and its own coordinatesq (p) and the generalized coordinates of all its low-order rigid bodies. See Appendix A for details. The recursive relationship can be simplified as: The kinetic energy T (p) and the gravitational potential energy V (p) g of the system can be derived from the sum of the kinetic and gravitational potential energy of all the objects in the system. The recursive relationship for energy is: The dynamic equation of a system composed of p − 1 rigid bodies can be written as: Then, the dynamic equation of the entire system composed of Equations (42) and (45) can be expressed as:

Dynamic Equation of the Rigid Flexible Multibody System
If the abovementioned flexible body p is discretized into m p rigid body elements, the generalized coordinate vector describing the i-th rigid body element in the global coordinate system and the transformation matrix of its local coordinate system to the global coordinate system can be expressed by the following equation: Then, the number of generalized coordinates describing the motion of RFE (p, i) of the flexible body p can be expressed as: n p,i = n p−1 + n p,1 + 6 = n p,1 + 6 (48) The generalized coordinate vector of the entire flexible body p is expressed by the following equation: where,q (p) = q(p,2) T · · ·q (p,m p ) T Then, the dynamic equation of the entire flexible body p can be expressed as: See Appendix A for details. Substituting the equations of motion of the flexible body p into dynamic Equation (45), the dynamic equations of the whole rigid-flexible multibody system can be obtained as follows: See Appendix A for details.

Field Test
In order to verify the correctness of the flexible multi-body dynamics modelling method of the belt conveyor in this paper, a DTII (A) 10040 type material belt conveyor of a cement plant was selected as an example for simulation calculation. The main parameters of the machine are as follows: The arrangement of the conveyor is as follows (Figure 6): The arrangement of the conveyor is as follows (Figure 6): Since the sensor cannot be installed in the belt, the indirect test method is adopted. Two pressure sensors with a range of 3 tons are installed under the upper roller about 10 m away from the machine head, as shown in Figure 7a. The height of the sensor is about 5 cm, and the dynamic tension in the cross section of the belt can be calculated by recording the signal of the pressure sensor. It can be calculated by formula Fd =2Fc/sinθ, where Fd is the dynamic tension in the cross section of the belt, Fc is the output pressure value of a single pressure sensor, and θ is the angle between the belt pushed up by the idler and the horizontal plane.
The running speed of the belt is measured by rolling a friction wheel with the belt without relative sliding. A disc is installed at the other end of the friction wheel shaft to rotate synchronously with the friction wheel, as shown in Figure 7b. The outer ring of the disc is engraved with a series of slots, and 121 pulses can be output through the photodiode when the disc rotates one circle. The instantaneous velocity of the belt can be calculated by the relationship between the pulse distance and the disc diameter. (2) Simulation model The multi-body dynamics simulation model of the machine is shown in Figure 8, in which the model includes the resistance of roller squeeze [23], the friction contact between the friction pulley and the conveyor belt [24] and the running resistance [25]. The model is established through the code we developed, and a variable topology band length solution is used. First, the conveyor belt is equally divided into N belt sections, and the belt sections are connected by an SDE. The length of the belt section is based on being able to meet the contact with the pulley. The belt sections that are not in contact with the pulley are set as a group of M, and fixed constraints are applied between the belt sections to form a long belt section. Before each simulation time step is carried out, if a certain belt section may be in contact with the pulley during this step, the constraint between the long belt sections will be dynamically lifted and converted into a normal belt section. Fixed constraints are Since the sensor cannot be installed in the belt, the indirect test method is adopted. Two pressure sensors with a range of 3 tons are installed under the upper roller about 10 m away from the machine head, as shown in Figure 7a. The height of the sensor is about 5 cm, and the dynamic tension in the cross section of the belt can be calculated by recording the signal of the pressure sensor. It can be calculated by formula F d = 2F c /sinθ, where F d is the dynamic tension in the cross section of the belt, F c is the output pressure value of a single pressure sensor, and θ is the angle between the belt pushed up by the idler and the horizontal plane. Since the sensor cannot be installed in the belt, the indirect test method is adopted. Two pressure sensors with a range of 3 tons are installed under the upper roller about 10 m away from the machine head, as shown in Figure 7a. The height of the sensor is about 5 cm, and the dynamic tension in the cross section of the belt can be calculated by recording the signal of the pressure sensor. It can be calculated by formula Fd =2Fc/sinθ, where Fd is the dynamic tension in the cross section of the belt, Fc is the output pressure value of a single pressure sensor, and θ is the angle between the belt pushed up by the idler and the horizontal plane.
The running speed of the belt is measured by rolling a friction wheel with the belt without relative sliding. A disc is installed at the other end of the friction wheel shaft to rotate synchronously with the friction wheel, as shown in Figure 7b. The outer ring of the disc is engraved with a series of slots, and 121 pulses can be output through the photodiode when the disc rotates one circle. The instantaneous velocity of the belt can be calculated by the relationship between the pulse distance and the disc diameter. (2) Simulation model The multi-body dynamics simulation model of the machine is shown in Figure 8, in which the model includes the resistance of roller squeeze [23], the friction contact between the friction pulley and the conveyor belt [24] and the running resistance [25]. The model is established through the code we developed, and a variable topology band length solution is used. First, the conveyor belt is equally divided into N belt sections, and the belt sections are connected by an SDE. The length of the belt section is based on being able to meet the contact with the pulley. The belt sections that are not in contact with the pulley are set as a group of M, and fixed constraints are applied between the belt sections to form a long belt section. Before each simulation time step is carried out, if a certain belt section may be in contact with the pulley during this step, the constraint between the long belt sections will be dynamically lifted and converted into a normal belt section. Fixed constraints are The running speed of the belt is measured by rolling a friction wheel with the belt without relative sliding. A disc is installed at the other end of the friction wheel shaft to rotate synchronously with the friction wheel, as shown in Figure 7b. The outer ring of the disc is engraved with a series of slots, and 121 pulses can be output through the photodiode when the disc rotates one circle. The instantaneous velocity of the belt can be calculated by the relationship between the pulse distance and the disc diameter.

Simulation Model
The multi-body dynamics simulation model of the machine is shown in Figure 8, in which the model includes the resistance of roller squeeze [23], the friction contact between the friction pulley and the conveyor belt [24] and the running resistance [25]. The model is established through the code we developed, and a variable topology band length solution is used. First, the conveyor belt is equally divided into N belt sections, and the belt sections are connected by an SDE. The length of the belt section is based on being able to meet the contact with the pulley. The belt sections that are not in contact with the pulley are set as a group of M, and fixed constraints are applied between the belt sections to form a long belt section. Before each simulation time step is carried out, if a certain belt section may be in contact with the pulley during this step, the constraint between the long belt sections will be dynamically lifted and converted into a normal belt section. Fixed constraints are added between ordinary belt sections that are out of contact to form long belt sections. The input parameters are shown in Table 1.

Model parameters and running resistance
Number of belt sections: 170 System degree of freedom: 309 Both the measured data and the simulation results show that the peak value of dynamic tension are significantly reduced 15 s after starting, and the dynamic tension tends to be stable. Therefore, this article only makes a comparative analysis of the data within 15 s after starting.

Verification
The comparison between the simulation results and the test data is shown in Figures 9 and 10. It can be seen that the changing waveforms and frequencies of the two are relatively consistent, indicating the practicality of the developed dynamic model. About 12 s after starting, the conveyor belt enters the normal working speed, and the fluctuation amplitude of the same starting tension tends to be stable. The simulation results show that the speed fluctuation of the conveyor belt is relatively gentle 0-12 s after the start of the conveyor, and the test results fluctuate obviously. This is mainly because the boundary conditions of the actual operation of the conveyor are much more complex than those given by the simulation. Overall, the dynamic analysis method of the flexible multi-body of conveyor belt is effective and accurate.  Both the measured data and the simulation results show that the peak value of dynamic tension are significantly reduced 15 s after starting, and the dynamic tension tends to be stable. Therefore, this article only makes a comparative analysis of the data within 15 s after starting.

Verification
The comparison between the simulation results and the test data is shown in Figures 9 and 10. It can be seen that the changing waveforms and frequencies of the two are relatively consistent, indicating the practicality of the developed dynamic model. About 12 s after starting, the conveyor belt enters the normal working speed, and the fluctuation amplitude of the same starting tension tends to be stable. The simulation results show that the speed fluctuation of the conveyor belt is relatively gentle 0-12 s after the start of the conveyor, and the test results fluctuate obviously. This is mainly because the boundary conditions of the actual operation of the conveyor are much more complex than those given by the simulation. Overall, the dynamic analysis method of the flexible multi-body of conveyor belt is effective and accurate.

Discussions
Since the dynamic tension of the belt conveyor reaches the maximum value in the starting stage, it is very important to study the factors affecting the dynamic tension in the starting process. There are many factors affecting the dynamic tension, including starting mode and time, driving mode, tension device and starting condition, etc. The ideal starting process should ensure that the acceleration curve changes gently without sudden change, and the peak value of acceleration is small. The ideal starting process should ensure that the acceleration curve changes smoothly, without abrupt changes, and the peak value of acceleration is small. US HARRISON [26] recommends the following starting speed and acceleration curve: where, v0 is the initial velocity (default 0), t is the time variable, T is the total acceleration time (staring time).
Nordell [26] gives the following starting speed curve:

Discussions
Since the dynamic tension of the belt conveyor reaches the maximum value in the starting stage, it is very important to study the factors affecting the dynamic tension in the starting process. There are many factors affecting the dynamic tension, including starting mode and time, driving mode, tension device and starting condition, etc. The ideal starting process should ensure that the acceleration curve changes gently without sudden change, and the peak value of acceleration is small. The ideal starting process should ensure that the acceleration curve changes smoothly, without abrupt changes, and the peak value of acceleration is small. US HARRISON [26] recommends the following starting speed and acceleration curve: where, v 0 is the initial velocity (default 0), t is the time variable, T is the total acceleration time (staring time). Nordell [26] gives the following starting speed curve: The starting curve of Nordell is a quadratic curve. In this paper, the following cubic starting speed curve is constructed by using STEP function: The common starting speed curve also has the equal acceleration curve: Figure 11 shows the speed and acceleration curves of four starting modes with 10 s starting time, and all of them reach the required speed of 2 m/s after the starting time is over. The common starting speed curve also has the equal acceleration curve: Figure 11 shows the speed and acceleration curves of four starting modes with 10 s starting time, and all of them reach the required speed of 2 m/s after the starting time is over.  Figure 12 shows the tension on the tight side of the conveyor belt near the drive drum calculated by the model. It can be seen that: the tension amplitude of the four starting modes tends to be the same. Due to the maximum instantaneous acceleration, the maximum tension of isoacceleration start appears earliest. The amplitude of tension fluctuation is the largest, but the maximum acceleration is the smallest. Therefore, in the case of short starting time, the maximum tension value is not different from the other three curves. The other three curves have little difference, among which Nordell curve has the largest maximum tension value and the smallest fluctuation.  Figure 12 shows the tension on the tight side of the conveyor belt near the drive drum calculated by the model. It can be seen that: the tension amplitude of the four starting modes tends to be the same. Due to the maximum instantaneous acceleration, the maximum tension of isoacceleration start appears earliest. The amplitude of tension fluctuation is the largest, but the maximum acceleration is the smallest. Therefore, in the case of short starting time, the maximum tension value is not different from the other three curves. The other three curves have little difference, among which Nordell curve has the largest maximum tension value and the smallest fluctuation.
The simulation results show that reasonable increase of starting time can reduce the acceleration peak value, thus effectively reducing the starting tension of the conveyor belt. Therefore, the true starting time is increased from 10 s to 30 s. Figure 13 shows the belt tension curve with a starting time of 30 s. It can be seen that when the peak acceleration decreases, the maximum starting tension of conveyor belt decreases from 110 kN to less than 80 kN, and the fluctuation amplitude of tension also decreases obviously. The overall shape of tension curve of 10 s and 30 s starting modes is basically consistent. It shows that for the large belt conveyor, there should be enough long starting time to ensure a small dynamic tension when starting. The simulation results show that reasonable increase of starting time can reduce the acceleration peak value, thus effectively reducing the starting tension of the conveyor belt. Therefore, the true starting time is increased from 10 s to 30 s. Figure 13 shows the belt tension curve with a starting time of 30 s. It can be seen that when the peak acceleration decreases, the maximum starting tension of conveyor belt decreases from 110 kN to less than 80 kN, and the fluctuation amplitude of tension also peak value, thus effectively reducing the starting tension of the conveyor belt. Therefore, the true starting time is increased from 10 s to 30 s. Figure 13 shows the belt tension curve with a starting time of 30 s. It can be seen that when the peak acceleration decreases, the maximum starting tension of conveyor belt decreases from 110 kN to less than 80 kN, and the fluctuation amplitude of tension also decreases obviously. The overall shape of tension curve of 10 s and 30 s starting modes is basically consistent. It shows that for the large belt conveyor, there should be enough long starting time to ensure a small dynamic tension when starting. Among the four starting curves, the tension of Harrison curve is the most ideal. The tension of STEP function curve is similar to it. Due to the belt sagging in static state, a delay of about 3~8 s is often added in the initial starting stage of large belt conveyor. Thus, the conveyor belt can be tensioned before the acceleration and torque increase. Each component can enter the running state at a lower speed, further reducing the peak value of starting tension and increasing the running stability. Such a starting mode is also called pre starting. In order to study the influence of pre starting on the dynamic characteristics of belt conveyor, this paper adds a 3 s delay section in four starting speed curves with 10 s starting time. Figure 14 shows the corresponding speed and acceleration curves. Among the four starting curves, the tension of Harrison curve is the most ideal. The tension of STEP function curve is similar to it. Due to the belt sagging in static state, a delay of about 3~8 s is often added in the initial starting stage of large belt conveyor. Thus, the conveyor belt can be tensioned before the acceleration and torque increase. Each component can enter the running state at a lower speed, further reducing the peak value of starting tension and increasing the running stability. Such a starting mode is also called pre starting. In order to study the influence of pre starting on the dynamic characteristics of belt conveyor, this paper adds a 3 s delay section in four starting speed curves with 10 s starting time. Figure 14 shows the corresponding speed and acceleration curves. It can be seen from Figure 15 that, when starting without delay for 10 s, the tension continuously rises to the maximum value within 2~3.5 s, and then oscillates down and tends to be stable. When the starting time is 10 s and the delay is 10 s, the tension begins to decrease gradually when it reaches 80 kN in about 2 s. After the delay is over, it rises slowly again until it reaches the maximum value, and then oscillates down and tends to be stable. No matter what starting mode, the maximum tension is about 95 kN after pre-starting, which is less than 105 kN without delay starting. It can be seen from Figure 15 that, when starting without delay for 10 s, the tension continuously rises to the maximum value within 2~3.5 s, and then oscillates down and tends to be stable. When the starting time is 10 s and the delay is 10 s, the tension begins to decrease gradually when it reaches 80 kN in about 2 s. After the delay is over, it rises slowly again until it reaches the maximum value, and then oscillates down and tends to be stable. No matter what starting mode, the maximum tension is about 95 kN after pre-starting, which is less than 105 kN without delay starting. It can be seen from Figure 15 that, when starting without delay for 10 s, the tension continuously rises to the maximum value within 2~3.5 s, and then oscillates down and tends to be stable. When the starting time is 10 s and the delay is 10 s, the tension begins to decrease gradually when it reaches 80 kN in about 2 s. After the delay is over, it rises slowly again until it reaches the maximum value, and then oscillates down and tends to be stable. No matter what starting mode, the maximum tension is about 95 kN after pre-starting, which is less than 105 kN without delay starting. It can be seen from Figure 16 that under the acceleration time of 30 s, the acceleration extremum under various starting modes is reduced by more than half. However, an acceleration time that is too long is not conducive to the improvement of transportation efficiency. Pre-start has little effect on acceleration, but can significantly reduce the extreme tension of a conveyor belt. The step function curve proposed in this paper obtains relatively low values, both in acceleration and tension. It can be seen from Figure 16 that under the acceleration time of 30 s, the acceleration extremum under various starting modes is reduced by more than half. However, an acceleration time that is too long is not conducive to the improvement of transportation efficiency. Pre-start has little effect on acceleration, but can significantly reduce the extreme tension of a conveyor belt. The step function curve proposed in this paper obtains relatively low values, both in acceleration and tension. From the discussion of parametric response of acceleration and tension, it can be concluded that: large belt conveyor should have enough long starting time to ensure small enough starting acceleration. The reasonable starting curve should meet the characteristics of minimum starting acceleration extreme value, gentle change of acceleration curve and no sudden change. It is an effective way to reduce tension extremum and fluctuation by adding delay sections to the start curve of a large belt conveyor and adopting step function curve starting.

Conclusions
The problems associated with the dynamic modelling of belt systems with tree system topological configurations is by the classical rigid finite body element method. The number of generalized coordinates of an open-loop system is equal to the degree of freedom of the system, so the number of system equations is the smallest. The results of dynamic simulation analysis were verified by field test. It proves the correctness of the rigid-flexible multi-body dynamic model and analysis method of the belt conveyor established in this paper.
Acceleration and tension are important evaluation indices of conveyor running stability and safety. How to get a smaller acceleration and tension on the basis of ensuring the transportation efficiency is the design optimization goal of the belt conveyor. Compared with the previous research, the method in this paper is more systematic and universal, and can easily apply belt attributes, idler From the discussion of parametric response of acceleration and tension, it can be concluded that: large belt conveyor should have enough long starting time to ensure small enough starting acceleration. The reasonable starting curve should meet the characteristics of minimum starting acceleration extreme value, gentle change of acceleration curve and no sudden change. It is an effective way to reduce tension extremum and fluctuation by adding delay sections to the start curve of a large belt conveyor and adopting step function curve starting.

Conclusions
The problems associated with the dynamic modelling of belt systems with tree system topological configurations is by the classical rigid finite body element method. The number of generalized coordinates of an open-loop system is equal to the degree of freedom of the system, so the number of system equations is the smallest. The results of dynamic simulation analysis were verified by field test.
It proves the correctness of the rigid-flexible multi-body dynamic model and analysis method of the belt conveyor established in this paper.
Acceleration and tension are important evaluation indices of conveyor running stability and safety. How to get a smaller acceleration and tension on the basis of ensuring the transportation efficiency is the design optimization goal of the belt conveyor. Compared with the previous research, the method in this paper is more systematic and universal, and can easily apply belt attributes, idler attributes, collision attributes and preload attributes. It provides a theoretical basis for the optimization design of conveyor.