Comparison of continuum damage models for nonlinear finite element analysis of timber under tension in parallel and perpendicular to grain directions

A reliable and efficient numerical modelling technique is essential to investigate the behaviour of timber and engineered timber products to promote their widespread use in construction. Wood is an anisotropic material and hence its mechanical properties largely depend on grain direction and type of loading i.e., material behaves differently under compression and tension. Material responses under tension parallel and perpendicular to the grain directions have been reported in the literature but the relevant progressive fracture behaviour has been ignored in typical numerical simulations, due to the complexities and uncertainties around modelling as well as lack of reliable test data. Fracture characteristics play a significant role in analysing crack initiation, propagation, and failure modes of timber so that its full potential can be utilised by knowing the post-elastic behaviour. This paper applies and compares four continuum damage mechanics based constitutive material models (MAT-22, MAT54/55, MAT-143 and MAT-261) available in the commercial finite element software LS-DYNA for simulating the post-elastic behaviour of general timber lamella products. Timber was modelled as both orthotropic and transversely isotropic material to simulate the fracture behaviour in tensile load cases. It is shown that the predicted fracture properties correlate well with experimental data. It was observed that all considered built-in continuum damage models in LS-DYNA are able to simulate the elastic response, but MAT-261, which was originally developed for modelling fibre reinforced composite materials, provides a simple yet reliable option for simulating fracture behaviour of timber.


