Elasto-plastic material model of green beech wood

Physically modelling the mechanical response of a tree by numerical simulation depends on having accurate data on the mechanical properties of green hardwood. Lacking such data, we developed and validated an orthotropic elasto-plastic (E–P) material model, based on the results of experiments performed on European beech ( Fagus sylvatica L . ) green wood, capable of including both the non-linearity and orthotropic properties of the material. We selected 655 clear samples with the special orthotropic structure of annual rings. All samples were prepared immediately after felling; their moisture content (MC) was 80% on average. The mechanical responses in normal directions and shear are represented by bi-linear stress–strain curves. The E–P model was validated by comparing the force– deflection response of three-point bending of green wood samples in a finite-element method (FEM) simulation (the average relative error was 4.6% for point-wise and 1.7% for integral-wise comparison). The output of this work was a consistent set of material constants for the E–P material model that is now available for the structural analysis of beech wood with MC above to fibre saturation point (FSP), especially green wood, subjected to relatively high loads (such that a plastic deformation appears) and that can very well predict a non-linear response above the proportional limits.


Introduction
The mechanical properties of wood are one of the key parameters in the mechanical analysis of natural wood structures such as trees [1].The moisture content (MC) of wood in standing trees exceeds the fibre saturation point (FSP), affecting its mechanical properties compared to construction wood [2].The mechanical properties of wood do not change dramatically above the FSP limit; however, for consistency, green wood properties should be measured above 50% MC [3].Various device-supported methods, like the pulling test [4], utilise green wood's mechanical properties to assess tree stability.These include assessing force and longitudinal strain using mechanical sensors [5] or optical techniques like digital image correlation (DIC) [6].Results of this tests are afterwards compared to catalogued green wood properties [7,8], detecting changes due to defects or root damages.
In recent years, numerical models based on the finiteelement method (FEM) have been established as a reliable instrument to predict the mechanical behaviour of trees.Jackson et al. [9] combine terrestrial laser scanning data and FEM to predict strain on the trunks of different trees and validate the model using strain data measured on the scanned trees.Vojackova et al. [10] present an FEM simulation of tree response to static loading and the interaction between stem and root-plate stiffness, which they validate by comparison with field-measured data.Middleton et al. [11] use different mechanical finiteelement (FE) models of tree forks to replicate rotations in the joints and validate them by bending mechanical testing of different Salix alba inosculations.
Material parameters for the structural analysis of trees are usually taken from mechanical tests of compression/ tension parallel to the grain or from bending tests [12,13].Wessolly and Erb [8] describe an approach to assessing tree stability by the pulling test method, where the longitudinal elastic modulus (E) in compression and the yield point (σ) in compression parallel to the grain are used to characterise the material.Lundström et al. [12,14] provide the material properties of living stems to predict the E and present an idealised stress-strain diagram.Niklas and Spatz [3] analyse the relationships among the elastic parameters of green wood from a taxonomically broad range of 161 species.Nevertheless, the biological nature of wood and its anisotropic structure make wood a material arduous to model [15,16].The mechanical behaviour of wood shows non-linear deformation in compression, mainly due to the collapsing of cells and the separation of fibres.In tensile loading and shear, the stress-strain behaviour is more linear with brittle failure [17].In simulating the mechanical response of woodbased composites through the generalised Hill's potential theory, Moses and Prion [17] introduce an artificial requirement: the equality of tensile σ and compressive σ to maintain model consistency.Furthermore, the σ is set artificially high and, however, this approach assumes a symmetrical response of wood to mechanical load which is not generally found, especially with high loads.Plastic behaviour above the σ is important for a reliable description of strain development in wooden structures [18] including trees.
The development of numerical methods allows more realistic models of wood and wood-based materials and structures to be used [19,20].Mackerle [21] published a comprehensive review of the general use of numerical models in wood analysis.The first general model of anisotropic plasticity was introduced by Hill [22].For the E-P modelling of mechanical response, 27 bi-linear material parameters and three Poisson's ratios are needed [19,23,24].The response of models is most often compared with the response in an experimental bending test, as this test combines tension, compression and shear, and bending is a very common mode of real loading in practice [25].Clouston and Lam [26] use a continuum mechanics approach to simulate the mechanical properties of wood strand composites by the implementation of an orthotropic E-P failure model.Moses and Prion [27] apply plasticity and the Weibull weakest link failure model to predict the properties of wood and wood-based composites to reach non-linear plastic behaviour in each of the orthogonal material directions.Yoshihara [28] proposes equations based on the beam theory to obtain critical stress intensity factors and assess the influence of the elastic properties of a double cantilever beam and three-point end-notched flexure.Wood compressive behaviour using the bi-linear material in fibre-reinforced glulam beams is modelled by Fiorelli and Dias [29].Oudjene and Khelifa [18] present a constitutive three-dimensional (3D) anisotropic plastic model for dried wood and use 3D FEM to describe the non-linear compressive pattern of wood.Hering et al. [30] show the moisturedependent stress-strain relationship of European beech wood under compression and describe the significant influence of MC on stress-strain behaviour for elastic and plastic domains.Hong et al. [31] describe a simplified non-linear orthotropic model of wood by reducing the number of parameters from 27 (for the general orthotropic bi-linear model) to six, with the experimental validation in compression [19], where E t is defined as a fixed small fraction (0.001) of the E. Using an E-P model, Milch et al. [23] describe the mechanical properties of dried wood by experimental and numerical methods by doing compression tests parallel and perpendicular to the grain in different directions, a three-point bending test and double shear joints with FEM.Šebek et al. [32] use an orthotropic non-linear model of beech wood based on split Hopkinson pressure and tensile bar experiments to reflect the anisotropic behaviour for tensile/compressive asymmetry.Tippner et al. [24] describe an orthotropic E-P model for oak wood at two different MC levels applicable for FEM.Although the E-P behaviour of wood has been investigated in the past, bi-linear numerical models of green wood are still lacking.The objective of this study is to determine the E-P material properties of a representative hardwood-the European beech (Fagus sylvatica L.), which is one of the most commonly occurring species in European forestry and arboriculture.In addition, a new numerical method is presented to identify the elastic and plastic zones of strain-stress curves.The sets of material constants are used in an FE model which is compared with the results of experimental tests in bending.This process produces a consistent set of material constants for an E-P material model that is available for the structural analysis of green beech wood where the wood is subjected to large deformations and plastic deformation takes place.

