Numerical simulation of transverse compression and densification of wood

Densification is commonly adopted to increase the mechanical performance of wood, but research on the micromechanical behaviour of the material during transverse compression is limited. Robust numerical models will enable better predictions of the performance of wood during compression and optimise the manufacturing process of densified wood minimising experimentation. The densification stress–strain response of wood after chemical treatment is reported via numerical simulations. A 3D finite element model of wood microstructure is studied under transverse compression using ABAQUS/Explicit software. A lower cellulose, hemicellulose and lignin content in the chemically treated wood is considered in the material parameters of the cell wall, and an ideal elastoplastic material model is used to represent the nonlinear stress–strain response. Parametric studies regarding the cell wall thickness, yield stress and chemical treatment are also considered. The numerical predictions agree well with microscopy studies of densified wood, and the nominal stress–strain curve obtained is similar to experimental findings under transverse compression as found in the literature. The cell wall thickness and yield stress are found to significantly affect the compressive stress–strain response of wood.


Introduction
The excellent mechanical properties of wood are determined by the microstructure of the wood cells (Chen et al. 2020;Toumpanaki et al. 2020). The wood cell wall is typically divided into five layers, and microfibrils of different layers are geometrically arranged in different orientations (Toumpanaki et al. 2020). So the wood cell wall is a natural fibre-reinforced composite material with a hierarchical organisation and has been demonstrated as a smart material in buildings (Eichhorn et al. 2022). Although fibre-reinforced composite materials can exhibit high-performance properties, the large amount of waste generated during the manufacture process, the difficulties in recycling, and their non-biodegradable characteristics mean that they are not always sustainable (Ramakrishna et al. 2001). Wood is characterised as being both a renewable and sustainable composite material. The material exhibits mechanical properties comparable to synthetic composites, and for its weight, commonly used materials like steel (Chen et al. 2020;Toumpanaki et al. 2020).
Wood cells in their most basic form comprise cell walls and lumens, and they are arranged in such way to form a honeycomb structure (Chen et al. 2020;Toumpanaki et al. 2020). The cell wall is a nano-scale composite material comprising cellulose, hemicellulose and lignin (Côté 1968;Isogai et al. 1989), where the lignin and hemicellulose are the "matrix" components, and cellulose is the "reinforcement". Wood cell walls are multi-layer composite material structures. Microscopically, their structure can be divided into three layers: the middle lamella (ML), primary wall (P), and secondary wall (S). The material properties of the secondary wall, which is divided into an outer layer (S1), a middle layer (S2) and an inner layer (S3), are anisotropic. The mechanical properties of wood are attributed to its cellular porous solid structure, and the molecular arrangement of its polymer components within the S2 layer plays an important role when the wood is loaded parallel to the cell grain direction (Côté 1968). The anisotropy of the microfibrils, mainly in the S2 layer, determines the orientationdependent mechanical properties of wood cells, and thereby the wood itself.
Modifying wood via heat treatment and transverse compression is a well-known method to improve its shear strength, surface hardness and durability. Since increasing the density of wood improves its mechanical properties, many attempts have been made to develop the densification process. Rautkari et al. (2013) compressed samples of Scots pine with different thicknesses, studying their density distribution and Brinell hardness after hydrothermal post-treatment (heated at 140 °C for 30 min and compressed with a 1000 N load for 15 s). They showed that the more compressed samples had a higher surface densification, resulting in an increase in the density of the surface layer and a significant increase in the elastic recovery of the densified wood.
However, a lack of dimensional stability is the main disadvantage of most thermomechanical treatments in densified wood, and studies have tried to address this problem either with resin impregnation or secondary treatment at high temperatures (Gabrielli and Kamke 2010;Schwarzkopf 2021;Gong et al. 2010). A recent study by Song et al. (2018) proposed a new approach for manufacturing a high-performance and dimensionally stable densified wood. Their densified wood changed the microstructure from the native state using chemical treatment (an aqueous mixture of NaOH and Na 2 SO 3 ) and a hot-pressing process at 100 °C (Song et al. 2018). Under the chemical treatment, the amount of hemicellulose and lignin in the cell wall was reduced, which caused the cell wall to become more readily compressible. After hot-pressing, the porous channel honeycomb structure became densified, and all the cell walls collapsed. The thickness of the wood was reduced by 80%, and the density was increased threefold (Song et al. 2018). Song et al. (2018) managed to achieve a tensile strength of the densified wood of up to ~ 600 MPa, which was 11 times that of untreated wood. The hardness was also found to be 30 times higher than that of the native undeformed wood (Song et al. 2018). These increases in mechanical properties were attributed to the higher degree of hydrogen bonding between the highly aligned cellulose fibres (Song et al. 2018).
Despite these great advances by Song et al. (2018), little research on the mechanistic understanding of wood compression and densification has been published. In the densification process, the stress transfer between the layers of the wood cell walls is critical. The material failure that occurs during deformation also affects the quality and magnitude of the mechanical properties of the densified wood (Salmén 2004).
Given the complexities of the wood structure, finite element (FE) modelling is a useful tool for understanding the stress transfer between layers, the stress distribution and stress-strain behaviour of wood microstructures during the compression and densification process. FE modelling can also be informative on how to optimise the densification procedure to avoid local stress concentrations. FE tools to model wood can be adopted at different scales ranging from micro-(cell wall layer) to macro-scale (wooden beam element).
Lukacevik et al. (2017) proposed a multi-surface failure criterion for brittle and ductile failures of earlywood and latewood to predict the failure mechanisms at the annual ring scale. With their proposed method, a new multi-surface failure criterion for clear wood was developed and validated by experiments under bi-axial loading. A new algorithm was presented which was able to describe both plastic failures and cracks. However, the numerical multiscale concept was initiated at the ring scale, and the layered structure of wood cell walls was not considered.
De Magistris and Salmén (2008) established a FE model of the wood microstructure, and numerically analysed the deformation of wood cell walls transverse to the fibre axis using the ABAQUS software. They created a representative volume model of wood microstructure based on a scanning electron microscope (SEM) image. The model consisted of twelve wood cells (each row consisted of four cells arranged in three rows). Each cell was divided by four layers of P, S1, S2 and S3, and each cell was connected by an M layer. The orthotropic elastic material parameters of each layer of cell wall were calculated using the rule of mixtures combined with a semiempirical model (Halpin-Tsai) (Halpin and Kardos 1976). The elastic properties of the dry cellulose, hemicellulose and lignin were taken from the work of Bergander and Salmén (2000a). Three load conditions of compression, shear and combined shear compression were applied to this microstructural model. The local and global buckling modes of the wood cell walls under different loading conditions and the stress concentration area of the cell wall deformation were studied. It was found that the cell wall buckled under a compressive load, and all the wood cells "bulged" in the same direction, with the earlywood showing a characteristic "S" shape. Under shear and compression loads, the cell wall exhibited a "brick" shape which resembles a compressed morphology as seen by light microscopy (De Magistris and Salmén 2008). However, the nonlinear material properties, the large deformation of the cell wall, and the stress distribution in the walls were not considered. Fortino et al. (2015) established a microstructural model of wood and reported the stress-strain curve for the cell wall deformation by means of a three-dimensional model under transverse load. They used the same theory as de Magistris and Salmèn (2008) to calculate the material parameters and studied the effect of cell geometry, temperature, and humidity on earlywood and latewood, in the stress distribution under transverse loads. It was found that the higher the temperature and humidity, the smaller the required load during the compression process. In addition, by means of the simulation results with different cell angles, it was found that as the cell angle decreases, the total elastic modulus rises, requiring a larger load for compression (Fortino et al. 2015). According to the experimental results, earlywood was found to be more prone to collapse under the same pressure and the stress gradient of the S2 layer for all models increased towards the S1 layer. The S1 and S3 layers also peeled off from the S2 layer due to the tensile stress (Fortino et al. 2015). This model was able to calculate the stress transfer between the cell wall layers and gave an explanation of the response of the secondary wall (S1, S2, and S3 layers) when the stress increased, highlighting its importance in the mechanics of wood. However, the simulated compression process was limited to elastic and small deformations only. Wang et al. (2019) established a FE model of a wood microstructure comprising twelve cells to study the relationship between the different cell wall layers and the fracture of the material under tension and shear loading. They analysed the stress concentrations and compared the crack position of a cell wall, as observed with a microscope. They discovered that the initial fracture position of the cell occurs at the edges of the S2 layer under longitudinal tension and at the S1/S2 interface under inplane shear deformation. Nairn (2006) used the material point method (MPM) to study the elastic and plastic properties of softwood and hardwood in both radial and tangential compression, also considering the effects of anatomy. In this work, the cell-wall material was treated as a transversely isotropic continuum and the numerical model ignored the effects of the microstructure of the cell wall. Plane strain models were established for numerical modelling. The cell wall material was assumed to be elastic-plastic without hardening, and a von Mises criterion was used as the material is isotropic in the plane strain model. However, the wood cell wall is known to consist of several layers, and each layer is orthotropic with a lamina structure. The anisotropic structure of the cell wall may influence the deformation characteristics of the wood during compression and densification.
All the previously published FE simulation studies considered only elastic properties and small deflections except Nairn's simulation study (Nairn 2006). His 2D elastoplastic model ignored the microstructural anisotropy of the cell wall, but the MPM method can model densification and large deformations without numerical difficulties and the true wood anatomy can be easily modelled based on SEM images. In the process of transforming wood into densified wood, all the cell walls will undergo very large deformations and collapse under large displacements. During this process, the plasticity of the cell walls must be considered. In the present study, a microstructural FE model is established to study the compression and densification process of wood, taking into account the plasticity of the material and the geometric nonlinearity of the cell wall during compression.

