Developing the orthotropic linear-elastic model for wood applications using the FE method

Abstract

The purpose of this work is to develop the three-dimensional (3D) finite element (FE) modeling approach for the linear mechanical behavior of the wood material. The developed framework consists of implementing the 3D constitutive equations using the linear elasticity theory. Wood is a complex, porous, fibrous, inhomogeneous, highly anisotropic material. Various wood materials are considered, such as poplar, spruce, and maple specimens to validate the applicability of the FE modeling. The framework is implemented in explicit code, written in Fortran language, and based on the PETSc library. To that purpose, using up-scaling methods, such as the homogenization technique that allows the prediction of the macroscopic property of materials, in our case the effective Young's modulus, is investigated. A numerical approach to control the equivalent micromechanical properties is presented using a representative elementary volume (REV) concept. Here, we investigate the convergence trend of material properties and structural geometries according to REV size. This framework aids considerably in predicting the mechanical properties of a given material microstructure. In order to determine the mechanical properties of the wood material in its anisotropic direction, the traction and compression loadings were performed using numerical tests. The FE modeling of some cases is presented for the final validation.

---

1 Introduction

Wood is a natural composite material with rich and excellent mechanical properties such as its low density.^1,2 Some wood is hard, containing fibrous tissue due to physical hardness or density. These fibers are robust in tension and emmersed in matrix lignin to resist compression. Wood is also a heterogeneous, anisotropic, and variable material with a complex structure which covers various length scales. This anisotropic property^3,4 explains the dependency of the mechanical properties on its material directions [radial (R), tangential (T), longitudinal (L)] and differs among species. This heterogeneity means that huge stiffness changeability is encountered in wood. In order to express this complex feature, multiscale frameworks have been developed on finite element simulations with computational homogenization aspects. This is used for generation in numerous applications,^5–7 such as building, construction, transportation, and landscaping timbers. In order to completely utilize the potentiality of the wood material, ameliorating knowledge of its mechanical behavior and constitutive frameworks is essential.

Over the last few years, wood has emerged as a promising material due to its significant and practical applications.^8,9 It has captivated many researchers due to its superior mechanical properties. They have provided much effort to investigate its large variability in these macroscopic material properties. The objective is to avert the requirement for essential safety factors derived from the highly unprofitable cost of timber members in these applications. Providing comprehensive knowledge about microstructural properties of individual wood specimens is required to correspond to them at the macroscale. The properties of these equivalent homogeneous materials are determined by employing homogenization techniques.^10,11 The link between these scales prepares accurate and compatible sets of macroscopic characteristics of wood in orthotropic directions.

There are two main sorts of wood cell, earlywood and latewood, also specified as tracheids, and they are devoted considerably to the volume of softwoods.^12,13 When favourable climatic conditions are met, there is a high growth rate of the earlywood cells that are laid down in the spring season.^14,15 However, the latewood cells are formed in the late summer, when the rate of growth decreases.^16,17 They provide mechanical strength with different characteristics. The cell walls of the earlywood cells are thin with a large open pore space, while those of the latewood are thick with a smaller pore space. Büyüksar et al.¹⁸ investigated the determination of the mechanical properties of earlywood and latewood sections in standard size samples with comparing the computed and measured mechanical properties.

Our primary motivation for developing the framework is to exploit 3D imaging techniques, including X-ray computer tomography (CT) and magnetic resonance imaging (MRI). These techniques offer many conveniences, resulting from significant improvements in high-resolution and reconstruction strategies. One technology that has benefited from such improvements is X-ray tomography. This technology is based on using different radiographs from a given object, viewed at various angles, to construct a bidimensional slice of that object. This paper will consider the conventional μCT imaging^19–21 to highlight many aspects of wood anatomy, in particular, a comprehensive characterization of the wood formation and its reconstruction of the geometric structure.

When determining the homogenization approach established on the finite elements method (FEM), the mechanical behavior of wood materials is usually expressed by using representative elementary volume (REV).^22,23 This modeling is aimed at analyzing the behavior of heterogeneous wood structures as it allows making the problem simple by examining a sub-volume of a total structure, without influencing the precision of the results.^24,25 In particular, the homogenized field variables of a material wood on the macroscopic scale depend on its size, making the effect of REV size questionable. The definition of REV has remained difficult in the literature for properties related to linear elastic models. The homogenization approach highlights the mechanical behavior obtained here since the micromechanical response at each sample point relies on the macroscopic deformation. In the beginning, the concept of REV analyses was presented by Hill (1963). According to a continuum mechanics viewpoint, the REV must contain an adequately large number of heterogeneities. Therefore, this work aims to investigate the convergence of the macroscopic elastic properties^26,27 by ensuring an appropriate choice of the REV size. This volume REV is subjected to specific meshing criteria to satisfy the convergence behavior. Some researchers have investigated various types of REV; for example, Gonzalez and Llorca,²⁸ studied a square REV, which is composed of a homogeneous and random dispersion of fibers included in the polymeric matrix. Previous work²⁹ was devoted to determining the adequate characteristics of a REV to perfectly model the mechanical behavior of a unidirectional glass fiber-reinforced polymer composite.

The present contribution is aimed at providing a simplified tool to determine the properties of the deformed microstructure with high computational efficiency based on the PETSc library,³⁰ widely used in various other packages.^31,32 This library is considered popular and suitable for determining solutions of partial differential equations (PDEs). It provides conveniently vast flexibility for users, that includes parallel linear and nonlinear solvers.^33,34 Many previous works^35–37 on deformation microstructure have been carried out in which the microstructure is influenced by phase compositions and processing parameters. In the present improvement, the global matrix associated with this discretization is stored in the sparse format supported by PETSc. This implementation focuses on the compressed sparse row (CSR) storage,^38–41 which defies scientific computing due to its efficiency allowing robust storage. The computational cost for this storage grows linearly depending on the system size, with some benefits. These benefits consist of reducing the memory overhead and avoiding storing zero elements. This benchmark problem illustrates the accuracy^42–44 of the proposed strategy to linear elasticity.

The rest of the paper is structured as follows. Section 2 presents a simple way to obtain wood morphology from tomographic images. Then, image-analysis techniques are used to give the morphology distributions and extract morphological details. This section describes briefly the implemented FEM framework for orthotropic elastic problems. Section 3 reports the numerical results of testing a wood structure with highlighting 3D mechanical analysis. To illustrate the newly developed homogenization method, numerical results for a one-two-phase wood species displaying the equivalent elastic behavior is presented. Finally, Section 4 concludes this work.

2 Materials and methods

This section is devoted to describing the main workflow for how to provide the binary representation of the wood morphology. The choice of studying these wood sample species is motivated by the porosity, anisotropy, and heterogeneity of the material characterized by some structural anomalies. All the specimens' preparation and conditioning will be detailed. For modeling the structures in anisotropic materials, the tool most often used is the general purpose FE tool. This section aims to describe mathematically this tool which is specifically designed for anisotropic materials to simulate small deformations.

2.1 Materials

In this study, the selection of the wood material was of the upmost importance, as the goal was to prepare samples containing only natural wood. The material wood samples involved in this paper were poplar (Populus euroamericana Koster), spruce (Picea abies) earlywood, spruce latewood, and maple specimens (Acer rubrum L.). The poplar specimen is from a forest located in Auménancourt-le-Petit (France). It was provided by the Huberlant sawmill in Cormicy (France), which cut the tree into boards and finally air-dried the wood. The spruce specimen is derived from a plantation in the Le Châtaignier forest in Riotord (France). These samples were provided from a single defect-free board with a horizontal grain angle. Fig. 1 depicts the scanned wood in three dimensional-space.

Fig. 1 The 3D scanned samples of poplar (left), spruce (middle), and maple (right).

2.2 Tomography and segmentation

The scanning lab tomography is used to obtain a stack of images to capture the structure of each morphology. X-ray nano-tomography is exploited as one tool available to wood anatomists to study the three-dimensional organization with a high-resolution image result. In order to carry out the segmentation of the 3D representation from these stack images by the reconstruction software, further processing is required using various software. The image processing is conducted by general software such as Avizo, Matlab, Python, etc. Fiji⁴⁵ is one of the practical software packages used in the current work to identify objects and study their properties based on the acquired images. Fiji is an extension of the ImageJ⁴⁶ and ICY⁴⁷ software. It is a software package that offers a wide range of plugins for enhancing image visibility. The purpose of this study is to provide quick image treatment by selecting the thresholding tool. The automatic thresholding method of Otsu⁴⁸ proposed by Fiji is one possible approach to get a binary image. The Fiji process menu offers a wide range of easy-to-use binary tools, such as a filter to remove noise or control a threshold value.^49,50 The aim is to identify and separate the present two phases (gas and solid) and generate the morphology's tetrahedral mesh.⁵¹