Materials and methods
A bi-linear orthotropic material model was used to simulate the E-P behaviour of green beech wood.The constitutive relation of the material loading and material response used a continuous curve consisting of two linear parts, whose slopes are referred to as E and the tangential modulus (E t ); the point where the two linear parts meet is referred to as σ.Thus, three material constants (E, E t and σ; G, G t and σ for shear, respectively) were used for each of the nine types of loading (normal tension and normal compression in three directions: longitudinal-L, radial-R and tangential-T; and shear in three planes: longitudinal-radial-LR, longitudinal-tangential-LT and radial-tangential-RT), giving a total of 27 constants for an orthotropic material.

Anisotropic material model in ANSYS mechanical APDL
The ANSYS ® Mechanical APDL FEM solver, release 19.2 (ANSYS Inc., USA) was used for the computation as a robust and widely accepted FEM tool.The ANSYS ANISO material model can handle both the anisotropy and the non-symmetric non-linearity of the material.However, this material model is based on the generalised Hill's potential theory in ANSYS ® and suffers from the strict requirements of this theory.More precisely, the E in tension and compression are supposed to be equal and the σ in tension and compression are required to satisfy the consistency equation Eq. (1.1) and the condition Eq. (1.2) on the M-tensor (see below) through the whole duration of loading.
The consistency condition (see ANSYS-help [33]) requires where σ +i and σ −i are yields in tension and compression (where i is one of L, R or T), respectively.
The condition on the M-tensor relates diagonal com- ponents of the tensor M = M ij defined by the relations: for a constant K and reads as (see ANSYS-help [33] for more details).