Introduction
During the last two decades, timber and timber composites have gained increasing popularity in the construction and building industry due to the inherent renewability of raw material and availability of engineered timber products that offer viable options to complement, reduce and, in some cases, replace other conventional materials such as concrete and steel. Engineered wood products such as cross-laminated timber (CLT), laminated veneer lumber (LVL) and glued laminated timber (glulam) are increasingly used in the construction sector (Li et al. 2020). It is widely known that timber is much stronger and stiffer when loads are applied parallel to the grain direction in comparison to perpendicular to the grain direction (Ross 2010). Failure due to tension in parallel and perpendicular directions could be catastrophic due to brittle cracking along the grain which is a relatively common cause of structural damage in timber construction (Danielsson 2013). Although timber design is predominantly dictated by its elastic behaviour, post-elastic response and fracture properties of timber are important aspects for its safe use. Fracture behaviour of timber, especially fracture energy release rate, is affected by a range of mechanical and material properties such as tensile or compressive strength, modulus of elasticity and shear modulus, density, moisture content, orientation of fibres as well as loading type and strain rate. Fracture in timber or in timber-based products is associated with crack initiation and the development of a fracture process zone consisting of multiple micro-cracks ahead of the crack tip. The energy to open and progress crack growth is referred to as fracture energy and represents the total amount of energy released during fracture.
To simulate fracture mechanisms, a wide range of constitutive laws has been proposed during the last two decades (Gharib et al. 2017). 3D constitutive laws based on FE modelling are popular for capturing both the local and 1 3 global response of timber under various types of loading. The development of suitable material models for timber has always been regarded challenging. Timber is relatively ductile under compression, while brittle failure is typically observed under tension. It is an orthotropic material which has different mechanical and fracture properties along longitudinal, radial, and tangential planes. However, properties along the radial and tangential planes are similar and hence, it can be simplified to transversely isotropic without compromising accuracy (Mackenzie-Helnwein et al. 2005). In aerospace applications, Continuum Damage Mechanics (CDM) is well established for simulating the mechanical response of fibre-reinforced composites (Reiner et al. 2019). Since fibre-reinforced composite materials are also orthotropic or transversely isotropic, CDM techniques developed for composites have the potential to be applied to timber. CDM is a mathematical concept, in which the mechanical properties of the damaged material within the fracture process zone are replaced by an equivalent material with associated damage behaviour while maintaining continuity (Reiner and Vaziri 2018a, b). The nonlinear degradation of CDM can be achieved by adjusting the element stiffness with relevant damage parameters when certain failure criteria are met. One of the main drawbacks of CDM is that these models can be highly mesh dependent, especially in modelling the softening response (Sandhaas et al. 2020). An advanced 3D material model for timber described by CDM in tension and shear, and a classic plasticity based model in compression was developed by the US Federal Highway Administration (FHWA) (Murray 2007). The material model was originally developed to simulate wooden guardrail posts against vehicle impact and was implemented as material card  in the commercial FE software LS-DYNA (Hallquist 2007).
A 3D constitutive material model based on CDM for three wood species (Spruce, Beech and Azobe) was recently developed by Sandhaas et al. (2020); the proposed model could capture nonlinear elastic failure modes in compression as well as brittle failure in tension and shear. The material model was developed through a userdefined subroutine and was implemented in the commercial FE software ABAQUS (Abaqus 2011). The failure response of wood was controlled by nine stress-based damage criteria. A simplified technique of modelling anisotropic behaviour of engineered wood based on fibrereinforced composite material was proposed by Valipour et al. (2014), in which timber core was treated as a matrix of a composite, and grains of timber were treated as fibre reinforcement. The material behaviour of matrix was simplified as isotropic, including a continuous piecewise failure envelop, while reinforcing fibres followed a uniaxial stress-strain relationship with specific strength values in tension and compression. Similarly, Gharib et al. (2017) developed an advanced 3D anisotropic constitutive model for fibrous materials with arbitrarily oriented grain directions. Specific failure modes and their related parameters were obtained through decomposition of single surface following Tsai-Wu yield function (Tsai and Wu 1971) in matrix and fibre direction. The uniaxial behaviour of timber in terms of hardening and bi-linear softening responses were achieved by utilizing CDM. The proposed material model was subsequently implemented in ABAQUS as a user defined subroutine.
Numerical modelling of timber or timber structures largely depends on the accuracy of material parameters, in particular modulus of elasticity, failure strength and fracture energies. Fracture energy is one of the key parameters to achieve a realistic post-elastic response. Fracture energy values for different materials are experimentally obtained using various fracture tests. Mode I fracture energy can be measured by various test configurations, such as, compact tension (CT), double cantilever beam (DCB), small scale notched specimen and single edge notch beam (SENB) tests (Boström 1994), more details are mentioned in Sect. 2.1. Stanzl-Tschegg et al. (1995) conducted mode-1 fracture tests on spruce (softwood species) in tangentiallongitudinal (TL) and radial-longitudinal (RL) plane. A new fracture splitting test method was developed to determine specific fracture energy G f, which ensures stable crack propagation (Reiterer et al. 2002). Dourado et al. (2015) characterised the mode I fracture of spruce timber by single-edge-notched beam loaded under three-point bending (SEN-TPB). A resistance-curve based on linear elastic fracture mechanics was proposed by applying data reduction schemes accounting for stress relief near the crack without monitoring crack length compared to classical methods.
Relevant literature shows that CDM has been used in FE modelling of timber but modelling of crack propagation and fracture behaviour of timber using CDM is scarce. Although efforts (Gharib et al. 2017;Sandhaas 2012;Sandhaas et al. 2020;Valipour et al. 2014) were made to develop new material models to simulate elastic-plastic and fracture behaviour of timber, most of these models were developed by user-defined subroutines or other in-house codes that are not available for general use. To illuminate such deficiencies, a rapid analysis technique with easily accessible material models is essential. LS-DYNA contains a large set of built-in material models for composite materials that could potentially be suitable for timber. In particular, it includes a CDM-based material model known as MAT-143 (Murray 2007) for wood and MAT-261 (Pinho et al. 2006) for fibre reinforced composites which have the potential for modelling post-elastic behaviour of timber. In this study, three different CDM based composite material models along with wood material model MAT-143 are addressed and compared with experimental test result.

Material properties of timber
Timber is an orthotropic material where mechanical properties can vary along its three orthogonal directions. It is described by three perpendicular planes, commonly known as longitudinal (L), radial (R) and tangential (T) as illustrated in Fig. 1a. Material strength is significantly higher in longitudinal direction, commonly known as parallel to the grain direction, than those in radial and tangential directions as shown in Fig. 1b. Orthotropic elastic behaviour of timber may be described by constitutive material laws including twelve constants, which are denoted by E LL , E RR , E TT , G LR , G RT , G TL , LR RT , LT , TR , RL and LR (as shown in Table 1),where E, G and ν represent modulus of elasticity, shear modulus and Poisson's ratio, respectively (Bodig and Goodman 1973). According to Hooke's law, material properties of timber are characterized by the compliance matrix shown in Eq. (1).
(1)  where, refers to strains, refers to stresses. The Poisson's ratios (ν RL, ν TL ) are very small, less than 0.1, thus difficult to measure and, including these ratios ν RL, ν TL and ν LT can be established from relation with modulus of elasticity MOE. The relationship between MOE and Poisson's ratio can be expressed as Therefore, twelve elastic constants are reduced to nine, i.e., three moduli of elasticity (E LL , E RR , E TT ), three shear moduli (G LR , G LT , G RT ) and three Poisson's ratios (ν LR , ν RT . ν LT ). Additionally, material properties of timber in radial and tangential directions are somewhat similar, and hence, timber may again be simplified to transversely isotropic material with five elastic constants such as E LL , E RR , LR , G LR and G RT considering To account for the variability of material properties, a thorough literature review was conducted as part of the current study to identify all elastic orthotropic parameters for different species of spruce grown in Europe as shown in Table 1. Table 1 shows that the modulus of elasticity parallel to the grain, E LL , is significantly higher than the other two moduli in radial and tangential direction. The shear modulus in the radial-tangential plane, G RT is considerably lower than the other two shear moduli, which often results in rolling shear failure (Li et al. 2021). Mechanical properties of timber also vary significantly with the type of loading. In tension, timber shows linear softening behaviour accompanied by brittle failure, whilst in compression timber is relatively ductile; details of stress-strain responses of timber are discussed in Sect. 3. Although all experimental specimens were reported to be produced from straight grain, i.e., without knots, cross grains, checks, or cracks, the reported strength values of spruce vary substantially as shown in Table 2; this type of variation is caused by density, orientation of growth ring, and moisture content (Dahl 2009). (2)

