Skip to content
BY 4.0 license Open Access Published online by De Gruyter August 4, 2022

Elasto-plastic material model of oak at two moisture content levels

Jan Tippner, Jaromír Milch, Václav Sebera and Martin Brabec
From the journal Holzforschung

Abstract

The mechanical properties of wood show a very high dependence on the moisture content (MC). A consideration of MC in numerical simulations increases the applicability of such prediction with respect to application and moisture states of the wood material. The goal of this work is to develop an accurate orthotropic elasto-plastic model for oak wood (Quercus robur L.) at two different MC levels applicable for finite element analysis (FEA). To achieve this goal, the following steps were carried out: (a) in-house standard specimens tests in compression, tension, and shear and in all three orthogonal directions, followed by three-point bending, where all specimens were conditioned to a 12 and 25.6% MC, prior to the mechanical test; (b) integration of all obtained material characteristics into the consistent numerical material models; (c) validation of the developed material models by comparing the numerically predicted values with the experimental ones; and (d) iterative calibration of the material models by adjusting the individual material characteristics to minimize error using a reference. Material models were successfully developed with the following mean relative errors: 5.2% for 12% MC and 5.8% for 25.6% MC, respectively. Both numerical material models consistently predicted the oak elasto-plastic response that can be easily integrated into any FEA.

1 Introduction

Wood is a complex material because it possesses a hierarchical structure and anisotropic physical and mechanical behaviors. In terms of its elastic mechanical behavior, wood is widely considered both orthotropic and anisotropic. However, once the elasto-plastic orthotropic and anisotropic description of wood behavior is of concern, the situation becomes more complicated even though a consistent theory background has been developed (Hill 1950; Tsai and Wu 1971). Knowing the wood elasto-plasticity, especially at various moisture content (MC) levels, is important since it allows a higher accuracy to be achieved for a numerical prediction of the behavior of wooden materials, wood connections, or new and historical timber structures in different or cyclic climates (Clouston and Lam 2002; Patton-Mallory et al. 1997; Yoshihara and Ohta 1992, 1994).

One method of coping with wood plasticity is to incrementally adjust the stiffness beyond the elastic limit with respect to the loading mode and anatomical direction as showed for Japanese ash in Tabiei and Wu (2000). Another convenient method to use when dealing with finite elements is to analyze the yielding stresses in every computed step and adjust the properties at the failure location. This approach applied on composites made of hardwood and with the help of orthotropic bilinear plasticity and the Weibull theory of weakest link was successfully presented by Moses and Prion (2004). However, recently, multi-surface plasticity (MSP) models were developed to describe wood behavior after the elastic limit, such as the constitutive model developed for spruce using biaxial tests in the longitudinal-radial (LR) plane by Mackenzie-Helnwein et al. (2003). Their MSP model utilized four surfaces representing four basic failure modes, which enabled the plasticity to be predicted using the hardening in radial compression and post-failure behavior in tension parallel to fiber. However, the use of an orthotropic single-surface formulation using a non-associative hardening/softening law may still be suitable for an analysis of both in-plane and out-of-plane shell problems, as shown by Mackenzie-Helnwein et al. (2005), and when applied to real applications (Fleischmann et al. 2007). A three-dimensional orthotropic MSP model for spruce was presented later by Schmidt and Kaliske (2006). The researchers employed the exponential function for post-failure behavior in tension using softening and in compression perpendicular to fibers using hardening. Oudjene and Khelifa (2009) created a constitutive law, including hardening associated with the densification phase, to analyze the compression of wood and dowel connections and found a good agreement with the experiments.

A post-elastic analysis of wood may be carried out using damage and strength approaches that consider the microstructure at various levels, as presented by Qing and Mishnaevsky (2011). More recently, the microstructure-based multisurface failure criterion that enables linking of spruce wooden cells with unit cell at the annual ring scale was developed by Lukacevic et al. (2017). The Tsai–Wu plasticity can be conveniently combined with fracture mechanic approaches, as shown by Danielsson and Gustafsson (2013), who utilized a cohesive zone model to predict the plastic behavior of spruce and pine wood when loaded perpendicular to the grain. Yoshihara (2014) used Hill’s plasticity theory (Hill 1950) to analyze the tangential strains during compression along the fiber direction for Sitka spruce and Japanese birch. His analysis demonstrated more nonlinearity in the relationship between longitudinal and tangential directions than in the plasticity theory. Pěnčík (2015) shows how to utilize the orthotropic plasticity of pine wood with various cylindrically oriented tree-rings in a four-point bending problem using finite elements, and Sirumbal-Zapata et al. (2018) reveals that plasticity may be approached using plasticity-damage constitutive models. Their model utilized plasticity yield surfaces with isotropic hardening and a derived closed form, which enabled them to predict wood behavior under cyclic loads with high strains that occur in connections. A combination of fracture and plasticity approaches to predict the collapse response of the two-dimensional (2D) cellular structure of wood can be conveniently applied to the cell wall level, as presented by Scheperboer et al. (2019). Recently, Benvenuti et al. (2020a) developed a novel orthotropic multi-surface model by coupling anisotropic damage and plasticity constitutive laws that allowed them to correctly capture brittle failure induced by strain localization and the possible occurrence of ductile plastic behavior. Although the effect of moisture was not studied, it could be implemented in this model by adding a new term to the stress potential. Their model was successfully validated against experimental embedment tests that simulated dowel-like connections between metal and wood. The same group of authors then applied that model to various tests in mode I and mixed modes of failure of both softwoods and hardwoods with a very high accuracy compared to the experimental data (Benvenuti et al. 2020b). On an engineering and more practical level, wood plasticity may be understood in terms of the strength design values, as shown by Marcolin et al. (2021), who developed an analytical plasticity model that linked strength values in the bending, compression, and tension of several hardwood Brazilian species. To employ the MSP model, one may omit finite elements and employ the material point method (MPM) instead, as shown by Adibaskoro et al. (2022), who utilized the MSP model from Schmidt and Kaliske (2006).