Basic geometry
To simulate the compression and densification of wood, a microstructural model of the wood for numerical simulation was set up as shown in Fig. 1a. An example of the microstructure of wood, as observed using microscopy, is shown in Fig. 1b for comparison. The model is divided into three rows, with four cells with different sizes in each row. These cells are set asymmetrically to represent the inherent variability among cell walls and better demonstrate the buckling deformation modes during compression. There are three larger cells and nine smaller cells in the model. The walls of each cell comprise four layers, namely P, S1, S2 and S3. All cells were connected with an M layer as shown in Fig. 1a. To describe the model deformation, the twelve cells in the model were numbered from 1 to 12. There were two sizes to the cells: one being 30 µm × 30 µm and the other 30 µm × 40 µm. The dimensions of the cell lumen were set to be representative of a softwood structure (Piermattei et al. 2020). To avoid stress concentrations, a chamfer radius of each cell was set to be 5 µm. The connecting layer M and the primary wall P layer were regarded as an oriented laminate. The material directions of the secondary wall S1, S2, S3 layers were determined by the microfibrillar arrangements (Fig. 1c). In this model, the microfibril angle of the S3 layer was set at −50°, and those of the S2 and S1 layers were, respectively, 10° and 70°. Finally, the fibre directions of P and M layers were set at 0° (Wang et al. 2019).
To compare the effects of cell wall thickness on the compression properties of the wood, two different wood cellular microstructural models with different total cell thicknesses were studied, namely 2.82 µm and 5.50 µm. The former could be representative of earlywood and the latter of latewood (Piermattei et al. 2020). However, the cell lumen size is more representative of earlywood. Their parameters are summarised in Supplementary Information (see Table S1, Supplementary Information). The heights of the wood models for earlywood and latewood were, respectively, 91.8 µm and 92.4 µm, and the widths were 132.3 µm and 133.0 µm. Because the fibre directions of the layers in the wall are different, it is not appropriate to simplify the models as either a plane stress or plane strain problem. Therefore, a 3D model was established, and the length of the two models was set to be 30 µm which is same as Salmén's work (Salmén 2004).

