Transient load simulation of forwarder rear frame

One difficulty in the design of the load bearing components of mobile machines is the transient and non-linear nature of the loads acting on them. A common method for tracking these loads is to use strain gauges and force transducers on a physical test machine. An alternative method for determining the transient loads by means of a mathematical model that intends to describe the response of a John Deere 1010E forwarder as it crosses a test track is utilized in this study. The model is based on finite element method and it is solved using explicit time integration and LS-DYNA® software. As a result of this study a model capable of replicating the real world with a reasonable accuracy was obtained. The forces acting on tires, which can be considered the most important results of this work, can be used as boundary conditions in consequent analyses, such as fatigue simulation.


Introduction
The purpose of this work is to build a simulation model that is capable of replicating the reaction forces measured from physical field test of a John Deere 1010E forwarder during June 2015.Simulation will reduce the need for physical testing which in turn will speed up the design process and possibly reduce the costs of machine development and the time spent on the design phase.
As forestry machines are driven in rough, constantly altering forest terrains, the forces acting on them are difficult to predict without physical testing.These forces depend on the stiffness properties of both the machine and the terrain, velocity of the machine, total mass of the machine and the geometry of the terrain.Comprehensive physical tests are expensive and they are unsuitable for testing new models early in the design phase.By using a simulation model that is capable of determining the forces exerted from the terrain to the machine and computing the response to these forces in desired operation conditions, new designs and effects of geometry modifications in them can be easily studied.
The frames of the forestry machines consist of welded sheet metal parts.Weld seams exposed to high stresses are prone to fatigue damage.By using the finite element method and the right loading conditions, the optimum design can be found.

Computational method
The Finite Element Method (FEM), in general, is a numerical method used for solving partial differential equations by discretizing the studied geometry using a finite amount of different types of elements [26 p. 2].FEM is used for complex geometries in situations where analytical solutions cannot be found.The method used in this simulation uses explicit time integration.The equation of motion to be solved is where the load vector is denoted with (), nodal acceleration vector (second time derivative of nodal displacement ) with ̈ and mass matrix with [].The load vector includes internal, external and damping forces, thus  = (, ̇, , ).Equation ( 1) represents the equilibrium equation for d'Alembert's principle, which states that external forces () acting on the domain must be equal to the sum of inertial ([]̈), damping ([]), and internal forces due to stiffness of the structure.Hourglass forces included in the load vector are unphysical and therefore the accumulated energy generated by these forces should be as low as possible as the solution proceeds.Due to the displacement dependency of internal force vector, Equation ( 1) is a non-linear ordinary differential equation, for which analytical solutions cannot be found and numerical methods must be used instead [13 pp. 611-612].Initial conditions ̈0, ̇0 and  0 are assumed to be known.Accelerations can be solved from Equation (1) as Now that accelerations at time  are known, the velocities at time  + ∆ 2 ⁄ can be calculated as And furthermore the nodal displacements  at time  + ∆ can be calculated using the central difference method employed in explicit time integration as This process is repeated until the desired time  is reached.Equations ( 2)-( 4) can be solved directly without iteration, hence the term explicit.
If a lumped mass matrix is used the equation of motion is uncoupled as opposed to the equations of implicit dynamics, and therefore the solution time for a single time step is smaller in explicit dynamics than it is in implicit.The inversion of a lumped (diagonal) matrix reduces to a simple division by a scalar in Equation (2), which is cheap to compute compared to the factorization of a consistent matrix [16 p. 24; 26 p. 648].However, due to the conditional stability of the explicit method (considered in the next section), the time step to be used is usually much smaller compared to the time step used in the implicit method.[2,8,9,13,16] As it was stated, the equations of the explicit time integration method are uncoupled, meaning that the results for each element are calculated separately.Therefore during one time step, the information from one side of the element can only travel to the other side of the element to keep the solution stable.The physical interpretation for this is the wave propagation in the structure (or highest natural frequency of the structure).The maximum time step ∆ that can be used in the analysis to ensure a stable solution is limited by Courant-Friedrichs-Lewy (CFL) condition as where  is the speed of sound in the material,  is the time step safety factor used in the analysis to increase stability (defaults to 0.9) and ℎ is the characteristic element dimension.The CFL condition restricts the wave from traveling more than the dimension ℎ in one time step.[2,5,8] The element that minimizes condition (5) dictates the time step to be used.Speed of sound also depends on stress state and damping properties of the material, but crude evaluation of the time step can be made even by neglecting their effect.Hexahedral/quad mesh is preferred over tetrahedral/triangular mesh to prevent the time step from getting too small while retaining solution accuracy.Smallest geometry details are de-featured and simplified to achieve reasonable computation times.
If the mesh contains only few small elements, automatic mass scaling can be used to prevent the time step from becoming too small.The speed of sound (elastic speed wave) in material is calculated as (for beam element) [13 p. 600] where  is the elastic modulus and  the density of the material,  is the volume and  is the mass of the element.By substituting (6) to (5), ∆ results as As can be seen, the minimum time step is directly proportional to the square root of the mass of the element.By scaling the mass of a single element the error in the total mass of the structure is small but the benefits from a greater time step might be significant.Mass scaling is used in this analysis.