The MC effect on post-elastic behavior is very high (Forest Products Laboratory 2021; Kollmann and Côte 1968; Tukiainen and Hughes 2016) as wood becomes much more compliant and has a greater proportion of plastic strain within the total strain. Hering et al. (2012) analyzed European beech wood’s elasto-plasticity in compression loading at four different MC levels using the Ramberg–Osgood approach and MSP. They obtained the MSP parameters using an inverse finite element analysis (FEA) and concluded that nonlinear material behavior played a major role in the development of moisture-induced stress. Pambou Nzienguia et al. (2017) studied the tropical hardwood species Okume (Aucoumea klaineana Pierre) and reported that variations in the MC (below and above 30%) negligibly influence the plasticity domain in compression parallel to fiber. However, despite the need for elasto-plastic material data for hardwoods at higher MCs, this information is largely missing in literature.

Wood plasticity at high strain rates is also of scientific interest. In this respect, Baumann et al. (2019) and Müller et al. (2020) used analogies from fiber composite strength principles to model three hardwoods (ash, beech, and birch) for an application in the automotive field using explicit FE codes. More recently, Šebek et al. (2020, 2021 developed and experimentally verified elasto-plastic damage models of European beech wood for dynamic and explicit FE analyses with the help of the split Hopkinson pressure bar.

Therefore, the aim of this study is to determine the material characteristics of English oak (Quercus robur L.) under two MC levels to develop an elasto-plastic material model with different behavior in compression and tension that is applicable for an FEA. The main hypotheses of the work are that: (1) increased MC close to the fiber saturation point influences both elastic and plastic strain characteristics of English oak wood; and (2) it is possible to obtain complex material with low relative error validated on a basis of more loading types. The partial steps to achieve the research aim are as follows: (1) to obtain elasto-plastic material characteristics through in-house experiments; (2) to evaluate the non-linear behavior of the elasto-plastic material model based on these characteristics; and (3) to optimize and calibrate the elasto-plastic material models by changing the material characteristics based on the minimization of relative differences between experimental and numerical datasets.

2 Materials and methods

2.1 Materials and experiments

English oak (Q. robur L.) was tested regarding compression and tension by adjusting the methodology of Brabec et al. (2015) and Milch et al. (2016) in the following ways: parallel (‖) and perpendicular (⊥) to the grain in the radial (R) and tangential (T) directions; in shear at longitudinal–tangential (LT), longitudinal–radiradial LR, and radial–tangential (RT) planes; and by three-point bending. Mechanical tests were performed on small, clear timber specimens corresponding to the standard BS 373 (1957; Figure 1). The specimens were conditioned in a climatic chamber at 20 °C and 65% relative humidity (RH) and at 20° and close to 100% RH until 12 and 25.6% equilibrium MC (EMC) was reached. The EMC and average oven-dry density of 738 and 754 kg m−3 were measured gravimetrically according to (ASTM D2395 2014).

Figure 1: 
Schemes of the experiments (adopted from Milch et al. 2016): (A) compression parallel to grain, (B) compression perpendicular to the grain, (C) 3-point bending, (D) tensile parallel to grain, (E) tensile perpendicular to grain, (F) shear at individual directions.

Figure 1:

Schemes of the experiments (adopted from Milch et al. 2016): (A) compression parallel to grain, (B) compression perpendicular to the grain, (C) 3-point bending, (D) tensile parallel to grain, (E) tensile perpendicular to grain, (F) shear at individual directions.

The mechanical tests were carried out using the universal testing machine Zwick Z050/TH 3A (Zwick Roell AG, Ulm, Germany) equipped with a 50 kN load cell. The experiment procedures were set and controlled by TestXpert v.11.02 (Zwick Roell AG, Ulm, Germany). The temperature was kept within a range of 20–22 °C during all mechanical tests. Normal air humidity (50–60%) in the testing room was maintained. To avoid moisture content losses, all specimens were tested immediately after removing from the climatic chamber, and, according to respected testing standards, the mechanical test were performed within 60 ± 30 s.

The compression and tension (‖ and ⊥, resp.) of the grain were tested using the specimens specified in Figure 1A–E and Table 1. The experimental test evaluation was derived from the Czech technical norms ČSN 49 0110–14 (ČSN 49 0110 1980; ČSN 49 0111 1992; ČSN 49 0112 1979; ČSN 49 0113 1992; ČSN 49 0114 1992). The compression strength (σ L,R,T,c ), normal elastic moduli (E L,R,T,c ) ‖ and ⊥ to the grain (Figure 1A and B), the tension strength (σ L,R,T,t ), normal elastic moduli (E L,R,T,t ) ‖ and ⊥ to the grain (Figure 1D and E) were calculated according to Equations (1) and (2):

(1) σ L , c ( t ) = F max R T ; σ R , c ( t ) = F yield L T ; σ T , c ( t ) = F yield L R

where F max is the maximum loading force (N), F yield is the loading force corresponding to the yield strength (N), T is the specimen dimension in the tangential direction (mm), R is the specimen dimension in the radial direction (mm), and L is the specimen dimension in the longitudinal direction (mm).

(2) E L , c ( t ) = ( F 40 F 10 ) l e x t R T ( u 40 u 10 ) E R , c ( t ) = ( F yield , 40 F yield , 10 ) l ext L T ( u yield , 40 u yield , 10 ) E T , c ( t ) = ( F yield , 40 F yield , 10 ) l ext L R ( u 40 u 10 )

where F 40 and F 10 are the forces at the 40 and 10% levels of F max (N), F yield,40 and F yield,10 are the forces at the 40 and 10% levels of F yield (N), l ext is the initial distance of the tracked points (mm), u 40 and u 10 are the relative changes of l ext at the force levels of F 40 and F 10 (mm), and u yield,40 and u yield,10 are the relative changes of l ext at the force levels of F yield,40 and F yield,10 (mm). When the specimens were compressed and tensed ‖ to the grain, the u 40 and u 10 were tracked by isolated points using “clip-on” extensometers (Zwick Roell AG, Ulm, Germany; Figure 1A). For the u yield,40 and u yield,10 in compression and tension ⊥ to the grain, the isolated points were tracked by the virtual extensometer (adjusted from Brabec et al. 2015) based on the digital image correlation (DIC) technique.

Table 1:

Dimensions (in mm) of individually tested specimens.

n L L 0 R T l ext l s Rad
Compression || to grain (A) 15 30 20 20 10
Compression ⊥ to grain in R direction (B) 15 30 20 20 R
Compression ⊥ to grain in T direction (B) 15 30 20 20 T
3-Point bending (C) 30 300 20 20 240
Tension || to grain (D) 15 350 20 4 50 62.5
Tension ⊥ to grain in R direction (E) 15 30 60 50 20 25 15
Tension ⊥ to grain in T direction (E) 15 30 60 20 50 25 15
Shear in LT plane (F) 15 30 50 30 20 1
Shear in LR plane (F) 15 30 50 20 30 1
Shear in RT plane (F) 15 50 20 4 T

2.2 Three-point static bending

Three-point static bending was carried out using the specimens specified in Figure 1C and Table 1, and all specimens were tested in a tangential direction. The deflection in bending was measured using a standard “clip-on” deflectometer (Zwick Roell AG, Ulm, Germany) that was clamped to the bottom surface at the middle point of the transverse dimension located at the midspan. The experimental test evaluation was derived from the Czech technical norms ČSN 49 0115 (ČSN 49 0115 1979) and ČSN 49 0116 (ČSN 49 0116 1986), with a modified span-to-depth ratio equal to 12. The bending strength (σ b ) was calculated from the F max, as shown in Equation (3):

(3) σ b = 3 F max l s 2 R T 2

where l s is the span of supports (mm).

The modulus of elasticity (E b ) was calculated using Equation (4) as follows:

(4) E b = ( F 40 F 10 ) l s 3 4 R T 3 ( u def , 40 u def , 10 )

where u def,40 and u def,10 are the deflections at force levels F 40 and F 10 (mm).

2.3 Shear tests

The shear tests were derived from Czech technical norm ČSN 49 0118 (ČSN 49 0118 1980). This approach was chosen in comparison to others (e.g., EN 789) since it provides a reliable investigation of shear strength for the purpose of the plastic response region definition. The LT, LR, and RT planes were conducted using the specimens specified in Figure 1F and Table 1. The LT and LR shear strains were calculated as engineering shear strains according to Equation (7) based on the relative rotation between the two (movable and static) isolated points, which were located 1 mm apart. The movable one was located near the shear plane and tracked by the virtual extensometer (Correlated Solutions Inc., Columbia, SC, USA). The movement of the static point near the fixed support was zero (Figure 1F). The shear strength (yield stress–τ ij ) was calculated as shown by ČSN 49 0118 (ČSN 49 0118 1980), i.e. Equation (5), and the shear elastic moduli (G ij ) were calculated from strains and force development following Czech technical standards e.g., (ČSN 49 0111 1992), as shown in Equation (6):

(5) τ L T = F max L T ; τ L R = F max L R ; τ R T = F max R T

(6) G i j = F 40 F 10 ϵ i j 40 ϵ i j 10

(7) ϵ i j = 2 arctg ( u ) l ext 100

where ε ij is the shear strain (−), ε ij,40 and ε ij,10 are shear strains at force levels F 40 and F 10 (−), and u is a displacement in the load direction of the tracked point (mm).

2.4 Optical measurements

The full-field optical technique utilizing the DIC was employed to measure the full-field strain and relative displacement of virtually tracked points, as shown by Brabec et al. (2015) and Milch et al. (2017). The acquisition set consisted of two cameras (AVT Stingray Copper F-504B, Allied Vision Technologies, Osnabrück, Germany, cell size of 3.45 μm and resolution of 2452 × 2056 pixels) equipped with lenses (Pentax C2514-M, Pentax Precision Co., Ltd., Tokyo, Japan, focal length of 25 mm) in the stereo vision three-dimensional (3D) configuration (Figure 5). The area of interest (AoI) on the specimen surface was covered by a high-contrast black and white speckle pattern using synthetic sprays. The AoI of the specimens was illuminated by two light sources by SobrietyCube 360 (Sobriety s.r.o., Kuřim, Czech Republic) and fitted with LED sensors (Luminus Phlatlight CSM-360, 90 W Luminus Devices Inc., Billerica, MA, USA). The images were captured at the acquisition interval of 0.25 s. The relative displacement of the tracked points obtained from the DIC was then used for the following: a) to obtain the load-displacement curves; b) to calculate the developed strains (relative changes of l ex) during testing; and c) to obtain the boundary conditions (BCs) for the subsequent FEA. The full-field strains (relative displacements) obtained from the DIC were used for the calculation of Poisson’s ratios. The calculation of the displacements from the images was performed by software Vic-3D v. 2010 (Correlated Solutions Inc., Columbia, SC, USA). Using non-contact video techniques to obtain Poisson’s ratio is comparable to standard contact techniques (clip-on extensometers and strain gages) or older non-contact techniques, such as Electronic Speckle Pattern Interferometry (ESPI), as presented by Kumpenza et al. (2018).