Fracture Properties of Timber
The fracture behaviour of timber is primarily influenced by micro-cracking, separation of front of crack tips, fibre bridging, strength and, most importantly, the total energy released in fracture processing. The fracture energy is defined as the work done by fracture and can be expressed as the amount of energy dissipated at the tip of the crack during the gradual development of the fracture zone of a unit area (Reiterer et al. 2002). To determine the fracture energy G f the area under the load-displacement curve is divided by cross-sectional area (A) of the fracture surface. Fracture energy may be calculated by: where, F is the load required for fracture, is the corresponding displacement, b is the width of the specimen and a is the crack length.
In anisotropic materials such as timber, the energy needed for crack propagation depends on the orientation of the crack plane as well as the method of load application. Considering three orthogonal directions and crack propagation, fracture energies in LR, RL, TL, LT, RT and TR planes generally differ, as illustrated in Fig. 2a (Danielsson 2013). Crack propagation systems are usually indicated using two letters, the first letter indicates the direction perpendicular to the crack plane, whilst the second letter refers to the direction of crack propagation. In-depth investigation of the fracture process in timber should typically comprise three different modes of crack propagation such as Mode-I, Mode-II and Mode-III as shown in Fig. 2b. These distinct fracture modes are defined as follows: in Mode-I, the crack surface deforms directly in opposite directions due to tension; in Mode-II, crack surfaces slide against each another due to in-plane shear and, in Mode-III, the crack surface is teared due to in-plane loading.
Mode-I and Mode-II are the typical scenarios of fracture progression in timber structures. Since behaviour of timber under tension is the primary focus of this paper, Mode-I fracture condition was chosen for detailed investigations.
(3)  Table 3 presents a summary of relevant studies on fracture energy of spruce when subjected to tension (Mode-I fracture). Fracture energy is one of the few available parameters required to predict the post-peak softening behaviour of timber under tension parallel and perpendicular to the grain direction.

Finite element modelling of timber
Single finite elements of timber were simulated under uniaxial tension to achieve a clear understanding in identifying differences in built-in material models currently available in LS-DYNA. Numerically obtained results from simple tensile tests were then compared with experimental data obtained from the literature. Finally, fracture energy values of timber were determined from numerical studies and compared with those reported in corresponding experimental investigations.

Material modelling
Brittle fracture behaviour of timber may be simulated using various FE approaches such as lattice model, classical continuum mechanics, fracture mechanics, cohesive zone modelling (CZM) and CDM (Sandhaas 2012). CDM is a nonlinear method and follows a simple and straightforward approach to model softening behaviour of a material. Continuum damage models are dictated by damage initiation criteria and damage evaluation laws defined in specific directions. Three-dimensional stress-strain constitutive laws are adopted by means of a damage parameter d , where material degradation or non-linear softening behaviour is achieved by transforming the stress tensor d using the relationship shown in Eq. (4). The damage index d ranges from 0 to 1 indicating no damage and full damage saturation, respectively (Sandhaas et al. 2020). Figure 3 illustrates the material degradation from initial stiffness E 0 to damaged stiffness E 0 (1-d) after activating damage initiation at strain 1 and stress T . With the initiation of failure, the relative stress components are multiplied by (1-d) for every time step to reduce the effective stress to zero. At the final point, stiffness of the element becomes zero and the element is removed from the model.
Material degradation can be categorised as sudden and gradual degradation models. In sudden degradation models, stiffness of the material immediately reaches zero after fulfilling the damage initiation criterion, following the path OAC in Fig. 3. One of these models follows Chang-Chang's failure formulation (Chang and Chang 1987b). Contrarily in gradual degradation models, material stiffness reduces to a series of lower stiffness values following Eq. (4) before reaching zero; a typical example of gradual degradation path is shown by OAB in Fig. 3. Several constitutive material models are available in   commercial FEM software packages that use both sudden and gradual degradation models. LS-DYNA contains several builtin orthotropic and transversely isotropic material models that are suitable for timber. In LS-DYNA, MAT-143 and MAT-261 are based on gradual degradation, whereas MAT-22 and MAT-54/55 follow sudden degradation. This paper compares the performance of these four built-in material models for timber when subjected to tension. The following subsections describe key aspects of the material models considered herein.

