Multi-scale Simulations of Dry Friction Using Network Simulation Method

The AFM, SFM, FFM tip physical scheme on a sample surface model, Fig. 1, assumes one rigid sliding body connected by three springs (one for each spatial direction) to the mobile microscope support. Each spring reflects the elastic interaction between the two surfaces during the contact [36, 38]. The sliding body involves inertial effect and the interaction with the sample surface.

Figure 1

The hypotheses currently assumed in FFM, AFM and SFM models are: (i) the microscope support is rigid and moves at a constant velocity, v_M, (ii) the cantilever and the tip show elastic behaviour, and (iii) interactional force between the tip and the sample surface is a function of the relative distances between the tip and the atoms forming the sample surface.

Several forces act on the tip in x-direction: inertial force represented by m · (d²x_t/dt²), where x_t is the absolute displacement of the tip of mass m; elastic force from the spring, represented by the term k_x · (x_M − x_t); and damping force represented by the term c_x · (dx_t/d_t). The coefficients of these forces are constants. Similar forces act in the other directions. In addition, the interaction force between the tip and the sample surface is implemented from a surface potential [36, 38]. Thus, the surface potential for NaF with an FFM tip, as used by Hölscher et al. [36] is:

V(xt,yt)=−V0⋅cos(2πax⋅xt)⋅cos(2πay⋅yt)$$\begin{array}{} \displaystyle V({x_t},{y_t}) = - {V_0} \cdot \cos \left( {\frac{{2\pi }}{{{a_x}}} \cdot {x_t}} \right) \cdot \cos \left( {\frac{{2\pi }}{{{a_y}}} \cdot {y_t}} \right) \end{array}$$(1)

where the potential amplitude, V₀, is 1eV , while a_x and a_y, the double of the lattice constant, are 4.62 for both parameters. Hölscher et al. [37] use a slight modification of Eq. 1 to study HOPG with a SFM tip:

VHOPG(xt,yt)=−V0⋅[2⋅cos(2πa⋅xt)⋅cos(2πa3⋅yt)+cos(4πa3⋅yt)]$$\begin{array}{} \displaystyle {V_{HOPG}}({x_t},{y_t}) = - {V_0} \cdot \left[ {2 \cdot \cos \left( {\frac{{2\pi }}{a} \cdot {x_t}} \right) \cdot \cos \left( {\frac{{2\pi }}{{a\sqrt 3 }} \cdot {y_t}} \right) + \cos \left( {\frac{{4\pi }}{{a\sqrt 3 }} \cdot {y_t}} \right)} \right] \end{array}$$(2)

where V₀ is 0.5eV and the lattice constant, a, is 2.46Å. Sasaki et al. [38] used the Lennard-Jones potential to study the graphite with an AFM tip:

VTS=Σt4ε⋅[(σr0i)12−(σr0i)6]$$\begin{array}{} \displaystyle {V_{TS}} = {\Sigma _t}4\varepsilon \cdot \left[ {{{\left( {\frac{\sigma }{{{r_{0i}}}}} \right)}^{12}} - {{\left( {\frac{\sigma }{{{r_{0i}}}}} \right)}^6}} \right] \end{array}$$(3)

where the strength of the potential, ε, is 0.87381 · 10⁻²eV , and the value of the distance between the tip and the surface atom at which the potential is zero, σ, is 2.4945Å. The distance between the tip and the i − th atom of the surface, r_0i, is calculated from the i − th atom position inside the lattice, which is defined by its constant, 2.46 Å, and the nearest distance between 2 carbon atoms, 1.424Å. Hölscher et al. used the same potential with the same sample surface and microscope [39], but σ is equal to 3.4Å.

As a result of the balance of the forces on the tip, Fig. 1, the following equations can be written:

{m⋅x¨t=kx⋅(xM−xt)−∂V∂xt−cx⋅x˙tm⋅y¨t=ky⋅(yM−yt)−∂V∂yt−cy⋅y˙tm⋅z¨t=kz⋅(zM−zt)−∂V∂zt−cz⋅z˙t$$\begin{array}{} \displaystyle \left\{ {\begin{array}{*{20}{c}} {m \cdot {{\ddot x}_t} = {k_x} \cdot \left( {{x_M} - {x_t}} \right) - \frac{{\partial V}}{{\partial {x_t}}} - {c_x} \cdot {{\dot x}_t}}\\ {m \cdot {{\ddot y}_t} = {k_y} \cdot \left( {{y_M} - {y_t}} \right) - \frac{{\partial V}}{{\partial {y_t}}} - {c_y} \cdot {{\dot y}_t}}\\ {m \cdot {{\ddot z}_t} = {k_z} \cdot \left( {{z_M} - {z_t}} \right) - \frac{{\partial V}}{{\partial {z_t}}} - {c_z} \cdot {{\dot z}_t}} \end{array}} \right. \end{array}$$(4)

where x_M, y_M and z_M are the microscope support displacements over the three coordinates. The movement of the microscope support is a set of trips in the x-direction. For each of these trips, the initial condition of tip displacement in the x-direction is null, and in the z-direction is null except for the AFM tip on graphite sample. In the y-direction, the initial condition for this displacement is the same as for the microscope support. The initial condition of tip velocity is null for each direction.

Several forces act on the block: inertial force, represented by m⋅d2xdt2$\begin{array}{} \displaystyle m \cdot \frac{{{d^2}x}}{{d{t^2}}} \end{array}$, where x is the absolute horizontal displacement of the block of mass m, and the elastic forces from the two springs, represented by the terms −k₁ · x and −k₂ · (x − y). Similar forces act on the crank: inertial force, represented by J⋅d2ydt2/L$\begin{array}{} \displaystyle J \cdot \frac{{{d^2}y}}{{d{t^2}}}/L \end{array}$, where y is the absolute vertical displacement of the horizontal arm tip of the crank with inertial moment of J, and the elastic forces from the two springs, represented by the terms −k₃ · y · L and −k₂ · (y − x) · L. The coefficients of these forces are constant. In addition, the interaction force between the block and the belt is implemented by [76, 77].

{|Ff|≤μS⋅FN=Fs;ifvrel=0Ff=−sgn(vrel)⋅μ⋅FN;ifvrel≠0$$\begin{array}{} \displaystyle \left\{ {\begin{array}{*{20}{c}} {\left| {{F_f}} \right| \le {\mu _S} \cdot {F_N} = {F_s};\;{\rm{if}}\;{{\rm{v}}_{{\rm{rel}}}} = {\rm{0}}}\\ {{F_f} = - sgn({v_{rel}}) \cdot \mu \cdot {F_N};\;{\rm{if}}\;{{\rm{v}}_{{\rm{rel}}}} \ne {\rm{0}}} \end{array}} \right. \end{array}$$(5)

In Eq. 5, the first equation represents the stick phase, when the block has the same velocity that the belt. In this case, it is not necessary to solve the dynamic equation to get information about the block velocity. The second equation represents the slip phase, with the block slipping, and the dynamic equation is used in order to obtain the acceleration. The friction force is governed by Coulomb’s Law, which involves step function.

From Eq. 5, it is possible to distinguish two physical phases: stick and slip. Stick is controlled by the following single equation

x˙=vdr$$\begin{array}{} \displaystyle \dot x = {v_{dr}} \end{array}$$(6)

whereas slip phase is controlled by a more complex equation. Thus, as a result of the balance of the forces on the block and on the crank, the following equations can be written

{m⋅x¨+k1⋅x+k2⋅(x−y)=Fk(vrel)J⋅y¨L+k3⋅y⋅L−k2⋅(x−y)⋅L=0$$\begin{array}{} \displaystyle \left\{ {\begin{array}{*{20}{c}} {m \cdot \ddot x + {k_1} \cdot x + {k_2} \cdot (x - y) = {F_k}({v_{rel}})}\\ {J \cdot \frac{{\ddot y}}{L} + {k_3} \cdot y \cdot L - {k_2} \cdot (x - y) \cdot L = 0} \end{array}} \right. \end{array}$$(7)

where F_k is a more complex version of F_f, which involves the relative velocity effect.

As other authors, we also consider k₂ = k₁ [56, 63–65]

The parametrized Eq. 7 is given by

{mk1⋅d2xdt2+2x−y=−μ⋅(−k3k1⋅y+mk1⋅g+mck1⋅g)⋅2π⋅arctan(ε⋅vrel)1+δ⋅|vrel|Jk1⋅L2⋅d2ydt2−x+k1+k3k1⋅y=0$$\begin{array}{} \displaystyle \left\{ {\begin{array}{*{20}{c}} {\frac{m}{{{k_1}}} \cdot \frac{{{d^2}x}}{{d{t^2}}} + 2x - y = - \mu \cdot \left( { - \frac{{{k_3}}}{{{k_1}}} \cdot y + \frac{m}{{{k_1}}} \cdot g + \frac{{{m_c}}}{{{k_1}}} \cdot g} \right) \cdot \frac{2}{\pi } \cdot \frac{{{\rm{arctan}}\left( {\varepsilon \cdot {v_{rel}}} \right)\;}}{{1 + \delta \cdot \left| {{v_{rel}}} \right|}}}\\ {\frac{J}{{{k_1} \cdot {L^2}}} \cdot \frac{{{d^2}y}}{{d{t^2}}} - x + \frac{{{k_1} + {k_3}}}{{{k_1}}} \cdot y = 0} \end{array}} \right. \end{array}$$(8)

In the Eq. 8, ε · v_rel is a nondimensional magnitude, and so ε has the dimension s/m.

The last term of the first equation of Eq. 8 represents F_k/k₁. The friction force magnitude is calculated from the normal reaction, an addition of weights and a spring force. The sign function is assessed by arctangent.

To simplify Eq. 8, it is necessary to define a nondimensional time, τ=k1/m$\begin{array}{} \displaystyle \tau = \sqrt {{k_1}/m} \end{array}$. Thus, the Eq. 8 can be written as

{d2xdτ2+2x−y=−μ⋅ (−k3k1⋅y+g')⋅2π⋅ arctan(ε'⋅v′rel)1+δ′⋅|v′rel|Jm⋅L2⋅d2ydτ2−x+k1+k3k1⋅y=0$$\begin{array}{} \displaystyle \left\{ \begin{array}{*{20}{c}} {\frac{{{d^2}x}}{{d{\tau ^2}}} + 2x - y = - \mu \cdot ( - \frac{{{k_3}}}{{{k_1}}} \cdot y + g') \cdot \frac{2}{\pi } \cdot \frac{{\arctan (\varepsilon ' \cdot {{v'}_{rel}})\,}}{{1 + \delta ' \cdot |{{v'}_{rel}}|}}} \\ {\frac{J}{{m \cdot {L^2}}} \cdot \frac{{{d^2}y}}{{d{\tau ^2}}} - x + \frac{{{k_1} + {k_3}}}{{{k_1}}} \cdot y = 0} \\ \end{array}\right. \end{array}$$(9)

where J/m · L² is the parameter ξ, and k₁ + k₃/k₁ is the nondimensional parameter κ. Eq. 9 can be written as

{d2xdτ2+2x−y=−μ⋅((1−κ)⋅y+g′)⋅2π⋅ arctan(ε'⋅v′rel)1+δ′⋅|v′rel|ξ⋅d2ydτ2−x+κ⋅y=0$$\begin{array}{} \displaystyle \{ \begin{array}{*{20}{c}} {\frac{{{d^2}x}}{{d{\tau ^2}}} + 2x - y = - \mu \cdot ((1 - \kappa ) \cdot y + g\prime ) \cdot \frac{2}{\pi } \cdot \frac{{\arctan (\varepsilon ' \cdot {v^\prime }_{rel})\,}}{{1 + \delta \prime \cdot |{v^\prime }_{rel}|}}} \\ {\xi \cdot \frac{{{d^2}y}}{{d{\tau ^2}}} - x + \kappa \cdot y = 0} \\ \end{array} \end{array}$$(10)

The initial horizontal displacement of the block is null, as is the vertical displacement of the horizontal arm tip of the crank. The initial velocity of the block is the belt velocity and the vertical velocity of the horizontal arm tip of the crank is null.

The initial conditions related to displacements, x_M − x_t = y_M − y_t = z_M − z_t = 0, and velocities, dx_t/dt = dy_t/dt = dz_t/dt = 0, are introduced in the specifications of the initial conditions of the capacitors (initial voltage value) and of the coils (initial current value), respectively. The whole network model is now run in a Duo Core, using the code PSpice [68,69] and very few rules are required since the model contains very few electrical devices.

The design of a reliable network model requires a formal equivalence between the equations of the model and those of the process, including the initial conditions. Since Eq. 4, the basic network model is related to this equation, to which initial conditions must be added. Basic rules for designing the model can be found in González-Fernández [1]. The first step is to choose the equivalence between mechanical and electric variables (different choices give rise to different networks). For the problem under consideration, the following equivalence is established: x_t (tip displacement in the x-direction) ≡ q(electric charge in the network), or dx_t/dt(tip velocity in the x-direction) ≡ dq/dt = i (electric current at the network). In the other directions, similar equivalences are applied.

Now, each term of the Eq. 4 is considered as a voltage in an electric component whose constitutive equation is defined by the mathematical expression of such a term. Hence, the equation itself can be considered as a voltage balance (second Kirchhoff’s law) of a network that contains as many electrical components connected in series as there are terms in the equation. The linear term of Eq. 4, d²x_t/dt², is associated to a coil (inductance) since the constitutive equation of this component is V_L = L(di/dt) = L(d²q/dt²), with L = 1H. The same electrical element is used for the other two directions. The non-linear terms cannot be implemented directly, and require the use of controlled current (or voltage) sources or auxiliary circuits. The first are devices whose output is defined by software as an arbitrary, continuous mathematical function of the dependent variables: voltages at any node of the network or the current in any component. This capability, which is beyond the scope of the classical electrical analogy that appears in most text books, makes the network method an interesting and efficient tool in the field of numerical computation.

To implement the term x_t, the value of this variable must be determined by integrating the current in the network model, dx_t/dt. To this end, an auxiliary circuit is implemented, Fig. 4(b), in the models of FFM on the NaF sample; Fig. 5(c) in the model of AFM on a graphite sample. The terms y_t and z_t are implemented in similar network models. F_1x is the controlled current source whose output is the current, dx_t/dt, that crosses the voltage generator (used as an ammeter), V_x. The voltage of this generator is zero so as not to disturb the network model of 4. The current of F_1x is integrated using a capacitor with C_x1 = C_y1 = C_x = 1F; in this way, the voltage through this capacitor,V_Cx1 = V_Cy1 = V_Cx = C_x−1 ∫ (dx_t/dt)dt, is simply the variable x_t (a resistor of very high value, R_INF, is needed for electrical continuity).

Figure 4

Figure 5

Once the variable x_t has been determined, a controlled voltage source is used to implement the term k_x · (x_M − x_t) in the model E₂ in the models of FFM on NaF sample. For the models of AFM on graphite, the controlled voltage source is E_Tx1. The variable x_M is obtained as the product of the microscope support velocity and time, and is provided by the auxiliary circuit shown in Fig.4(c). The term related to ∂V∂xt$\begin{array}{} \displaystyle \frac{{\partial V}}{{\partial {x_t}}} \end{array}$ is implemented by a pair of controlled voltage sources, E₁ and E_1a, for FFM on the NaF sample. For the models of AFM on graphite sample, the controlled voltage source is E_Tx2. When the potential used is a Lennard-Jones potential, an auxiliary circuit is implemented, Fig. 5(b), for the model of AFM on graphite sample. The terms related to ∂V∂yt$\begin{array}{} \displaystyle \frac{{\partial V}}{{\partial {y_t}}} \end{array}$ and ∂V∂zt$\begin{array}{} \displaystyle \frac{{\partial V}}{{\partial {z_t}}} \end{array}$ are implemented in a similar way. The last term of the first equation in Eq. 4 is implemented by a resistor with a resistance of c_x in the models of FFM on NaF sample. The last term of the other equations in Eq. 4 is implemented in similar form. The initial conditions related to displacements, x_M − x_t = y_M − y_t = z_M − z_t = 0, and velocities, dx_t/dt = dy_t/dt = dz_t/dt = 0, are introduced in the specifications of the initial conditions of the capacitors (initial voltage value) and of the coils (initial current value), respectively. The whole network model is now run in a Duo Core, using the code PSpice [68, 69] and very few rules are required since the model contains very few electrical devices.

The design of a reliable network model requires the formal equivalence between the equations of the model and those of the process, including the initial conditions. The basic network model is related to Eq. 10 to which initial conditions must be added. Basic rules for designing the model can be found in González-Fernández [1]. The first step is to choose the equivalence between mechanical and electrical variables (different choices give rise to different networks). For the problem under consideration, the following equivalences are established: x (horizontal displacement of the block) ≡ q₁ (electric charge in the network), y (vertical displacement of the horizontal arm tip of the crank) ≡ q₂ (electric charge in the network), or dx/dτ (velocity of the block) ≡ dq₁/dt = i₁ (electric current at the network), dy/dτ (vertical velocity of the horizontal arm tip of the crank) ≡ dq₂/dt = i₂ (electric current at the network).

Now, each term of Eq. 10 is considered as a voltage in an electric component whose constitutive equation is defined by the mathematical expression of such a term. Hence, the equation itself can be considered as a voltage balance (Kirchhoff’s second law) of a network that contains as many electrical components connected in series as there are terms in the equation. The linear terms of Eq. 10, d²x/dτ² and d²y/dτ², are associated to two coils (inductances) since the constitutive equation of this component is V_L = L(di/dt) = L(d²q/dt²), with L = 1H. Thus, the coils used in the network model have inductances of L₁ = 1H and L_1A = ξ = J/mL². The non-linear terms cannot be implemented directly, and require the use of either controlled current (or voltage) sources or auxiliary circuits. The first are devices whose output is defined by software as an arbitrary, continuous mathematical function of the dependent variables: voltages at any node of the network or currents in any component. This capability, which is beyond the scope of the classical electrical analogy that appears in most text books, makes the network method an interesting and efficient tool in the field of numerical computation.

To implement the terms x and y, the values of these variables must be determined by integrating the currents, dx/dτ and dy/dτ, in the network model. To this aim, two auxiliary circuits are implemented, Fig. 6(b) and Fig. 6(h). F₁ and F_1A are controlled current sources whose outputs are the currents, dx/dτ and dy/dτ, which cross the voltage generators (used as ammeters), V₁ and V_1A respectively. The voltage of these generators is zero so as not to disturb the network model of Eq. 10. The current of F₁ is integrated using a capacitor with C₁ = 1F; in this way, the voltage through this capacitor, VC1=C1−1∫(dx/dτ)dτ$\begin{array}{} \displaystyle {V_{C1}} = C_1^{ - 1}\int {(dx/d\tau )} d\tau \end{array}$, is simply the variable x (a resistor of very high value, R_INF1, is needed as a requirement of electrical continuity). In the same way, the current of F_1A is integrated using a capacitor with C_1A = 1F, the voltage through which is the variable y (another resistor of very high value, R_INF1A, is needed as a requirement of electrical continuity).

Figure 6

Once these variables have been determined, a controlled voltage source is used to implement the term 2 · x − y for the model; this is E₁. In the same way, a controlled voltage source is used to implement the term −x + κ · y for the model; this is E_1A. The last term of the first equation in the system described in Eq. 10, μ⋅((1−κ)⋅y+g')⋅2/π⋅arctan(ε′⋅vrel′)/(1+δ′·|vrel′|)$\begin{array}{} \displaystyle \mu \cdot ((1 - \kappa ) \cdot y + g{\rm{'}}) \cdot 2/\pi \cdot \arctan (\varepsilon \prime \cdot v_{rel}^\prime )/(1 + \delta \prime \cdot\left| {v_{rel}^\prime } \right|) \end{array}$, is also implemented by a controlled voltage source, E₂, and requires the variable dx/dτ. This is a problem because this variable is also the current that crosses this controlled voltage source. To solve this problem an auxiliary circuit represented in Fig. 6(f) provides the voltage used in E₂. To help this auxiliary circuit, another one is needed, Fig. 6(e), which calculates 2/π⋅arctan(ε′⋅vrel′)/(1+δ′·|vrel′|)$\begin{array}{} \displaystyle 2/\pi \cdot \arctan (\varepsilon \prime \cdot v_{rel}^\prime )/(1 + \delta \prime \cdot\left| {v_{rel}^\prime } \right|) \end{array}$.

The conditions that select the stick or the slip phase are implemented by two voltage controlled switches, S₁ and S₂, following the flowchart represented in Fig. 10. To implement this, it is necessary to obtain the voltages that appear in it by using the auxiliary circuits represented in Figs. 6(c), 6(d) and 6(e). The voltage provided by controlled voltage source E_control, in Fig. 6(c), is a scalar associated to the slip phase when the voltage of E_Elas, in Fig. 3(d), is bigger than the voltage of E_FRS, Fig. 6(e). Otherwise, there are two possibilities: if the current that crosses L₁ is not equal to the belt velocity, then the voltage of E_control will provide the same value of this scalar of the last case; if not, another check will be necessary. When the voltage in the coil is not equal to zero, then this controlled voltage source provides the same value of this scalar, which is the slip phase. If this is not the case, stick phase arises. The voltage of controlled voltage source E_Elas, in Fig. 6(d), provides the absolute value of the force in the horizontal springs, while controlled voltage source E_FRS, in Fig. 6(e), provides the maximum friction force.

Figure 10

Thus, if S₁ is closed (stick phase), then the current will be determined by the Eq. 6. On the other hand, if S₁ is open S₂ will be rapidly closed (slip phase) and the current will be determined by the Eq. 7.

The initial conditions of the problem are implemented by the initial conditions of the capacitors and coils. The full network model is run on a Core Duo processor, using PSpice code [69, 78, 79].

The Network Simulation Method has been applied to sample surfaces and tips for which data are available: FFM tip on NaF [36], SFM tip on HOPG [37], and AFM tip on graphite [39]. The first sample has a cubic crystal net and the others have a hexagonal crystal net, which tests the potential of the network model. The characteristics of the tip, the test conditions and the sample size for the different tests are:

the FFM tip on NaF sample is characterized by m_x = m_y = 10⁻⁸kg, c_x = c_y = 10⁻³N · s/m and k_x = k_y = 10N/m. The microscope support moves with a constant velocity of 400 Å/s, and the results are depicted with a scan range of 20x20Å and scan angle of 0 [36]. The proposed model solves the Eq. 4 using the potential defined in Eq. 1.

Eq. 4 is solved in each test for the potential functions defined in the following equations:

Eq. 1 in the FFM tip on NaF sample model

In the last case, for graphite sample, only short computational times are necessary for the software to find the right values of the parameters not supplied by Sasaki et al. [38]. Fig. 7(b) shows the components of the tip elastic force, F_x and F_y, in each position on the NaF sample surface, whose crystal net is shown in Fig. 7(a). The correspondence between both depicted images enables us to use this software as a tool to interpret those images from the microscope. Thus, this software is able to manage different theoretical crystal nets, even with crystallographic defects, providing the corresponding theoretical elastic forces, which can be compared with those from the microscope. The results of the simulation agree with those obtained by Hölscher et al. [36].

Figure 7

Fig. 8(b) shows the components of the tip elastic force, F_x and F_y, in each position on the HOPG sample surface, whose crystal net is shown in Fig.8(a). The results of the simulation agree with those obtained by Hölscher et al. [37]. Besides, Fig. 9(b) shows the x-direction component of the tip elastic force, F_x, in each position on the graphite sample surface, whose crystal net is shown in Fig. 9(a). In this case, a 3D graph is used to show the details better, to be more precise the stick - slip mechanism. The lines of the graph are traced at different y-coordinates, separated by a fixed step. When one of these lines crosses the sample surface at a critical position, the instabilities appear, Fig. 9(b). These critical positions have a relation with the minimal threshold of the distance between the tip and the atom on the surface. The results of the simulation agree with those obtained by Sasaki et al. [38].

Figure 8

Figure 9

The model used in this numerical simulation follows the mechanism shown in Fig. 11, and is similar to that used by Awrejcewicz et al [56]. Its dimensions and properties are:

crank arm length, L = 0.105m

Figure 11

Let us assumed that the inertial moment of the crank, J, is equal to m_c · L²/3. The dimensionless parameters, from the data above, are: ξ = 1/3 · m_c/m = 0.412 and κ = 1.5. The value of the modified gravity acceleration is: g^′ = g · m/k₁ = 0.01778m.

The parameters that characterize the friction are:

dimensionless accuracy parameter, ε^′ = 3.5121e7

The static friction coefficient, μ_s, is calculated from a measure using a mass on the belt. A friction force of 1.6 N is obtained from a mass of 0.277kg, which means a value of μ_s equal to 0.59. The modified belt velocity used, v^′_dr, is equal to -0.002788 m.

Fig. 12 shows the phase plane for the belt velocity used for a time span from 100 to 1,000. Stick phase can be recognized as a horizontal line with ordinate V₁ = −0.002788m, Fig. 12. In addition, Fig. 13 shows the frequency spectrum of the signal X₁ for a time span from 100 to 1,000.

Figure 12

Figure 13

A waveband arises centred on the value of 0.143433 Hz, a value that has no relation with belt velocity. Fig. 14 shows the Poincaré map using a period of 1/0.143433 for a time span from 0 to 3,200.

Figure 14

Fig. 14 shows the chaotic behaviour of the system. To quantify the chaos, an approximation to the maximum Lyapunov exponent is calculated from a time series of the solution by the software OpenTSTool [80], developed by Göttingen University. This approximate value is equal to 0.0514, which suggests chaotic behavior of the solution.