2.5 The Poisson’s ratios

Elastic strains calculated by Green-Lagrangian strain tensor evaluated by optical full-field measurement, served for the Poisson’s ratios determination. The estimation of Poisson’s ratios was based on the active, i.e., compressive (parallel to loading), and passive, i.e., tensile (perpendicular to loading), elastic strains in each correlated point, as shown in Equation (8).

(8) v A P = ϵ P ϵ A

where ν AP is Poisson’s ratio in the XY plane (A stands for active and P stands for passive), ε A is the strain parallel to loading (active), and ε P is the strain perpendicular to loading (passive).

The final resulting values were obtained by averaging Poisson’s ratios over the whole AoI, where the extremes of Poisson’s ratios (lower than zero and higher than one) obtained in certain regions within the AoI were excluded from the averaging process. As Hearmon’s (1948) notations were used, Poisson’s ratios for LR, RT, and LT are the major ones.

2.6 Finite element analysis

To validate and subsequently calibrate the input material characteristics in the elasto-plastic range, the FEA of the compression and bending tests based on 3D solid models was performed. The finite element models of laboratory-sized clear wood specimens were built in the ANSYS v. 14.5 computing system using the ANSYS Parametric Design Language (ANSYS Inc., USA). The 3D FE analysis of the compression and bending tests precisely reflected the physical experiments in the specimens’ dimensions, geometry of support and loading heads, and applied load. The wood was modeled as a 3D orthotropic elasto-plastic material with non-linear isotropic hardening, where different behaviors in compression and tension (yield stress and tangent moduli) were considered. Figure 2 depicts the workflow for the experimental and numerical analyses.

