Modelling and simulation of an harvester head mechanism

– This paper presents an approach for modelling and simulating a harvester head mechanism in the delimbing process of hardwood stands. An original model of the mechanical transmission is developed which allows, in the early phases of design, sizing of actuators and structural elements but also the validation and optimization of the various stages of the logging operation. The considered type of trees is characterized by ﬂexuous trunks and big branches. Hence, the current harvester heads are not enough eﬃcient to process such trees. New designs of harvester head are required to overcome these problems. In order to predict the process performances associated to several novel architectures of harvester heads, in preliminary design stage, a modelling of the process is developed. It is based on parameterized tree models with realistic geometry. The global harvesting head model makes interact original sub-models representative of basics functions as feeding trunk into the head and cutting branches by knife. These models identiﬁed experimentally, allow a good prediction of the interactions between the diﬀerent parts of the machine and trees. The modelling and simulation of harvester head mechanisms are implemented in a numerical multibody environment. Simulation results of the global model are in a good agreement with the observed behaviour of the real system and represent a design support tool to develop innovative concepts of harvester head.


Introduction
This paper presents researches carried-out within the ECOMEF project (Eco-design of mechanized equipment for hardwood harvesting) framework. This project aims to develop a harvester head more adapted to process and fell broad-leaved trees, in order to face the lack of motor manual workers and the low rate of mechanization in hardwood logging operations [1]. The processing operation is composed of delimbing trunks and bucking them into logs [2] (Fig. 1a). Due to crooked trunks, bigger branches with sharper angles and hardwood, broad-leaved trees are indeed harder than conifers to feed and delimb. Therefore, specific tools to improve the energy efficiency and the productivity of the operation are required [3]. This work focuses particularly on the modelling and simulation of the single-grip harvesters used in the cut-tolength method (CTL), like the Kesla 25 RH (Kesla Oyj 2013) harvester head (Fig. 1). The processing operation is composed of delimbing trunks and bucking them into logs (Ćuprić and Bajić 2009).
In the existing harvester head mechanisms (Fig. 1b), the action of the mobile knives (2 & 2 and 5 & 5 ) cona Corresponding author: vincent.gagnol@sigma-clermont.fr strains the trunk, by surrounding and tackling it against the support (0) and the fixed knife (1), to translate relatively to the support under the action of the feeding rollers (3 & 3 ). Thanks to the straight shape of the trunks in conifers, this motion is rapid and continuous and knives spacing self-adapts to the variations of the trunk diameter. The delimbing is produced by the impact of the knives' sharp edges against the branches. However, the processing of crooked trunks may fail during hardwood harvesting. Therefore, the ECOMEF project has addressed the challenging task of designing new mechanisms inside which crooked trunks can slide while ensuring power transmission and avoiding locking phenomena: bitting, over-center locking. In order to predict, from the preliminary design stages, the performances of these mechanisms, in terms of energy efficiency and productivity, an approach based on process modelling and simulation is implemented. Process simulations should enable testing and comparing several designs of harvester mechanisms.
A top-down functional analysis is presented in Figure 2 to describe the functional interaction inside the harvester head. Four functions have been identified: modelling the trunk gripping (function A1), feeding (function A2), delimbing (function A3) and felling and bucking  (function A4). In this representation, a closed loop interaction can be observed between feeding and delimbing model function in order to achieve the global function of harvesting trees. The presented paper is focused on the modelling of this closed loop to ensure the overall harvesting trees function. Given the proposed functional decomposition, three models are necessary to be developed: a driving model of the trunk inside the harvester head and a specific cutting model representative of the impact of the delimbing knives against the branches. It is also necessary to model the specific geometry of crooked trunks. The presentation of the different models and their integration into a global multibody system model organise the paper. Hence, the Section 2 presents the modelling of the building block, used to describe the model of driving the trunk in the harvester head mechanism both analytically and experimentally. The Section 3 is dedicated to the cutting model development. Two approaches are conducted simultaneously: the development of analytical approaches and the development of experimental test campaigns for the registration and / or identification of model parameters. In Section 4, an existing method for geometric modelling of trees is enriched in order to integrate trunk flexuosity. Finally, Section 5 presents the integration of the elementary modeling blocks into a global numerical multibody environment and simulations results.

