Modeling of modulus and strength in void-containing clay platelet/cellulose nanocomposites by unit cell approach

Abstract

Clay platelets/cellulose nanofibril nanocomposites are recyclable engineering materials of interest for sustainable development. There is substantial experimental data for mechanical properties, but modeling efforts are scarce. Here, a conceptual unit cell with voids was used in finite element modeling. Predictions for the in-plane modulus taking voids into account were more accurate than the classical rule of mixtures. Simulations also reveal that the cellulosic matrix undergoes shear and tensile deformation with inclined fracture located near the end of platelets for low clay content while predicting brittle tensile failure for high clay content. The results from the unit cell approach provide improved understanding of experimental observations, supporting the strive to better understand mechanisms of deformation and fracture.

Graphical Abstract

Keywords:

• Artificial nacre
• brick-and-mortar
• porosity
• FEM
• nanocellulose
• biocomposites

1. Introduction

Cellulose nanocomposites from cellulose nanofibrils (CNF) are of interest for sustainable development as they are from renewable resources and have good mechanical properties [1]. Unlike synthetic thermoset composites, cellulosic composites may show advantages for recycling and biodegradation in line with the immediate need for sustainable solutions [2]. Also, the in-plane modulus for CNF as a matrix is much higher than for an isotropic polymer. Transparent nacre-inspired nanocomposites are of particular interest. Such composites are normally thin films. CNF nanofibrils are mixed with clay platelets in water and filtered and dried. Capillary effects and strong physical interparticle interactions lead to the formation of a fairly dense and uniform CNF matrix with, e.g. montmorillonite (MTM) natural clay as a nanoscale reinforcement phase [3]. The MTM inclusions have platelet shape with a ‘diameter’ of roughly 100–400 nm and height of 1–15 nm (depending on agglomeration tendency). Platelets often show high in-plane orientation. This class of nanocomposites can show a brick-and-mortar microstructure at high clay content [4], and stiffness and strength can be tailored. Interfaces between the constituents are often weak as there are no primary chemical bonds; the only interactions are of secondary nature (hydrogen bonds, van der Waals interactions etc.). In composites, weak interfaces promote higher strain-to-failure and improved toughness as cracks propagate through the weak interfaces, and the overall toughness of high inorganic content composites is relatively high [5, 6].

The stiffness-strength-toughness balance of such composites is fundamentally important [7]. Modeling at multiple length-scales is interesting for nanocomposites, which can be atomistic [8, 9] or continuum mechanics [10, 11]. Here, we limit our discussion to the continuum-based approach and compare results to experimental test data of macroscopic properties (e.g. mechanical properties measured from tensile testing), which is suitable for materials design purposes. Composite modulus and strength for materials with in-plane oriented platelet reinforcements have been targeted for nacre-like composites [12–15]. For the in-plane composite modulus, Christensen [16], proposed a model with perfect platelet-matrix adhesion, which is a variant of the classical rule of mixtures. Also, shear-lag-based models were extensively explored [17–19]. The shear-lag models were also extended to predict the composite strength by considering interfacial failure controlled by the interfacial shear strength. However, all mentioned models are one-dimensional simplifications, and lack accurate details because the composite geometry leads to a very complex stress state.

Finite element (FE) simulations associated with a two- or three-dimensional unit cell (composite geometry) provide more detailed predictions. Gao and co-workers [20, 21] did FE simulations for biological nacre with protein matrix and calcium carbonate platelets and studied the effect of nano-inclusions on effective properties. Kim et al. [22] explored computational homogenization, using an FE-based model [10], to study printed platelet-staggered samples to validate analytical modeling and obtained good agreement with experimental measurements. For biocomposites, a complex three-dimensional unit cell generated from X-ray microcomputed tomography data was coupled to computational homogenization where elastic composite properties were well predicted [23]. All simulations essentially used the FE method (discretizing the geometry, solving the governing equations with boundary conditions by numerical integration), but what is specific is the unit cell. The proper composite geometry could also predict composite strength using FE damage models [24–26].

A unit cell is needed for CNF/MTM nanocomposites and the current ones for nacre composites do not account for porosity. Neat CNF films have some porosity [27] and in CNF/MTM the high MTM content compositions have significant nanoporosity [28, 29] since platelets in tactiods (stacks of MTM platelets) will have different lengths with irregular platelet edges. The discrepancy between modulus predictions from classical micromechanics (e.g. the rule of mixtures) based on no voids and experimental data is noticeable [28, 29].