Hill's theory for anisotropic material with one dominant direction
Conditions (1.1) and (1.2) present a kind of artificial restrictions which must be met during the FEM simulation.A closer investigation reveals that even though the conditions from Hill's theory are symmetric, the condition (1.1) is highly sensitive to the properties in R and T directions and (1.2) is highly fragile.We briefly describe (1.1) the restrictions coming from these conditions and how to overcome them.Equation (1.1) can be rewritten in the form Thus, the harmonic mean of yield moduli in tension equals that in compression.However, the harmonic mean is less intuitive than arithmetic or geometric means.Large changes in dominant values have minimal impact, while small changes in lesser values are substantial.This reflects behaviour of solid wood, where mechanical properties in L directions are dominant to those in the R and T directions.Condition (1.1) is sensitive to the mechanical properties in R and T direction and less sensitive to L direction.
The left-hand side of inequality Eq. (1.2) is a bi-linear form with matrix: This matrix describes a cone in a three-dimensional space and (1.2) claims that the point (M 11 , M 22 , M 33 ) is inside this cone.The material properties of wood imply that this point is always near the boundary of the cone.Really, the anisotropy of wood implies that the yields product σ −L σ +L is much higher than the yield products σ +R σ −R and σ +T σ −T .Thus, M 11 is small compared with M 22 and M 33 .Consequently, the point (M 11 , M 22 , M 33 ) is close to the plane X = 0 which is a tangent plane to the cone defined by matrix A .This reveals that the M-tensor condition Eq. (1.2) is not robust against perturbations of the quantities M ii , i ∈ {1, 2, 3} , since a small change in these quantities may shift the point outside the cone prescribed by the condition (1.2).
To overcome the difficulties from Hill's theory, σ values which satisfy both Eqs.(1.1) and (1.2) were used and the E t was redefined to be close to zero (following Hong et al. [19]).This ensured that Eqs.(1.1) and (1.2) remained valid after loading and kept more realistic values for all σ (equal in compression/tension, in correspondence with Moses and Prion [17]).
The consistent and validated E-P material model was developed in three main steps.The first step was to identify the mechanical properties of green wood for each loading mode independently.Subsequently, the constants for the E-P material model were established.The output from this step was 27 constants referred to as a bi-linear material model (B-MM).In the second step, the data from all normal loading modes were processed simultaneously and adjusted to the Hill's theory and conditions Eqs. (1.1) and (1.2).The output of this step is referred to as a computational material model (C-MM).Finally, a FEM simulation of three-point bending was compared with the experimental data.This step produced information for fine-tuning the constants in the computational material model.Using these constants, we improved an agreement between the FE model and the experiment in three-point bending.The output of the last step is referred to as an adjusted computational material model (AC-MM).

Timber selection and sample crafting
Four European beech trees were chosen from a natural beech forest stand (GPS 49.2753611 N, 16.7176978 E; 485 msl) in the university enterprise Masaryk Forest in Křtiny (Brno, CZ) and felled.The trees' average height was 32 m, the breast height diameter was 45 cm and the age of the trees was 95 years.One-metre-long logs were taken from a height of approximately 10 m on each trunk to make the samples.Due to the visual assessment and the high of trees, the part of the log was chosen as a representative source of material without defect (decay, cracks, wood structure deviation) with a common range of density representing the beech wood.The timber was processed within 4 days after felling, one log per day.For mechanical testing, ten types of the specially orthotropic samples (Fig. 1) were crafted.To keep the MC as close as possible to the  2) parallel to the grain; C tension tests perpendicular to the grain; D three-point bending test; E tension tests parallel to the grain.For brevity, the loading is omitted in pictures A3, A4, B2 and C2 and the dimensions are omitted in A1, A2, B1 and C1.In these cases, the values are identical in all corresponding tests (adapted from Tippner et al. [24]) original MC in standing trees, the test samples were crafted immediately after log processing; the MC did not fall below 50% during crafting.This was achieved by storing all samples over water in sealed boxes to ensure the high relative humidity of the air within (99%).