Wood material model: MAT-143
Material model MAT-143 (Murray 2007) in LS-DYNA was developed for wood considering it as a transversely isotropic material. Parallel direction represents longitudinal (L) axis and perpendicular direction refers to radial (R) or tangential (T) axis of timber. Hashin's formulations (Hashin 1980) were used to define different strength criteria. Since this study investigates timber under tension parallel and perpendicular direction, only brittle tensile failure of Hashin's model has been included in Eqs. (5) and (6). In parallel to the grain direction, damage initiates when Eq. (5) is satisfied, whilst Eq. (6) must be true to initiate damage in perpendicular to the grain direction.
where LL , LR , LT , RR , TT , RT are the stress components for orthotropic material in relevant planes; f Lt and f LT are tensile and shear strength in parallel to the grain direction, respectively. f Rt and f TR are tensile and shear strength in perpendicular to the grain direction, respectively. Damage evolution follows separate rules for parallel and perpendicular to grain directions. Failure in parallel direction is often considered to be catastrophic since timber cannot carry any load after failure in this case. If damage occurs in parallel mode, all six-stress components are degraded gradually. However, if damage occurs in perpendicular direction, only perpendicular stress component is degraded. Based on the following assumptions, MAT-143 is represented by the degradation models shown in Eqs. (7)-(10).
where B, C, and D are softening parameters to be obtained from experimental data. The parameter d L(max) and d RL(max) is the maximum damage in parallel and perpendicular direction. It ranges from 0 to 1, where 0 means no damage and 1 means maximum (saturated) damage. Here G fLT and G fTL represent fracture energy in parallel and perpendicular direction and G f 0 denotes initial fracture energy. Fracture energy G f plays an important role in progressive damage modelling and should be applied uniformly to all elements, which is implemented by using separate mesh-size regulations for MAT-143 when used for parallel and perpendicular models (Murray 2007).

Continuum damage model: MAT-261
To evaluate damage in unidirectional composite material, a CDM based material model MAT-261 was developed by Pinho et al. (2006). Two separate failure criteria are used in this model; one for the fibres under tension, and the other is based on Mohr-Coulomb criteria for the matrix surrounding fibres (Manan et al. 2020). Equations (11)-(15) show the failure criteria for fibres and matrix under tension.
where, is the angle of fracture surface measured in through thickness direction. Damage in MAT-261 follows a linear progression law, in which softening behaviour occurs with the gradual reduction in effective stress. Damage evaluation laws are shown in Eqs. (16) and (17).
where, f XY is the material strength in XY direction, and 0 is corresponding initial strain.
The final damage parameter of MAT-261 is governed by the failure strain 2 , which is directly related to the fracture toughness Γ and the characteristic length L of the element. Fracture toughness Γ can be calculated from the fracture energy G f following Eqs. (18) and (19).