Here, a representative two-dimensional unit cell with voids (void content experimentally verified) is proposed, to address the problem of lower mechanical property data, such as elastic modulus and strain-to-failure, for nanocomposites with high MTM content. The new porous unit cell was used in the FE framework to predict in-plane composite modulus, whose values were compared to the rule-of-mixture (no porosity) commonly used for the studied nanocomposites. Similarly, the composite strength was evaluated with the FE method, but using the phase-field damage model to promote matrix rupture in order to evaluate the strain-to-failure. Predictions are compared with experimental data for two material systems: CNF matrix with MTM platelets by Medina et al. [29] and TEMPO-oxidized cellulose nanofiber (TOCN) matrix with well-dispersed MTM platelets by Li et al. [30]. Pore content was characterized in both studies.

2. Material and methods

2.1. Micromechanical model geometry

The idealized material microstructure contains in-plane oriented MTM platelets in a homogeneous matrix representing random-in-the-plane cellulose nanofibers (CNF) fused into a continuous phase. Figure 1 (from ref. [30]) shows a nanostructural TEM micrograph of a cross-section of CNF/MTM with 10 vol% of MTM. In the context of polymer matrix nanocomposites, the micrograph confirms a comparably homogeneous platelet distribution.

Figure 1. Transmission electron microscopy image from the cross-section of a nanocomposite of CNF matrix with a volume fraction of 10% MTM platelets (Image by Dr. Ogawa, CERMAV, Grenoble).

For the composition in Figure 1, porosity was estimated from nanocomposite density and found to be less than 5%. Since porosity increases with MTM volume fraction, voids are related to MTM edges and agglomerated stacks of MTM (cf. Ref. [30]). This hypothesis is supported by the fact that horizontal platelets in Figure 1 have different lengths and are longer than expected. An agglomerated stack of several 1 nm thick MTM platelets will then have a ‘jagged’ edge with nanoscale pores in between.

The idealized two-dimensional geometry of clay in the present micromechanical composite model is depicted in Figure 2. The voids are positioned at the ends of the platelets. The blue dashed line indicates the representative unit cell of the composite. The small 1–2 coordinate systems represent the projection of the CNF matrix as if the 1-axis were properties along the nanofiber bundles and the 2-axis were transverse. Using the dimension of the unit cell, the volume fraction of the platelets VF_p can be determined as (1) $${\mathrm{VF}}_{\mathrm{p}}= \frac{h_{\mathrm{p}}{l}_{\mathrm{p}}}{(h_{\mathrm{p}}+{2h}_{\mathrm{m}})(l_{\mathrm{p}}+d_1)}=\frac{{h_{\mathrm{p}}}^2\rho }{(h_{\mathrm{p}}+{2h}_{\mathrm{m}})(l_{\mathrm{p}}+d_1)}$$(1) where $h_{\mathrm{p}}$ is the thickness of the platelet, $l_{\mathrm{p}}$ is the length of the platelet, $h_{\mathrm{m}}$ is half of the thickness of matrix sections between the overlapped platelets, and $d_1$ is the horizontal distance between the platelets. Equation (1)(1) $${\mathrm{VF}}_{\mathrm{p}}= \frac{h_{\mathrm{p}}{l}_{\mathrm{p}}}{(h_{\mathrm{p}}+{2h}_{\mathrm{m}})(l_{\mathrm{p}}+d_1)}=\frac{{h_{\mathrm{p}}}^2\rho }{(h_{\mathrm{p}}+{2h}_{\mathrm{m}})(l_{\mathrm{p}}+d_1)}$$(1) accounts for the aspect ratio of the platelet $\rho =l_{\mathrm{p}}/h_{\mathrm{p}}$ as well. In the same fashion, porosity (2) $$P= \frac{h_{\mathrm{p}}{d}_1}{(h_{\mathrm{p}}+{2h}_{\mathrm{m}})(l_{\mathrm{p}}+d_1)}$$(2) can also be expressed in terms of the dimensions in Figure 2.

Figure 2. Illustration of the geometry of the composite with its dimensions.

The geometry allows studying the effects from volume fraction of the platelet, the aspect ratio of the platelet and porosity on composite modulus. To reduce complexity, limitations related to imperfect interfaces between CNF and MTM platelets are not considered here although they would influence results in real materials.

2.2. Composite modulus: computational homogenization

