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),
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]
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):
where
and
are the position vectors of point
A and point
B, respectively;
is the position of
B relative to
A; and
is the penetration velocity of
B into
A, as expressed by (4):
where
and
are the velocity vectors of points
A and
B, respectively;
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),
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.
- -
is the damping-enabling factor, defined as in (6)
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 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).
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).
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),
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):
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)
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 (tk∗n) and the previously determined normal vectors (nk∗n, 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
k∗th tooth with respect to the bushing center (
F1) are computed as (12)
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),
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):
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)
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):
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).
- 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
from (16) into Equation (1)
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).
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).
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),
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).
The penetration of the roller (
r) into the U channel (
u) is computed as in (23),
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),
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),
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).
- 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).
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),
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).
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).
- 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]
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],
where
c is the cohesion of the terrain, expressed in kPa;
γs is the weight density of the terrain, expressed in kg/m
3; and
n,
k′
c, 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).
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/m
2 and kN/m
3, 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).
where
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
xg –
zg 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 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 and , 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),
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).
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).
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),
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):
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):
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),
gOw, is used to calculate the sinkage velocity (
0) as in (44),
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),
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
,
are the sinkage and sinkage velocity functions of the contact angle as expressed by (47).
- 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),
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),
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),
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)
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),
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):
where
l1T is the length of the first arm of the track frame (segment
T1T2 in
Figure 2), and R
i 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)
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.