Mechanical testing
The wood samples were tested in compression and tension parallel and perpendicular to the grain in the R and T directions; in shear at the LT and LR planes, and by three-point bending in the T direction-the notation accords with Hearmon [34], as presented in Tippner et al. [24].
The experimental setup (i.e. the dimensions of the samples and the course of the tests) was adopted from the British standard BS 373 [35] and the Czech standards for compression [36][37][38], shear [39] and tension [40,41].A three-point bending experiment described in standard ČSN 49 0115 [42] was used for model validation.The dimensions of the samples and directions of load are described in Fig. 1.
All samples were destructively tested on a Zwick Z050/ TH 3A universal testing machine (Zwick Roell AG, Ulm, Germany) equipped with a 50 kN load cell (Fig. 2).The experimental procedures were controlled using TestXpert v.11.02 (Zwick Roell AG, Ulm, Germany).The temperature was kept within a range of 20-22 °C during all mechanical tests.To avoid MC losses, all samples were processed immediately after removal from the sealed boxes.
The surface of each sample was covered with a black sprayed speckle pattern and the displacement was measured using DIC for all measurements except tension parallel to the grain.A pair of AVT Stingray Copper F-504 CCD cameras (Allied Vision, Stadtroda, Germany) with a cell size of 3.45 μm and resolution of 2452 × 2056 pixels, equipped with Pentax C2514-M lenses (Pentax Precision Co., Ltd., Tokyo, Japan) was used to record the experiment (Fig. 2).The captured images were further evaluated in Mercury DIC software (Sobriety Ltd., Brno, Czech Republic), in which the relative displacement of virtual tracked points on the samples was monitored.The location of the virtual point marker was set in the central part of the samples according to Brabec et al. [43].
The displacement during tests in tension parallel to the grain was measured using "clip-on" mechanical extensometers (Zwick Roell AG, Ulm, Germany) (Fig. 1E).Data obtained from the universal testing machine and DIC were synchronised in the MATLAB 2022b (Math-Works, Inc., Portola Valley, California, USA) environment.The shear characteristics in the RT plane were not tested experimentally and were obtained from Eqs. (2.1) and (2.2).
The dimensions and weight of the samples were measured before performing mechanical tests.The dimensions and weights of the compression and bending samples were used to determine MC and density (these samples were selected for their simple geometry, which allowed reliable volume measurement).MC was assessed by the gravimetric method according to Czech standard [44].
In total, 655 samples were tested.Due to the careful selection of defect-clear and specially orthotropic samples, the number of samples was not the same in each group as shown in Table 1.

Establishing constants for the B-MM
The parameters for the B-MM (σ, E, E t , G and G t ) were identified from the experimental stress-strain curves.Each curve from the strain-stress diagrams obtained in the experiments was fitted by a bi-linear function (piece-wise linear function consisting of two linear  parts) in MATLAB R2020b using the lsqcurvefit function, which allows a user-defined non-linear function to be fitted to the experimental data.This produced the optimal bi-linear fit for each normal loading mode (L, R and T) in tension and compression and for shear in two planes (LR and LT).The mechanical parameters from the normal load were additionally marked with " + " for tension and "-" for compression.Figure 3 shows sample data with average bi-linear function in experimental stress-strain diagrams.
The G and σ for the shear RT plane (σ RT ) were not experimentally obtained.Two equations were used from Bachtiar et al. [45] for the calculation of these parameters.
The G RT was calculated by Eq. (2.1): where K is an empirically determined constant to adjust the amplitude of the sinusoidal function for the RT plane ( K = 0.2 ), E −R and E −T are moduli of elasticity in com- pression and µ R and µ T are the Poisson's ratios.
The σ RT was calculated by (2.1) where σ −R and σ −T are yields in compression.
To assess sample deformation optically, virtual trackers were placed both parallel and perpendicular to the load direction.This approach provided both active (compressive) and passive (tensile) strains at the centre of the sample during compression tests in the elastic region.Poisson's ratios were calculated from these strains using Eq. ( 3): where μ represents the Poisson's ratio in the XY plane, ε i is the strain parallel to the loading direction (active), while ε j refers to the strain perpendicular to the loading direction (passive).
In accordance with the notation of direction in Hearmon [34], LR, LT and RT were considered the major Poisson's ratios.Since it was not possible to establish the exact deformations of the sample surface of the compression samples for Trees 2-4, only Poisson's ratios obtained from Tree 1 were used further. (2.2)