Material model and parameters
In this study, the layered structure of the cell wall is considered, and the material is also assumed to be ideally elastoplastic.
According to the microstructures of the S1, S2, S3 and P layers, each of these layers is orthotropic and the material parameters of each layer were calculated by the rule of mixtures. The M layer is considered as an isotropic material. For an in-plane random long fibre lamina within the P layer, the global elastic properties were determined according to Hull and Clyne (1996). The elastic properties calculated for the three directions of the untreated wood cell layers are given in Table S2 (Supplementary Information). The elastic properties of cellulose, hemicellulose and lignin were derived from previous approaches (Bergander and Salmén 2000b;Thuvander et al. 2002) and are summarised in Table S3 (Supplementary Information). The volume content and weight ratio of every polymer component in the wood layers were taken from Bergander and Salmén (2000b) as summarised in Table S4 (Supplementary  Information).
As reported by Song et al. (2018) before the wood is compressed, it is treated with chemical reagents. To make the wood cell wall readily compressible, the cellulose, hemicellulose and lignin are partially removed from each layer during the chemical treatment process. Therefore, the material parameters of wood after chemical treatment of each layer are different to those of the untreated wood. Song et al. (2018) tested the content of compressed wood after chemical and high temperature treatment; the lignin content was reduced by 45%, the cellulose content was reduced by 12%, and the hemicellulose content was reduced by 74%. So, a wood structure with 45% of the lignin removed was also used in the present work. Table S5 (Supplementary Information) gives the contents of cellulose, hemicellulose, and lignin of wood after chemical treatment. DW-x refers to the densified wood with a certain amount (x) of lignin removal and subsequent densification.
So, according to these reductions in cell wall components, we can estimate the material parameters of each layer of the cell wall after chemical treatment from the material parameters of untreated wood. The weight ratio W C and volume ratio V C of wood after chemical treatment of the cellulose, hemicellulose and lignin are listed in Table S6 (Supplementary Information).
The elastic properties of each layer were determined by using the equations: where E 1 is the elastic modulus along the microfibril growth direction, E 2 is the transverse elastic modulus, G 12 is the shear modulus, and υ 12 is Poisson's ratio. V f represents the fibre volume fraction, and V mL and V mH represent the volume fractions of lignin and hemicellulose, respectively. In addition, G m and G f are the shear moduli of the matrix and cellulose fibrils, respectively. η 2 is the stress partitioning parameter, and η 12 is the stress distribution parameter which is determined by For an in-plane random long fibre lamina, like the P layer, the global elastic properties were determined according to Hull and Clyne (1996) using the equations: Using Eqs. 1-9 and the values from Table S3 to S6 (Supplementary Information), the parameters of the different layers of a cell wall are shown in Table S7 and more details on the modelling approach are also given in Supplementary Information. An ideal elastoplastic model for the cells is adopted in the present numerical simulation. So, under elastoplastic conditions, the yield stress is a main factor in determining the stress-strain relationships. To understand the effect of the cell wall yield strength assumption in the stress-strain curve, a parametric study was carried out with yield stresses in the range of 30-50 MPa and 100-300 MPa; these values are within the range of yield stresses considered in the work of Nairn for loblolly pine (Nairn 2006). Yield stress values in a similar range were subsequently adopted by Aimene and Nairn (2015) for poplar, where hyperelastic-plastic material parameters were used to simulate the large-deformation of transverse wood compression.
According to the experimental results of Song et al. (2018), the density of untreated wood is 0.43 g cm −3 , and the density of densified wood with 45% lignin removal is 1.3 g cm −3 . The density of the densified wood is primarily the outcome of the reduced volume due to the completely compressed cell lumens and the reduced weight of the chemically treated cell walls. A density of 1.5 g cm −3 was used for the wood in the present model (Table 1).

