Micromechanical modelling of creep in wooden materials

(1)

MAT-VET-F 21006

Examensarbete 15 hp June 2021

Micromechanical modelling of creep in wooden materials

Kevin Coleman

Oskar Falkeström

Malin Nilsson

(2)

Teknisk- naturvetenskaplig fakultet UTH-enheten

Besöksadress:

Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0

Postadress:

Box 536 751 21 Uppsala

Telefon:

018 – 471 30 03

Telefax:

018 – 471 30 00

Hemsida:

http://www.teknat.uu.se/student

Abstract

Micromechanical modelling of creep in wooden materials

Kevin Coleman, Oskar Falkeström, Malin Nilsson

Wood is a complex organic orthotropic viscoelastic material with a cellular structure. When stressed, wood will deform over time through a process called creep. Creep affects all wooden structure and can be difficult, time-consuming and expensive to measure.

For this thesis, a simple computer model of the wooden

microstructure was developed. The hypothesis was that the modelled microstructure would display similar elastic and viscoelastic

properties as the macroscopic material.

The model was designed by finding research with cell geometries of coniferous trees measured. The model considered late- and earlywood geometries as well as growth rings. Rays were ignored as they only composed 5-10% of the material.

By applying a finite element method, the heterogeneous late- and earlywood cells could be homogenized by sequentially loading the strain vector and calculating the average stress.

The computer model produced stiff but acceptable values for the elastic properties. Using the standard linear solid method to model viscoelasticity, the computer model assembled creep curves comparable to experimental results.

With the model sufficiently validated, parametric studies on the cell geometry showed that the elastic and viscoelastic properties changed greatly with cell shape. An unconventional RVE was also tested and shown to give identical result to the standard RVE.

Although not perfect, the model can to a certain degree predict the elastic and viscoelastic characteristics for wood given its

cellular geometry. Inaccuracies were thought to be caused by assumptions and approximations when building the model.

(3)

Populärvetenskaplig sammanfattning

Den 14 juli 1902 kollapsade San Marcos klocktorn i Venedig, 729 år efter att den byggdes.

Detta var en stor tragedi för Venedigs befolkning som beslöt sig för att återuppbygga tornet på precis samma sätt som innan. Klocktornet i San Marco är därmed dömd att i framtiden kollapsa igen [1]. I Sverige bedrivs forsking för att något liknande inte ska hända den svenska nationalklenoden Vasaskeppet. Vasaskeppet, likt San Marcos klocktorn, lider av något som kallas för krypning. Krypning är en konstant försvagning av material under lång tid. Jämför med hur en plastpåse långsamt kan slitas sönder av innehållets vikt. På samma sätt kryper alla material med tiden. Frågan är alltså inte om krypning går att helt och hållet stoppa, utan till vilken grad det kan motverkas.

Detta är inte bara viktigt för redan existerande strukturer som Vasaskeppet, utan även för framtida konstruktioner. För vad är meningen med dagens fantastiska konstruktioner och arkitekturer om inte eftervärlden får ta del av dem? Denna studie försöker besvara frågan hur krypning kan motverkas i barrved som material.

Trä som byggnadsmaterial har varit viktig för mänskligheten sedan stenåldern. Eftersom trä både är hållfast och förnybart så är det fortfarande relevant idag. Trä är ett så kallat viskoelastiskt material. Detta innebär att de egenskaper som bestämmer hur materialet töjs kan förändras med tiden. Denna egenskap är själva grunden till varför krypning sker i trämaterial. Linjärelastisitet å andra sidan, är ett begrepp som innebär att de krafter som påverkar materialet är kopplade till materialets töjning på samma sätt hela tiden.

Trä kan beskrivas på olika sätt. På den yttersta, så kallade makroskopiska nivån, består träd av årsringar. De ljusa partierna av årsringarna växer tidigt på säsongen och kallas därmed för vårved. De mörka partierna växer senare på säsongen och kallas för sommarved. Studien bygger på hypotesen att trämaterialets beståndsdelar på en mikroskopisk nivå, i form av träceller, kan beskriva och kontrollera hur det makroskopiska trämateri- alet beter sig. Ett av dessa beteenden var till vilken grad trämaterialet kan motstå att deformeras, så kallad styvhet.

Figure 1: Honungsbin som bygger en vaxkaka av bivax [2]

På ännu mindre skala består träcellerna av långa fibrer som sträcker sig genom hela trädet. Tillsammans bildar fibrerna ett stort nät, likt en vaxkaka skapad av bin, se figur 1. Det är dessa nät av fibrer som ger trädet sin styvhet vilket möjliggör för trädet att stå upp. Träcellerna i sig består av så kallade mikrofibriller av cellulosa. Träcellerna

(4)

har egenskaper såsom tjocklek på cellväggen, bredd/höjd och mikrofibrillernas vinkel i fibercellerna. Dessa egenskaper skiljer sig åt mellan olika trädtyper. Hypotesen testades genom att en datormodell byggdes upp över en träcell. Datormodellen byggdes upp som den minsta möjliga beståndsdelen av denna cellstruktur som kan beskriva hela materialet.

Först tillverkades en modell som beskriver trädet som linjärelastiskt, eftersom detta är det simplaste fallet. Nämnda egenskaper hos träcellen, samt mängden vårved/sommarved ändrades sedan i denna datormodell för att jämföras med kända värden från olika typer av träd. På sätt kunde modellen verifieras. Efter detta utfördes en så kallad parameterstudie där den linjärelastiska träcellens olika egenskaper ändrades för att se hur materialet i stort reagerade på detta. Det kunde då observeras att exempelvis mikrofibrillvinkeln, mängden sommarved och cellens relation mellan bredd och höjd påverkade trädets styvhet.

Sedan användes de tidigare resultaten för att bygga en viskoelastisk datormodell som därmed tog med hur trädet reagerar på krafter över tid. En parameterstudie utfördes även för att se huruvida krypning förbättras eller förvärras med ändringar hos träcellens egenskaper. Det kunde exempelvis konstateras att krypning i en riktning kan motverkas genom att minska mikrofibrillens vinkel. Om cellens bredd eller höjd ökas, minskar krypning i respektive riktning. Generellt kunde det visas att ju styvare materialet blir, desto mer kan materialet motverka krypning.