Modification of constants for C-MM
The 27 constants describing the B-MM could not be used directly in a FEM solver, since the restrictions of the generalised Hill's theory had to be met.They had to be modified artificially to fulfil the conditions of Eqs.(1.1) and (1.2).The generalised Hill's theory does not distinguish between the E in tension and compression either and this reduces the number of material constants in the model to 24.The stability and robustness of the computations require setting the E t in tension, compression and G t close to zero.This simplifies the definition of these nine constants and reduces the number of constants to 15.More precisely, there remain six shear constants (G and σ in each of the three planes) and nine constants related to tension and compression (for each direction, we have a σ in tension, σ in compression and the E in tension and compression).Only nine of these constants are independent (all E, G and three shear σ).The remaining six σ in tension and compression must satisfy Eqs.(1.1) and (1.2).
The B-MM obtained was replaced by the C-MM which is as close as possible to the original bi-linear curve but has negligible E t and G t .The same lsquarefit function was used in the computation, but we used an artificial upper bound 1.0 MPa for each E t and G t .See Fig. 4 for a graphical illustration of both E-P material models.
The σ for normal loading must satisfy Eqs.(1.1) and (1.2) and thus all constants related to these loading modes had to be modified.We merged tension and compression for each of the three modes into a single diagram, which produced three constitutive laws; see the red curves in the three diagrams on the right-hand part of Fig. 4. To find the bi-linear constitutive law satisfying Eqs.(1.1) and (1.2), we looked for the best continuous piece-wise linear approximation with equal E in tension and compression and with an E t smaller than 1.0 MPa.The optimality criterion was the total sum of the residua across all three diagrams.The computations were performed in the MATLAB R2021b environment.The constitutive law for the resulting computational model is shown in Fig. 4 as a blue curve.
The fixed value of 0.4 MPa was used for all values of E t and G t in the final simulation.This constant was determined by computational testing as the value which does not break the FEM solution.
Because of the requirements of the C-MM, the values of material constants obtained are different from the typical values observed in experiments.Recall that these values should not be interpreted as physical constants of wood.They have artificial values which are a compromise between the restrictions of the generalised Hill's theory and the mechanical properties of wood.