Boundary conditions and contact interactions
To simulate the full compression and densification process, two rigid plates with sizes of 200 µm × 50 µm × 1.0 µm were placed on the upper and lower surfaces of the model as shown in Fig. 2 Further details on the finite element (FE) model are contained within Supplementary Information, including the definitions of the geometries and orientations within the fibre cell walls (Fig. S1), and the simulation of the compression and densification process including the control of the displacement and internal and kinetic energies (Fig. S2).

Validation of numerical simulation
The deformation of the wood cells at a displacement of -25 µm is shown in Fig. 3a as derived from the present numerical analysis for Model 1. In the initial stage of the compression process, the stress does not reach the yield stress, and the cell walls undergo elastic deformation. It is observed that plastic hinges occur at the corners of the cell walls when plastic deformation takes place, but the strains in other zones  (2008), as shown in Fig. 3b. Moreover, the shape of the compressed cells is similar to images of wood compressed transversely, as observed using light microscopy (Fig. 3c). This similarity to some extent validates the modelling approach. The shape of the modelled wood cells when compressed to a densified form is shown in Fig. 4a. During the compression of untreated wood, all the cell walls and lumens collapse. The densified state of the present model is very close to the actual densified wood microstructure as shown in Fig. 4b. This indicates that a reasonable simulation of the compression and densification process of untreated wood is achieved.

