Numerical modelling of wood under combined loading of compressionperpendicular to the grain and rolling shear

Numerical modeling is an efficient tool for experimental validation and for gaining a deeper understanding of complex material phenomena, especially when causal relationships are overlaid by material variability. Wood is such a highly orthotropic and complex material, which in engineering problems however is considered as macro-homogeneous. The aim of this study is to numerically investigate stress and strain states of wood in the radial-tangential plane and the influence of the orthotropic material behavior on the structural response. Model validation is based on experiments performed on clear wood of Norway spruce ( Picea abies ) by using a biaxial test setup. Three material models were used, namely Hill ’ s plasticity model, the Hoffman criterion and a novel quadratic multi-surface (QMS) criterion. After validation on the local material scale, the models were applied to the engineering problem of compression perpendicular to the grain for studying the effect of the unloaded length. As a novel part, the influence of the annual ring structure on the local material behavior and the global elasto-plastic force – displacement behavior of wood under compression perpendicular to the grain were numerically investigated. Hill ’ s failure criterion was found to be the least suitable at both length scales, local material behavior and global structural response. The Hoffman and the QMS criteria showed quite good agreement with the biaxial experiments in terms of force – displacement relations and strain distributions for different loading situations, especially for combinations with radial compression, while there was less agreement with experiments for the behavior of combinations with tangential compression. Application of these material models to compression perpendicular to the grain for studying the unloaded length effect yielded similar trends as observed in structural tests. A reasonable and similar force – displacement response by Hoffman and QMS criteria was observed, while Hill ’ s model yielded significantly overestimated force carrying capacity. Differences in force-– displacement response for different loading situations were well in line with literature findings and the influence of the annual ring curvature on the overall force – displacement behavior could be quantified.