For computational homogenization, the model is based on small deformations and the platelet-matrix and platelet-platelet bonds are assumed to be perfect. The composite stress (3) $$\overline{\mathit{\sigma }}=\frac{1}{V}{\int }_V\mathbf{\sigma }\mathrm{ }\mathrm{d}V$$(3) is calculated as the average stress in the volume of unit cell $V.$ The localization problem is given as follows: (4a) $$\nabla ·\mathbf{\sigma }=0\text{ in }V$$(4a) with the constitutive law as (4b) $$\mathbf{\sigma }=\mathbb{C}\left(\mathbf{x}\right)\mathbf{\epsilon } \mathrm{for}\mathrm{ }\forall \mathbf{x}\in V,$$(4b) where $\mathbb{C}(\mathbf{x})$ is the heterogeneous elasticity tensor as a function of the space variable $\mathbf{x}$ in microscopic scale. Equation (4)(4a) $$\nabla ·\mathbf{\sigma }=0\text{ in }V$$(4a) was resolved via the finite element method with Comsol Multiphysics solver along with the displacement continuity (periodic boundary conditions) [10] (5) $${\mathbf{u}}_{\mathbf{dst}}\text{-}\mathrm{ }{\mathbf{u}}_{\mathbf{src}}=\mathrm{ }\overline{\mathit{\epsilon }}\left({\mathbf{x}}_{\mathbf{dst}}\text{-}\mathrm{ }{\mathbf{x}}_{\mathbf{src}}\right) \text{ for }\forall \mathbf{x}\in V.$$(5)

In Comsol Multiphysics, this is implemented using Boundary Pair subnode to the outer boundaries of the unit cell. The boundaries source (src) and destination (dst) are shown in Figure 3a. The composite strain $\overline{\mathit{\epsilon }}$ was prescribed throughout the unit cell in Eq. (5)(5) $${\mathbf{u}}_{\mathbf{dst}}\text{-}\mathrm{ }{\mathbf{u}}_{\mathbf{src}}=\mathrm{ }\overline{\mathit{\epsilon }}\left({\mathbf{x}}_{\mathbf{dst}}\text{-}\mathrm{ }{\mathbf{x}}_{\mathbf{src}}\right) \text{ for }\forall \mathbf{x}\in V.$$(5) . In this study, the main interest in the composite modulus in the longitudinal direction $E_L$ and therefore the need in $\overline{\mathit{\epsilon }}$ is only for composite strain component ${\overline{\epsilon }}_L.$ The unit cell was discretized into small quadrilateral elements with quadratic interpolation functions. Both material components in the composite are considered linear elastic. The platelets are assumed to be orthotropic with elastic modulus in the longitudinal direction of 139 GPa, elastic modulus in the transverse direction of 174 GPa, and a Poisson’s ratio of 0.34 [31]. For the CNF matrix, the in-plane properties of the homogeneous matrix based on neat CNF film data are $E_1$= 15–17 GPa [29, 30] and ${\nu }_{12}=0.33$ [32], and the out-of-plane properties perpendicular to the fibrils are estimated as $E_2= {0.18E}_1$ [33]. The shear modulus $G_{12}$ is 2.5 GPa [34]. Once computed from Eqs. (3)(3) $$\overline{\mathit{\sigma }}=\frac{1}{V}{\int }_V\mathbf{\sigma }\mathrm{ }\mathrm{d}V$$(3) , (4), and (5)(5) $${\mathbf{u}}_{\mathbf{dst}}\text{-}\mathrm{ }{\mathbf{u}}_{\mathbf{src}}=\mathrm{ }\overline{\mathit{\epsilon }}\left({\mathbf{x}}_{\mathbf{dst}}\text{-}\mathrm{ }{\mathbf{x}}_{\mathbf{src}}\right) \text{ for }\forall \mathbf{x}\in V.$$(5) , the composite modulus in the longitudinal direction (horizontal, in-plane) is estimated by dividing the average stress in the unit cell by the constant composite strain $E_L={\overline{\sigma }}_L/{\overline{\epsilon }}_L.$ From classical micromechanics, the composite modulus in the longitudinal direction can be estimated by the rule of mixtures as [16] (6) $$E_L=\mathrm{ }{E}_{\mathrm{p}}{\mathrm{VF}}_{\mathrm{p}}+\mathrm{ }{E}_{\mathrm{m}}{\mathrm{VF}}_{\mathrm{m}}$$(6) where E is the elastic modulus, VF is the volume fraction and the subscripts p and m stand for platelet and matrix, respectively.

Figure 3. Geometry and boundary conditions for the modeling: (a) boundary source and destination for in-plane modulus study and (b) ¼ of the unit cell to estimate the composite strength based on cellulose matrix tensile fracture assumption.

2.3. Composite strength: phase-field damage and cohesive zone models

Strength criteria could be based on various failure mechanisms such as fracture of platelets, interfacial sliding and rupture of the matrix. Here, we account for interfacial sliding and assume that the composite fracture takes place when the matrix ruptures; if platelet fracture was the critical failure, composite strength would be much higher than experimental data. We investigate how two models integrated into the finite element framework can provide guidance when using the proposed unity cell. One way to study the contribution of the weak interface and pull-out mechanisms simultaneously is to numerically promote the matrix fracture by the phase field damage model in combination with cohesive zones at the interfaces for platelet sliding [24, 25]. The phase field damage model is a powerful tool to predict damage initiation and propagation in the matrix of which its variational form contains the fracture energy and a damage variable, which degrades the strain energy density [35]. Cohesive zone models (CZM) are well-established techniques that also describe crack initiation and propagation assuming traction forces in front of the crack tip. In a practical way, however, it is normally used along predefined boundaries instead. In composite material analysis, CZM are commonly used to model delamination mechanisms at interfaces.