Figure 2: 
Flowchart for experimental and numerical work (adopted from Milch et al. 2016).

Figure 2:

Flowchart for experimental and numerical work (adopted from Milch et al. 2016).

The elasto-plastic orthotropy was modeled using the TB, ANISO function in APDL, which is based on the generalized anisotropic Hill plasticity theory (Hill 1950). Hill’s criterion is a generalized form of the von Mises criterion, which accurately describes the anisotropy of the materials. The FE mesh was built using the 3D hexahedral 20-node structural solid element (SOLID95) with plasticity capabilities. The element size was set to 1 mm in the compression tests (12,000 of FEs) and 2 mm in the three-point bending test (15,000 of FEs). In the compression tests, the virtual specimens were placed between two areas (one fixed; the other free in the L, R, and T load directions) with rigid target contact surfaces. In the three-point bending test, the support and loading heads were modeled as rigid surfaces with the same settings as the compression tests. The friction coefficient (μ) was set to 0.13, and the contact algorithm was defined as the augmented Lagrange method, with the normal contact stiffness factor 1.0 for all analyses.

The MRE in Equation (9) was calculated by comparing the product of the optimization loop (set of material properties) with the experimental dataset by considering material parameters. Therefore, the MRE represents an arithmetic average of certain relative errors computed from single material properties.

(9) MRE = 100 | E exp E num | / E exp

3 Results and discussion

Fundamental orthotropic elasto-plastic material models of English oak were obtained from in-house experiments and assembled into two sets of material models for the studied EMC levels. Each set, which comprised 27 characteristics, was the result of the iterative procedure (Figure 2) and is listed in Tables 2 and 3. The values of presented characteristics differed from publications providing sets of material properties (Forest Products Laboratory 2021; Hearmon 1948; Kollmann and Côte 1968). The modification of measured values was the point of an iterative procedure, influenced by theory limitations and a model validated response. The final sets of all characteristics represented the necessary inputs for the FE analyses using the bilinear orthotropic approach, which was convenient for both the proprietary computing system used in the presented work (ANSYS®) and other systems. Therefore, it is necessary to consider presented characteristics as a tool for numerical study, not as a reference of material properties in the sense of measured values.

Table 2:

Set of final material characteristics of English oak at EMC = 12.0%.

Data English oak – 12% EMC
ρ av  = 738 (kg/m3)
Elastic
L
R
T
E i [MPa] 14 669 656 494
LR RT LT
G ij (MPa) 1900 127 2015
ν ij (−) 0.037 0.326 0.028

Plastic L R T

σ i,Te (MPa) 100.013698041898a 10.13 9.83
E i,Te,Tan (MPa) 1.0001 0.1013 0.0983
σ i,Co (MPa) 62.00 11.30 9.45
E i,Co,Tan (MPa) 1.24 1.13 5.884

LR RT LT

τ ij (MPa) 6.91 4.75 6.14
G ij, Tan (MPa) 0.1382 0.0950 0.1228

  1. ρ av is average density, E i is normal elastic modulus, G ij is shear elastic modulus, ν ij is Poisson’s ratio, σ i,Te is yield stress in tension, E i,Te,Tan is normal tangent modulus in tension, σ i,Co is yield stress in compression, E i Co, Tan is normal tangent modulus in compression, τ ij is shear yield stress, G ij,Tan is shear tangent modulus. Indices i and j characterize anatomical directions L, R, and T. aThis value must be defined with all decimal points to meet conditions in ANSYS® Mechanical APDL for TB, ANISO function that is based on Hill’s plasticity.

Table 3:

Set of final material characteristics of English oak at EMC = 25.6%.

Data English oak – 26.6% EMC
ρ av  = 754 (kg/m3)
Elastic
L
R
T
E i (MPa) 9250 158 130
LR RT LT
G ij (MPa) 858 112 905
ν ij (−) 0.37 0.31 0.54

Plastic L R T

σ i,Te (MPa) 60.671234882694a 5.01 5.44
E i,Te,Tan (MPa) 0.001 0.001 0.001
σ i,Co (MPa) 27.20 6.65 4.70
E i,Co,Tan (MPa) 0.0002 0.0002 0.0002

LR RT LT

τ ij (MPa) 3.5 4.7 4.7
G ij,Tan (MPa) 0.001 0.001 0.001

  1. ρ av is average density, E i is normal elastic modulus, G ij is shear elastic modulus, ν ij is Poisson’s ratio, σ i,Te is yield stress in tension, E i,Te,Tan is normal tangent modulus in tension, σ i,Co is yield stress in compression, E i,Co,Tan is normal tangent modulus in compression, τ ij is shear yield stress, G ij,Tan is shear tangent modulus. Indices i and j characterize anatomical directions L, R, and T. aThis value must be defined with all decimal points to meet conditions in ANSYS® Mechanical APDL for TB, ANISO function that is based on Hill’s plasticity.