Composite damage models MAT-22 and MAT-54/55
Composite damage model MAT-22 is applicable to both shell and solid elements. This model was developed based on Chang-Chang's failure criteria (Chang and Chang 1987b) and can effectively simulate linear elastic behaviour with brittle failure in tension and shear. MAT-22 only requires strength parameters, there is no need for fracture energy or failure strain to be input. MAT-54/55, which is an updated version of MAT-22, is built by combining both Chang-Chang's and Tsai and Wu (1971) failure models. MAT-54/55 has definite progressive damage criteria as maximum strain failure for tensile fibre in parallel and perpendicular directions.
Tensile failure in parallel to the grain is triggered when Eq. (20) is satisfied.
where β is a shear stress weighting factor.
In perpendicular to the grain direction, tensile damage is initiated when Eq. (21) is met.
A summary of important characteristics for the afore mentioned material models is presented in Table 4.
It is worth noting that all material models reduce to the same maximum strength criterion for damage initiation in pure tensile loadings, and hence the initial elastic response is identical. Damage evaluation rules are, however, different in each of the considered material models. MAT-143 and MAT-261 follow gradual degradation using different damage evaluation rules. For MAT-143 damage index d 143 depends on fracture energy G and constants B, C and D that should be obtained from experiments as shown in Eqs. (9) and (10). In the case of MAT-261, damage index d 261 depends  Chang and Chang (1987a) Solid and Shell on failure strain, which is a function of fracture toughness as shown in Eqs. (16) and (17). Mechanical properties of timber can vary significantly as they depend on several factors such as grain orientation, growth, geographical location, age of the tree, temperature, and moisture content (Gharib et al. 2017). Material input data used in FE models developed as part of the current study were taken directly from the corresponding reported literature to account for the inevitable inherent uncertainty and considerable variability in timber properties. Material properties used for the following FE simulations were categorised in four different material IDs as 1, 2, 3 and 4 which were taken from four different sources reported in the literature as shown in Table 5. Here only orthotropic material properties are presented which are used in MAT-22, MAT-54/55 and MAT-261. MAT-143 is transversely isotropic thus only five elastic constants E LL , E RR , LR , G LR , G RT are used from Table 5. However, strength and fracture energy parameters are identical in orthotropic and transversely isotropic material.

Single element simulation
Single element simulation is an efficient option to explore different material models as it becomes difficult to directly evaluate the effect of material models in full scale analyses (Reiner and Vaziri 2018a, b). Although MAT-143 was developed for wood, other material models used in the current study such as MAT-261, MAT-22 and MAT-54/55 were originally developed for fibre reinforced composite materials. Hence, it is important to evaluate the behaviour of all material models under identical loading conditions prior to developing full scale models. To evaluate the fundamental similarities and differences of each of the considered material models, investigation of a single three-dimensional solid element under tensile loading in both parallel and perpendicular to grain directions was carried out. An eight-node solid element with unit (1 mm) length and width as shown in Fig. 4a and b was used for simulation where the prescribed displacement was applied parallel and perpendicular to the timber grain direction. The degrees of freedom in vertical direction at the base of the sample were constrained, and the displacement was applied on the top surface, as shown in Fig. 4. Material properties used in this model were taken from Material ID-1 (Table 5) and obtained results are shown in Fig. 5.
As discussed in Sect. 3.1, MAT-143 and MAT-261 use damage parameters based on fracture energy values and fracture toughness, respectively. In MAT-261, the damage index d 261 follows a linear progression law, which is established from fracture toughness, material strength and length of element as shown in Eqs. (16) and (17). MAT-143 also depends on its relevant damage index d 143 , which is formulated not only from the fracture energy but also from post-elastic softening curve fitting parameters B, C and D as shown in Eqs. (7)-(10). According to FHWA (Murray 2007) guidelines, higher value of D produces flatter peaks while larger value of C produces severe softening. For the current single element modelling and all other simulations, values for B, C and D of 30, 10 and 30 respectively, were taken from softening parameters reported for southern yellow pine by FHWA (Murray 2007). To ensure the same  Fig. 5a and b that MAT-261, MAT-143 and MAT-55 produced the same amount of energy as determined by the area under the stress-strain curves. Since all initiation criteria for materials presented above reduce to a maximum stress criterion in single element models under pure tensile loads, the elastic response and strength predictions were identical for all considered material models. It is worth noting that MAT-55 and MAT-22 exhibited different post-elastic responses as the damage evaluation techniques for those materials are different. Although MAT-55 showed a sudden degradation mode, failure initiates at the maximum stress and completely fails when it reaches the maximum strain. MAT-22, on the other hand, showed a sudden failure after reaching the ultimate stress. The magnitude of fracture energy for MAT-22 is lower than the other material models used in this current research.