Även om denna studie kom fram till användbara resultat efterfrågas vidare forskning.

Detta för att kunna motverka ytterligare försvagning av äldre konstruktioner såsom Vasaskeppet, samt för att motverka krypning av framtida byggnationer.

(5)

Acknowledgements

We are grateful for being offered to take part in this research group and would like to express our gratitude to individuals that made this thesis project possible.

A special thank you to PhD student Rhodel Bengtsson for supervising this project and for his contributions in supplying relevant books and compendiums, lectures, weekly check-up meetings and input on the report. We thank you and wish you all the best with your continued work towards your PhD.

Furthermore, we would like to thank Professor Kristofer Gamstedt for including us in the project, warmly welcoming us to the research group and allowing us to present our results to the Applied Mechanics division at Uppsala University.

Moreover, we want to show our gratitude to associate Professor Mahmoud Mousavi and Researcher Reza Afshar giving us guidance and valuable input throughout the project.

We also want to thank our project mentor, assistant Professor Ocean Cheung, for answering course-related questions and continuously reviewing our progress.

Last but not least, we would like to thank the course coordinators at the Department of Materials Science and Engineering at Uppsala University: Peter Svedlindh, Tomas Edvinsson, Germán Salazar Alvarez, Martin Sjödin and Ken Welch for enabling our cross-program thesis work.

Kevin Coleman Oskar Falkeström Malin Nilsson

(6)

1 Introduction

Wood as a material has been used by humanity since the stone age [3]. Its abundance and versatility has made wood a good building material commonly used all throughout human civilization. Wooden artefacts at least 5000 years old have been found in the pyramids. The word material is itself derived from the Latin word for lumber, materia.

Compared to many other popular building materials, wood is also a renewable resource [3].

Although prevalent, wood is actually a highly complex material. The structure of wood looks very different at different sizes, from small wood cells to the planks used when building. Being an organic material, the mechanical properties of wood can vary between trees within the same species depending on growth conditions. Furthermore, the strength between different tree species can be significant. For example, oak can be more than 10 times stronger than balsa [4]. Accurately describing the mechanical properties of "wood"

can therefore be tricky and, often, requires extensive experimentation.

Wood is also a, so called, viscous material. This means that wood will deform with time.

A wooden structure will slowly sag until it collapses. This deformation, called creep, is also a complex phenomena. Understanding creep is necessary to prevent damage to wooden structures or buildings. A relevant example is the ship of Vasa, a Swedish 17th century warship. Vasa capsized in the Stockholm port during her maiden voyage [5]. The ship is now on display at the Vasa Museum. The ship is however suffering from creep deformation. Research on the ship’s deformation has been done by gathering plenty of data and measurements on the ship [6]. It is commonplace that conducting creep experiments are difficult, time consuming and expensive.

Figure 2: Sagging wooden house close to collapsing.

1.1 Objective

The underlying idea behind this thesis work was to develop and test a simple mechanical model of wood. Since wood is a cellular material the idea was to model wood at the cellular level. The hypothesis was that the mechanical properties of the cellular material would be similar to the wooden material. The differences between tree species could then partly be described by their differing cell geometries.

Accurately modelling the cell geometries required reading and learning of: Mechanics of materials, numerical methods and wood cell structures. The basic principle was to find

(8)

species with a measured cell geometry, mechanical properties and creep behaviour. Then for the model, numerically computing the mechanics of materials equations, extracting values, and comparing them to the measured properties.

The project was divided into two stages. The model would first be tested for the station- ary case to ensure that the linear elastic mechanical properties were sufficiently close.

Subsequently the time dependence was added to investigate the linear viscoelastic creep behaviour.

1.2 Goals

The goal was to develop a model that gave reasonable predictions for a tree species elastic and viscoelastic behaviour given its cell geometry. A prediction was considered reasonable if the elastic properties were similar to the experimental value. Furthermore, the elastic properties should preferably show reasonable relation to each other e.g. the model should accurately rank the stiffness of the three orthogonal direction. Simulated creep curves were considered acceptable if they looked similar to the experimental curves.

With the model established, the aim was to conduct parameter studies on the cell geometry and see how the elastic and creep behaviours were affected. The goal was to understand which mechanical properties were affected by the geometric parameter, and to what extent.

The project was considered complete when the following had been achieved:

• A working model with implemented linear- elastic and viscoelastic behaviour

• A linear elastic parameter study

• A linear viscoelastic parameter study

(9)

2 Theory

2.1 Mechanics of Materials

Mechanics of materials is a field which aims to understand and calculate deformations and failure of loaded structures. Mechanics of materials is an important field in engineering as the knowledge of material strengths and behaviours is vital when constructing any type of structure. This chapter will present the basic principles of mechanics of materials used in this thesis project.

As introduction, the words and concepts of strain, stress and displacement need to be understood. Stress, denoted σ, is the force applied over an area. Stress and applied force are thus highly interconnected. However, stresses are often preferable to forces since stresses are area independent. An equal force applied on different sized surfaces yields a highly varied deformation. The stress-force relation is described in section 2.1.1.

The displacement of a particle in a material is intuitively the distance it is moved during deformation. Displacement is often represented as a vector field u where the components of u equal the dimensions of possible deformation (in 3D denoted u, v & w). Each point {x1,x2,x3} on the undeformed body is mapped on the deformed body through the vector field u. This vector field is shown in equation 1.

upu, v, wq “

»

— –

upx1, x₂, x₃q vpx1, x₂, x₃q wpx1, x₂, x₃q

fi ffi

fl (1)

Strains, denoted εij, describe the relative deformation of the material. One unit of strain would therefore mean a doubling of length, or equivalently that the deformed object is twice the length of the undeformed object. The two indices represent the surface normal direction i and the strain direction j. Commonly a strain matrix is built with all nine three dimensional strain components, equation 2.

ε “

»

— –

ε₁₁ ε₁₂ ε₁₃ ε₂₁ ε₂₂ ε₂₃ ε₃₁ ε₃₂ ε₃₃

fi ffi

fl (2)

Stress, strain and displacement are, as one could expect, not independent from each other.

While the stress-strain relation is discussed in section 2.1.2 the strain-displacement relation can be derived through a theoretical deformation of a line segment with undeformed length δl0. Two cases are imagined: one with uniaxial strain, figure 3a, and one with shear strain, figure 3b.

