Multi-Surfaced Elasto-Plastic Wood Material Model in Material Point Method

Wood is a naturally occurring material widely used for construction. Due to its natural origin, wood properties vary and its behaviour is complex. This paper shows an implementation of a multi-surface elasto-plastic constitutive material model for wood into a custom explicit material point method code. The constitutive model chosen is one proposed by Schmidt & Kaliske with minor modifications to ensure better internal consistency. The model parameters are chosen based on literature data for spruce. The paper presents two Convected Particle Domain Interpolation Material Point Method simulations of experiments, both performed with the previously established model parameters. The first simulation replicates a compression test of a spruce specimen perpen- dicular to grain direction, carried out at the Department of Civil Engineering, Aalto University. The second simulation replicates an experiment from literature, in which a spruce specimen with knots is tensioned until failure. The numerical simulations successfully replicate the experimental outcomes qualitatively in terms of the deformation and load-displacement curves. Simulations of the three knotted specimens under tension, with introduced slight variation in wood grain direction, replicate different failure patterns with a similar failure load, resembling the behaviour of natural wooden structural elements. Additionally, one of the obtained failure pat- terns replicates that of the experiment well.


Wood material
Wood has been used as a construction material for centuries. With growing concern on improving sustainability and lowering the environmental impact, wood has become even more appealing due to its renewability and often low carbon footprint of the overall production process. Yet, wood behaviour is complex as it is best approximated by an orthotropic material model where the mechanical properties such as strength and stiffness vary depending on grain and loading directions. Furthermore, wood is relatively brittle in tension as opposed to the significant ductility in compression. Finally, wood behaviour is affected by, among others, temperature, moisture content, and heterogeneity.
Until the mid-1980s, wood testing comprised mainly of uniaxial loading. The mixed loading scenarios were introduced by e.g. Hemmer (Hemmer, 1985), Ehlbeck (Ehlbeck, 1986), and Spengler (Spengler, 1986), and later, Eberhardsteiner et al. (Eberhardsteiner et al., 1999), Eberhardsteiner (Eberhardsteiner, 2002), and many others evaluated wood strength under various biaxial loading conditions and showed that the tensile loading beyond the wood tensile strength results in an instantaneous loss of strength at the macroscopic level. Yet, at the microscopic level, Faessel et al. (1999) discovered the gradual strain softening behaviour during tensile failure. Subsequently, Mackenzie-Helnwein et al. (2003) suggested a connection between the energy dissipated by this mechanism and the mode I fracture energy, both for the longitudinal and radial directions. Further fracture tests on all possible orthogonal crack propagation directions were conducted by King et al. (1999), who comprehensively experimented on fracture toughness of pines. Other similar studies for spruce and pines are also available, i.e. (Aicher, 1994;Stanzl-Tschegg et al., 1995;Reiterer et al., 2002;Dourado et al., 2008;Tukiainen and Hughes, 2013;Tukiainen and Hughes, 2016). On the other hand, compression in the transverse direction beyond the compressive strength results in a significant plateauing followed by an increase in stiffness, a phenomenon called densification. Compression in the longitudinal direction ultimately leads to an inelastic response -initial softening due to grain buckling is followed by a longer plateau and densification afterwards (see Adalian and Morlier, 2001;Adalian and Morlier, 2002 for the experimental results). Therefore, while the tensile failure results in a quasi-brittle behaviour, the compressive failure exhibits a significant strain-hardening and ductility.
Numerical modelling of the anisotropic material failure is challenging under multi-axial loading because of the interaction between the stress components in the orthogonal directions. The historical effort by Tsai and Wu (1971) to model such material utilizes an elliptical failure surface, which at some point was tried on wood (Eberhardsteiner, 2002). However, as wood failure modes are easy to identify, a multi-surface approach, in which each yield surface is associated with a failure mode is an attractive approach with a clear physical interpretation. Such multi-surface plasticity constitutive model for wood is conceptually similar to the models of Tresca and Mohr-Coulomb, both featuring six surfaces representing all permutations of 2 out of 3 principal stresses. Discussions on Tresca and Mohr-Coulomb are widely available in many solid mechanics or computational plasticity textbooks, e.g. (Simo and Hughes, 1998;de Souza Neto et al., 2011;David and Zdravkovic, 1999). For wood, the multi-surface approach has been used e.g. by Mackenzie-Helnwein et al. (2001). The additional yield surfaces allow for defining different nonlinear behaviour for the radial and longitudinal failures. Mackenzie-Helnwein et al. (2003) then implemented more surfaces to account for failure phenomena in the radial and the longitudinal compression and tension. Kaliske (2006, 2009) developed the multi-surface concept further for a three-dimensional case. They introduced two additional surfaces for tension and compression in the tangential direction, and another surface exclusively for the shear failure. The equations were generalized, practically allowing the same set of expressions to describe all seven surfaces, differentiated by different parameters.
Considering the wood inherent orthotropic behaviour and the different failure behaviours when subjected to loading, the material model by Kaliske (2006, 2009) seems to be the most capable for the task. Unfortunately, the numerical implementation quickly revealed that the model produced an inconsistent stress increment for an elasto-plastic deformation when different units were used. In this paper, we identify the dimensional inconsistency in the model formulation of Kaliske (2006, 2009) (Appendix A) and present a necessary modification, following the established principles of the material model (Appendix B). The proposed material model also utilises an additional algorithm to tackle the issue of severe softening, potentially arising from the highly-brittle failures (Appendix C). In this paper, we assume a constant moisture content and temperature and take wood as a homogenous material. Simulating the effects of changing temperature and moisture remains an open possibility for future studies.