The wood samples are acquired by X-ray nano-tomography, EasyTomXL Ultra 150–160, RX Solutions. 2D-projection images are realized by measuring the X-ray attenuation originated from a selection. The system is equipped with a nano-focus source with different detector types for both fast and in situ (flat panel detector) and high-resolution scanning (CCD camera). The scanning resolution is on the order of 1 μm and takes at least four hours. Fig. 1 depicts the full scan of the wood samples, which have a diameter of about 1 mm, and the radius and height of this cylinder are represented by 2.82 and 3.11 millimeters, respectively. The image file sizes are extremely large, which forces an adequate ROI (Region of Interest) to be chosen within the sample. This sub-volume is a parallelepiped whose center and dimensions should be selected. Each parallelepiped sample is treated with ImageJ-Python software to make 3D volumic reconstruction possible. The post-processing of 3D reconstructed wood specimens is presented in Fig. 2, which is shared into six diagrams. The top three diagrams Fig. 2b–d show that each sample is a randomized phase two study. The first is associated with the gaseous phase, and the second is for the solid phase (matrix). Since the free water occupied by the gaseous phase does not interact with the solid one, the porosities do not lead to mechanical variations. Then, the displacement field basis of all the solid matrix nodes is considered to neglect the gas distribution in the pore space. However, diagrams Fig. 2f–h depict the extracted solid phase from the same 3D reconstructed wood specimens. For this reason, only the mesh of the solid phase has been implemented in the finite element software.

Fig. 2 The extracted solid phase from the 3D reconstructed wood specimens scanned by X-ray nano-tomography presented in Fig. 1 to model the material microstructure. Each specimen having the same dimensions of 128 × 128 × 128 voxels is treated with ImageJ-Python software.

2.3 Elastic properties for the model material

We assume that the mechanical properties used in the model are considered constant and independent of the porosity. This hypothesis is justified in simplifying the modeling work. The aim is to find better parameters existing in the literature for which the model is convergent and affects the numerical results. The ultrastructural parameters reported for wood materials are a collection of different values of engineering elastic constants. These values are computed from molecular models^52,53 or predicted from the behavior of similar materials.^54,55 Neagu and Gamstedt⁵⁶ focus mainly on the knowledge of the wood fiber structure by providing these hygroelastic properties from some wood samples. These properties are important information used essentially in the FE calculations.⁵⁷ The mechanical parameters to identify the solid phase for these species, poplar, spruce, and maple samples, are described as follows. Note that the hygroelastic properties of the cell wall layers computed with the recursive algorithm of Sutcu⁵⁸ must be given as input to the mechanical model presented in terms of physical units Gigapascal (GPa). The longitudinal Young's moduli E_L of wood specimens is equal to 64.0 GPa, while the cell wall and lumen are obtained from the radial and tangential Young's modulus E_R and E_T using the same value of 6.08 GPa, i.e., E_L = 64.0 GPa, E_R = E_T = 6.08 GPa. Due to the transversely isotropic elastic properties, the Poisson ratios ν_TL and ν_RL, and the shear moduli G_TL and G_RL, are the same, i.e., ν_TL = ν_RL = 0.03, G_TL = G_RL = 2.54. The Poisson ratio ν_RT and the shear moduli G_RT, for transverse strain in the tangential direction (T) when stress is applied in the radial direction (R), are 0.42 GPa and 2.14 GPa, respectively, formally, ν_RT = 0.42 GPa, G_RT = 2.14 GPa.

2.4 Up-scaling procedure

One of the important challenges in the FEM framework is enforcing the boundary conditions. The penalty method has been proposed^59,60 to control this problem. This method offers an approximation for boundary conditions. The benefit of this method is that there are no further global constraint equations and only one unknown, which is the displacement field.^61,62 Note that the penalty parameter involved in the implementation has to be larger, exceeding the other elements of the stiffness matrix.⁶³ Hadjicharalambous et al.⁶⁴ showed the high accuracy and good convergence from the present penalty method with some concrete numerical examples. They tested this convergence by increasing the degrees of freedom and analyzing the global system matrix. This choice leads to good convergence according to a refined mesh. Consider a 3D material microstructure in the form of a cubic shape, subjected to uniform displacement on its one particular face following the orthotropic direction as shown in each diagram of Fig. 3. In the following, the penalty method is applied and validated under three computational homogenization scenarios in wood materials. Each scenario means that a loading is applied taking into account the boundary conditions for uniaxial loading cases. When the displacement is applied at one face, for instance the upper face normal to the longitudinal direction (L), as shown in the left diagram, all the lateral and lower faces stay blocked to capture much-averaged values. This means that the penalty method is operated to apply a null displacement at these faces, and also the targeted displacement. To validate the purpose of the homogenization approach, this strategy is continued with imposing the associate displacements in the transverse plane (R,T), as illustrated in the middle and right diagrams.

Fig. 3 A 3D cubic material microstructure subjected to a uniform displacement on each particular face.

Let Ω be the domain of a body and ∂Ω its boundary disjointly defined by the surfaces Ω^u and Ω^t, which are related to Dirichlet's and Neumann's boundary conditions. This means that the boundary ∂Ω is subdivided into two parts Ω^u and Ω^t, i.e., ∂Ω = Ω^u ∪ Ω^t, Ω^u ∩ Ω^t = ∅. The mixed fundamental boundary value problem of the linear elastic model aims to determine the distribution of stress σ, strain ε and displacement u throughout the body, when it has constrained displacements ū defined on Ω^u and is loaded by an external system of distributed surface and body forces with densities denoted by t on Ω^t and f in Ω, respectively. In what follows, we will introduce the governing equations for solving the problem of linear elasticity. These equations are used to solve the elasticity problem on the REV model including the balance equation interfaced with boundary conditions for non-homogeneous materials. The REV occupies the computational domain. Ω is denoted as the REV domain and V is the volume of the REV. In the model, displacements and tensions are supposed to be continuous across interfaces, meaning that all the phases presented in the REV model are superbly bonded. Thus, the constitutive equations modeling the linear elastic problem can be formulated as,

In which σ is the stress tensor and f is the body force vector. The vector t represents the tension load on part of the static boundary Ω^t prescribed by the value ‖t‖ and n is the exterior unit normal on the boundary Ω^t. It is well established that the mechanical properties of wood materials depend on their orientation. Mechanics of wood materials usually describe longitudinal and transverse directions related to their structure, which are described in Fig. 3. The orthotropic wood models consider different material behavior for each orthogonally related direction. Consequently, more parameters are requested to report material behavior fully compared to isotropic materials. Unlike isotropic materials that are parameterized by a unique Young's modulus and Poisson's ratio, orthotropic materials have three various Young's moduli E₁, E₂, and E₃, one for each orthogonal direction, and six Poisson ratios ν_ij, for i ≠ j, only three of which are independent. In order to keep the elasticity tensor matrix C symmetric, the Poisson's ratios have to satisfyAssuming that the physical variables belong to the microscale problem, that is, are related to the geometry and material of the volume REV. C and B denote the matrices of elasticity and partial derivatives of the interpolation functions. The stain ε and stress σ tensors could be expressed in terms of the displacement u as,

ε = Bu. (3)

σ = CBu. (4)

The stiffness matrix C, the fourth-order tensor of orthotropic wood materials, shows elastic constants, which can be presented in the following form: σ = Cε,where the term Δ_c is defined as,C is Hooke's matrix associated with the stiffness form of linear elastic Hooke's law. The stiffness operator is represented by a 6 × 6 symmetric matrix whose elements are expressed in equality (5). In our case of a fully anisotropic elastic material, 9 independent terms are required to fill the matrix C. V⁰_Ω is denoted by a set of continuous and sufficiently regular functions, zero-valued in Ω^u. The first expression of eqn (1) is multiplied by a test function v lying in the space V⁰_Ω of kinematically admissible displacements at zero, and next integrated by parts. In this context, the resolution of the elasticity problem consists of determining the displacement field u corresponding to the solution of the following variational problem,Based on the established micromechanical modeling, the homogeneous finite element simulation of REV models is carried out. The REV models of increased voxel size with different volume fractions are established, and their effective stiffnesses are calculated by using the numerical homogenization calculation technique. By combining the numerical homogenization computation technique, the effective stiffness property was predicted. To investigate the mechanical effects from the homogenized behavior, uniaxial tests are considered. This consists of doing uniaxial tensile tests. Fig. 3 depicts the schematic diagram of building the loading method adopted in the tensile test process. A constant displacement is applied to the sides of the specimen following the direction loading. The macroscopic stress and strain tensors [small sigma, Greek, macron] and [small epsilon, Greek, macron] of the macroscopic homogeneous REV model can be computed by averaging the microscopic variables or local stress and strain σ and ε over the entire volume of REV, and expressed according to micro-macro relation,