Field test
The test track is located at the John Deere Forestry Tampere facility.The simulated part of the test track consists of semicircular bumps that are arranged as shown in Figure 1.
Bumps are constructed of steel plates with a thickness of 12 mm.The machine is fitted with customized wheel hubs that include force transducers.The radial force reactions measured in each rear frame tire during the drive around the test track with a 12 ton load are shown in Figure 2. The machine crosses the studied portion of the test during the first 35 seconds (Figure 3 shows the same results as Figure 2, but with horizontal axis limited).Another bump track with a different bump arrangement is crossed at a time period ranging from approximately 45 to 75 seconds.The right front (RF, see Figure 4) tire of the rear frame is the first one to hit the bump (no. 1 in Figure 1).This is shown as a peak in the reaction force (approximately 83 kN).The second high peak is the result of left front tire (LF) hitting bump no. 2 (approximately 89 kN).It can be seen that the forces acting on the rear tires (RR and LR) are lower than those acting on the front tires.

Model
The details of modeling for all components in the calculation model are presented in this chapter.The justifications for the chosen element types, simplifications compared to the physical real life machine and for the material models are also given.Since this is a dynamic analysis, it is important that the mass in each component is close to the mass of the real life component it represents in order to achieve correct inertial behavior.The densities of the simulation model components are modified to match weights to measured values shown in Table 1.Mesh and component positions at initial configuration are shown in Figure 5.The time step is set to 2 microseconds.The minimum allowed element dimension for each element type using this time step size is calculated based on the material properties using Equation (6).Some elements in the vicinity of complex geometry details are smaller than the allowed minimum sizes presented in Table 2. Mass scaling will be used in these elements.Belytschko-Lin-Tsay shell elements are used for modeling the bump track and all other steel parts in the structure.The boundary conditions shown in Figure 6 are rigid which might lead to overly stiff ground response, but studying the stiffness of the ground was beyond the scope of this work.
Mass damping (LS-DYNA® keyword *DAMPING_PART_MASS [12]) with a magnitude resulting as 3 % of the critical damping of the lowest natural frequency of the bump track is used in the analysis.As the rear frame consists of moderately thin sheet metal plates assembled by welding it is suitable to be modeled using shell elements.Some thick plates in the vicinity of the rear axle violate the applicability range of the Mindlin-Reissner plate theory, which might lead to errors in the results for these parts although no excessive deformation is expected in them (see [19]).By using shell elements in all parts they can be directly joined together in the nodal level.If solid and shell parts are joined in the nodal level, shell rotations are left undetermined since a node of solid element has only translational degrees of freedom.To "weld" (to constrain all necessary DOFs) shell and solid parts together requires the use of some other constraint method than nodal connections (e.g.contacts).Consequently, to keep the model as simple as possible while still retaining an acceptable level of accuracy, all parts of the rear frame are meshed using the Belytschko-Lin-Tsay shell elements.Since all connecting parts are joined via shared nodes, possible reduced stiffness caused by welds with cross-sectional area smaller than of the plate's is ignored, i.e. the welds are assumed to have a stiffness equal to the base material.
It is assumed that the heat from welding does not cause warpage to the plates or any residual stresses.No misalignment of the plates during assembly is considered.That is, the geometry is assumed to be ideal; as designed.This is expected to have a very small influence to the dynamic behavior of the structure.Residuals or imperfections are not usually accounted for in these types of analyses [11,21].Initial imperfections (wavelike deformations in plate fields and plate misalignments) affect the buckling resistance of the structure, but the loads used in this simulation are assumed to be much lower than those resulting to buckling (global loss of load bearing capacity).
The chosen density for the element grid is a compromise between accuracy, required CFL time step, and computation time.The results should approach the continuum solution, and therefore mesh independency, as the element size approaches zero.However by decreasing element size, the CFL time step also decreases and the total number of DOFs for the model increases, both of which increase the total computation time.Sufficient element size was studied using undamped modal analysis (using explicit analysis to study mesh dependency is evaluated to be computationally too expensive).
Modal analysis solves eigenvalues for linearized form of Equation ( 1) when damping and external force terms equal to zero.These eigenvalues represent the modal frequencies (or eigenfrequencies, natural frequencies) of the structure.Modal shape vectors can be calculated based on modal frequencies [15].Modal shape vectors describe the shape in which the structure sinusoidally oscillates around its equilibrium position.The solver used in calculating the modal solution is ANSYS® Mechanical PCG (Preconditioned Conjugate Gradient) Lanczos Eigensolver for which the Belytschko-Lin-Tsay shell element is not supported.The SHELL181 element with reduced integration is used instead since it is closest to Belytschko-Lin-Tsay within elements supported by PCG Lanczos Eigensolver [1].Modal analysis is conducted with different element sizes to obtain eigenfrequencies as a function of element density.The results shown in Figure 8 are considered to be mesh independent if eigenfrequencies do not change (within a tolerance) with further mesh refinement.If this tolerance is chosen to be 3 %, sufficient global element size is 25 mm, resulting to 21 710 elements in the model of the rear frame when studying first three non-rigid eigenmodes shown in Figure 7.The chosen element size is highlighted with red in Figure 8.The correct modeling of tire stiffness behavior is critical in this analysis since the forwarder does not have a separate suspension system.Tires are therefore a major contributor in the total dynamic behavior of the machine.Tire models used in some heavy machine simulations [11,21] have been coarse since the stiffness behavior of the machine has been dominated by the suspension system (springs and shock absorbers).
The tire model used in the simulation of the entire machine is the one that is evaluated to have the best accuracy to computational efficiency ratio.This ratio is evaluated by simulating a single tire and varying critical parameters.
The mesh dependency of the results was studied by compressing the tire against a rigid wall using different element densities in the model.The results in Table 3 show that the deflection is mesh dependent even when using the finest studied mesh (even though a converging trend can be observed).The computation time for the fine mesh model is more than ten times longer due to the higher number of elements and smaller

