Next Article in Journal
Robotic Manipulator’s Expressive Movements Control Using Kinematic Redundancy
Next Article in Special Issue
An Adaptive Control Method for the Distribution Valve of a Digital Pump
Previous Article in Journal
Estimation of Unmeasurable Vibration of a Rotating Machine Using Kalman Filter
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Multibody Modeling of a New Wheel/Track Reconfigurable Locomotion System for a Small Farming Vehicle

1
Dipartimento di Ingegneria dell’Innovazione, Università del Salento, 73100 Lecce, Italy
2
Graduate School of Life Science and Systems Engineering, Kyushu Institute of Technology, Kitakyushu 808-196, Japan
*
Author to whom correspondence should be addressed.
Machines 2022, 10(12), 1117; https://doi.org/10.3390/machines10121117
Submission received: 24 October 2022 / Revised: 18 November 2022 / Accepted: 21 November 2022 / Published: 24 November 2022
(This article belongs to the Special Issue Mechatronic Systems: Developments and Applications)

Abstract

:
Tracks and wheels are the most widely used running gear for the locomotion of agricultural vehicles. The main difference between the two systems is the contact area with the ground and, consequently, the pressure distribution. Evaluating the pressure distribution on the ground is important because soil damage and vehicle performance depend on it. This analysis is especially difficult for tracked vehicles, owing to their complexity compared with wheeled systems. In this paper, we describe a multibody model of a flexible track to evaluate the pressure distribution upon contact with soft terrain. The track considered in this study is part of a reconfigurable locomotion system of a small farming vehicle, which can vary the pressure distribution by switching from a wheeled vehicle to a half-tracked vehicle. The aim of such a vehicle is to minimize soil damage in addition to optimizing its performance. The model is used to characterize this vehicle and evaluate the pressure distribution with varying characteristic parameters, such as track tension, the position of the vehicle’s center of gravity, the weight distribution on the track itself, and the stiffness of the suspension system.

1. Introduction

Tracked systems are among the most common running gear used for various off-road vehicles, from military applications to the agricultural field. Tracked vehicles can be divided into two main categories: rigid and flexible tracks [1]. In the first case, the track belt structure consists of a solid chain track made of steel plates, whereas in the second case, the belt is made of rubber reinforced with steel wires along the longitudinal length. The track belt is forced onto the ground by a series of wheels called rollers, and it is moved by an additional motored wheel called a sprocket, the position of which depends on the system layout. Although the two types of tracks (rigid and flexible) lead to differing performances of the vehicles, a common aspect is the capability to deform upon contact with the soil and adapt to it to some extent. An optimal design of these systems requires taking into account this flexibility. This makes the analysis of tracked vehicles challenging compared with wheeled vehicles. The vehicle’s performance depends on the forces exerted on the ground as a function of the normal and shear stress at the contact patch. Accordingly, it is important to have a reliable tool to quantify the stress distribution at the point of contact between the running gear and soil in order to analyze and design tracked vehicles.
Several methods have been proposed to carry out such analysis. The first adopted strategy uses empirical methods [2], whereby simple soil measurements and some main vehicle characteristics are used to assess the performance. Another strategy involves using analytical models, among which the method developed by Bekker is the most famous [2]. These strategies have the drawback of relying on simplified assumptions, and in the case of analytical models, the pressure distribution at the contact point with the ground is assumed beforehand. Common pressure distributions are constant, trapezoidal, and multipeak sinusoidal. In this last case, the pressure peaks are related to the rollers. The real pressure distribution is unlikely to be of the functions mentioned above. It is affected by the position of the vehicle’s center of gravity, the tension of the track belt, the geometry, the disposition of the rollers, and the suspension system forcing the rollers onto the belt. Whereas parametric models of increasing complexity have been proposed to analyze tracked systems [3,4,5,6], the high computational capacity of modern computers enabled the development of numerical methods to quantify the vehicle performance. An example is NTVPM [2], which can model the interaction between the track and soil in detail.
Multibody systems have proven to be promising tools for the modeling of tracks. In [7], a multibody model was developed, assuming that the track can be approximated as a set of imaginary wheels. In [8], a multibody model was developed, employing a nonlinear finite-element representation that can simulate both flexible rubber tracks and rigid tracks.
In this work, we propose a step-by-step implementation of a multibody model to study the pressure distribution of a reconfigurable wheeled/tracked vehicle on soft terrain. The vehicle consists of a rear-wheel axle and a wheel/track reconfigurable front axle, which can adjust the contact area with the ground. A tracked configuration can better distribute pressure and causes less compaction than a wheel. On the other hand, tracks cause more soil distortion [9] and are less efficient during steering maneuvers than wheels. The vehicle is intended to adjust its contact area based on the soil conditions to minimize soil damage while considering the rolling resistance and traction performance. To design such a system, it is important to quantify the pressure distribution at the contact point with the soil; this is the purpose of the multibody model developed in this work. Although a similar task could be accomplished by relying on dedicated plugins of commercial software [10], adopting a more diffuse and general-purpose tool such as MATLAB® enables the development of a model from scratch, facilitating a deeper understanding and a better compromise between accuracy and computational effort. Moreover, a step-by-step implementation permits customization of the model based on future research needs.
In this paper, we describe the strategy pursued in modeling the system consisting of half the robot frame, the reconfigurable front axle, the rear wheeled axle, and the ground.

2. Materials and Methods