In which Ω_s is the solid phase of all the microstructure REV. The homogenized strain [small epsilon, Greek, macron]_Th is computed exactly from the imposed displacement field. This choice is motivated by taking into consideration the repartition of pores, particularly their volumes. The macroscopic properties of a wood material, are partially controlled by the material properties contrasting between the solid and porosity phases. The macroscopic strain yields a homogeneous deformation on the REV, which is imposed by a displacement. The result [small epsilon, Greek, macron]_Th is presented for the average support displacement u_Direchlet over the support length L_s, defined as,

Let C̄ be the homogenized effective stiffness tensor. In reality, C̄ should be considered as a 6 × 6 matrix. Since three computed homogenized material constants of the shear moduli are not involved in our computational homogenization procedure, the effective elasticity matrix C̄ will be of 3 × 3 elements. The numerical modeling effort is performed using three series of computational homogenization analyses. Then, the homogenization method for heterogeneous media is carried out meaningfully to estimate the three effective elastic properties. The effective stress–strain relation for the averaged stress and strain can be written as,

[small sigma, Greek, macron] = C̄[small epsilon, Greek, macron]_Th. (11)

[small sigma, Greek, macron]_ij denotes the volume average of the corresponding stress component in REV. Under radial loading, the first line of the averaged stiffness matrix C̄ can be calculated as,Similarly, from the tangential and longitudinal loadings, we can completely fill the matrix C̄,By inverting the constitutive matrix defined by (13), the compliance matrix S is obtained. Then the effective elastic constants of REV can be calculated according to the matrix S. In particular, the effective Young's moduli E^h₁, E^h₂, and E^h₃ are exactly the inverse of the diagonal elements of S̄,The strictly lower triangular part from matrix C̄ is used to obtain the three effective Poisson's ratios ν^h₁₂, ν^h₁₃, and ν^h₃₂. However, the normalized effective Young's moduli ν^h₂₁, ν^h₃₁, and ν^h₂₃ are deduced from the strictly upper triangular part. Due to the symmetry of matrix C̄, these ratios satisfy the proportionality equation defined by (2). Formally,

ν^h₁₂ = −C₂₁E₁, ν^h₃₁ = −C₁₃E₁, ν^h₃₂ = −C₃₂E₃. (15)

The up-scaling procedure aims to train material wood models for the homogenization responses. The microscale REVs are then subjected to three various loadings, and the responses are recorded and used to validate the macroscopic data-driven model. For completeness, we give a summary of this procedure and the applications to computational anisotropic mechanics problems. Our computational process often involves the following steps:

• First, the computational domain is decomposed into small tetrahedral elements of four points, and each point has three degrees of freedom.

• Then, the local stiffness matrix for each element is of size 12 × 12. Each local matrix is presented by three unidimensional vectors satisfying the optimization procedure to avoid storing the nonzero entries of the matrix. The local element stiffness matrix can be determined in terms of the discrete matrix B and C given in equalities (3)–(5).

• The assembly process consists of taking into account all the elementary contributions. Each local system determines each contribution to build a global linear system.

• The applied displacement on the boundary is considered by updating the coefficient of the global matrix, which is presented in a sparse format. These boundary conditions are incorporated in the FE formulation via the invoked penalty method. Finally, the global linear system is solved with PETSc in parallel using a KSP (Krylov subspace accelerator) solver to determine the unknown displacement field.

• The application of homogenization for stress–strain behavior is demonstrated. The macroscopic elastic properties are determined via micro-to-macro scale transition analysis. In particular, the macroscopic stress [small sigma, Greek, macron] and strain [small epsilon, Greek, macron] tensors are estimated in a discrete assembly viarelations (8) and (9).

The stages mentioned in Fig. 4 summarize the resulting computational homogenization leading to the upscaling of the microscopic scale towards a macroscopic one for specific loading.^65,66 During the loading process restricted to three anisotropic directions, three homogenization schemes are grouped and factorized^67,68via equality (11). The explicit expressions of the effective elastic properties are computed using equalities (14) and (15) based on the matrix C̄.

Fig. 4 The computation process based on the homogenization scheme employing a FE discretization with displacement.

3 Results and discussion

This section aims to demonstrate the capability of a finite deformation approach for analyzing elastic orthotropic materials. For modeling the structures in anisotropic materials, the tool most often used is the general purpose FE tool described in Section 2.4. This standard approach can be used for different categories of materials and their properties. This model is validated for the wood specimens. A realistic loading is applied at the microstructure volume during the deformation process. The deformed REV is exhibited and analyzed in terms of the realistic microstructure. As a simple rule based on the presented curve shape, we used the range of deformation [0%,3%] for all mechanical tests.

3.1 Numerical test

This part is devoted to illustrating the ability of the finite element method (FEM) framework to function. First, this framework predicts the mechanical properties, particularly the nodal displacements in the standard orthotropic directions for the poplar, maple, and spruce wood specimens. This framework is validated against some uniaxial mechanical testing systems responding to the tensile loadings. Fig. 5 depicts the displacement field at each test following the orthotropic direction. The material is assumed to be elastic and homogeneous, and is determined by mechanical properties, including elastic properties described in Section 2.3. The Dirichlet boundary conditions have been defined by applying a horizontal displacement on the top lateral face following the radial direction (R), and by blocking the bottom lateral face, as shown in Fig. 5a. There is a null displacement on the four rest faces. The strain ε₁₁ in the R-direction is computed according to the strain-displacement relation (3). From this, and according to the stress–strain equality, the component of stress σ₁₁ is computed. Next, the stiffness matrix S is utilized to evaluate the transverse strains. Then, the elementary displacement field u is determined in terms of strain ε, using an integral operator given by equality (3). The component of the resulting displacement u₁₁ is illustrated in Fig. 5a. Similarly, the same situation occurs for the tangential u₂₂ and longitudinal u₃₃ displacements as shown respectively in Fig. 5b and c. The displacements u₁₁, u₂₂, and u₃₃ in their orthotropic directions grow linearly with each direction. The model demonstrates that it is consistent with respect to the expected calculation proving the Poisson effect consequently. Each displacement varies progressively from the lower value, which is presented by the bottom color of the scale to the higher value given by its top where the force is applied.

Fig. 5 A series of mechanical tests demonstrating the validation of the orthotropic linear-elastic model for the poplar (left), maple (middle), and spruce (right) wood specimens. Each nodal displacement is generated for the uniaxial tensile test.

From a morphological point of view, the deformation intrinsically characterizes the new state of the object. It is necessary to compare the transformation and shape between the original object and the deformed one. For instance, in such fields as neurosurgery,^69,70 the deformed object permits predicting a tumor's growth. Furthermore, the behavior of forward deformation models is investigated in ref. 71 and 72 to predict the variations in human tissue characteristics in response to the movement of bones. The current results are reported and used to verify the proposed method by interpreting the deformation. The main objective is to analyze deeply the deformed paths with the simulation tests in three-dimensional viewing. Only one typical test of these simulations is performed using the uniaxial tensile test. Several influencing factors should be considered when testing configurations are realized within the framework model. These tests contribute to evaluating the accuracy of the deformed field from the original microstructure.

Fig. 6 depicts the microscopic deformation subjected to the orthotropic tensile loadings. Each mechanical test presents a response under uniaxial loading. The loadings follow the radial (R), tangential (T), and longitudinal (L) directions, corresponding respectively to the poplar (see Fig. 6a), maple (see Fig. 6b), and spruce (see Fig. 6c) wood specimens. The deformation provides information about the new structure and its visualization is an important step. The analyzed deformation demonstrates that no enlargement occurs in the transversal plane. The profile in the deformed state differs a little from the original microstructure presented in the initial state since the lateral faces are blocked. The enlargement is analyzed during each test following the uniaxial orthotropic direction normal to this transversal plane. The 3D observations bring a qualitative overview of the poplar and maple microstructures, augmenting the shape and size of the porosity. The poplar's and maple's pores elongate in the principal direction without contracting perpendicularly. However, for the spruce specimen loaded in the longitudinal direction, the distribution of pores becomes uniform. This deformed structure has provided the natural response of this expected elastic behavior. The deformation of the specimen becomes considerable, and its damage happens early when an important displacement field is imposed.

