Fracture predictions in impact three-point bending test of European beech

Hardwood has become widespread in European forests. The strongest factor is climate change and damage to coni-fers by the bark beetle. The effort to study hardwoods grows with increasing volume of applications. Therefore, European beech wood was investigated under two impact loads in two material directions, resulting in four unique combinations supplemented by the measurement of the friction coefficient. Then, it was computationally simulated to reproduce the cracking, while the material model reflected the orthotropic behaviour in elasticity, plasticity and failure. The model was coded using the user subroutine in Abaqus to initiate and propagate the crack using the element deletion. The resulting reaction forces were in good agreement with those from the experiments. Cracking was numerically simulated in three of four cases as experimentally observed, however, upon larger deflections. Therefore, the model is applicable for further investigations.


Introduction
Wood is a widespread natural composite, which has been used as a construction material for centuries because it is renewable, biologically degradable, environmentally friendly, strong, lightweight, easy to manufacture, electrically resistant, absorbing noise and aesthetic.Therefore, models capable of predicting deformation and failure under operational or random loading are needed.
Bending is one of the common tests of wood.Jansson [1] investigated failure modes and stresses using static and impact bending tests showing decreasing impact bending strength with decreasing time to failure for the Spruce Pine Fir (SPF).Yoshihara et al. [2] used a static three-point bending test to determine the shear modulus with the help of a correction function and a modified Timoshenko beam for six wood species: Sitka spruce, western hemlock, akamatsu, yellow poplar, shioji and balsa.Yoshihara et al. [3] conducted a similar study for asymmetric four-point bending tests.Yoshihara and Oka [4] carried out the compression bending test to determine the elastic modulus, proportional limit and bending strength on specimens of Japanese fir of various lengthto-thickness ratios.Then they compared the results with the conventional bending test to confirm that the correct flexural properties can be obtained from the compression bending test for a large length-to-thickness ratio.Kubojima et al. [5] used the impact bending test to estimate the elastic modulus for Japanese cedar, hondo spruce, hiba arborvitae, Japanese red pine, paulownia, Manchurian ash and Japanese evergreen oak.Polocoşer et al. [6] investigated the effect of low-velocity impact on failure stresses and stiffness using the bending test consisting of a pendulum and specimens having a length of 650 mm, width of 50 mm and thicknesses of 20, 30 and 40 mm to find that failure was significantly different compared to quasi-static tests for beech, larch and pine.Polocoşer et al. [7] found higher energy absorption in pine and spruce reinforced with E-glass on the surface loaded in tension on a pendulum compared to specimens without reinforcement.The factor was 1.4 for pine and 2.5 for spruce.Jacques et al. [8] carried out full-scale testing on individual light frame lumbers made from SPF.Then, the four-point bending tests with strain rates of 6 × 10 -6 to 4 × 10 -1 s -1 served to propose the stress-strain relationship applicable to the modelling of blast loadings.Olmedo et al. [9] followed the dynamic response of fresh stems under impact loading using the Mouton-Charpy pendulum.Then they analysed the load-bearing capacity of wooden constructions made of felled trees serving as protection against falling rocks.Brancheriau et al. [10] proposed an analytical relationship, which converts the elastic modulus estimated from the three-point bending test to the elastic modulus appearing in the fourpoint bending test.The results were experimentally validated on specimens from six wood species with different densities.Different elastic moduli of spruce and oak in three-and four-point bending tests were also analysed by Babiak et al. [11].Gaff et al. [12] examined the effect of thermal modification of European oak and Norway spruce using the Charpy pendulum to find that changes in basic chemical components and impact bending strength were less affected for oak than for spruce.The effect of thermal modification on impact bending strength was also investigated by Hassan Vand and Tippner [13] to find that thermal modification caused a decrease in deflection and maximum longitudinal strain up to approximately 50% according to digital image correlation (DIC) for three-point bending tests of five wood species: ash, beech, larch, oak and spruce.Hassan Vand et al. [14] also tested the effect of moisture content on the behaviour of three wood species of beech, oak and spruce.They evaluated the work required to initiate the crack and break the specimens and found, using DIC, that the maximum deflection and longitudinal tensile strain increased with increasing moisture content.
The minimum of the literature deals with the computational modelling of wood under impact.Therefore, this work focuses on the dynamic behaviour of European beech (Fagus sylvatica L.) under the three-point bending test.All tests were carried out using a drop-weight impact testing machine with two hammers to achieve various impact velocities.Then, numerical simulations followed to develop a material model with good fracture predictability compared to experiments.It should be noted that many works [15,16] use the orthotropic yield criterion according to Hill [17] for the description of plastic behaviour.However, some material parameters (three out of six) for the equivalent stress can be negative due to significantly distinct yield stresses in respective wood directions [18][19][20].The negative material parameters can then cause the equivalent stress to be undefined.Therefore, the orthotropic yield criterion according to Barlat et al. [21] will be used to overcome this issue.