Introduction
The natural cellular and fibrous material wood can be considered as a cylindrically or polar orthotropic material due to its cylindrically shaped, circular annual ring structure and differences in material properties in three principal directions: longitudinal, radial and tangential.The orientation dependence of mechanical properties of wood thus originates from the micro-structure of the material and shows spatial fluctuations at different length scales [1].In engineering applications, however, wood is considered as a macro-homogeneous continuum and, as regards its application in timber structures, transversely isotropic material behavior is commonly assumed.The weak cross-grain directions have great importance in engineering applications and must be carefully considered in the design and detailing of timber structures [2].Considering differences in the two transverse material directions, the radial and tangential directions, are important for the understanding of corresponding phenomena and for improving engineering design rules.Due to the differences in stiffness and strength in the radial and tangential directions and the curvature of annual rings, a combined and complex stress state arises already under uniaxial compression or tension and becomes even more complex in case of loading under an angle to the principal material directions.This calls for a multi-axial stress analysis.
The aim of this study is to numerically investigate the material behavior of wood in the radial-tangential plane, to validate a material model by comparison with test data, and to apply the model to the engineering problem of partly compressed wood perpendicular to the grain.Model validation rests on comparison of the experimentally observed mechanical behavior of Norway spruce (Picea abies) clear wood under compression perpendicular to the grain, rolling shear and combined compression perpendicular to the grain with rolling shear [3], and on numerical modeling considering elasto-plastic material behavior.As a novel part of the presented work, the influence of the annual ring structure on the global elasto-plastic force-displacement behavior of wood under compression perpendicular to the grain will be numerically investigated.This effect is difficult to assess in experiments since it is overlaid by the variability in material properties.Experimental findings that form the basic data-set for the phenomenon to be modeled herein are discussed in the following, in combination with previous experimental and numerical studies available in the literature.
Several research works [4][5][6] investigated variations of elastic stiffness and strength in radial and tangential directions.Hoffmeyer et al. [5] conducted experimental tests on five different types of Norway Spruce structural timber specimens and studied the dependence of stiffness on the annual ring orientation with respect to the applied loading direction.The transverse stiffness was found lowest when the annual ring structure was orientated at 45 • to the loading direction.The radial stiffness was found to be about two times higher than the tangential stiffness, which is in agreement with Hall [7], Madsen et al. [8] and Siimes and Liiri [9].A recent study by Akter and Bader [3] on Norway spruce clear wood found radial stiffness to be about 1.5 times higher than tangential stiffness, in case of both tensile or compressive loading.Farruggia and Perré [4] reported higher anisotropy ratios of between 3 to 4 from tensile tests on Norway spruce juvenile earlywood at the microscopic level.These results indicate lower anisotropy ratios for higher wood densities.Regarding strength in the transverse crosssectional plane of wood, Hoffmeyer et al. [5] and Madsen et al. [8] concluded that strength did not show much dependency on annual ring orientation, which however is in contradiction to Hall [7], Kennedy [10] (cited in Brandner [11]) and Akter and Bader [3], who found radial compressive strength to be higher than tangential strength.The latter experimental work presented a biaxial test setup that allowed to quantify stiffness and strength for stress combinations in the radial-tangential plane of wood.A positive effect of compressive stresses on the rolling shear strength was found and the test data was compared to different failure criteria.Tests were conducted for the combination of radial compressive and tangential compressive stresses with rolling shear stress.
Regarding timber engineering design situations, the effect of boundary conditions (the so-called effect of the unloaded length) in compression of wood perpendicular to the grain is of great importance.It describes the apparent increased strength of a partly compressed compared to a fully compressed wooden element perpendicular to the grain.A huge variation of design strengths became obvious in experimental tests [12][13][14][15].In engineering design standards, boundary effects are considered by the so-called k c,90 factor (with values from 1 to 1.75, see e.g.EN 1995-1-1 [16] and Blass and Görlacher [13]).Brandner [11] and Leijten [17] recently reviewed design models for spruce loaded in compression perpendicular to the grain and discussed stress distribution for the above-mentioned effect in glued-laminated and cross-laminated engineered wood-based products.
In many cases, the experimentally determined material responses are compared with numerical models for validation purposes and for gaining better insight into this complex material behavior of wood.However, numerical modeling of the highly orthotropic, hierarchically organized and porous material wood is a challenging task for the engineering community, especially regarding its highly ductile behavior under compression loading and brittleness under shear and tensile loading.Therefore, numerical models of wood are mostly based on specific approaches depending on the problem area.Commonly, elastic models [5,6,18] in plane-stress conditions are in use while comparing with experimental elastic behavior and stress distribution at low load level.Brittle behavior can be modeled by adopting fracture mechanics or adding a cohesive zone model [19,20].Classical flow theory of plasticity considering mainly Hill [21], Hoffman [22], or Tsai-Wu [23] failure criteria can be used to represent the ductile behavior of wood under compression.
Shipsha and Berglund [6] conducted a comparative study of strain distributions in wooden boards of Norway Spruce, subjected to transverse compression.Two types of boards were investigated, namely one from the pith area and another one from the periphery region of the stem.The size of the specimens was 30 × 20 × 15 mm 3 (l × h × b).Loading of the specimens was kept within the elastic limit and very small strains.Therefore, elastic material behavior could be considered in the numerical study.Polar orthotropy was considered for defining local material orientations in the model by prescribing the distance from the pith and the position of the specimen in the board.Radii of 0 and 60 mm represented the center and periphery specimen, respectively.The numerical results for the center specimen identified a very non-uniform strain distribution (with a factor of three for differences in local strain values) in comparison to the periphery specimen, which showed a more homogeneous strain distribution, as an effect of the smaller curvature of the annual ring structure.
Hoffmeyer et al. [5] and Wallin [24] reported pre-mature failure of glued laminated timber (GLT) due to the development of tensile stresses perpendicular to the grain when the specimen was under compressive force.Differences in radial and tangential elastic stiffness as well as annual ring curvature caused such failure.For validation purposes and for gaining better insight into stress distributions, an elastic numerical study was performed on solid wood and GLT.Material anisotropy in radial and tangential directions and the annual ring curvature effect were considered by assigning the location of the pith for each individual lamella in GLT.Severe non-homogeneous stress distributions were found from numerical models for both structural timber (tangential direction) and GLT specimens, which is in agreement with the study by Aicher et al. [18].High tangential tensile stress and shear stress were observed in GLT and solid timber under radial compression.This confirmed the reason for the early failure of GLT before reaching its compressive strength.
Nonlinear material behavior was not only observed under compressive loading of wood but was also found for shear loading.Dahl and Malo [25] investigated corresponding nonlinearities for spruce wood loaded in the three orthotropic shear planes of the material and compared experimental data to results from Finite Element modeling.Material parameters for linear, bi-linear and Voce-type (exponential formulation) modeling were determined from experiments and an elasto-plastic model based on Hill's criterion and a bi-linear stress-strain relationship was developed.The latter yielded about 5-30% higher values of rolling shear modulus than the experimentally determined values.
Bogensperger et al. [26] conducted a numerical analysis of solid timber blocks and cross-laminated timber (CLT) cubes.Numerical results were compared with corresponding experiments done by Augustin et al. [27] and Salzmann [28] with the aim of calculating the so-called k c,90 factor, a factor used in EN 1995-1-1 [16] to account for various stress distribution effects in compression perpendicular to the grain.A 3-D finite element model was used and cylindrical orthotropic elasticity was assumed when assigning local material directions to the wood.Material plasticity was considered in radial direction only with isotropic hardening to define the evolution of plastic strains by assuming linear hardening and small strain theory.The model correlated well with experimental results as regards the influence of parameters for the compression perpendicular to the grain strength and corresponding k c,90 values.
Alternative modeling approaches that are of relevance for the abovedescribed design situations are models based on continuum damage S.T. Akter et al. mechanics [29,30], the material point method [31] and hybrid models [32][33][34].In continuum damage mechanics, a modified stiffness matrix taking into account material damage is used, where the stiffness parameters are reduced after failure initiation.A hybrid foundation model was developed by Toussaint [34] using elastic-plastic orthotropic beam elements in combination with a crushable foam model for the coupling of beam elements oriented parallel to the grain.This model allows capturing both plasticity as related to the crushing of wood cells and brittle failure in tension perpendicular to the grain and shear.The interaction of stresses perpendicular to the grain and shear is crucial in this case.Another type of hybrid model was applied to study compression perpendicular to the grain in timber structures by Persson [33], where an overlay technique was used.The approach included combining a standard elastic orthotropic material model and an isotropic plastic crushable foam model by overlaying geometrically coinciding elements.
The macroscopic phenomenon of compression perpendicular to the grain is a consequence of the microstructure of the material, which however is only implicitly considered in continuum models for clear wood.Several micro-mechanical models for wood are available in literature [1,[35][36][37] and have led to a better understanding of the deformation behavior of the cellular material.These models are based on representative volume elements at different microscopic levels, often in a hierarchical manner at several length scales of the material.
Multi-surface plasticity models [38,39] have been developed for the prediction of the non-linear and plastic behavior of wood at different length scale.Li et al. [36] developed a multi-surface plasticity model for wood by implementing nonlinear material behavior in the micromechanical models of two unit cells, namely for earlywood and latewood.Periodic boundary conditions represented the repetition of the cell walls.The model was validated by comparison of the combination of multiple model predicted failure surfaces with experimental data from biaxial tests on Norway spruce clear wood, reported by Eberhardsteiner [40].The promising outcome of the comparison study showed the effectiveness of the multi-scale plasticity model.However, the multiscale plasticity model is not easy to follow and not readily available for engineers.
Herein, a comparatively simpler material model is verified, calibrated and applied.The model can represent the overall mechanical behavior of wood in the transverse plane for all loading cases except tensile loading, such that the model can be used in application examples for, e.g., solid timber, GLT or CLT to determine the local stress distributions and predict the structural response under compressive loading.The calibration of the model is based on an experimental work on defect free, dog-bone shaped, clear wood specimens of Norway spruce (Picea abies) loaded in tension or compression in the R-or T-direction, loaded in rolling shear, and loaded by tensile and compressive forces in the R-or T-direction in combination with rolling shear [3].In the experiments, loading was applied until failure or up to approximately 30% nominal strain, whereas the model was used for below 20% nominal strain.The following sections of this paper present the considered failure criteria with a description of the corresponding FE-modeling (Section 2), results in terms of comparisons between numerical modeling and experiments for material testing and the effect of the unloaded length (Section 3), and finally, a general discussion and main conclusions drawn from this work (Section 4).

Elasticity
In this work, an orthotropic elastic wood material behavior with three symmetry planes, namely longitudinal (L), radial (R) and tangential (T), is assumed.The stresses in the elastic material can be calculated from the stress-strain relationship in the form of the generalized Hooke's law, which reads σ = D ε, (1) where, D is the elasticity tensor in matrix notation, which is assumed to hold nine independent elastic material parameters, namely three Young's moduli E L , E R and E T , three shear moduli G LR , G LT and G RT , as well as three independent Poisson's ratios ν LR .ν LT and ν RT .The inverse of the elasticity tensor is the compliance tensor C.This allows to express the Hooke's law in matrix notation as ε = C σ. (2) In Voigt notation, it can be written as Due to the fact that C is symmetric, the following relations apply so that only nine independent elastic material parameters remain.

Strain decomposition
The theory presented herein is formulated with small strain theory, and assuming that the total strain ε is the sum of the elastic ε e and the plastic strain ε p .In rate form, this relationship reads as from which the plastic strain rate is given as where a dot denotes the time rate, which is commonly used in theory of plasticity.Using Eq. ( 2), plastic strain can be expressed in finite increment form as

Yield criteria
In this work different failure criteria will be used.Thus, stress based failure criteria according to are adopted where f is the failure criterion, which includes material parameters related to strength.Orthotropic material behavior of wood with different strength values in tension and compression is considered in this work.The three different failure criteria considered herein, namely Hill's, Hoffman's, and the quadratic multi-surface (QMS) failure criterion, are illustrated in Fig. 1 and discussed in the following.

Hill's criterion.
Hill's orthotropic yield criterion [21] is a generalized formulation of the classical von Mises criterion for anisotropic materials, which assumes that hydrostatic stress does not contribute to yielding.The criterion can be expressed in matrix notation as where P is a fourth-order tensor expressed as a six by six matrix holding the material parameters and assumes that yielding only depends on the deviatoric part of the stress tensor, where e is the hydrostatic part of the stress tensor.Eq. ( 9) can thus be rewritten using the deviatoric part of the stress tensor as This reduces the number of independent material parameters in Hill's criterion to six, which leads to the following definition of the P matrix which yields Hill's criterion The material paramters F, G, H, L, M, N characterize the material orthotropy.These material parameters can be determined from six uniaxial tests in tension or compression and in shear.Denoting uniaxial normal strength f i and shear strength f ij and relating these to the three principal material directions (planes) of the orthotropic wood material, the parameters of the material strength matrix, P, are defined as In component form, Hill's criterion reads as where f R , f T and f L are the three normal strength values and f RT , f RL and f TL are the three shear strength values.This failure criterion is not capable of accounting for any possible difference in tensile and compressive strength.The criterion is often implemented in commercial Finite Element codes and might be easy to apply.It should however be noted that it is limited to a certain degree of anisotropy.This is a consequence of the assumption that yielding does not depend on the hydrostatic part of the stress tensor.For large degrees of anisotropy, this could lead to unrealistic stress states in wood, since the failure criterion may constitute an open failure surface in the deviatoric stress space [41][42][43].In order to avoid this, certain inequality conditions need to be satisfied [43].This is however not of relevance for the herein considered interaction of normal stresses perpendicular to the grain and rolling shear stress.

Hoffman's criterion.
The Hoffman [22] failure criterion can be interpreted as a further extension of the Hill criterion.It considers the difference in tensile and compressive strength in three principal directions, which makes this criterion more relevant to wood.The criterion, furthermore, accounts for the effect of hydrostatic stress in yielding.Eq. ( 9) is extended by a linear term and a second-order tensor q.In matrix-vector notation, the criterion then reads as The components of the material strength matrix P using the notation of Eq. ( 12) are given as Akter et al. where, the indices t and c denote tensile and compressive strength respectively in the corresponding directions.q in Eq. ( 21) cannot involve components that affect shear stresses, and consequently q T = [q R q T q L 0 0 0] In component form the Hoffman criterion is defined as It should be noted that the off-diagonal terms of P in Eq. ( 12) were neglected in the implementation of the Hoffman criterion in this work.This means that the failure surface, the six-dimensional ellipsoid, is oriented with its principal axes along the main axes of the normal stress expressed in the main material directions.Therefore, the material matrix P reads as 2.2.2.3.Quadratic multi-surface (QMS) failure criterion.The quadratic multi-surface (QMS) failure criterion is based on a quadratic interaction criterion without coupling terms between the individual stress components.It is defined as a multi-surface criterion for each combination of stress components in the stress space with the failure criterion.This criterion has several advantages compared to the previously discussed single surface failure criteria such as its applicability to any extent of material anisotropy since the failure criterion is composed of several ellipsoid surfaces.The inequality conditions required for the components of P in Eq. ( 12) for ensuring a closed convex surface of Hill's and Hoffman's criteria do not apply here.Another advantage relates to the predicted development of plastic strain at uniaxial loading when associated plasticity is assumed.For a failure criterion that includes offdiagonal terms (see Eq. ( 12)), uniaxial loading then results in a Poisson effect that might be non-physical, or at least clearly in contrast to the expected behavior [44].The formulation of the criterion reads as The material matrix P includes material parameters and is defined as follows, The QMS criterion allows to take into account the difference in tensile and compressive stress by setting the normal strength values f R , f T and f L in Eq. ( 34) equal to f R,t , f T,t and f L,t or f R,c , f T,c and f L,c , depending on whether the corresponding normal stresses σ RR , σ TT and σ LL are greater or smaller than zero.The ellipsoid defined by Eq. ( 33) is thus constructed by different sub-surfaces, corresponding to the various combinations of positive and negative values of the normal stresses, and its design assures continuity of the failure surface and its normal.In component form, the criterion can be written as

Plastic flow rule and strain hardening
Associated plasticity was assumed and thus where εp is the plastic strain rate and λ is the plastic multiplier.In case of hardening of the material, the yield function is re-written as where σ y is the current normalized yield strength, which is equal to 1 in case of perfect plasticity.Hardening was defined by a piecewise exponential function using λ as an internal variable.It was assumed that the hardening could be described accurately enough with two such functions, such that after some plastic straining, a densification could be simulated (see Serrano [44]).Consequently, σ y was defined as The reason for this choice is that the σ-ε p relation will become piecewise linear, for the QMS failure criterion introduced here.The abovementioned Hoffman and QMS-criteria were implemented as user subroutines for the commercial FE-software ABAQUS [45].In that implementation, the approach as described in Schellekens and De Borst [46] for integration of elasto-plastic equations and calculation of the elastoplastic tangent stiffness was followed at large.See also Serrano [44].The use of Hill's criterion was in this work based on the built-in functionality of the software.Hardening was for that criterion defined so as to correspond to Eq. (38), by specifying yield stress and plastic strain values in tabular format.

Modeling of uniaxial and biaxial loading cases in the RT-plane
The commercial finite element software Abaqus [45] was used for numerical modeling of the setup and specimens used in the experimental study by Akter and Bader [3].In that study, experiments were carried out for pure compression, pure shear and multiple combinations of compressive or tensile normal stress with rolling shear stress.Dog-bone shaped specimens were used in the tests, and these specimens were cut such that the annual ring orientation differed between the two series, namely for radial compression (R) and tangential compression (T), which were also used for rolling shear loading and corresponding stress combinations.Outer dimensions of the specimen were 50 × 60 × 20  The failure criteria discussed in Subsection 2.2.2 were applied as yield criteria in elasto-plastic material modeling for studying stress interactions and global force-displacement behavior of the abovementioned material test-setup for the R-T plane of wood.Numerical modeling was performed for loading corresponding to pure compression, pure shear and corresponding to three different combinations of compression and rolling shear loading namely SC-30, SC-45 and SC-60, defining loading paths with a combination of compression and rolling shear loading and with displacement ratios of horizontal to vertical displacement of 0.577, 1.00 and 1.73, respectively.
The 3-D geometry was modeled, but assuming symmetry in the thickness direction (longitudinal direction of the wood) in order to reduce the computational cost, see Fig. 3.In a first step, a Cartesian coordinate system was used to define the local material orientation, as the curvature of the annual rings was limited in the specimens due to the position of the specimens in the outer part of the log.
To accurately model the fixing to the testing machine by the load introduction system, the boundary conditions at the top and bottom parts of the specimen were applied not only to the surface of the specimen but also to the interior parts, until a depth of 2 mm.This was meant to represent the mechanical grips of the test setup, see Fig. 3. Numerical models under compression, shear and biaxial loading were analysed for total displacements of 5 mm, 2 mm and 5 mm, respectively.
The elastic material properties and strength values used in the FE modeling are given in Tables 1 and 2. These were determined from the experimental study presented in Akter and Bader [3].The experimental work was limited to material parameters in the R-T plane.The remaining material parameters in Tables 1 and 2 were assumed based on literature data.It should be noted that the latter parameters insignificantly affect the results presented in this work.Hill's failure criterion differs significantly in terms of failure stress, since the difference in tensile and compressive strength was not taken into account and the tensile strength was considered to be equal to the compressive strength.The differences between the Hoffman criterion and the QMS criterion are small, in terms of failure stress.Note that the normal of the failure surface differs.Hardening was defined as in Eq. ( 37), where material orthotropy was considered.Linear hardening was defined by the exponential functions in Eqs.(38) and (39), where the parameters k 1 and k 2 were set equal to 1.0 and 1.5, respectively, and λ 1 was set equal to 0.4.
Mesh dependency was checked for one model of the test setup loaded in radial compression by using linear and quadratic elements for a mesh size variation of 0.5 mm to 5 mm and the Hoffman yield criterion, see Table 3.The specimen was loaded up to a displacement of 5 mm and the variation of the corresponding vertical reaction force was studied.Variations in maximum force were mainly found for linear elements, while quadratic elements yielded insignificant differences for different mesh densities.Coarser and linear mesh resulted in higher reaction force, which decreased with decreasing element size and higher order elements.Based on the mesh dependency study, 8-node linear brick elements (C3D8) of 1 mm size were used, which was found to give a reasonable balance between accuracy, spatial resolution and computational cost.
The influence of the annual ring curvature on the results of the modeling was investigated by comparing the results of the Hoffman model for compressive and shear loading, by changing the coordinate system from Cartesian to cylindrical and defining the pith location.Three different pith locations, namely Cyl-1, Cyl-2, and Cyl-3 (see Fig. 2) were considered for the investigation of the annual ring curvature effect on stiffness and load carrying capacity.For specimens under radial compression, Cyl-1 and Cyl-2 define pith locations at the center position of the bottom edge, 50 mm below the bottom edge center position respectively.For specimens under tangential compression Cyl-1 and Cyl-2 define pith locations at 30 mm to the right side from the center of the specimen, at 80 mm to the right side from the center of the specimen, respectively.Cyl-3 defines the pith location at 50 mm below and 50 mm to the right side of the bottom right corner, which results in annual rings orientated at 45 • with the direction of compression and shear force for both types of specimens for the above mentioned loading cases.

Numerical study of the unloaded length
The material models described in Subsection 2.2.2 and assessed with material testing as described in Subsection 2.3, were used in a FE-model to predict the force-displacement behavior in compression perpendicular to the grain loading of Norway spruce solid wood for fully and partially loaded situations.The latter was used to study the influence of the so-called unloaded length in structural timber.Dimensions of the specimens were chosen according to the experimental study performed in Kathem et al. [15].A qualitative comparison of model predictions with experimental data from that study will be done, in order to assess the model's abilities to predict the influence of the unloaded length.Moreover, the model allowed to assess the influence of the annual ring pattern on the overall force-displacement response.
The cross-sectional size of the structural wooden element was 95 ×  45 mm 2 , while the length in the longitudinal direction was 500 mm for the partially loaded case and 90 mm for the fully loaded case, respectively, see Fig. 4. Load was applied either over the entire length of the specimen or at the mid-or at the edge position of the specimen with a steel bar measuring 90 × 40 × 95 mm 3 (l × h × w).The models were denoted as Mid and Edge for partially loaded cases and as Full for the fully loaded compression test.The 3-D geometry was modeled and constrained displacement boundary conditions on the lower edge of the wooden element were assumed.Loading was applied through a prescribed vertical, but constrained horizontal displacement of the steel bar.8-node linear brick elements (C3D8) of 1 mm size were used, which was found to be a reasonable balance between accuracy, spatial resolution and computational cost.
The influence of the annual ring curvature on the load carrying capacity and stiffness in all afore-mentioned loading cases was investigated by changing the location of the pith.Three different pith locations, similar to Cyl-1, Cyl-2, and Cyl-3 were considered, see Fig. 2.Only Hoffman's yield criterion was used for this numerical study.

Modeling of uniaxial and biaxial loading cases in the RT-plane
The outcome of the numerical modeling was compared with experiments for material testing in terms of force-displacement curves.Additionally, shear stress-normal stress relationships are given in order to compare the experimentally determined stress interaction with the numerical predictions.Results are illustrated in Figs.5-9.It should be noted that the stiffness parameters that were used in modeling, see Table 1, were calculated from the unloading path of the experimental tests by Akter and Bader [3].
Material behavior under uniaxial compression and shear in the two different loading directions, radial and tangential, is discussed first.As regards radial compression and shear in the RT-plane, Fig. 5 shows that all the numerical models predict stiffness properties very close to the unloading stiffness, since material parameters were evaluated from       these parts of the experimental curve.The loading stiffness however is slightly overestimated, as well as the non-linear transition from the elastic to the plastic region.Similar observations are made for the compressive behavior in the tangential direction, see Fig. 6.The shear stiffness was slightly lower in the T-specimens as compared to the Rspecimens, which is however not considered in the model, where an average value was used.The overestimation of the loading stiffness is a consequence of the determination of material properties from the unloading path, which is considered to represent the elastic material behavior.The experimentally determined loading stiffness however is usually lower due to nonlinearities because of the load transfer between the test equipment and specimens, which is not considered herein.Similar observations and modeling approaches for the corresponding nonlinearities were shown in research on timber joints, e.g., [47][48][49].
Strength was experimentally determined as the stress corresponding to 1% nominal compressive strain.In combination with the above discussed stiffness parameters, this does obviously not lead to the same stress-strain point, i.e., the transition point from linear elastic to the plastic behavior, in numerical modeling, see Figs. 5(a) and 6(a).The overall experimentally observed behavior is however well captured for some of the models and loading cases.Hoffman's and the QMS models yielded very similar load-displacement behavior as observed in experiments for compressive loading in the radial direction.The hardening behavior used in combination with the Hoffman and QMS failure criteria gave good agreement with the experiments.Impressively, the model with the QMS yield criterion followed the experimental hardening slope very closely for this loading case.Contrarily, using Hill's yield criterion resulted in a considerably larger compressive force compared to experiments, especially at the elastic to plastic transition zone.This indicates a different development of plastic strains after the yield limit.
The behavior under tangential compression is in general less accurate in terms of predicting the experimentally recorded behavior.The experiments yielded a stress peak from which the stress decreased up to a certain strain level before it started to increase again [3].It should be mentioned that the experimental investigation in [3] showed that the tangential compressive strength was dependent on specimen shape and height, and thus, may not reflect the real material strength.This complex response is difficult to model with a rather limited set of hardening parameters (k 1 , k 2 and λ 1 ) in the material model.The hardening parameters used were aimed at mimicking the behavior of radial compression.
A perfectly brittle or progressively brittle failure, depending on the orientation of the annual rings with respect to the shear loading direction, was observed under shear loading.For the T-specimens this resulted in a more non-linear behavior and a slightly higher strength than for the R-specimens, see Figs. 5(b) and 6(b).This difference is however not accounted for in the considered orthotropic material behavior.Moreover, brittle failure is not included in the elasto-plastic Next, the behavior of the test specimens under biaxial loading is discussed.The comparison of simulations with experiments is limited to the interaction of compressive loading with shear loading, since tensile loading and combination with shear loading yielded brittle material failure not accounted for in the material model.A very clear difference between the material models became obvious.
Nominal engineering stresses and strains were calculated from the numerical models to compare with experimental data for the different stress states.Nominal normal stresses (σ RR , σ TT ) were calculated by dividing the vertical force obtained from the models by the initial minimum cross-sectional area of the specimens, which was 300 mm 2 .
Nominal normal strains (ε RR , ε TT ) were obtained from vertical displacement divided by the initial unsupported height of the specimen.Horizontal/shear force from models was divided by the above mentioned minimum cross-section, at the center of the specimens, to calculate the shear stress, τ RT , while shear strain, γ RT was determined by the horizontally applied displacement divided by the unsupported height of the specimen.The initial unsupported height of the specimen was 30 mm, excluding the supported height of the gripped plate.Fig. 7 shows radial compressive and shear forces as a function of the corresponding displacement, Fig. 8 shows tangential compressive and shear forces that developed for different prescribed biaxial displacement paths, and Fig. 9 shows a combination of these two force directions in terms of nominal compressive and shear stresses in a stress plot.For the R-specimen (Fig. 7), Hoffman and QMS models were well in line with experiments, especially for the SC-60 and the SC-45 paths, except for some deviation at the later stage of the post-elastic stress evolution direction.This is also seen in the stress plot (Fig. 9(a)), where the Hoffman and QMS models are rather close to the experimental curves.For SC-30 paths, however, these two models yielded higher compressive force and lower shear force as compared to the experiments, which lead to a different stress path in the stress plot.Results calculated with Hill's model did not follow the experimental curves, but resulted in considerably higher compressive forces and lower shear forces, including even a shear softening behavior, for all interaction paths.Thus, the stress paths in Fig. 9(a) show a different evolution of plastic behavior.The latter is a consequence of a stress re-distribution which is different to the Hoffman and QMS models.Note that small strain plasticity was considered in all of the numerical models.However, the models were compared with the experiments at comparably high strain levels up to 16%.The deviation in the evolution of plasticity for combined loading cases at higher strain levels as noticed in FE models could be an effect of considering small strain plasticity.
The experimental material behavior as regards tangential compressive and shear forces for different prescribed biaxial displacement paths and corresponding simulation results are shown in Figs.8(b) and 9(b).Since the material models already gave differences for uniaxial tangential compression, simulation results showed big differences compared to the experimental force-displacement curves.This is also visible in the stress plot, see Fig. 9(b) for nominal stress interaction.Similar behavior as under radial compression was predicted by the model.The Hill yield model resulted in considerably higher compressive forces and a shear force softening, which led to stress paths that considerably differed to the ones predicted by the Hoffman and the QMS yield model.In the work presented here, to keep the model simple, the same hardening behavior for both radial and tangential directions was considered which, as seen from comparing the numerical results with the test results, seems to be a limitation.It is difficult to compare the simulations with the experiments, but it appears that the evolution of the plastic behavior is closer to the Hill model predicted re-distribution to compressive stresses, rather than a re-distribution and increase of shear stresses as predicted by Hoffman and QMS yield models.
The experimental test data allows also to compare the strain fields at the surface of the specimens, which were measured by means of an optical digital image correlation system, with model results.The comparison is shown for uniaxial compression in radial and tangential direction, shear loading in the RT-plane, as well as for a SC-45 biaxial loading state with radial compression and for a SC-30 biaxial loading state with tangential compression in Figs. 10 and 11.Corresponding strain fields were calculated and compared for the Hill, Hoffman and QMS material models.Strain distributions in the experiments were found to depend on local material situations i.e., the strength of earlywood and latewood layers.The assumption of a homogeneous material, as was used herein, was not capable of capturing the differences in strains in earlywood and latewood layers.However at higher stress levels, as the earlywood layers were compressed (collapsed), strain concentrations become less pronounced and the difference between earlywood and latewood could possibly be disregarded in a homogeneous approach.
Under shear loading, a comparatively wider spread strain field was observed in the experiments, Figs.10(b) and 11(b) even though failure occurred in the notched part of the specimen.The numerical models yielded a uniform band in the notched part with stress concentration in the lower part of the notch.The models for combined loading develop strain concentrations at the upper left and lower right corners of the notch, which can be explained in terms of the superposition of compressive strain in the vertical direction due to the horizontal and vertical displacement of the specimen.In the experimentally measured strains, this was not observed.
Differences in strain fields between the applied yield models became obvious.Under compressive loading, Hoffman's model yielded the most uniform strain field in both R-and T-specimens.Hill's model yielded similar strain concentrations as observed in experiments for tangential compression, see Akter and Bader [3].Under biaxial loading, Hill's model resulted in high edge concentrations in comparison to Hoffman's and QMS models.
Finally, the influence of the annual ring curvature on stiffness and strength in the RT-plane was studied by means of the simulations.Force-displacement relationships using three different cylindrical coordinate systems Cyl-1 to Cyl-3, with different pith locations were compared with experiments and simulation results with a Cartesian coordinate system.Figs. 12 and 13 show the results from specimens under radial and tangential compression, respectively.The general observation was that elastic stiffness and compressive load carrying capacity decrease with decreasing pith distance from the center of the  specimen, whereas the behavior was vice versa for shear testing.From Figs. 12 and 13, it is clearly visible that both stiffness and load carrying capacity for specimens under radial compression showed higher dependency on pith locations than specimens under tangential compression.Under compressive loading in the radial direction, Cyl-2 i.e., when the pith was at 50 mm below from the bottom edge center position, yielded similar force-displacement relations as the model with Cartesian coordinates.The pith was located far enough to produce almost straight annual rings.Cyl-1 position, where the pith was at the bottom edge center position, yielded slightly lower elastic stiffness and compressive force than the Cartesian coordinate system.Cyl-3, i.e., with annual rings orientated at 45 • , yielded the lowest elastic stiffness and compressive force, which is in agreement with previous findings, see e.g.Hoffmeyer et al. [5].For shear loading on R-T-specimens, the load carrying capacity increased when the pith was located closer to the specimen center.The highest shear stiffness and capacity were found for Cyl-3.This result is in agreement with previous findings, e.g., by Aicher and Dill-Langer [50] and Ehrhart and Brandner [51].Similar to the behavior in compressive testing, Cartesian and Cyl-2 models yielded almost the same shear stiffness and strength, when the pith was located far away enough so that there was no influence on the behavior of the annual ring curvature.Specimens under tangential compression showed a very small difference in elastic stiffness and compressive strength between Cartesian and cylindrical coordinate systems for different pith locations.A similar but smaller increasing trend as for testing in the radial direction was found for shear stiffness and strength with close pith locations to the center of the specimen, see Fig. 13.

Numerical study of the unloaded length
The comparatively simple elasto-plastic material modeling approach was motivated by the aim of its application to the engineering problem of compression perpendicular to the grain (CPG).Figs. 14 and 15 show simulation results in terms of force-displacement relationships for three different loading cases (fully compressed and partially compressed at mid and edge loading positions) in order to quantify the effect of the unloaded length.The length effect in CPG strength is clearly predicted by the models.Model results in Fig. 14 show the force-displacement relationships for radial compression and annual ring pattern Cyl-2.Central or mid position of the load introduction bar (Mid) resulted in higher load carrying capacity than the edge position (Edge), which in turn yielded higher load carrying capacity than the fully compressed (Full) wooden member.The reason for the increased capacity is the lateral load dispersion in the case of protruding material for mid and edge position of  the loading element, see e.g.Madsen et al. [8], Augustin et al. [27], Bogensperger et al. [26], and Kathem et al. [15].Loading at the mid position results in the highest capacity due to the load being dispersed at both sides of the loaded length, along the fiber direction.This trend of increased capacity was predicted by all material models.Similar elastic behavior was observed for all models.However, the absolute values of the compressive force beyond the elastic limit were different.Hill's model yielded considerably higher and unrealistic forces as compared to Hoffman's and QMS models, which yielded very similar force--displacement curves, as was also found for the material testing results presented in Subsection 3.1.
The influence of the annual ring curvature, i.e., of the location of the pith, on stiffness and strength in CPG for compression in the radial and tangential material orientations was studied by using Hoffman's model considering the similar responses by Hoffman's and QMS models at both material testing and application levels.Results are shown in Fig. 15.Moreover, different material models were applied to compression in the radial direction and a comparison of modulus of elasticity and strength is summarized in Tables 4 and 5.
For the Cyl-2 pith location, Hoffman's yield criterion resulted in a MOE, E c,90 , of 686 N/mm 2 , 803 N/mm 2 and 924 N/mm 2 for Full, Edge and Mid position of the loading element, respectively.Corresponding strengths, f c,90 were found to be 4.16 N/mm 2 , 4.87 N/mm 2 and 5.40 N/ mm 2 for Full, Edge and Mid positions of the loading element, respectively.Thus, Hoffman's model predicted an increase of the stiffness E c,90 of 17% and 35% for Edge and Mid positions, respectively as compared to the Full loaded specimen.Moreover, the increase in strength f c,90 was equal to 17% and 30% for Edge and Mid positions as compared to the Full loaded specimen.The model predicted increase is slightly smaller than the increase found in experiments of Kathem et al. [15], who reported 26% increased strength for Mid position as compared to the Edge position.The model however even allowed to assess the influence of the annual ring orientation.Taking the latter into account, the corresponding strength increase was found to be 11.7%, 10.9%, and 12.2% for Cyl-1, Cyl-2, and Cyl-3 pith locations, respectively.A more pronounced influence of the annual ring orientation was found for MOE.For the fully loaded specimen (Full), the Cyl-2 pith location yielded a two times higher stiffness than the Cyl-3 pith location.Corresponding strength was found to be 15.6% higher for Cyl-2 than for Cyl-3, cf.Tables 4 and 5.

Conclusions
The material behavior of wood under transverse uniaxial loading and a combination of compressive stresses with rolling shear stress in the R-T plane was studied by FE modeling considering material elasto-plasticity.Three material models were used, including a built-in model for the used FE-software, namely Hill's plasticity model.A user-defined subroutine was used to implement the Hoffman criterion, as well as a novel quadratic multi-surface (QMS) criterion.A biaxial test setup for identification of material parameters and the engineering problem of compression perpendicular to the grain in structural wooden elements were simulated.
Hill's failure criterion was found to be the least suitable in the present study, since the cellular material wood has different tensile and compressive strengths in the three principal directions, something that cannot be modeled with Hill's criterion.Moreover, Hill's failure criterion yielded unrealistic stress re-distributions, which resulted in far too high reaction forces when modeling experimental setups.The Hoffman and the QMS criteria showed quite good agreement with the biaxial experiments in terms of force displacement relations and strain distributions for different loading situations, especially for combinations with radial compression.Differences in material response in the R-and Tdirections were noticed in the experiments in terms of stiffness, strength and hardening behavior.Consequently, consideration of only a small number of hardening parameters in the Hoffman and QMS models was not sufficient to reflect the different material response under compression, shear and combined stress interaction for both material orientations.A deviation in the evolution of plasticity for combined loading cases at higher strain levels was noticed in the FE-models, which could be the effect of not considering large strain plasticity.The influence of the annual ring curvature or pith location on elastic and shear stiffness and strength was found to be more pronounced in the case of radial compression as compared to tangential compression.The highest shear stiffness and shear strength and lowest elastic stiffness and force under compression were noticed in the specimen with the annual rings oriented at 45 • with the loading directions.
Application of the material models to compression perpendicular to the grain for studying the unloaded length effect, yielded similar trends as observed in structural tests.Differences in force-displacement response for different loading situations were well in line with literature findings and the influence of the annual ring curvature on the overall force-displacement behavior could be quantified.
In conclusion, consideration of different hardening parameters for plastic material behavior in the tangential direction and brittle shear failure could be relevant to consider for further development of the concept in order to achieve better correlation with experiments on clear wood.However, the application of Hoffman's material model in predicting the compression perpendicular to the grain behavior for different loading situations indicates the suitability of the models in predicting the response at the structural level, which are usually not exposed to very large strain.The presented findings will help to develop efficient wood-based products and design safe structural elements.

Table 4
Comparison of elastic moduli from different material models: influence of the unloaded length effect and annual ring orientation on radial compressive stiffness.
S.T. Akter et al.mm3 (b × h × l), see Fig.2(a).Experiments were realized in a MTS biaxial test frame, which was equipped with two servo-hydraulic actuators with a capacity of 100 and 50 kN (MTS Model 661.20F) for the vertical and horizontal orientation, respectively.

Fig. 4 .
Fig. 4. Specimens and loading details for the study of the unloaded length effect; (a) Mid position, (b) Edge position and (c) Full specimen.

Fig. 5 .Fig. 6 .
Fig. 5. Comparison of force vs. displacement from experiments and numerical modeling of R-specimens; (a) compression test and (b) shear test.

Fig. 7 .
Fig. 7. Comparison of (a) compressive force vs. displacement and (b) shear force vs displacement, from experiments and numerical modeling of combined stress interaction tests on R-specimens.

Fig. 8 .
Fig. 8.Comparison of (a) compressive force vs. displacement and (b) shear force vs displacement, from experiments and numerical modeling of combined stress interaction tests on T-specimens.

Fig. 9 .
Fig. 9. Comparison of shear stress vs. normal stress from experiments and numerical modeling of combined stress interaction tests on (a) R-specimens and (b) T-specimens.

Fig. 12 .
Fig. 12.Comparison of force vs. displacement from experiments and numerical modeling of R-specimen by considering Cartesian and cylindrical coordinate systems with different pith locations (a) compression test and (b) shear test.

Fig. 13 .
Fig. 13.Comparison of force vs. displacement from experiments and numerical modeling of T-specimen by considering Cartesian and cylindrical coordinate systems with different pith locations (a) compression test and (b) shear test.

Fig. 14 .
Fig.14.Force-displacement plots from unloaded length effect tests of solid wood (radial compression) for different loading positions and cyl-2 pith location as predicted with different material models.

Fig. 15 .
Fig. 15.Force-displacement behavior as predicted by Hoffman's material model for unloaded length effect test of solid wood for different pith locations; (a) radial compression and in (b) tangential compression.

Table 1
Elastic material properties for FE modeling.

Table 2
Material strength for FE modeling.

Table 3
Mesh dependency for modelling radial compression with Hoffman's model.