(10)

x x ` δlo

upx, x2, x3q

upx ` lo, x2, x3q

x1

x2

x3

du e fin ă 3

(a) Uniaxial strain of line element.

x x ` δlo

vpx, x2, x3q

vpx ` lo, x2, x3q

x1

x2

x3

γ

δl₀ ` δl0 Bu Bx1

δl₀ ^Bv Bx1

(b) Shear strain of line element lying in x1

direction.

Figure 3: Deformations of line element with undeformed length δl₀.

The engineer strain of the line element for uniaxal strain is defined as the unit length increase. Shear strain is calculated by the angle created from a shear stress. As shear stress are two-dimensional in nature, for a multidimensional body one must also account for a line element lying in a separate direction and its sheer angle. The results of the two elements are however analogous. The shear strain is then half of the two angles γ.

ε₁₁“ upx ` δl0, x₂, x₃q ´ upx, x2, x₃q δl₀

δl0Ñ0

ùùùùñ Bu Bx1

(3)

γ_x₁ « tanpγx1q “

Bv Bx1

1 ` _Bx^Bu

1

« Bv Bx1

ε₁₂“ 1

2`γ_x₁ ` γx2

˘

(4)

Performing the same calculation for all components in the ε-matrix yields component expressions for small deformations seen in equation 5. Note that the matrix is symmetric.

ε₁₁“ Bu Bx1

, ε₂₂ “ Bv Bx2

, ε₃₃ “ Bw Bx3

ε12“ 1 2

´Bu Bx2

` Bv Bx1

¯

, ε13 “ 1 2

´Bu Bx3

` Bw Bx1

¯

, ε23“ 1 2

´ Bv Bx3

` Bw Bx2

¯ (5)

Combining strain and displacement notation from equations 1 & 2 a matrix representation of equation 5 can be written. This is shown in equation 6.

ε “ 1

2p∇u ` ∇u^|q (6)

2.1.1 Cauchy’s Equilibrium Equation

The central governing equation when modelling a 3D body, as done in this thesis, is the equilibrium equation. The underlining principles for the equation is Newton’s second law of motion, that the total force acting on a body is equal to the change of momentum [7].

(11)

volume forces (or body forces) that arise within the body. These are forces per volume that, for example, occur due to gravity and will in this thesis be denoted fv. Then there are the boundary forces acting on the surface of the body. For a three dimensional body, each boundary can have forces with components in three directions. Recalling from section 2.1, force is preferably described with force per area (stress) and the convenient way to represent this is with the stress tensor. Coupled with a surface’s normal vector, the force per area on that surface is received. A representation can be seen in figure 4 and equation 7.

σ₂₁

σ22

σ₂₃ σ₃₁

σ₃₂ σ₃₃

σ₁₁

σ₁₂ σ₁₃

x₁ x2

x₃

Figure 4: All possible stress components acting on a three dimensional body.

σ ¨ n ”

»

— –

σ₁₁ σ₁₂ σ₁₃ σ₂₁ σ₂₂ σ₂₃ σ₃₁ σ₃₂ σ₃₃

fi ffi fl ¨

»

— –

n₁ n₂ n₃

fi ffi

fl (7)

The arbitrary body can be divided into cubes with small volume δV and surface δS. The total force on this body is gained by summing all the cubes after multiplying the volume forces and equation 7 by δV and δS respectively. Taking the limit yields:

ÿ

BV

pσ ¨ nqδS ` ÿ

V

fvδV “ 0ùùùùùñ^{δS,δV Ñ0} ż

BV

σ ¨ n dS ` ż

V

fv dV “ 0 (8) Equation 8 is then described as a volume integral by applying the divergence theorem.

ż

BV

σ ¨ n dS ` ż

V

f_v dV “ ż

V

∇ ¨ σ dV ` ż

V

f_v dV “ 0 ùñ

ż

V

∇ ¨ σ ` fv dV “ 0

ùñ ∇ ¨ σ ` fv “ 0 (9)

Equation 9 is the equilibrium equation for a resting body [8]. The equation consists of three partial differential equation and is sometimes written as a set of equation rather

(12)

Bσ11

Bx1

`Bσ12

Bx2

`Bσ13

Bx3

` f_v¹ “ 0 Bσ21

Bx1

`Bσ22

Bx2

`Bσ23

Bx3

` f_v² “ 0 Bσ31

Bx1

`Bσ32

Bx2

`Bσ33

Bx3

` f_v³ “ 0

2.1.2 Linear Elasticity and Hooke’s Law

An elastic material is by definition capable to restore itself to its original shape after an applied stress is removed. If the strains in a material are small and proportional to the stress, a linear relationship can be described. These are the conditions of the ’Linear Elastic model’. Linear elastic dependence of stress and strain is illustrated in figure 5.

These conditions hold for a certain material until the elastic limit of that material is reached.

σ

ε

unload load

Figure 5: The graph illustrates a linear elastic dependence.

There are different material definitions based on the intrinsic symmetries of the material properties. For example, an isotropic material has the same material properties in every direction. An anisotropic material on the other hand has no directional dependence and no symmetry. Materials such as wood and other composites are often part of a third material definition called orthotropic materials. These materials have different properties in its perpendicular directions. Because of the three orthogonal planes in the orthotropic material, it will have a symmetry that can reduce the many independent elastic constants to only nine constants. Consequently, a shear stress (index 4,5 and 6) in one direction will act the same but negatively in the opposite direction. The relations can be written as:

σ₁ “ C1xε_x, σ₂ “ C2xε_x, σ₃ “ C3xε_x

σ₄ “ C4xε_x, σ₅ “ C5xε_x, σ₆ “ C6xε_x (10) In the opposite direction an induced strain, -ε , creates the following stresses:

(13)

Furthermore, to simplify the complete stiffness matrix the following symmetries can be used:

ε₁ ñ σ1 “ σ¹₁ ε2 ñ σ2 “ σ¹₂ ε₃ ñ σ3 “ σ¹₃

(12)

ε₄ ñ σ4 “ ´σ₄¹ , C₁₄ “ ... “ C34 “ 0 ε₅ ñ σ5 “ ´σ₅¹ , C₁₅ “ ... “ C45 “ 0 ε₆ ñ σ6 “ ´σ₆¹ , C₁₆ “ ... “ C56 “ 0