Simulation of different coupon level tests
Several different experimental studies from the literature were selected to evaluate the performance of considered Significantly different experimental cases such as tension parallel and perpendicular to the grain, progressive fracture tests in Mode-I using single edge notch and micro wedge splitting tests were considered in the current study to evaluate the relative performance of the considered material models. This section presents FE modelling techniques used to simulate those loading and geometric conditions, followed by mesh sensitivity analyses for each case. Figures 6 and 7 show experimental set ups used by Karagiannis (2017) for tensile test specimens in parallel and perpendicular to the grain Fig. 6 Tensile test in parallel to the grain direction Specimen-1 (S1) Fig. 7 Tensile tests in perpendicular to grain direction, Specimen-2 (S2) directions. One edge of the sample was fully clamped whilst the displacement was applied at the free end at a rate of 1 mm/min. The single edge notched beam (SEN-TPB) test involves a relatively simple mechanical setup, which was used by Dourado et al. (2015) for Norway spruce. 1 mm wide precrack notch of length a 0 was fabricated at the centre line of the beam as shown in Fig. 8. Material properties for FE simulation were taken from Dourado et al. (2015), which is listed as Material ID-4 in Table 5. Figure 8a shows geometric dimensions of the beam and Fig. 8b shows the corresponding FE model developed in the current study. Load at the pin was applied until the beam ruptured completely along the crack.
The third case considered for FE modelling was the micro wedge splitting test also used for Mode I fracture testing. Figure 9a shows the test arrangement used by Frühmann et al. (2003) to conduct the micro-wedge splitting fracture test on spruce and beech, whilst Fig. 9b shows the FE model developed as part of the current study to evaluate performance of different material models in simulating the reported mechanical behaviour.

Mesh sensitivity analysis
Most CDM based material models are highly dependent on element size and mesh regularity, particularly those used for failure modelling. Mesh uniformity also is important to ensure uniform distribution of fracture energy to all elements. It is worth noting that MAT-143 and MAT-261 can distribute fracture energy uniformly to all elements (Hallquist 2007), whilst MAT-22 and MAT-54/55 do not require fracture energy as material input data. Although the considered four material models have different damage evaluation parameters, all of them are CDM based showing similar patterns in elastic response and failure initiation. Hence, mesh sensitivity analyses were conducted only using MAT-261 assuming that other material models will converge in a similar fashion. Mesh sensitivity analyses were conducted for all considered loading cases i.e., tension (S1 and S2), SEN-TPB (S3) and micro-wedge splitting (S4). Load and boundary conditions were applied according to the relevant experiment as illustrated in Sect. 3.3. Results shown in Fig. 10 demonstrate that higher mesh density produced more accurate results but was computationally expensive. Figure 10 also shows differences in load-deformation response in considered cases in the current study due to change in element size; Table 6 lists the optimal size of elements as obtained from FE results.
Specimen-1 and-2 were subjected to tension in parallel and perpendicular to the grain directions resulting in uniaxial stress and uniaxial strain conditions. Hence, maximum stress vs number of elements for S1 and S2 were evaluated for the mesh convergence study, and it is obvious from Fig. 10 that at 4 mm element size, both Fig. 10 Mesh sensitivity analysis using MAT-261 in four different simulation cases S1 and S2 have converged. In Samples 3 and 4, model geometry and loading were complex due to the presence of the notch used to investigate progressive fracture. Load-deformation responses obtained from FE models of S3 and S4 were compared against those reported from experiments to study mesh sensitivity. Numerical analysis of micro-edge splitting test (S3) showed high mesh sensitivity, and even developed unstable crack propagation due to inadequate meshing. Hence, finer mesh with 0.5 mm element size was used, which eventually produced stable load-displacement response in the FE model developed to simulate micro-edge splitting test as shown in Fig. 10c. In SEN-TPB test, i.e., S4, 2 mm element size produced good agreement with test results as shown in Fig. 10d. Table 6 lists mesh convergence results and the proposed element sizes for all loading cases considered in the current study.