Material Point Method
In this work, numerical simulations utilize Convected Particle Domain Interpolation (CPDI) (Sadeghirad et al., 2011) Material Point Method (MPM). In Material Point Method, originally introduced by Sulsky et al. (1994), the computational domain is discretized by particles cast over a background grid. Particles store all kinematics and internal variables, while the equations of motion are solved at the background grid. Interpolation functions, similar to shape functions in Finite Element Method (FEM), map the variables between the particles and the grid. The particles move during the computations, while the background grid remains undeformed between timesteps.
Compared to FEM, MPM performs better when involving problems with large deformations, as large element distortions in FEM may result in ill-conditioned Jacobian matrices (see e.g. Więckowski, 2004;Beuth et al., 2011). Nonetheless, the original MPM formulation leads to errors when particle crosses cell boundaries (known as cell crossing noise). The problem is significantly reduced by introducing domains to the material points, which was first proposed by Bardenhagen and Kober (2004) in Generalized Interpolation Material Point Method. Sadeghirad et al. (2011) later proposed the Convected Particle Domain Interpolation, which featured basis function approximation with deformable parallelogramshaped particle domain (parallelepiped for 3D). The approximation usually introduces insignificant amount of error, although the inaccuracy might turn substantial in some special cases. Nonetheless, CPDI shows superior performance over MPM with static particle domain (Kamojjala et al., 2015). While shear-deformability may cause severe particle domain distortion, such occurrence does not degrade accuracy as element distortion does to FEM. Regardless of particle domain, the weak form of the equation of motion still integrates over the undeformed cells of the background grid.
This paper chooses CPDI MPM with Euler-forward time integration as the discretization method, which implementation follows the description in (Sadeghirad et al., 2011). After initial discretization of the material into particles with square domains, the algorithm of CPDI involves 3 phases in each time step: particle-to-node mapping, nodal calculation, and particle update from nodal solution. In the first phase, CPDI basis function maps particle kinematics and internal variables into nodal values. For example, nodal mass, momentum, and internal force calculations follow: where p = particle number, i = node number, n p = number of particles, t = time, p → t i = nodal momentum, m p = particle mass, v → t p = particle velocity, and ϕ app ip = CPDI-approximated grid basis function described as (see Sadeghirad et al., 2011 for derivation): where x → p 1 to x → p 4 denotes the positions of the 4 corners of the particle domain, and S i = standard grid basis function (i.e. standard MPM basis function).
Mapping other variables such as nodal acceleration, body force, and traction follow similar way (see Sadeghirad et al., 2011;Sulsky et al., 1994 for more details) Nodal solution for the next time step is then calculated in the second phase. For instance, nodal momentum update (see Sadeghirad et al., 2011;Sulsky et al., 1994 for other variables) is calculated as: where f → int i and f → ext i are nodal internal and external forces, and Δt = time-step size.
Finally, particle update from nodal solution takes place in the third phase. The update of e.g. particle position follows: where x → t p = particle position, and n i = number of nodes. Node-toparticle mapping of other variables follow similar procedure depending on the governing equations (see Sadeghirad et al., 2011;Sulsky et al., 1994 for the mapping procedures of other variables).
Afterwards, we go back to the first phase to solve the next time step. Throughout the simulation, particle domain changes shape as parallelogram, following material deformation according to: where p = particle number, r → t 1,p and r → t 1,p are the two vectors defining the parallelogram at the current configuration, r → 0 1,p and r → 0 1,p are those at initial configuration, and F t p is deformation gradient of current configuration. Such approach naturally allows domain elongation, shearing, and rigid-body rotation.
MPM simulation involving wood material has been relatively numerous (e.g. Matsumoto and Nairn, 2012;Nairn, 2015;Nairn, 2016;Aimene and Nairn, 2015). Nonetheless, unlike the mentioned contributions, the simulations presented in this paper utilise CPDI Material Point Method, which is more accurate numerically. The paper reports on the implementation of the wood material model into a custom dynamic explicit CPDI MPM code following the formulation described in (Sadeghirad et al., 2011). Our solver is yet to include implicit time integration, as MPM is most beneficial for fully dynamic cases while current state of the art implementations of MPM with explicit time integration are more robust. Additionaly, our future works require the implementation of the wood material model in conjunction with dynamic fracture. Unfortunately, suitable dynamic validation cases for wood are scarce, therefore our simulations attempt to replicate quasi-static cases instead. The first validation test case involves a novel asymmetric compression experiment of spruce perpendicular to the grain direction. We also demonstrate the capability to reproduce the tensile failure of spruce beam with knots in our second validation test case.