Material
The material studied was European beech belonging to hardwoods.The specimens were made from tree grown near Brno (Czech Republic).The wood was considered orthotropic with a direction parallel with three primary (height) growth called longitudinal (L), a direction perpendicular to the tree primary growth passing through the pith of the tree called radial (R) and a direction perpendicular to tree growth and simultaneously as a tangent to annual rings called tangential (T) as depicted in Fig. 1.All specimens were stored in a climatic chamber at a temperature of 20 °C and a relative humidity of 65% to produce a uniform equilibrium moisture content of 12%.

Three-point impact bending test
All specimens were manufactured with fibres oriented along the length of 300 mm having a square cross section with an edge of 20 mm.Before testing, specimens were weighed in analytical balance with 1 mg readability, which resulted in an arithmetic mean of density of 727 kg/m 3 with a standard deviation of 33 kg/m 3 .The symmetric three-point bending was carried out in two material directions R and T (Fig. 1), respectively.The span-to-depth ratio was 12, resulting in a support span of 240 mm.Both the supports and the hammer had a radius of 15 mm and a surface roughness Ra of 0.2 μm.
The tests were performed on the drop-weight impact testing machine DPFest 400 from Labortech (right Fig. 2) at room temperature.Two different hammers were dropped from two initial heights to obtain two impact velocities as given in Table 1 to provide a diverse material for numerical simulations.The testing machine automatically measured the true impact velocities (Table 1).For each of these two conditions, 6 specimens were tested in the R and T directions, resulting in 4 datasets.The hammer displacement was measured with 0.01 mm precision, while the force was measured using a piezoelectric force transducer CFT + 50 kN from HBM attached to the hammer.
The tests were recorded with Fastcam SA-X2 from Photron.The high-speed camera had a cell size of 20 μm and was equipped with a Nikon Micro-Nikkor G lens with a focal length of 105 mm and a Nikon Z TC-2.0 × .It was placed approximately 0.9 m from the lateral specimen surface, which was parallel to the camera sensor.Additional lighting was ensured by two MultiLed QT standalone lamps (left Fig. 2).The field of view was fitted to a central part of lateral specimen's surface with a hammer.Images with a resolution of 1024 × 672 px were captured with a frame rate of 20,000 fps.The image postprocessing was performed in Vic-2D 2010 software from Correlated Solutions.In order to determine the conversion factor, a simple calibration was done using the one-dimensional scale defined by the known distance.
The displacement field was calculated from the field of 3 × 3 pt.The Lagrange strain field was determined with the lowest possible strain filter size of 5 × 5 pt so that the maximum possible spatial resolution was achieved.In

Measurement of friction coefficient
Tribological experiments were performed to collect input for numerical modelling.The tests were carried out using the Universal Mechanical Tester (UMT) TriboLab from Bruker (Fig. 3).The pin of 5 mm in diameter was made of AISI 52100 alloy steel having a surface roughness Ra of 0.2 μm, the same as the supports and the hammer in the three-point impact bending test.It was forced against wooden specimens by 196 N, which was measured by a dual force sensor DFH-100 with a sampling frequency of 100 Hz and corresponded to a contact pressure of 10 MPa.This contact pressure was chosen as it reached tens of MPa in the following numerical simulations.The pin oscillated with a frequency of 2 Hz and amplitude of 5 mm for 120 s.Six specimens with dimensions of 40 × 20 × 10 mm were tested.These specimens were manufactured from the same tree as those for the three-point impact bending test and stored in the climatic chamber to reach the same uniform equilibrium moisture content of 12%.The pin was forced against three LR and three LT surfaces of 40 × 20 mm moving in the R and L directions, respectively.Therefore, the measurement was within one annual ring for the L direction and approximately 5 annual rings for the R direction with an average annual ring width of 2.5 mm.
The pin trajectory was 10 mm, but the results were processed in MATLAB R2024a for a 6-mm-long portion of the specimen in order to omit the dead ends where the pin decelerated and accelerated.Then, the considered friction was kinetic and area for evaluation approximately 50 mm 2 .Also, the evaluation was not conducted for the beginnings of the tests where a ploughing was present.