Nominal stress-strain relation
To analyse the mechanical characteristics of the microstructural model during compression, nominal stress σ nom and nominal strain ε nom are defined as: where F is the contact force between the model and the upper or bottom rigid plate during compression, A is the initial contact area between them, d y is the vertical displacement specified on the upper rigid plate, and H is the initial height of the model. The relationship between the nominal stress and strain is discussed in this section.

Influence of yield stress of cell wall
The typical transverse compression nominal stress-nominal strain curve of the untreated wood of model 1 is shown in Fig. 5a (Huang et al. 2020), in which the dotted lines show the method used to extract modulus of elasticity (E), yield stress and strain for the onset of yield (σ y and ε y ), and the onset of densification. The location of σ y and ε y corresponds to the yield point of untreated wood, and after the plastic region the curve represents the onset of densification. It is noted that the macro yield stress defined by the nominal stress-strain curve is the nominal stress. The influence of the yield stress of the cell wall, σ 0 (see equations S1 and S2, Supplementary Information), in the nominal stress-strain relation is discussed in this section. Comparing the typical compression curve for wood (Fig. 5a), and the simulated results (Fig. 5b), they appear similar in the sense that both curves increase during the elastic region and then change slope at the yield point. However, the FE model is only representative of the complete structure of wood, and only contains twelve cells, so the curve appears more sensitive to local buckling failures during deformation. After the plastic deformation region, both curves enter the densification region and present an upward trend where the wood changes to a densified state.
In Fig. 6, the bulk stress-strain relationships of six models with different yield stress values and a cell wall thickness of 2.82 μm are compared with the experimental results for white spruce, aspen and white ash. It is observed that yield stress values of 30 MPa, 40 MPa and 50 MPa result in nominal yield stresses of 4 MPa, 5 MPa and 6 MPa, respectively. The nominal compressive stress-strain curves are similar to experimental data for native softwood such as white spruce earlywood (Tabarsa and Chui 2000). The transverse compressive strength of softwoods has been reported  (Song et al. 2018) to be in the range 3-6 MPa (Tabarsa and Chui 2001). The numerical results with a yield stress of 50 MPa agree better with the experimental values of aspen which is characterised as a soft hardwood. The numerical nominal stress-strain results of the models with a yield stress value in the range of 200-300 MPa correlate better with experimental values of white ash as reported by Aimene and Nairn (2015). A yield stress of 100 MPa results in a nominal transverse compressive strength of 10 MPa that is representative of a hardwood species. For hardwoods, the nominal yield stress of European oak has been reported as 10.4 MPa (Ljungdahl et al. 2006) and white ash as 10.9 MPa Chui 2001, 2000).
It should be noted that the cell geometry, size, and density in the numerical simulations are representative of a softwood structure and firm conclusions for hardwood   Tabarsa and Chui (2001) species cannot be derived. Moreover, the numerical results are based on chemically treated cell walls and the comparison with the bulk compressive yield strength of native wood specimens is indicative. It is unknown to what extent the chemical treatment can affect the yield stress of the cell wall and if the additional hydrogen bonding reported by Song et al. (2018) takes place after or before the densification process. An increased degree of hydrogen bonding will increase the yield stress in the cell walls. To the best of the authors' knowledge, there is lack of experimental stress-strain data for chemically treated wood under transverse compression in the literature which limits firmer conclusions. The experimental nominal compressive stress-strain plots in Fig. 6 show a plateau at the yield point, whereas small decreases after the yield point are observed in the numerical results. This is attributed to the micro-scale of the numerical model that cannot fully represent the results of macroscopic wood compression, and it is more sensitive to local buckling of the cell walls.
Comparing the three curves of the numerical model with yield stresses of 100 MPa, 200 MPa and 300 MPa as shown in Fig. 6, in the initial stage of compression (compression displacement is from 0 to 5 µm-linear elastic range), the three curves are close to each other exhibiting the same elastic compressive stiffness because the elastic parameters of the three models are the same. The slopes of the curves for the models with 100, 200 and 300 MPa yield stresses are 263, 266 and 274 MPa, respectively. The slope in the initial linear elastic range of the experimental results for white ash yields a modulus of 260 MPa which is similar to the numerical findings. A higher compliance is expected in chemically treated wood and the initial elastic slope of treated wood should be smaller than that of untreated wood. This will be further discussed in the sect. Stress-strain relations of untreated and chemically treated wood. The model with a yield stress of 100 MPa is the first to undergo plastic deformation. Comparing the nominal stress values of the three curves, it shows that the higher the yield stress of the wood cell walls, the higher the nominal yield stress during the compression and densification process. Comparing the three curves of the models with yield stresses of 200 MPa and 300 MPa, it was found that their features are similar, and yet undergo volatile changes after the initial yielding. However, if the yield stress of the cell wall is 100 MPa, its stress-strain curve changes smoothly, and there are no sudden drops during compression after the first yielding, which is attributed to more elements being plastically deformed and may lead to less damage of the cells during the manufacture of the densified wood. Figure 7 shows the compression results of wood under three yield stresses. From the finite element results, it is found that when the yield stress is 100 MPa, 13% of the elements have equivalent plastic strain, when the yield stress is 200 MPa, 15% of the elements have equivalent plastic strain, and 16% of the elements have equivalent plastic strain when the yield stress is 300 MPa. It shows that the larger the material yield stress, the more elements undergo plastic deformation during the compression process.
The influence of yield stress on the wood cell wall plasticisation and deformation was also studied. Earlywood, following simulated chemical treatment, was selected to study the nominal stress-strain relations in the cases with yield stresses of 100 MPa, 200 MPa and 300 MPa (Models 1, 2 and 3). The densified states of the models with different yield stresses obtained by FE simulation are shown in Fig. 7. In the densified wood compression experiments by Song et al. (2018), untreated wood was compressed to 20% of the thickness. For the numerical models with 100 MPa, 200 MPa and 300 MPa yield stress, the models are all compressed to 21% of their original thickness, so close to the experimental approach by Song et al. (2018).