Fig. 6 Typical small deformation for the poplar (left), maple (middle), and spruce (right) wood specimens in each traction testing presented in Fig. 5. The original microstructure is presented in mesh by the orange color, while the deformed microstructure is presented by the blue color.

3.2 Convergence test

Choosing a representative elementary volume (REV) is a crucial step in modeling to correctly describe the material's microstructure. There is not a unique and required concept of the REV for a heterogeneous material with complex microstructures. That might explain the existence of different definitions of the REV, see ref. 73 for a review. Many approaches have been investigated in the literature to decide this size. According to Hill,⁷⁴ the REV must contain sufficient heterogeneity and its size must be much larger than the characteristic length of the microscale and small enough concerning the macroscopic scale of the material. A statistical framework based on the convergence of the effective properties is investigated in ref. 75 and 76 to determine the size of the volume REV. In this contribution, our approach is made by evaluating the visible characteristics of the material by executing various numerical computations on different sizes of the microstructure.

The evaluation of the REV size in computational homogenization is approached by treating mechanical response quantities in terms of mesh type to utilize their properties. A novel approach to determining the convergence of effective properties as a function of domain size is explained as follows. This approach is often recommended in material science, since it converges faster with increasing the number of tetrahedral cells. Its accuracy is demonstrated through three main meshing tests, which are carried out for each wood specimen with a fixed REV size. To achieve higher accuracy and convergence rates of the standard FE method, especially in the presence of complex meshes,^77–79 three types of mesh are investigated. In this contribution, we study the convergence and stability from the mesh research until the mechanical effective quantities become stable. The flowchart of the proposed procedure is shown in Fig. 7. Given a wood sample, we employ multiple mesh types via a progressive refinement strategy, avoiding morphology and topology modification. From a more practical standpoint, the idea touched upon here is considered in the smoothing stabilization and the FE discretization. According to the wood specimen and its REV size, the adaptive mesh is regulated following the refinement process to capture reasonable convergent macroscopic values. Once the suitable computed mesh has converged, the proportionality rule from a given sample resolution is made to determine the convergent mesh without running any simulation.

Fig. 7 Flowchart to ensure the convergence and set-up scheduling the REV size.

The numerical tests are selected by augmenting the grid resolution, presented in Fig. 8 to show the convergence of the numerical method. The pores are maintained to be sufficiently smoothly meshed. The first is obtained with the coarse mesh; the second is realized with the medium mesh; and finally, the last is acquired with the refined mesh. A set of information concerning the three meshes is given in Table 1 for the spruce and poplar specimens. The microstructure of the intermediate spruce species is displayed below. This states that the initial mesh is around 228 791 volumic cells, and the final mesh reaches around 540 729. It also displays the effective elastic modulus in the radial and tangential directions. Since the poplar and maple have the same characteristics of vessels, the procedure testing is repeated sufficiently for the poplar wood with REV size, which is 385 × 385 × 6 voxels.

Fig. 8 An enlargement of coarse and refined mesh zones for the intermediate spruce specimen with dimensions 192 × 192 × 6.

Table 1 The effect of the meshing process on the spruce and poplar samples

Specimen	Mesh type	Cell number	Voxels	Volume fraction	E ^EffT (MPa)	E ^EffR (MPa)
Spruce	Coarse	228.791	192 × 192 × 6	0.367	987	1050
	Medium	311.938	192 × 192 × 6	0.367	900	987
	Refine	540.729	192 × 192 × 6	0.361	862	932
Poplar	Coarse	482.399	380 × 345 × 9	0.307	139	255
	Medium	756.203	380 × 345 × 9	0.307	106	199
	Refine	1.622.184	380 × 345 × 9	0.303	96	222

---

In summary, the above study introduced an adaptive meshing framework that was based on a refinement technique to guarantee convergence to corresponding values. The convergence criteria of the two specimens were estimated following the steps of the invoked strategy. First, these specimens were homogenized in terms of a three-dimensional mesh of the REV size. Three levels of refinements are investigated to converge to the final solution. During the refinement process, the mesh varies by conserving the same morphology and the volume fraction of the aggregates. Then, in the second step, the effective mechanical values were compared. The variation of aggregate meshing from coarse to medium impacts insignificantly until stabilizing at the converged values. As a result, it can be seen that the coarse mesh is sufficient in terms of the existing domain size used in the calculations. The key idea is to inspire from this study to the next one to choose according to the REV size, the convergent mesh. The next study follows further development, particularly complex block calculation with the largest REV, which is much more time-consuming.

In this part, a convergence test was handled to determine the REV sizes whose values satisfy the stability of the effective mechanical property. The computation is carried out in the case of the standard wood directions. The studied REVs are considered within the full 3D image of each sample. P. Perré et al.⁸⁰ have used a similar approach to conduct a convergence test regarding numerical experiments concerning thermal conductivity. In our case of spruce, we will demonstrate that the established REV does not impact the convergence behavior. Consider the user input center (x₀,y₀,z₀) to select the REV geometry. This REV is investigated in parallelepiped form with dimensions (2Δx,2Δy,2Δz), given as follows:

At each test, the volume of the REV is simultaneously increased in the longitudinal direction following the same step Δz = 5, but fixed in the transverse plane (R,T), i.e., Δx = Δy = 96. The information concerning the mesh, such as the computed solid fractions and the volumeric cell numbers, are presented in Table 2. The predicted macroscopic Young's moduli are given in the same table. The convergence is noticed in the middle volumes as the REV size is a bit important. Fig. 9 depicts the effective Young's modulus convergence following the longitudinal direction as the REV size rises. This effective elastic property is unstable on the first volumes. After the sixth volume, it becomes stable since these volumes can be considered representative. This successful convergence test permits us to consider this REV size reasonably well for spruce. Notably, the transverse values of spruce evolve in a quasi-similar way. However, since poplar is presented here as one example of a dual-porosity organization (vessel lumens and fiber lumens), its convergence is open to question.

Table 2 The explanation of the longitudinal effective mechanical properties of the spruce sample in terms of the REV size and information from the mesh

No. REV	Size in voxel	Cell number	Volume fraction	E ^EffR (MPa)	E ^EffT (MPa)
1	192 × 192 × 10	75.263	0.563	2082	1556
2	192 × 192 × 20	81.339	0.563	2189	1706
3	192 × 192 × 30	101.753	0.475	1736	1362
4	192 × 192 × 40	109.242	0.477	1775	1398
5	192 × 192 × 50	115.370	0.479	1804	1427
6	192 × 192 × 60	120.401	0.483	1844	1466
7	192 × 192 × 70	134.473	0.478	1810	1432
8	192 × 192 × 80	151.955	0.479	1814	1431
9	192 × 192 × 90	183.737	0.479	1820	1437
10	192 × 192 × 100	223.411	0.479	1820	1433
11	192 × 192 × 110	274.278	0.476	1806	1414
12	192 × 192 × 120	335.900	0.472	1784	1390
13	192 × 192 × 130	377.461	0.472	1775	1378
14	192 × 192 × 140	445.745	0.472	1782	1377
15	192 × 192 × 150	604.740	0.470	1763	1350
16	192 × 192 × 160	532.384	0.474	1796	1381
17	192 × 192 × 170	557.400	0.489	1831	1453
18	192 × 192 × 180	598.147	0.487	1795	1409
19	192 × 192 × 190	650.010	0.489	1836	1452
20	192 × 192 × 200	699.557	0.486	1790	1399

---

Fig. 9 Convergence of the macroscopic Young's modulus in the case of the spruce wood. This convergence is aimed in the longitudinal direction (L) in terms of REV defined by (16).

We have just assured the convergence behavior of spruce along the longitudinal direction (L). The aim is devoted to surveying the effects of the selected sizes of the REV on the macroscopic mechanical properties of some wood species by carrying out a series of homogenization analyses. The focus will be on the convergence of the porosity from spruce, poplar, and maple wood by increasing the REV size in the transverse plane (R,T). For that purpose, the REV has to be representative of porosity. A configuration on each considered REV for these species with sufficiently suitable size is given by,