Numerical simulations
The calculations were performed in Abaqus/Explicit commercial code based on the explicit formulation of the finite element method.The material model was implemented using the VUMAT user subroutine, as it is not a standard one.The failure was modelled by deleting elements that reached critical damage, which is a simple technique that does not require remeshing, contrary to the node separation method.

Model of material
The hammer and supports were made of steel, which is much stiffer than the tested wood specimens.Therefore, the steel parts were modelled as rigid, saving some computational time.On the contrary, wood was modelled as a homogeneous orthotropic continuum.The heterogeneity was neglected to save some computational time again through material model simplification, including its calibration.
The orthotropic elasticity was described by the generalised Hooke's law as where ε i , E i and σ i are the normal strain, elastic modulus and normal stress, respectively, ε ij , ν ij , G ij and σ ij are the tensorial shear strain, Poisson's ratio, shear modulus and shear stress for i, j = L, R, T , respectively.All elastic con- stants (Table 2) were taken from [20], where the disintegration of the same wood was carried out under various strain rates in various material directions.
The orthotropic plasticity was described by the model of Barlat et al. [21], who proposed the yield condition as (1)  where σ y is the yield stress, which is dependent on the equivalent plastic strain (Fig. 4), while σ B is the equiva- lent stress according to Barlat et al. [21] as (3) where m is the exponent (influencing the shape of the yield surface), K 1 , K 2 and K 3 are the principal values of the linearly transformed deviatoric stress tensor where a , b , c , f , g , and h are the material parameters, which were fitted (Table 3; Fig. 5) so that the yield surface corresponded to that in [20], where the model of Hill [17] was used with a negative material parameter (H < 0) .A change in the yield criterion was sought due to the elimination of a negative material parameter resulting from a significant difference in the yield stresses in the L direction compared to others (R and T in Fig. 5), which can be problematic in some cases.The yield criterion was taken from [20], where it was not necessary to model different yield stresses in tension and compression in finite elements.Then, the yield criterion was recalibrated for European beech, still neglecting the strength differential effect.However, more sophisticated yield criteria should be considered when significant differences in tensile and compressive yield stresses were observed.
(  The orthotropic cumulative damage was described by the damage parameter as: where ε f i , ε D i and ε p i are the fracture strain, plastic strain for a given loading path and plastic strain, respectively.As mentioned above, the element is removed when the critical value is reached, specifically when max (D i ) = 1 .Furthermore, fracture strains are dependent on stress triaxiality (apart from the rate dependence introduced later) to incorporate the tension/compression failure asymmetry (which is simpler than the fracture models in [15,[22][23][24]) where σ m is the mean stress and σ is the equivalent stress according to von Mises The dependence of fracture strains on stress triaxiality was based on [20] and was further refined using the trial and error method using tens of numerical simulations so that the results matched experimental observation in impact three-point bending tests.Fracture strains dependent on the stress triaxiality, ε f i (η) , are shown in Fig. 6 for the reference equivalent plastic strain rate of 1 s -1 .
Finally, fracture strains were additionally dependent on the equivalent plastic strain rate as experiments were conducted at various impact velocities.Contrary to the linear rate dependence in [20] for significantly greater strain rates, an exponential dependence on the equivalent plastic strain rate based on the equation proposed by Johnson and Cook [25] was recalibrated to fit the experiments as follows (using ε f i (η) that has already been calibrated previously for the reference equivalent plastic strain rate of 1 s -1 ) where C is the material parameter, calibrated as 0.1 after several additional numerical simulations for reference (5) equivalent plastic strain rate ε0 = 1 s −1 .Finally, εp is the equivalent plastic strain rate where εp is the plastic strain rate tensor.

Model of geometry and boundary conditions
As mentioned above, the advantage of the element deletion technique is its simple implementation.However, it has its drawbacks as it depends on the size of the element, which is usually kept as small as possible to realistically propagate the crack using one or two elements.Therefore, the elements had a size of 0.1 mm in areas of wood failure as well as in contact regions (Fig. 7).Other areas were meshed with elements of a size of 1.82 mm.Only a 0.1-mm-width specimen was modelled as a plane strain to represent the inner layer of wood instead of using 20 mm to save considerable computational time.Geometry was meshed with 8-node linear brick finite elements with reduced integration and hourglass control (labelled C3D8R in Abaqus).The mesh consisted of 128,774 nodes and 63,850 elements of which 4 were additionally 6-node linear wedge finite elements with reduced integration and hourglass control (labelled C3D6R in Abaqus) near the contact regions where the fine mesh changed into coarse mesh (highlighted in red in Fig. 7).
As mentioned earlier, the hammer and supports were modelled as rigid, therefore, as surfaces which come into contact with the specimen.The width of the rigid  6 Fracture strain for all material directions and the reference equivalent plastic strain rate of 1 s -1 bodies, 0.3 mm, was greater than the width of the specimen, 0.1 mm, and were centred with each other.The hammer and supports were meshed with 4-node bilinear quadrilateral rigid finite elements (labelled R3D4 in Abaqus) having a size of 0.1 mm, therefore, 3 elements per width.One support, modelled the same as a hammer, had 1892 nodes and 1416 elements.Moreover, an element was added to the reference point of the hammer (highlighted by a blue pentagram in Fig. 7) in order to incorporate the volume load (labelled MASS in Abaqus), equivalent to the hammer weight.
Displacements and rotations of reference points of the supports were restricted.The same applies to the hammer except for the vertical displacement.The initial velocity was prescribed in this direction, corresponding to the average true impact velocity of the experiment (Table 1).The specimen had constrained displacements only in the width direction to simulate the plane strain condition (the geometry of the specimen was discretised by 1 element along the width).The acceleration due to gravity of 9.807 m/s 2 was applied on the whole geometry in the direction of impact.
In addition to the definition of the contact between the specimen, hammer and supports, self-contact was applied to all element edges in the area of expected failure (where the elements with 0.1 mm edges were).It ensured the contact of new free surfaces, which were not initially present in the geometry, emerging after the element deletion that simulated the cracking.The normal behaviour of the contacts was set as 'hard' to allow any pressure when the surfaces are in contact, while the friction coefficient of 0.15 obtained from the experiments in the previous section was set in the tangential direction.

Results and discussion
Figures 8 and 9 summarise all the force responses against the hammer displacements obtained from the impact testing machine for hammers that weighed 9.05 and 4.55 kg, respectively.The DIC did not serve to obtain the displacements.The dotted lines highlight the tests in which the specimens remained intact.The corresponding specimens bended forward and back (the hammer returned approximately to the reference position corresponding to the displacement of 0 mm-it is zero deflection).The dashed-dotted lines highlight the tests in which the specimens partially broke (did not break through the whole height).The corresponding specimens exhibited significant springback, therefore, the hammer returned approximately to the reference position.The solid lines highlight the tests in which the specimens failed (broke through the entire height).The corresponding specimens remained bent after the test (all in Fig. 8).Moreover, the moment of crack initiation is highlighted by a circle.The average value of the hammer displacement corresponding to the moment of crack initiation is highlighted by a thick vertical dashed grey line.This line is not plotted in left Fig. 9, where only 2 out of 6 specimens partially broke and none failed for the hammer weighing 4.55 kg loading in the R direction.Only 1 specimen failed for the hammer that weighed 4.55 kg in the T direction (right Fig. 9).Other specimens partially broke, while the specimen corresponding to Test 6 in right Fig. 9 was cracked along almost the whole height under these conditions, therefore, the hammer started returning later.
The arithmetic mean friction coefficient of 0.15 with a standard deviation of 0.005 was measured irrespective of the material direction.Then, the experimentally and computationally obtained results are compared in Figs. 10 and 11.The predicted forces were in good Fig. 7 Finite element mesh with details agreement with the experiments before the maximum deflection.Forces oscillated more in numerical simulations than in impact three-point bending tests.However, this is an intrinsic feature of explicit dynamics.Cracking was predicted later in numerical simulation than in experiments for a 9.05 kg hammer in Fig. 10  was achieved for the R direction with a 4.55 kg hammer (left Fig. 9) computationally to reproduce the experiments (left Fig. 11).However, partial cracking was not achieved computationally for the T direction with a 4.55 kg hammer as in the experiments (right Fig. 9).The elastic strain energy was so high in numerical simulations for that case that led to the prediction of crack propagation along the whole height, contrary to the impact three-point bending tests where specimens broke just partially (right Fig. 9).The homogeneity of the model also contributed to the lack of prediction of crack arrest, which occurred in experiments after some delamination because, among others, the wood is strongly inhomogeneous.This should be taken into account by modelling the geometry based on data from computed tomography specifically for each specimen, which might not be practical for industrial applications.
The predicted contours of the damage parameter are displayed in Figs. 12 and 13 after the impact three-point bending tests compared to experiments approximately in scale.However, it should be noted that the predicted deflections are greater (approximately 1.5 mm), as is evident from Figs. directions, respectively.The predicted cracks were vertical for the T direction regardless of the weight of the hammer.The crack bifurcated for 9.05 kg hammer and then unified and propagated in a single path after a while (left Fig. 12), whereas no crack initiated for the R direction and 4.55 kg hammer (left Fig. 13).Moreover, right Fig. 13 shows the test, where the highlighted significant zigzag cracking stopped after a while, which was not reproduced by the numerical simulation as discussed above.It could stop because, among others, the strain energy transformed into kinetic energy of large pieces chipped off, which lowered the crack driving force.Finally, predictability was assessed based on the DIC in Fig. 14, where the total horizontal Lagrange strain is compared to the numerical simulation, which provided a true (logarithmic) measure.It was evaluated for a 4.55 kg hammer loading the specimen in the T direction with a displacement of 4 mm.It should be noted that the DIC was evaluated on the face (free surface under plane stress condition), but the computation corresponds to the condition inside the specimen (plane strain).Regardless, the match is quite good, with corresponding two localised regions under the hammer that resulted from contact (the two strain measures are close for small strains).In general, the contours from DIC are not as smooth as from computations due to actual wood heterogeneity.The extreme tensile strains were in the bottom part of the bent specimen in the L direction, as expected, while the top part was loaded in compression along the L direction as well as in the impact direction (R or T, respectively) in the vicinity of contact with the hammer.There, both the compressive strain in the L direction and equivalent plastic strain were approximately twice greater compared to the bottom part (Fig. 14).Still, the failure initiated from the bottom part of the specimen in the finite element simulations due to the correctly calibrated fracture model, which accounted for the tension/compression failure asymmetry apart from orthotropic failure behaviour.

Conclusions
The European beech was studied under impact threepoint bending.Two impact velocities and two material directions were experimentally tested.The specimens completely failed under a 9.05 kg hammer in both directions, partially broke under 4.55 kg in the T direction, and almost not cracked under 4.55 kg in the R direction.All tests were modelled computationally, for which the measurement of friction coefficient was carried out, considering orthotropic behaviour of wood in elasticity, plasticity and fracture.Then, the reaction forces were in good correspondence with the experimental observation.However, the cracking was generally predicted later than in reality (under greater deflection).The geometry was modelled as a homogeneous continuum, which caused unsatisfactory reproduction of crack propagation, which was realised through the element deletion technique within the explicit finite element method.This could potentially be solved by techniques that have roots in classical fracture mechanics rather than using phenomenological or empirical criteria.

Fig. 1
Fig. 1 Loading in the R (top) and T directions (bottom) with all material directions during impact three-point bending tests (all dimensions in mm)

Fig. 2
Fig. 2 High-speed camera with accessories (left) and drop-weight impact testing machine (right)

Fig. 3
Fig. 3 UMT TriboLab used for measurement of friction coefficient

Fig. 4 Fig. 5
Fig.4 Flow curve of the European beech wood[20] Fig.6 Fracture strain for all material directions and the reference equivalent plastic strain rate of 1 s -1

Fig. 9 Fig. 10
Fig. 9 Force-displacement responses for a 4.55 kg hammer in the R (left) and T directions (right), respectively (for interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article)

Fig. 11 Fig. 12
Fig. 11 Comparison of force-displacement responses from numerical simulations and experiments for a 4.55 kg hammer in the R (left) and T (right) directions, respectively

Fig. 13 Fig. 14
Fig. 13 Experiments (top with highlighted zigzag cracking) compared to predicted contours of the damage parameter (bottom) for the 4.55 kg hammer and the R (left) and T (right) directions (approximately in scale)

Table 1
Configuration of the impact three-point bending tests

Table 2
[20]tic constants for European beech wood[20] . No cracking Force-displacement responses for a 9.05 kg hammer in the R (left) and T directions (right), respectively (for interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article)