FE model setup
The geometry of the FE model was a 20 × 20 × 300 mm rectangular volume with three rigid body supports in Fig. 4 Stress-strain diagrams for Tree 1 with a comparison of the B-MM by analyses of the experiments and the C-MM used in FEM simulations correspondence with the experimental setup.The model was meshed by ANSYS SOLID95 element type (quadratic 20-node 3D structural solid which allows a generalised Hill's material model, in terms of the ANSYS ANISO feature).The element size of regular sweep mesh 0.003 m was used as a good compromise between detail and the speed of the computations.
The boundary conditions in the form of 0.012 m displacement of the middle support and zero displacements of two rigid supports were applied.Due to the use of the same testing machine and testing of wood with high MC, the parameters for the FE model were set based on the model by Tippner et al. [24].Consequently, the contact pairs between the solid volume and the rigid supports were surface-to-surface with TARGE170 and CONTA174 element types with a friction coefficient of 0.13.According to the experimental setup, the probe for deflection was in the centre of the bottom area, opposite the middle loading point.The solution was divided into 100 steps to ensure the convergence of the solution and recording steps of loading/response.The FEM simulation produced a force-deflection curve which was compared with the experimental one.

FE model analysis and validation
It is natural that some of the material parameters from the C-MM used in the FE model are less significant for the resulting force-deflection curve.The Ansys probabilistic design system (PDS) was used to detect the most significant parameters and to identify correlations among all parameters.
A collection of 1000 computations were evaluated and analysed to obtain statistically significant output from the PDS analysis.The same configuration of the FE model was used.However, this high number of computations necessitated the reduction of computation time which was achieved by increasing the element size to 0.004 m.
The PDS sets random values using the Monte Carlo method for input parameters, runs the simulations and processes the outputs.When enough simulations are performed, a correlation between the output parameters and input data is revealed.However, a completely random setting of parameters is impossible in our model, since there is zero probability that the conditions from Hill's theory are met for random values of parameters.This issue was resolved by recalculating the constants for the R direction.The four σ from the C-MM (tension and compression in L and T) had random values in the PDS analysis and the remaining two σ (tension and compression σ in the R direction) had values derived from formulas Eqs.(1.1) and (1.2).The other values, such as elastic constants and shear parameters, had random values prescribed by the PDS setting.The input range of parameters was defined with a Gaussian distribution by an average value of a parameter, and standard deviation matched the coefficient of variance 0.2 in all cases, which corresponds to the common variability in the mechanical properties of wood.
The following output parameters were considered in the PDS: the deflection and the reaction force at the end of the loading, the deflections at selected values of reaction force (20%, 40% and 75% of the force at the end of loading), the reaction force for selected deformations (20%, 40% and 75% of the deflection at the end of loading), the parameters of the linear approximation of the elastic part of the curve, the parameters of the linear approximation of the plastic part of the curve and the integral of the force-deflection curve with respect to deformation.The Python library Pandas was used to find the correlations suitable for consideration in model verification.These correlations were used to fine-tune the initial phase (the slope of the elastic part) and the terminal phase (the height of the tangential part) of the curve in the bending force-deflection diagram.This resulted in a final set of 27 constants referred to as AC-MM.
The average curves of force-deflection from the experiments were compared with the corresponding curves obtained by running the simulation in ANSYS.Since there were minor inaccuracies in the width and height of the wooden samples, all the force-deflection curves from the experiment were preprocessed and the curves adjusted to a 20 mm × 20 mm cross-section.Thus, the resulting curves are free of inaccuracies caused by minor deviations in cross-section.

Results and discussion
The verified E-P material model of green beech wood was created based on experimentally determined mechanical parameters sorted into a set of 27 material characteristics with their measured variability.Because of the theory limitation of the bi-linear-orthotopic approach for FEM solver (ANSYS ® ), the experimental material characteristics were processed to get constants for the E-P material model.The FE model using optimised AC-MM shows close similarity with experimental curves from three-point bending.
The values of MC and the density of samples in the test time are shown in Table 2.There are the values of compression and bending samples and the values summarising both groups as a more extensive dataset for wood coming from one tree.These data show an 80% MC on average with a standard deviation of 17%.This MC is well above the FSP and, according to Niklas and Spatz [3] should correspond to the MC in the standing tree.Table 2 also shows the values of the basic density with an average of 555 ± 47 kg/m 3 .The MC variability exceeded 10% and averaged 21%; for the basic density, the variability slightly exceeded 10% in some cases.

Material constants for the B-MM
The material constants of the B-MM from processing experimental data for each tree are in Table 3.The experiments in compression R and T were performed for Tree 1 only and the corresponding values are missing from the table.The way we substituted these missing data for the C-MM and adjusted AC-MM is described below.
Figure 3 shows the experimental stress-strain curves with a fitted bi-linear function that captures the average values of B-MM.These parameters for individual samples were determined by fitting a bi-linear function to each stress-strain curve separately, which proved to be an effective tool for evaluating the mechanical properties of wood.Our approach appeared to be more efficient than the generalised approach derived from BS 373 [35] and ČSN 49 0111 [37] used in studies [23,43,46], where the linear part of the stress-strain diagram is specified by section between 10 and 40% or 30-60% of the applied force.
Some of the values of B-MM in Table 3 (namely E t ) reveal a high variability.This comes from the fact that we were dealing with a biological material using a simplified B-MM and the material was under a high load.The parts of the stress-strain diagrams (Fig. 3) where the material fails (from the point of view of the theory of elasticity) are described by E t .This part shows considerable dispersion, therefore, the E t has high variability and significantly low values.This is also visible in the results of studies by Hong et al. [31], Milch et al. [23] and Tippner [24].Especially for the numerically small quantities (such as E t,−L ), we can easily obtain a larger standard deviation than the average value.This, however, does not affect the applicability of the material model, since it is a natural consequence of the fact that values of E t are artificially low in order to keep the possibility of solving the ANSYS model within Hill's theory.

C-MM and FEM simulation
The graphical representation of the experimental data for Tree 1, the B-MM based on these data and the C-MM reflecting Hill's theory assumptions are shown in Fig. 4. Note that the requirements of Hill's theory have a substantial influence on the fit quality.However, even though some of the material values have been adjusted using artificial modifications, there is still agreement between the model and the experimental data.We believe that this is achieved using simultaneous fit across all loading modes and the global approach reflects the real properties of wood better.
The missing data from the experiments with compressions R and T, which failed for Trees 2, 3 and 4, were proportionally replaced by the data for Tree 1.The constant of proportionality was determined by the σ −L , since this number is a significant characteristic in the stress-strain diagram and is available for all four trees.
The PDS analysis was performed for random values of material constants and the results are graphically represented in Fig. 5.This analysis proved a strong correlation between E L and the slope of the force-deflection curve in the elastic part ( r = 0.99 ).The other moduli E R and E T are not correlated with the slope ( r = −0.02and r = 0.02 ).Among the correlations between the material constants and the characteristics describing the terminal phase of the curve, the strongest correlation was detected between σ −L and the force required to get deformation was 8 mm (F 8 ) ( r = 0.89 ).Note also that there is no cor- relation either between E L and F 8 ( r = 0.18 ) nor between σ −L and the elastic slope of the force-deflection curve ( r = 0.02 ).The regressions between the elastic slope and E L and between σ −L and F 8 were used to refine the values of the constants E L and σ −L to get the curve from the FE model as close as possible to the average of the experimental data.

FE model validation
Figure 6 shows the experimental curves of force-deflection for each tree, the average of the experimental curves,      [18] and Pěnčík [47], the curves were checked for point-wise and interval correspondence.For point-wise correspondence, we used the relative difference of FEM simulation from the experimental data.The initial interval (up to 2 mm) was neglected since the values are small and the relative error does not yield useful information.A maximum of the relative error was estimated on the rest of the forcedeflection curves, similarly to Milch et al. [23] and Tippner et al. [24].The values of maximum relative error are 5.87%, 5.60%, 4.21% and 2.56% for individual trees with an average of 4.56%.The integrals of the force-deflection curves over the interval from 0 to 9 mm were used for interval comparison of the curves.The relative errors between the experimental and theoretical curves were 3.56%, 1.24%, 0.88% and 1.03% for individual trees with an average of 1.68%.Three-point bending is a more complex evaluation of the E-P material model than the uni-axial experiment in Hong et al. [19].Despite this, we still have both qualitative and partial quantitative compliance between the FEM simulation and the experimental result.It has been shown that compliance can be improved by adjusting the constants which govern the shape of the force-deflection curve.More precisely, using the correlation from PDS analysis, we were able to adjust the value σ −L to get the terminal part of the computed curve close to the experiment.This output with adjusted data (AC-MM) is shown as a blue curve in Fig. 6.The constants for the AC-MM models are summarised in Table 3 in parentheses.

Conclusions
The material properties of green wood are crucial in developing reliable models of trees under static and dynamic load and for predicting tree failure.We performed experiments on 655 samples to determine the main mechanical characteristics of beech green wood.The data obtained were used to develop an E-P material model, which was verified using a three-point bending test.Comparison of force-deflection curves from the experiment and from the FE model with C-MM show sufficient accuracy (maximum relative error 4.56%, relative error 1.68% for interval comparison).
The experiments resulted in a set of mechanical characteristics with natural variability derived from destructive tests with uni-axial load in compression, tension and shear.The correspondence between the C-MM and the experiment proves that the material characteristics obtained from the experiments with uni-axial loading can be used to build a complex E-P material model.It was verified that the simplified E-P material model consisting of a bi-linear constitutive law for each loading type can be used as a starting point to build a model suitable for numerical simulations in the scope of the generalised Hill's theory.
Important material properties governing a particular case of three-point bending were determined by statistically processing 1,000 numerical simulations in Fig. 6 Force-deflection curves for trees and for an average model.The figures for single trees show the experimental curves, the average experimental curve, the FE model with data from the C-MM and the curve from the FE model with AC-MM according to the elastic slope and force for a deformation of 8 mm.The last picture shows all experimental data as a single set with averages of single trees, the total average of all data and both E-P material models (raw and adjusted) probabilistic analysis.The regression model based on these simulations was successfully used to improve the correspondence between the output of numerical simulation and the real experiment.
The main goal of the research, to build a complex E-P material model which allows the handling of the elastic and plastic deformation of wood, was achieved.For European beech (Fagus sylvatica L.) green wood (80% MC and 555 kg/m 3 basic density-Table 2.) without any defects, we produced a set of validated constants (Table 3) for direct use in structural analyses in a complex loading scenario.These results reveal the real mechanical characteristics of green beech wood, but above all offer a set of constants that can be used primarily in tree biomechanics-for example, in the numerical simulation of static or dynamic loading.However, the use of the E-P model is not limited to the tree biomechanics field, it can be also used in FEM simulations of beech wood with MC above FSP in the wood processing industry applications such as the assessment or processing of beech wood.
Our work proved that an orthotropic E-P material model based on experimental values can describe green wood under load.The use of this E-P material model for beech green wood in tree biomechanics or raw material processing can bring numerical simulations in FEM software even closer to reality.
The procedure for model building we tested is applicable for the development of relevant E-P material models of other wood species under changing conditions (MC, decay).Due to the high natural variability of wood, however, other issues arise, such as the E-P material model of wood with defects caused by rot, knots, etc.Moreover, the applicability of the E-P model can be assessed in wood processing by testing both green wood and wood with MC suitable for structural use, both obtained from the same tree.Testing the developed E-P model in more realistic larger scale cylindrical orthotropy might be the next step in model expansion.Another possible development of this work is to try to simplify the material model.Since the properties in the L direction are very different from those in the R and T directions, it could be useful to examine if it is possible to simplify the C-MM by not distinguishing the latter directions and considering the longitudinal and non-longitudinal directions only.

Fig. 1
Fig. 1 Schemes and dimensions (mm) of testing samples.Samples for A compression tests (1), (3) parallel to the grain (2), (4) perpendicular to the grain; B shear tests (1), (2) parallel to the grain; C tension tests perpendicular to the grain; D three-point bending test; E tension tests parallel to the grain.For brevity, the loading is omitted in pictures A3, A4, B2 and C2 and the dimensions are omitted in A1, A2, B1 and C1.In these cases, the values are identical in all corresponding tests (adapted from Tippner et al.[24])

Fig. 2
Fig. 2 Setup for experimental mechanical testing.Universal testing machine with DIC set of stereo cameras and lights

Fig. 3
Fig. 3 Stress-strain curves merged into a single common diagram.Experimental curves are in grey.The B-MM has been fitted through every single curve and the final B-MM is determined as an average model among the curves considered The last column contains average values for the beech model based on all trees.(The values of σ +R and σ −R in the average model are not averages of single trees but have been calculated from the assumptions of Hill's theory.)For FEM simulation, great accuracy is required since the conditions from Hill's theory are checked with high precision.The C-MM does not distinguish between E +L and E −L and the common values are written as values for tension.The values in parentheses are adjusted values from regression formulas (AC-MM) obtained from PDS analysis the force-deflection curve from the simulation running for the material constants obtained by processing the experimental data (C-MM) as described above and the force-deflection curve for material constants adjusted using regression formulas (AC-MM).The figure shows that we got qualitative agreement between the experimental data and the data from the FEM simulation for all trees.For Trees 1 and 3, we also got quantitative agreement.For Tree 4 the FE model with C-MM (the black curve) is stiffer than the experimental curve (the red curve) and for Tree 2 the FE model with C-MM is less stiff in the elastic part and stiffer in the plastic part of the diagram.So that the agreement of the experimental curves and the FE model curves is not assessed only visually, as stated by Oudjene and Khelifa

Fig. 5
Fig. 5 Heat map with correlation coefficients between E-P material models based on C-MM and output characteristics.Only selected output parameters are shown (elastic and plastic slope, plastic intercept, force at the end of the experiment, integral of force concerning deformation and force required for deformation 8 mm)

Table 1
Numbers of tested samples from the wood of Trees 1-4.The groups of samples are described in Fig.1

Table 2
Values of MC and basic density of compression and bending samples.The average of both types of samples is in the last row.The values in the table are mean and standard deviation

Table 3
Experimental material constants of B-MM and constants of C-MM and AC-MM for FEM simulations for each tree E t,−i Plastic modulus for compression in direction i , where i is one of L, R, and T σ +i Yield in tension in direction i , where i is one of L, R, and T σ −i Yield in compression in direction i , where i is one of L, R, and T G ij Shear elastic modulus in ij plane, where ij is one of LR, LT, and RT G t,ij Shear tangential modulus in ij plane, where ij is one of LR, LT, and RT σ ij Shear yield in ij plane, where ij is one of LR, LT, and RT µ ij Major Poisson ratio for ij plane, where ij is one of LR, LT, and RT F 8Force required for deformation 8 mm in the middle of the bottom side of the beam F END Force at the end of the experiment FE/FEM Finite-element/finite-element method/finite-element analysis +i Elastic modulus in direction i for tension, where i is one of L, R, and T E −i Elastic modulus in direction i for compression, where i is one of L, R, and T E iCommon value of E +i andE −i in FEM model E t,+iPlastic modulus for tension in direction i , where i is one of L, R, and T