The step following the longitudinal direction is fixed for each wood sample, i.e., Δz = 3. Normally, the geometrical center (x₀,y₀,z₀) of the REV is fixed, but it is modified in some cases to capture many vessels. This concerns the smallest element size for poplar and maple samples. After each test, the volume of the REV is concurrently raised in the radial and tangential directions following the uniform steps Δx = Δy = 32. Table 3 reports the identified values of the macroscopic properties via the computational homogenization of the spruce, poplar, and maple species along the radial and tangential directions: this table summarizes the main results of this work. From this set of values, we obtained a data set that gave each sample the macroscopic Young's modulus as a function of the REV size for the various solid fractions obtained by the mesher. Fig. 10 depicts the convergence trends associated with changing REV sizes for the considered samples. As in previous convergence, the homogenized elasticity property is sensitive in the first volume of this series size. It can be seen from this figure that all the homogenized values become close to each other as the REV size increases. This result agrees with the strong and successful convergence test. This convergence is noted in the last volumes as the size of the REV is sufficiently large. Since the shape of vessels for both poplar and maple species is extended in the radial direction, the macroscopic transverse properties are distinguished between these species. This difference is much more significant compared to the spruce wood.

Table 3 The macroscopic Young's modulus and solid fraction of spruce, poplar, and maple in the orthotropic directions as a function of the REV size

Specimen	No. REV	Size in voxels	Cell number	Volume fraction	E ^EffT (MPa)	E ^EffR (MPa)
Spruce	1	192 × 192 × 6	70.811	0.255	331	975
	2	256 × 256 × 6	130.103	0.278	169	381
	3	320 × 320 × 6	198.435	0.298	310	543
	4	384 × 384 × 6	286.267	0.311	203	498
	5	448 × 448 × 6	376.485	0.332	238	445
	6	512 × 490 × 6	497.632	0.320	271	425
	7	576 × 576 × 6	606.485	0.329	286	532
	8	640 × 640 × 6	732.556	0.325	297	578
	9	704 × 704 × 6	870.257	0.320	289	608
	10	768 × 768 × 6	1.039.558	0.318	290	571
Poplar	1	192 × 192 × 6	133.723	0.284	353	482
	2	256 × 256 × 6	247.039	0.301	168	585
	3	320 × 320 × 6	401.630	0.314	204	660
	4	384 × 384 × 6	597.862	0.298	220	562
	5	448 × 448 × 6	838.734	0.344	278	827
	6	512 × 512 × 6	1.049.153	0.335	324	773
	7	576 × 576 × 6	1.278.623	0.319	364	816
	8	640 × 640 × 6	1.590.835	0.322	343	827
	9	704 × 704 × 6	3.176.169	0.346	410	828
	10	768 × 768 × 6	3.653.953	0.345	266	799
Maple	1	192 × 192 × 6	1.609.65	0.313	487	916
	2	256 × 256 × 6	246.450	0.310	583	658
	3	320 × 320 × 6	418.806	0.307	465	864
	4	384 × 384 × 6	556.532	0.337	443	620
	5	448 × 448 × 6	919.130	0.339	523	953
	6	512 × 512 × 6	1.143.100	0.340	527	972
	7	576 × 576 × 6	1.470.580	0.334	479	935
	8	640 × 640 × 6	1.849.768	0.340	522	976
	9	704 × 704 × 6	2.178.739	0.346	500	938
	10	768 × 768 × 6	2.587.940	0.345	492	905

---

Fig. 10 Convergence of the macroscopic Young's modulus in the case of the spruce, poplar, and maple species. This convergence is conducted in the transverse plane (R,T) in terms of REV defined in relation (17).

In light of Table 3, the predicted simulations were performed in the two transverse directions of wood (radial and tangential). The predicted macroscopic values are higher in the radial direction than in the tangential direction for all samples, including the intermediate spruce species. The cellular morphology used here for the intermediate wood differs from other species. However, its solid fraction corresponds more closely to the poplar and maple ones, consequently justifying the interesting outputs. It is clearly seen that the macroscopic mechanical properties depend strongly on the porosity fraction and morphological species. It is explained physically by the dominated solid components on the model at different spatial scales. Note that the convergence has reached a rapid level for the spruce compared to other species.

The mechanical behavior of two wood tissues, poplar and maple specimens, was simulated using the FE method in both the radial and tangential directions. The FE modeling and analytical analysis of energy stress in the microstructure were performed for validation and confirmed the convergence results. This analysis also contributes to the mechanical response, as described in many works.^81,82 Sretenovic et al.⁸¹ investigated the stress transfer from the matrix material to reinforce discontinuous fibers for a wide range of composite materials. For the same purpose, the mechanism of stress transfer has been studied in ref. 82, when the wood is embedded in a polymer matrix to provide a general understanding of its mechanical properties and microscale deformation process. Fig. 11 depicts the transfer properties of the wood species. To clarify this transfer, we adopt a strategy that assumes two REV types can be chosen. The figure is shared into two lines in which the minimal REV presented in the top diagrams is of 192 × 192 × 6 voxels, while the maximal REV presented in the bottom diagrams is of 640 × 640 × 6 voxels. We can see in the top diagrams that the energy grew radially from the left side and gradually extended to the right side of the sample. The same behavior was observed in the bottom diagrams. The energy spreads from the top to the bottom, following the tangential direction. We report that the energy transfer stays valid tangentially for the smallest domain size that complies with the minimum REV size, and radially for well-distributed pores across the microstructure with the maximum REV size. Note that the size of the REV varies spatially and depends on the presented quantity. A better largely quantitative REV size compared to the dimensions of the pores satisfies convergence. This occurs for the selected specimens since they have considerable variability in vessel dimensions compared to the spruce species. The presented simulations are consistent with our visual observation, mainly, there is a relationship between the domain size and convergence.

Fig. 11 Snapshots of radial and tangential transfers of some wood species under tension.

The framework's straightforward application is illustrated in Fig. 11, which shows the formation and development of a stress transfer zone across the microstructure. In all simulation tests, a great variety of high stresses is observed. The red color scales the higher values, whereas the blue color presents the lower values. At first sight, the value of this stress transfer depends on the loading direction. Again, following this analysis, the simulated stress transfer resulting from the FE model demonstrated a good agreement. The transfer starts forming energy zones, increasing this energy from the contact zone, where the stress force is enforced, until the end zone, where it is free. Under these radial (respectively tangential) loading conditions, the transfer propagates so that the radial (respectively tangential) stress always has higher values. In this case, when the radial component is excessive, the tangential component is low, and vice versa. In the same situation, as the load is increased, the radial transfer will be much more important, and the tangential transfer will be reduced significantly.

In what follows, we will carry out several numerical simulations to achieve the validation of the prediction approach by the framework. The homogenization approaches discussed below yield the estimated macroscopic elastic properties from the microscopic structure of heterogeneous materials. Indeed, the microscopic characterization at each tetrahedral cell scale is a necessary step that allows the anticipation of the material's macroscopic behavior. Generally, the homogenized model regarded as a two-phase model made up of a gas brick embedded in a connected solid phase provides the effective transverse modulus E^h expressed as follows:

This novel parallel-series model consisting of a series model and a parallel model is introduced, to compute the averaged mechanical constants in the transverse plane from the averaged stiffness matrix C̄ introduced previously by the relationship (11). This model given by (18) is formulated in terms of E^g_Eff, E^s_Eff, α_g, and α_s which are the effective Young's modulus and volume fraction for gas and solid, respectively, while the subscripts g and s indicate the phase for gas and solid states, respectively. Since the choice of a homogenization scheme to model the wood material depends on its microstructure, the poplar, spruce, and maple wood specimens are investigated. In addition to the intermediate spruce specimen treated through this paper, the earlywood and latewood are considered here to present much variability among wood species. The aim is to enrich this investigation by varying the solid fraction. The analysis of the dimensionless equivalent elasticity leads to illustration of the advantages of the proposed specimens. The present homogenization technique provides the transition of earlywood and latewood information. In particular, its minimal and maximal solid fractions indicate the capability of determining the transition using these specimens. The literature^83–85 offers several contributions concerning the distinction between earlywood and latewood regarding their homogenized properties. From these specimens, a homogenization strategy of only the solid phase is proposed to describe the mechanical and physical properties of macroscopic structures with heterogeneous microstructures. This proposed homogenization is based on the assumption that the micro-stress/micro-strain distribution is homogeneous in the solid matrix. We refer to the macroscopic and microscopic stress measures [small sigma, Greek, macron] and σ, and their associated work-conjugate strains [small epsilon, Greek, macron] and ε in Section 2.4.

