50 years Annals of Forest Science
 Original Paper
 Open Access
 Published:
New modelling approaches to predict wood properties from its cellular structure: imagebased representation and meshless methods
Annals of Forest Science volume 73, pages 147–162 (2016)
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 socalled 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 nonlinear 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 ultrastructural 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 Xray 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 structureproperty relationships in the transverse plane of wood (radialtangential) are more complicated. In this plane, solid components act both in parallel and in series at different spatial scales (multilayered 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 (SanchezPalencia 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):

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.
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 coexist. 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; SedighiGilani 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 multilayered 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 Xray μ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 imagebased 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 imagebased 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 dualporosity 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 sublayers of the cell wall to be represented in singleporosity tissues). Each method is then briefly presented, and the last section presents results and discussion.
2 Imagebased 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 imagebased 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).
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 ringporous 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 dualscale 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.
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 socalled 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 singleparticle distribution function obeys the following equation:
where \( {\overrightarrow{c}}_i \) is the discrete particle velocity vector and the righthand side term represents the collision term, the socalled single relaxation time. The set of vectors \( {\overrightarrow{c}}_i \) are the possible velocities a particle can move from a lattice site to a nearestneighbouring site at each time step t + δt. For simplicity, we used a “single relaxation time” form for the collision term:
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 BhatnagarGrossKrook (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 MaxwellBoltzmann 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} ∀ i = {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:
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.
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:
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):
The standard 2D diffusion equation for the normalised mass or temperature field,
was solved in a 2D rectangular porous medium, [x_{1},x_{n}]×[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: Ψ(x, y, 0) = 0. Along the xdirection, Dirichlet boundary conditions were applied: Ψ(x_{1}, y, t) = 1 and Ψ(x_{n}, y, t) = 0. Adiabatic conditions were set along the ydirection: φ(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 ChapmanEnskog expansion for the LBM equation (Eq. 1 including Eq. 2). The relationship between macro and mesoscales gives
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:
The normalised flux has the following form:
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):
where Ψ_{1} and Ψ_{n} are fictive variables obtained from balance equations to ensure the desired Dirichlet value. For example, for the plane (x = x_{1}), we use
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
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 [x_{1}, x_{n}] × [y_{1}, y_{m}] is, therefore, divided into two subsets: the cell lumens, Σ_{1}, and the solid phase, Σ_{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 steadystate 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 noncontinuous field defined by Dirac functions (capital letters are used for Lagrangian values):
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, fournode, 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:
Internal forces are mapped on the grid using boundary conditions and the stress level at each point p (step 3 of Fig. 7):
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:
At the grid nodes, the weak form of the momentum equation reduces to
The velocities at the grid nodes are computed from the point velocities using a leastsquares approach with point masses as weighting factors:
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):
From the nodal information, the velocities and positions are updated back to material points using the element shape functions (step 5 of Fig. 7):
The strain increment of point p is computed from the element shape functions:
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).
Such a paradigm is very different from classical continuous matter mechanics, in which only shortrange 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 integrodifferential equation:
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 socalled force state on both the \( {\overrightarrow{r}}^{\prime } \) and \( \overrightarrow{r} \) positions (Silling 2000; Silling et al. 2007). Among all possible statebased peridynamic models, in ordinary statebased 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 bondbased 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 fibrereinforced 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:
Here, we consider only linear models as a function of the bond relative elongation s, as follows:
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 nonlocal approach. Both the elastic modulus, E, and the surface energy, G, can be deduced from the microscale parameters c and s_{0} with the help of analytical expressions (Ha and Bobaru 2011) as follows
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 \( {\overrightarrow{r}}_i \), linked by springs above closest neighbours. The equation of motion becomes:
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 wellknown velocityVerlet integration algorithm (Allen and Tildesley 1986).
To take heterogeneous mechanical properties into account, both k and s_{0} 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 k_{ij} = 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 s_{0} becomes s_{oij} = s_{0}(φ_{i}, φ_{j}) in the same way. As a consequence, a stiffness k_{αβ} and a critical elongation s_{0αβ} can be attributed to each {α, β} phasephase couple. Bulk phase elastic moduli E_{α} can be deduced from k_{αα} values following Eq. (27), bulk phase fracture energies G_{α} can be deduced from s_{0αα} following Eq. (28) and interface fracture energies G_{αβ} can be deduced from s_{0αβ} 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 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 placed in series or in parallel, respectively). The respective solutions read as
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).
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.
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 stressstrain 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.
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.
4.3 Fragmentation
Peridynamic studies of quasistatic crack growth are quite rare (Kilic and Madenci 2010). To perform quasistatic 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 16core 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: 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 prenotch 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 prenotch 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.
5 Conclusion
This paper presents a comprehensive strategy to predict different wood properties from anatomical images. The starting point is an imagebased 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 postprocessing.
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/i2012121177
Agwai A, Guven I, Madenci E (2011) Crack propagation in multilayer thinfilm 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 woodwater relationship due to heat treatment assessed on microsamples 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/S0718221X2014005000015
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/17426596/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/09650393/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 Xray 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 largescale digital imagebased finite element analyses. J Biomech 33:255–259. doi:10.1016/S00219290(99)001414
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/s1070400993463
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 3D 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
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 ChapmanEnskog 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/S00457949(98)003319
Hu W, Ha YD, Bobaru F (2012) Peridynamic model for dynamic fracture in unidirectional fiberreinforced 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 quasistatic 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/S00457949(02)004066
Lallemand P, Luo LS (2003) Lattice Boltzmann method for moving boundaries. J Comput Phys 184:406–421. doi:10.1016/S00219991(02)000220
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 hygroelastic properties of wood fibres. J Mater Sci 42:10254–10274. doi:10.1007/s1085300611999
Perré P (2001) Wood as a multiscale 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 imagebased 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 woodwater relations on microsamples. 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/S13858947(01)002704
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 bioproducts 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 microsamples in the three material directions. Wood Sci Tech 47:481–498. doi:10.1007/s002260120511x
Ruiz G, Oritz M, Pandolfi A (2000) Threedimensional simulation of the dynamic Brazilian tests on concrete cylinders. Int J Num Methods Eng 48:963–994
Salmen L (2004) Micromechanical understanding of the cellwall structure. C R Biolog 327:873–880. doi:10.1016/j.crvi.2004.03.010
SanchezPalencia E (1980) Nonhomogeneous 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
SedighiGilani M, Navi P (2007) Experimental observations and micromechanical modeling of successivedamaging phenomenon in wood cells’ tensile behavior. Wood Sci Tech 41:69–85. doi:10.1007/s0022600600945
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 longrange forces. J Mech Phys Solids 48:175–209. doi:10.1016/S00225096(99)000290
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/s1065900791251
Succi S (2001) The lattice Boltzmann equation for fluid dynamics and beyond. Clarendon Press Oxford, New York, p 308. doi:10.1063/1.1537916
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 historydependent materials. Comp Methods Appl Mech Eng 118:179–196. doi:10.1016/00457825(94)901120
Suquet PM (1985) Element of homogenization for inelastic solid mechanics in homogenization techniques for composite media. Lecture Notes in Physics 272 SanchezPalencia and Zaoui Ed SpringerVerlag
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/i2007102019
Ulrich D, van Rietbergen B, Weinans H, Rüegsegger P (1998) Finite element analysis of trabecular bone structure: a comparison of imagebased 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 multiaxial dynamic loading. J Mech Phys Solids 49:2768. doi:10.1016/S00225096(00)000284
Acknowledgments
The authors would like to thank the French Agence Nationale de la Recherche (ANR, project ANR 10 – HABISOL00503) for its participation in founding this research work.
Author information
Authors and Affiliations
Corresponding author
Additional information
Handling Editor: Erwin Dreyer
Contribution of the coauthors
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.
About this article
Cite this article
Perré, P., Almeida, G., Ayouz, M. et al. New modelling approaches to predict wood properties from its cellular structure: imagebased representation and meshless methods. Annals of Forest Science 73, 147–162 (2016). https://doi.org/10.1007/s1359501505190
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s1359501505190
Keywords
 Lattice Boltzmann Method
 Material Point Method
 Peridynamic
 Real anatomical structure
 Mechanical properties
 Stiffness
 Fracture
 Diffusion
 Thermal conduction
 Collapse