Stress-strain relations of untreated and chemically treated wood
The cell wall is rich in lignin and hemicellulose rendering the densification process of its natural state difficult. To make the cell wall easier to collapse during the manufacture of densified wood, chemical reagents are used to treat the material to remove part of the lignin and hemicelluloses.
In the sect. Material model and parameters, the elastic parameters of the wood after chemical treatment were reported. We note that we do not actually model either chemical or heat treatment of wood, and only the effects of the changes in cellulose, hemicellulose and lignin content on mechanical properties are considered in this study. To compare the influence of the chemical treatment of the wood on the compression process, the stress-strain curves of the untreated wood before and after simulated chemical treatment are shown in Fig. 8 (Models 1 and 4). It is noted that the elastic parameters only affect the compression stress and stiffness below 5% strain. The nominal compressive yield stress and stiffness of the untreated wood is higher than that of the wood after chemical treatment. However, when the material reaches the yield stress, the subsequent plastic deformation of the two models is almost equal, and the abrupt changes in load during the compression process take place at similar strain values irrespective of the presence or not of chemical treatment. It is noted that yield stress is the main influencing factor of the densification process under large deformations.

Mechanical behaviour during compression
Focusing on the wood after chemical treatment, and the sample with a 2.82 µm cell wall thickness, with a 100 MPa yield stress (Model 1), the mechanical behaviour of the wood during compression was analysed step-by-step. The deformation states corresponding to the typical 9 points as shown in Fig. 8 are shown in Fig. 9. State 1 is the initial state of the model. From state 1 to 3, the slope of the stress-strain curve shows an upward trend. During this phase, the model undergoes elastic deformation, but in state 2, some elements of the cell walls undergo plastic deformation. When the strain reaches 5% (state 3), the stress-strain curve changed from an upward to a downward trajectory. This indicates that the model undergoes yielding in stage 3, and large deformation occurs up to a pronounced increase in the compression load. The yield point of the curve is 10 MPa, but the maximum stress of the model reaches 40 MPa at complete densification. Therefore, the yield stress of the cell wall and the nominal yield stress of the whole model (wood) are different, i.e., the cell wall reaches the yield limit before the cell wall collapses.
The nominal stress-strain curve from state 3 to state 4 shows a downward trajectory. During this process, the cells of both the upper (row 1) and bottom row (row 3) incline to the right side. Due to the buckling of the cell walls, the bottom surface of the model "sinks" and the reaction force of the bottom surface decreases, so the curve shows a downward trajectory. With an increase in the deformation, the stress-strain curve shows an upward trajectory from state 4 to state 5. During this process, the larger deformation of cells in the bottom row takes place, and the inner walls (S3 layer) of those cells undergo self-contact, which leads to an increase in the stress until state 5 is reached. The process from state 5 to state 6 and the process from state 6 to state 7 are similar, with further consolidation of the structure. State 6 reflects the collapse of cells in row 1, indicating that the inner cell wall undergoes Fig. 8 Nominal stress-strain curve for treated and untreated wood (yield stress 100 MPa, cell wall thickness 2.82 μm). Annotated with different stages (1-9) of deformation self-contact, and state 7 reflects the collapse of cells in row 2. It is worthwhile noting that in the initial state when the model begins to compress, the areas where buckling deformation and cell wall collapse first occur are in cells 9 to 12 (see Fig. 1a). So, the lower cells first reach the yield stress and undergo plastic deformation. The stress-strain curve at state 8 shows the densification process, and state 9 is the final state of the densified wood.