A series of mechanical uniaxial traction and compression tests of invoked material samples are carried out to determine the homogenized Young's moduli. During these tests, the lateral faces remain blocked to capture much-averaged values. However, the face in which the force is applied stays free for the mechanical uniaxial testing. The mean homogenized Young's moduli E^h_R, E^h_T, and E^h_L are determined following respectively the radial, tangential, and longitudinal directions. In this study, the REV microstructures are presented with large voxel size. The real brick specimens in three-dimensional space are chosen in such a way as to cover different values of solid fractions. The significant reference of volume is subject to good precision of the macroscopic components. Following the procedure established above, the effective elastic properties of single-phase microstructures without pores are obtained. Fig. 12 provides these results for various microstructure realizations that are presented in terms of volume fractions. First, it can legibly be seen that the REVs with greater volumic fraction have bigger radial and tangential Young's moduli. The effective Young's modulus E^h_L has much higher values compared to homogenized Young's moduli E^h_R and E^h_T, as shown in Fig. 12a and b. These radial and tangential averaged values have a declined tendency of effective material properties due to the influence of material property uncertainties. The radial Young's moduli E^h_R are higher than the tangential Young's moduli E^h_T for all the wood species. Note that the macroscopic value is very close to the series model concerning the spruce earlywood. However, the value becomes close to the parallel model for the spruce latewood. The anisotropic ratio is inverted if the porosity is dilated slightly for these specimens. In this situation, it becomes the opposite due to the morphological reasons and connected phase. Mainly, it is caused by the flattened cells. This means the solid fraction is augmented by around 0.2 for these specimens. Then, the effective elastic modulus is lower in the radial direction than in the tangential direction. These results confirm those obtained by Perré et al.⁸⁶ for thermal/mass diffusivity and mechanical properties of the real wood tissue structure.

Fig. 12 The passage from the microscopic scale to the macroscopic one is realized thanks to the homogenization approach. The effective material properties of wood species in their orthotropic directions are presented in each diagram. The solid line refers to the parallel model.

The consistency of the survey estimated macroscopic values following the longitudinal direction is examined as follows. Fig. 12c and d illustrate the effective macroscopic values E^h_L corresponding approximately to the fundamental parallel model. If these values are exactly the same as the parallel model, then it means that the obtained result can be false. Nevertheless, this approximation is a great sign of result conformity. Generally, these results show consistency concerning our macroscopic numerical predictions of the homogenization aspect. It confirms the capacity framework to forecast the homogenized elastic properties using the FE method.

3.3 Heterogeneity overview of wood

This work aims to develop a realistic heterogeneous micromechanical model of wood by taking into account its microstructures in two regions and applying this model to analyze the effect of wood microstructure on the homogenization aspect. In this regard, the numerical homogenization technique is used to estimate the effective physical properties of dual-phase heterogeneous wood materials.^87,88 The mechanical computational leads bring out the influence of the porosity on the physical properties. This is achieved using two different homogeneous materials. The global microstructure and composition are examined as heterogeneous, where each phase is homogeneous with uniform characteristics. The physical characteristics affecting the solid structure are preserved, as described in part 2.3. Recalling that cellulose and hemicellulose are the two fundamental constituents of wood at the ultrastructural scale, these constituents presented by the second phase are treated as an equivalent material and subjected to 2% for the solid elastic properties. Note that many research works^89–92 have explored the relationship between microstructure with heterogeneous multiple layers and mechanical properties in cellulose-based materials. For instance, a multiscale mechanical framework for analyzing the elastic properties of densified wood, considering its chemical component, is developed in ref. 89 and 90 using analytical and FE approaches. Additionally, the aim of the investigations in ref. 91 and 92 is to study the mechanical properties of cellulose interactions on the anisotropic structure in food and medical applications. Due to the complex imbrication of structural parameters on wood elasticity, our feature of the 2% chain composition leads to simplifying and highlighting the generality of the heterogeneity procedure.

Previously, the classical homogeneous model provided a clear and comprehensive overview of the importance and impact of treating the local wood properties. The heterogeneity approach is useful for understanding which structural parameters control the material properties, such as the realistic stiffness and failure analysis. The mechanical behavior of the anisotropy and heterogeneity (Fig. 13) was simulated using our FE framework in both the radial and tangential directions. Running heterogeneous wood materials is time-consuming and consists of tensile and compression tests at small deformations scale. The investigated materials, such as the poplar and maple species, have regularities in the cell shape and large sizes. Because meshes made with interconnecting pores are required via porous meshes, the FE framework is extended so that the effect of material heterogeneity is explored within the possibilities offered by this model. The effect of wood microstructures on its elastic properties was studied within the systematic computational modeling. In our studies, this model will be used to predict the wood response while mechanical loading is applied. The same figure shows two exhibit diagrams covering the stress tensor at each material point under tensile loading conditions applied to the REV. This contribution establishes a transfer of tensorial stress data, radially and tangentially, respectively, in the left and right diagrams. This result is consistent with the energy transfer illustrated in Fig. 11. The ultimate stress with higher values evolves radially when a tensile force is applied in this direction, and vice versa while following another direction. The improved framework contributes through the porous structure a weak mechanical energy during the deformation. From a physical point of view, one must expect to obtain a smoother response, rather than this non-variability. This is explained by attributing the low mechanical properties and other characteristics of the pore structure.

Fig. 13 (a) and (b) Design of respectively the radial and tangential tensile stresses for maple and poplar specimens via heterogeneous FE modeling.

We believe that the present framework can provide further investigations on complex heterogeneous wood materials. More numerical simulations are demanded to provide information on the basic constituents. The homogenization-based multi-scale techniques are adopted to physically treat the intricate mechanisms of wood at the macroscopic level. The results provided by the homogeneous and heterogeneous wood materials via the FE method are performed on a suitable REV.⁹³ A series-parallel mixture model via relation (18) is proposed to predict the overall properties. As shown in Fig. 13, the two elastic materials are connected together in parallel to form the investigated wood material. To further understand the effect of cellulose components on the mechanical properties of wood, we have compared the average tensile Young's modulus along the anisotropic direction (Fig. 14). To this aim, recall that the overall Young's modulus is the average of the moduli of the constituents weighted by the volume fraction of each phase.^94–96 The assumed porosity highly influences the prediction of effective properties in the case of heterogeneous against the homogeneous framework.^97–99 First, it can clearly be seen that the REVs with greater porosities have smaller predictive effective elastic properties that must remain in the physical ranges. At the same time, the radial Young's moduli E^h_R are greater than the tangential Young's moduli E^h_T for all heterogeneous samples. This remark is expanded for the spruce specimen. As expected, the radial and tangential heterogeneous homogenized Young's moduli E^h_R and E^h_T are markedly greater than the homogeneous ones. The wood cellulose effect demonstrates superior performance over homogeneous, so practically, higher values indicate better mechanical properties of heterogeneous materials. We conclude that the results obtained via the presented methodology with the FE tool yield simulation of any wood microstructure with great accuracy.^100–102

Fig. 14 An upscaling strategy, such as the homogenization approach, is used to obtain the effective material properties of wood species over all phases compared with the homogeneous model at the macroscale.

4 Discussion and concluding remarks

The wood elements are still known for multi-scale and anisotropic aspects. The scientific interest of this work consists of understanding the wood microstructures utilizing the X-ray tomography technique. This technique is based on scan parameters with specific environmental conditions to access their real morphology. The poplar, spruce, and maple wood specimens with representative dimensions are scanned using the tomograph and analyzed. Next, image processing is carried out to provide a 3D digital morphology of these samples. The aim is to use this representation to generate a 3D unstructured mesh that is the necessary input data for the finite element method (FEM) framework.

In the scope of computer science, we have shown that the FEM is a numerical procedure that aims to describe the physical structure in a three-dimensional space, respecting the mathematical approach. This framework of FEM was consistent with the structure topologies of different phases species. The mechanical model of the 3D real structure has been further validated, through a numerical procedure. In particular, the numerical results highlight the effect of wood's anisotropy at the microscopic scale. The higher mechanical properties of these tested samples are presented by showing the stress–strain curves. The summarized results for each sample provide a very realistic microstructure deformation.

