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

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

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).
Finally, for each material point of the solid phase, the local anisotropy directions are generated from the two nearest cell contours in such a way as to mimic how the living cell generated the secondary cell wall (Fig. 3). In the final grid, each point carries its own information, such as material type, local anisotropic directions, density, etc.
An important feature was added to the MeshPore software for the specific needs of the present paper. In the case of fragmentation, it is important to distinguish the middle lamella from the secondary cell wall layer. In this case, the position of the middle lamella was iteratively generated from the lumen contour using iterative increments along a direction perpendicular to the contour while accounting for a minimum gap between two neighbouring internal contours. To correctly represent the middle lamella at the intersection of three cells, the contour evolution locally stops when a maximum curvature is attained. Once all contours are stabilised, two "solid" types (secondary wall and middle lamella) are allocated to all points located inside the cell wall (Fig. 4). Figure 5 depicts the three tissue morphologies tested in the present work. The most refined meshes were generated for spruce (earlywood and latewood), in which the middle lamella was distinguished from the secondary cell wall. This feature was required to model fracture. The third case is a cellular arrangement of poplar, a ring-porous hardwood species. This example is interesting due to the double porosity: fibre and longitudinal parenchyma lumen and vessel lumen. Using a medium refinement at the cell wall level, this dual-scale porosity can be represented using a number of material points comparable to those needed in the case of spruce, where the refinement is used to separate the middle lamella and the secondary cell wall. 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; and, therefore, has interesting properties, such as simplicity of development, suitability for parallel computing and flexibility in geometrical shape .
In the standard LBM, the single-particle distribution function obeys the following equation: where c ! i is the discrete particle velocity vector and the righthand side term represents the collision term, the so-called single relaxation time. The set of vectors 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 eq i x ! ; t À Á is the equilibrium distribution function at site x ! and time t. The above equation is well known as the Bhatnagar-Gross-Krook (BGK) approximation for the collision operator Ω i x ! ; t À Á . This means that after a collision, the distribution function tends to an equilibrium Maxwell-Boltzmann distribution.
In the present work, a standard D2Q9 model was adopted, i.e. a 2D computational domain (x, y) involving nine discrete velocities, c i ∀i={0,1,…8}. These velocities allow the particles to jump to neighbouring nodes in one lattice time step

Contour 1
Contour 2 Point of the grid   6). They are defined as follows: where the lattice speed, c, is related to the lattice spacing c ¼ δx δ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 x-direction, 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 Chapman-Enskog  expansion for the LBM equation (Eq. 1 including Eq. 2). The relationship between macro-and meso-scales 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 wheref k andf 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 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 sub-sets: 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 steady-state regime. The equivalent diffusivity is then computed as the average flux divided by the macroscopic gradient imposed by the Dirichlet boundary conditions.