Influence of cell wall thickness on stress-strain relation
Using the same elastoplastic material parameters, wood microstructural models with different cell wall thicknesses were used to study the influence of this parameter on the stress-strain behaviour. These data are shown in Fig. 10. In the simulation of the compression process by ABAQUS, the automatic increment time step is used. When the model is compressed to a very dense state and can hardly compress any further, the time step deceases to very small increments approaching zero. The simulation then stops. The simulation of the model with smaller wall thickness stops Comparing the stress-strain curves of the two models, it is found that the changing trends and turning points of the two curves are almost the same, and the cell walls collapse at similar strain values. Yet, the smaller increases in compression load after the first yield, defined as states 4-5 and 6-7 in the sect. Mechanical behaviour during compression, seem to occur at smaller strain values in the model with the highest thickness. This is attributed to the greater self-contact in the thicker cell walls at a given displacement due to the resulting decrease in the cell lumen size. The model using a larger cell wall thickness has larger stresses associated with it due to the higher geometric stiffness. In addition, when the curve has not reached the yield stress, the slope of the nominal stress-strain curve is also affected by the thickness of the cell wall. It can be seen from Fig. 10, the thicker the cell wall the higher the linear elastic slope of the curve. The bulk yield stress is higher in the case of thick cell walls and is attained at higher strain values, indicating that higher energy and longer times are required during the compression process. This correlates well with SEM studies (Felhofer et al. 2020) of the densification process in Scots pine (Pinus sylvestris L.) that showed that the deformation starts first in earlywood due to the weaker thinner cell walls.