Overview
This section reviews the wood material model by Schmidt and Kaliske (2006). Several minor modifications are introduced into the implementation of the material model, elaborated in Appendices A-C. Nonetheless, the principles stay true to the original model.
The model assumes orthotropic linear elasticity in radial (R), tangential (T), and longitudinal (L) (see Fig. 1), and seven yield criteria: for tensile and compressive failures in the three orthogonal directions and for a failure in shear. A unique hardening function accompanies each yield surface. The hardening function for each tensile failure surface simulates the brittle response (Eberhardsteiner et al., 1999;Eberhardsteiner, 2002), with the total plastic work taken equal to the fracture energy of the relevant loading direction (Aicher, 1994;Dourado et al., 2008). Meanwhile, the hardening function for the compression failure simulates plateauing and densification as observed in experiments (Adalian and Morlier, 2001;Adalian and Morlier, 2002).
Due to the orthotropic nature of wood, Young's moduli, shear moduli, and Poisson's ratios vary for the R, T, and L directions, which are reflected in the elastic constitutive tensor (see Lekhnitskii et al., 1964). For simplicity, this paper utilises Voigt's notation with order the same as that in (Schmidt and Kaliske, 2006). For example, the order of components of stress tensor σ and strain tensor ε are as follows:

Yield criteria
The seven yield criteria are given by a generic equation: where q m is the hardening function for failure mode m; a (m) and b (m) are 2nd and 4th order tensors defining the yield surface for failure mode m, given by Schmidt and Kaliske (2006). In Voight's notation, the 2nd order tensors a m are given by: Fig. 1. Orthotropic directions of wood.

Table 1
Components of b tensor.
The 4th order tensors b (m) are defined as diagonal matrices in Voigt's notation. For 1⩽m⩽7, the components of b (m) are given in Table 1 where f tr , f tt , and f tl are the tensile strengths in the radial, the tangential, and the longitudinal direction, respectibely; f cr , f ct , and f cl are the compressive strengths in the radial, the tangential, and the longitudinal direction, respectively; f vrt , f vtl , and f vrl are the shear strengths in the radial-tangential, the tangential-longitudinal, and the radiallongitudinal direction, respectively.

Tensile hardening functions
For every m, the hardening function q m from Eq. (7) is independently controlled by its own hardening parameter α m , 1⩽m⩽7: The model assumes an associated flow rule: where λ is a plastic multiplier, which is not dimensionless (see Appendix A for the full derivation).
Hardening for the tensile yield criterion simulates a brittle failure. The hardening function is then defined as: where: where l c is the characteristic length, f m stands for the tensile strength associated with surface m, G (f) m corresponds to the fracture energy associated with surface m failure, and κ m is the relative residual strength. The correct value of the relative residual strength is zero, but in the simulation it may be taken as a value close to zero for numerical reasons. The idealization of the tensile failure in Eqs. (12) and (13) is shown in Fig. 2 (see Appendix C for stress-strain curve derivation in up to equation C5).
The characteristic length ℓ c is represented as the element size in the numerical discretization, and cannot be larger than a critical value:

Compression hardening functions
Unlike tension, compression failure exhibits a significant amount of ductility. Experiments by Morlier (2001, 2002) shows a significant plateauing after a compression failure, followed by a rapid increase in the compressive strength due to densification. Fig. 3 shows the typical stress-strain curve of a wood sample under compression in the L and T directions. To take it into account, the hardening function is described as follows: where α m,d and α m,max marks the start and end of the densification phase.
The value of κ m should be selected to be approximately 1 in order to simulate a plateau. For a declining plateau (L compression; m = 6), κ m should be taken less than 1, which results in an idealized stress-strain curve shown in Fig. 4. Otherwise, for an inclined plateau (R and T compressions and shear; m = 2,4,7) κ m is set to be larger than 1, resulting in the curve shown in Fig. 5. The shear failure follows the same hardening function described in Eq. (15) with greater than 1κ m .