Tension in parallel and perpendicular to the grain directions
Timber shows linear brittle behaviour under tension in both parallel and perpendicular to the grain directions. As a natural bio-composite material, considerable variations are often experienced in strength, MOE and failure strain during experimental investigations (Sirumbal-Zapata et al. 2018). It is evident from relevant literature that tensile strength, MOE, and failure strain in parallel to grain direction are considerably higher than those in perpendicular to grain direction; Table 7 lists the experimental values for these parameters reported in the literature. Test results showed a comparable mean with the characteristic value of C30 strength class from EN 338 (Standard 1995). Average values of MOE and material strength were taken from reported experimental results and used as input parameters in the current study.
Fracture energy values for spruce were taken from Daudeville (1999) which also is a C30 strength class of EN 338 (Standard 1995). Detailed material properties are presented in Table 5 as material ID-3. Figure 11 shows a typical contour plot of maximum principal stress and damage pattern developed in S1 and S2 samples when subjected to tension parallel and perpendicular to the grain. The support and loading conditions are applied as described in the corresponding experiments, which is   Figure 12 compares numerically obtained stress vs strain responses under tension in parallel and perpendicular to grain directions to those reported by Karagiannis (2017).
Detailed comparisons between experimental data and the corresponding FE results are shown in Table 7 which illustrates that all considered material models indicated good agreement with the experimental results for timber under tension in both parallel and perpendicular to the grain directions. The average value of modulus of elasticity in parallel direction E LL and E RR from experimental tests were recorded as 10,494 MPa and 181.1 MPa. Similarly ultimate strength in parallel direction LL and perpendicular direction RR are 36.2 MPa and 1.55 MPa. Table 7 illustrates that all considered material models indicate good agreement with these average experimental measurements for timber under tension in both parallel and perpendicular to the grain directions. The ratio of the ultimate stress obtained from finite element analysis and from average test data in parallel and perpendicular direction are in the range of 0.99-0.94. It is worth noting that simplified material models such as MAT -22 and  MAT-55, which do not require fracture energy data as input, can produce reliable predictions under tensile loading.
As part of the current study, numerical simulations for higher strength spruce were also conducted and numerically obtained results were compared with those reported by Franke (2008) and Poulsen (1998). Uniaxial tension tests in parallel to the grain direction were conducted by Franke (2008) on clear timber dog bone shaped specimens of 2 mm width and 20 mm length. On the other hand, uniaxial tension tests perpendicular to the grain were carried out by Poulsen (1998) on a 20 mm wide and 20 mm long clear wood sample. Experiments were conducted following the European standard EN 408 (EN 2003) with 12-13% moisture content. From the test the average failure strength in parallel and perpendicular to the fibre direction was recorded as 61 MPa and 3.5 MPa. Material ID 2 (as shown in Table 5), which was reported from the test, was used in the current study to develop corresponding FE models. Figure 13 compares stress-strain response of tested samples against those obtained numerically for all considered material models. Figure 13a shows that there was considerable variability in parallel to the grain test response, but all material models produced an almost identical response, which is in line with findings from previous comparisons presented here. In perpendicular to the grain direction, test responses were more consistent as shown in Fig. 13b. In FE simulations, MAT 22 failed to capture the complete stress-strain response in perpendicular to the grain direction as the failure criteria are solely defined by strength, whereas other material models MAT-143, MAT-261 and MAT-55 captured the post elastic zone more accurately as these material models contain post-peak softening options in form of fracture energy or failure strain as input data. This clearly demonstrates the significance of including post-peak softening capabilities to provide reliable FE models for simulating damage evolution in timber.

Progressive fracture tests
Among various testing methods, the Mode-I fracture evaluation by single-edge-notched beam (SEN-TPB) testing is widely used for stable crack propagation and simple loading arrangements. Micro-wedge splitting tests for Mode-I fracture, on the other hand, can produce unstable crack propagation and complex loading conditions (Frühmann et al. 2003).
Both test types were modelled as part of the current study to evaluate the effect of material models in FE simulations under different loading arrangements. From the comparison of predicted and experimentally measured fracture energy values presented in Tables 8 and 9, it is obvious that fracture energy values are significantly higher in micro-wedge splitting tests than SEN-TPB tests considering the same material.  Figure 14a shows the failure mode and deformed shape of the single edge notched beam obtained from FE simulation using MAT-261. The maximum displacement of the notched beam was recorded as 3.41 mm before ultimate rupture, i.e., specimen broken into two pieces. Fracture energies from FE were calculated using Eq.
(3) and the values are presented in Table 8. The area under the load vs displacement curve was calculated using the "trapz" function available in MATLAB. Calculation of fracture energy requires corresponding crack length, which was determined from FE models. Detailed load-displacement data at the mid-span obtained for all types of material models are illustrated in Fig. 14b and were compared against experimental measurements. The ratios of FE results and experiments on failure strength are 1.03 for MAT-143, MAT-261, MAT-55 whereas MAT-22 yields a comparatively lower ratio of 0.87. As the fracture energy is calculated from the area under the curve of load-displacement, it also followed a consistent FE and  (Frühmann et al. 2003) and numerical results obtained for micro-wedge splitting tests are presented in Fig. 15 and Table 9. It was observed that in micro-wedge splitting test, results obtained using MAT-143 were not consistent with other material models used in the current study. Considerable local failure was observed in the FE simulation using MAT-143, which resulted in the sample to incline on one side as shown in Fig. 15b. As a result, crack propagation was unstable, which is obvious from the stepwise post-peak softening response in Fig. 15a. This ultimately produced significantly higher fracture energy (903.23 N/mm) when compared to those obtained from other material models. Figure 15b depicts timber failure at the steel pin support for MAT-143 and Fig. 15c shows typical failure observed for MAT-261, MAT-22 and MAT-55. MAT-22 performed relatively better in simulating micro-wedge splitting tests in comparison to SEN-TPB tests as the former process results in rapid softening. Fracture process zone and crack tips from the experiment and those captured from the FE model using MAT-261 are shown in Fig. 16. As explained in Sect. 3.1, post-peak responses of the considered material models can be divided into two groups such as sudden degradation and gradual degradation. This fundamental criterion is obvious in FE simulation results for both notched beam and micro-wedge fracture tests. As damage initiation criteria for MAT-143, MAT-261 and MAT-54/55 in tension reduce to the same threshold, they produced almost identical responses until the onset of failure where the maximum failure loads obtained from three material models were very close, i.e., 106.32 N, 105.54 N and 106.04 N, respectively. FE models involving these three material models showed similar damage progression before final rupture and released consistent amount of fracture energy showing good agreement with average experimental data with differences of 16%, 18% and 19% using MAT-143, MAT-261 and MAT-54/55, respectively. MAT-22, on the other hand, produced a lower peak force of 89.9 N and a considerably smaller amount of fracture energy when compared against other material models. This finding is in line with observations for single element simulations and tensile tests outlined in previous sections. This clearly shows that MAT-22 may be used to simulate the mechanical response of timber up to the ultimate strength but is not suitable for modelling post-peak behaviour as it cannot capture the fracture response realistically due to the nonexistence of progressive damage evolution.