Stress distribution in cell wall during compression
In the FE model, the cell wall is regarded as a laminate composed of S1, S2, S3 and P layers. During compression of the model, the stresses of the layers are 1 3 discontinuous due to material mismatch of the layers. The distributions of the stress in the fibre direction of each layer in cell 10 at typical states during compression are shown in Fig. 11, from which it is seen that the stress distributions of different cell layers are obviously different due to the differences of material parameters and fibre orientations of the layers. At the beginning of compression, the distributions of stress in the fibre direction first occur within the vertical cell wall areas and stress concentrations are observed close to the chamfered areas. Then, the high stress transmits to the cell wall in the vertical direction and buckling occurs when the material reaches the yield stress (elasto-plastic buckling). The early plasticisation of the cell corners alters the local boundary conditions of the cell walls facilitating buckling failure modes. As the cells collapse, the stress gradually increases in all parts of the cell wall, and finally, the stress concentrates on the upper surface cell wall due to the increased contact surface area. It is assumed here that the material properties remain constant within the cell wall and the same yield strength applies in all cell wall layers. However, a recent study by Felhofer et al. (2020) showed that there is a rearrangement of wood polymer components during lateral compression with highly lignified regions occurring in the middle lamella, at the cell corners, and in compressed regions. It was postulated that lignin interacts with hemicellulose via electrostatic bonding and can therefore easily migrate upon mechanical loading. The compressed regions were stiff regions of densely packed cellulose microfibrils and with high lignin accumulation. This seems to be a natural reinforcing mechanism of wood during lateral compression to minimise stress concentrations at the cell corners as indicated in Fig. 11.
Initially, the stress first increases in the P layer, the outermost material of the cell wall, and then gradually transfers to the S3 layer. Although the stress of the P layer increases first, the S3 layer reaches the stress yield limit first because the thickness of the S3 layer is too small and the geometric stiffness is the smallest. Then, the materials of the P layer and S1 layer reach a higher stress during the buckling deformation process. The thickness of the S2 layer material is the largest, so the higher stress occurs in the cell wall in the collapsed state. At the same moment of the compression process, the S3 layer exhibits higher S11 stresses, and the S2 layer relatively smaller stresses. The thicker the layer, the smaller the stress. During the compression process, the S3 layer is compressed first, followed by the S1 and P layers. The thickest layer S2 layer supports the cell structure during the compression process preventing it from collapse. Therefore, the S2 layer clearly plays a very important role in the cell wall during compression.

Conclusion
Simplified microstructural FE models of wood before and after chemical treatment have been developed to simulate the compression and densification process of the material. The nominal bulk compressive stress-strain relationships of the models during the compression have been obtained, and the deformation modes and stress variations during the compression process have been analysed. The deformation process and nominal stress-strain relationships of the models during the compression and densification process obtained by the FE simulation are consistent with those observed by experiments and microscopy studies; this validates the use of the FE modelling method. The presented simulation method could therefore be used to model the densification manufacturing process of different species for different chemical treatments and with a range of lignin, hemicellulose, cellulose contents and microfibril angles. A higher cell wall yield stress results in increased bulk compressive strengths during the densification process and higher energy requirements in the manufacturing process. It is proven, through numerical simulation, that the yield stress is the main factor affecting the stress-nominal strain curve. The thicker the cell wall, the larger the geometric stiffness, so a larger compression load is required to deform the material. Therefore, during the densification process earlywood cell walls deform first as has been observed with microscopy in the literature. By comparing nominal bulk compressive stress-strain curves of wood before and after chemical treatment it is found that the degradation of the elastic modulus properties only affects the compression stress up to 5% strain. The subsequent plastic deformation of the two models after the first yielding is similar. However, the numerical simulation did not consider the change of yield stress due to resulting reductions in lignin and hemicellulose content. It is expected that chemical treatment affects the cell wall yield stresses, and more research is needed to derive firm conclusions on the stress-strain response of chemically modified wood during the densification process. The stress distribution in each layer of an indicative cell was analysed with increasing displacement. Stress concentrations were first observed at the cell wall corners that led to yielding and elasto-plastic buckling of the cell walls.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.