# New modelling approaches to predict wood properties from its cellular structure: image-based representation and meshless methods

Patrick Perre, Giana Almeida, Mehdi Ayouz, Xavier Frank

► To cite this version:

Patrick Perre, Giana Almeida, Mehdi Ayouz, Xavier Frank. New modelling approaches to predict wood properties from its cellular structure: image-based representation and meshless methods. Annals of Forest Science, 2016, 73 (1), pp.147-162. 10.1007/s13595-015-0519-0. hal-01237586

HAL Id: hal-01237586

https://centralesupelec.hal.science/hal-01237586

Submitted on 27 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# New modelling approaches to predict wood properties from its cellular structure: image-based representation and meshless methods

Patrick Perré<fnref n="1" /> · Giana Almeida<fnref n="2" /> · Mehdi Ayouz<fnref n="1" /> · Xavier Frank<fnref n="3" />

Received: 2 December 2014 / Accepted: 24 August 2015

© INRA and Springer-Verlag France 2015. This article is published with open access at Springerlink.com

## Abstract

*Key message* The real tissue structure, including local anisotropy directions, is defined from anatomical images of wood. Using this digital representation, thermal/mass diffusivity and mechanical properties (stiffness, large deformation, rupture) are successfully predicted for any anatomical pattern using suitable meshless methods.

*Introduction* Wood, an engineering material of biological origin, presents a huge variability among and within species. Understanding structure/property relationships in wood would allow engineers to control and benefit from this variability. Several decades of studies in this domain have emphasised

the need to account simultaneously for the phase properties and the phase morphology in order to be able to predict wood properties from its anatomical features. This work is focused on the possibilities offered by meshless computational methods to perform upscaling in wood using actual tissue morphologies obtained by microscopic images.

*Methods* After a section devoted to the representation step, the digital representation of wood anatomy by image processing and grid generation, the papers focuses on three meshless methods applied to predict different macroscopic properties in the transverse plane of wood (spruce earlywood, spruce late-wood and poplar): Lattice Boltzmann Method (LBM) allows thermal conductivity and mass diffusivity to be predicted, Material Point Method (MPM) deals with rigidity and compression at large deformations and peridynamic method is used to predict the fracture pathway in the cellular arrangement.

*Results* This work proves that the macroscopic properties can be predicted with quite good accuracy using only the cellular structure and published data regarding the cell wall properties. A whole set of results is presented and commented, including the anisotropic ratios between radial and tangential directions.

**Keywords** Lattice Boltzmann Method · Material Point Method · Peridynamic · Real anatomical structure · Mechanical properties · Stiffness · Fracture · Diffusion · Thermal conduction · Collapse

## 1 A short review of structure/properties relationships in wood

As this paper is published in the 50th anniversary issue of the *Annals of Forest Science*, we reviewed literature starting from

Handling Editor: Erwin Dreyer

**Contribution of the co-authors** Patrick Perré: coordinating the research project, contributing in writing the paper, MPM modelling

Giana Almeida: contributing in writing the paper, wood anatomical structure

Mehdi Ayouz: contributing in writing the paper, Lattice Boltzmann modelling

Xavier Frank: contributing in writing the paper, peridynamic modelling

<figure><img src="image_2.png" /></figure>

Patrick Perré
patrick.perre@ecp.fr

Giana Almeida
giana.perre@agroparistech.fr

Mehdi Ayouz
mehdi.ayouz@ecp.fr

Xavier Frank
xavier.frank@supagro.inra.fr

<footnote n="1">1 LGPM, Ecole Centrale Paris, Grande Voie des Vignes, 92 290 Châtenay-Malabry, France</footnote>

<footnote n="2">2 UMR 1145, Agro ParisTech, INRA, 1 Av. des Olympiades, F-91744 Massy, France</footnote>

<footnote n="3">3 UMR 1208 IATE, INRA-CIRAD-UM-Supagro, 2 place Pierre Viala, 34060 Montpellier, France</footnote>
the 1960s. This choice was motivated by the fact that most of
the research dedicated to this topic began to appear during this
period. Over the years, the investigation of the structures/
properties in wood was motivated by two major facts:

- The huge variability of wood properties, within species or between species

- The intuition, provided by studying wood anatomy, that the anatomical pattern is able to explain, at least partly, this variability

When wood is used as a structural material, two properties are of primary matter: longitudinal stiffness and transversal shrinkage. This is why most research efforts have focused on these properties. In this case, the term “property” refers to the so-called macroscopic property, as defined using solid wood samples, with typical sizes of several centimeters in each direction, possibly some tens of centimetres in the longitudinal direction. Studies showed that density, which represents the quantity of lignocellulosic matter embedded in the wood, is highly variable (ranging from 100 to 1200 kg m⁻³ among species) and is likely to account for most of the variability in the properties of wood. Therefore, it is not surprising that the first attempts to predict wood properties were in the form of linear or non-linear correlations dependent on density. For example, this strategy works nicely for longitudinal stiffness and hardness (Kollmann and Côté 1968; Bosshard 1984). The hidden, and coarse, assumption made in this simple approach is that all phases of wood are in parallel and aligned along the longitudinal direction. Using this simple upscaling strategy, the macroscopic property is simply a weighted average of the microscopic property over all phases of the heterogeneous medium. This explains why poor correlations are obtained for certain properties, such as transverse shrinkage. In addition, even if a rather good correlation is obtained, the residual variability is too large to accurately predict wood properties. This means that, for a given sample, the deviation from the general correlation might be large in terms of relative error.

To improve the knowledge of wood properties and, more
specifically, to elaborate models to explain the dramatic vari-
ability observed in its properties, microscopic features (at the
anatomical and ultra-structural levels) have to be considered.
The task is not easy, as many spatial scales contribute to mac-
roscopic behaviours. It is now well established that three spa-
tial scales are particularly relevant: the cell wall level, namely
through the microfibril angle (MFA); the cellular structure,
which explains the tissue properties; and the anatomical pat-
tern, in which all the anatomical tissues are organised in pro-
portion and in space. A multiscale approach would be ideal to
account for all these wood features. In practice, key factors
were gradually introduced to explain the deviations observed
from the statistical correlations.

For example, regarding the stiffness or shrinkage values in
the longitudinal direction, the MFA in the secondary cell wall
was proposed as an explanatory parameter decades ago
(Harris and Meylan 1965; Meylan and Probine 1969). These
findings were a major improvement in the understanding of
the longitudinal behaviour of wood. The determination of the
MFA by X-ray diffraction (Cave 1966) played an important
part in this progress. It is noteworthy to mention that such a
clear influence of the MFA on the longitudinal properties of
wood is a rare example where a factor at a low spatial scale
(nanoscale) has a straightforward effect on a macroscopic
property: scaling in material sciences is usually more complex
and involves successive upscaling steps. This relative simplic-
ity, due to the fact that all solid components act in series in the
longitudinal directions, allowed analytical models to be pro-
posed in the same period (Barber 1968).

The understanding and prediction of structure-property re-
lationships in the transverse plane of wood (radial-tangential)
are more complicated. In this plane, solid components act both
in parallel and in series at different spatial scales (multi-lay-
ered cell wall, cellular morphology, anatomical pattern).

Again, observations and measurements came before
modelling. Several scientists tried to use anatomical features
as input parameters in statistical explanations. For example,
the occurrence of ray cells (Barkas 1941; Boutelje 1962;
Kelsey 1963; Keller and Thiercelin 1975; Guitard and El
Amri 1987) and the shape of cells (Mariaux and Narboni
1978; Masseran and Mariaux 1985) were tested as possible
explanations of shrinkage variability. In 1989, Mariaux ob-
served that the transverse anisotropy of tissues depends on
the mean elongation of the cell, but that shrinkage was not
isotropic for “isotropic” cells (same mean diameter in both the
radial and tangential directions).

In the meantime, theoretical works were proposed to ex-
plain transverse properties from the cellular structure (Barber
and Meylan 1964; Gillis 1972; Koponen et al. 1991; Gibson
and Ashby 1988). These works are based on analytical models
and assume that the cellular structure is represented by a
unique tracheid.

Indeed, earlier works pointed out the need to account for
the real morphology of the cellular structure for a successful
prediction of transverse properties (Farruggia 1998; Holmberg
et al. 1999; Perré 2001; Nairn 2006; Abbasi 2013). To do this,
computational approaches, which take advantage of advances
made in applied mathematics and mechanics regarding scaling
approaches (Sanchez-Palencia 1980; Suquet 1985), need to be
applied to wood science.

Upscaling, such as the homogenisation of periodic struc-
tures, is a deterministic approach that includes several steps
(Fig. 1):

1. Representation: choice of the representative elementary volume (REV), also called the Unit Cell. This REV

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 1** Upscaling in material sciences: principle</figcaption></figure>

should be defined in a suitable way for subsequent calculation (finite element mesh, collection of material points...).

2. Characterisation of the properties for each phase of the unit cell.

3. Solution: the theoretical formulation (i.e. the homogenisation of periodic media) has to be solved using a suitable computational method.

4. Validation: the predicted macroscopic properties should be tested against experimental data.

5. Localisation: this step is not mandatory, but it allows the local (microscopic) fields (shrinkage, strain, stress, temperature, etc.) to be computed inside the REV under the macroscopic conditions applied to the product.

Figure 1 indicates that a macroscopic property depends on both the local properties of the different phases of the REV (2) and their spatial organisation (morphology) (1). For example, in the transverse plane, the anisotropy of tissue stiffness is mostly explained by the cellular structure (morphology) (Farruggia 1998; Perré 2001), while shrinkage, namely the difference between normal wood and reaction wood, strongly involves the cell wall behaviour in the transverse plane (local property) (Watanabe and Norimoto 1996; Perré and Huber 2007) and, eventually, by the alternation of earlywood and latewood (Lanvermann 2014; Perré and Turner 2002, 2008). In the case of a biological product such as wood or lignocellulosic materials, steps 1 and 2 are particularly difficult. In addition, they are strongly entwined, yet the quality of the prediction depends mostly on the qualities of these two first steps. Whatever the target scale of the scaling approach (the macroscopic scale), the different phases of the unit cell have to be well defined, both in shape and in values, at a smaller scale. For example, if the goal is to obtain the properties of the cellular structure, the representation step consists of defining the size of the REV and the cell morphology inside this volume. Then, computing the macroscopic properties requires the cell wall behaviour to be used as input data (step 2). Regarding this step, products of biological origin are different from other materials in the sense that the constituents of the unit cell do

not exist alone and are, therefore, very difficult to characterise. Indeed, two strategies co-exist. The first one consists of direct characterisation. In this case, the size of the sample has to be sufficiently reduced so that it is representative of the local scale. A number of researchers have proposed this approach for mechanical properties (Bergander and Salmén 2000; Sedighi-Gilani and Navi 2007; Farruggia and Perré 2000; Perré et al. 2013) and for shrinkage (Perré and Huber 2007; Perré 2007; Almeida et al. 2009; Almeida et al. 2014). In the second strategy, the local properties are deduced from previous scaling approaches (Holmberg et al. 1999; Hofstetter et al. 2005; Neagu and Gamstedt 2007). If the phase morphology is correctly represented, a third possibility exists: a complete scaling approach is performed, but, instead of using this approach to predict the macroscopic properties, an inverse analysis allows local properties to be deduced from macroscopic measurements (arrow 6 of Fig. 1). This approach was applied successfully to demonstrate an equivalent stiffness value of the cell wall in the transverse plane (Farruggia 1998).

As this paper is devoted to the behaviour of wood tissues, the local scale is, therefore, the cell wall, whereas the macroscopic scale will be that of the tissue (a representative subset of cells). It is well known that the cell wall is multi-layered and presents a specific microfibril orientation and macromolecular organisation, which explains its heterogeneous and anisotropic properties (Salmen 2004; Gierlinger and Schwanninger 2006).

Several options exist regarding the modelling approach. The simplest one would assume the cell wall to be homogeneous and isotropic. At this level, one has to be aware that the assumptions made regarding the solid phase of the unit cell have dramatic effects on steps 1 and 2. The representation step should provide a geometrical description relevant to the assumptions. If needed, heterogeneities and material directions must be generated, ideally using real anatomical structures, which makes image processing difficult. In this sense, modern imaging tools, such as environmental electron scanning microscopy (ESEM), confocal scanning laser microscopy (CLSM), confocal Raman microscopy and computed X-ray µ-tomography, are of great interest (Perré 2011). It is now quite common to acquire chemical images of wood sections (Gierlinger and Schwanninger 2006; Perré 2011; Cabrolier 2012).

Once the image exists, the representation step consists of building a digital representation of the image. This is commonly performed by generating a finite element (FE) mesh. In this case, the computational part will just solve the basic problems of the FE mesh using FE theory. The simplest way to build a FE mesh consists of generating a Q4 element from each pixel of the image (in 2D) or a Q8 element from each voxel of the image (3D). More sophisticated approaches define the geometry from the image (contours in 2D and surface in 3D) and mesh the volume from the geometry. This strategy
is quite well developed in materials sciences (Ulrich et al. 1998; Charras and Guldberg 2000; Kwon et al. 2003). However, currently, a whole family of meshless methods are available that facilitate the representation step (Belytschko et al. 1996; Frank and Perré 2010). This paper is devoted to this new approach that is able to take advantage of modern imaging tools in defining the real pore morphology. In the following section, the mechanical and transfer properties (mass and thermal diffusion) of wood tissue were computed using image-based representation and solved with various meshless methods. The choice of the method depends on the behaviour of interest:

- The Lattice Boltzmann Method (LBM) was chosen for thermal and mass diffusion. This is an elegant method in which both the value and the flux are known at each grid point (Succi 2001; Mohamad 2011).

- The Material Point Method (MPM) was chosen for the mechanical behaviour. This method predicts stiffness, but also analyses the behaviour of the cellular structure in large deformations (Sulsky et al. 1994; Sulsky and Schreyer 2004; Bardenhagen and Kober 2004).

- Finally, fracture in the cellular structure was simulated using the peridynamic method (Silling 2000; Silling et al. 2007; Silling and Askari 2005).

The next section of this paper is devoted to the image-based representation of cellular structures. Three tissues will serve as input files in Section 4: spruce (Picea abies) earlywood, spruce latewood and poplar (Populus tremula x alba). Poplar is proposed here as one example of a dual-porosity organisation (vessel lumens and fibre lumens), which requires ca. 200, 000 points to accurately represent the actual morphology (the same number of points allows the sub-layers of the cell wall to be represented in single-porosity tissues). Each method is then briefly presented, and the last section presents results and discussion.

## 2 Image-based representation

Whatever the meshless method used, the digital representation consists of a collection of points. Throughout this work, this point collection is generated directly from anatomical images of the cellular tissue. The vectorial image processing facilities provided by the custom software MeshPore allows this image-based grid generation to be efficiently performed (Perré 2005). The image is first segmented to generate contours of the lumens: closed chains of vectors for internal lumens and open chains for border lumens. A regular grid of points is then generated at any desired point density: the grid refinement is, therefore, completely independent of the pixel size of

the anatomical image. The contours are then used to determine whether a point belongs to the solid phase or to the voids (Fig. 2). The type “void” is allocated to each point located inside a closed chain (lumen) or inside a border chain, while the type “cell wall” is allocated to all other points. In the case of MPM or peridynamic modeling, the voids are discarded because their mechanical contribution is negligible. This is not the case when the pores are filled with water, whose negative pressure is likely to provoke a collapse of the cellular structure (Perré et al. 2012).

Finally, for each material point of the solid phase, the local anisotropy directions are generated from the two nearest cell contours in such a way as to mimic how the living cell generated the secondary cell wall (Fig. 3). In the final grid, each point carries its own information, such as material type, local anisotropic directions, density, etc.

An important feature was added to the MeshPore software for the specific needs of the present paper. In the case of fragmentation, it is important to distinguish the middle lamella from the secondary cell wall layer. In this case, the position of the middle lamella was iteratively generated from the lumen contour using iterative increments along a direction perpendicular to the contour while accounting for a minimum gap between two neighbouring internal contours. To correctly represent the middle lamella at the intersection of three cells, the contour evolution locally stops when a maximum curvature is attained. Once all contours are stabilised, two “solid” types (secondary wall and middle lamella) are allocated to all points located inside the cell wall (Fig. 4). Figure 5 depicts the three tissue morphologies tested in the present work. The most refined meshes were generated for spruce (earlywood and latewood), in which the middle lamella was distinguished from the secondary cell wall. This feature was required to model fracture. The third case is a cellular arrangement of poplar, a ring-porous hardwood species. This example is interesting due to the double porosity: fibre and longitudinal parenchyma lumen and vessel lumen. Using a medium refinement at the cell wall level, this dual-scale porosity can be represented using a number of material points comparable to those needed in the case of spruce, where the refinement is used to separate the middle lamella and the secondary cell wall.

<figure><img src="image_6.png" /><figcaption>**Fig. 2** The vectorial image processing used to build a grid of material points from an anatomical image of wood. From left to right: ESEM image, vectorial contours and an example of a coarse grid</figcaption></figure>

<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 3 Using the neighbouring contours, an algorithm was developed to propagate the material directions to each material point of the solid phase (left)</figcaption></figure>

# 3 The meshless methods used in this work

## 3.1 Lattice Boltzmann Method

The LBM was first developed to solve the macroscopic momentum equation for viscous flows (Succi 2001). Based on velocity distributions on a regular lattice, this method

simulates the macroscopic behaviour as an emerging property of the discrete movement of particles (propagation and collision). Indeed, an asymptotic development of the lattice rules allows the macroscopic set of equations to be theoretically derived from the discrete rules. Following this strategy, LBM progressively became a general numerical method to solve any kind of partial differential equations. It belongs to the family of so-called meshless methods (Belytschko et al. 1996; Frank and Perré 2010) and, therefore, has interesting properties, such as simplicity of development, suitability for parallel computing and flexibility in geometrical shape (Frank et al. 2010).

In the standard LBM, the single-particle distribution function obeys the following equation:

$$f_i(\vec{x} + \vec{c}_i \delta t, t + \delta t) - f_i(\vec{x}, t) = \Omega_i(\vec{x}, t) \quad (1)$$

where $\vec{c}_i$ is the discrete particle velocity vector and the right-hand side term represents the collision term, the so-called single relaxation time. The set of vectors $\vec{c}_i$ are the possible velocities a particle can move from a lattice site to a nearest-neighbouring site at each time step $t+\delta t$. For simplicity, we used a “single relaxation time” form for the collision term:

$$\Omega_i(\vec{x}, t) = -\frac{1}{\tau} [f_i(\vec{x}, t) - f_i^{eq}(\vec{x}, t)], \quad (2)$$

where $\tau$ is the dimensionless relaxation time and $f_i^{eq}(\vec{x}, t)$ is the equilibrium distribution function at site $\vec{x}$ and time $t$. The above equation is well known as the Bhatnagar-Gross-Krook (BGK) approximation for the collision operator $\Omega_i(\vec{x}, t)$. This means that after a collision, the distribution function tends to an equilibrium Maxwell-Boltzmann distribution.

In the present work, a standard D2Q9 model was adopted, i.e. a 2D computational domain $(x, y)$ involving nine discrete velocities, $c_i \forall i=\{0, 1, \dots, 8\}$. These velocities allow the particles to jump to neighbouring nodes in one lattice time step

<figure><img src="image_9.png" /><figcaption>Fig. 4 Generation of internal contours and allocation of two solid phase types (blue=lumens, red=secondary wall and yellow=middle lamella)</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 5** Three examples of wood tissues used in the following sections (the number of points is given for mechanical simulation and, hence, it does not include the lumens)</figcaption></figure>

(Fig. 6). They are defined as follows:

$$ \vec{c}'_i = \begin{cases} (0,0)c & i=0 \\ (\pm 1,0)c, (0,\pm 1)c & i=1,2,3,4 \\ (\pm 1,\pm 1)c & i=5,6,7,8 \end{cases} \quad (3) $$

where the lattice speed, c, is related to the lattice spacing $c = \frac{\delta x}{\delta t}$. $\delta x(=\delta y)$ is the lattice increment and $\delta t$ is the lattice time step.

For a diffusion problem, it is appropriate to consider the equilibrium distribution function as constant, where no macroscopic velocity is involved. This can be formulated as follows:

$$ f_i^{\text{eq}}(x, y, t) = \omega_i \Psi(x, y, t), \quad (4) $$

where $\omega_i$ are the weighting factors and $\Psi(x, y, t)$ is the normalised macroscopic variable (temperature field $\theta$ or mass $X$ (see Section 4), such as $(x,y,t) \in [0,1]$. In the D2Q9 model, the weighting factors are given as (Mohamad 2011):

$$ \omega_i = \begin{cases} 4/9 & i=0 \\ 1/9 & i=1,2,3,4 \\ 1/36 & i=5,6,7,8 \end{cases} \quad (5) $$

The standard 2D diffusion equation for the normalised mass or temperature field,

$$ \frac{\partial \Psi}{\partial t} = D \left[ \frac{\partial^2 \Psi}{\partial x^2} + \frac{\partial^2 \Psi}{\partial y^2} \right], \quad (6) $$

was solved in a 2D rectangular porous medium, $[x_1,x_n] \times [y_1,y_m]$, using the LB method. In the above equation, $D$ (assumed constant) is the macroscopic mass or thermal diffusion coefficient. Initially, a constant field is assumed: $\Psi(x,y,0)=0$. Along the x-direction, Dirichlet boundary conditions were applied: $\Psi(x_1,y,t)=1$ and $\Psi(x_n,y,t)=0$. Adiabatic conditions were set along the y-direction: $\varphi(x,y,t)=0$ at $y=y_1$ and $y=y_n$.

In the literature (Higuera 1990), it has been proved that the diffusion equation can be derived using the Chapman-Enskog

<figure><img src="image_11.png" /><figcaption>**Fig. 6** The main steps of the LBM approach in the case of the D2Q9 model (the thickness of the arrows represents the values of the distribution function $f_i$)</figcaption></figure>
expansion for the LBM equation (Eq. 1 including Eq. 2). The relationship between macro- and meso-scales gives

$$D = \frac{1}{3} \left( \tau - \frac{\delta t}{2} \right) c^2, \qquad (7)$$

where $\tau$ has a dimension of time, expressed in seconds. The macroscopic normalised quantities, $\Psi$ and related flux $\phi$, are obtained through moment summations in the velocity space:

$$\Psi(x, y, t) = \sum_{i=0}^{8} f_i(x, y, t) \qquad (8)$$

$$\vec{\phi}(x, y, t) = \sum_{i=0}^{8} \vec{c}_i f_i(x, y, t) \qquad (9)$$

The normalised flux has the following form:

$$\vec{\Phi}(x, y, t) = \frac{\left(\tau - \frac{1}{2}\right)}{\tau} \vec{\phi}(x, y, t) \qquad (10)$$

In LBM formalism, Eqs. (1) and (2) consist of two steps, collision and streaming. Then, boundary conditions must be applied to complete one time step (Fig. 6). On the domain boundaries, the components of the distribution functions for velocities going outwards are known from the streaming process. On the contrary, they are unknown for velocity components getting towards the domain. For example, the distribution functions ($f_1, f_5, f_8$) and ($f_3, f_6, f_7$) associated to velocities ($c_1, c_5, c_8$) and ($c_3, c_6, c_7$) are known at the lattice faces $x=x_1$ and $x=x_n$, respectively. On the contrary, ($f_3, f_6, f_7$) at $x=x_1$ and ($f_1, f_5, f_8$) at $x=x_n$ are unknown. For more details, see Bao et al. (2008) and references therein.

The implemented Dirichlet conditions can be expressed for the unknown components of the distribution function as follows (Bao et al. 2008):

$$\begin{aligned} f_i(x_1, y, t) &= \omega_i \Psi_1, & \forall i &= \{1, 5, 8\} \\ f_i(x_n, y, t) &= \omega_i \Psi_n, & \forall i &= \{3, 6, 7\} \end{aligned} \qquad (11)$$

where $\Psi_1$ and $\Psi_n$ are fictive variables obtained from balance equations to ensure the desired Dirichlet value. For example, for the plane ($x=x_1$), we use

$$\begin{aligned} \Psi_1 &= \frac{1}{\sum_i \omega_i} \left[ \Psi(x_1, y, t) - \left( \sum_{k=0}^{8} \tilde{f}_k(x_1, y, t) - \sum_i \tilde{f}_k(x_1, y, t) \right) \right] \qquad (12) \\ \forall i &= \{1, 5, 8\} \end{aligned}$$

where $\tilde{f}_k$ and $\tilde{f}_i$ represent the distribution functions after collision and streaming, and $k$ an index running over all velocities. Notice that $i$ runs only over the unknown velocities indexes.

For adiabatic boundaries, the unknown distribution functions ($f_2, f_5, f_6$) at $y=y_1$ and ($f_4, f_7, f_8$) at $y=y_m$ can be obtained

using as moving boundary conditions (see Lallemand and Luo 2003 for more details) such as

$$\begin{aligned} f_2(x, y_1, t) &= f_4(x, y_1, t), f_5(x, y_1, t) = f_8(x, y_1, t), \\ f_6(x, y_1, t) &= f_7(x, y_1, t), f_4(x, y_m, t) = f_2(x, y_m, t), \\ f_7(x, y_m, t) &= f_6(x, y_m, t), f_8(x, y_m, t) = f_5(x, y_m, t) \end{aligned} \qquad (13)$$

As our goal is to compute the equivalent macroscopic properties of the REV, the method has to be applied to heterogeneous domains. To do this, the LBM equation (Eq. 1), including the BGK approximation (Eq. 2), can be used by distinguishing two relaxation times, $\tau_1$ and $\tau_2$. The values of these two relaxation times are related to the two diffusion coefficients of the respective phase, as defined in Eq. (7). The discrete particle velocities are represented in a computational domain, such as earlywood of spruce from Fig. 2. The rectangular computational domain [$x_1, x_n$] [$y_1, y_m$] is, therefore, divided into two sub-sets: the cell lumens, $\Sigma_1$, and the solid phase, $\Sigma_2$, which are related to the macroscopic diffusion coefficients $D_1$ and $D_2$, respectively.

Numerically, within each time step, the main steps of the calculation can be summarised as follows (Fig. 6):

1. Calculate the equilibrium density function and perform the collision procedure at each node.

2. Stream the density distribution populations.

3. Apply boundary conditions for the density distribution function.

A Fortran95 LBM code was developed at LGPM (CentraleSupelec) to model the thermal and mass diffusion in heterogeneous medium. This code is flexible and can take account for different morphologies. To compute the equivalent property of the REV, this set of steps is iterated until equilibrium is reached. In general, a large number of iterations ($10^5-10^6$) is required to obtain the steady-state regime. The equivalent diffusivity is then computed as the average flux divided by the macroscopic gradient imposed by the Dirichlet boundary conditions.

## 3.2 Material Point Method

The Material Point Method (MPM) is derived from the Particle in Cell (PIC) method originally developed in the 1960s by Harlow in computational fluid mechanics. The MPM has been successfully applied to computational solid mechanics problems (Sulsky et al. 1994; Sulsky and Schreyer 2004; Bardenhagen and Kober 2004; Guilkey et al. 2006). One of the most important features of the MPM is that it allows a straightforward discretisation of complex material shapes, including direct discretisation from 2D or 3D images,
as well as an efficient and robust handling of contacts between material surfaces.

In the MPM, the domain of interest is defined by a collection of material points $p, p=1... N_p$. As a direct consequence of this principle, each material point carries its properties, such as position, velocity, acceleration, strain and stress, in a Lagrangian way (step 1 of Fig. 7). The particle mass, $m_p$, is computed from the medium density and the initial point distribution. In the classical version of MPM, no volume is affected to material points: therefore, the resulting density field is approximated by a non-continuous field defined by Dirac functions (capital letters are used for Lagrangian values):

$$ \rho(\vec{x}, t) = \sum_{p=1}^{N_p} m_p \delta(\vec{x} - \vec{X}_p(t)) \quad (14) $$

The $m_p$ values are constant over time, which ensures mass conservation. The weak form of the momentum equation is solved with the aid of a background grid. This grid consists of simple finite elements, usually isoparametric, four-node, quadrilateral elements in two dimensions. This simple grid only has to include all the material points and is reset at each time increment, thereby avoiding the classical FE problems that result from large deformations. The algorithm allows all quantities at material points to be updated stepwise. To do so, the first step is to map the material points on nodes $i, i=1... N_g$ of the FE mesh (step 2 of Fig. 7). Superscript $k$ refers to time step $k$ and superscript $L$ refers to the end of the Lagrangian step. Classical FE shape functions, $N_i$, are used for this purpose:

$$ \begin{aligned} S_{pi}^k &= N_i (\vec{X}_p^k) \\ G_{pi}^k &= \nabla|_{X_p^k} N_i \end{aligned} \quad (15) $$

Internal forces are mapped on the grid using boundary conditions and the stress level at each point $p$ (step 3 of Fig. 7):

<figure><img src="image_12.png" /><figcaption>Fig. 7 Principle of the Material Point Method (adapted from Sulsky et al. 1994)</figcaption></figure>

$$ \begin{aligned} \vec{f}_i^k &= \vec{f}_i^{\text{int},k} + \vec{f}_i^{\text{ext},k} \\ \text{where} \quad \vec{f}_i^{\text{int},k} &= -\sum_{p=1}^{N_p} m_p \left( \vec{G}_{ip}^{T,k} \sigma_p^k \right) \end{aligned} \quad (16) $$

In Eq. (16), $\overline{\sigma}^S$ is the specific stress tensor at point $p$. The external forces include the boundary conditions, which can be applied either to points or to nodes.

Similarly, the mass matrix is computed by mapping the point masses on the grid:

$$ m_{ij}^k = \sum_{p=1}^{N_p} m_p (S_{ip}^{T,k} S_{pj}^k) \quad (17) $$

At the grid nodes, the weak form of the momentum equation reduces to:

$$ \sum_{j=1}^{N_g} m_{ij}^k \vec{a}_j^k = \vec{f}_i^k \quad (18) $$

The velocities at the grid nodes are computed from the point velocities using a least-squares approach with point masses as weighting factors:

$$ \sum_{j=1}^{N_g} m_{ij}^k \vec{v}_j^k = \sum_{p=1}^{N_p} m_p \vec{v}_p^k N_i(X_p^k) \quad (19) $$

The momentum Eq. (18) allows the nodal velocities and positions to be updated at the end of the Lagrangian phase (step 4 of Fig. 7):

$$ \begin{aligned} \vec{v}_i^L &= \vec{v}_i^k + \Delta t \vec{a}_i^k \\ \vec{X}_i^L &= \vec{X}_i^k + \Delta t \vec{v}_i^L \end{aligned} \quad (20) $$

From the nodal information, the velocities and positions are updated back to material points using the element shape functions (step 5 of Fig. 7):

$$ \begin{aligned} \vec{V}_p^L &= \vec{V}_p^k + \Delta t \sum_{i=1}^{N_g} \vec{a}_i^k N_i(X_p^k) \\ \vec{X}_p^L &= \vec{X}_p^k + \Delta t \sum_{i=1}^{N_g} \vec{v}_i^k N_i(X_p^k) \end{aligned} \quad (21) $$

The strain increment of point $p$ is computed from the element shape functions:

$$ \Delta \epsilon_p = \frac{\Delta t}{2} \sum_{i=1}^{N_g} \left\{ \vec{G}_{ip}^k \vec{v}_i^L + \left( \vec{G}_{ip}^k \vec{v}_i^L \right)^T \right\} \quad (22) $$

The stress tensor at each material point is then simply updated by applying the constitutive equations. The updated stress
tensor, point position and velocity are now available to proceed further in time. The grid is usually reset at its original shape before proceeding to the next time step (step 6 of Fig. 7). MPM simulations have been proposed for wood, elastic and plastic constitutive equations (Nairn 2006), but with a quite poor morphological description (no local anisotropy nor a refined cell wall description) compared with the present work.

The MPM code developed at LGPM (CentraleSupelec), *MPM_Pore* accounts for large displacements and for anisotropic materials. *MPM_Pore* is written in Fortran 95 and parallelised by domain decomposition using Message Passing Interface (MPI) instructions.

### 3.3 Peridynamic approach

The growth of cracks in heterogeneous materials is of crucial interest in many fundamental and industrial domains, and has been extensively studied using various numerical approaches, such as the finite element method (Zavattieri et al. 2001; De la Osa et al. 2009; Sancho et al. 2007; Ruiz et al. 2000; Itakura et al. 2005), deformable lattice methods (Kitsunezaki 2013) and lattice element methods (Topin et al. 2007; Affes et al. 2012), to cite a few. Recently, the peridynamic approach emerged as an alternative method, based on integral equations, rather than partial differential equations (Silling 2000).

In the peridynamic framework, each material point is involved in interactions with points in a limited neighbourhood $H_r = \{\vec{r}' \mid \|\vec{r}' - \vec{r}'\| < \delta\}$, where $\vec{r}'$ is a point position at $t=0$, $\vec{r}'$ are neighbouring positions and $\delta$ is a cutoff distance called the horizon (Fig. 8).

Such a paradigm is very different from classical continuous matter mechanics, in which only short-range forces are involved. A bond is defined as $\vec{\xi} = \vec{r}' - \vec{r}$ when $\vec{r}'$ is in the neighbourhood of $\vec{r}$, $\vec{u}(\vec{r}, t)$ is the displacement and $\vec{\eta} = \vec{u}(\vec{r}', t) - \vec{u}(\vec{r}, t)$ is the deformation of the bond $\vec{\xi}$

<figure><img src="image_13.png" /><figcaption>**Fig. 8** Schematic view of the peridynamic approach: initial state ($t=0$) and deformed state ($t>0$). Definition of horizon $\delta$, neighbourhood $H_r$ of $\vec{r}$, bond $\vec{\xi}$, displacement $\vec{u}(\vec{r}, t)$ and bond deformation $\vec{\eta}$</figcaption></figure>

(Fig. 8). The force balance leads to the following integro-differential equation:

$$ \rho(\vec{r}) \ddot{\vec{u}}(\vec{r}, t) = \int_{H_r} \vec{f}(\vec{u}(\vec{r}', t) - \vec{u}(\vec{r}, t), \vec{r}' - \vec{r}) d\vec{r}' + \vec{b}(\vec{r}, t) \quad (23) $$

where $\rho$ is the density, $\vec{f}$ is the internal force term and $\vec{b}(\vec{r}, t)$ is a volume force, such as gravity. The internal force, $\vec{f}$, depends on the so-called force state on both the $\vec{r}'$ and $\vec{r}$ positions (Silling 2000; Silling et al. 2007). Among all possible state-based peridynamic models, in ordinary state-based models, forces between points are exerted in the same direction as the deformed bond $\vec{\xi} + \vec{\eta}$ (Hu et al. 2012). In such a case, the intensity of $\vec{f}(\vec{u}(\vec{r}', t) - \vec{u}(\vec{r}, t), \vec{r}' - \vec{r})$ can still depend on all bonds of both $\vec{r}'$ and $\vec{r}$. A simpler way is to use bond-based peridynamic approaches (Silling and Askari 2005; Silling et al. 2007), which leads to valuable results in various fields, such electronic devices durability (Agwai et al. 2011) and fibre-reinforced composites fractures (Hu et al. 2012), concrete failure mechanics (Gerstle et al. 2007), polycrystal fracture (Askari et al. 2008) and nanoscale rupture mechanics (Bobaru 2007; Celik et al. 2011). In such models, opposite forces of the same magnitude are exerted through a deformed bond, $\vec{\xi} + \vec{\eta}$, which has the same direction as the bond and depends only on the bond deformation. As a consequence, Eq. (23) can be rewritten in a simpler way:

$$ \rho \ddot{\vec{u}}(\vec{r}, t) = \int_{H_r} f(\|\vec{\eta} + \vec{\xi}\|, \|\vec{\xi}\|) \frac{\vec{\eta} + \vec{\xi}}{\|\vec{\eta} + \vec{\xi}\|} d\vec{r}' + \vec{b}(\vec{r}, t) \quad (24) $$

Here, we consider only linear models as a function of the bond relative elongation $s$, as follows:

$$ f(\|\vec{\eta} + \vec{\xi}\|, \|\vec{\xi}\|) = cs \quad (25) $$

$$ s = \frac{\|\vec{\eta} + \vec{\xi}\| - \|\vec{\xi}\|}{\|\vec{\xi}\|}, \quad (26) $$

where $c$ is the micromodulus (Ha and Bobaru 2011). To take crack growth into account, a critical elongation, $s_0$, is defined: when this critical value is reached, the bond is definitely broken, and the bond force exerted is set to 0 (Fig. 9).

Cracks appear as damaged zones with finite thicknesses, as peridynamics is a non-local approach. Both the elastic modulus, $E$, and the surface energy, $G$, can be deduced from the
<figure><img src="image_14.png" /><figcaption>Fig. 9 Bond force as a function of bond relative elongation, definition of critical elongation</figcaption></figure>

microscale parameters c and s₀ with the help of analytical expressions (Ha and Bobaru 2011) as follows

$$E = c\pi\delta^3(1-v)/6 \qquad (27)$$

$$G = 9E\delta s_0^2 / 4\pi \qquad (28)$$

The material is discretised with the help of a regular Cartesian grid of material points. Concretely, the system involved in the simulations is made of masses, $m_i$, initially on grid node positions $\vec{r}_i$, linked by springs above closest neighbours. The equation of motion becomes:

$$m_i \ddot{\vec{u}}_i(t) = \sum_{\vec{r}_j \in H} ks(\vec{\eta}_{ij}, \vec{\xi}_{ij}) \frac{\vec{\eta}_{ij} + \vec{\xi}_{ij}}{\|\vec{\eta}_{ij} + \vec{\xi}_{ij}\|} + \vec{b}_i(t) \quad (29)$$

where $\vec{u}_i(t) = \vec{u}(\vec{r}_i, t)$, $\vec{\xi}_{ij} = \vec{r}_j - \vec{r}_i$ is the bond between points initially located on the $\vec{r}_i$ and $\vec{r}_j$ positions, $\vec{\eta}_{ij} = \vec{u}_j - \vec{u}_i$ is the deformation of $\vec{\xi}_{ij}$, $k=c(\Delta x)^4$ is a stiffness coefficient, $\Delta x$ being grid step, and $\vec{b}_i(t)$ is the body force on the material point initially located on $\vec{r}_i$. Clearly, such an approach is similar to molecular dynamics (Seleson et al. 2009), and Eq. (29) can be integrated with the help of the well-known velocity-Verlet integration algorithm (Allen and Tildesley 1986).

To take heterogeneous mechanical properties into account, both k and s₀ depend upon the phase (i.e. cell wall or middle lamella) of points $\vec{r}_i$ and $\vec{r}_j$. Concretely, if $\varphi(\vec{r}_i)$ is a phase index, k becomes $k_{ij}=k(\varphi_{i}\varphi_{j})$ in Eq. (29), where $\varphi_i = \varphi(\vec{r}_i)$ and $\varphi_j = \varphi(\vec{r}_j)$. Critical elongation s₀ becomes $s_{oij}=s_0(\varphi_i, \varphi_j)$ in the same way. As a consequence, a stiffness $k_{\alpha\beta}$ and a critical elongation $s_{0\alpha\beta}$ can be attributed to each {$\alpha, \beta$} phase-phase couple. Bulk phase elastic moduli $E_\alpha$ can be deduced from $k_{\alpha\alpha}$

values following Eq. (27), bulk phase fracture energies $G_\alpha$ can be deduced from $s_{0\alpha\alpha}$ following Eq. (28) and interface fracture energies $G_{\alpha\beta}$ can be deduced from $s_{0\alpha\beta}$ when $a \neq \beta$.

# 4 Results

## 4.1 Thermal and mass diffusivity

The LBM was used to compute the equivalent thermal conductivity and the mass diffusivity of the cellular structure of spruce, which was obtained in Section 2. The simulation was performed in the two transverse directions of wood (radial and tangential), which means a flux along the horizontal and vertical directions, respectively, in the images depicted in Fig. 5.

The solid fraction $\varepsilon_s$ was equal to 0.40 and 0.78 for earlywood and latewood, respectively. The thermal conductivities of the solid and the gas were set to $\lambda_s=1$ W m⁻¹ K⁻¹ and $\lambda_{air}=0.023$ W m⁻¹ K⁻¹, respectively. For mass diffusion, we used the dimensionless mass diffusivity f, which accounts for the diffusion resistance relative to binary diffusion in air (Perré and Turner 2001b). Therefore, this value ranges from 0 to 1. The values of f used for the solid and gaseous phases were $f_s=0.004$ and $f_{air}=1$, respectively. These values are representative of bound water diffusion in wood at a moisture content value of ca. 12 % (Siau 1984). A full LBM computational run on a grid of about 200×200 points requires a couple of hours (2 to 3) using an Intel I7 processor at 3.4 GHz. Note that 3–10 h were required for simulating mass diffusion because the connected phase (the solid) phase has a low diffusivity in this case.

The computed equivalent values were compared to the lowest and highest possible values (the heterogeneous phases

<figure><img src="image_15.png" /><figcaption>Fig. 10 The local field obtained at steady state in the case of thermal conductivity. Temperature values are in colour and heat flux is shown as arrows</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. 11** Dimensionless equivalent conductivity, $\lambda_{eq}$, versus solid fraction, $\varepsilon_s$, for radial (R) and tangential (T) directions of spruce early- and latewood. Series (dashed line) and parallel (solid line) theoretical models are represented</figcaption></figure>

placed in series or in parallel, respectively). The respective solutions read as

$$ \Gamma_{\text{series}} = \frac{1}{\frac{1-\varepsilon_s}{\Gamma_{\text{air}}} + \frac{\varepsilon_s}{\Gamma_s}} \quad (30) $$

$$ \Gamma_{\text{parallel}} = (1-\varepsilon_s)\Gamma_{\text{air}} + \varepsilon_s\Gamma_s \quad (31) $$

where $\Gamma$ represents $\lambda$ or $f$ for thermal conduction and mass diffusion, respectively.

Figure 10 depicts an example of a local field obtained in the case of thermal conduction. In this case, it is obvious that the

heat flux takes advantage of the connected and conductive phase (the solid phase), avoiding the inclusions of low conductivity (the lumens).

Because this morphology is favourable to heat conduction, it is not surprising to observe that the predicted macroscopic values are quite close to the parallel model (Fig. 11). As the cell walls are rather aligned along the radial direction in earlywood, a consequence of cambial cell division in trees, the thermal conductivity is higher in the radial direction than in the tangential direction. In contrast, the tracheids are flattened in latewood, which forms lumens having a large tangential extension. These inclusions, with a low conductivity, block

<figure><img src="image_17.png" /><figcaption>**Fig. 12** Dimensionless equivalent diffusivity $f$ versus solid fraction, $\varepsilon_s$, for radial (R) and tangential (T) directions of spruce early- and latewood. Series (dashed line) and parallel (solid line) theoretical models are represented. The inset depicts the radial and tangential diffusivity for the latewood on an enlarged scale</figcaption></figure>
the thermal flux mainly in the radial direction, which explains the inverted anisotropic ratio.

In the case of mass diffusivity, the situation is opposite because the connected phase, the solid phase, has a low diffusivity. As a consequence, the macroscopic value is very close to the series model (Fig. 12). Again, for the same morphological reasons, the anisotropic ratio is inverted between earlywood and latewood. These results confirm those obtained by Perré and Turner (2001b) for a simplified geometry of cellular structure (Perré and Turner 2001a).

## 4.2 Mechanical behaviour

The mechanical behaviour of the three tissues (Fig. 5) was simulated using MPM in both the radial and tangential directions. One run basically consists of a compression test with large deformations. This method was implemented in a custom code written in FORTRAN95 and parallelised using Message Passing Interface (MPI) routines. A full computational run on a grid of about 200,000 points requires a couple of hours (2 to 3) using 4 cores on an Intel I7 processor at 3 GHz. To allow the REV to pave the plane, the lateral faces of the domain are forced to stay straight, and an iterative algorithm was derived to keep the average lateral force at zero. This strategy assumes that the REV has two planes of symmetry (Farruggia 1998). The local cell wall properties were deduced from the paper by Neagu and Gamstedt (2007). A blocked rotation was applied to the values of this work to account for an AMF of ca. 20°. We finally used 8 GPa along the tangential direction (parallel to the lumen contour) and 6.4 GPa across the cell wall.

Table 1 summarises the stiffness values obtained by the MPM simulation. The stiffness was obtained by a linear regression of the stress-strain curve in the elastic domain. As a simple rule based on the curve shape, we used the range of deformation [0 %, 3 %] for all tests. Importantly, one can observe that by accounting for the cellular morphology, the macroscopic properties are in rather good agreement with those in the literature (Farruggia and Perré 2000; Perré et al. 2013). This agreement is quite perfect for poplar, where the approach is really predictive, both in terms of its value and

<table><caption>**Table 1** Summary of the stiffness computed by MPM modelling for the three tissues selected in the present work</caption><thead><tr><th rowspan="2">Tissue</th><th rowspan="2">Solid fraction (%)</th><th colspan="2">Stiffness (MPa)</th></tr><tr><th>Tangential</th><th>Radial</th></tr></thead><tbody><tr><td>Spruce, earlywood</td><td>40</td><td>416</td><td>1090</td></tr><tr><td>Spruce, latewood</td><td>78</td><td>3470</td><td>2430</td></tr><tr><td>Poplar</td><td>38</td><td>364</td><td>811</td></tr></tbody></table>

anisotropic ratio. In the case of spruce earlywood, the prediction quality is very good as well (the cellular morphology used here corresponds more closely to that of intermediate wood, as a solid fraction of 40 % would give a density of 600 kg m⁻³). In the case of latewood, no experimental data are really available in the very last layer of latewood, due to experimental issues. Therefore, the model should be considered as a way to generate new data. Among the interesting outputs is the inversion of the anisotropy: due to the flattened cells, the tangential stiffness is higher than the radial stiffness.

Figures 13 and 14 depict the main results of the numerical test performed for spruce earlywood and poplar, respectively. Thanks to the ability of MPM to deal with large deformations and contacts, the model has the remarkable capacity to predict the entire compression curve: initial elastic behaviour, compression plateau and final densification. Spruce depicts a progressive collapse of entire tracheid layers, starting with the weakest ones, while in poplar, vessels act as weak elements allowing the cellular structure to deform before the cell lumen eventually collapses.

However, one should emphasise that the stress levels and the domain of elasticity are higher than expected. Several effects are likely to explain these differences, including boundary conditions, cell wall behaviour (assumed to be elastic here) and refinement of the background grid. All these effects will be presented in detail in a future paper that will focus on MPM.

<figure><img src="image_18.png" /><figcaption>**Fig. 13** Stress-strain curves obtained for spruce earlywood in the radial and tangential directions. Computational images of the cellular deformation for the radial test were selected at strain levels of 5, 12, 25 and 48 %</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig. 14** Stress-strain curves obtained for poplar in the radial and tangential directions. Images of the cellular deformation for the radial test were selected at strain levels of 10, 20 and 42 %</figcaption></figure>

## 4.3 Fragmentation

Peridynamic studies of quasi-static crack growth are quite rare (Kilic and Madenci 2010). To perform quasi-static numerical tests, an equilibrium search approach is required. To perform damping of mechanical energy in the material and reach an equilibrium state, we applied a viscous dissipation as a body force to each material point. After each strain increment, critical bonds are broken and a new equilibrium search is performed, if needed, until a stable mechanical equilibrium is reached (Fig. 15).

This algorithm was implemented in parallel in a custom code written in FORTRAN 95 and using Message Passing Interface (MPI) routines. A typical computing time is 8 h on a 16-core server.

In the present application, we used a horizon, $\delta=4\Delta x$, where $\Delta x$ is the grid step, which is a good compromise

<figure><img src="image_20.png" /><figcaption>**Fig. 15** Algorithm used in the present study to perform quasi-static tests. After each strain increment, successive equilibrium searches and bond break steps are performed until an equilibrium search without a critical bond is reached</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 16** Snapshots of successive steps of a peridynamic simulation of spruce earlywood failure under tension (colour online). The crack is in black. Colours are related to tensile stress: green when stress is close to 0, blue below 0 and red above 0</figcaption></figure>

between precision and numerical weight. Two phases have to be distinguished: the cell wall phase and middle lamella. The tenacity, $K = \sqrt{EG}$, is defined for each phase: $K_{cw}$ for the cell wall and $K_{ml}$ for lamella. Here, we fixed the ratio between tenacities as $K_{cw}/K_{ml} = 4$. A pre-notch was performed on one side of the spruce latewood sample and the tensile fracture test was performed. The same test was performed for the earlywood and latewood of spruce.

As we can see in Fig. 16, a crack grew from the pre-notch and gradually extended to the other side of the sample. The crack remained inside the middle lamella, and cell walls were not damaged during the failure test, as expected according to the tenacity ratio. The same behaviour was observed in the case of spruce latewood (Fig. 17). In both cases, the crack remained within the most fragile phase, the middle lamella.

Such a test clearly shows that the peridynamic approach is a convenient framework to simulate crack growth in wood, despite the complex structure and mechanical properties inside this material.

<figure><img src="image_22.png" /><figcaption>**Fig. 17** Snapshots of successive steps of a peridynamic simulation of spruce latewood failure under tension (colour online). The crack is in black. Colours are related to tensile stress: green when stress is close to 0, blue for negative values and red for positive values</figcaption></figure>
## 5 Conclusion

This paper presents a comprehensive strategy to predict different wood properties from anatomical images. The starting point is an image-based representation that is able to benefit from any present and future imaging tools, either in 2D or in 3D. This representation step accounts for the local anisotropy and heterogeneity of the cell wall. Then, several meshless methods are proposed to compute the properties of wood tissues from its cellular structure. The emphasis here was on meshless methods, including the Lattice Boltzmann Method, Material Point Method and peridynamic method, which are able to account quite easily for any complex geometries and to predict thermal and mass diffusivities, stiffness and fracture, respectively.

A selected set of computational results proved the predictive ability of this modelling strategy and its potential to predict properties that would be difficult, if not impossible, to measure.

Further studies are in progress to extend this modelling approach in 3D and to extract general trends by comparing various anatomical patterns. The extension from 2D to 3D is today perfectly possible, namely thanks to µ-tomography available at synchrotron facilities. However, it requires much higher computational resources for three major steps: image processing, the computational part itself and post-processing.

**Acknowledgments** The authors would like to thank the French Agence Nationale de la Recherche (ANR, project ANR 10 - HABISOL-005-03) for its participation in founding this research work.

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

## References

Abbasi AR (2013) Multiscale poroelastic model: bridging the gap from cellular to macroscopic scale (Doctoral dissertation, Diss., Eidgenössische Technische Hochschule ETH Zürich, Nr. 20821, 2013)

Affes R, Delenne JY, Monerie Y, Radjaï F, Topin V (2012) Tensile strength and fracture of cemented granular aggregates. Eur Phys J E 35:117. doi:10.1140/epje/i2012-12117-7

Agwai A, Guven I, Madenci E (2011) Crack propagation in multilayer thin-film structures of electronic packages using the peridynamic theory. MicroelecReliab 51:2298-2305. doi:10.1016/j.microrel. 2011.05.011

Allen MP, Tildesley DJ (1986) Computer simulation of liquids. Oxford University Press, American Mathematical Society, 57(195):442-444. doi:10.2307/2938686

Almeida G, Brito J, Perré P (2009) Changes in wood-water relationship due to heat treatment assessed on micro-samples of three Eucalyptus species. Holzforchung 63:80-88. doi:10.1515/HF.2009.026

Almeida G, Huber F, Perré P (2014) Free shrinkage of wood determined at the cellular level using an environmental scanning electron microscope. Maderas Cienc Tecnol 16:187-198. doi:10.4067/S0718- 221X2014005000015

Askari E, Bobaru F, Lehoucq RB, Parks ML, Silling SA, Weckner O (2008) Peridynamics for multiscale materials modelling. J Phys Conf Ser 125:012078. doi:10.1088/1742-6596/125/1/012078

Bao J, Yuan P, Schaefer L (2008) A mass conserving boundary condition for the lattice Boltzmann equation method. J Comput Phys 227: 8472-8487. doi:10.1016/j.jcp.2008.06.003

Barber NF (1968) A theoretical model of shrinking wood. Holzforschung 22:97-103

Barber NF, Meylan BA (1964) The anisotropic shrinkage of wood: a theoretical model. Holzforschung 18:146-156

Bardenhagen S, Kober E (2004) The generalized interpolation material point method. Comput Model Eng Sci 5:477-495. doi:10.3970/cmes.2004.005.477

Barkas WA (1941) The influence of ray cells on the shrinkage of wood. Trans Faraday Soc 37:535-548

Belytschko T, Krongauz Y, Organ D, Fleming M, Krysl P (1996) Meshless methods: an overview and recent developments. Comput Method Appl Mech Eng 139:3-48

Bergander A, Salmén L (2000) The transverse elastic modulus of the native wood fibre wall. J Pulp Pap Sci 26:234-238

Bobaru F (2007) Influence of van der Waals forces on increasing the strength and toughness in dynamic fracture of nanofibre networks: a peridynamic approach. Modelling Simul Mater Sci Eng 15:397-417. doi:10.1088/0965-0393/15/5/002

Bosshard HH (1984) Holzkunde Band 2 Zür Biologie, Physik und Chemie des Holzes Birkhäuser Verlag Basel

Boutelje J (1962) The relationship of structure to transverse anisotropy in wood with reference to shrinkage and elasticity. Holzforshung 16: 33-46

Cabrolier P (2012) Caractérisation des propriétés structurales et mécaniques des composantes pariétales du bois à l'échelle du tissu. PhD manuscript 254 pages AgroParisTech Nancy

Cave ID (1966) Theory of X-ray measurement of microfibril angle in wood. Forest Prod J 16:37-42

Celik E, Guven I, Madenci E (2011) Simulations of nanowire bend tests for extracting mechanical properties. Theor Appl Fract Mech 55: 185-191. doi:10.1016/j.tafmec.2011.07.002

Charras GT, Guldberg RE (2000) Improving the local solution accuracy of large-scale digital image-based finite element analyses. J Biomech 33:255-259. doi:10.1016/S0021-9290(99)00141-4

De la Osa MR, Estevez R, Olagnon C, Chevalier J, Vignoud L, Tallaron C (2009) Cohesive zone model and slow crack growth in ceramic polycrystals. Int J Fract 158:157-167. doi:10.1007/s10704-009- 9346-3

Farruggia F (1998) Détermination du comportement élastique d'un ensemble de fibres de bois à partir de son organisation cellulaire et d'essais mécaniques sous microscope. PhD thesis ENGREF Nancy

Farruggia F, Perré P (2000) Microscopic tensile tests in the transverse plane of earlywood and latewood parts of spruce. Wood Sci Tech 34:65-82. doi:10.1007/s002260000034

Frank X, Perré P (2010) The potential of meshless methods to address physical and mechanical phenomena involved during drying at the pore level. Drying Technol 28:932-943. doi:10.1080/07373937. 2010.497077

Frank X, Almeida G, Perré P (2010) Multiphase flow in the vascular system of wood: from microscopic excursion to 3-D Lattice Boltzmann experiments. J Multiphase flow 36:599-607. doi:10. 1016/j.ijmultiphaseflow.2010.04.006

Gerstle W, Sau N, Silling SA (2007) Peridynamic modelling of concrete structures. Nuclear Eng Design 237:1250-1258. doi:10.1016/j. nucengdes.2006.10.002

<figure><img src="image_23.png" /></figure>
Gibson LJ, Ashby MF (1988) Cellular solids, structure and properties. Oxford: Pergamon Press, p 357. doi:10.1002/adv.1989.060090207

Gierlinger N, Schwanninger M (2006) Chemical imaging of poplar wood cell walls by confocal Raman microscopy. Plant Physiol 40:1246–1254. doi:10.1104/pp. 105.066993.1246

Gillis PP (1972) Orthotropic elastic constant of wood. Wood Sci Tech 6: 138–156

Guilkey JE, Hoying JB, Weiss JA (2006) Computational modeling of multicellular constructs with the material point method. J Biomech 39:2074–2086. doi:10.1016/j.jbiomech.2005.06.017

Guitard D, El Amri F (1987) La fraction volumique en rayons ligneux comme paramètre explicatif de la variabilité de l'anisotropie élastique du matériau bois. ARBOLOR Actes du 2nd colloque Sciences et Industrie du Bois 405–412

Ha YD, Bobaru F (2011) Characteristics of dynamic brittle fracture captured with peridynamics. Eng Fract Mech 78:1156–1168. doi:10.1016/j.engfracmech.2010.11.020

Harris JM, Meylan BA (1965) The influence of microfibril angle on longitudinal and tangential shrinkage in *Pinus radiata*. Holzforschung 19:144–153

Higuera F (1990) Lattice gas method based on the Chapman-Enskog expansion. Phys Fluids A 2, 1049. doi:10.1063/1.857645

Hofstetter K, Hellmich C, Eberhardsteiner J (2005) Development and experimental validation of a continuum micromechanics model for the elasticity of wood. Eur J Mech A/Solid 24:1030–1053. doi:10.1016/j.euromechsol.2005.05.006

Holmberg S, Persson K, Petersson H (1999) Nonlinear mechanical behaviour and analysis of wood and fibre material. Comput Struct 72: 459–480. doi:10.1016/S0045-7949(98)00331-9

Hu W, Ha YD, Bobaru F (2012) Peridynamic model for dynamic fracture in unidirectional fiber-reinforced composites. Comput Meth Appl Mech Eng 217–220:247–261. doi:10.1016/j.cma.2012.01.016

Itakura M, Kaburaki H, Arakawa C (2005) Branching mechanism of intergranular crack propagation in three dimensions. Phys Rev E 71:055102(R). doi:10.1103/PhysRevE.71.055102

Keller R, Thiercelin F (1975) Influence des gros rayons ligneux sur quelques propriétés du bois de hêtre. Ann Sci Forest 32:113–129

Kelsey K (1963) A critical review of the relationship between the shrinkage and structure of wood, Division of Forest products technological paper n° 28 CSIRO Melbourne

Kilic B, Madenci E (2010) An adaptive dynamic relaxation method for quasi-static simulations using the peridynamic theory. Theor Appl Fract Mech 53:194–204. doi:10.1016/j.tafmec.2010.08.001

Kitsunezaki S (2013) Cracking condition of cohesionless porous materials in drying processes. Phys Rev E 87:052805. doi:10.1103/PhysRevE.87.052805

Kollmann FP, Côté WA (1968) Principles of wood science and technology, vol 1. Springer, Solid Wood

Koponen S, Toratti T, Kanerva P (1991) Modelling elastic and shrinkage properties of wood based on cell structure. Wood Sci Technol 25: 25–32

Kwon G, Chae SW, Lee KJ (2003) Automatic generation of tetrahedral meshes from medical images. Comput Struct 81:765–775. doi:10.1016/S0045-7949(02)00406-6

Lallemand P, Luo LS (2003) Lattice Boltzmann method for moving boundaries. J Comput Phys 184:406–421. doi:10.1016/S0021-9991(02)00022-0

Lanvermann C (2014) Sorption and swelling within growth rings of Norway spruce and implications on the macroscopic scale (Doctoral dissertation, Diss., Eidgenössische Technische Hochschule ETH Zürich, Nr. 21761, 2014)

Mariaux A (1989) La section transversale de fibre observée avant et après séchage sur bois massif. Bois et Forêts de Tropiques 221:65–76

Mariaux A, Narboni P (1978) Anisotropie du retrait et structure du bois: essai d'approche statistique. Bois et Forets de Tropiques 178:36–44

Masseran C, Mariaux A (1985) Anisotropie de retrait et structure du bois. Recherche de l'influence des caractères morphologiques transverse des fibres. Bois et Forêts de Tropiques 209:35–47

Meylan BA, Probine MC (1969) Microfibril angle as a parameter in timber quality assessment. Forest Prod J 19:30–34

Mohamad AA (2011) Lattice Boltzmann method: fundamentals and engineering applications with computer codes. Springer, London

Nairn JA (2006) Numerical simulations of transverse compression and densification in wood. Wood and Fiber Sci 38:576–591

Neagu RC, Gamstedt EK (2007) Modelling of effects of ultrastructural morphology on the hydroelastic properties of wood fibres. J Mater Sci 42:10254–10274. doi:10.1007/s10853-006-1199-9

Perré P (2001) Wood as a multi-scale porous medium: observation, experiment, and modelling, keynote lecture. 1st International conference of the European Society for wood mechanics Lausanne Switzerland 403–422

Perré P (2005) MeshPore: a software able to apply image-based meshing techniques to anisotropic and heterogeneous porous media. Drying Technol J 23:1993–2006. doi:10.1080/07373930500210432

Perré P (2007) Experimental device for the accurate determination of wood-water relations on micro-samples. Holzforschung 61:419–429. doi:10.1515/HF.2007.075

Perré P (2011) A review of modern computational and experimental tools relevant to the field of drying. Drying Technol 29:1529–1541. doi:10.1080/07373937.2011.580872

Perré P, Huber F (2007) Measurement of free shrinkage at the tissue level using an optical microscope with an immersion objective: results obtained for Douglas fir (*Pseudotsuga menziesii*) and spruce (*Picea abies*). Ann For Sci 64:255–265. doi:10.1051/forest:2007003

Perré P, Turner I (2001a) Determination of the material property variations across the growth ring of softwood for use in a heterogeneous drying model. Part I: capillary pressure, tracheid model and absolute permeability. Holzforschung 55:318–323. doi:10.1515/HF.2001.052

Perré P, Turner I (2001b) Determination of the material property variations across the growth ring of softwood for use in a heterogeneous drying model. Part II: use of homogenisation to predict bound water diffusivity and thermal conductivity. Holzforschung 55:417–425. doi:10.1515/HF.2001.069

Perré P, Turner IW (2002) A heterogeneous wood drying computational model that accounts for material property variation across growth rings. Chem Eng J 86:117–131. doi:10.1016/S1385-8947(01)00270-4

Perré P, Turner IW (2008) A mesoscopic drying model applied to the growth rings of softwood: mesh generation and simulation results, Maderas. Cienc Tecnol 10:251–274

Perré P, Almeida G, Frank X (2012) MPM modelling of the cellular collapse of bio-products due to capillary forces, proceedings of the 18th Int. Drying Symposium 5 pages Xiamen China

Perré P, Dinh T, Assor C, Frank X, Pilate G (2013) Stiffness of normal, opposite, and tension poplar wood determined using micro-samples in the three material directions. Wood Sci Tech 47:481–498. doi:10.1007/s00226-012-0511-x

Ruiz G, Oritz M, Pandolfi A (2000) Three-dimensional simulation of the dynamic Brazilian tests on concrete cylinders. Int J Num Methods Eng 48:963–994

Salmen L (2004) Micromechanical understanding of the cell-wall structure. C R Biolog 327:873–880. doi:10.1016/j.crvi.2004.03.010

Sanchez-Palencia E (1980) Non-homogeneous media and vibration theory. Lecture Notes in Physics 127 Springer Verlag

Sancho JM, Planas J, Cendón DA, Reyes E, Gálvez JC (2007) An embedded crack model for finite element analysis of concrete fracture. Eng Fract Mech 74:75–86. doi:10.1016/j.engfracmech.2006.01.015

<figure><img src="image_24.png" /></figure>
Sedighi-Gilani M, Navi P (2007) Experimental observations and micromechanical modeling of successive-damaging phenomenon in wood cells' tensile behavior. Wood Sci Tech 41:69-85. doi:10.1007/s00226-006-0094-5

Seleson P, Parks ML, Gunzburger M, Lehoucq RB (2009) Peridynamics as an upscaling of molecular dynamics. Multiscale Model Simul 8: 204-227. doi:10.1137/09074807X

Siau JF (1984) Transport processes in wood. Springer, New York, p 245

Silling SA (2000) Reformulation of elasticity theory for discontinuities and long-range forces. J Mech Phys Solids 48:175-209. doi:10.1016/S0022-5096(99)00029-0

Silling SA, Askari E (2005) A meshfree method based on the peridynamic model of solid mechanics. Comput Struct 83:1526-1534. doi:10.1016/j.compstruc.2004.11.026

Silling SA, Epton M, Weckner O, Xu J, Askari E (2007) Peridynamic states and constitutive modeling. J Elasticity 88:151-184. doi:10.1007/s10659-007-9125-1

Succi S (2001) The lattice Boltzmann equation for fluid dynamics and beyond. Clarendon Press Oxford, New York, p 308. doi:10.1063/1.537916

Sulsky D, Schreyer L (2004) MPM simulation of dynamic material failure with a decohesion constitutive model. Eur J

Mech A/Solid 23:423-445. doi:10.1016/j.euromechsol.2004.02.007

Sulsky D, Chen Z, Schreyer H (1994) A particle method for history-dependent materials. Comp Methods Appl Mech Eng 118:179-196. doi:10.1016/0045-7825(94)90112-0

Suquet PM (1985) Element of homogenization for inelastic solid mechanics in homogenization techniques for composite media. Lecture Notes in Physics 272 Sanchez-Palencia and Zaoui Ed Springer-Verlag

Topin V, Delenne JY, Radjaï F, Brendel L, Mabille F (2007) Strength and failure of cemented granular matter. Eur Phys J E 23:413-429. doi:10.1140/epje/i2007-10201-9

Ulrich D, van Rietbergen B, Weinans H, Rüegsegger P (1998) Finite element analysis of trabecular bone structure: a comparison of image-based meshing techniques. J Biomech 31:1187-1192

Watanabe U, Norimoto M (1996) Shrinkage and elasticity of normal and compression woods in conifers. Mokuzai Gakkaishi 42:651-658

Zavattieri PD, Raghuram PV, Espinosa HD (2001) A computational model of ceramic microstructures subjected to multi-axial dynamic loading. J Mech Phys Solids 49:2768. doi:10.1016/S0022-5096(00)00028-4

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>