Here, the nanofibril matrices are ductile by nanofibril network deformation mechanisms [36]. For this work, an elastic-plastic version of the phase field model was chosen (see detail in Appendix A) to account for plasticity before fracture. The platelet properties are the same as in previous simulations and it is assumed that the strength of the platelet is always higher than platelet stress. For the composite strength, the geometry is simplified as one-fourth of the unit cell as shown in Figure 3b. Rollers are applied on top and bottom to constrain the opening and promote interface sliding only. At the interface, a bilinear cohesive zone model (CZM) is then set [24], where the shear strength and fracture toughness are the main inputs. The shear strength is defined as 30 MPa, determined in Appendix B, and the interfacial fracture toughness is 1 J/m² predicted by in situ tests on MTM/cellulosic interfaces [37]. The calculations were implemented in COMSOL Multiphysics following the model implementation in ref. [38] and the variational form in [39].

2.4. Previous experimental results

This section revisits previous experiments to determine the dimensions of the unit cell. Both Medina et al. [29] and Li et al. [30] prepared their nanocomposites by mixing CNF and MTM suspensions, filtering through a microporous membrane and drying the mixture to form a nanocomposite film. Li et al. [30] improved sample preparation for better MTM dispersion and increased in-plane Young’s modulus. Figure 4a shows modulus as a function of MTM content for both materials. The moduli show a steady increase of up to ≈ 40 vol% of MTM content. Then, the modulus decreases with MTM content in Medina et al. [29], whereas Li et al. [30] report a stable modulus of ≈ 40 GPa also at higher MTM content. The drop in modulus in ref [29] is linked to increased porosity at higher MTM content, see Figure 4b. In ref [30], the MTM platelets show lower thickness, see Figure 4c (based on WAXD data), and this correlates with reduced porosity.

Figure 4. Experimental data: (a) in-plane composite modulus from tensile tests, (b) total nanocomposite porosity, (c) thickness of MTM platelets by Scherrer size from wide angle X-ray diffraction, data from Li et al. [30], Medina et al. [29].

Using the experimental data, the size of the unit cell can be determined as: (i) The lengths of the platelets measured by atomic force microscopy are $l_{\mathrm{p}}$= 100 nm for Medina et al. and $l_{\mathrm{p}}$= 400 nm for Li et al. [29, 30], (ii) thickness of the platelets $h_{\mathrm{p}}$ can be directly measured by Scherrer size from WAXD in Figure 4c, and (iii) the dimensions $d_1$ and $h_{\mathrm{m}}$ are estimated from Eqs. (1)(1) $${\mathrm{VF}}_{\mathrm{p}}= \frac{h_{\mathrm{p}}{l}_{\mathrm{p}}}{(h_{\mathrm{p}}+{2h}_{\mathrm{m}})(l_{\mathrm{p}}+d_1)}=\frac{{h_{\mathrm{p}}}^2\rho }{(h_{\mathrm{p}}+{2h}_{\mathrm{m}})(l_{\mathrm{p}}+d_1)}$$(1) and (2)(2) $$P= \frac{h_{\mathrm{p}}{d}_1}{(h_{\mathrm{p}}+{2h}_{\mathrm{m}})(l_{\mathrm{p}}+d_1)}$$(2) supported by data in Figure 4b and the given MTM platelet volume fraction [29, 30].

3. Results and discussion

3.1. Composite modulus: the volume fraction, platelet aspect ratio and porosity