The stiffness matrix can thus be written:

»

—

— –

σ₁ σ₂ σ₃ σ₄ σ₅ σ₆

fi ffi ffi ffi ffi ffi ffi ffi ffi fl

“

»

—

— –

C₁₁ C₁₂ C₁₃ 0 0 0 C₂₂ C₂₃ 0 0 0

C₃₃ 0 0 0

C₄₄ 0 0 C₅₅ 0 C₆₆

»

—

— –

ε₁ ε₂ ε₃ ε₄ ε₅ ε₆

(13)

Equation 13 describes the relationship between strain and stress. This relationship is known as Hooke’s Law, and can also be written on the following form.

σ “ C ¨ ε (14)

Where C is a material constant called stiffness. This law can be applied to solve various homogeneous problems involving these two properties. By knowing some of the relations, the rest can be evaluated.

Inverting the matrix in equation 13 yields the compliance matrix. The elements in the orthogonal compliance matrix are known and shown in equation 15. Note that the indices 4,5 & 6 have been replaced with indices representing the plane shear strain or stress. The indices 1,2 & 3 thus represents the three orthogonal directions. Shear strains/stresses are thereby defined by two directions and uniaxial strains/stresses by one, as shown in their respective indices.

»

—

— –

ε₁ ε₂ ε₃ 2ε₁₂ 2ε₁₃ 2ε

“

»

—

— –

1

E1 ´^ν_E²¹

2 ´^ν_E³¹

3 0 0 0

´^ν_E¹²₁ _E¹₂ ´^ν_E³²₃ 0 0 0

´^ν_E¹³₁ ´^ν_E²³₂ _E¹₃ 0 0 0

1

G12 0 0

1

G13 0

1

»

—

— –

σ₁ σ₂ σ₃ σ₁₂ σ₁₃ σ

(15)

(14)

In Equation 15, three new parameters are now obtained, shear modulus Gij, elastic modulus Eij and Poisson’s ratio νij [9]. The modulus of elasticity, or Young’s modulus, measures the tensile stiffness of a solid material, where a higher value makes the material stiffer. If the material is affected by shear stresses and it is proportional to the shear strain, it will result in a slope which is equivalent to the shear modulus for the elastic region in a stress-strain graph. Poisson’s ratio for elastically anisotropic materials works in a different way compared to isotropic materials. When a material is being exposed to a tensile stress, it gives rise to an elongation. This elongation will in turn give rise to constrictions in the lateral directions that is perpendicular to an applied stress. For isotropic materials that have an applied stress in its z-direction, the parameter is described as the ratio between the lateral and axial strains. But for an elastically anisotropic material, it is the crystallographic direction that is the determining factor for the elastic behavior [10].

2.1.3 Linear viscoelasticity and the standard linear solid model

The elastic description of a material discussed in the previous section models the stress- strain relation as linearly proportional. Another possible material description is a viscous model. Viscosity η, more commonly used when describing fluid or gas, is the internal resistance of flowing [11] Concretely, the stress-strain relationship of viscous materials linearly relates the stress to the time derivative of the strain η 9ε “ σ. Loading and unloading might therefore yield different strains depending on how the stress was applied, this is seen in figure 6.

σ

ε

Hysteresis loop

load

unload

Figure 6: Illustration of a viscoelastic dependence.

Even though solid materials do not behave like fluids, some materials tend to show time-dependent behaviours. This would suggest that a degree of viscosity exists in said material. A suitable choice is to combine the elastic and viscous material model to a viscoelastic model. A viscoelastic model displays both elastic and viscous properties.

A common model is the standard linear solid model (SLS). This model combines linear elastic springs and linear viscous dashpots (a dampener), forming a linear viscoelastic model.

(15)

E2

E₁

η₂

σ σ

Figure 7: Stanard linear solid model for viscoelastic materials.

The lower branch in figure 7 is called a Maxwell branch. Adding further Maxwell branches in parallel gives rise to the more advanced generalized Maxwell model. In this thesis however, the simpler SLS model will suffice. The combined stress-strain relation of one Maxwell branch (figure 7 without E1) is easily derived from Hooke’s law and the viscous behaviour of a dashpot, shown in equation 16. The derivation uses the Laplace transform [12] as it turns the set of differential equations to a set of algebraic equations. Note that in equation 16, ε1 and ε2 denoted the strain over the spring and dashpot respectively.

$

’&

’%

ε_tot “ ε1` ε2

ε₁ “ _E¹σ ε9₂ “ ¹_ησ

Laplace

ùùùùùùùñ

T ransf orm

$

’&

’%

˜

ε_tot “ ˜ε₁ ` ˜ε₂

˜

ε₁ “ _E¹σ˜ s˜ε₂ “ _η¹σ˜ ñ ˜σ ` η

Es˜σ “ ηs˜ε_tot ù^Inverseùùùùñ

Laplace σ ` η

Eσ “ η 9ε9 tot

(16)

Deriving the stress-strain relation for the entire model in figure 7 can be done in a similar way. Some relations for the stress/strain can be determined from the figure.

Firstly, the strain on both the upper and lower branch has to equal each other along with both equalling the strain on the entire model. Secondly, the stress from both branches additively form the total stress. Stress-strain relations follows earlier discussion and the relation for the Maxwell branch is seen in equation 16.

(16)

$

’’

’&

’’

’%

σ_tot “ σ1` σ2

ε_tot “ ε1 “ ε2

σ₁ “ E1ε₁ σ₂`_E^η²

2σ9₂ “ η2ε9

Laplace

ùùùùùùñ

Transform

$

’’

’&

’’

’%

˜

σ_tot “ ˜σ₁` ˜σ₂

˜

ε_tot “ ˜ε₁ “ ˜ε₂

˜

σ₁ “ E1ε˜₁

˜ σ₂` _E^η²

2s˜σ₂ “ η2s˜ε₂ ñ ˜σ_tot “ E1ε˜_tot` η2

1 ` _E^η²

2sε˜_tot ñ σ˜_tot` η2

E₂s˜σ_tot “ E1`1 ` η2

E₂s˘ ˜ε_tot` η2s˜ε_tot ñ σ `˜ η₂

E₂s˜σ “ E1ε˜_tot`η₂pE1` E2q E₂ s˜ε_tot

Inverse

ùùùùùñ

Laplace σ_tot` η2

E₂σ9_tot “ E1ε_tot`ηpE1` E2q

E₂ ε9_tot (17)

2.1.4 Creep

Consequently from the linear viscoelastic model, solving the stress-strain constitutive equation derived in equation 17 requires solving a time-dependent differential equation.

This means that the strain will be time-dependent. When loading a linear viscoelastic material there will first be an instantaneous elastic strain followed by a time-dependent viscous strain. A viscoelastic structure experiencing a continuous stress will thus contin- ually deform over time. This time-dependent strain is called creep.

Creep usually follow a specific curve, a creep curve. In Figure 8 the three different stages of a creep curve can be seen. The three stages are the primary, secondary and the tertiary creep where rapture occurs immediately after. A creep tests is a tensile test with constant load. The primary creep shows the elastic region and in the secondary creep (which is called steady state) the creep is nearly constant. This stage takes the longest time. The last and tertiary creep leads to ductile fracture and occurs before the rupture where the graph ends [13].

Often, due to the difficulty of simulating tertiary creep, only the primary and secondary stages creep are considered. This will be done in this thesis project.

(17)

ε

t ptimeq

Primary creep

Secondary creep

Tertiary creep

Rupture

Instantaneous deformation

Figure 8: Illustration of a creep curve and its three stages.

A viscoelastic structure under a fixed strain will over time have its material stress reduced.

This process is called relaxation. Defining the relaxation modulus Eptq and the creep compliance Jptq for a constant stress σ0 and strain ε0 respectively can be seen in equation 18. For a constant unit strain and stress the relaxation modulus and the creep compliance are each others inverse.

εptq “ σ0J ptq, σptq “ ε0Eptq (18) One way of understanding viscoelastic properties is by doing a relaxation test. A constant strain is applied and from the measured, or calculated, relaxation curve the relaxation modulus can be extracted. For the SLS viscoelastic model derived i section 2.1.3 the relaxation modulus can be calculated from the stress-strain relation in equation 17 when the applied strain is constant, εptq “ ε0.

εptq “ ε0 ùñ 9εptq “ 0 (19)

With the condition presented in equation 19 the SLS stress-strain differential equation, equation 17, becomes a standard linear ordinary differential equation and is easily solved with a general solution [14].

σ `9 E₂

η₂σ “ E₁E₂ η₂ ε₀ ùñ σptq “ e^´

ş_E2

η2dt` ż

e^E2^η2^tE₁E₂ η2

ε₀ dt ` C˘

ùñ σptq “ e^´

E2

η2tÈ₁ε₀eÊ2^η2^t` ε0C˘ ùñ σptq “È₁` Ce^´

E2 η2t˘ loooooooomoooooooon

Eptq

ε₀

(20)

(18)

Solutions to the differential equation can be seen in equation 20. Furthermore, at t “ 0 the material should display only elastic characteristics. Due to this requirement, the sum of E1 and C needs to equal the material’s elastic properties. This is shown in equation 21. Luckily, linear elastic properties are thoroughly discussed in section 2.1.2.

Ept “ 0q “ E1` C “ Celastic (21)

Going from one dimensional linear viscoelastic material to a three dimension linear viscoelastic material simply involves replacing the scalar parameters E1 & C with matrices so that they together form the linear elastic matrix in equation 13. Also, conventionally the exponent _E^η²₂ is replace by a scalar value τ. This new scalar parameter is called the relaxation time as it governs the rate at which the material relaxes.

(19)

2.2 Numerical methods for mechanics of materials

2.2.1 Finite element formulation

The finite element method, or FEM, is a method used in approximating solutions to partial differential equations (PDE). FEM is a powerful numerical tool in many different fields and mechanics of materials is no exception.

The basic idea behind a FEM approximation is to discretize the physical domain into small elements with nodes connecting the elements to each other. The elements’ shape and number of nodes can be varied depending on the geometry of the physical domain.

This discretization process is called meshing and generates a mesh of elements and nodes.

In each element, the PDE solution is approximated in every node with a chosen set of shape functions. These shape functions are often simple polynomial which are easy to work with. Commonly these shape functions are piecewise functions and the requirement on differentiability imposed by the PDE must be softened. This is done by turning the PDE to a, so called, weak form. This is usually done by integrating the PDE over the domain and is thereby sometimes called "integral form".

n₁ n₂

n₃

n₄ K

Ω

Figure 9: Physical domain and example tetrahedron element with surrounding nodes.

The governing PDE is the equilibrium equation, introduced in section 2.1.1 equation 9.

The weak formulation of this equation is derived by integrating over the domain after multiplying the equilibrium equation by a test function v in a suitable function space.

Depending on the specific boundary conditions (BC), what constitutes as a suitable function space may vary. The equilibrium equation multiplied by the test function and the defined function space V is seen in equation 22. Translating the mathematical language used when defining the vector space, V is the space of all functions v. The functions v have three components, each part of a Hilbert room defined on the domain Ω. All functions v have to be bounded and follow any possible BC on the domain’s boundary.

p∇ ¨ σqv ` fVv “ 0

v P V “ tv P rHpΩqs³ : ||v||² ` ||∇v||² ă 8 : v|BΩ“BC u (22)

(20)

Integrating over the domain and applying the divergence theorem:

ż

Ω

p∇ ¨ σqv dV ` ż

Ω

f_vv dV “ 0

ùñ ż

Ω

p

3

ÿ

j,i“1

pBσij

Bxj

qvi dV ` ż

Ω

f_vv dV “ 0

ùñ

3

ÿ

j,i“1

ż

BΩ

pσij ¨ njqvi dS ´

3

ÿ

j,i“1

ż

Ω

σ_ij Bvi

Bxj

` ż

Ω

f_vv dV “ 0

ùñ

3

ÿ

j,i“1

ż

Ω

σ_ijBvi

Bxj

“ ż

BΩ

pσ ¨ nqv dS looooooomooooooon

Boundary term

` ż

Ω

f_vv dV

The left term needs further work. The contraction operation is thus introduced, seen in equation 23. The PDE simplifies as follows:

A : B ”

3

ÿ

i,j“1

A_ijB_ij (23)

ż

Ω

σ : ∇v “ ż

BΩ

pσ ¨ nqv dS ` ż

Ω

f_vv dV (24)

Last step involves rewriting the leftmost term by splitting the gradient of v (a 3x3 matrix) into a symmetric and anti-symmetric matrix:

∇v “ 1

2p∇v ` ∇v^|q ` 1

2p∇v ´ ∇v^|q (25)

When substituting the gradient in equation 24 the contraction operator between the symmetrical σ and the anti-symmetrical part of ∇v becomes zero. This is due to their symmetrical properties and that the anti-symmetrical matrix has a diagonal of zeros.

Furthermore, recalling the definition of the engineer strain, section 2.1 & equation 6, the symmetric matrix is the engineer strain for a deformation v.

ż

Ω

σ : ∇v dV “ ż

Ω

σ : 1

2p∇v ` ∇v^|q looooooomooooooon

“εpvq

` σ : 1

2p∇v ´ ∇v^|q looooooooomooooooooon

“0

dV “ ż

Ω

σ : εpvq dV (26)

The final weak formulation of the equilibrium equation is gained by substituting equation 26 in equation 24. The resulting PDE is presented in equation 27.

(21)

2.2.2 General solution for linear elastic case

Solving equation 27 with FEM requires that the Ω domain gets discretized into elements and nodes. The specific geometry of the elements can be varied, but for this derivation the elements will be tetrahedrons with four nodes. Each node can be displaced in three dimensions, figure 10 shows a tetrahedron element K and its four nodes.

n₁ n₂

n₃

n₄

K

w₂ v2

u₂ u₁ v₁

w₁

u₃

v₃ w₃

u₄

v₄ w₄

Figure 10: A tetrahedron element with three dimensions of nodal displacement.

Regarding the shape functions, different implementations may choose different sets of shape functions. Keeping the derivation general, the functions class is denoted Q. The unspecified set of shape functions will be denoted tϕiu⁴₁. Defining a suitable function space where the approximate solution uh is to be found:

u_h P Vh “ tv P rQpKqs³ : v|BK “ BCu Vh “ spanptϕiuq

Some changes in notation for the stress and engineer strain allows representing Hooke’s law, section 2.1.2, with a matrix relationship. The specifics of the matrix C was presented in equation 15.

ε “

”

ε₁ ε₂ ε₃ 2ε₁₂ 2ε₁₃ 2ε₂₃ ı|

σ “

”

σ1 σ2 σ3 σ12 σ13 σ23

ı|

σ “ Cε (28)

ùñ ε : σ “ ε^|σ “ ε^|Cε (29)

The element displacement is the sum of the nodal displacements, superposition of displacement allowed due to linear deformation [15]. Each node can be displaced in three

(22)

dimensions and the displacement in each direction is approximated by the shape functions. The ansatz for the element displacement thus becomes:

u_h “

»

— – u v w

fi ffi fl “

”_ϕ₁ ₀ _{0 ϕ}₂ ₀ _{0 ϕ}₃ ₀ _{0 ϕ}₄ ₀ ₀

0 ϕ1 0 0 ϕ2 0 0 ϕ3 0 0 ϕ4 0 0 0 ϕ1 0 0 ϕ2 0 0 ϕ3 0 0 ϕ4

ı loooooooooooooooooooomoooooooooooooooooooon

ϕ

»

—

— –

u1

v1

w1

u2

v2

w2

u3

v3

w3

u4

v4

w4

fi ffi ffi ffi ffi ffi fl loomoon

d

(30)

The strain-displacement relation from equation 6 can also be rewritten as a matrix following the introduced notation. The ansatz for element displacement thus relates to the engineers strain by the differential matrix L. Together with the shape function matrix ϕ the strain matrix, B, is formed.

ε “

»

—

— –

B

Bx1 0 0

0 ^B

Bx2 0

0 0 ^B

Bx3

0 ^B

Bx3

B Bx2

B

Bx3 0 ^B

Bx1

B Bx3

B

Bx1 0

fi ffi ffi ffi ffi ffi ffi ffi ffi fl looooooooomooooooooon

L

u_h “ Lϕ loomoon

B

d “ Bd (31)

The finite element formulation is formed by combining the weak form PDE (equation 27), the stress-strain relationship (equation 28) and the strain-displacement approximation (equation 31). ke is the stiffness matrix for the element (not to be confused with stiffness matrix C from section 2.1.2) and fe is the load vector. The stiffness matrix for the whole domain is simply the sum of all elements’ stiffness matrices. Note that σ ¨ n has been replaced with the traction vector t, the boundary force vector.

` ż

K

B^|CB dV loooooooomoooooooon

ke

˘ d “ ż

BK

ϕ^|t dS ` ż

K

ϕ^|fv dV looooooooooooooomooooooooooooooon

fe

ðñ ked “ fe

(32)

K_Ω “

Ne

ÿ

e“1

k_e (33)

2.2.3 Homogenization and Representative Volume Element

(23)

a mesh able to distinguish between material stiffness on the microscopic level would be incomprehensibly small and computationally unfeasible. A more realistic approach is by creating a similar homogeneous material which display the same material characteristics on the macroscopic level as the heterogeneous material. Finding these characteristics to model the heterogeneous material as homogeneous is done through a process called homogenization.

Homogenization

Figure 11: Basic principle of homogenization shown with a composite material.

The homogenization process for a material in which the heterogeneity occur in the microscopic level might still require an enormous amount of computations. Luckily each material heterogeneity is not unique, rather, the microscopic heterogeneities are often highly repetitive. One can often find a unit cell from which the entire material is built.

A better way of homogenizing is by utilizing the repetitiveness of the unit cell and only homogenizing the small unit cell. This unit cell is called representative volume element or RVE. The criteria for the RVE is that it must be representative of the whole material when repeated and preferably as small as possible. A basic illustration of choosing a suitable RVE for a heterogeneous material can be seen in figure 12.

Valid but needlessly large RVE.

Would harm computational ease

Invalid RVE, does not correctly represent the entire

structure

Suitable RVE, smallest volume still

able to represent the whole

Figure 12: Visualisation of choosing an appropriate RVE.

With a suitable RVE the homogenization process starts. The goal is finding an average of the microscopic stiffness so that the macroscopic stiffness matrix C is a constant in the macroscopic environment, thereby making it usable in equation 28. Notation wise, index with lowercase m is the microscopic behaviour and uppercase M is the macroscopic behaviour.

Similarly to the finite element formulation for the macroscopic level, section 2.2.1, the

(24)

however, additional approximations and boundary conditions appear. Firstly, the body forces are so small they can be neglected reducing the equilibrium equation. Secondly, since the unit cell is meant to behave identically throughout the material and are built atop each other the boundaries need to behave symmetrically. A periodic boundary condition is therefore applied to the equilibrium PDE, requiring that the displacement be symmetrical and the traction vector anti-symmetrical.

$

’&

’%

∇ ¨ σ “ 0 on Ω

u symmetric on BΩ

t “ σ ¨ n anti symmetric on BΩ

(34)

The basic principle of retrieving CM is by either finding the volume average of microscopic stress or strain εm & σm. The two method are therefore logically called stress/strain averaging. In order for the microscopic average to equal the macroscopic value the, so called, Hill-Mandel Condition needs to be satisfied. The proof and validation that the periodic boundary condition satisfy the Hill-Mandel Condition as well as the condition itself will be left out of this report, but can be further understood by reading Yvonnet’s book on homogenization [16].

In this thesis the average strain method will be used as it is better fit for the displacement- based finite element framework [17]. Acquiring CM is done by imposing a macroscopic strain on the boundary and solving the equation 34 in a similar way as described in section 2.2.1. When the solution is computed the volume average of stress, equation 35, is calculated. This average is by the Hill-Mandel lemma [16] the same as the macroscopic stress σM. By doing multiple calculation and iteratively loading components of the imposed strain εM the stiffness matrix CM can be constructed and the heterogeneous body homogenized.

σ_M “ 1 V_Ω

ż

Ωm

σ_m dΩ (35)

(25)

2.3 Wooden materials

Wooden materials come in different shapes and forms. The two major wood types are softwoods and hardwoods. Softwood is wood or lumber from a coniferous tree such as pine, spruce and larch trees. Softwoods contain no vessels and are characterized by its non-porosity. Hardwood on the other hand comes from dicot trees. Oak and maple trees are examples of dicot trees [18]. This study will focus on softwoods and thus coniferous trees alone.

2.3.1 Wood properties and cell structure

Figure 13: A visualization of wood on different scales [17].

On a macroscopic scale wooden materials are distinguished by their repeating sequences of dark and light bands called growth rings as visualized in the leftmost picture of figure 13. Growth rings are a consequence of the trees growth process. The light bands are called earlywood (EW) and the dark bands are called latewood (LW). The names are explanatory as earlywood is produced out of cells growing early in the season and latewood is produced later in the season [18]. The percentage of LW and EW can differ within a certain species as the factor is linked to seasonal growth and thus depend on geographical location. The EW and LW transition of a timber is given in the second picture from the left in figure 13.

If growth ring curvature is neglected, wood can be assumed to be ortothropic. The wood will thus have three orthogonal planes as shown in figure 14. The axes are called longitudinal/axial, tangential and radial. Wood is viscoelastic in nature. At small strains however, wood behaves linear elastically in all three directions, thus the elastic modulus plays a significant role for the characteristics of the structure. Stiffness and strength is greatest in the axial/longitudinal direction, i.e. the elastic modulus is largest in these directions. In the tangential and radial direction, the elastic moduli are of the same order of magnitude, however some differences can be observed which will be discussed later [4].

(26)

Figure 14: The three orthogonal axes of a tree; longitudinal/axial, tangential and radial [4]

Further down in scale wood is a cellular material made out of anisotropic wood cells called fibres or tracheids. The tracheids are elongated, pointed and constituates the bulk of the material [19]. The tracheids lie in a honeycomb-like hexagonal array as shown in the third picture from the left in figure 13. The figure shows the tracheid matrix which elongates in the longitudinal/axial direction. As a reference, figure 15 shows a microscopic image of a typical honeycomb cell structure in a softwood. The more spacious parts are EW and the denser parts are LW.

Apart from the tracheids, there are also horizontal fibre cells called rays which trans- port sap across the grain. They are pointed along the radial axis [18]. Rays compose approximately 5-10% of softwoods [20].

(27)

in nature [19] and are organized in fibrils made out of crystalline cellulose [18]. The fibrous cellulose can explain some of the anisotropic characteristics of wood. This means that the mechanical properties differ if a force is applied along or across fibres. Cell shape also affects the anisotropy. Outside and in between the cell walls is a matrix of hemicellulose and lignin which binds the cell together [4]. The cell wall can thus be approximated as a fiber reinforced composite [19].

The tracheid cell wall can be looked at in even higher detail as shown in the first picture from the left in figure 16. The tracheid structure is composed of a compound middle lamella (CML) consisting of the primary wall layer (P) and the Middle Lamella (ML) as well as secondary cell wall layers called S1, S2 and S3 respectively. The main differences between these layers is the cell wall thickness and the distribution of the microfibrils in the cell wall. ML is primarily made out of ligning and its main purpose is to bind the tracheids together [22]. In the P-layer which is the outermost layer of the cell wall, the microfibrils are packed loosely and at random. The S1 and S3 layers in the secondary layer of the wall are very thin. Research has shown that the S2 layer makes up about 80%

of the cell wall thickness [23]. The cell wall properties of the S2 layer should consequently have the most influence over the physical properties of the wood cells. The dark lines in the S2 layers represent the cellulose microfibrils and their direction. A closer look on the microfibrils are given in the second picture from the left in figure 16.

Figure 16: Cell wall layers and the microfibrils in the S2 layer [17].

2.3.2 Modelling wood structure

The established method when modelling wood structures is to extract the macroscopic properties from the microscopic cell structure. The macroscopic mechanical properties can thus be produced by cellular properties such as cell shape and relative density. Cell shape properties can for example be those shown in figure 17 which explain the shape

(28)

thickness, T is the length in the tangential direction and R is the length in the radial direction. Another important property is the so called microfibrilar angle (MFA). The MFA describes the angle of the microfibrils in the cell wall compared to the tracheids elongation in the longitudinal/axial direction. The second picture from the left in figure 16 shows the MFA as an angle θ in the S2 level. These cellular properties has shown to have direct effect on the macroscopic mechanical properties. For example, the MFA has a direct negative influence on the elastic modulus in the axial direction. This is also true for the moisture content of wood. An increase in moisture content leads to a reduced stiffness [24]. On the other hand, according to previous research, external factors such as age, air moisture or air temperature play less of a role. Another important factor is that the macroscopic properties of wood such as density differ significantly between different tree species. Cell wall properties are on the other hand more similar, which can be utilized when creating a model with the purpose of simulating different wood types [19]. As can be visualized in figure 15 most of the cells in the hexagonal array has shapes closely resembling rectangles. Modelling could thus be greatly simplified with this geometric assumption. The cell shape properties in figure 17 can easily be translated to the properties of a rectangle. This approximation has been implemented in previous research. However one observed effect is that a rectangular model leads to increased elastic moduli compared to a hexagonal model [25].

Figure 17: Geometric parameters of a hexagonal tracheid [17].

On a cellular level, EW and LW has different properties. LW is characterized by thick

(29)

elongated in the radial direction, will increase the radial stiffness. Rays are higher in density than tracheids, and as density increases the stiffness, the radial elastic modulus is larger than the tangential. According to Gibson’s analytical model, elastic modulus in the radial and tangential directions approximates to ER« 1.5ET. These will however still be significantly smaller than elastic modulus in the axial direction. This means the cell is stiffer against axial deformation compared to bending. One reason for this is that the characteristics of the honeycomb structure only are distinguished in the radial-tangential plane. Hexagonal, prismatic cells are intrinsically stiffer along the prism axis. The wood structure is thus by nature stiffer in the out of plane direction. Another reason for the axial stiffness is that the microfibrils of cellulose in the cell walls are placed closer together in the axial direction [19].

(30)

3 Method

3.1 RVE modelling of wood cell structure

For the purpose of conducting parameter studies on an arbitrary softwood an RVE model was built in Comsol Multiphysics. Comsol is a commercial software used for building geometries and simulating physical phenomenons formulated by differential equations.

In Comsol, the equations are mostly solved with the finite elements method (FEM).

Broadly, the final model consists of two geometries representing individual cells of EW and LW respectively as well as a third geometry representing the entire wood structure.

Examples of the three geometries are given in figures 18a, 18b and 18c.

When constructing the model, the EW and LW RVEs were built first. The individual tracheid cells in the honeycomb structure were approximated to be rectangular rather than hexagonal. This approximation was motivated by studying images such as figure 19a and 19b where many of the individual cells resemble the geometry of a rectangle. This approximation greatly simplifies the model which is necessary for both computational and administrative reasons when conducting the parameter studies. To account for the entire honeycomb geometry, the RVE’s were built as in figures 18a and 18b where the extensions of the rectangle are pointed in the tangential direction. These two RVE geometries can be compared to the image of a real cell demonstrated in figure 19c. As can be visualized from figures 19a-19c, an RVE built from a single cell can if duplicated represent the entire cell structure. However, an equivalent and more computationally convenient RVE was discovered after the linear elastic case study which will be discusses in section 3.3.

(a) Model of earlywood cell (b) Model of latewood cell

(c) Model of EW and LW distribution

Figure 18: Images of the simulated cellular model

(31)

(a) Honeycomb cell structure (b) Several tracheid cells

(c) A single tracheid cell

Figure 19: A microscopic image of a Cell structure for Norwegian Spruce on different scales [26]

The tracheid dimensions as well as MFA differ between tree species and the RVE geometries were thus adjusted based on which species was studies. The considered parameters were ’cell wall thickness’, ’radial diameter’, ’tangential diameter’ and MFA. In table 1 these parameters are specified for EW and LW respectively, for the species in which the model is validated for. The tracheid dimensions for Norwegian Spruce is well documented [27]. The values for pulpwood is used for Norwegian Spruce. Pulpwood is found far out on the tree log, where the curvature of the growth rings is not as prevalent. As discussed in 2.3, the ortothropic characteristic of wood is more evident if curvature can be neglected.

The dimensions of a European Larch tracheid [21] and a Scots Pine tracheid [28] are also given from previous research.

An MFA of 10° was used for Norwegian spruce based on previous research [29]. For European larch, the MFA has in average been shown for Normal wood to be 18° for EW and 5° for LW [30]. The model assumes 25% LW and 75% EW which as a weighted average results in an MFA of 14.75° for European Larch. For Scots Pine the MFA is chosen as 13° based on previous research for an annual ring far from the pith [31].

(32)

Table 1: Tracheid dimensions in µm and MFA in degrees, specified for EW and LW for different species of tree.

Norwegian spruce European Larch Scots Pine

EW LW EW LW EW LW

Cell wall thickness 2.13 3.88 3.60 7.68 2.0 3.5 Radial length 31.75 23.90 52.75 22.65 32.6 23.7 Tangential length 30.74 30.05 26.38 30.7 29.2

MFA 10.0 14.75 13.0

For the actual implementation of the MFA as seen in figure 13, local coordinate systems were added to each side consisting of the same color in figure 20. The rotation matrices in equation 36 express how the local coordinate systems relate the direction of the mi- crofibrills to the global coordinate system. The coloured indices in equation 36 reference figure 20.

Figure 20: Colour representation of RVE sides with same local coordinate systems

R_Blue “

„ cosp´MFAq 0 sinp´MFAq

0 1 0

´ sinp´MFAq 0 cosp´MFAq



, R_Red “

„ cospMFAq 0 sinpMFAq

0 1 0

´ sinpMFAq 0 cospMFAq



R_Green “

„cosp´MFAq ´ sinp´MFAq 0 sinp´MFAq cosp´MFAq 0

0 0 1



, R_{Y ellow} “

„cospMFAq ´ sinpMFAq 0 sinpMFAq cospMFAq 0

0 0 1

 (36)

Micromechanical modelling of creep in wooden materials

Examensarbete 15 hp June 2021

Micromechanical modelling of creep in wooden materials

Kevin Coleman

Oskar Falkeström

Malin Nilsson

Abstract

Micromechanical modelling of creep in wooden materials

Populärvetenskaplig sammanfattning

Acknowledgements

Contents

1 Introduction

1.1 Objective

1.2 Goals

2 Theory

2.1 Mechanics of Materials

2.2 Numerical methods for mechanics of materials

K

2.3 Wooden materials

3 Method

3.1 RVE modelling of wood cell structure