Re-derivation of the hardening parameter evolution
The original material model assumes the amount of plastic work throughout uniaxial tensile failure to be equal to the fracture energy: Based on the generalized Eq. (16), we show in Appendix B that the increment of the hardening parameter on the active yield surface is linked with the scalar plastic multiplier by the relationship: where f t is the tensile strength for the appropriate failure surface. In further calculations, Eq. (17) replaces Eq. (11), allowing for a consistently dimensionless hardening parameter α while maintaining the intended original hardening mechanism. Therefore, the model used in this paper follows the hardening rule: Three sets of Eq. (18) exist, one for each tensile failure surface and relevant when said failure surface is active in an elasto-plastic deformation.

Material parameters
The material model requires experimental results to select elastic parameters, material strengths (tension, compression, shear), fracture energy, and compression load-displacement curves, in all orthogonal directions for each wood material. Despite the modifications, the presented material model requires the same number of parameters as the previous material model (Schmidt and Kaliske, 2006), which is 41 parameters. Such a high number of parameters is necessary due to repetitions for all 3 orthogonal directions and 7 failure surfaces, as each surface requires its own hardening parameter. The first 9 parameters determine the orthotropic elastic behaviour, which values are taken as follows (Schmidt and Kaliske, 2009): The next 9 parameters involve material strengths in terms of tension, compression, and shear in each direction (Schmidt and Kaliske, 2006;Miksic et al., 2013): Another 3 parameters determine fracture energy for tensile failures in the 3 orthogonal directions, with physical meaning in the context of the proposed material model explained in subsection 2.3 and further in Appendix B.1.
The remaining 20 parameters determine the densification curves in each direction, matching the stress-strain curves from experiments in T, L and R directions. Subsection 2.4 discusses the physical meanings behind the parameters. The taken parameters in the calculations correspond to those for Norwegian Spruce at 12% moisture content (Mackenzie-Helnwein et al., 2003), aimed to match experimental result in (Adalian and Morlier, 2001;Adalian and Morlier, 2002).
Additionally, wood is a natural material, with significant spatial variability. This study has not implemented any variability of material properties, beyond the change of the grain direction around the knots. Introducing random variation of material parameters and performing large number of Monte-Carlo simulations would likely result in a better representation of real wood behaviour. However, getting data to justify the parameters variation would likely be challenging, hence the introduced parameters random distributions would likely be judged as arbitrary.

MPM implementation
Verification cases presented in this paper utilizes the in-house dynamic-explicit CPDI MPM code built at Aalto University programmed with the proposed material model. For the MPM simulations, total strain increment is defined as (Sulsky et al., 1994): which assumes a small strain increment between the timesteps while representing a finite strain for a large number of timesteps. We assume total strain is split additively into elastic strain and plastic strain: where ε e = elastic strain, and ε p = plastic strain. Naturally, the following relationship applies: where D = the elastic orthotropic stiffness tensor, while σ = stress tensor. The stress increment is a function f of material parameters P , strain increment ε , and hardening parameter α To calculate stress increment σ during elasto-plastic deformation, we use of 2nd order explicit stress integration (derivation in Appendix D), with implicit-explicit hybrid integration for when snapback due to rapid softening occurs (derivation in Appendix C).
In a given time step total strain update bases on the velocity gradient and leads to logarithmic strain over the course of the simulation: where Additionally, the rotation of the material points is calculated via polar decomposition, the same way as e.g. (Biswajit, 2006;Yarahmadi et al., 2013).