Material Point Method
The Material Point Method (MPM) is derived from the Particle in Cell (PIC) method originally developed in the 1960s by Harlow in computational fluid mechanics. The MPM has been successfully applied to computational solid mechanics problems (Sulsky et al. 1994;Sulsky and Schreyer 2004;Bardenhagen and Kober 2004;Guilkey et al. 2006). One of the most important features of the MPM is that it allows a straightforward discretisation of complex material shapes, including direct discretisation from 2D or 3D images, as well as an efficient and robust handling of contacts between material surfaces.
In the MPM, the domain of interest is defined by a collection of material points p, p=1… N p . As a direct consequence of this principle, each material point carries its properties, such as position, velocity, acceleration, strain and stress, in a Lagrangian way (step 1 of Fig. 7). The particle mass, m p , is computed from the medium density and the initial point distribution. In the classical version of MPM, no volume is affected to material points: therefore, the resulting density field is approximated by a non-continuous field defined by Dirac functions (capital letters are used for Lagrangian values): The m p values are constant over time, which ensures mass conservation. The weak form of the momentum equation is solved with the aid of a background grid. This grid consists of simple finite elements, usually isoparametric, four-node, quadrilateral elements in two dimensions. This simple grid only has to include all the material points and is reset at each time increment, thereby avoiding the classical FE problems that result from large deformations. The algorithm allows all quantities at material points to be updated stepwise. To do so, the first step is to map the material points on nodes i, i=1… N g of the FE mesh (step 2 of Fig. 7). Superscript k refers to time step k and superscript L refers to the end of the Lagrangian step. Classical FE shape functions, N i , are used for this purpose: 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), σ 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 least-squares 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 Fig. 7 Principle of the Material Point Method (adapted from Sulsky et al. 1994) 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.

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 where r ! is a point position at t=0, r ! 0 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 short-range forces are involved. A bond is defined as ξ is the deformation of the bond ξ ! (Fig. 8). The force balance leads to the following integrodifferential equation: where ρ is the density, f ! is the internal force term and b ! r ! ; t À Á is a volume force, such as gravity. The internal force, f ! , depends on the so-called force state on both the r ! 0 and 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 ξ ! þ η ! (Hu et al. 2012). In such a case, the intensity of f !
still depend on all bonds of both r ! 0 and 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, ξ ! þ η ! , 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 non-local approach. Both the elastic modulus, E, and the surface energy, G, can be deduced from the Fig. 8 Schematic view of the peridynamic approach: initial state (t=0) and deformed state (t>0). Definition of horizon δ, neighbourhood H r of r ! , bond ξ ! , displacement u ! r ! ; t À Á and bond deformation η ! 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 r ! i , linked by springs above closest neighbours. The equation of motion becomes: i is the bond between points initially located on the r ! i and r is a stiffness coefficient, Δx being grid step, and b ! i t ð Þ is the body force on the material point initially located on r ! i . Clearly, such an approach is similar to molecular dynamics (Seleson et al. 2009), and Eq. (29) can be integrated with the help of the well-known velocity-Verlet integration algorithm (Allen and Tildesley 1986).
To take heterogeneous mechanical properties into account, both k and s 0 depend upon the phase (i.e. cell wall or middle lamella) of points r 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 {α, β} 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 s 0αα following Eq. (28) and interface fracture energies G αβ can be deduced from s 0αβ when a≠β.

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 =1W m −1 K −1 and λ air = 0.023W 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 Fig. 9 Bond force as a function of bond relative elongation, definition of critical elongation Fig. 10 The local field obtained at steady state in the case of thermal conductivity. Temperature values are in colour and heat flux is shown as arrows 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 Fig. 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 Fig. 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 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).

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. 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.  Fig. 13 Stress-strain curves obtained for spruce earlywood in the radial and tangential directions. Computational images of the cellular deformation for the radial test were selected at strain levels of 5, 12, 25 and 48 %

Fragmentation
Peridynamic studies of quasi-static crack growth are quite rare (Kilic and Madenci 2010). To perform quasi-static numerical tests, an equilibrium search approach is required. To perform damping of mechanical energy in the material and reach an equilibrium state, we applied a viscous dissipation as a body force to each material point. After each strain increment, critical bonds are broken and a new equilibrium search is performed, if needed, until a stable mechanical equilibrium is reached (Fig. 15). This algorithm was implemented in parallel in a custom code written in FORTRAN 95 and using Message Passing Interface (MPI) routines. A typical computing time is 8 h on a 16-core server.
In the present application, we used a horizon, δ=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 ¼ ffiffiffiffiffiffiffi EG p , is defined for each phase: K cw for the cell wall and K ml for lamella. Here, we fixed the ratio between tenacities as K cw /K ml = 4. A pre-notch was performed on one side of the spruce latewood sample and the tensile fracture test was performed. The same test was performed for the earlywood and latewood of spruce.
As we can see in Fig. 16, a crack grew from the pre-notch and gradually extended to the other side of the sample. The crack remained inside the middle lamella, and cell walls were not damaged during the failure test, as expected according to the tenacity ratio. The same behaviour was observed in the case of spruce latewood (Fig. 17). In both cases, the crack remained within the most fragile phase, the middle lamella.
Such a test clearly shows that the peridynamic approach is a convenient framework to simulate crack growth in wood, despite the complex structure and mechanical properties inside this material. Fig. 16 Snapshots of successive steps of a peridynamic simulation of spruce earlywood failure under tension (colour online). The crack is in black. Colours are related to tensile stress: green when stress is close to 0, blue below 0 and red above 0

Stress (MPa)
Fig. 14 Stress-strain curves obtained for poplar in the radial and tangential directions. Images of the cellular deformation for the radial test were selected at strain levels of 10, 20 and 42 % New modelling approaches to predict wood properties

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.