Solving the linear system in parallel was implemented to exemplify the improvement of the computational efficiency and speed of the code on the computer architecture. Indeed, the post-processing of mechanical modelling of these microstructures can be very time consuming. The coefficient of the global matrix is presented in a sparse format, which makes this method more accurate in performing small deformation scales. We would like to show that the parallel computing using Message Passing Interface (MPI), and Compute Unified Device Architecture (CUDA), can greatly increase the computational efficiency of the modelling framework. The presented work will be extended for multi-scale simulation of 3D fiber-reinforced composite modelling. Then we will investigate the influence of fiber distribution on the simulation, in which multiple fibers with different orientations are presented. We will carry out experimental tests using a universal testing machine. The studied specimens with fibers will be placed in this machine according to the standard of these structures.

Author contributions

Tarik Chakkour: conceptualization; methodology; validation; formal analysis; investigation; data curation management; writing – original draft; visualization. Patrick Perré: conceptualization; methodology; formal analysis; investigation; resources; data curation management; writing – review and editing; supervision; project administration; funding acquisition.

Consent to participate

All authors contributed to this work. All authors read and approved the final manuscript.

Data availability

Data will be made available on request.

Conflicts of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

First, the authors address many thanks to the anonymous reviewers for their helpful and valuable comments that have greatly improved the paper. This study was carried out in the Centre Européen de Biotechnologie et de Bioéconomie (CEBB), supported by the Région Grand Est, Département de la Marne, Greater Reims and the European Union. In particular, the authors would like to thank the Département de la Marne, Greater Reims, Région Grand Est and the European Union along with the fund (FEDER Grand Est 2021–2027) for their financial support of the Chair of Biotechnology of CentraleSupélec.

Notes and references