Compression test
A compression test of a wood specimen with the grain orientation perpendicular to the compression direction is carried out with a Universal Testing Machine (UTM), which provides a displacementcontrolled compression, along with the output of the acting compressive force and displacement throughout time. The specimens used in the compression test have their growth ring centre located around one edge. Fig. 6a shows the orientation of the test specimen's wood grains. The specimens are situated on a hemispherical knee joint to avoid inducing bending. The knee joint is adjusted and then locked after a preload is applied. The loading rate is set at 0.5 mm/minute, so that failure according to 1% offset method happens at around 300 s (5 min) of loading (Finnish Standard Association, 2012).
Three specimens are tested, shown in Fig. 6b to d. The loading processes of specimens 1 and 2 are video-recorded. The recording device, unfortunately, failed while recording the loading of specimen 3, though the load and displacement data are still recorded by the UTM.
The key results of the experiment with specimen 1 are shown in Fig. 7. At 8 mm vertical displacement, a curved shape is observed. At 12 mm displacement, small shear cracks start to appear at the bottom-right corner of the specimen, getting more visible at 17 mm displacement. At the end of the loading process, a large crack formed on the right side of the specimen.
The experiment is replicated with CPDI MPM under 2D plane strain assumption. For specimen 1, the simulation was carried out using two mesh densities, i.e. spatial resolutions. Fig. 8a, 8b, and 8c show the schematic of the 2D representation and the MPM initial configuration for two mesh densities (the red lines represent the orientation of the material's tangential direction). The particles are spread evenly with 4 particle-per-cell arrangement at the initial configuration. The coarser setup involves 176 particles occupying 44 cells, while the finer mesh density includes 704 particles within 176 cells at the start of the simulation. The wood sample is fixed at the bottom. The steel plate moves  with the velocity of 1 mm/s, with the fixed velocity boundary initially set at the grid nodes at the bottom of the steel plate. The thickness of the steel plate is exactly the same as the target displacement of 20 mm, hence at the end of the simulation, the top of the steel plate coincides with the location of the fixed-velocity boundary. Such arrangement is necessary as the used MPM code is yet to feature displacement control at the particles. Fig. 8b and 8c represent the steel plate as blue circles without a red line and show the positioning of the steel plate in the MPM initial configuration.
The choice of timestep follows the requirement of Courant stability condition for Euler-forward explicit time integration: where Δt = time-step size, Δt cr = critical time-step size, Δx = cell size, and c = the speed of elastic wave propagation calculated as: where E = elastic modulus, and ρ = material density.
The time-step size chosen in the presented simulations are half of the critical value, i.e. Δt = 0.5(Δx/c).
We define the steel material for the steel plate with elastic modulus of E s = 200 GPa, Poisson's ratio of ν = 0.3, and density of ρ s = 7600 kg/m 3 . Such assumption causes the steel material to be the determining critical time-step size. Anyhow, the resulting critical timestep size would be prohibitive to our simulation in terms of computational cost, therefore we introduce mass scaling into our simulation. Artifically increasing material density through mass scaling allows quasi-static simulation to run at higher critical time-step limit. An increase by a factor β increases the critical time-step size by ̅̅ ̅ β √ . Ceccato et al. (2016) implemented the concept in MPM with quasi-static case, which reported insignificant deviation with mass scaling factor up to 400. MPM simulations in e.g. (Martinelli and Galavi, 2021;Mishra et al., 2019) employ the same approach, while other examples besides MPM include the applications in e.g. (Kim et al., 2003;van den Boogaard and Huétink, 2006;Wang et al., 2007). Our simulations in this paper uses mass scaling factor of 100. The arrangement resolves into Δt = 1.479 × 10 −4 s as the chosen time-step size for the denser mesh simulation and Δt = 2.958 × 10 −4 s for the coarser counterpart. Fig. 9a and b display the result of the MPM simulation of specimen 1 with a coarser mesh, and Fig. 10a and b with a denser mesh. A comparison between the resulting load-displacement curves is shown in Fig. 11. Calculations with both mesh densities show similar deformedshape outlines and load-displacement curves. However, as expected, the result with the finer mesh resembles the experimental result better, therefore is taken for further discussion. Additionally, the simulations of specimens 2 and 3 use the finer mesh density only.
At 8 mm of displacement, the outline of the specimen starts to curve similar to the experimental result. The plastic deformations start to  concentrate at the top and the bottom right corners. After further loading, the maximum plastic shear strain in Fig. 10b shows a shearband-like pattern forming from the centre of the annual ring diagonally towards the top and bottom right corner of the specimen which   goes diagonally from the top left towards the bottom right. For a better comparison, the MPM result is overlaid on top of the experimental results in Fig. 12. The numerical and the experimental results are in very good agreement up to around 8 mm of displacement. Some deviation starts to appear at around 12 mm displacement. At 17 mm the deviation of the MPM result is significant, while the large crack that starts to form in the experiment is not replicated. This is likely because the current CPDI-MPM implementation does not include discontinuity or fracture. Similar behaviour is observed in specimen 2, as shown in Fig. 13. Nonetheless, Fig. 14 shows a good qualitative agreement between the load-displacement curves produced by the experiments and the MPM simulations. The simulation aims for qualitative replication, though quantitative comparison shows that the displacement at yield is captured accurately, while the load at yield are underpredicted by 11%

Tensile test with knots
At a structural element level, the presence of imperfections such as knots adds to the complexity of the behaviour of wood material.  and  investigated the effect of the knot clusters in structural elements in tensile tests. The specimens tested were 44x126x4000 mm 3 spruce beams with a knot or a knot cluster, while the observed area was 126x150 mm 2 for each beam. The study concentrated mainly on the strain distribution around the knot, which was acquired via digital image correlation. Several specimens were tensioned up to the failure with crack patterns recorded. MPM simulation attempts to replicate one of the crack patterns in the form of a smeared crack indicated by the accumulation of the plastic strain. The material model uses the same material parameters as those described in subsection 3.1.
Let us consider a tensile test in the y direction. There are three aspects considered in representing the knots in a simulation . First, most of the tensile stress in y-direction flows around the knots, which makes the knots themselves structurally insignificant in transmitting that load. Second, the knots are compressed in the x direction and, in this case, provide significant resistance. Last, the grain direction becomes irregular around the knots, commonly showing a pattern where the wood grains are seemingly going around those knots. Thus, in MPM the knots are represented by an elastic material with elastic modulus E y close to zero and E x equal to 1.5 times E r , while the   simulation setup takes into account the variation of the wood grain orientation around it. Fig. 15 shows the MPM model representing the tensioned beam and illustrates the observed area of the particular specimen to be replicated (see  for further details). There are two oval regions in the clear wood which are replaced by the anisotropic elastic material described earlier to simulate the two knots in the specimen. The orientation of the material points representing the wood material is adjusted to simulate the variation of the grain direction in the clear wood around the knots, with the longitudinal direction indicated in the figures by the lines plotted over the material points. The MPM model does not simulate the whole 4000 mm length of the beam. Instead, two 125 mm zones are added at the top and bottom of the observed area of 150 × 125 mm 2 to minimise the influence of the supports. Other approaches to replicate the distribution of the grain direction in the MPM model are shown in Fig. 16. The grain direction in the immediate area around the top-left knot goes around in a more circular manner than in the second case (Fig. 16a), while in the third case (Fig. 16b) this approximation also applies to the bottom-right knot.
The properties of the steel material attached to the top of the wood specimen, boundary condition, and calculation of time-step size are the same as those in the previous test case (subsection 4.1). The difference in cell size results in time-step size of Δt = 2.4650 × 10 −4 s. The different geometry necessitates adjustment of loading rate to 5 mm/s. The tensile failure path of the experiment is illustrated in Fig. 17. From the left side, the crack goes through the top knot, then down, and through/over the bottom knot before finishing off with a horizontal crack to the right side. Crack as a result of the tensile failure is indicated in the MPM simulation by a band of high plastic strain concentration. This particular pattern is replicated well by the MPM simulation with the 1st approach of the grain direction distribution modelling, shown in Fig. 18a. However, a slight change of the grain direction distribution leads to significantly different crack patterns. The different modelling approaches of what is principally the same grain direction distribution results in very different crack patterns, shown in Fig. 18b and c. Nevertheless, this phenomenon is also observed in the actual experiments. As a natural material, wood features an inherent variability throughout a specimen or element. This variability results in qualitatively-varying possible failure modes, some of which can be realized by experiment replication, while the failure loads remain similar. The MPM simulation replicates this behaviour: three failure patterns are shown (one of which matches that of the original experiment), while the failure loads remain constant at around 100 kN, as shown in Fig. 19.
The actual failure load of the lab specimen is 140 kN, higher than the result from the simulation. Although a qualitative replication is achieved, the quantitative agreement could be perhaps reached with material parameters obtained from laboratory tests of the specific wood used in the test, instead of taken from the literature.

Conclusions and future work
This paper presents simulations of wood failure with a CPDI Material Point Method. First, we show an implementation of a constitutive material model for wood by Schmidt and Kaliske (2006), augmented with a modified relation between the hardening parameter and the corresponding scalar plastic multiplier. The modification solves the dimensional inconsistency of the original model and provides better robustness in the case of severe softening.
We simulated two sets of experiments of Norway Spruce. In the experiment with compression perpendicular to the grain direction, the simulation has shown a good qualitative agreement with the experimental deformation until the point where a large crack opening starts to occur in the experiment. The load-displacement curves are well replicated before the yield, while the yield-loads are underpredicted by 11%. In the simulations of the tensile test, where a specimen with knots was loaded in the direction parallel to the grain, the simulation achieved good qualitative agreement, with the failure load within 29% of the peak load from the experiment.
Quantitative differences between the simulation and the experiment in the two cases most likely arise from the values of material properties utilised in the simulationsthey were taken from literature and are not fully accurate for the woods in the experiments. In future studies, testing of the mechanical properties of the same wood as used in the experiments, as well as incorporating the natural variability of wood material in a probabilistic manner may improve the quantitative agreement between the simulation and experimental results. We also work on adding a capability to model discrete cracks. Additionally, simulation of the timber beams with knots under tension could utilise a contact law between the knots and the clear wood, therefore allowing separation under tension while load transfer under compression remains possible.
Future works within the same project, aimed towards wood cracking simulation at high temperature, will further discuss the aforementioned material property changes due to moisture and temperature, alongside the effect of thermal degradation, i.e. charring. Said degradation poses a particular challenge to the use of timber in construction. Charring produces a protective barrier preventing further degradation, yet the formation of cracks and fissures partially counter this effect. The capability of simulating such complex processes will lead to better understanding and improving the safety of timber structures. Fig. 19. Load-displacement curves of the 1st (a), 2nd (b), and 3rd (c) tensile test modelling approaches.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements
This work was supported by the Academy of Finland (Grant Number 297030).

Appendix A. Dimension Analysis
In the original model formulation, the evolution of state variables, for the given computation point loaded with a certain strain increment, vary depending on the chosen units of the material parameters. A simple dimensional analysis presented in this appendix shows the inconsistency within the formulations of the material model.
To cover the base units, L, M, and T represent dimensions of length, mass, and time respectively. For the derived units, F and E represent force and energy respectively, which are equivalent to: First, we identify the dimension of α m from Eqs. (11), (12), and (13). The dimensions of Eq. (13) are as follows: an indication that η is dimensionless. Since ηα m has to be dimensionless for exp( − ηα m ) in Eq. (12), α m is also dimensionless.
Secondly, we identify the dimension of α m from Eqs. (7), (10), and (11). The dimensions of the variables in Eq. (10) are as follows: Plastic strain ε p is dimensionless, therefore its first derivative against time ε p is in the dimension of T −1 . Also by inspecting Eq. (7), F m is dimensionless. Substituting these back to Eq. A(4) gives us: The relationship between λ and α from Eq. (11) allows the calculation of the dimension of α m : which concludes in α m having the dimension of F⋅L −2 , inconsistent with Eq. A(3).

Appendix B:. Tension Hardening Modification
To solve the dimensional inconsistency (see Appendix A), we revisit the intended amount of energy dissipation for tensile failure.

B.1. Tensile failure energy dissipation
The amount of the dissipated energy per unit area D under axial tension until failure is as follows: where ℓ c is the characteristic length and ε p is the plastic strain.
During plasticity, taking Eq. (7) with F = 0 for the uniaxial tension, we get: Therefore, the change of elastic strain during plasticity is: Expanding dF = 0 results in (see Eq. C(1) for 3D case, or see e.g. (Simo and Hughes, 1998): Substituting Eq. B(3) to equation B (5): Substituting Eqs. B(2) and B(7) to equation B(1) gives: By using integration by parts, equation B(8) becomes: By going on for arbitrary N times, equation B(8) becomes: Many λ(α) functions can fit into equation B12, including polynomial functions of any order. The original model takes 1 st order polynomial (linear) function of λ(α), therefore the 2 nd derivative of λ(α) and higher vanishes from equation B12. Additionally, λ(α) being 1 st order polynomial implies constant 1 st derivative. therefore: From the original model, the amount of dissipated energy is taken as equal to the fracture energy 1 in the same loading direction.
The decay rate η is defined in Eq. (13), thus Eq. (B14) becomes: 1 Note that the model only takes fracture energy as the amount of energy dissipation instead of incorporating Griffith theory of fracture. See e.g. (Zehnder et al., 2013) for the definition of fracture energy and its role in linear elastic fracture mechanics.

B.2. Critical Characteristic Length
The concept of critical characteristic length is to prevent snapback due to sudden softening (Carpinteri, 1989). For this particular case, let us first rearrange equation B (7): To prevent dα/dε from going towards infinity or negative at the most extreme case of α = 0, the following needs to be satisfied: Inequality B21 establishes the critical characteristic length, the same expression as reported in (Carpinteri, 1989;Mackenzie-Helnwein et al., 2001;Schmidt and Kaliske, 2009).

C.1. Explicit stress integration
The previous material model (Schmidt and Kaliske, 2006) recommends the use of implicit stress integration described in (Simo and Hughes, 1998). Due to the independent hardening, an implicit stress integration may converge towards multiple solutions that cannot be proven incorrect. Therefore, this paper proposes implementation with explicit stress integration (unless under special circumstances), which always resolves into a unique solution. The explicit integration employs the Modified Euler algorithm with error control (Sloan, 1987).
The explicit stress integration algorithm follows the algorithm given in (Sloan, 1987), which in principle keeps the stress state on the yield surface by keeping dF = 0 (see dF = 0 expansion in e.g. (Simo and Hughes, 1998): Because of the independent hardening, derivatives on the last term of equation where [D] = elastic constitutive matrix; T = dimensionless time, which is 0 at the beginning of the current elastoplastic timestep and 1 at the end. The method follows the Modified Euler integration with error control (Sloan, 1987) and drift correction (Sołowski and Sloan, 2016). C1. Stress-strain curve with the incorrect Δλ < 0 due to severe softening.

C.2. Severe Softening Special Case
Let us first establish an analytical stress-strain curve during plasticity via parametric function. By integrating equation B(5) and combining equation B(2), we have: Let us then revisit equation C(4). Large enough ∂q/∂α may cause the Δλ < 0, clearly an incorrect result. Fig. C1 compares numerical solution to that from the analytical expression in equation C(5). The analytical path exhibits a snap-back behaviour as observed in (Carpinteri, 1989). In the case of highly brittle material, this problem may not be avoidable. Considering 1 st order explicit stress integration fails in this condition, this paper introduces 2 nd order explicit stress integration, allowing for the correct identification of Δλ > 0 . One-step implicit stress correction then follows the stress integration, guaranteed to converge towards state variables with Δλ > 0 by adjusting the initial guess.

C.2.1. Yield surface second derivative
During severe softening, 2 nd order approximation of Δλ is as follows: The derivation of equation C6 is written in Appendix D.
There are two solutions of Δλ, and the one with positive value is taken. With this Δλ from equation C6, a new stress state {σ n+0.5 } is calculated. Although the equation results in positive Δλ, the resulting drift is considerably large. A one-step implicit stress integration then follows, nudged towards the correct solution of state variables with Δλ > 0 by having more favourable initial guess via λ-iteration.

C.2.2. λ-iteration
An implicit method is utilized for correcting the current stress state {σ n+0.5 }, which is implemented after equation C6. The implicit correction will result in a new stress state and an additional Δλ . Due to the severe softening, the implicit iteration may converge into a solution with Δλ > 0 (the correct solution), or with Δλ < 0 (the incorrect solution).
To assure the convergence towards the correct solution, the algorithm shifts the initial guess so that the initial stress state at the beginning of the implicit solving process is inside the yield surface. This paper utilises a simple iterative procedure to calculate the appropriate Δλ and the stress state for the initial guess, whichfor the ease of referencingwill be called λ-iteration.
In the beginning, we assume that the whole applied strain increment in the current time step {Δε e } contributes towards the plastic strain increment. It is not possible to be exact, as the direction of plastic strain increment is determined by the direction of the yield surface. Nonetheless, the plastic strain multiplier Δλ can be estimated in the following way: With this Δλ, the initial guess of the stress state and the hardening parameter for the implicit solver are: If F > 0 (meaning { σ guess } falls outside the yield surface), Δλ is multiplied by two and Eqs. (C8) -(C10) are revisited with the new Δλ . This process is repeated until Δλ results in F < 0, in which case the most recent { σ guess } and α guess is passed on to the implicit solver as the initial guess. Fig. C2 shows λ-iteration algorithm in a more procedural manner.
Afterwards, a one-step implicit stress integration is carried out with the closest point projection (Simo and Hughes, 1998). Be noted that { σ guess } is not a valid stress state that satisfies the plasticity rules laid out in the previous sections, but rather just an initial guess for the implicit solver. The trial stress for this implicit stress integration is still {σ n+0.5 } .
Following the proposed method allows the material model to work within reasonable stress increment output, useful when the material is highly brittle. However, the energy dissipation may no longer be accurate. If one desires to recover the accurate energy dissipation, keeping mesh size below ℓ c with some safety factor is advised. Fig. C3 shows the resulting stress-strain path when the algorithm is followed.

Appendix D. 2 nd order Taylor approximation of explicit stress integration
In our 2 nd order approximation, the aim is to satisfy: . C2. Flowchart of the λ-iteration algorithm for finding initial guess to avoid convergence towards the wrong solution from the implicit solver.

∂λ ∂T
ΔT + 1 2 Since ΔT is very small, the 1 st term of equation D15 (the one with ΔT 2 multiplier) is relatively much smaller, thus neglegible. Also, the 3 rd term of equation D15 does not have ΔT as multiplier in it, which will cause problem with the error control algorithm within the explicit stress integration procedure (see Sloan, 1987). Therefore, the 3 rd term is also ignored for calculating Δλ at this stage. The resulting error is dealt with by the drift correction scheme employed afterwards. Equation D15 is then reduced into: where Q and R are defined in Eqs. (D17) and (D18) respectively. His main research themes are the numerical fire modelling and simulation, thermal radiation, pyrolysis modelling, and fire risk analyses. Currently, he is the PI for the Academy of Finland research project on the fire safety of wooden materials.