A visual and numerical comparison between the FE predicted (red curves) and experimentally measured (black curves) mechanical behavior of oak wood at an EMC equal to 12.0% in compression and bending is shown in Figure 3. Figure 3A shows the compression parallel to the grain, and the FE predicted response lies in the experimental ranges. Figure 3B and C show the numerically predicted compression in the radial and tangential directions, respectively, as a bilinear response that fits the experimental data. Even though the transversal plastic moduli might be higher in a reality for both observed MC’s, the study showed the numerical prediction based on such small values which exhibited quite low errors. For the sample of MC = 25.6%, the transversal plastic behavior was nearly perfectly plastic as plastic moduli are nearly zero when compared to elastic moduli. Hong et al. (2016) also present application of small values of transversal plastic moduli (as 1% of elastic moduli) with satisfying agreement between FEA and experiments. Tabiei and Wu (2000) show that low plastic tangent moduli a nearly perfect plastic plateau lead to a reliable prediction of a transversal behavior of Japanese ash at low strain rates. The last stage of the radial compression (densification) is not visible in Figure 3B because the maximal strains were kept below 6%, which is too low for densification to occur. Nairn (2006) shows that transversal densification often begins after 10% of strain, depending on species and cell wall yield strength. The compression yield values were found higher in the radial direction compared to the tangential one, which represents a typical sign for hardwood species due to radially oriented pith rays that function as a reinforcement. This corresponds not only to experimental findings but also numerical predictions using the MPM applied on real microscopical morphology of wood (Nairn 2006).

Figure 3: 
Force-displacement/deflection curves of compression parallel to the grain (A) and perpendicular to the grain in R (B), T (C) directions and bending (D). The red curves represent the FE predictions, and the black curves represent the experimental elasto-plastic deformation behavior of oak at 12.0% EMC. The blue curve (E) characterizes the MRE during the calibration procedure; the red crosses show the MRE, and table (F) shows relative errors of particular parameters in the individual calibration steps (S1–S5).

Figure 3:

Force-displacement/deflection curves of compression parallel to the grain (A) and perpendicular to the grain in R (B), T (C) directions and bending (D). The red curves represent the FE predictions, and the black curves represent the experimental elasto-plastic deformation behavior of oak at 12.0% EMC. The blue curve (E) characterizes the MRE during the calibration procedure; the red crosses show the MRE, and table (F) shows relative errors of particular parameters in the individual calibration steps (S1–S5).

Figure 3D shows a bending prediction, which represents the interaction between optimized bilinear models in normal and shear directions. When the maximal bending force is within the limit, the deflection at maximal force is rather conservative as it shows lower values than the experiments. There are two main reasons for this: (1) the material reaches plastic domain in bending due to plastic strain first developed in compression parallel to fiber; and (2) the experiments showed higher deflections than FE models due to the surface imperfections and roughness of both the experimental device and specimens. This was also related to higher relative errors of bending characteristics (E b and σ b ) in the optimization procedure (Figure 3), which somehow searches for a trade-off between all particular characteristics to be minimized in terms of MRE.

Figure 3E shows a progression of the MRE during the optimization of the material model, starting with 56.6% and reaching its final value of 5.2%. The final MRE can be considered low due to the high natural variability of wood that was, among others, clearly visible in all experiments and the number of material parameters that entered the optimization loop (Tables 2 and 3). The main material parameters entering the optimization loop (see Figure 3) showed that the transversal elastic moduli and transverse yield strengths underwent the biggest adjustments in the process to lower the global MRE. For instance, the E R adjusted its MRE from 230% in the first iteration to 2.7% in the final one. In the global MRE, we may claim that the calibrated numerical material model faithfully represented the mechanical behavior of oak in the elastic and plastic ranges at an MC of 12%.

Like the 12% MC, the optimization result of the numerical material model for English oak at an EMC of 25.6% is shown in Table 3. A higher MC of wood lowered most elastic and strength parameters in a statistically significant manner and was generally assumed when more MC was present in the wood cell walls. E L , E R , and E T at an MC of 25.6% were 63, 24, and 26% of the ones at an MC of 12%, respectively. For yield stresses, the values at an MC of 25.6% were 60 and 44% of the 12% MC oak properties for stress parallel to the fiber in tension and compression, respectively. A similar decrease was seen for both shear elastic moduli and shear yield stresses. Even though MRE is a mean value and a general metric to assess the error of our material models, it is convenient because it takes into account eight different material parameters in various anatomical directions, and it includes uniaxial (compression) and combined (bending) modes of loading. In addition, many studies have assessed the error of numerical simulation or material model only by using of visual comparison of loading diagrams of a certain loading mode (Benvenuti et al. 2020a, 2020b; Oudjene and Khelifa 2009; Pěnčík 2015). However, using material parameters obtained only from standard uniaxial tests does not always provide satisfying simulation results, especially when a multiaxial stress state occurs (Hong et al. 2016). Therefore, this work, following the study from Milch et al. (2016), focused on finding a set of material models using the MRE because it directly reflects the compression and bending phenomena.

As shown in Figure 4E, the numerical material model of oak with 25.6% EMC before optimization exhibited a significantly better prediction than the model for 12% EMC before optimization (Figure 3E). Nevertheless, the final numerical material model for English oak at 25.6% MC has relative errors in a range of 1–10%, which is a slightly broader range than for oak at 12% MC, which is 2–8%.

Figure 4: 
Force-displacement/deflection curves of compression parallel to the grain (A) and perpendicular to the grain in R (B), T (C) directions and bending (D). The red curves represent the FE predictions, and the black curves represent the experimental elasto-plastic deformation behavior of oak at 25.6% EMC. The blue curve (E) characterizes the MRE during the calibration procedure; the red crosses show the MRE, and table (F) shows relative errors of particular parameters in the individual calibration steps (S1–S5).

Figure 4:

Force-displacement/deflection curves of compression parallel to the grain (A) and perpendicular to the grain in R (B), T (C) directions and bending (D). The red curves represent the FE predictions, and the black curves represent the experimental elasto-plastic deformation behavior of oak at 25.6% EMC. The blue curve (E) characterizes the MRE during the calibration procedure; the red crosses show the MRE, and table (F) shows relative errors of particular parameters in the individual calibration steps (S1–S5).

Figure 5: 
Experimental test configuration: stereo vision optical system (3D) of the tensile test perpendicular to the grain.

Figure 5:

Experimental test configuration: stereo vision optical system (3D) of the tensile test perpendicular to the grain.