Trunk feeding model
Cutting the branches off the trunk is usually achieved by feeding the tree at 3 to 7 m.s −1 using feed rollers, and creating an impact of the branches against delimbing knives. Cutting and feeding models require to be integrated in the dynamic modelling of the delimbing process. Most of the harvesting heads use rollers for trunk feeding. As explained in Reference [10], concentric gripping mechanisms are not very suitable for small stems, whereas lateral gripping mechanisms facilitates passing of crooked parts of the trunk but could lack pressing force when processing large trunks. Several hybrid solutions adapting the kinematics of the gripping motion to the diameter of the trunk, in order to be efficient in the whole range of considered diameters have been proposed [10].

Analytical approach
The aim is to develop a model to calculate the maximum drive forces transmitted to the trunk by the rollers, in order to size the various components of the lumberjack head. A geometric definition of the elliptical spikes is proposed (Fig. 3), from which equivalent surfaces are defined corresponding to the projected areas of portions of the roll which entered the trunk in the gripping direction (perpendicular to the fibers) and advance (in the sense of fibers). A relationship is established between the projected surfaces and maximum forces transmitted through the maximum allowable yield strength in each of the directions. In order to integrate these models in a more comprehensive global delimbing simulation (detailed in Sect. 4.2), an equivalent coefficient of adhesion c adh representative of the roller / trunk contact is established. The adopted approach to compute feeding forces is similar to that used by Klamecki [11] to calculate the friction coefficient due to the creation of a furrow by a conical asperity travelling through the chip (Fig. 4a). Klamecki expresses the transverse and frontal forces on the asperity using the simple relation F = σ · S, F being the force and S the projected area in the considered direction. N A represents the number of conical asperity in charge.
For the transverse force T , the considered surface At is half a disk of diameter d (maximum penetration of the cone into the chip, and the stress linking force and area is the yield strength of the chip material. For the frontal force N , the Af area represents the projection of the cone in a radial plane, perpendicular to the force (triangle of the height of x and of the base of d). The considered stress is associated with the strength of the chip material.
The originality of the proposed models is to define a gripping and a feeding models representative of the contact between flexuous trunk and training roller equipped with spikes. The proposed models are adaptable to the wide variety of existing geometry and configurations. First, an analytical model allows the calculation of the penetration depth and the projected area normal to the feeding direction, parallel to fibres. Then, the maximum stress admissible by a given wood species in the gripping direction, perpendicular to fibers, needs to be experimen- Klamecki model of force generated by a conical asperity [11].
tally determined, in order to link the measured gripping force with the projected area normal to this direction. Finally, the maximum stress admissible in this direction is also determined experimentally and enables calculating the maximum transmissible feeding force.
Still being one of the most used types of stud on feed rollers, the elliptical stud or steel spike could be geometrically obtained by two successive oblique cuts in a rod (Fig. 5). The proposed model allows calculating the areas when using elliptical studs. To extend the study to other shapes of studs, derived analytical formulations or a CAD-based calculation method are proposed to help engineers to optimize the feeding process: maximizing the maximum transmissible feeding force by adapting the roller geometry as well as the shape and the pattern of the studs. The Figure 5 shows the development of equivalent frontal surfaces Af and cross surface At that are established from the proposed analytical formulation for the entire roll range of the manufacturer Mense Oy (2013). The geometrical designation used for the studs corresponds to the dimensions Φ S × h S in millimeters unit.
The calculation of surfaces Af and At, considered as portions of ellipses, is necessary for the determination of maximum allowable stresses and then transmissible force in each direction. Detailed formulas presented in Equation (1) to (2) allow the calculation of these surfaces through the expression of the parameters a, a and b shown in Figure 5, based on the geometric parameters definition of an elliptical spike.
To calculate the frontal area Af, the elliptical portion is considered that expressed in Equation (1), to which a rectangular portion can be added depending on the height of the spike (if h S − h a f , as presented in Figure 6a (b). In order to consider each front surface in a plane normal to the feed direction rather than a radial plane of the drive roller, a projection θ r taking into account the rotation of the roller and each spike row is necessary.
a f a f . The calculation of the transverse surface At requires the determination of clearance on each spike. Considering a staggered arrangement of spikes between two successive rows (Fig. 5), it is possible to define every free height according to the clearance of a reference spike, defined as the central spike of a row consisting of an odd number of spikes.
as presented in Figure 6b (b).

Experimental calibration
Two experimental test bench are designed to characterise experimentally the values of the maximum stress admissible by a given wood species in the gripping direction (perpendicular to fibres) and in the feeding direction (parallel to fibres). These measured admissible stresses respectively denoted σ ⊥ and σ // allow tuning the proposed models in gripping and feeding direction. In order to integrate the high variability of the measured maximum permissible stresses, a probabilistic approach is developed. Then, an equivalent sticking coefficient representative of   the roller/trunk contact is then determined from the proposed models.
The tests presented concern the experimental characterization of the maximal stresses perpendicular to the fibers σ ⊥ and in the direction of the fibers σ // . Figure 7 shows the test bench in the clamping direction designed and built specifically and a standard curve of the clamping force versus displacement is presented. Experimental characterization are carried out on a traction/compression machine (maximum force 250 kN), the upper jaw portion is provided with a roller and the lower jaw of a support for maintaining the trunk. In order to refine the determination of the penetration depth, the crushing of the trunk under the effect of clamping is considered. The positioning of the sensor (inductive sensor plunging core) for measuring the crushing of the trunk is shown in Figure 7. The parameters tested and the associated levels are detailed in Table 1.
The average maximum allowable stress is calculated on the basis of the equivalent surface At resulting from the   stress perpendicular to the fibers σ ⊥ average value is evaluated to 22.14 MPa.
A second test bench, designed and built specifically for testing in the feeding direction is presented in Figure 8. It allows the simultaneous application of a clamping force of the roller against the trunk, adjustable and measurable, and a driving torque, applied by means of a hydraulic cylinder causing the rotation of the roller. The clamping forces Fs, the maximum transmissible drive effort Fe max and the maximum allowable stresses in the direction of the fibers are evaluated during each test. The tested parameters are similar to those shown in Table 1. The clamping force is the value recorded at the end of the roll against the trunk clamping phase (just before training).
The drive force Fe max is the maximum thrust force of the trunk, calculated by dividing the drive torque of the roller by the distance separating the axis of rotation of the roller and the center of the front surface. The maximum allowable stress in the direction of the fibers is σ // determined for each test by dividing the maximum drive effort Fe max by the frontal surface At before training. The depth of penetration required for calculation of frontal surface is directly measured and readjusted as a function of the deformation of the trunk by means of inductive displacement sensors.
As in the clamping direction, a wide variability of experimental results are observed in the training direction which required applying a probabilistic approach to determined the value of the stress. The average value of σ // = 24.14 MPa is chosen as the reference mean value. This value is very close to the compressive strength in the direction of the fibers which is 24.4 MPa for fresh wood of Quercus michauxii referenced by Kretschmann (2010).

Equivalent adhesion coefficient
The proposed approach allows modelling the driving by obstacle of a trunk by feeding rollers equipped with spikes. In order to integrate this model in a more comprehensive simulation (detailed in Sect. 4.2), an equivalent coefficient of adhesion cadh is established representative of such training. Indeed, the software ADAMSTM to be used for the dynamic simulation of scrolling trunk offers by default the definition of coefficients of friction (static and dynamic) for sliding contacts. The clamping and feeding forces experimentally measured associated with the analytical model of gripping and feeding allows determining the equivalent coefficient of adhesion c adh for standard rollers configurations (Fig. 10).
The Figure 11 illustrates the calculated variation of the equivalent coefficient of adhesion c adh depending on the type of spikes, trunk diameter and roll diameter. This spikes design analysis towards equivalent coefficient of adhesion is a design support tool to optimize roller training efficiency.

Analytical formulation
Unlike many cutting studies based on the type and the section of a chip for several machining operations, it is not possible in the studied case to define a chip, since the "chip" is the whole branch. The lack of cutting model representative of delimbing operation in literature confirmed the need to develop an original predictive cutting force model. The geometry of the cutting situation is defined on the basis of McKenzie [5] approach. In this notation, two angles are used to define cutting process. The first one is formed by the cutting edge with the grain direction, the second one by the cutting speed direction with the grain direction. Delimbing can thus be considered as a 90 • -90 • orthogonal cutting process, sometimes called cross-cutting. In the same way, chipping can follow the same definition. The work from [6] is adapted to define the cutting situation, as presented in Figure 12.
Such a definition of cutting angles conserve the usual relation α + β + γ = 90. Based on image correlation analysis, performed in a preliminary study, the principles for the definition of a reduced model have been defined. The calculation of cutting force (Eq. (3)) is based on three components and specific geometric parameters presented in Figure 13. All these components depend on the penetration depth p. Constants k Sx , k Sz and k L have been identified experimentally by using two types of test benches for slow and rapid delimbing [7]. They depend on the humidity ratio T H of the wood. (3) underlines the importance of the projected surface A zeq (Fig. 12) that should be reduced to decrease the cutting force along the feeding direction z. This design rule has been used to develop an innovative ribbed knife published in Refference [8] and patented [9]. A xeq , A Zeq , L eq are determined analytically. For the calculation of the equivalent surfaces, each branch is represented by an elliptical base cylinder whose axes are Φ // and Φ ⊥ . The L eq formulation is: The equivalent surface A xeq orresponds to the projection in the cutting plane of the intersection between the branch and a L eq cutting edge length portion. This intersection is shown in Figure 12a, and its projection in the section plane is shown in Figure 12b. The definition of a branch by an elliptical cylinder basis allows describing the equivalent surface A xeq as a portion of an ellipse (Fig. 14). By considering half of an ellipse centered on the origin point O and whose semi-major axis a and semi-minor axis b, another portion may be delimited by defining a vertical line with abscissa x = a (Fig. 14). The equivalent surface A Zeq is modeled by a rectangle of length L eq and whose height corresponds to the thickness of the blade at a distance p of the cutting edge (with p the depth of penetration).
Experimental validation is documented in Reference [7]. It allows determining through a probabilistic approach the coefficient of Equation (3) as presented in Table 2.   By defining

Models integration and delimbing simulation
In order to simulate the global process, feeding and cutting models, presented in the previous sections, are integrated in a multi-body model of the harvester head mechanism. In these simulations, cutting and feeding forces are generated as function of the relative motions of the trunk relatively to the knives and the rollers. Phenomena affecting process efficiency, such as over-centre locking and bitting, can be reproduced in order to compare several design solutions for a given geometry of the tree. A geometric model of crooked trees is developed in order to integrate trunk flexuosity.

Crooked trunks geometric modelling
Hardwood species are characterized by the sinuosity of their trunks. In this section we focus on the geometric description of this type of trees. In order to evaluate and compare by simulation the performances of different harvesting head designs, it is useful to generate realistic tree models. Therefore, a number of representative trunks will then allow the comparison of simulations while avoiding the morphological changes observed during tests in operating conditions. This work is carried out before the simulation step.
The geometric trees model can be adapted to several cases, in which the contact can be modeled by planar or 3D kinematic constraints, or more complex dynamic studies with contact handling based on a collision detection algorithm. Generated trunks must incorporate representative breaks and discontinuities of more or less crooked trees as well as the first level of branches in order to simulate delimbing operations. The complexity of vision-based reconstruction methods as well as the difficulty of selecting representative samples of the great morphological diversity of trees led us to move towards a parameterized model, especially the model of Weber and Penn [4]. This type of model offers an important versatility regarding the considered simulation context. It offers the following possibilities: -Recursive model: each child stem inherit its properties from its parent. -More than 80 parameters defining the general shape, and each level of stem. Spatial developments and aleas in the profile curve of each trunk are introduced in the Weber and Penn's model in order to make it capable of generating irregular crooked (flexuous) trunks. The morphological characteristics specific to species considered are also integrated by establishing correlations between several key parameters definition of trees. Introduction of additional variables enables adjusting the degree of flexuosity (Fig. 15).

Dynamic simulation of an harvester head during delimbing
The approach adopted to integrate feeding model consisted in using equivalent sticking coefficient for the contact modelling between the trunk and the rollers. This coefficient has been determined experimentally. Nevertheless, analytic calculation of surfaces offers the possibility to determine such a coefficient for different dimensions of the rollers. Indeed, this approach requires the numerical coupling and co-simulation between ADAMS r , MATLAB r , OpenSCAD and Blender based on Boolean intersection evaluations and the use of analytical expressions of the equivalent surfaces used for force estimations.
In the multibody model presented (Fig. 16a), the solids are assumed to be undeformable and only the rigid body motions are studied. Based on the kinematic model of a harvester head presented in Figure 16b, various subassemblies have been identified for the development of a model set of articulated frame: -A sub-assembly "upper knives" composed of two mo-45 bile upper knives 2+2, the fixed knife and a part of the frame called top frame 1 and finally the mobile knives clamping mechanism (not shown in the cinematic diagrams); -A sub-assembly "drive roller" consisting of the central part of the chassis 1, as well as two rollers 3 + 3 , the roller support arm 4 + 4 and the associated clamping mechanism (not shown); -A sub-assembly "lower knives" compound bottom frame 1", the clamping mechanism (not shown) and two lower mobile knives 5 + 5 . Process simulations have been performed on a standard trunk in a one-piece frame (constituting the reference architecture) as well as on articulated frame, in order to compare these two architectures. The curves generated for the different parameters validate the proper operation of the simulation. The proposed delimbing model allows predicting the evolution of the branch cutting effort along the trunk advance. Its use requires knowledge of the evolution of frontal and transverse equivalent areas A xeq and A zeq and the equivalent length Leq in order to be able to determine the cutting force. The differents components necessary for its implementation in ADAMS TM are developed for a simple case characterized by a trunk and a perpendicular branch, both cylindrical and a rectilinear movement of the trunk. The cutting forces relative to the delimbing simulation for the reference architecture are presented in Figure 17a. This helps locate the contact with each branch on a time scale.
During the simulation, the roller engines are controlled in speed. It induces torque peaks (Fig. 17c) when the knives enter in contact with the branches that are more important than what the system is actually able to provide (some observed peaks exceed 5000 Nm, while the maximum values provided from the manufacturers is 2500 Nm per roller). Indeed, assuming an average delimbing force of 28 kN and roller diameter of 350 mm, one can determine a driving torque of about 2500 Nm per roller. This behavior is also observed in Figure 16b on the power peaks reaching 130 kW (against 75 to 100 kW total power announced). It is important to note that the modeled system will never block in a branch by lack of power. The user must therefore be extra attentive to power consumed. The total power drawn in Figure 16b shows an efficient value 41 kW. This is in agreement with the data values for the harvester head Kesla 25RH Kesla Oyj (2013), since a required power of 75 kW to 100 kW is announced by the manufacturer.

Conclusions
In this paper, an original modelling approach for simulating the delimbing process of hardwood stands has been presented. An existing method for generating realistic tree geometric model has been developed. Parameterized trees with irregular and crooked trunks provided with branches have been generated.
In order to improve the realism of the operation, the developed models of feeding trunk into the head and of cutting branches by knife have been integrated into a dynamic multibody simulation implemented on ADAMS r software. The integration of the equivalent coefficient of adhesion resulting from the feeding model, and cutting forces whose evolution depends on the contact surfaces between the trunk and knives and the equivalent engaged tool length (calculated based on the depth of penetration) have been introduced in the proposed closed loop delimbing model. Its setting allows many dimensional variations of the head as well as the modification of the kinematics architecture. The scrolling of a standard trunk in a reference architecture of harvester head has been simulated and compared to manufacturers data. Hence, the proposed model appears as an efficient tool for predicting and comparing the performances of several novel harvester head mechanisms from preliminary design stages. Future work will focus on enhancing simulation realism in order to reproduce precisely bitting and over-center locking phenomena and optimize the energetic efficiency of the proposed design solutions.