Linear woodcutting of European beech: experiments and computations

Hardwood species are becoming increasingly important with the growing need for a diversity of forests that have recently been facing global temperature changes or conifer pests. This further leads to the growth of its potential as a building material that may originate from sustainable production. As hardwoods need to be properly processed, the article deals with the disintegration of European beech. The influence of wood grain direction, uncut chip thickness and cutting speed on the cutting force magnitudes was experimentally investigated using the device with a rotating arm of approximately 4 m in diameter. Then, the disintegration process was modelled using the finite element method in Abaqus/Explicit. The developed material model consisting of orthotropic elasticity and plasticity with rate-independent and rate-dependent tensile–compressive failure asymmetry was implemented through the user subroutine, while the crack initiation and propagation were realized using the element deletion technique. The computationally predicted average values of cutting forces and chip shapes were, except for a few tests, in good agreement with the experiments. It means that the model may be used for further investigation, such as the influence of tool wear.


Introduction
There is a broad spectrum of applications for hardwoods, whose presence is increasing in European forests. Hardwood is utilized in wood-based composites, construction of buildings, furniture industry and others. The material has to be partitioned or treated in the vast majority of wood processes. Therefore, there is an effort to improve the mechanical disintegration techniques of hardwoods within the industry. Numerical modelling may decrease the number of physical tests, and therefore, the process may be optimized in an effective way that may help to protect the environment. Better surface quality may be achieved or tool wear may be decreased with energy intensity. Nevertheless, it should be noted that the computations are not intended to completely replace the experiments. Kivimaa (1950) was one of the first researchers to study the influence of fundamental parameters (wood density, moisture content, temperature, tool geometry and cutting speed) on the cutting force. Kivimaa (1950) found that the wood grain direction is most influential, which was confirmed in this study. As in this work, Franz (1956) focussed on linear wood machining and followed the depth of cut, rake angle and wood feeding influence on the force response and surface quality, and further analysed cutting parallel to wood grains and derived an analytical equation for controlling the surface quality based on optimal cutting angles. Circular wood cutting was analytically investigated by McKenzie (1961), who measured the forces on an adapted lathe. An analytical approach was also sought by Axelsson et al. (1993) and Porankiewicz et al. (2011), who focussed on the dependence of various process parameters (cutting speed, moisture content and density) and analysed it using statistical data. All the mentioned references chose an analytical approach, which is in contrast to this work that utilizes a numerical approach in the sense of the finite element method, which might be broadly applicable. Costes and Larricq (2002) conducted the numerical analysis (as in this work) of the surface roughness of milled beech with a numerically controlled routing machine using optics. The tests were carried out for constant chip thickness, while the cutting speed ranged from 3 to 62 m·s -1 . The results showed that increasing cutting speed leads to better surface quality, especially for dry wood cut parallel to wood grains. Similar research was later carried out for beech, oak and spruce by means of a sliding mitre saw with manual feeding (Kminiak and Gaff 2015). Hlásková et al. (2015) presented a model to determine the tooth cutting edge positions and cutting speed directions using the fracture toughness and shear yield stresses. The method was applied to beech wood using a sash gang saw on a circular sawing machine. Eyma et al. (2004) studied the most important factors influencing the force needed for the chip separation. Research was carried out for 14 species of wood in the peripheral milling in the fibre direction. The beech was among the wood studied as in the present work. Krenke et al. (2017Krenke et al. ( , 2018 used a modified pendulum in order to achieve conditions close to the linear cutting of wood. The spruce was investigated and a cutting speed of approximately 7 m·s -1 was reached. Then, the transfer function was used in order to filter the parasitic effects caused by the damping of the system and by the knife vibrations, while the function was based on the falling ball impulse. Finally, Pfeiffer et al. (2015) analysed the chip formation for an orthogonal cutting of beech and Douglas fir using the Chardin pendulum. This paragraph briefly summarizes works dealing mainly with linear and rotary machining processes. The pendulum-based device seems to be suitable for high-speed cutting close to the linear one. Therefore, such a device with a rotor that produces cutting speeds of up to 100 m·s -1 was developed with further details in Dvoracek et al. (2021a, b).
The finite element method (FEM) has been widely used for simulations of chip formation during machining, but mostly for metallic materials (Carroll III and Strenkowski 1988;Marusich and Ortiz 1995;Hua and Shivpuri 2004;Calamaz et al. 2008;Mahnama and Movahhedy 2010;Jagadesh and Samuel 2015;Paresi et al. 2021). FEM has been utilized for wood cutting simulations rarely. Most of the models have been based on theoretical, empirical and statistical data. Uhmeier and Persson (1997) created a two-dimensional finite element model (the one presented within this work is three-dimensional) for the prediction of the shape and thickness of the chip. The material model consisted of orthotropic elasticity, Hill plasticity (Hill 1948) and failure based on the fracture stresses (normal and shear) and energies in tension and shear. The results were compared with the experiments in the literature (Uhmeier and Persson 1997). The failure was realized through cohesive zones. This approach requires the use of cohesive elements where the fracture is expected. This is a limitation compared to the present study, where the fracture might initiate in any finite element where the criterion is reached. Nairn (2016) used the material point method-similar to FEM-for the simulation of chip formation for Douglas fir wood. Orthotropic elasticity and Hill plasticity (Hill 1948) were used, while fracture characteristics were implemented using a thin layer between the specimen and the presumed chip on the basis of the cohesive zone. The model was first verified towards the analytical method, and then the influence of tool sharpness, grinding angle, friction and position of the chip breaker and base plate on the chip formation was studied. Pichler et al. (2018) designed and assembled a small machine based on a rotating arm for the investigation of the cutting processes of spruce and European beech. The knife holder was equipped with a strain gauge that measures the bending stresses that occur during the disintegration process. Then, the process was simulated using FEM within Abaqus/Explicit. The fracture was realized through element deletion as in the present work. Nevertheless, the simplification was that the orthotropic behaviour was simulated by a variation of isotropic material parameters with respect to the cutting and wood grain directions. Therefore, such a model is not entirely complex as in the present study. Meng et al. (2019) conducted experiments on the small cutting machine with a circular saw blade for mulberry, while the cutting forces were measured, and the surface quality of branches was investigated under various conditions (rotation speed, number of saw teeth, diameter of branches and cutting speed). Then, the simulations were performed within ANSYS/LS-DYNA, employing the orthotropic elasticity but giving no details about the plasticity and damage. Aboussafy and Guilbault (2021) presented a two-dimensional model for the prediction of the cutting forces that arise in linear woodcutting. The finite element model based on fracture mechanics (containing the orthotropic elasticity with Hill plasticity) needed to contain an initial crack that propagated due 1 3 to the constant remeshing during the simulation. The direction of crack propagation was based on stress intensity factors determined on the basis of the displacement extrapolation method. Computations were done for Douglas fir wood, and the chosen parameters like depth of cut, rake angle and cutting speed were followed. This study presents a different approach. In the present fracture analyses, no initial crack was required, nor the remeshing, which is at a high computational cost.
This article presents linear machining of European beech wood. The experiments were carried out on a large device with a rotating arm that has a radius of approximately 2 m, reaching a radial velocity of up to 100 m·s -1 at its end, which is rarely seen in the literature. The influence of depth of cut, cutting speed and wood grain direction (grain angle) was experimentally investigated while also numerically simulated using a model built within FEM in Abaqus/Explicit. The model covered orthotropic elasticity, plasticity as well as fracture. In addition, the tensile-compressive failure asymmetry was incorporated using stress triaxiality so that the force responses corresponded well to the measurements. As written in the previous paragraph, there are only a few studies dealing with woodcutting using the finite element method in such complexity. The presented model is readily applicable to other areas, such as crashworthiness studies (Müller et al. 2020).

Test device
The device for cutting force measurements consisted of a 4050-mm-diameter rotor arm. The length of the arm was chosen in order to allow the study of linear cutting processes using state-of-the-art cutting speeds. Although the process is not completely linear, the arm length makes it possible to assume. The greatest deviation from straightness was 0.62 mm along the cut length of 100 mm, which made the angular deviation proximately 0.7° (on 50 mm length as it was symmetric). The test rig is based on the pendulum principle, i.e. the specimen moves in a circle and the knife is stationary. Therefore, the cutting speed corresponds to the feeding speed. The specimen is attached to one end of the arm, while the counterweight that keeps the assembly in balance is on the opposite side. Modal and structural analyses were conducted in the design stage so that the structure had an optimal stiffness to weight ratio as it rotates at high speeds. In areas that vibrated above the norm, special reinforcements were made. To validate the proper dynamic behaviour of the structure, several accelerometers were mounted to the rotor arm. The device can achieve cutting speeds up to 100 m·s -1 , so the inner structure was enclosed for safety reasons (Fig. 1). The walls were made from 4-mm-thick steel plates covered with aluminium foam panels from the inside that served as impact absorbers in the unlikely event of an accident. A low-speed set-up mode is applied in the case of operations around the processing zone. Additionally, various safety sensors were implemented (Fig. 2) to stop the operation when anything moved or approached the zone around the arm or the knife holder. The steel knife holder was topologically optimized in order to have a high stiffness and thus a natural frequency (Dvoracek et al. 2021b). For the setup used in the study, a LEITZ turn blade knife (HW-05:30 × 12 × 1.5 mm) was attached to the knife holder at an angle of 60°. The knife had a clearance angle of 20°, a wedge (sharpness) angle of 57° and a rake angle of 13°. The knife holder was equipped with a pre-tensioned force sensor. Then, the assembly was attached to a BAZUS Kinematic PK-2 slide table moving only in the horizontal (radial) direction. The slide table moved forward a set distance in order to achieve the required (uncut) chip thickness (depth of cut) prior to the cutting and returned to the reference position after the operation. Force was measured with a three-dimensional piezoelectric Kistler 9027C force sensor and conditioned using a charge amplifier Kistler 5080A and 3 × 5067. The force was sampled with a highspeed data acquisition unit HBM Gen3t and GN815. Finally, the process was controlled by the National Instruments CompactRIO-9045 real-time system with an inhouse algorithm implemented within LabVIEW 2018 (Dvoracek et al. 2021b).
The experiments were also captured by an iX Cameras i-SPEED 726 high-speed monochromatic acquisition set with a frame rate of 200,000 fps. The exposure time of 13 ns required a Dedotec COOLH laser stroboscope as a light source. The recorded images served for visual comparison of the chip shapes with the simulated ones.

Material and samples
Wood is a natural composite with a complex structure that contains a large number of chemical compounds such as cellulose, hemicellulose, lignin and so forth. However, wood is considered an orthotropic homogeneous continuum in the vast majority of cases, so its mechanical characteristics are distinct and independent in three perpendicular directions. The direction parallel to the grain growth (tree axis) is designated as longitudinal (L), perpendicular to the grain pointing to the tree axis is designated as radial (R) and perpendicular to the grain, while tangent to the trunk is designated as tangential (T) direction. Schematic illustration of material directions with respect to grain growth is seen in Fig. 3.
The wood studied was European beech (Fagus sylvatica, L.) grown in forests near Brno (Czech Republic). The samples had a thickness of 20 mm, a width of 50 mm and a length (cutting direction) of 100 mm. Once manufactured, the specimens were stored under standard climatic conditions of 20 °C and 65% of relative humidity for Fig. 3 Illustration of the material directions with respect to grain growth one month prior to testing, resulting in an equilibrium moisture content of approximately 12%. The specimens were measured and weighed before testing, resulting in an average density of 660 kg·m -3 with a standard deviation of 16 kg·m -3 . The surfaces of the specimens to be cut were rounded using a small feed value of 0.05 mm in order to guarantee a smooth surface and, therefore, constant chip thicknesses during measurements (Fig. 4c). It should be noted that one sample was cut multiple times. To isolate one cut from another and eliminate the influence of unwanted factors, smoothing of the surface was performed between each cut (similarly after a change in the cutting speed, which causes the arm to shorten or lengthen due to centrifugal forces).

Measurements without rate dependency
The tests were carried out at a cutting speed of 90 m·s -1 in order to study the influence of wood grain direction and depth of cut. The chips were cut parallel to the grain in the longitudinal-radial (LR) plane and perpendicular to the grain in the radial-tangential (RT) plane (Fig. 4a, b). Two wood grain directions (LR and RT) and four depths of cut (0.1, 0.2, 0.3 and 0.4 mm) were considered, which made 8 unique conditions. Each condition was repeated 8 times. Therefore, a total of 64 runs were made.
The measured force responses were filtered using the transfer function, which was determined based on modal tests, so that parasitic effects, such as damping of the system and knife vibrations, caused by the initial impact of the specimen with the knife, were eliminated (Dvoracek et al. 2021a, b). The remaining oscillations could be attributed to the complex natural structure at the macroscopic level (mainly alternating different density of early and late wood zones of annual rings and relatively large soft radial rays) and at the microscopic level The force responses of linear woodcutting for the LR and RT wood grain directions, the cutting speed of 90 m·s -1 and the depth of cut of 0.1 mm are shown in Fig. 5. Other tests exhibited a similar scatter in results. Data were re-sampled in MATLAB R2019a in order to obtain the average force response for each condition (Fig. 6). The average force response increased with increasing depth of cut. Furthermore, the responses for LR wood grain direction were approximately twice lower than those for RT wood grain direction. Furthermore, the variability of force progression increased with increasing depth of cut. Finally, the force responses served for the estimation of the average cutting forces based on the trapezoidal integration, which is used for a further comparison with FEM.

Rate-independent model
Computations were performed within the Abaqus/Explicit commercial finite element code. Explicit time integration scheme is suitable for dynamic simulations containing discontinuities-cracks. These were initiated and propagated using the element deletion technique. The material model described below was not available in the software. Therefore, it was implemented using the Vectorized User MATerial (VUMAT) user subroutine. The elastic-plastic behaviour implemented within VUMAT was carefully checked towards Abaqus/Explicit using a cubic specimen subjected to various loadings. The relative errors in strains, stresses and forces were below 1%.
The knife was considered a rigid body. Therefore, it did not vibrate within the simulations. As it vibrated within experiments, the signal was filtrated as discussed in Subsection Measurements without rate dependency.
The sample (European beech) was modelled as a homogeneous orthotropic continuum. Orthotropic elasticity contained 9 elastic constants and related strains and stresses in the following form where i , i and E i are the normal strains, the normal stresses and the Young's moduli for i = L, R and T , respectively, while ij , ij , ij and G ij are the tensorial shear strains, shear stresses, Poisson's ratios and shear moduli for i, j = L, R and T , respectively. All the elastic constants given in Table 1 were taken from Šebek et al. (2020) where the same material was tested under similar conditions. Orthotropic plasticity was described by Hill's quadratic yield criterion (Hill 1948). The corresponding yield condition is written as where H is the Hill's equivalent stress and y is the actual reference yield stress, which was a linear function (Fig. 7) of the equivalent plastic strain p (cumulative) defined as where ̇ p is the plastic strain rate tensor and t is the time. Finally, Hill's equivalent stress reads (Hill 1948) where F , G , H , L , M and N are the material parameters, which may be expressed using the yield stresses in individual normal directions ̃i for i = L, R and T , the yield stresses in shear ̃i j for i, j = L, R and T , the initial reference yield stress 0 (Fig. 7) and the initial reference yield stress in shear 0 as follows

Fig. 7 Reference flow curve with initial reference yield stress
The relationship between the reference yield stress and the reference yield stress in shear is given by the von Mises yield criterion as The material parameters for Hill's plasticity were taken from Šebek et al. (2020), as in the case of elasticity, where the wood was tested in order to study the behaviour to allow for simulations in this article. Parameters are summarized in Table 2 where H < 0 , which was caused by a different behaviour in the L direction when compared to other directions ̃L >>̃R ∼̃T . Finally, the isotropic hardening was assumed with an associated flow rule.
The approach to fracture was chosen simpler-uncoupled-when compared to Šebek et al. (2020, 2021). Elastic-plastic behaviour is not influenced by damage in uncoupled models, contrary to coupled ones. The coupled approach exhibits the stress softening that leads to the loss of ellipticity, which results in the localization of field variables. This can be treated by adopting the non-local approach, which is unfortunately at a high computational cost, therefore, omitted in the present study. The evolutions of damage parameters in individual material directions L, R and T were governed by where D i , f i and p i are the plastic strains for the given loadings, fracture strains and plastic strains for i = L, R and T , respectively. The crack initiation (linked to the elimination of the first element) was driven by the condition while the gradual deletion of the surrounding elements then corresponded to the crack propagation.
Fracture strains for individual directions depended on the stress triaxiality, which had to be determined first as where m is the mean stress and is the equivalent stress according to the von Mises yield criterion. Stress triaxiality ranges from −∞ to ∞ and secured tensile-compressive failure asymmetry. The fracture strains taken from Šebek et al. (2020, 2021) were slightly recalibrated by the trial-and-error method so that the average cutting forces corresponded to the experiments as much as possible. Around hundred runs were needed for the final dependency given in Fig. 8 and analytically described for the studied European beech as It should be noted that = −1∕3 corresponds to the uniaxial compression and = 1∕3 to the uniaxial tension in theory and that the given fracture strains correspond to the element size (given further) that was fixed to minimize the mesh dependency.
The knife was modelled as rigid (because the filtration was applied to the experimental results in order to eliminate the vibration of the knife) using a surface that came into contact with the sample. The very fine mesh was needed in the chip surroundings as discussed later. Therefore, the width of the specimen was modelled as 0.02 mm instead of 20 mm in order to save computational time. No effects were expected to occur along the width in the numerical simulations since nothing of the kind was observed during the investigation either. The geometries of the knife and sample are given in Fig. 9. Dimensions are provided by Dvoracek et al. (2021b).
The knife was discretized with 4-node bilinear quadrilateral R3D4 finite elements of the size of 0.018 mm in the process zone, 0.7 mm away from the knife tip towards the chip breaker (Fig. 9). Only one layer of elements was placed along the width of 0.02 mm (as in the case of the specimen). The element size was 0.5 mm in the remaining parts of the knife. The process zone at a distance of twice the depth of cut away from the cut edge ( Fig. 9) was discretized with 0.02-mm-sized C3D8R 8-node linear brick finite elements with reduced integration and viscous hourglass control. Behind that zone, there was an area of 0.1-mm-sized elements 1.2 mm away from the cut edge (Fig. 9). The remaining volume was filled with a coarse mesh of 2-mm-sized elements.
Rigid body motion was prevented in the case of the knife-all degrees of freedom were disabled. The initial angular velocity was prescribed to the specimen. The axis of rotation was 2025 mm right to the tip of the knife (Fig. 9). The coupling was then applied to the right surface of the specimen (Fig. 9). The degrees of freedom of the respective nodes of that surface were coupled to the axis of rotation to which the required angular velocity (the same as the initial one and corresponding to the respective cutting speed) was entered. Finally, the displacements were constrained in the width direction in order to simulate the plane strain state, as only a thin layer of elements was modelled. Otherwise, the computational cost would have been too high. The plane stress state was present only on free surfaces and its close vicinity, so the results were not biased much, as a large part of the 20 mm width was under plane strain condition.
The friction coefficient of 0.45 was considered for the whole model based on previous works by Šebek et al. (2020, 2021). It should be noted that the derived fracture strains also depend on that value that is influenced by many factors. Nevertheless, the current work does not focus on the effect of friction, since only one moisture content was tested under one (room) temperature.

Comparison of experiments and computations without rate dependency
Comparison of the average cutting forces from experiments and numerical simulations (both obtained as described in Subsection Measurements without rate dependency) is given in Table 3. The deviation with respect to the experiments is also included.   The shapes of the chips predicted by the computations are compared to the experiments (high-speed camera recordings) in Figs. 12, 13 for LR and RT wood grain directions, 90 m·s -1 cutting speed and various depths of cut. The chips were split into many small pieces in both experiments and numerical simulations.

Measurements with rate dependency
The LR wood grain direction with 0.1 mm depth of cut was chosen in order to investigate the influence of the cutting speed, while four speeds of 20 m·s -1 , 40 m·s -1 , 60 m·s -1 and 80 m·s -1 were analysed. Each of them was repeated 8 times, so there were 32 runs. The scatter of the results is similar to the one in Subsection Measurements without rate dependency (Fig. 5). Moreover, the scatter decreased slightly with decreasing cutting speed. The same procedure as in Subsection Measurements without rate dependency served for the determination of average force responses, which then served for the estimation of the average cutting forces discussed later. It is obvious that the force responses plotted in Fig. 14 decreased with decreasing cutting speed.

Rate-dependent model
The model presented in the Subsection Rate-independent model was rate independent and not capable of describing the behaviour under various strain rates-cutting speeds. Therefore, the dependency on the plastic strain rate was introduced into the fracture criterion. Other parts of the model remain the same as in the Subsection Rate-independent model. The original fracture strains presented in the Subsection Rate-independent model (  (Fig. 16). The rate-independent numerical simulations with cutting speeds of 90 m·s -1 and 20 m·s -1 in LR wood grain direction served as an initial guess. The average strain rates were extracted from the loaded zones in the mentioned simulations. Based on these results, the material parameters A and B were set so that the fracture strains were similar to those in Subsection Rateindependent model for 90 m·s -1 and lower for 20 m·s -1 . To verify the validity of the model, one case from the Section Influence of wood grain direction and uncut chip thickness on cutting force (LR wood grain direction, 90 m·s -1 cutting speed and 0.1 depth of cut) was calculated to find that the average cutting force was similar to the model without rate dependency (only more pronounced oscillations of the force occurred) (compare Tables 3, 4). It should be noted that the model did not account for any damping.

Comparison of experiments and computations with rate dependency
Average force responses from experiments and numerical simulations are compared in Fig. 16 for LR wood grain direction, 0.1 mm depth of cut and various cutting speeds. The computations corresponded well to the experiments, which is documented by the deviations in Table 4.
The shapes of the chips predicted by the computations are compared to the experiments (high-speed camera recordings) in Fig. 17 for LR wood grain direction, 0.1 mm depth of cut and various cutting speeds. The lower the cutting speed, the more the wood tends to form the spiral chip. There are many small pieces of chips for 60 m·s -1 cutting speed, semi-spiral chips for 40 m·s -1 cutting speed and a spiral one for 20 m·s -1 cutting speed (Fig. 17). The model did not predict spiral chip formation, which still resulted in fine pieces of chip smaller than those predicted by the model without strain rate dependency-compare with Figs. 12, 13. However, a Fig. 15 Rate dependency of fracture strain in L direction different set of material parameters allowed for the modelling of spiral chip formation (Fig. 18). Nevertheless, it was still not exactly as in the experiment, so these parameters are not given here. Moreover, it led to a much worse overall prediction of cutting forces, which were prioritized over the chip shapes in the present article.
The rate-dependent material model produced much more pronounced force oscillations, which could be caused by many factors. The first knife-specimen contact was an impact, causing an abrupt increase in plastic strain rate and the propagation of stress waves, among others. Then, the plastic strain rate was highly inhomogeneous in the process zone. The knife-specimen contact was interrupted when the chip split into small pieces, causing shocks that appeared in the form of  the force oscillations. This was confirmed by the last computation, where there was a quite continuous-spiral-chip. Then, the force oscillations were negligible, even smaller than those in the respective experiment.

Conclusion
Linear hardwood disintegration was examined at cutting speeds of up to 90 m·s -1 . Specimens prepared from European beech rotated with a rotor arm approximately 2 m in length. The cutting process was recorded by a piezoelectric force sensor and a high-speed camera in order to study the effect of process parameters (wood grain direction, depth of cut and cutting speed). These examinations served as a basis for the development of a finite element model. It was found that cutting perpendicular to the grain growth (RT wood grain direction) resulted in approximately twice the forces of cutting parallel to the fibre direction (LR wood grain direction). Naturally, the larger the depth of cut and the cutting speed, the greater the force responses. A finite element model, rate independent and dependent, was developed and implemented for force prediction. It resulted in good agreement with the experiments, except for the smallest and the greatest depths of cut, the LR wood grain direction and the cutting speed of 90 m·s -1 . In addition to the force responses, the formation of the chips was followed. It was captured satisfactorily by the computational model, which predicted variously sized small pieces of the chip for all cases, which was against the experiment only for the lowest cutting speed (20 m·s -1 ). Nevertheless, it was overcome by recalibrating the model, which also additionally revealed that the parasitic force oscillations did not appear in numerical simulations with the spiral chip. It should be noted that the recalibrated model only showcased the prediction capabilities, and it resulted in a much higher force response than in the respective experiment. Finally, the model built in Abaqus/Explicit may be used for an optimization of the cutting process or even its further investigation, like the influence of the tool wear. Fig. 18 Shape of chip from experiment (left) and computation (right, different set of material parameters not given in this article) for LR wood grain direction, 0.1 mm depth of cut and cutting speed of 20 m·s -1 (time of 2.5 ms)