Moreover, when different behavior in the compression and tension of wood was considered, it provided a more realistic prediction of bending, especially beyond the proportional limit. As the bending behavior combines compression, tension, and shear, such a prediction improvement was expected. However, this improvement was at the expense of the less precise (overvalued) material characteristics of the tension behavior of wood loaded perpendicularly to the grain (in both R and T directions). However, wood is loaded by tension perpendicular to the grain only in a few specific and rare loading cases compared to bending, which is the most common one in practice.

The numerical material models represent datasets that may be easily implemented in FE software packages, which could improve predictions for both dried wood and wood with an MC close to the fiber saturation point. The applications of material model of nearly green wood goes into the processing of the raw material as well as to the numerical simulations of tree biomechanics, which utilizes the FEA frequently (Dupuy et al. 2007; Rahardjo et al. 2014; Vojáčková et al. 2019; Yang et al. 2014).

4 Conclusions

The work presents a broad range of experiments that provided suitable datasets to model the orthotropic mechanical behavior of oak wood at two different MC levels. The datasets represent the elasto-plastic material models for oak wood, including full orthotropic elasticity and plasticity with non-linear isotropic hardening, where different compression and tensile characteristics, such as yield stress and tangent moduli, are considered. Validation based on the complex bending loading mode showed a very reliable agreement with nonlinear experimental outputs. However, the robust and straightforward bi-linear representation at the level of uniaxial and shear stress-strain response was used, which may be the limitation for applications just in modes close to uniaxial or shear loading. The detailed verification of the presented material model in these loading modes can be the objective of future studies. Further work should focus on calibrating the material models in terms of transversal yield strengths that are still somewhat higher than in a reality for a sake of model numerical consistency given by Hill’s plasticity theory. For the applications in tension and shear modes of loading, the addition of fracture damage criteria makes a natural next step in the model development. The assembled fundamental elasto-plastic material models were validated and subsequently successfully calibrated based on the measured data obtained from in-house experiments on clear small wood specimens. This fact determines the current application area. The influence of defects (such as knots, drying cracks, and wave grain) requires more extensive studies. The predicted load-displacement/deflection curves representing the mechanical behavior of oak wood showed an excellent agreement with the experimental data, especially when the natural variability of wood properties is considered. The MREs of both models were 5.2 and 5.8% for 12% MC and 25.6% MC, respectively. Both numerical material models represent the reliable prediction of the English oak elasto-plastic response, including the ultimate strength, which can be easily integrated into any FEA. The models may be used to analyze variety of wooden elements and constructions or, in the case of moist material models, for an analysis of green wood at the level of tree behavior in plant biomechanics research.


Corresponding author: Jan Tippner, Department of Wood Science and Technology, Faculty of Forestry and Wood Technology, Mendel University in Brno, Zemědělská 3, 613 00 Brno, Czech Republic, E-mail:

Funding source: Ministry of Education Youth and Sports in Czech Republic