1. J. Wang, X. Han, W. Wu, X. Wang, L. Ding, Y. Wang, S. Li, J. Hu, W. Yang and C. Zhang, et al. , Int. J. Biol. Macromol., 2023, 231, 123343.
2. M. Chutturi, S. Gillela, S. M. Yadav, E. S. Wibowo, K. Sihag, S. M. Rangppa, P. Bhuyar, S. Siengchin, P. Antov and L. Kristak, et al. , Sci. Total Environ., 2023, 864, 161067.
3. S. M. Koch, C. Goldhahn, F. J. Muller, W. Yan, C. Pilz-Allen, C. M. Bidan, B. Ciabattoni, L. Stricker, P. Fratzl and T. Keplinger, et al. , Mater. Today Bio, 2023, 22, 100772.
4. Y. Chen, L. Zhang, C. Mei, Y. Li, G. Duan, S. Agarwal, A. Greiner, C. Ma and S. Jiang, ACS Appl. Mater. Interfaces, 2020, 12, 35513–35522.
5. S. He, X. Zhao, E. Q. Wang, G. S. Chen, P.-Y. Chen and L. Hu, Annu. Rev. Mater. Res., 2023, 53.
6. P. Niemz and M. Dunky, Biobased Adhesives: Sources, Characteristics and Applications, 2023, pp. 621–658.
7. T. Chakkour, Eng. Math. Lett., 2018, 4.
8. J. Wang, D. Zhang and F. Chu, Adv. Mater., 2021, 33, 2001135.
9. Y. Zuo, J. Feng, T.-O. Soyol-Erdene, Z. Wei, T. Hu, Y. Zhang and W. Tang, Chem. Eng. J., 2023, 463, 142332.
10. F. Aldakheel, E. S. Elsayed, T. I. Zohdi and P. Wriggers, Comput. Mech., 2023, 1–17.
11. R. Yazdanparast and R. Rafiee, Eng. Comput., 2023, 39, 373–386.
12. J. Liu, W. Liu, R. Tan, W. Huang, B. Wang, S. Zhang and M. Zhang, Drying Technol., 2023, 41, 746–753.
13. A. Liszka, R. Wightman, D. Latowski, M. Bourdon, K. B. Krogh and J. J. Lyczakowski, Front. Plant Sci., 2023, 14, 1283093.
14. G. K. Zelenov, L. V. Belokopytova, E. A. Babushkina, D. F. Zhirnova, B. Yang, X. Peng, J. Liu, G. A. Sitnikov and E. A. Vaganov, Forests, 2024, 15, 249.
15. X. Zhang, H. Liu and T. Rademacher, Sci. Total Environ., 2024, 912, 169165.
16. J. J. Camarero, X. Serra-Maluquer, M. Colangelo, A. Gazol and M. Pizarro, IAWA J., 2023, 1, 1–14.
17. A. Slupianek, E. Myskow, A. Kasprowicz-Maluski, A. Dolzblasz, R. Zytkowiak, M. Turzanska and K. Sokolowska, J. Exp. Bot., 2023, erad469.
18. Ü. Büyüksar, N. As and T. Dündar, BioResources, 2017, 12, 4004–4012.
19. W. Li, Z. Zhang, C. Mei, P. Kibleur, J. Van Acker and J. Van Den Bulcke, Wood Mater. Sci. Eng., 2022, 1–10.
20. R. Lehnebach, M. Campioli, J. Gricar, P. Prislan, B. Mariën, H. Beeckman and J. Van den Bulcke, Front. Plant Sci., 2021, 12.
21. J. L. Paris and F. A. Kamke, Int. J. Adhes. Adhes., 2015, 61, 71–80.
22. X. Zhang, B. Xiao, H. Andrä and Z. Ma, Compos. Struct., 2016, 137, 18–32.
23. S. Dai and P. Cunningham, Compos. Struct., 2016, 142, 298–312.
24. A. Catapano and M. Montemurro, Compos. Struct., 2014, 118, 664–676.
25. K. Spanos, S. Georgantzinos and N. Anifantis, Compos. Struct., 2015, 132, 536–544.
26. Z. Wang, T. Wang, W. Wang and Z. Zhang, Acta Geotechnica, 2023, 18, 1883–1900.
27. S. Aney and A. Rege, Math. Mech. Solids, 2023, 28, 1624–1634.
28. C. González and J. LLorca, Compos. Sci. Technol., 2007, 67, 2795–2806.
29. L. Riaño, L. Belec and Y. Joliff, Compos. Struct., 2016, 154, 11–16.
30. S. Balay, S. Abhyankar, M. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W. Gropp, et al., 2019, https://ora.ox.ac.uk/objects/uuid:fa2b9e7c-1c58-429c-90fd-f780a3c3dc7d.
31. Y. Saad, SPARSKIT: a basic tool kit for sparse matrix computations, 1994.
32. S. L. Wood, K. Jacobson, W. T. Jones and W. K. Anderson, AIAA Scitech 2020 Forum, 2020, p. 0317.
33. M. Othman, Bull. Am. Phys. Soc., 2024 , https://meetings.aps.org/Meeting/APR24/Session/J16.3.
34. W.-L. Feng, J. Zhong, T. Chen and X.-F. Yuan, J. Phys.: Conf. Ser., 2023, 012012.
35. Z. Wang, T. A. Palmer and A. M. Beese, Acta Mater., 2016, 110, 226–235.
36. M. Z. Ahmed, K. Chadha, S. Reddy, D. Shahriari, P. P. Bhattacharjee and M. Jahazi, Metall. Mater. Trans. A, 2020, 51, 6406–6420.
37. H. Yue, Y. Chen, X. Wang and F. Kong, J. Alloys Compd., 2018, 750, 617–625.
38. R. T. Mills, M. F. Adams, S. Balay, J. Brown, A. Dener, M. Knepley, S. E. Kruger, H. Morgan, T. Munson and K. Rupp, et al. , Parallel Comput., 2021, 108, 102831.
39. J. D. Trotter, J. Langguth and X. Cai, Parallel Comput., 2023, 118, 103051.
40. T. Chakkour, Front. Comput. Sci., 2024, 5.
41. T. Chakkour, Front. Comput. Sci., 2024, 6.
42. V. Hapla, M. G. Knepley, M. Afanasiev, C. Boehm, M. van Driel, L. Krischer and A. Fichtner, SIAM J. Sci. Comput., 2021, 43, C127–C153.
43. P. Jolivet, J. E. Roman and S. Zampini, Comput. Math. Appl., 2021, 84, 277–295.
44. M. A. Salazar de Troya and D. A. Tortorelli, Struct. Multidiscipl. Optim., 2020, 62, 2467–2479.
45. J. Schindelin, I. Arganda-Carreras, E. Frise, V. Kaynig, M. Longair, T. Pietzsch, S. Preibisch, C. Rueden, S. Saalfeld and B. Schmid, et al. , Nat. Methods, 2012, 9, 676–682.
46. C. A. Schneider, W. S. Rasband and K. W. Eliceiri, Nat. Methods, 2012, 9, 671–675.
47. F. De Chaumont, S. Dallongeville, N. Chenouard, N. Hervé, S. Pop, T. Provoost, V. Meas-Yedid, P. Pankajakshan, T. Lecomte and Y. Le Montagner, et al. , Nat. Methods, 2012, 9, 690–696.
48. N. Otsu, IEEE Trans. Syst. Man Cybern., 1979, 9, 62–66.
49. C. A. Beretta, S. Liu, A. Stegemann, Z. Gan, L. Wang, L. L. Tan and R. Kuner, Cells, 2023, 12, 704.
50. E. Kugler, E.-M. Breitenbach and R. MacDonald, Curr. Protoc., 2023, 3, e654.
51. T. Chakkour, Oxford Open Mater. Sci., 2024, 4.
52. C. Zhang, M. Chen, S. Keten, B. Coasne, D. Derome and J. Carmeliet, Sci. Adv., 2021, 7, eabi8919.
53. C. Zhang, M. Chen, B. Coasne, S. Keten, D. Derome and J. Carmeliet, Composites, Part B, 2022, 228, 109449.
54. A. Alkadri, D. Jullien, O. Arnould, E. Rosenkrantz, P. Langbour, L. Hovasse and J. Gril, Wood Sci. Technol., 2020, 54, 1269–1297.
55. P. Grönquist, P. Panchadcharam, D. Wood, A. Menges, M. Rüggeberg and F. K. Wittel, R. Soc. Open Sci., 2020, 7, 192210.
56. R. C. Neagu and E. K. Gamstedt, J. Mater. Sci., 2007, 42, 10254–10274.
57. N. Charupeng and P. Kunthong, Int. J. Comput. Mater. Sci. Eng., 2022, 11, 2150027.
58. M. Stucu, Int. J. Solids Struct., 1992, 29, 197–213.
59. A. Düster, J. Parvizian, Z. Yang and E. Rank, Comput. Methods Appl. Mech. Eng., 2008, 197, 3768–3782.
60. X. Shen, Y. Xue, Q. Yang and S. Zhu, Adv. Appl. Math. Mech., 2022, 14, 365–385.
61. Z. Yu, Y. Cui, Q. Zhang, J. Liu and Y. Qin, Mech Mach Theory, 2022, 174, 104906.
62. B. Zheng, T. Li, H. Qi, L. Gao, X. Liu and L. Yuan, Int. J. Mech. Sci., 2022, 223, 107282.
63. W. Wang, Y. Wu, H. Wu, C. Yang and Q. Feng, Soil Dyn. Earthq. Eng., 2021, 140, 106420.
64. T. Zhu and S. Atluri, Comput. Mech., 1998, 21, 211–222.
65. M. Königsberger, M. Lukacevic and J. Füssl, Mater. Struct., 2023, 56, 13.
66. M. Königsberger, V. Senk, M. Lukacevic, M. Wimmer and J. Füssl, Composites, Part B, 2024, 281, 111571.
67. Z. Guo, X. Liu, L. Huang, S. Adhikari and X. Liang, Compos. Struct., 2024, 118391.
68. J.-F. Ganghoffer, A. Wazne and H. Reda, Mech. Res. Commun., 2023, 130, 104114.
69. M. Tonutti, G. Gras and G.-Z. Yang, Artif. Intell. Med., 2017, 80, 39–47.
70. G. M. Della Pepa, G. Menna, V. Stifano, A. M. Pezzullo, A. M. Auricchio, A. Rapisarda, V. M. Caccavella, G. La Rocca, G. Sabatino and E. Marchese, et al. , Neurosurg. Focus, 2021, 50, E15.
71. N. H. Hart, S. Nimphius, T. Rantalainen, A. Ireland, A. Siafarikas and R. Newton, J. Musculoskelet. Neuronal. Interact., 2017, 17, 114.
72. P. Zioupos, H. O. Kirchner and H. Peterlik, Bone, 2020, 131, 115176.
73. I. Gitman, H. Askes and L. Sluys, Eng. Fract. Mech., 2007, 74, 2518–2534.
74. R. Hill, J. Mech. Phys. Solids, 1963, 11, 357–372.
75. D. Savvas, G. Stefanou and M. Papadrakakis, Comput. Methods Appl. Mech. Eng., 2016, 305, 340–358.
76. X. Tu, A. Shahba, J. Shen and S. Ghosh, Int. J. Plast., 2019, 115, 268–292.
77. S. Bernaud, I. Ramière, G. Latu and B. Michel, Ann. Nucl. Energy, 2024, 205, 110577.
78. F. Freddi and L. Mingazzi, Appl. Eng. Sci., 2023, 14, 100127.
79. G. Liu, W. Huang, N. Leng, P. He, X. Li, M. Lin, Z. Lian, Y. Wang, J. Chen and W. Cai, Bioengineering, 2024, 11, 519.
80. B. Mazian, E.-H. Quenjel and P. Perré, Int. J. Heat Mass Transfer, 2023, 217, 124700.
81. A. Sretenovic, U. Müller and W. Gindl, Composites, Part A, 2006, 37, 1406–1412.
82. E. Jungstedt, C. Montanari, S. Östlund and L. Berglund, Composites, Part A, 2020, 133, 105853.
83. L. Melzerová, L. Kucková, T. Janda and M. Šejnoha, et al. , Wood Res., 2016, 61, 861–870.
84. M. Šejnoha, L. Kucková, J. Vorel, J. Sykora and W. De Wilde, Int. J. Comput. Methods Exp. Meas., 2019, 7, 167–180.
85. A. Rafsanjani, D. Derome, F. K. Wittel and J. Carmeliet, Compos. Sci. Technol., 2012, 72, 744–751.
86. P. Perré, G. Almeida, M. Ayouz and X. Frank, Ann. For. Sci., 2016, 73, 147–162.
87. Y. Wang, G. Han, X. Liu, Y. Ren and H. Jiang, Compos. Struct., 2024, 328, 117724.
88. C. Organista, R. Tang, Z. Shi, K. Jefimovs, D. Josell, L. Romano, S. Spindler, P. Kibleur, B. Blykers and M. Stampanoni, et al. , Sci. Rep., 2024, 14, 384.
89. R. Song, F. Deng, X. Liang, J. Song, S. Shen and T. Li, J. Mech. Phys. Solids, 2024, 191, 105761.
90. N.-T. Phan, F. Auslender, J. Gril and R. Moutou Pitti, Wood Sci. Technol., 2024, 1–29.
91. S.-Q. Chen, P. Lopez-Sanchez, D. Mikkelsen, M. Martinez-Sanz, Z. Li, S. Zhang, E. P. Gilbert, L. Li and M. J. Gidley, Food Hydrocolloids, 2023, 136, 108283.
92. X. Lu, H. Jiao, Y. Shi, Y. Li, H. Zhang, Y. Fu, J. Guo, Q. Wang, X. Liu and M. Zhou, et al. , Carbohydr. Polym., 2023, 308, 120669.
93. S. Okamoto, N. Akiyama, Y. Araki, K. Aoki and M. Inayama, J. Wood Sci., 2021, 67, 1–25.
94. S. Sugiman, S. Salman and M. Maryudi, Mater. Today Commun., 2020, 24, 101360.
95. Y. K. Galadima, W. Xia, E. Oterkus and S. Oterkus, Eng. Comput., 2023, 39, 461–487.
96. C. Naili, I. Doghri, T. Kanit, M. Sukiman, A. Aissa-Berraies and A. Imad, Compos. Sci. Technol., 2020, 187, 107942.
97. R. Sharma, V. K. Mittal and V. Gupta, J. Inst. Eng. (India): Ser. D, 2024, 105, 1–20.
98. H. A. Madkhali, M. Nawaz and S. O. Alharbi, Ain Shams Eng. J., 2024, 15, 102449.
99. S. Yakupov, H. Kiyamov, N. Yakupov and I. Mukhamedova, Case Stud. Constr. Mater., 2023, 19, e02360.
100. J. Tao, P. Tahmasebi, M. A. Kader, D. Feng, M. Sahimi, P. D. Evans and M. Saadatfar, Mater. Des., 2023, 228, 111812.
101. C.-W. Chang, C.-Y. Lin, H.-W. Kao and F.-C. Chang, Wood Sci. Technol., 2023, 57, 75–92.
102. R. Bengtsson, M. Mousavi, R. Afshar and E. K. Gamstedt, Mech. Mater., 2023, 179, 104586.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4ma00554f

---