Conclusion
The current study was designed to evaluate the suitability of built-in CDM based material models for the numerical analysis of progressive damage evolution in timber using the commercial finite element software LS-DYNA. Four different types of material models, i.e., MAT-143, MAT-261, MAT-54/55 and MAT-22 were used for simulating different types of loading cases. The following are the key outcomes of the current study: 1. When timber is subjected to standard tensile tests, either in parallel or perpendicular to the grain direction, all considered material models MAT-143 MAT-261, MAT-54/55 and MAT-22 showed good agreement with test results up to the ultimate load. This initial investigation demonstrates that despite being developed to model composite materials, MAT-261, MAT-54/55 and MAT-22 may be used to model timber as an orthotropic or transversely isotropic material under tension to investigate elastic and strength properties. 2. MAT-143 MAT-261, MAT-54/55 are well suited to represent progressive fracture behaviour and can capture the overall nonlinear response of timber under progressive tensile loading as those models require post-peak criteria including fracture energy, failure strain or other damage parameters as input to capture brittle or gradual softening of timber. In contrast, MAT-22 showed some limitations in simulating progressive fracture in timber as there is no effective damage evolution mechanism to produce a gradual softening response. Therefore, MAT-22 may be used for simulating elastic response of timber only. 3. MAT-143 performed well in simulating progressive fracture in timber, as expected. However, as MAT-143 is a transversely isotropic material model, material axes combinations of longitudinal (L) to radial (R) or longitudinal (L) to tangential (T) showed good agreement with experiment but radial (R) to tangential (T) combination are not yet proven to be effective as MAT-261. Additionally, MAT-143 is dependent on softening parameters B, C and D, which are dependent on experimental data, and hence may vary significantly for different species or even for the same timber species. As a result, MAT-143 requires a larger amount of input data to be obtained from experimental tests compared to other material models; this can be difficult and expensive considering timber as a natural bio-composite material possessing inherent imperfections. 4. MAT-261, which was originally developed for fibre reinforced composite materials, was found to be effective for simulating tensile and fracture behaviour of timber. MAT-261 did not show any issues on any combination of timber material axis, i.e., radial (R) to tangential (T) or longitudinal (L) to radial (R) or tangential (T). This orthotropic CDM based material model offers a relatively simple yet reliable technique for simulating fracture behaviour of timber.
These efficient material models can be potentially used for rapid analyses of various large-scale timber structures such as cross-laminated timber, glue-laminated timber, and other engineered timber products. To describe the full sequence of failure modes, future work will extend the presented study to include damage evolution in structures subjected to compressive loadings and interlaminar failure in form of delamination.
Funding Open Access funding enabled and organized by CAUL and its Member Institutions.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.