A parametric study was carried out varying $h_{\mathrm{p}},$ $h_{\mathrm{m}}$ and $d_1$ while keeping $l_{\mathrm{p}}$ = 200 nm to investigate effects of the volume fraction of the platelets, aspect ratio of the platelet, and porosity on the in-plane composite modulus. Figure 5a shows the computed $E_{\mathrm{L}}$ as a function of VF_p for three configurations: $d_1$ = 40, 80 and 120 nm with a constant MTM aspect ratio of 50. The voids have negligible influence on $E_{\mathrm{L}}$ at very low platelet content, but a significant divergence in the $E_{\mathrm{L}}$ predictions can be seen at VF_p = 13 vol%. The $E_{\mathrm{L}}$ differ by 15% at VF_p = 13 vol% by changing from small voids ($d_1$=40) to large voids ($d_1$=120). The predictions would diverge even more for higher platelet content. The $E_{\mathrm{L}}$ versus aspect ratio was also studied in three configurations: $d_1$ = 40, 80 and 120 nm while $h_{\mathrm{m}}$ = 20 nm as shown in Figure 5b. In this case, porosity is driven in two ways, e.g. along with the platelet thickness in the abscissa, and through the three curves by $d_1.$ The same trend of low $E_{\mathrm{L}}$ for high porosity was found where the green curve ($d_1$ = 120 nm) exhibited lower $E_{\mathrm{L}}$ than the other two. In addition, for all curves, $E_{\mathrm{L}}$ increases to aspect ratio $\rho$ ≈ 15, and then steadily drops until $\rho$ = 50. This behaviour is due to the effect of platelet volume fraction via the parameter $h_{\mathrm{p}}$ on porosity via the parameter $d_1.$ The last graph in Figure 5c demonstrates how the porosity directly affects $E_{\mathrm{L}}$ for aspect ratios $\rho$ = 16.7, 25 and 50 while $h_{\mathrm{m}}$ = 20 nm. Here, porosity was increased by varying $d_1$ for all cases and the $E_{\mathrm{L}}$ decreases towards high porosity levels. The example illustrated in Figure 5 is not trivial to visualize because the VF_p, $\rho$ and P are simultaneously counteracting each other for $E_{\mathrm{L}}.$ Interestingly, the example shows how to optimize a system when the size of the platelets is controlled (e.g. identifying the maximum value of $E_{\mathrm{L}}$ at $\rho$ ≈ 15), which fits the materials design purpose in this study. Coupling the proposed geometry with FE-based computational homogenization turns out to be a useful material design tool for optimization studies such as those carried out by Barthelat [14], but adding porosity effects, as introduced in the unit cell.

Figure 5. Parametric study from the FE model showing how the composite modulus in the longitudinal direction is affected by the platelet volume fraction, platelet aspect ratio, and porosity. The length of the platelet is $l_{\mathrm{p}}$ = 200 nm, while varying $h_{\mathrm{p}},$ $h_{\mathrm{m}}$ and $d_1$ in different configurations. Note the info on fixed quantities in each graph.

3.2. Composite modulus: case studies

To investigate whether the results are reasonable, a new set of simulations were performed where the geometric parameters of the unit cell – the dimensions $h_{\mathrm{p}},$ $h_{\mathrm{m}},$ $d_1$ and $l_{\mathrm{p}}$ – were estimated from experiments [29, 30]. The size of the nanoplatelets is known ($h_{\mathrm{p}}$ and $l_{\mathrm{p}}$) from the Scherrer size and AFM measurements. The remaining quantities $h_{\mathrm{m}}$ and $d_1$ are calculated by the Eqs. (1)(1) $${\mathrm{VF}}_{\mathrm{p}}= \frac{h_{\mathrm{p}}{l}_{\mathrm{p}}}{(h_{\mathrm{p}}+{2h}_{\mathrm{m}})(l_{\mathrm{p}}+d_1)}=\frac{{h_{\mathrm{p}}}^2\rho }{(h_{\mathrm{p}}+{2h}_{\mathrm{m}})(l_{\mathrm{p}}+d_1)}$$(1) and (2)(2) $$P= \frac{h_{\mathrm{p}}{d}_1}{(h_{\mathrm{p}}+{2h}_{\mathrm{m}})(l_{\mathrm{p}}+d_1)}$$(2) given that the Porosity and volume fraction of the platelets are also known. The parameters in the three sets of simulations for each material system were deliberately chosen to match the experimental modulus drops for Medina et al. [29] while the modulus keeps increasing for Li et al. [30]. The input values for all dimensions are tabulated in Table 1.

Table 1. Input values for geometrical microstructure parameters (see Figure 2) on which the predictions in Figures 6 and 7 are based.

Materials	MTM platelets [Vol%]	$l_{\mathrm{p}}$ [nm]	$h_{\mathrm{p}}$ [nm]	$h_{\mathrm{m}}$ [nm]	$d_1$ [nm]	$P$ [Vol%]
	7	100	3.0	15	40	3
CNF/MTM (Medina et al.)	36	100	3.0	2.0	20	11
	46	100	5.0	2.0	20	17
	10	400	2.5	10	50	2
TOCN/MTM (Li et al.)	20	400	2.5	4	70	5
	40	400	2.5	1.5	50	9