A 3D model of the system in the simulation environment is shown in Figure 1. A detailed model of the tracks and the reconfiguration mechanism was combined with a simplified body representing the chassis and the rear-axle wheel.
A functional sketch of the model is shown in Figure 2, in which the main components are represented:
Ground, indicated with the letter G, represents the soil;
Frame, indicated with the letter F, represents the vehicle chassis;
Track module, indicated with the letter T, is a set of bodies assembled to form the system shown in Figure 3. The assembly in Figure 4 is hinged to the vehicle chassis (F) through a revolute joint located at point T₁;
Shock absorber, indicated with the letter S (blue line). According to one of the requirements discussed above, this element is inserted to endow the mechanism with the capability of letting the track better adapt to bumpy, uneven ground. The bottom end of the shock absorber is hinged to the track frame (T) through a revolute joint;
Rocker, indicated with the letter R, represents the key element of the switching mechanism, deputed to transfer the motion from the actuator to the track frame. It is a ternary link that is often used in passive [11] and active adjustable suspension systems [12,13]. It is hinged to the vehicle chassis (F) through a revolute joint at point R₁. The upper end of the shock absorber (S) is connected to the rocker through a revolute joint located at point R₂;
Actuator, indicated with the letter A (red line), is the element that drives the mechanism. One end of the actuator is hinged to the frame (F) through a revolute joint located at point A₁, whereas the other end is linked to the rocker (R) through a revolute joint at point R₃; and
Wheel, indicated with the letter W, represents the wheel of the robot’s rear axle.
The track module is a set of bodies assembled together. This module comprises the sprocket; the idler; the track frame; the undercarriage; and the chain, which constitutes the track belt. In order to capture the functioning of this module, the aforementioned elements, especially the links making up the chain, cannot be rigidly connected to one another, as they experience multiple contacts, the implementation of which is detailed in Section 2. The shock absorber is composed of two main bodies—the cylinder and the piston—connected by a prismatic joint. The multibody model of the actuator is composed of three blocks: the main body, the leadscrew, and the plunger (also called piston or rod. The relative motion between the actuator body and the plunger is dictated by a prismatic joint, allowing only translation along the axis of the leadscrew; the position along the longitudinal axis of the cylinder varies from 0 mm when the actuator is fully retracted to 100 mm (equal to the stroke) when the actuator is fully extended. The wheel is modeled as a rigid body composed of two parts, the tire and rim, which are clearly visible on the right side of Figure 1.
The functioning of the proposed switching mechanism is extremely simple; when the actuator (A) retracts (distance between points A₁ and R₃ decreases), it causes the rocker (R) to rotate about point R₁ in a counterclockwise direction; consequently, the track module (T) is forced through the shock absorber (S) linked to the rocker (R) to rotate about point T₁ in a counterclockwise direction, thus lifting the idler section of the track from the ground. Conversely, if the actuator (A) extends (the distance between points A₁ and R₃ increases), the rotations of the rocker R and the track module T occur in a clockwise direction, and the reconfigurable section of the track is lowered to the ground. The motion imposed by the actuator is transferred to the track module in an effective yet simple way, thus providing for the reconfiguration of the system.

2.1. Track Module

The track module shown in Figure 3 was designed, realized, and tested at the Ishii Laboratory of the Kyushu Institute of Technology [14]. A roller chain track was chosen among the various possible mechanical solutions to build a custom track. The belt consists of an attachment chain (also known as a conveyor chain) with built-in flanges on which treads can be easily attached to form a track.
As shown in Figure 3, the roller chain track is made up of a chain engaging with a toothed drive sprocket on one side (hereafter referred to as a sprocket) and a toothed idler sprocket on the other (hereafter referred to as an idler). The sprocket moves the track along, and the idler keeps the track straight. The sprocket has 33 teeth and an outer diameter of 281 mm, whereas the idler has 15 teeth and an outer diameter of 135 mm, with pitch diameters of 267.21 mm and 122.17 mm, respectively.
On each one of the chain attachments, an aluminum U channel is fastened by two bolts, as shown in Figure 5.
The links of the track belt that gradually enter the lower branch, specifically in the section separating the sprocket from the idler, must be pressed against the ground to effectively redistribute the weight of the vehicle.
This important function is accomplished by the undercarriage shown in Figure 6. It is composed of two rows carrying two rollers each. The shafts about which the rollers are free to rotate are supported by two L–shaped pieces, which are connected on one end to the idler shaft and on the other end to the arm between the sprocket and idler shafts.
The arm connecting the sprocket and idler shaft is represented in Figure 7, hereafter referred to as the track frame. It is an ensemble of three telescopic squared aluminum frames: a central frame, a solid block, and a frame that can slide inside the other two, which are hollow and connected by a hexagonal nut. This nut acts as a tensioner for the track belt; one of the ball joints at the end of the nut is right-threaded, whereas the other has a left thread so that the rotation of the nut can increase or decrease the distance between the idler and sprocket, thus tensioning the chain.

2.1.1. Track Frame

The track frame shown in Figure 7 is the element that connects the sprocket shaft and the idler shaft. It is composed of two telescopic frames that can slide relative to each other; a tensioner nut with a right thread on one side and a left thread on the other side provides a simple way to tension the chain that constitutes the track belt. On the idler side of the track frame, the shock absorber has to be mounted to allow the switching mechanism raise and lower part of the track; therefore, as shown in Figure 8, the attachment bracket of the tensioning system is equipped with an additional hole to house the lower eyelet of the shock absorber. To endow the model with the possibility of specifying an arbitrary tension for the track chain, thus improving its fidelity, the two parts composing the track frame that rotate about the sprocket shaft and carry the idler shaft are connected by a prismatic joint with base (Ts1) and follower (Ts2) (see Figure 9). This joint may be assigned a null displacement or a certain force. In the first case, a null displacement between Ts1 and Ts2 is imposed, and the chain is loose, i.e., the tension is zero; in the second case, an arbitrary force (Fts) oriented from Ts1 to Ts2 is applied to the joint, and consequently, the track chain experiences a tension (Tch) as expressed by (1),
T c h = F t s 2 cos α
where α is the chain angle shown in Figure 9. In light of such implementation, it is possible to simulate the model with differing values of chain tension, which are usually expressed as fractions of the vehicle weight [1].

2.1.2. Track Belt

The track belt comprises a chain with attachments (see Figure 10). Although the chain is composed of articulated links forming a closed kinematic loop, the belt can be modeled by connecting all the links to form a unique rigid body made up of the chain, the sprocket, the idler, the undercarriage, and the track frame. Such a modeling strategy, although simple and fast, does not capture the flexibility of the track belt, thus leading to an unrealistic, almost uniform pressure distribution over the ground contact area, which is equivalent to assuming a rigid footing, as is often the case in the literature on rolling resistance evaluation. To fulfill the goal of assessing how the pressure originated by the track–terrain interaction distributes over the contact patch, it was necessary to model the chain more accurately.
The chain shows the alternation of external links (26 pin links such as that on the left in Figure 11a) and internal links (26 bushing links such as that on the right in Figure 11a). Each link is connected to the two adjacent links using revolute joints to permit relative rotations, as shown in Figure 11a. The total of 52 links are divided into 13 groups of 4. Figure 11a,b indicate the contact surfaces, highlighted in different colors:
  • In red, the surfaces of the bushings of the inner link that come into contact with the two gears of the track module;
  • In green, the portions of the upper surface of U channels that come in contact with the undercarriage rollers; and
  • In blue, the surface of the rubber lug that comes in contact with the ground.
  • The contacts mentioned above were analyzed and implemented in the proposed model.

2.2. Contact Model

As previously mentioned, the bodies of the multibody model, especially those comprising the track module shown in Figure 11, in addition to coming into contact with the ground (the rubber lugs), collide with each other, so a precise and accurate representation of impact or contact is needed [15,16,17]. The contact of the track links with the gears, and the track module’s rollers are modeled using a linear spring–damper model. According to this method, the local deformations of the contacting bodies in the contact region are simplified, i.e., allowing the bodies to penetrate each other to a limited extent. The amount of penetration, denoted by δ, is equal to the sum of local deformations undergone by the two contacting bodies. The penetration causes a pair of resistive contact forces to act on the two bodies in opposite directions.
The contact force can be decomposed into the normal and tangential components, which act along unit vectors n and t. The normal contact force that A exerts on B (FnAB) shown in Figure 12 is computed as if it were generated by a linear spring–damper system placed between points A and B so that [18]
F n A B = σ   δ B A · K · δ B A + γ δ ˙ B A D δ ˙ B A n A
where:
-
K is the contact spring stiffness and represents the resistance opposing the geometric penetration of the two colliding bodies. The larger the value of the spring stiffness, the harder the contact. The value K = 106 N/m was set based on an iterative tuning;
-
D is the damping coefficient representing the resistance of the contact damper opposed to motion when the bodies penetrate. The larger the damping coefficient value, the more energy is lost during the collision and the faster the contact vibrations are dampened. The value D = 103 Ns/m was chosen based on iterative tuning; and
-
δBA is the penetration depth of B into A, defined by (3):
δ B A = r b r a · n = r B A · n
where r A and r B are the position vectors of point A and point B, respectively; r B A is the position of B relative to A; and δ ˙ B A is the penetration velocity of B into A, as expressed by (4):
δ ˙ B A = v b v a · n = v B A · n
where v A and v B are the velocity vectors of points A and B, respectively; v B A is the velocity of B relative to A (i.e., as demonstrated by a reference frame attached to A);
-
σ(δBA) is the contact force activation factor, defined as in (5),
σ ( δ B A ) = 0                                                                                                         i f   δ B A 0 t a n h ( ln 199 2 · δ T R W · δ B A                                       i f δ B A < 0  
where δTRW is the transition region width, which specifies the region over which the spring–damper force increases to its full value. The hyperbolic tangent that appears in the definition of the function σ (δ) smoothly ramps up the force, whereas the penetration depth moves through the transition region (0 < δδTRW). When δ = δTRW, at the end of the transition region, the contact force reaches a value equal to 99% of the full stiffness and damping. Similarly, on the rebound, the elastic and viscous components of the force are smoothly decreased back to zero. Reducing the transition region width makes the onset of contact sharper and the model more accurate; however, the solver is forced to take smaller time steps to deal with a discontinuous problem, and the simulation speed decreases. The value δTRW = 0.1 mm was chosen, as it proved to be a good compromise between speed and accuracy.
-
γ δ ˙ B A is the damping-enabling factor, defined as in (6)
γ ( δ ˙ B A ) = 0                                                                                                         i f   δ ˙ B A 0   1                                                                                                         i f   δ ˙ B A < 0  
which deactivates the damping component of the contact force on the rebound.
-
nA is the unit normal vector pointing outwards from the surface of A at the point of contact.
The force that B exerts on A (FnBA) is equal and opposite to FnAB, i.e., FnBA = −FnAB. Owing to contact friction, the tangential component was neglected in the present analysis. It does not play a considerable role with respect to the switching phases and causes the solver to incur an increased computational burden. This is the reason why the subscript n is hereafter dropped from the symbol of the contact force.

2.2.1. Contact between Track Links and Gears

The alternation of the internal links (bushing links) and external links (pin links) shown in Figure 11a forms the roller chain track of the robot. It is evident that to simulate the switching phases, it is necessary to model the contact between the chain bushing/rollers and the teeth of both the sprocket and the idler [19,20,21]. First, an assumption was made about the motion of the links, i.e., that they are constrained through a planar joint to perform a planar motion on the plane of the two gears so that the model accounts for the contact between the bushing cylinder (shown in red in Figure 11) and the surface of the sprocket teeth, whereas any side contact is excluded. A reference frame is attached to the sprocket as follows:
  • The origin is placed at the center of the sprocket denoted by the symbol Os;
  • Axes xs and ys are configured as in Figure 13; and
  • The zs axis is perpendicular to the plane of the illustration in Figure 13.
The profile of the surface of each tooth of the sprocket is quite complex and would require detailed geometric features, such as an x–y spline curve, to be faithfully represented.
In the present analysis, to simplify the model and reduce the computation time, for each tooth of the sprocket, six points were selected relative to the sprocket center in such a way as to construct an approximated piecewise linear tooth profile sufficiently close to the real profile, as shown in Figure 14. The coordinates (xsi,1, ysi,1) expressed in the reference system (Os; xs, ys) of the point Pi,1 i = 1, …, 6 used to approximate the first tooth are summarized in Table 1. The six points identify five lines that represent as many surfaces of possible contact; the nth of these lines, with n = 1, …, 5, has a tangent unit vector given by (7).
t n , 1 = x n + 1 , i s x n , i s x n + 1 , i s x n , i s 2 + y n + 1 , i s y n , i s 2 · x s + y n + 1 , i s y n , i s x n + 1 , i s x n , i s 2 + y n + 1 , i s y n , i s 2 · y s
where xs and ys are the unit vectors of the xs and ys axes, respectively. The unit vector of the tooth profile can be obtained through a 90 degree clockwise rotation of tn, i.e., in (8).
n n , 1 = y n + 1 , i s y n , i s x n + 1 , i s x n , i s 2 + y n + 1 , i s y n , i s 2 · x s x n + 1 , i s x n , i s x n + 1 , i s x n , i s 2 + y n + 1 , i s y n , i s 2 · y s
The first tooth profile that is the set composed of the points Pi,1, i = 1 …, 6 and the vectors tn,1 and nn,1, n = 1, …, 5 is replicated around the zs axis 33 times according to the number of sprocket teeth (see Figure 13), yielding 6× (33 teeth) = 198 points, 5× (33 teeth) = 165 tangent vectors, and as many normal vectors. The ith point of the profile of the kth sprocket tooth, with i = 1, …, 6, and k = 1, …, 33, is denoted by the symbol Pi,k. Similarly, the vectors tangent and normal to the nth profile segment of the kth tooth, with n = 1, …, 5, are indicated by tn,1 and nn,1, respectively. Because the coordinates of the characteristic points of the tooth profile remain constant for the sprocket reference system, the tangent and normal vectors, given by Equations (7) and (8), respectively, are computed before the simulation and stored for later repeated calculations to reduce the running time. Figure 15 shows the geometry of the contact between a chain bushing and one of the sprocket teeth. Each inner link of the chain has two bushings (highlighted in red in Figure 11) that come into contact with the toothed surface of the sprocket; as shown in Figure 11b, a reference frame is attached to the single bushing, with the origin placed at the center point (F1) and axes xb, yb, and zb, with the latter perpendicular to the plane of the diagram.
The algorithm that calculates forces and torques originated by such interactions is articulated through the following steps:
1.
The distance between the sprocket center (Os) and the bushing center (F1) is calculated by (9),
d b = x F 1 s 2 + y F 1 s 2
where xsF1 and ysF1 are the x y coordinates of point F1 with respect to the reference system (Os; xs, ys, zs).
2.
A check is performed to verify whether the bushing in question can collide with the sprocket. Contact can take place if and only if the distance (db) in (9) is less than or equal to a prespecified proximity range, i.e., (10):
d b δ p r o x s = R s + P c h
The proximity range (δsprox) was set equal to the sum of the outer sprocket radius (Rs) and the pitch of the chain (Pch). Inequality (10) is a necessary but not sufficient condition for the sprocket and the bushing to come into contact; if it is verified, the algorithm continues through the next steps; otherwise, the algorithm does not continue. The check in question avoids needless calculations for links that are excessively distant from the sprocket.
3.
The sprocket tooth with which the considered bushing is likely to collide is identified based on the angle (ψ) between the vector from Os to F1 and the xs axis of the sprocket reference frame, i.e., (11)
Ψ = arctan ( y F 1 s x F 1 s )
Knowing the value of ψ, it is possible to select the index (k) of the tooth candidate for contact, along with relevant information, such the tangent vectors (tkn) and the previously determined normal vectors (nkn, n = 1, …, 5) (see the Equations (7) and (8)). Thanks to this step, the calculations exposed in the next points are performed only for the candidate tooth, thus preventing the algorithm from cycling over all the sprocket teeth.
4.
The positions of the profile points of the kth tooth with respect to the bushing center (F1) are computed as (12)
L i , k * = x F 1 s x i , k * s · x s + y F 1 s y i , k * s · y s ,                                             i = 1 , , 6
where xsi,k and ysi,k are the x and y coordinate, respectively, of the point (Pi,k) expressed in the sprocket reference system (Os; xs, ys).
5.
For each segment composing the k*th tooth profile, a check is performed to verify whether contact could occur with the circle representing the bushing. Contact can occur between the bushing circle and the profile line segment if the orthogonal projection of the center of the circle (F1) onto the line of the profile falls between the two ends of the segment (see Figure 15). Mathematically, this translates into the following inequality (13),
( t n , k * · L n , k * ) ( t n , k * · L n + 1 , k * ) 0                                             n = 1 , , 5
where represents dot product, and the vectors Ln,k are given by (12). The profile segments for which condition (13) is true are those probably in contact with the bushing, whereas the others attribute no contribution to the interaction in question. Here, for clarity and brevity of the explanation of the following steps, we assume that only the n*th line of the tooth profile satisfies the inequality (13); the code specifically written for the task at hand is capable of dealing with completely generic circumstances in which the bushing contacts more than one side of the tooth profile at the same time.
6.
The penetration of the bushing (b) into the sprocket (s) is defined in (14):
δ b s = d R b
where Rb is the bushing outer radius, and d is the distance from the center of the bushing (point F1) to the n*th line of the k*th tooth profile given by (15)
d = P n * + 1 , k * P n * , k * L n * , k * P n * + 1 , k * P n * , k *
where Pn*,k* = xsn*,k*xs + ysn*,k*ys and Pn*+1,k* = xs n*+1,k*xs + ysn*+1,k*ys are the position vectors of the n*th and the (n* + 1)th point of the sprocket k*th tooth profile, respectively, expressed in the reference frame (Os; xs, ys).
7.
The velocity with which the bushing penetrates into the sprocket is calculated as in (16):
δ ˙ b s = v b s d ω b s t n * , k * · n n * , k *
where vsb is the velocity vector of the bushing center (F1) as demonstrated by the sprocket reference system (Os; xs, ys); and ωsb is the angular velocity of the frame (F1, xb, yb) attached to the bushing, as demonstrated by the sprocket reference system (Os; xs, ys). The quantity in parentheses in Equation (16) represents the velocity of the point of contact (Pc) shown in Figure 15, the position vector of which in the coordinate system (Os; xs, ys) is given by (17).
P c = x F 1 s x s + y F 1 s y s d n n * , k *  
8.
The force exerted by the sprocket tooth on the chain bushing (Fsb) is computed (18) by substituting the penetration δbs from (14) and the penetration velocity δ ˙ b s from (16) into Equation (1)
F s b s = σ δ b s K δ b s + γ δ ˙ b s D δ ˙ b s n n * , k *
where the superscript s indicates that the force vector components calculated through Equation (18) are expressed in the reference system integral with the sprocket; therefore, to apply the force to the bushing, it is necessary to express its components in the coordinate system attached to the bushing itself. This can be achieved using the relative rotation matrix [R]sb, as in (19).
F s b s = R x x R x y R x z R y x R y y R y z R z x R z y R z z s b T F s b s · x s F s b s · y s F s b s · z s
Fbsb is then applied to the bushing.
9.
The force exerted by the chain bushing on the sprocket tooth is determined by simply inverting the sign in (19), i.e., Fsbs = −Fssb. The force (Fsbs) also causes a torque (Msbs) with respect to the sprocket center (Os) (arm denoted by an in Figure 15), given by (20).
M b s s = P c F b s s
Both the force (Fssb) and the moment (Msbs) are applied to the sprocket.
At each time step, the implemented model carries out the calculations related to the contact of the two gears of the track module with each of the 52 chain bushings. Much of the algorithm outlined above, except the first two points, is implemented through a specially written MATLAB® function. The contact of the chain links with the idler is dealt with using a procedure similar to that described for the sprocket.

2.2.2. Contact between Track Links and Rollers

The rollers placed in the undercarriage of the track module shown in Figure 2 are essential to evenly distribute the vehicle’s weight over the ground contact area. Therefore the contacts between the rollers and every link composing the chain must be modeled to obtain reliable results. Because the chain links are constrained to perform a planar motion in the plane of the two gears, the contact in question ultimately reduces to the problem of the interaction between a finite line (representing the upper surface of the U channel highlighted in red in Figure 11) and a circle (representing the outer diameter of the roller), as shown in Figure 16. A reference frame is placed on the U channel of each chain link, with the origin positioned at the midpoint (F3) of the contact surface, the xu axis and yu axis parallel and normal, respectively, to the contact surface, and the zu axis perpendicular to the plane of the illustration in Figure 16. A coordinate system is also associated with the roller, originating at the center point (Or) and axes xr, yr, and zr. The x and y coordinates of the roller center’s reference frame (F3; xu, yu, zu) are denoted by the symbols xuOr and yuOr, respectively. Because the total length of the contact line is 22 mm, three possible contact situations may arise:
1.
xuOr > 11 mm, in which case the projection of the roller center onto the line of contact falls “beyond” the second edge (E2), as represented in Figure 16b. Contact takes place at point E2, and the penetration occurs along the direction of the following unit vector (21),
n u = x O r u x E 2 u d x u + y O r u d y u
where xuE2 = 11 mm is the x coordinate of point E2 in the U-channel frame (F3; xu, yu, zu), and d is the distance from the edge point (E2) to the roller center (Or) given by (22).
d = x O r u x E 2 u 2 + y O r u 2
The penetration of the roller (r) into the U channel (u) is computed as in (23),
δ r u = d R r
where Rr = 50 mm is the outer radius of the roller. The penetration velocity of the roller into the U channel is given by (24),
δ ˙ r u = v O r u d ω r u z u Λ n u · n u
where vuOr is the velocity vector of the roller center (Or), and ωur is the angular velocity of the roller from the perspective of the reference frame attached to the U channel (F3; xu, yu, zu).
The quantity in parentheses in Equation (24) is the velocity vector of the point of contact. The contact force that the U channel exerts on the roller can be computed by substituting (23) and (24) into Equation (2), as in (25),
F u r u = σ δ r u K δ r u + γ δ ˙ r u D δ ˙ r u n u
where the superscript u indicates that the components of the vector are expressed in the reference frame (F3; xu, yu, zu), which is why (25) must be premultiplied by the transpose of the relative rotation matrix ([R]ur) before being applied to the roller. Accordingly, the U channel undergoes a force (Furu = −Fuur), which also causes a torque about point F3, given by (26).
M r u u = x E 2 u x u F r u u
2.
−11 mm ≤ xuOr ≤ 11 mm, in which case the projection of the roller center onto the line of contact falls between the edges E1 and E2, as represented in Figure 16a. The penetration of the roller into the U channel is defined in (27).
δ r u = y O r u R r
where Rr = 50 mm is the outer radius of the roller. The velocity with which the roller penetrates the U channel is given by (28),
δ ˙ r u = v O r u + y O r u ω r u x u · y u
where the quantity in parentheses is the velocity vector of the contact point Pc. The contact force that the U channel exerts on the roller can be computed by substituting (27) and (28) into Equation (2), as in (29).
F u r u = σ δ r u K δ ˙ r u + γ δ ˙ r u D δ ˙ r u y u
where the superscript u indicates that the components of the vector are expressed in the reference frame (F3; xu, yu, zu), which is why (29) must be premultiplied by the transpose of the relative rotation matrix ([R]ur) before being applied to the roller (19). Accordingly, the U channel undergoes a force (Furu = −Fuur), which also causes a torque about point F3, given by (30).
M r u u = x O r u x u F r u u
3.
xuOr < −11 mm, i.e., the projection of the roller center onto the line of contact, falls “behind” the first edge (E1). This third case is analogous to the first, the only difference being that the line of action of the force passes through point E1 rather than E2; hence, no further explanation is needed.
The interaction between every chain link and the rollers, governed by the three cases listed above, is implemented using a specially written MATLAB® function.

2.3. Ground Modeling

2.3.1. Soil Model

To estimate the sinkage and the normal pressure distribution at the vehicle–terrain interface, a model capable of describing the soil behavior when compressed in a normal direction is needed. M. G. Bekker developed the most famous and widely used model, which combines on-field measurements taken with a bevameter (a measuring tool designed by Bekker himself) with an analytical formulation to provide the following pressure–sinkage relationship (31), (also known as the Bekker equation) [22]
p = k c b + k ϕ y n
where p is the normal pressure exerted by soil, expressed in kPa; y is the vertical sinkage coordinate, expressed in m; b is the smaller dimension of the contact patch, expressed in m; and n, kc, and kϕ are the pressure–sinkage parameters.
Equation (31) is essentially an empirical equation, according to which the parameters kc and kϕ have variable dimensions, depending on the value of the exponent n (equal to kN/mn+1 and kN/mn+2, respectively). The exponent n is dimensionless.
Based on a more rigorous approach to soil mechanics and experimental evidence, Reece proposed the following pressure–sinkage relationship (32) for homogeneous soil [23],
p = c k c + γ s b k Φ ( y b ) n
where c is the cohesion of the terrain, expressed in kPa; γs is the weight density of the terrain, expressed in kg/m3; and n, kc, and kϕ are the pressure–sinkage parameters for the Reece equation. The constants c, γs, k′c, and k′ϕ in Equation (32) can be combined as k’′c =ck′c and k’′ϕ=γs k′ϕ to yield (33).
p = k c + b k Φ ( y b ) n
The parameters n, k″c, and kϕ are derived from a procedure for fitting pressure–sinkage data with the Reece Equation (33); their mean values for various mineral terrains are tabulated in [2]. Because the parameters k″c and kϕ, unlike what happens with Equation (31), have constant dimensions (kN/m2 and kN/m3, respectively). Equation (33) was chosen to characterize the pressure–sinkage relationship of the modeled soil. Furthermore, a viscous linear term was introduced in Equation (33) to account for energy dissipation, yielding (34).
p = k c + b k Φ ( y b ) n + D g y ˙
where y ˙ is the sinkage velocity (expressed in m/s), and Dg is the soil damping coefficient [24], the value of which was set to 500 kPa·s/m.

2.3.2. Contact between Track Links and Ground

The soil is represented as a brick body fixed to the world frame.
As shown in Figure 17, a coordinate system is associated with the ground body, with an origin at point Og and:
  • The xg axis lying on the upper plane of the brick solid representing the soil and oriented from left to right (fixed horizontal direction);
  • The yg axis normal to the upper plane of the brick solid representing the soil and oriented upwards (fixed vertical direction); and
  • The zg axis emerging from the plane of the diagram in Figure 17 so that the right-hand rule is satisfied.
The dimensions of the body representing the ground are 1.5 m, 0.1 m, and 0.7 m in the xg, yg, and zg directions, respectively. The vertical sinkages of the links that make up the track belt are evaluated with respect to the undeformed upper surface of the ground solid, i.e., the xgzg plane of the soil reference system. As described, each track link has a rubber lug meant to come in contact and exchange forces with the ground. To quantify the sinkage of the lower surfaces of the lugs, two points (F4 and F5) are attached to each one of these elements, as represented in the diagram in Figure 17. A reference frame is placed at the midpoint (F6) of segment F4F5, with:
  • the xl axis parallel to the contact surface (shown in blue in Figure 17) and oriented from point F4 to point F5;
  • The yl axis normal to the contact surface and oriented towards the inside of the rubber lug; and
  • The zl axis emerging from the plane of the diagram in Figure 17 so that the right-hand rule is satisfied.
The reference frame (F6; xl, yl, zl) is fixed to the rubber lug.
The mathematical model of deformable soil defined by the non-linear elastic relationship (34) was integrated into the MATLAB® Simscape™ environment to correctly simulate the contact and the behavior of the track on agricultural terrain. In particular, at each time step of the simulation of the multibody model:
1.
The vertical sinkages ygF4 and ygF5, together with the sinkage velocities y ˙ F 4 g and y ˙ F 5 g , of the points F4 and F5, respectively, are evaluated with respect to the ground reference system (Og; xg, yg, zg);
2.
The normal pressures p4 and p5, acting at points F4 and F5, respectively, are calculated according to Equation (34), yielding (35) and (36),
p 4 = σ 4 k c + b k Φ ( y F 4 g b l ) n + γ 4 D g y F 4 g ˙
p 5 = σ 5 k c + b k Φ ( y F 5 g b l ) n + γ 5 D g y F 5 g ˙
where bl is the length of the segment F4F5, which represents the smaller dimension of the contact patch, equal to 25 mm; σ4 and σ5 are the contact-enabling factors, defined as in (37); and γ4 and γ5 are the damping-enabling factors, defined as in (38).
σ 4 = 1                 i f   y F 4 g < 0   0                   i f   y F 4 g 0             σ 5 = 1                 i f   y F 5 g < 0   0                   i f   y F 5 g 0    
γ 4 = 1                 i f   y F 4 g < 0   0                   i f   y F 4 g 0         γ 5 = 1                 i f   y F 5 g < 0   0                   i f   y F 5 g 0
The factors σi, i = 4; 5, in (37) ensure that the pressure is zero when no contact occurs between the considered lug and the ground. Similarly, the factors γi, i = 4; 5, in (38) are essential to activating the damping term of Equation (34) only when the point of the rubber lug, either F4 or F5, is moving towards the inside of the ground body;
3.
The force (Fgl) and the moment (Mgl) exerted by the ground on the rubber lug are calculated assuming a linear distribution of the pressure over the contact area (a simplification that has proven reasonable), as expressed in (39).
p x l = p 5 p 4 b x l + p 4 x F 5 l p 5 x F 4 l b
where xlF4 = −bl/2= −12.5 mm and xlF5 = bl/2= 12.5 mm are the (constant) coordinates of points F4 and F5, expressed in the reference frame (F6; xl; yl; zl) attached to the lug. If the pressure (39) is supposed to be constant along the direction normal to the plane of the diagram in Figure 17, integrating over the contact area yields (40) and (41),
F g l = w l x F 4 l x F 5 l p x l d x l y l = w l b l 2 p 4 + p 5 y l
M g l = w l x F 4 l x F 5 l p x l x l d x l z l = w l b l 2 12 p 5 p 4 z l
where wl is the width of the track, equal to 160 mm, i.e., the dimension of the rubber lug along the direction of the zl axis; and yl and zl are unit vectors of the yl zl axis of the reference frame attached to the lug, respectively. The force lies along the direction normal to the contact surface, whereas the torque acts in the longitudinal plane of the track module, as represented in Figure 17. The previously calculated force (Fgl) and moment (Mgl) are applied to the considered rubber lug. The steps listed above are repeated for every link composing the track belt in each simulation time step. The ground experiences the force (Flg = −Fgl) and the moment (Mlg = −Mgl) as a consequence of the contact with the single lug; however, these produce no effect, as the soil is rigidly connected to the fixed world frame.

2.3.3. Contact between Rear Wheels and Ground

The purpose of modeling half the robot necessarily requires considering the wheeled rear axle’s presence and the reconfigurable module on the front. Because the vehicle does not move during the reconfiguration phases, the wheel is rigidly attached to the frame. The method adopted for the implementation of the contact between the wheel and the ground is that proposed by Bekker [25], which is part of the mechanics of a rigid wheel over unprepared terrain and relies on three major assumptions:
  • The terrain reaction at all points on the contact patch is purely radial;
  • The terrain reaction is equal to the normal pressure beneath a horizontal sinkage plate at the same depth in a pressure–sinkage test, as expressed by Equation (35); and
  • No shear stress acts on the wheel–terrain interface.
Figure 18 represents a diagram of the simplified wheel–soil interaction model described here. The center point of the wheel body is denoted by the symbol Ow; it is the origin of the right-hand reference frame integral with the wheel, with xw, yw, and zw axes. The right-hand coordinate system (Og; xg, yg, zg) is associated with the ground body. To implement the force exchanged between the wheel and, consequently, the vehicle frame and the body representing the ground, the following steps are taken:
1.
The maximum sinkage (y0) is calculated as in (42):
y 0 = σ w R w y O w g
where ygOw is the y coordinate of the wheel center (Ow) in the ground reference (Og; xg, yg, zg), Rw is the outer wheel radius (equal to 0.17 m), and σw is the contact-enabling factor, defined as in (43):
σ w = 1                 i f   y O w g R w   0                   i f   y O w g > R w  
This ensures that the maximum sinkage is zero when there is no contact between the wheel and soil. Furthermore, the y component of the wheel center velocity expressed in the ground reference frame (Og; xg, yg, zg), y ˙ gOw, is used to calculate the sinkage velocity ( y ˙ 0) as in (44),
y 0 ˙ = σ w γ w y ˙ O w g )
where γw, similarly to (38), is a factor needed to activate damping only when the wheel moves towards the inside of the ground body.
2.
The normal pressure distribution over the contact area is computed according to Equation (34) by (45),
p θ w = k c + b w k Φ ( y θ w b w ) n + D g y ˙ θ w
where θw is the contact angle of the wheel, ranging between zero and the maximum value (θ0) given by (46); bw is the width of the wheel measured in the direction normal to the plane of the diagram in Figure 18, equal to 0.16 m; and y θ w , y ˙ θ w are the sinkage and sinkage velocity functions of the contact angle as expressed by (47).
θ o = arccos ( R w y 0 R w )
y ( θ w ) = R w c o s θ w 1 + y 0 y ˙ ( θ w ) = c o s θ w y 0 ˙
3.
The force exerted by the ground on the wheel (Fgw) is determined by taking the integral of the pressure distribution (45) over the contact area, as in (48),
F g w = 2 b w R w 0 θ 0 p θ w cos θ w d θ w y w
where yw is the unit vector of the yw axis of the reference frame attached to the wheel. In deriving the force (48), the pressure (p) is assumed to be uniform along the direction of the zw axis, i.e., normal to the plane of the diagram in Figure 18. The integral in (48) is solved numerically at each simulation time step. The force (Fgw) is then applied to the wheel. Because the shear stress at the wheel–terrain interface is supposed to be null, and the pressure (45) is symmetrical with respect to the wheel centerline, no torque contribution is included in this case.
Once again, the reaction force (Fwg = −Fgw) of the wheel on the ground has no effect.

3. Results and Discussion

The normal pressure distribution over the contact area between the track module and the ground is one of the most relevant, if not the principal, outcomes of the analysis and simulation of the implemented model (see Figure 19 and Figure 20).
The reasons why it is important to predict how pressure distributes at the track–ground interface are diverse. First, the resistance opposed by the ground to the motion of the vehicle due to terrain compaction (Rc) is proportional to the work done in compacting the terrain [25], as shown in (49),
W o r k = b 0 l 0 y x p y x d y d x
where b is the width of the track contact area (equal to the width of the rut made by the track), l is the length of the track contact area (equal to the length of the rut made by the track), x is the distance from the front of the contact area, and y(x) is the sinkage of the track in the point of coordinate x along the contact area. Because Rcl = Work, Equation (49) yields (50),
R c = b l 0 l 0 y x p y x d y d x
From which it is evident that the motion resistance depends on the pressure distribution between the track and the ground. Usually, the uniform contact pressure is assumed in the first place so that the integral in (50) can be calculated and the resistance can be assessed. Other pressure distributions, such as trapezoidal or sinusoidal distributions, can be used to obtain Rc from (50); however, all of these lead to first approximation estimates. As stated in [14], assuming a more uniform pressure distribution than the real distribution leads to underestimated results, as in this case, as numerically verified by applying equation (50) to calculate the rolling resistance both with the pressure distribution obtained from the simulation and a uniform distribution deduced by dividing a quarter of the vehicle weight by the track contact area (see Table 2).
The shearing of the terrain produces the tractive effort developed by the track or wheel. The torque applied to the tire or the track sprocket initiates a shearing action over the vehicle running gear–soil interface.
To predict vehicle thrust and associated slip, the relation between shear stress and the shear displacement for the considered terrain is needed, as highlighted by (51)
F = b 0 l τ j d x
where F is the tractive effort, and τ and j are the shear stress and displacement, respectively. With the sole objective of qualitatively demonstrating the influence of normal pressure distribution on the generation of tractive effort, among the three most observed types of shear stress–shear displacement semiempirical relationship, we select the exponential function of the form proposed by Janosi and Hanamoto [26] so that Equation (51) becomes (52),
F = b 0 l c + p t a n ϕ 1 e j k d x
where c and ϕ are the cohesion and angle of internal shearing resistance of the terrain, respectively, and K is the shear deformation modulus, representing a measure of the magnitude of the required shear displacement.
Equation (52) demonstrates how the tractive effort is affected, among other factors, by the normal pressure distribution (p) over the contact area. Although traction is not subject to the present work, it is notably influenced by the pressure distribution, the shape of which is therefore a matter of absolute interest.
Another reason for being interested in assessing the pressure distribution over the contact area is dictated by the peculiarity of the proposed system itself. Because wheel/track reconfiguration plays a leading role, it is required to monitor the development of contact pressure in response to what the actuator is controlled to do. In other words, it is important to predict the shape of the pressure curve, particularly its peak value, as the target imposed by the control system varies. This aspect is closely associated with the ultimate purpose for which the system is controlled, as the goal is to achieve a track configuration in which the peak pressure and, consequently, soil compaction, are minimized. In Figure 21, a side view of the front module captured during the track mode interval of the simulation cycle (from t2 = 20 s up to t3 = 25 s) is shown. As evidenced by the deflection of the chain due to the applied load and consequent deformation of the ground, the accurate modeling detailed in Section 2 has made it possible to achieve a plausible sinkage profile of the track. The rigid links making up the chain that constitutes the track belt can be divided into two groups: those that are in contact both with the terrain and the roadwheels (with this name indifferently denoting the sprocket, the idler, or the rollers) and the others in contact with the terrain only. The shape of the track segments in contact with the roadwheels stems from the shape of the roadwheels themselves.
In contrast, the sinkage of the track in contact with the terrain only is determined by the roadwheel spacing and the chain tension, as well as the assumed pressure–sinkage relationship. The observed sinkage profile resembles the normal pressure graphed in Figure 21; as expected, the distribution is far from uniform, as it presents with peaks under the sprocket, the two rows of rollers, and the idler. The lowest pressure value occurs where the spacing is greatest, i.e., in the section between the sprocket and the first row of rollers.
The system must be capable of moving from an initial posture to a final assigned posture, so the extremal points of the motion path must be defined. The wheel mode is achieved when the idler wheel is sufficiently distant from the ground, i.e., when the inclination angle of the track frame reaches a value of ϑT,wheel = 10°. Regarding the track mode, it is difficult, if not impossible, to know a priori what position the mechanism has to assume to produce an optimal pressure distribution over the contact area. Therefore, the definition of the inclination angle (ϑT,track) corresponding to the track configuration is more formal than substantial, as there is no way to select a single, specific value; instead a range in the proximity of a “geometric” reference position must be examined (53):
θ T , t r a c k = sin 1 R i R s l 1 T δ θ T
where l1T is the length of the first arm of the track frame (segment T1T2 in Figure 2), and Ri and Rs are the outer radii of the idler and the sprocket, respectively. The first term in Equation (53) corresponds to the posture in which, supposing an ideal, perfectly flat, rigid ground, the lower branch of the track belt would be horizontal, and the entire track area would contact the soil. The second term in Equation (53) (δϑT), later referred to as control target deviation, can take either positive or negative values close to zero; its effect is demonstrated by the pressure curves in Figure 21:
  • The blue curve corresponds to a control target deviation of δϑT = −0.5°, which yields ϑT,track = −11.61°. In this case, the actuator is controlled to reach an overall length of 338.8 mm (the plunger is driven to travel 82.8 mm along the available stroke).
  • The orange curve corresponds to a control target deviation of δϑT = 0°, yielding ϑT,track = −12.11°. In this case, the actuator is controlled to reach an overall length of 341 mm (the plunger is driven to travel 85.0 mm along the available stroke).
  • The yellow curve corresponds to a control target deviation of δϑT = 0.5°, yielding ϑT,track = −12.61°. In this case, the actuator is controlled to reach an overall length of 343.3 mm (the plunger is driven to travel 87.3 mm along the available stroke).
Observation of the three curves mentioned above indicates the significance of the effect of the control deviation (δϑT) on the pressure distribution. According to Equation (53), the higher δϑT is, the harder the control system forces the actuator to push the track frame against the ground; as a consequence, the peak of pressure beneath the idler (and partly that beneath the second row of rollers) increases at the expense of that under the sprocket. This is easily understood by considering that as δϑT increases from −0.5° to 0.5°, as shown in Figure 21 (from the blue curve to the yellow curve), the actuator drives the rocker slightly further and causes the shock absorber to exert an increased compressing action on the track frame.
As an extreme consequence of this reasoning, if δϑT exceeds a threshold limit, an undesired situation may occur in which the peak pressure beneath the sprocket vanishes; this means that the sprocket is detached from the ground and the vehicle front axle rests on the idler only. If δϑT were too low instead, the part of the track belt beneath the rollers and the idler would barely touch the ground, providing insufficient contact pressure. Accordingly a value or a narrow range of values can be assigned to δϑT to generate a distribution characterized by minimum peak pressure. The highlighted relationship between δϑT and the peak of pressure beneath the idler is difficult to assess, as it vitally depends on knowledge of the exact geometry of the parts involved, as well as the exact value of all the elastic parameters. Therefore, a small mistake in any procedural step could cause the reconfiguration to fail, leading to a lack of contact or the generation of excessive force. The switching system, especially during wheel-to-track reconfiguration must exert a definite force on the external environment rather than reaching a certain position. Therefore, the task for which the system is intended should be approached with a force–control or impedance–control strategy [27,28]. Nonetheless, the position control scheme adopted here is a possible way to approach the question, even if only for the practicality of implementing position feedback in a first embodiment of the system, as opposed to the necessity of obtaining a measured and controlled force. With the implemented model, it is possible to evaluate the effect of variation of the control target, together with other important parameters of the system, on the contact pressure peak between the track and the ground.

3.1. Effect of Chain Tension

In Figure 22, the effect of the track belt tension on the contact peak pressure is examined for a control target deviation (δϑT) in the range of −1° to 1°. Three curves are shown:
  • The blue curve corresponds to null tension, i.e., the chain is loose;
  • The orange curve corresponds to a tension equal to 25% of the vehicle weight, i.e., 0.981 kN; and
  • The yellow curve corresponds to a tension equal to 50% of the vehicle weight, i.e., 1.962 kN.
As the track belt tension increases, the contact peak pressure beneath the track decreases. This is beneficial effect, as the tension stiffens the system and therefore helps to redistribute the pressure more evenly. From the point of view of the control target deviation (δϑT), the three curves show no substantial difference.

3.2. Effect of Shock Absorber Spring Rate and Preload

Figure 23 shows the effect of the stiffness of the shock absorber spring on the contact peak pressure for a control target deviation (δϑT) in the range of −1° to 1°. Three curves are shown:
  • The blue curve corresponds to a shock absorber with a spring rate equal to 40 N/mm;
  • The orange curve corresponds to a shock absorber with a spring rate equal to 100 N/mm; and
  • The yellow curve corresponds to a shock absorber with a spring rate of 160 N/mm.
The above-listed alternatives cover a wide range of elasticity, from a soft spring to the hardest one. All three curves exhibit a minimum point corresponding to different δϑT values. The blue curve reaches the minimum at δϑT = 0.875°, whereas for the orange and the yellow curves, the minimum peak pressure is obtained at δϑT = −0.375° and δϑT = −0.5°, respectively. Such differentiations are easily understood by considering that the shock absorber and the ground can be thought of as two springs in series, with stiffnesses of ks and kg, respectively, so the overall stiffness is defined as in (54)
k = k s k g k s + k g = k g 1 + k g k s
Equation (52) shows that if a certain force is to be developed on a given terrain, the softer the spring (the lower ks), the further the control system has to push the idler to achieve a minimum peak pressure (the higher the δϑT). An comparable result is obtained by varying the preload of the spring, as shown in Figure 24.
Figure 25 shows the effect of the position of the center of gravity of the vehicle frame on the contact peak pressure.
Such effect is important because it is reasonable to expect that the distribution of mass on the vehicle (specifically the payload) will likely change during normal operation. The control target deviation (δϑT) is given for the x axis and varies from −1° to 1°, whereas for the y axis, it is the peak of the pressure distribution between track and ground. Three curves are shown:
  • The blue curve corresponds to a frame with the center of gravity displaced by 115 mm towards the front axle with respect to the middle point of the wheelbase;
  • The orange curve corresponds to a frame with the center of gravity in the middle of the vehicle wheelbase; and
  • The yellow curve corresponds to a frame with a center of gravity 115 mm towards the rear axle with respect to the middle point of the wheelbase.
The vehicle wheelbase is 768 mm. The distance separating the frame center of mass from the middle of the wheelbase for the blue and yellow curves is set as equal to 15% of the wheelbase. First, a substantial increase in the minimum peak pressure is observed when passing from the yellow to the blue curve (from about 27 kPa to 43 kPa); this is easily explained by the increased fraction of the vehicle weight supported by the front wheel/track reconfigurable axle. For the same reason, as the center moves towards the front axle, the control target deviation necessary to achieve the minimum peak pressure increases: equal to −0.75° for the yellow curve, −0.375° for the orange curve, and 0.125° for the blue curve.

4. Conclusions

The development of a multibody model to study the pressure distribution of a reconfigurable wheeled/tracked vehicle on soft terrain is reported herein. The main aim of this research was the construction of a reconfigurable wheel/track locomotion system to be installed on board a grape-transporting robot to switch between a half-tracked configuration and a wheeled configuration depending on the soil type and conditions. The objective was to realize a running gear capable of actively adapting to the environment to attain improved performance while preserving the porous structure of the soil, which is often damaged by the intensive use of heavy machinery, to the greatest extent possible.
The scope of the presented multibody model is to simulate the wheel-to-track and track-to-wheel switching phases. As detailed in Section 2, the goal in question was achieved through the accurate modeling of the multiple interactions between each of the links of the chain that constitute the track belt and the various elements that are part of the system, i.e., the sprocket wheel, the idler wheel, the undercarriage rollers, the soil, etc.. The chain that constitutes the track is modeled as groups of rigid links constrained by revolute joints to account for the flexibility of the track itself. A specifically designed algorithm was implemented through MATLAB® functions to compute the contact forces applied to the track elements. The model used to characterize the soil behavior when compressed in a normal direction is based on Reece’s non-linear elastic pressure–sinkage relationship. Furthermore, the contact between the rear wheel and the ground was accounted for to ensure the possibility of assessing the effect of the position of the center of mass of the vehicle frame on the pressure distribution over the ground contact patch.
As mentioned in Section 3, the proposed multibody model allows for realistic simulation of the vehicle reconfiguration from half-tracked mode to wheeled mode and vice versa and to obtain useful information regarding both the controlled actuation unit and the contact of the track module with the ground. The main outcome is the normal pressure distribution between the track module and the ground; when in tracked configuration, it presents a non-uniform shape with peaks directly beneath the sprocket wheel, the idler wheel, and the rollers. In particular, the pressure distribution in track mode is strongly influenced by the action exerted on the system by the controlled actuator; the further the plunger is pushed forward, the more the pressure peak under the idler increases. Numerous simulations of the system prove that there exists a value, or at least a range of values, of the control target that yields the optimal ground pressure distribution characterized by the minimum peak. Herein, we demonstrated that such an optimum range is highly affected by the spring rate and preload of the shock absorber, as well as the position of the vehicle center of mass, with a much more minimal influence of the chain tension.

Author Contributions

Conceptualization A.G. and E.d.M.; methodology, A.G.; software, A.G. validation, A.G., E.d.M., N.I.G. and K.I.; formal analysis, A.G. and E.d.M.; investigation, A.G. and E.d.M.; resources, N.I.G. and K.I.; data curation, A.G.; writing—original draft preparation, A.G. and N.I.G.; writing—review and editing, A.G., E.d.M. and N.I.G.; visualization, A.G. and E.d.M.; supervision, N.I.G. and K.I.; project administration, E.d.M.; funding acquisition, N.I.G. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Muro, T.; O’Brien, J. Terramechanics: Land Locomotion Mechanics; A. A. Balkema Publishers: London, U.K., 2005. [Google Scholar]
  2. Wong, J.Y. Terramechanics and Off-Road Vehicle Engineering: Terrain Behaviour, Off-Road Vehicle Performance and Design; Butterworth-Heinemann: Oxford, UK, 2010. [Google Scholar]
  3. Nicolini, A.; Mocera, F.; Somà, A. Multibody simulation of a tracked vehicle with deformable ground contact model. Proc. Inst. Mech. Eng. Part K J. Multi-Body Dyn. 2019, 233, 152–162. [Google Scholar] [CrossRef]
  4. Slättengren, J. Utilization of ADAMS to Predict Tracked Vehicle Performance. SAE Transactions. 2000; Volume 109, pp. 1–7. JSTOR. Available online: http://www.jstor.org/stable/44650730 (accessed on 1 January 2022).
  5. Solis, J.M.; Longoria, R.G. Modeling track–terrain interaction for transient robotic vehicle maneuvers. J. Terramech. 2008, 45, 65–78. [Google Scholar] [CrossRef]
  6. Mocera, F.; Somà, A.; Nicolini, A. Grousers Effect in Tracked Vehicle Multibody Dynamics with Deformable Terrain Contact Model. Appl. Sci. 2020, 10, 6581. [Google Scholar] [CrossRef]
  7. Mahalingam, I.; Padmanabhan, C. A novel alternate multi-body model for the longitudinal and ride dynamics of a tracked vehicle. Veh. Syst. Dyn. 2021, 59, 433–457. [Google Scholar] [CrossRef]
  8. Ma, Z.D.; Perkins, N.C. A Super-Element of Track-Wheel-Terrain Interaction for Dynamic Simulation of Tracked Vehicles. Multibody Syst Dyn. 2006, 15, 347–368. [Google Scholar] [CrossRef]
  9. Lamandé, M.; Greve, M.H.; Schjønning, P. Risk assessment of soil compaction in Europe—Rubber tracks or wheels on machinery. Catena 2018, 167, 353–362. [Google Scholar] [CrossRef]
  10. Grazioso, A.; Galati, R.; Mantriota, G.; Reina, G. Multibody Simulation of a Novel Tracked Robot with Innovative Passive Suspension. In Mechanisms and Machine Science; Niola, V., Gasparetto, A., Quaglia, G., Carbone, G., Eds.; Advances in Italian Mechanism Science. IFToMM Italy 2022; Springer: Cham, Switzerland, 2022; Volume 122. [Google Scholar] [CrossRef]
  11. Grazioso, A.; Di Maria, E.; Margheriti, S.; Gallone, R.; Prete, G. Modellazione della Sospensione Formula SAE SRT15 in MSC Adams/Car. November 2014. Available online: https://www.mscsoftware.com/it/page/modellazione-della-sospensioneformula-sae-srt15-msc-adamscar (accessed on 1 January 2022).
  12. Visconte, C.; Cavallone, P.; Carbonari, L.; Botta, A.; Quaglia, G. Design of a Mechanism with Embedded Suspension to Reconfigure the Agri_q Locomotion Layout. Robotics 2021, 10, 15. [Google Scholar] [CrossRef]
  13. Jiang, H.; Xu, G.; Zeng, W.; Gao, F.; Chong, K. Lateral Stability of a Mobile Robot Utilizing an Active Adjustable Suspension. Appl. Sci. 2019, 9, 4410. [Google Scholar] [CrossRef] [Green Version]
  14. di Maria, E. Running Gear Size Selection of a Wheel/Track Reconfigurable Grape Transporting Vehicle by FEM Analysis. Ph.D. Thesis, Graduate School of Life Science and Systems Engineering, Kyushu Institute of Technology, Fukuoka, Japan, 2021. [Google Scholar]
  15. Cheli, F.; Pennestrì, E. Cinematica e Dinamica dei Sistemi Multibody. Casa Editrice Ambrosiana: Milano, Italy, 2009; Volume 1. [Google Scholar]
  16. Flores, P.; Lankarani, H.M. Contact Force Models for Multibody Dynamics; Springer: Berlin/Heidelberg, Germany, 2016. [Google Scholar] [CrossRef]
  17. Lankarani, H.M.; Nikravesh, P.E. Continuous contact force models for impact analysis in multi-body systems. Nonlinear Dyn. 1994, 5, 193–207. [Google Scholar] [CrossRef]
  18. Nikravesh, P.E. Planar Multibody Dynamics: Formulation, Programming with MATLAB®, and Applications; CRC Press: Boca Raton, FL, USA, 2018. [Google Scholar]
  19. Omar, M.A. Chain Drive Simulation Using Spatial Multibody Dynamics. Adv. Mech. Eng. 2014, 6, 378030. [Google Scholar] [CrossRef] [Green Version]
  20. Pedersen, S. Simulation and Analysis of Roller Chain Drive Systems. Ph.D. Thesis, Department of Mechanical Engineering, Technical University of Denmark, Lyngby, Denmark, 2004. [Google Scholar]
  21. Pedersen, S. Model of contact between rollers and sprockets in chain-drive systems. Arch. Appl. Mech. 2005, 74, 489–508. [Google Scholar] [CrossRef]
  22. Bekker, M.G. Introduction to Terrain-Vehicle Systems; University of Michigan Press: Ann Arbor, MI, USA, 1969. [Google Scholar]
  23. Reece, A.R. Principles of Soil-Vehicle Mechanics. Proc. Inst. Mech. Eng. Automob. Div. 1965, 180, 45–66. [Google Scholar] [CrossRef]
  24. Rubinstein, D.; Coppock, J.L. A detailed single-link track model for multi-body dynamic simulation of crawlers. J. Terramechanics 2007, 44, 355–364. [Google Scholar] [CrossRef]
  25. Wong, J.Y. Theory of Ground Vehicles; John Wiley & Sons: Hoboken, NJ, USA, 2008. [Google Scholar]
  26. Janosi, Z.; Hanamoto, B. The analytical determination of drawbar-pull as a function of slip for tracked vehicles in deformable soils. In Proceedings of the 1st International Conference on the Mechanics of Soil-Vehicle Systems, Torino, Italy, 12–16 June 1961; pp. 707–736. [Google Scholar]
  27. Villani, L.; De Schutter, J. Force Control. In Springer Handbook of Robotics; Siciliano, B., Khatib, O., Eds.; Springer Handbooks; Springer: Cham, Switzerland, 2016. [Google Scholar] [CrossRef]
  28. De Schutter, J.; Van Brussel, H. Compliant Robot Motion II. A Control Approach Based on External Control Loops. Int. J. Robot. Res. 1988, 7, 18–33. [Google Scholar] [CrossRef]
Figure 1. 3D view of the multibody model of the farming robot.
Figure 1. 3D view of the multibody model of the farming robot.
Machines 10 01117 g001
Figure 2. Scheme of the main components of the model.
Figure 2. Scheme of the main components of the model.
Machines 10 01117 g002
Figure 3. Track module.
Figure 3. Track module.
Machines 10 01117 g003
Figure 4. Detail the track belt (left) and the track link (right) components.
Figure 4. Detail the track belt (left) and the track link (right) components.
Machines 10 01117 g004
Figure 5. Components of the track module and final assembly.
Figure 5. Components of the track module and final assembly.
Machines 10 01117 g005
Figure 6. The undercarriage of the track module: rollers and L–shaped supports.
Figure 6. The undercarriage of the track module: rollers and L–shaped supports.
Machines 10 01117 g006
Figure 7. Track frame connecting the sprocket shaft and the idler shaft.
Figure 7. Track frame connecting the sprocket shaft and the idler shaft.
Machines 10 01117 g007
Figure 8. Track frame connecting the sprocket shaft (in blue) and the idler shaft.
Figure 8. Track frame connecting the sprocket shaft (in blue) and the idler shaft.
Machines 10 01117 g008
Figure 9. Scheme of the track belt tensioner.
Figure 9. Scheme of the track belt tensioner.
Machines 10 01117 g009
Figure 10. Detailed photo of the track belt.
Figure 10. Detailed photo of the track belt.
Machines 10 01117 g010
Figure 11. Links of the chain that constitutes the track belt: (a) outer and inner links, (b) contact surfaces.
Figure 11. Links of the chain that constitutes the track belt: (a) outer and inner links, (b) contact surfaces.
Machines 10 01117 g011
Figure 12. Linear spring–damper for frictionless contact modeling.
Figure 12. Linear spring–damper for frictionless contact modeling.
Machines 10 01117 g012
Figure 13. Sprocket gear.
Figure 13. Sprocket gear.
Machines 10 01117 g013
Figure 14. Approximation of the sprocket tooth profile using five straight-line segments.
Figure 14. Approximation of the sprocket tooth profile using five straight-line segments.
Machines 10 01117 g014
Figure 15. Illustration of the contact between a chain bushing and a sprocket tooth.
Figure 15. Illustration of the contact between a chain bushing and a sprocket tooth.
Machines 10 01117 g015
Figure 16. Illustration of the contact between a sprocket tooth and a roller pin. (a) the projection of the roller center onto the line of contact falls between the edges E1 and E2. (b) the projection of the roller center onto the line of contact falls “beyond” the second edge E2.
Figure 16. Illustration of the contact between a sprocket tooth and a roller pin. (a) the projection of the roller center onto the line of contact falls between the edges E1 and E2. (b) the projection of the roller center onto the line of contact falls “beyond” the second edge E2.
Machines 10 01117 g016
Figure 17. Illustration of the contact between rubber lugs and the ground.
Figure 17. Illustration of the contact between rubber lugs and the ground.
Machines 10 01117 g017
Figure 18. Illustration of the contact between the rear wheel and the ground.
Figure 18. Illustration of the contact between the rear wheel and the ground.
Machines 10 01117 g018
Figure 19. Normal pressure distribution (left y axis) and sinkage (right y axis) of the reconfigurable module in track mode.
Figure 19. Normal pressure distribution (left y axis) and sinkage (right y axis) of the reconfigurable module in track mode.
Machines 10 01117 g019
Figure 20. Normal pressure distribution (left y axis) and sinkage (right y axis) of the reconfigurable module in wheel mode.
Figure 20. Normal pressure distribution (left y axis) and sinkage (right y axis) of the reconfigurable module in wheel mode.
Machines 10 01117 g020
Figure 21. Normal pressure distribution under the reconfigurable module in track mode: effect of control target deviation (δϑT).
Figure 21. Normal pressure distribution under the reconfigurable module in track mode: effect of control target deviation (δϑT).
Machines 10 01117 g021
Figure 22. Peak pressure beneath the track module as a function of δϑT: chain tension effect.
Figure 22. Peak pressure beneath the track module as a function of δϑT: chain tension effect.
Machines 10 01117 g022
Figure 23. Peak pressure beneath the track module as a function of δϑT: spring rate effect.
Figure 23. Peak pressure beneath the track module as a function of δϑT: spring rate effect.
Machines 10 01117 g023
Figure 24. Peak pressure beneath the track module as a function of δϑT: spring preload effect.
Figure 24. Peak pressure beneath the track module as a function of δϑT: spring preload effect.
Machines 10 01117 g024
Figure 25. Peak pressure beneath the track module as a function of δϑT: center of gravity position effect.
Figure 25. Peak pressure beneath the track module as a function of δϑT: center of gravity position effect.
Machines 10 01117 g025
Table 1. Coordinates of the profile points of the first sprocket tooth.
Table 1. Coordinates of the profile points of the first sprocket tooth.
Pointx (mm)y (mm)
P1,113.4140.0
P2,112.0140.1
P3,15.2125.6
P4,1−5.2125.6
P5,1−12.0140.1
P6,1−13.4140.0
Table 2. Motion resistance for three soil types: values calculated with a uniform pressure and with simulated pressure.
Table 2. Motion resistance for three soil types: values calculated with a uniform pressure and with simulated pressure.
Soil Type (Moisture Content %)Rolling Resistance (N)
Uniform PressureSimulated Pressure
Grenville loam (18.2%) 5.869.28
Upland sandy loam (44.3%)3.585.10
North Gower clayey loam (52.0%)11.4216.83
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Grazioso, A.; di Maria, E.; Giannoccaro, N.I.; Ishii, K. Multibody Modeling of a New Wheel/Track Reconfigurable Locomotion System for a Small Farming Vehicle. Machines 2022, 10, 1117. https://doi.org/10.3390/machines10121117

AMA Style

Grazioso A, di Maria E, Giannoccaro NI, Ishii K. Multibody Modeling of a New Wheel/Track Reconfigurable Locomotion System for a Small Farming Vehicle. Machines. 2022; 10(12):1117. https://doi.org/10.3390/machines10121117

Chicago/Turabian Style

Grazioso, Andrea, Enrico di Maria, Nicola Ivan Giannoccaro, and Kazuo Ishii. 2022. "Multibody Modeling of a New Wheel/Track Reconfigurable Locomotion System for a Small Farming Vehicle" Machines 10, no. 12: 1117. https://doi.org/10.3390/machines10121117

APA Style

Grazioso, A., di Maria, E., Giannoccaro, N. I., & Ishii, K. (2022). Multibody Modeling of a New Wheel/Track Reconfigurable Locomotion System for a Small Farming Vehicle. Machines, 10(12), 1117. https://doi.org/10.3390/machines10121117

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop