Skip to main content

Volume 73 Supplement 1

50 years Annals of Forest Science

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

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 latewood 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.

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 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−3 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 variability 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 macroscopic behaviours. It is now well established that three spatial 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 pattern, in which all the anatomical tissues are organised in proportion 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 simplicity, due to the fact that all solid components act in series in the longitudinal directions, allowed analytical models to be proposed in the same period (Barber 1968).

The understanding and prediction of structure-property relationships 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-layered 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 observed 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 explain 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 structures, is a deterministic approach that includes several steps (Fig. 1):

Fig. 1
figure 1

Upscaling in material sciences: principle

  1. 1.

    Representation: choice of the representative elementary volume (REV), also called the Unit Cell. This REV should be defined in a suitable way for subsequent calculation (finite element mesh, collection of material points…).

  2. 2.

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

  3. 3.

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

  4. 4.

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

  5. 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 modelling, 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).

Fig. 2
figure 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

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.

Fig. 3
figure 3

Using the neighbouring contours, an algorithm was developed to propagate the material directions to each material point of the solid phase (left)

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.

Fig. 4
figure 4

Generation of internal contours and allocation of two solid phase types (blue = lumens, red = secondary wall and yellow = middle lamella)

Fig. 5
figure 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)

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\left(\overrightarrow{x}+{\overrightarrow{c}}_i\updelta t,t+\updelta t\right)-{f}_i\left(\overrightarrow{x},t\right)={\varOmega}_i\left(\overrightarrow{x},t\right) $$
(1)

where \( {\overrightarrow{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 \( {\overrightarrow{c}}_i \) are the possible velocities a particle can move from a lattice site to a nearest-neighbouring site at each time step t + δt. For simplicity, we used a “single relaxation time” form for the collision term:

$$ {\varOmega}_i\left(\overrightarrow{x},t\right)=-\frac{1}{\tau}\left[{f}_i\left(\overrightarrow{x},t\right)-{f}_i^{\mathrm{eq}}\left(\overrightarrow{x},t\right)\right], $$
(2)

where τ is the dimensionless relaxation time and \( {f}_i^{\mathrm{eq}}\left(\overrightarrow{x},t\right) \) is the equilibrium distribution function at site \( \overrightarrow{x} \) and time t. The above equation is well known as the Bhatnagar-Gross-Krook (BGK) approximation for the collision operator \( {\varOmega}_i\left(\overrightarrow{x},t\right) \). 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, cii = {0, 1, … 8}. These velocities allow the particles to jump to neighbouring nodes in one lattice time step (Fig. 6). They are defined as follows:

$$ {\overrightarrow{c}}_i\left\{\begin{array}{ll}\left(0,0\right)c\hfill & i=0\hfill \\ {}\left(\pm 1,0\right)c,\left(0,\pm 1\right)c\hfill & i=1,2,3,4\hfill \\ {}\left(\pm 1,\pm 1\right)c\hfill & i=5,6,7,8\hfill \end{array}\right. $$
(3)

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

Fig. 6
figure 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 fi)

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^{\mathrm{eq}}\left(x,y,t\right)={\omega}_i\varPsi \left(x,y,t\right), $$
(4)

where ωi are the weighting factors and Ψ(x, y, t) is the normalised macroscopic variable (temperature field θ or mass X (see Section 4), such as (x, y, t)  [0, 1]. In the D2Q9 model, the weighting factors are given as (Mohamad 2011):

$$ {\omega}_i\left\{\begin{array}{ll}4/9\hfill & i=0\hfill \\ {}1/9\hfill & i=1,2,3,4\hfill \\ {}1/36\hfill & i=5,6,7,8\hfill \end{array}\right. $$
(5)

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

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

was solved in a 2D rectangular porous medium, [x1,xn]×[y1,ym], 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: Ψ(x, y, 0) = 0. Along the x-direction, Dirichlet boundary conditions were applied: Ψ(x1, y, t) = 1 and Ψ(xn, y, t) = 0. Adiabatic conditions were set along the y-direction: φ(x, y, t) = 0 at y = y1 and y = yn.

In the literature (Higuera 1990), it has been proved that the diffusion equation can be derived using the Chapman-Enskog 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{\updelta t}{2}\right){c}^2, $$
(7)

where τ has a dimension of time, expressed in seconds. The macroscopic normalised quantities, Ψ and related flux φ, are obtained through moment summations in the velocity space:

$$ \varPsi \left(x,y,t\right)={\displaystyle \sum_{i=0}^8{f}_i\left(x,y,t\right)} $$
(8)
$$ \overrightarrow{\varphi}\left(x,y,t\right)={\displaystyle \sum_{i=0}^8{\overrightarrow{c}}_i{f}_i\left(x,y,t\right)} $$
(9)

The normalised flux has the following form:

$$ \overrightarrow{\varPhi}\left(x,y,t\right)=\frac{\left(\tau -\frac{1}{2}\right)}{\tau}\overrightarrow{\varphi}\left(x,y,t\right) $$
(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 (f1, f5, f8) and (f3, f6, f7) associated to velocities (c1, c5, c8) and (c3, c6, c7) are known at the lattice faces x = x1 and x = xn, respectively. On the contrary, (f3, f6, f7) at x = x1 and (f1, f5, f8) at x = xn 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{array}{ll}{f}_i\left({x}_1,y,t\right)={\omega}_i{\varPsi}_1\hfill & \forall i=\left\{1,5,8\right\}\hfill \\ {}{f}_i\left({x}_n,y,t\right)={\omega}_i{\varPsi}_n\hfill & \forall i=\left\{3,6,7\right\}\hfill \end{array} $$
(11)

where Ψ1 and Ψn are fictive variables obtained from balance equations to ensure the desired Dirichlet value. For example, for the plane (x = x1), we use

$$ {\varPsi}_1=\frac{1}{{\displaystyle {\sum}_i{\omega}_i}}\left[\varPsi \left({x}_1,y,t\right)-\left({\displaystyle \sum_{k=0}^8{\tilde{f}}_k\left({x}_1,y,t\right)-}{\displaystyle \sum_i{\tilde{f}}_k\left({x}_1,y,t\right)}\right)\right]\kern1em {\forall}_i=\left\{1,5,8\right\} $$
(12)

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 (f2, f5, f6) at y = y1 and (f4, f7, f8) at y = ym can be obtained using as moving boundary conditions (see Lallemand and Luo 2003 for more details) such as

$$ \begin{array}{l}{f}_2\left(x,{y}_1,t\right)={f}_4\left(x,{y}_1,t\right),{f}_5\left(x,{y}_1,t\right)={f}_8\left(x,{y}_1,t\right),{f}_6\left(x,{y}_1,t\right)={f}_7\left(x,{y}_1,t\right)\hfill \\ {},{f}_4\left(x,{y}_m,t\right)={f}_2\left(x,{y}_m,t\right),{f}_7\left(x,{y}_m,t\right)={f}_6\left(x,{y}_m,t\right),{f}_8\left(x,{y}_m,t\right)={f}_5\left(x,{y}_m,t\right)\hfill \end{array} $$
(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, τ1 and τ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 [x1, xn] × [y1, ym] is, therefore, divided into two sub-sets: the cell lumens, Σ1, and the solid phase, Σ2, which are related to the macroscopic diffusion coefficients D1 and D2, respectively.

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

  1. 1.

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

  2. 2.

    Stream the density distribution populations.

  3. 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 (105–106) 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… Np. 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, mp, 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 \left(\overrightarrow{x},t\right)={\displaystyle \sum_{p=1}^{N_p}{m}_p\updelta \left(\overrightarrow{x}-{\overrightarrow{X}}_p(t)\right)} $$
(14)
Fig. 7
figure 7

Principle of the Material Point Method (adapted from Sulsky et al. 1994)

The mp 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… Ng 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, Ni, are used for this purpose:

$$ \begin{array}{l}{S}_{pi}^k={N}_i\left({\overrightarrow{X}}_p^k\right)\hfill \\ {}{G}_{pi}^{-k}={\left.\nabla \right|}_{X_p^k}{N}_i\hfill \end{array} $$
(15)

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

$$ \begin{array}{l}{\overrightarrow{f}}_i^k={\overrightarrow{f}}_i^{\mathrm{int},k}+{\overrightarrow{f}}_i^{\mathrm{ext},k}\hfill \\ {}\mathrm{where}\kern1em {\overrightarrow{f}}_i^{\mathrm{int},k}=-{\displaystyle \sum_{p=1}^{N_p}{m}_p\left({\overline{G}}_{ip}^{T,k}{\overline{\overline{\sigma}}}_p^{S,k}\right)}\hfill \end{array} $$
(16)

In Eq. (16), \( {\overline{\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={\displaystyle \sum_{p=1}^{N_p}{m}_p\left({S}_{ip}^{T,k}{S}_{pj}^k\right)} $$
(17)

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

$$ {\displaystyle \sum_{j=1}^{N_g}{m}_{ij}^k{\overrightarrow{a}}_j^k={\overrightarrow{f}}_i^k} $$
(18)

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

$$ {\displaystyle \sum_{j=1}^{N_g}{m}_{ij}^k{\overrightarrow{v}}_j^k={\displaystyle \sum_{p=1}^{N_p}{m}_p{\overrightarrow{V}}_p^k{N}_i\left({X}_p^k\right)}} $$
(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{array}{l}{\overrightarrow{v}}_i^L={\overrightarrow{v}}_i^k+\varDelta t\kern0.5em {\overrightarrow{a}}_i^k\hfill \\ {}{\overrightarrow{x}}_i^L={\overrightarrow{x}}_i^k+\varDelta t\kern0.5em {\overrightarrow{v}}_i^L\hfill \end{array} $$
(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{array}{l}{\overrightarrow{V}}_p^L={\overrightarrow{V}}_p^k+\varDelta t{\displaystyle \sum_{i=1}^{N_g}{\overrightarrow{a}}_i^k{N}_i\left({X}_p^k\right)}\hfill \\ {}{\overrightarrow{X}}_p^L={\overrightarrow{X}}_p^k+\varDelta t{\displaystyle \sum_{i=1}^{N_g}{\overrightarrow{v}}_i^L{N}_i\left({X}_p^k\right)}\hfill \end{array} $$
(21)

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

$$ \varDelta {\overline{\overline{\varepsilon}}}_p=\frac{\varDelta t}{2}{\displaystyle \sum_{i=1}^{N_g}\left\{{\overline{G}}_{ip}^k{\overrightarrow{v}}_i^L+{\left({\overline{G}}_{ip}^k{\overrightarrow{v}}_i^L\right)}^T\right\}} $$
(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=\left\{\left.{\overrightarrow{r}}^{\prime}\right|\left\Vert \overrightarrow{r}-{\overrightarrow{r}}^{\prime}\right\Vert <\updelta \right\} \), where \( \overrightarrow{r} \) is a point position at t = 0, \( {\overrightarrow{r}}^{\prime } \) are neighbouring positions and δ is a cutoff distance called the horizon (Fig. 8).

Fig. 8
figure 8

Schematic view of the peridynamic approach: initial state (t = 0) and deformed state (t > 0). Definition of horizon δ, neighbourhood Hr of \( \overrightarrow{r} \), bond \( \overrightarrow{\xi} \), displacement \( \overrightarrow{u}\left(\overrightarrow{r},t\right) \) and bond deformation \( \overrightarrow{\eta} \)

Such a paradigm is very different from classical continuous matter mechanics, in which only short-range forces are involved. A bond is defined as \( \overrightarrow{\xi}={\overrightarrow{r}}^{\prime }-\overrightarrow{r} \) when \( {\overrightarrow{r}}^{\prime } \) is in the neighbourhood of \( \overrightarrow{r} \), \( \overrightarrow{u}\left(\overrightarrow{r},t\right) \) is the displacement and \( \overrightarrow{\eta}=\overrightarrow{u}\left({\overrightarrow{r}}^{\prime },t\right)-\overrightarrow{u}\left(\overrightarrow{r},t\right) \) is the deformation of the bond \( \overrightarrow{\xi} \) (Fig. 8). The force balance leads to the following integro-differential equation:

$$ \rho \left(\overrightarrow{r}\right)\overset{..}{\overrightarrow{u}}\left(\overrightarrow{r},t\right)={\displaystyle \underset{H_r}{\int}\overrightarrow{f}\left(\overrightarrow{u}\left({\overrightarrow{r}}^{\prime },t\right)-\overrightarrow{u}\left(\overrightarrow{r},t\right),{\overrightarrow{r}}^{\prime }-\overrightarrow{r}\right)d{\overrightarrow{r}}^{\prime }+\overrightarrow{b}\left(\overrightarrow{r},t\right)} $$
(23)

where ρ is the density, \( \overrightarrow{f} \) is the internal force term and \( \overrightarrow{b}\left(\overrightarrow{r},t\right) \) is a volume force, such as gravity. The internal force, \( \overrightarrow{f} \), depends on the so-called force state on both the \( {\overrightarrow{r}}^{\prime } \) and \( \overrightarrow{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 \( \overrightarrow{\xi}+\overrightarrow{\eta} \) (Hu et al. 2012). In such a case, the intensity of \( \overrightarrow{f}\left(\overrightarrow{u}\left({\overrightarrow{r}}^{\prime },t\right)-\overrightarrow{u}\left(\overrightarrow{r},t\right),{\overrightarrow{r}}^{\prime }-\overrightarrow{r}\right) \) can still depend on all bonds of both \( {\overrightarrow{r}}^{\prime } \) and \( \overrightarrow{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, \( \overrightarrow{\xi}+\overrightarrow{\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 \overset{..}{\overrightarrow{u}}\left(\overrightarrow{r},t\right)={\displaystyle \underset{H_r}{\int }f\left(\left\Vert \overrightarrow{\eta}+\overrightarrow{\xi}\right\Vert, \left\Vert \overrightarrow{\xi}\right\Vert \right)\frac{\overrightarrow{\eta}+\overrightarrow{\xi}}{\left\Vert \overrightarrow{\eta}+\overrightarrow{\xi}\right\Vert }}d\overrightarrow{r}+\overrightarrow{b}\left(\overrightarrow{r},t\right) $$
(24)

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

$$ f\left(\left\Vert \overrightarrow{\eta}+\overrightarrow{\xi}\right\Vert, \left\Vert \overrightarrow{\xi}\right\Vert \right)=cs $$
(25)
$$ s=\frac{\left\Vert \overrightarrow{\eta}+\overrightarrow{\xi}\right\Vert -\left\Vert \overrightarrow{\xi}\right\Vert }{\left\Vert \overrightarrow{\xi}\right\Vert }, $$
(26)

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

Fig. 9
figure 9

Bond force as a function of bond relative elongation, definition of critical elongation

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 microscale parameters c and s0 with the help of analytical expressions (Ha and Bobaru 2011) as follows

$$ E=c\pi {\updelta}^3\left(1-v\right)/6 $$
(27)
$$ G=9E\updelta {s}_0^2/4\pi $$
(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, mi, initially on grid node positions \( {\overrightarrow{r}}_i \), linked by springs above closest neighbours. The equation of motion becomes:

$$ {m}_i{\overset{..}{\overrightarrow{u}}}_i(t)={\displaystyle \sum_{{\overrightarrow{r}}_j\in {H}_{{\overrightarrow{r}}_i}}ks\left({\overrightarrow{\eta}}_{ij},{\overrightarrow{\xi}}_{ij}\right)\frac{{\overrightarrow{\eta}}_{ij}+{\overrightarrow{\xi}}_{ij}}{\left\Vert {\overrightarrow{\eta}}_{ij}+{\overrightarrow{\xi}}_{ij}\right\Vert }}+{\overrightarrow{b}}_i(t) $$
(29)

where \( {\overrightarrow{u}}_i(t)=\overrightarrow{u}\left({\overrightarrow{r}}_i,t\right) \), \( {\overrightarrow{\xi}}_{ij}={\overrightarrow{r}}_j-{\overrightarrow{r}}_i \) is the bond between points initially located on the \( {\overrightarrow{r}}_i \) and \( {\overrightarrow{r}}_j \) positions, \( {\overrightarrow{\eta}}_{ij}={\overrightarrow{u}}_j-{\overrightarrow{u}}_i \) is the deformation of \( {\overrightarrow{\xi}}_{ij} \), k = c(Δx)4 is a stiffness coefficient, Δx being grid step, and \( {\overrightarrow{b}}_i(t) \) is the body force on the material point initially located on \( {\overrightarrow{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 s0 depend upon the phase (i.e. cell wall or middle lamella) of points \( {\overrightarrow{r}}_i \) and \( {\overrightarrow{r}}_j \). Concretely, if \( \varphi \left({\overrightarrow{r}}_i\right) \) is a phase index, k becomes kij = k(φi, φj) in Eq. (29), where \( {\varphi}_i=\varphi \left({\overrightarrow{r}}_i\right) \) and \( {\varphi}_j=\varphi \left({\overrightarrow{r}}_j\right) \). Critical elongation s0 becomes soij = s0(φi, φj) in the same way. As a consequence, a stiffness kαβ and a critical elongation s0αβ can be attributed to each {α, β} phase-phase couple. Bulk phase elastic moduli Eα can be deduced from kαα values following Eq. (27), bulk phase fracture energies Gα can be deduced from s0αα following Eq. (28) and interface fracture energies Gαβ can be deduced from s0αβ when a ≠ β.

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 ε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 λs = 1 W m− 1 K− 1 and λair = 0.023 W m− 1 K− 1, 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 fs = 0.004 and fair = 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 placed in series or in parallel, respectively). The respective solutions read as

$$ {\varGamma}_{\mathrm{s}\mathrm{eries}}=\frac{1}{\frac{1-{\varepsilon}_{\mathrm{s}}}{\varGamma_{\mathrm{air}}}+\frac{\varepsilon_{\mathrm{s}}}{\varGamma_{\mathrm{s}}}} $$
(30)
$$ {\varGamma}_{\mathrm{parallel}}=\left(1-{\varepsilon}_{\mathrm{s}}\right){\varGamma}_{\mathrm{air}}+{\varepsilon}_{\mathrm{s}}{\varGamma}_{\mathrm{s}} $$
(31)

where Γ represents λ 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).

Fig. 10
figure 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

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 the thermal flux mainly in the radial direction, which explains the inverted anisotropic ratio.

Fig. 11
figure 11

Dimensionless equivalent conductivity, λeq, versus solid fraction, εs, for radial (R) and tangential (T) directions of spruce early- and latewood. Series (dashed line) and parallel (solid line) theoretical models are represented

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).

Fig. 12
figure 12

Dimensionless equivalent diffusivity f versus solid fraction, ε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

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 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−3). 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.

Table 1 Summary of the stiffness computed by MPM modelling for the three tissues selected in the present work

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.

Fig. 13
figure 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 %

Fig. 14
figure 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 %

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.

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).

Fig. 15
figure 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

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, δ = 4Δx, where Δx is the grid step, which is a good compromise 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: Kcw for the cell wall and Kml for lamella. Here, we fixed the ratio between tenacities as Kcw/Kml = 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.

Fig. 16
figure 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

Fig. 17
figure 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

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.

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.

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

    Article  CAS  Google Scholar 

  • 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

    Google Scholar 

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

    Article  Google Scholar 

  • 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

    Article  CAS  Google Scholar 

  • 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

    CAS  Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Article  Google Scholar 

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

    Article  Google Scholar 

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

    Article  Google Scholar 

  • 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

    Google Scholar 

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

    Article  CAS  Google Scholar 

  • 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

    Article  Google Scholar 

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

    CAS  Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Article  CAS  Google Scholar 

  • 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

    Google Scholar 

  • 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

    Article  CAS  Google Scholar 

  • 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

    Article  CAS  Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Article  CAS  Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Article  CAS  Google Scholar 

  • 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

    Article  CAS  Google Scholar 

  • 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

    Article  Google Scholar 

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

    Article  Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Article  Google Scholar 

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

    Article  Google Scholar 

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

    Article  CAS  Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Article  Google Scholar 

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

    Article  Google Scholar 

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

    Book  Google Scholar 

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

    Article  Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Google Scholar 

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

    Google Scholar 

  • 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

    Google Scholar 

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

    Google Scholar 

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

    Book  Google Scholar 

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

    CAS  Google Scholar 

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

    Article  CAS  Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Google Scholar 

  • 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

    Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Article  Google Scholar 

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

    Article  CAS  Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Article  CAS  Google Scholar 

  • 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

    Article  Google Scholar 

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

    Book  Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Article  Google Scholar 

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

    Article  Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Article  Google Scholar 

  • 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

    Article  CAS  Google Scholar 

  • 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

    Article  CAS  Google Scholar 

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

    Google Scholar 

  • 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

    Article  Google Scholar 

Download references

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.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Patrick Perré.

Additional information

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

Rights and permissions

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.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Perré, P., Almeida, G., Ayouz, M. et al. New modelling approaches to predict wood properties from its cellular structure: image-based representation and meshless methods. Annals of Forest Science 73, 147–162 (2016). https://doi.org/10.1007/s13595-015-0519-0

Download citation

  • Received:

  • Accepted:

  • Published:

  • Issue Date:

  • DOI: https://doi.org/10.1007/s13595-015-0519-0

Keywords