Figure 6 shows the predictions of $E_{\mathrm{L}}$ versus VF_p, where the same trend is found between simulations (solid line) and experiments (dashed line). The predictions are slightly higher than the experimental data, but the values are comparable. As seen in the experiments, the modulus decreases in Medina et al. and it correlates with increased platelet thickness from 3 to 5 nm. From the porous material modeling perspective, this lowers the stress in the platelets (c.f. shear lag assumptions, e.g. Jackson et al. [17] and Kotha et al. [18]) and the lower stress in the platelets obviously leads to lower composite modulus when computing Eq. (3)(3) $$\overline{\mathit{\sigma }}=\frac{1}{V}{\int }_V\mathbf{\sigma }\mathrm{ }\mathrm{d}V$$(3) . The analysis is complex but the lower stress in platelets is primarily because of associated porosity, see Table 1. The aspect ratio (l_p/h_p) is also lowered, but this leads to a weaker effect. The other material system (Li et al. [28]) has a lower aspect ratio of the platelets but a very low porosity, which provides a higher composite modulus than for Medina et al. [29]. The modulus drop does not happen in the Li et al. TOCN/MTM system because the thickness of the platelets is constant (2.5 nm) for the three simulations. The graph also contains calculated values by the rule of mixtures from Eq. (6)(6) $$E_L=\mathrm{ }{E}_{\mathrm{p}}{\mathrm{VF}}_{\mathrm{p}}+\mathrm{ }{E}_{\mathrm{m}}{\mathrm{VF}}_{\mathrm{m}}$$(6) (in green) for the sake of comparison. Interestingly, for low volume fraction, the rule of mixtures predicts the composite modulus well, but it deviates for high platelet content. The rule of mixtures does not account for aspect ratio and porosity effects, which are important in this problem as shown in Figure 5; High platelet content leads to relatively thick tactoids and large voids.

Figure 6. In-plane composite modulus predictions as well as experimental values for two CNF/MTM material systems. The solid lines are the FE simulations, which follow the same trend as the experiments (dashed lines). the rule of mixtures in green overestimates the modulus for high MTM content.

3.3. Composite strength and toughness: case studies

Here, six simulations for strength analysis were performed where the size of the unit cells are listed in Table 1. The predicted stress-strain curves shown in Figure 7a,b were plotted so that the composite stress was computed by Eq. (3)(3) $$\overline{\mathit{\sigma }}=\frac{1}{V}{\int }_V\mathbf{\sigma }\mathrm{ }\mathrm{d}V$$(3) and the tensile strain by applied displacement over half length of the unit cell. The simulations were executed until complete rupture of the matrices (i.e. when the phase field damage variable reached a unit value). The values of composite stress for the two nanocomposites are comparable to experimental observations in [29, 30]. Strain-to-failure was decreased from matrix strain concentration effects as the MTM content was increased, in agreement with experimental observations [29, 30]. Figure 7c shows fields of the damage evolution in the matrices from a few load steps to complete matrix rupture for TOCN/MTM and CNF/MTM nanocomposites with lower MTM content. The fields reveal where the fracture occurs as well as crack orientation. For the CNF/MTM nanocomposite with 7 vol% MTM, the matrix rupture occurs near the ends of the platelets and is inclined, which indicates shearing in the matrix before fracture. For CNF/MTM nanocomposites with 10 vol% MTM and all other remaining examples, the matrix rupture occurred in the middle matrix segment region connecting the two platelets (c.f. dimension $d_1$ in Figure 2). The crack is straight and perpendicular to platelet orientation. Values of ‘matrix thickness’ $h_{\mathrm{m}}\mathrm{ }\ge$ 15 nm result in matrix failure due to shearing and $h_{\mathrm{m}}\mathrm{ }\le$ 10 nm result in matrix failure due to tension. Higher $h_{\mathrm{m}}$ also lead to higher strain-to-failure in the composites. Note that the CNF matrix is modeled as orthotropic elastic and isotropic plastic, allowing for matrix ductility.

Figure 7. Stress-strain curves (a–b) based on one-fourth of the unit cell and computed from the nanocomposite data in Table 1. The fields (c) illustrate the nature of matrix deformation and fracture during increased loading.

The dimension $h_{\mathrm{m}}$ is related to the volume fraction of MTM. High values of $h_{\mathrm{m}}$ means low MTM content in nanocomposites. High $h_{\mathrm{m}}$ allows nanofibril sliding, which promotes global plasticity and high strain-to-failure for the nanocomposites [36]. Lower $h_{\mathrm{m}}$ means high MTM content and higher local strain (and stress) in the CNF matrix at a given global composite strain. The strain to failure of the nanocomposite is thus reduced.

4. Concluding remarks

Modulus and stress-strain behavior of cellulose/clay (CNF/MTM) nanocomposites were modeled based on a two-dimensional unit cell with voids and imperfect cellulose/clay interfaces, in order to analyze reasons for lower nanocomposite mechanical properties found in previous publications. Predictions for in-plane modulus $E_{\mathrm{L}}$ for these nacre-mimetic nanocomposites follow the experimental trend for nanocomposites with different MTM contents. Since nanocomposites have imperfect structures, classical micromechanical models without voids, such as the rule of mixtures, cannot provide good estimate values, especially for high MTM content.

Stress-strain curves of the nanocomposites with different MTM contents were also predicted based on imperfect cellulose/clay interfaces (CZM) and a damage development of the cellulose matrices (the phase-field damage model), including cellulose matrix plasticity and voids. Stress and strain levels were well predicted including the decreased strain to failure with increased clay content. Local strain concentration effects and damage development were magnified with increased clay content, and the nature and location of matrix and interface damage could be identified.

Recent experimental work on CNF/MTM nanocomposites has shown that the extent of MTM agglomeration (tactoid formation, stacking of platelets) can be very limited [30], even at MTM content as high as 40 vol%. Yet modulus and strength are much lower than expected from micromechanics. The present investigation suggests that porosity, imperfect interfaces and strain concentration effects are important reasons. The CNF/MTM interface problem in particular needs further work, both in terms of materials design to improve nanocomposite properties but also modeling efforts based on nanostructure data. Since continuum mechanics models have limitations for nanocomposites, they need to be combined with other approaches, e.g. molecular dynamics analysis of interface problems.

Acknowledgements

The authors acknowledge financial support by Knut and Alice Wallenberg (KAW) foundation through the Wallenberg Wood Science Center. Discussions with Dr. Lengwan Li on his sample preparation and experimental results of MTM/TOCN nanocomposites are also acknowledged.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was funded by Knut och Alice Wallenbergs Stiftelse.

Notes on contributors

Marcus Vinícius Tavares da Costa

Marcus Vinícius Tavares da Costa received a PhD degree in Engineering Physics with specialization in Applied Mechanics from Uppsala University, Sweden, in 2020. After, he was a postdoctoral researcher at Royal Institute of Technology (KTH) investigating the strength-toughness balance of novel cellulose nanocomposites by exploring microstructures, weak interfaces, anisotropy, and lamella thickness. He is currently an assistant professor at Karlstad University, Sweden.

Lars A. Berglund

Lars A. Berglund is a professor in Wood composites at Royal Institute of Technology (KTH). He has been a visiting researcher at Stanford University, Cornell University and Kyoto University. His research focus is environmentally friendly composites based on nanocellulose and wood. He has published more than 300 journal papers, obtained around 10 patents, examined more than 30 PhD’s and is a member of the Royal Swedish Academy of EngineeringSciences.He justfinished anERCAdvancedGrant project on wood nanotechnologies. He has more than 35,000 citations and an H-index of 93 at Google Scholar.

Appendix A.

Material parameters for CNF and TOCN matrices

The chosen elastic-plastic version of phase field models that govern only the matrix domain $\mathrm{\Omega }$ is defined as [39] (A.1) $$\mathrm{\Pi }(\mathit{u},\mathrm{ }\varphi )={\int }_{\mathrm{\Omega }}{\left(1\text{-}\mathrm{ }\varphi \right)}^2\left(\frac{1}{2}\mathrm{ }{\mathit{\epsilon }}^{\mathrm{e}}\left(\mathit{u}\right):\mathit{\sigma }+\mathrm{ }\frac{1}{2}H{\gamma }^2+\mathrm{ }{\sigma }_{\mathrm{Y}}\gamma \right)d\mathrm{\Omega }+\mathrm{ }{\int }_{\mathrm{\Omega }}{G}_{\mathrm{c}}\left(\frac{{\varphi }^2}{2\mathcal{l}}\mathrm{ }+\mathrm{ }\frac{\mathcal{l}}{2}\nabla \varphi ·\mathrm{ }\nabla \varphi \right)d\mathrm{\Omega }\mathrm{ }$$(A.1) where $\varphi$ is the damage variable, ${\mathit{\epsilon }}^{\mathrm{e}}$ is the elastic strain tensor, $\mathit{\sigma }$ is the stress tensor, $H$ is the effective hardening modulus, ${\sigma }_{\mathrm{Y}}$ is the initial yield stress and $\gamma$ is equivalent plastic shear strain. When the damage variable $\varphi$ reaches a unit value, fracture happens. Otherwise, the material is intact. The second integral accounts for the fracture process [35, 40], where $G_{\mathrm{c}}$ is fracture energy and $\mathcal{l}$ is the length scale parameter that governs the size of the fracture zone. The properties $G_{\mathrm{c}}$ and $\mathcal{l}$ have to be calibrated. As mentioned in Section 2.3, this version of the phase-field damage model was chosen because the ductility of the CNF and TOCN matrices originated from the nanofibrils network. In particular, the sliding of nanofibrils has a major contribution (above 60%) to plasticity in cellulose microfibrils [36]. For this damage mechanism, the fracture toughness of sliding is 0.2 J/m² [8] and this value was therefore used in Eq. (A.1)(A.1) $$\mathrm{\Pi }(\mathit{u},\mathrm{ }\varphi )={\int }_{\mathrm{\Omega }}{\left(1\text{-}\mathrm{ }\varphi \right)}^2\left(\frac{1}{2}\mathrm{ }{\mathit{\epsilon }}^{\mathrm{e}}\left(\mathit{u}\right):\mathit{\sigma }+\mathrm{ }\frac{1}{2}H{\gamma }^2+\mathrm{ }{\sigma }_{\mathrm{Y}}\gamma \right)d\mathrm{\Omega }+\mathrm{ }{\int }_{\mathrm{\Omega }}{G}_{\mathrm{c}}\left(\frac{{\varphi }^2}{2\mathcal{l}}\mathrm{ }+\mathrm{ }\frac{\mathcal{l}}{2}\nabla \varphi ·\mathrm{ }\nabla \varphi \right)d\mathrm{\Omega }\mathrm{ }$$(A.1) throughout this work. The model in Eq. (A1) is essentially calibrated from the stress-strain curves on the neat matrices reported in the main references [29, 30]. The calibration is performed as the 2D plane stress approach in the same geometries used in tensile testing, e.g. stripes where the gauge lengths GL were set in the optical extensometers. The geometry as one-fourth of stripes is used to reduce the calculation time. Symmetry boundary conditions are set on the left and bottom while prescribed displacement on top as shown in Figure A1a. The phase-field damage model is then set in a small domain of the geometry to localize the fracture and avoid unnecessary calculation time. The in-build COMSOL plasticity material model as isotropic hardening is used for the remaining regions. Higher-order quadrilateral elements were used with mesh refinement in the phase field domain.