Award Identifier / Grant number: ERC CZ #LL1909

  1. Author contributions: All the authors have accepted responsibility for the entire content of this submitted manuscript and approved submission.

  2. Research funding: This work was created at the Research Center Josef Ressel, Brno-Útěchov, Mendel University in Brno and funded by the project of the Ministry of Education Youth and Sports in Czech Republic (grant number ERC CZ #LL1909).

  3. Conflict of interest statement: The authors declare no conflicts of interest regarding this article.

References

Adibaskoro, T., Sołowski, W., and Hostikka, S. (2022). Multi-surfaced elasto-plastic wood material model in material point method. Int. J. Solid Struct. 236–237: 111333, https://doi.org/10.1016/j.ijsolstr.2021.111333.Search in Google Scholar

ASTM D2395 (2014). Standard test methods for density and specific gravity (relative density) of wood and wood-based materials. Philadelphia, PA, USA: American Society for Testing and Materials.Search in Google Scholar

Brabec, M., Tippner, J., Sebera, V., Milch, J., and Rademacher, P. (2015). Standard and non-standard deformation behaviour of European beech and Norway spruce during compression. Holzforschung 69: 1107–1116, doi:https://doi.org/10.1515/hf-2014-0231.Search in Google Scholar

BS 373 (1957). British standard: methods of testing small clear specimens of timber. London: British Standard Institution.Search in Google Scholar

Baumann, G., Hartmann, S., Müller, U., Kurzböck, C., and Feist, F. (2019). Comparison of the two material models 58, 143 in LS Dyna for modelling solid birch wood. In: 12th European LS-DYNA conference 2019, Koblenz, Germany.Search in Google Scholar

Benvenuti, E., Orlando, N., Gebhardt, C., and Kaliske, M. (2020a). An orthotropic multi-surface damage-plasticity FE-formulation for wood. Part I: constitutive model. Comput. Struct. 240: 106350, https://doi.org/10.1016/j.compstruc.2020.106350.Search in Google Scholar

Benvenuti, E., Orlando, N., Gebhardt, C., and Kaliske, M. (2020b). An orthotropic multi-surface damage-plasticity FE-formulation for wood. Part II: numerical applications. Comput. Struct. 240: 106351, https://doi.org/10.1016/j.compstruc.2020.106351.Search in Google Scholar

Clouston, P.L. and Lam, F.A. (2002). Stochastic plasticity approach to strength modeling of strand-based wood composites. Compos. Sci. Technol. 62: 1381–1395, https://doi.org/10.1016/s0266-3538(02)00086-6.Search in Google Scholar

ČSN 49 0110 (1980). Drevo. Medza pevnosti v tlaku v smere vlákien [in Slovak].Search in Google Scholar

ČSN 49 0111 (1992). Skúšky vlastností rastlého dreva. Metóda zisťovania modulu pružnosti v tlaku pozdĺž vlákien [in Slovak].Search in Google Scholar

ČSN 49 0112 (1979). Drevo. Tlak naprieč vlákien [in Slovak].Search in Google Scholar

ČSN 49 0113 (1992). Skúšky vlastností rastlého dreva. Metóda zisťovania pevnosti v ťahu pozdĺž vlákien [in Slovak].Search in Google Scholar

ČSN 49 0114 (1992). Skúšky vlastností rastlého dreva. Metóda zisťovania pevnosti v ťahu naprieč vlákien [in Slovak].Search in Google Scholar

ČSN 49 0115 (1979). Drevo. Zisťovanie medze pevnosti v statickom ohybe [in Slovak].Search in Google Scholar

ČSN 49 0116 (1986). Drevo. Metóda zisťovania modulu pružnosti pri statickom ohybe [in Slovak].Search in Google Scholar

ČSN 49 0118 (1980). Drevo. Medza pevnosti v šmyku v smere vlákien [in Slovak].Search in Google Scholar

Danielsson, H. and Gustafsson, P.J. (2013). A three-dimensional plasticity model for perpendicular to grain cohesive fracture in wood. Eng. Fract. Mech. 98: 137–152, https://doi.org/10.1016/j.engfracmech.2012.12.008.Search in Google Scholar

Dupuy, L.X., Fourcaud, T., Lac, P., and Stokes, A.A. (2007). Generic 3D finite element model of tree anchorage integrating soil mechanics and real root system architecture. Am. J. Bot. 94: 1506–1514, https://doi.org/10.3732/ajb.94.9.1506.Search in Google Scholar PubMed

Fleischmann, M., Krenn, H., Eberhardsteiner, J., and Schickhofer, G. (2007). Numerische Berechnung von Holzkonstruktionen unter Verwendung eines orthotropen elasto-plastischen Werkstoffmodells. Holz Roh-Werkst. 65: 301–313, https://doi.org/10.1007/s00107-007-0185-5.Search in Google Scholar

Forest Products Laboratory (2021). Wood handbook: wood as an engineering material. General Technical Report FPL-GTR-282. Madison, WI: U.S. Department of Agriculture, Forest Service, Forest Products Laboratory, p. 543.Search in Google Scholar

Hearmon, R.F.S. (1948). The elasticity of wood and plywood. Special Report no. 7. London: HM Stationery Office.Search in Google Scholar

Hering, S., Saft, S., Resch, E., Niemz, P., and Kaliske, M. (2012). Characterisation of moisture-dependent plasticity of beech wood and its application to a multi-surface plasticity model. Holzforschung 66: 373–380, https://doi.org/10.1515/hf.2011.162.Search in Google Scholar

Hill, R. (1950). The mathematical theory of plasticity. Oxford: The Oxford Engineering Science Series, pp. 97–114.Search in Google Scholar

Hong, J.P., Lee, J.J., Yeo, H., Kim, C.K., Pang, S.J., and Oh, J.K. (2016). Parametric study on the capability of three-dimensional finite element analysis (3D-FEA) of compressive behaviour of Douglas fir. Holzforschung 70: 539–546, https://doi.org/10.1515/hf-2015-0151.Search in Google Scholar

Kollmann, F.F. and Côte, W.A. (1968). Principles of wood science and technology I. Solid wood. In: Principles of wood science and technology. New York: Springer-Verlag.10.1007/978-3-642-87928-9Search in Google Scholar

Kumpenza, C., Matz, P., Halbauer, P., Grabner, M., Steiner, G., Feist, F., and Müller, U. (2018). Measuring Poisson’s ratio: mechanical characterization of spruce wood by means of non-contact optical gauging techniques. Wood Sci. Technol. 52: 1451–1471, https://doi.org/10.1007/s00226-018-1045-7.Search in Google Scholar

Lukacevic, M., Lederer, W., and Füssl, J. (2017). A microstructure-based multisurface failure criterion for the description of brittle and ductile failure mechanisms of clear wood. Eng. Fract. Mech. 176: 83–99, https://doi.org/10.1016/j.engfracmech.2017.02.020.Search in Google Scholar

Mackenzie-Helnwein, P., Eberhardsteiner, J., and Mang, H.A. (2003). A multi-surface plasticity model for clear wood and its application to the finite element analysis of structural details. Comput. Mech. 31: 204–218, https://doi.org/10.1007/s00466-003-0423-6.Search in Google Scholar

Mackenzie-Helnwein, P., Müllner, H.W., Eberhardsteiner, J., and Mang, H.A. (2005). Analysis of layered wooden shells using an orthotropic elasto-plastic model for multi-axial loading of clear spruce wood. Comput. Methods Appl. Mech. Eng. 194: 2661–2685, https://doi.org/10.1016/j.cma.2004.07.051.Search in Google Scholar

Marcolin, L.A., Moritani, F.Y., Rodegheri, P.M., and Rocco Lahr, F.A. (2021). Properties relationship evaluation and plasticity analytical model approach for Brazilian tropical species. Eur. J. Wood Prod. 79: 477–485, https://doi.org/10.1007/s00107-020-01605-x.Search in Google Scholar

Milch, J., Tippner, J., Sebera, V., and Brabec, M. (2016). Determination of the elasto-plastic material characteristics of Norway spruce and European beech wood by experimental and numerical analyses. Holzforschung 70: 1081–1092, https://doi.org/10.1515/hf-2015-0267.Search in Google Scholar

Milch, J., Brabec, M., Sebera, V., and Tippner, J. (2017). Verification of the elastic material characteristics of Norway spruce and European beech in the field of shear behaviour by means of digital image correlation (DIC) for finite element analysis (FEA). Holzforschung 71: 405–414, https://doi.org/10.1515/hf-2016-0170.Search in Google Scholar

Moses, D.M. and Prion, H.G. (2004). Stress and failure analysis of wood composites: a new model. Compos. B Eng. 35: 251–261, https://doi.org/10.1016/j.compositesb.2003.10.002.Search in Google Scholar

Müller, U., Jost, T., Kurzböck, C., Stadlmann, A., Wagner, W., Kirschbichler, S., Baumann, G., Pramreiter, M., and Feist, F. (2020). Crash simulation of wood and composite wood for future automotive engineering. Wood Mater. Sci. Eng. 15: 312–324, https://doi.org/10.1080/17480272.2019.1665581.Search in Google Scholar

Nairn, J.A. (2006). Numerical simulations of transverse compression and densification in Wood. Wood Fiber Sci. 38: 576–591.Search in Google Scholar

Oudjene, M. and Khelifa, M. (2009). Elasto-plastic constitutive law for wood behaviour under compressive loadings. Construct. Build. Mater. 23: 3359–3366, https://doi.org/10.1016/j.conbuildmat.2009.06.034.Search in Google Scholar

Pambou Nziengui, C.F., Ikogou, S., and Moutou Pitti, R. (2018). Impact of cyclic compressive loading and moisture content on the mechanical behavior of Aucoumea klaineana Pierre. Wood Mater. Sci. Eng. 13: 190–196, https://doi.org/10.1080/17480272.2017.1307281.Search in Google Scholar

Patton-Mallory, M., Cramer, S.M., Smith, F.W., and Pellicane, P.J. (1997). Nonlinear material models for analysis of bolted wood connections. J. Struct. Eng. 123: 1063–1070, https://doi.org/10.1061/(asce)0733-9445(1997)123:8(1063).10.1061/(ASCE)0733-9445(1997)123:8(1063)Search in Google Scholar

Pěnčík, J. (2015). Tests of wooden specimens from Scots pine (Pinus sylvestris) with the help of anisotropic plasticity material model. Drvna Ind 66: 27–33, https://doi.org/10.5552/drind.2015.1362.Search in Google Scholar

Qing, H. and Mishnaevsky, L. (2011). A 3D multilevel model of damage and strength of wood: analysis of microstructural effects. Mech. Mater. 43: 487–495, https://doi.org/10.1016/j.mechmat.2011.05.007.Search in Google Scholar

Rahardjo, H., Harnas, F.R., Indrawan, I.G.B., Leong, E.C., Tan, P.Y., Fong, Y.K., and Ow, L.F. (2014). Understanding the stability of Samanea saman trees through tree pulling, analytical calculations and numerical models. Urban For. Urban Green. 13: 355–364, https://doi.org/10.1016/j.ufug.2013.12.002.Search in Google Scholar

Scheperboer, I.C., Suiker, A.S.J., Luimes, R.A., Bosco, E., and Jorissen, A.J.M. (2019). Collapse response of two-dimensional cellular solids by plasticity and cracking: application to wood. Int. J. Fract. 219: 221–244, https://doi.org/10.1007/s10704-019-00392-8.Search in Google Scholar

Schmidt, J. and Kaliske, M. (2006). Zur dreidimensionalen Materialmodellierung von Fichtenholz mittels eines Mehrflächen-Plastizitätsmodells. Holz Roh-Werkst. 64: 393–402, https://doi.org/10.1007/s00107-006-0102-3.Search in Google Scholar

Šebek, F., Kubik, P., Brabec, M., and Tippner, J. (2020). Modelling of impact behaviour of European beech subjected to split Hopkinson pressure bar test. Compos. Struct. 245: 112330, https://doi.org/10.1016/j.compstruct.2020.112330.Search in Google Scholar

Šebek, F., Kubik, P., Tippner, J., and Brabec, M. (2021). Orthotropic elastic-plastic-damage model of beech wood based on split Hopkinson pressure and tensile bar experiments. Int. J. Impact Eng. 157: 103975, https://doi.org/10.1016/j.ijimpeng.2021.103975.Search in Google Scholar

Sirumbal-Zapata, L.F., Málaga-Chuquitaype, C., and Elghazouli, A.Y. (2018). A three-dimensional plasticity-damage constitutive Mmodel for timber under cyclic loads. Comput. Struct. 195: 47–63, https://doi.org/10.1016/j.compstruc.2017.09.010.Search in Google Scholar

Tabiei, A. and Wu, J. (2000). Three-dimensional nonlinear orthotropic finite element material model for wood. Compos. Struct. 50: 143–149, https://doi.org/10.1016/s0263-8223(00)00089-1.Search in Google Scholar

Tsai, S.W. and Wu, E.M. (1971). A general theory of strength for anisotropic materials. J. Compos. Mater. 5: 58–80, https://doi.org/10.1177/002199837100500106.Search in Google Scholar

Tukiainen, P. and Hughes, M. (2016). The effect of temperature and moisture content on the fracture behaviour of spruce and birch. Holzforschung 70: 369–376, https://doi.org/10.1515/hf-2015-0017.Search in Google Scholar

Vojáčková, B., Tippner, J., Horáček, P., Praus, L., Sebera, V., and Brabec, M. (2019). Numerical analysis of branch mechanical response to loading. Arboric. Urban For. 45: 120–131, https://doi.org/10.48044/jauf.2019.011.Search in Google Scholar

Yang, M., Défossez, P., Danjon, F., and Fourcaud, T. (2014). Tree stability under wind: simulating uprooting with root breakage using a finite element method. Ann. Bot. 114: 695–709, https://doi.org/10.1093/aob/mcu122.Search in Google Scholar PubMed PubMed Central

Yoshihara, H. (2014). Plasticity analysis of the strain in the tangential direction of solid wood subjected to compression load in the longitudinal direction. Bioresources 9: 1097–1110, https://doi.org/10.15376/biores.9.1.1097-1110.Search in Google Scholar

Yoshihara, H. and Ohta, M. (1992). Stress-strain relationship of wood in the plastic region. I. Examination of the applicability of plasticity theories. Mokuzai Gakkaishi 38: 759–763.Search in Google Scholar

Yoshihara, H. and Ohta, M. (1994). Stress-strain relationship of wood in the plastic region. II. Formulation of the equivalent stress-equivalent plastic strain relationship. Mokuzai Gakkaishi 40: 263–267.Search in Google Scholar

Received: 2022-01-16
Accepted: 2022-07-13
Published Online: 2022-08-04

© 2022 the author(s), published by De Gruyter, Berlin/Boston

This work is licensed under the Creative Commons Attribution 4.0 International License.

Scroll Up Arrow