Number of elements [thousands]
Torsion mode 1st bending mode 2nd bending mode CFL time step compared to coarse mesh model.Therefore the use of high mesh density is not practical.Since the tire cross-section layup details are unknown, they have to be modified to iteratively find the correct stiffness for the tire compared to experimental results.This allows the use of relative coarse mesh since the artificial stiffness introduced by it can be compensated through altering cross-section properties.Rubber is used as a matrix material in the composite layup (Figure 10) and it is assumed to possess isotropic material properties.Since the material behaves in a nonlinear manner (in the elastic region), simple constitutive relations (e.g.Hooke's Law) cannot represent the stress-strain relation correctly.In addition to this, Hooke's law is numerically incapable of modeling incompressible ( = 0.5) materials like rubber.Rubber stiffness properties also exhibit strong relation to temperature and strain rate which further complicate the implementation of universally correctly behaving material model [3,7].Despite these non-linearities, studies by Li et al. [11] and Szurgott et al. [21] have used a linear material model for rubber due to the fact that hyperelasticity is not supported for shell elements with layered composite properties (material properties for linear rubber are shown in Table 2).In this study the tire cross sections are modeled using shell elements and therefore hyperelastic relations cannot be used.The stiffness behavior of the composite layup is dominated by the reinforcements and therefore the error caused by the use of a linear model for rubber should not be excessive.
Because none of the material models discussed above are able to simulate hysteretic energy dissipation, some other form of damping must be included to prevent excessive non-physical vibrations.Relative damping (LS-DYNA® keyword *DAMPING_RELATIVE [12]) is used for this purpose due to its ability in only damping vibrations relative to rigid body motions.By using this method the rigid motions themselves are not damped which is a necessary property in simulating a rolling tire.This method for damping is used in recent analyses by Reid et al. [14] and Shiraishi et al. [17].The rigid body for which the damping will be relative to in each tire is the rim.
A correct damping magnitude for rubber is evaluated to be 5 % of critical damping.Relative damping requires two different inputs; the frequency to be damped (FREQ) and the fraction of critical damping (CDAMP).The selected frequency is recommended to be close to the lowest eigenfrequency [12 p. 1027].
Due to material symmetry, only two material constants are needed for isotropic materials.For orthotropic material, however, there exist three orthogonal symmetry planes, and a total of nine independent material constants need to be defined [10 p. 85].The reinforcement layer in the tire is modeled using an orthotropic linear material model.The stiffness of the reinforcement layer in fiber direction is assumed to be dominated by steel wires.In direction perpendicular to the fiber, the stiffness is dominated by rubber matrix (the stiffness is set two magnitudes larger in comparison to rubber properties in Table 2 to avoid tire stability issues encountered using lower stiffness).As this is a 3D analysis, stiffness in direction normal to the laminate surface must also be defined.In addition to the elastic moduli, shear moduli must also be defined for three main directions as there is no similar connection as for isotropic material.The shear moduli for both directions orthogonal to the fiber are assumed to be dominated by steel properties and therefore the values to be used are adopted from isotropic steel.The shear modulus in the fiber direction is estimated to be dominated by isotropic rubber properties.It was observed, however, that setting the modulus equal to the shear modulus calculated for linear rubber resulted as excessive shear deformations and tire instability.Therefore, instead of choosing the value for rubber, a value of steel was also set for the fiber direction.Six Poisson's ratios for the orthotropic layer need to be defined.The major ratios in this case are  12 ,  13 and  32 .The minor ratios  21 ,  31 and  23 can be calculated based on major ratios as [10 p. 95].The most dominant value of this layer is  1 .Other values should not crucially influence the results (other than tire stability).The above material properties are summarized in Table 4.To include the effect of pressure change during compression, LS-DYNA® *AIRBAG_SIMPLE_PRESSURE_VOLUME (this feature is called airbag for historical reasons; it was initially developed for automotive industry to model airbags) model is used.This model computes the volume change of the enclosed "airbag" (in this case the volume of the tire) which is inversely proportional to the applied pressure [12 p. 142].
The layup of shell elements in the tire cross section consists of a layer of isotropic linear elastic and homogeneous rubber with reinforcement layers in different orientations as shown in Figure 10 and Table 5.The reinforcement layers are thought to be comprised of steel wires surrounded by rubber matrix.It must be emphasized that the displayed layup is not necessarily similar to the layup of the real tire; it is rather the result of iterative process.Table 5. Properties of tire layers.Fiber orientations shown with respect to default orientation (Figure 9).

Layer Material Fiber orientation [˚] Layer thickness [mm]
A Linear rubber (see page 86 It was found through preliminary analyses that the under-integrated Belytschko-Lin-Tsay shell was unstable in modeling the tire, and therefore it cannot be used.A more robust fully integrated shell element (LS-DYNA® type 16 [13 pp. 169-173]) is used instead.This element is also reported to be used in studies by Barsotti [4] and Shokouhfar et al. [18].Due to full integration, hourglass modes will not be present in the tire with the exception of transverse shear modes that will be damped using hourglass control type 8 with a coefficient of 0.1 [6; 12 p. 1506].This method of hourglass damping was also used by Barsotti [4].To further improve stability of rotating components, second order objective stress update (based on Jaumann's stress rate) is activated (*CONTROL_ACCURACY keyword flag OSU=1 (Objective Stress Updates)) as was also done by Barsotti [4] and Reid et al. [14].
A series of static radii results for different compressive forces and inflation pressures were obtained from Nokian Tyres [22].The uncompressed static radius was however not provided; it was read from a CAD (Computer Aided Design)-model to be 615 mm.Vertical tire deflections shown in Figure 11 with polynomial trend lines are calculated based on these static radii.
To obtain correct cross section properties, an iterative process of comparing the simulation results to the measurements [22] and then modifying the cross-section properties was conducted.The final results (with optimized layer thicknesses shown in Table 5) are shown in Figure 12.The cross section properties are adjusted to obtain good agreement with nominal pressure (500 kPa) and moderate loads.The results show overly soft behavior with underinflated tire (400 kPa) when compared to the measured results.This might be caused by severe bending moments experienced by tire (support from internal pressure is lower) elements connected to the rim.Although comparing simulation model results to a dynamic test results (such as the cleat test) would have also been important, no data were available for this.The damping properties of the tire would have been acquired from dynamic test results.
The friction coefficients in the contact between tires and the test track (LS-DYNA® keyword *CONTACT_AUTOMATIC_SURFACE_TO_SURFACE [12]) vary with respect to moisture level and other ambient conditions.Therefore the coefficients are only estimates.The friction is calculated for tire-track contact as where   is the relative velocity between contacting surfaces.

Results
All results from the simulation are recorded at intervals matching to the sample rate of measured quantities.Therefore, the phenomena that have a timescale shorter than this are not captured.
The computing time required to simulate the 12 ton load case (duration 5 seconds) was approximately 10 hours using MPP (Massively Parallel Processing) LS-DYNA with 20 processor cores at 2.8 GHz.Other load cases were cheaper to compute.
The simulated and measured results are plotted to same graphs to be compared for validation.The radial forces of right hand side tires are also decomposed in vertical (normal to ground) and horizontal (parallel to driving direction) components to better understand the radial results.
Immediately it can be seen that the simulated results are more "restless" on high frequencies (small short duration oscillations) compared to measured results.This suggests that the damping values used on high frequencies should be higher on the simulation model.
As can be noted from Figure 14, the horizontal force peaks on right hand side tires (also holds true for left hand side tires) are opposite in sign, meaning that the resultant horizontal force of the tire pair, shown in Figure 15 (sum of RF and RR forces), does not have oscillations as high as individual tires do.The first horizontal peak at RF tire occurs at  ≈ 0.7 .This is, caused by the tire hitting bump no. 1.An opposite reaction for this peak with a smaller magnitude can be seen in RR tire.
The second force peak occurs at  ≈ 2.7 , when both right-hand-side tires are at bump no. 1 (see Figure 13).These peaks are larger in simulation model than they are in reality.It is estimated that this is caused by the constant rotational velocity applied on the rims.In the physical machine the tires are not forced to rotate at exactly same velocity at all times.The use of a differential and a torque limiter might improve results in this area.The friction between the test track and tires also plays an important role in the horizontal forces; differences can be explained to some level by a possibly different tire grip on field test and simulation.One more possible explanation for the second force peak is the stiffness of the fixed boundary condition used for the test track (see Figure 6); in reality the ground response might be softer, resulting as lower force peak.
With the exception of difference between second force peak results, the simulated horizontal forces correlate well with the measured values.Vertical forces shown in Figure 16 also experience a peak force when RF tire collides to bump no. 1.The machine rear frame tilts around its longitudal axis as right hand side tires cross the bump (1  <  < 2 ).This tilting causes the rear frame CoG to shift towards the left hand side tires.This can be seen as reduced vertical reactions on the right hand side.Similarly the reactions increase on the left hand side.A high peak is observed at  ≈ 2.7 .It is estimated to be caused by same reasons as the horizontal peak (see page 90).Figures [17][18] show the radial reaction forces acting on the tires (see [19] for results of all tires).Radial force is the resultant of horizontal and vertical force vectors.The correlation is reasonable with the exception of the force peaks discussed earlier.Some temporal deviation between simulated and measured force peaks are observed at  > 3 .This can be interpreted as the physical machine slowing down since the measured peaks occur later than the simulated one.The simulated peak of LF tire hitting bump no. 2 is higher than the measured value.
Figure 19 shows the time history of maximum von Mises stresses for the rear frame.It can be seen that most of the highest stresses are caused by forces due to weight of the cargo.Stress peak is reached just before applying initial velocity ( ≈ −0.05 ).The high damping at this phase however dissipates some of the energy and the maximum stress decreases to a "quasi-static" state (0.15  <  < 0.45 ) as the machine rear frame is moving on a flat surface.Just before the first impact (rear frame RF tire to bump 1, see Figure 3 and Figure 1) the stress level plunges and rises again at impact ( ≈ 0.7 ).The stress states of these two points (illustrated in Figure 19) are shown in Figure 20.The most critical areas are the load space mounting points and the vicinity of the rear axle.

Summary
The goal of creating a simulation model capable of replicating the forces measured from the physical test was achieved with reasonable accuracy.The simulated tire forces correlated well with measured values with only one major deviance at ~2.7 seconds (see Figures 14 and 16 for explanation).The correct evaluation of the largest force peaks is important since they contribute most to the fatigue life of the machine components.The measured and simulated force peak results and their comparison is shown in Table 6.
For example, the difference between simulated and measured values of RF tire (see Figure 4) hitting bump no. 1 (see Figure 1) is approximately 5 % (shown in Figure 17).The forces on LF tire, when hitting bump no. 2, are overestimated.The stiffness properties of the tires were found to be a dominant factor in the dynamic behavior of the entire machine.Test results for tire dynamic properties were unavailable and therefore the validity of the used tire model in terms of dynamic properties could not be evaluated.In addition to tire stiffness, the precision of modeling the torque transmitted to tires via the power train was found to be a major contributor in the resulting tire forces.Tire friction depends on ambient conditions, thus the values of friction coefficients have to be more or less "guessed".A sensitivity analysis (see [19]) showed that the magnitudes of the force peaks are sensitive to the used values of friction coefficients.
The stress levels of the rear frame and load space were observed to be governed by the forces exerted from the weight of the cargo.The rear frame stresses had larger oscillations due to the tire forces compared to the stresses on the load space.Since shell elements with coarse discretization were used, the results of local stress concentrations are not accurate.

Figure 1 .
Figure 1.Bump track general view.Driving direction shown by red arrow.Bumps are numbered to be referenced later in the results.

Figure 2 .
Figure 2. Rear frame wheel radial reaction forces measured at the test track with 12 ton load.

Figure 3 .
Figure 3. Results zoomed to crossing of the first bump track.

Figure 5 .
Figure 5. Mesh for the entire simulation model has approximately 151 thousand elements.The initial position of the machine is shown.

Figure 6 .
Figure 6.All DOFs are fixed at the highlighted edges (blue) of bump track.The structure has RHS tubes as reinforcement at these edges even though they are not visible in Figure 1.
a) Torsion mode b) 1 st bending mode c) 2 nd bending mode Figure7.First three non-rigid eigenmodes of rear frame.