Figure A1. (a) Boundary conditions and material models and (b) FE predictions overlapped on the experimental results.

Figure A1b shows the calibrated stress-strain curves of each matrix overlapped with the experimental results. The curves were fitted well for both matrices where material parameters from Eq. (A1) are reported in Table A1. In phase field model calibration, the strain-to-failure was adjusted from the parameter $\mathcal{l}$ reported also in Table A1.

Table A1. Input values for prediction in Figure A1.

Materials	Width [mm]	GL* [mm]	Thickness [µm]	E [GPa]	H [GPa]	${\sigma }_{\mathrm{Y}}$ [MPa]	$G_{\mathrm{c}}$** [J/m^2]	$\mathcal{l}$ [nm]
TOCN (Li et al.)	5	25	30	18	4.2	120	0.2	1.4
CNF (Medina et al.)	4	10	40	12	3	115	0.2	1.1

*Gauge length from the references.

** From reference [Hou et al.].

Appendix B.

Composite strength: effect of the interfacial shear strength

The bilinear cohesive zone modeling needs the interfacial shear strength as an input. The sensibility of this material parameter is verified in this section on the stress-strain curves where composite in-plane stress is calculated by Eq. (3)(3) $$\overline{\mathit{\sigma }}=\frac{1}{V}{\int }_V\mathbf{\sigma }\mathrm{ }\mathrm{d}V$$(3) . Here, arbitrary dimensions were chosen as $h_{\mathrm{p}}$ = 4 nm, $h_{\mathrm{m}}$ = 20 nm, $d_1$= 40 nm and $l_{\mathrm{p}}$= 200 nm as an example. In addition, we assume the interfacial fracture toughness as 1 J/m² from previous work [37].

Figure B1 shows the curves from three sets of values: 5, 15 and 30 MPa. The first value of 5 MPa was measured for clay/cellulosic interfaces in a previous study [37]. The curve from this value provides a slope deviation after the initial curve (see the arrow), which is not seen in the experimental observations [29, 30]. In the numerical setting, this curve deviation is originated from platelet sliding and poor stress transfer. For the higher values of 15 to 30 MPa, the transition resembles the experimental data, which indicates, in turn, an unexpectedly stronger interfacial strength since the system has a chemically weak interaction at the nanofibrils/clay interface without a binder. During the composite assembly, the high aspect ratio platelets act as a nucleate site for the agglomeration of nanofibril bundles [41]. The nanofibrils placed and compacted on the rough surfaces of the platelets promote an interlocking mechanism during tension [42]. In this work, instead of explicitly setting the roughness of platelets in the geometry, we increase the interfacial shear strength up to 30 MPa in the cohesive zone modeling because adding such details in the modeling would increase dramatically the computation time.

Figure B1. Stress-strain curves from the one-fourth of the unit cell computed from three distinct values of interfacial shear strength.