Figure 8 .
Figure 8. Converging behavior of natural modes shown in Figure 7.

Figure 9 .
Figure 9. Default fiber orientation is tangential to the red line.

Figure 10 .
Figure 10.Layers of tire laminate (layer thicknesses are out of scale for better visualization).

Figure 11 .
Figure 11.Measured static radii and trend lines crossing at origin.

Figure 12 .
Figure 12.Comparison of simulation results to measured data show satisfactory agreement.

Figure 13 .
Figure 13.The right-hand-side tires of the rear frame crossing bump no. 1 at  ≈ . .

Figure 14 .
Figure 14.The horizontal forces of the right-hand-side tires of the rear-frame.

Figure 15 .
Figure 15.Sum of horizontal forces of right hand side tires (RF and RR in Figure 4).

Figure 16 .
Figure 16.The vertical forces of the right-hand-side tires of the rear-frame.

Figure 17 .
Figure 17.RF tire radial reactions.Simulated radial peak force of RF tire hitting bump no. 1 is 5 % larger than the measured one.

Figure 18 .
Figure 18.LF tire radial reactions.Simulated radial peak force of LF tire hitting bump no.2is 27 % larger than the measured one.

Figure 19 .
Figure 19.Maximum von Mises stress time history of the rear frame.

Figure 20 .
Figure 20.Von Mises stresses of the rear frame before first impact and at the impact (see Figure19for notation).

Table 3 .
Results of mesh convergence study.

Table 6 .
Comparison of measured and simulated radial force peaks [kN].