Multi-Scale Analysis and Testing of Tensile Behavior in Polymers with Randomly Oriented and Agglomerated Cellulose Nanofibers

Abstract

Cellulose nanofiber (CNF) has been accepted as a valid nanofiller that can improve the mechanical properties of composite materials by mechanical and chemical methods. The purpose of this work is to numerically and experimentally evaluate the mechanical behavior of CNF-reinforced polymer composites under tensile loading. Finite element analysis (FEA) was conducted using a model for the representative volume element of CNF/epoxy composites to determine the effective Young’s modulus and the stress state within the composites. The possible random orientation of the CNFs was considered in the finite element model. Tensile tests were also conducted on the CNF/epoxy composites to identify the effect of CNFs on their tensile behavior. The numerical findings were then correlated with the test results. The present randomly oriented CNF/epoxy composite model provides a means for exploring the property interactions across different length scales.

1. Introduction

Cellulose nanofiber (CNF) has been accepted as a valid nanofiller that can improve the mechanical properties of composite materials by both mechanical and chemical methods [1,2,3,4]. Nishino et al. [5] determined that the elastic modulus of natural cellulose in the longitudinal direction was 138 GPa. Blefdzki et al. [6] proposed that the ultimate tensile strength of cellulose is estimated to be 17.8 GPa, which is seven times higher than that of steel. Sain and Oksman [7] studied the mechanical properties of CNF to design CNF-reinforced composite materials. Chirayil et al. [8] and Ansari et al. [9] showed that the mechanical properties of CNF-dispersed polymer composites increased compared with those of pure polymer. However, it was reported that an agglomeration of CNF in the matrix, poor adhesion, and water uptake decreased the mechanical properties. Kurita et al. [10] investigated the flexural properties of the CNF-dispersed epoxy layer inserted into woven glass fiber-reinforced polymer (GFRP) composite laminates. They indicated that the flexural strength increases, whereas the flexural modulus is almost constant when inserting two CNF/epoxy interface layers. Xie et al. [11] prepared extremely low CNF-reinforced epoxy matrix composites and tried to understand the strengthening mechanism of the CNF/epoxy composite through a three-point flexural test and finite element analysis (FEA). They proposed the possibility that CNF behaves as an auxiliary agent to enhance the structure of epoxy molecules and not as a reinforcing fiber in the epoxy matrix.

Prediction of the stress and strain requires a detailed FEA. The finite element method is regarded as an important tool in the design and analysis of composite materials, and some finite element works on the modeling of composite materials have been presented by our group [12,13,14].

In this paper, the tensile behavior of CNF-reinforced polymer composites is investigated. A finite element model that considered the randomly oriented CNF in the polymer was developed, and FEA was performed for the composites with randomly oriented CNFs. Tensile tests were then carried out on the CNF/epoxy composites. The stress–strain curve and Young’s modulus for the CNF-reinforced polymer composites were then discussed in detail.

2. Numerical Procedure

2.1. Definition of Representative Volume Element

In this study, ANSYS Multiphysics code (ver. 11.0) was used to perform a three-dimensional FEA. Figure 1 shows the representative volume element (RVE) of the CNF-reinforced polymer composites. As shown in Figure 1a, the RVE was created such that the CNF was embedded in the center of the polymer matrix. The Cartesian coordinate system o-$xyz$ was defined such that the $x$-axis was in the longitudinal direction, and the $y$- and $z$-axes were in the transverse plane of the RVE. $L^{\mathrm{f}}$ and $L_{}^{\mathrm{m}}$ represent the length of the CNF and matrix, respectively. The superscript $\mathrm{f}$ represents the fiber, and $\mathrm{m}$ denotes the matrix. $d_{}^{\mathrm{f}}$ is the diameter of the CNF. $W_{}^{\mathrm{m}}$ and $T_{}^{\mathrm{m}}$ are the width and thickness of the matrix, respectively. Here, $T_{}^{\mathrm{m}}$ was set to be equal to $W_{}^{\mathrm{m}}$ so that the cross-section $yz$-plane is square. The CNF aspect ratio $L_{}^{\mathrm{f}}$/$d_{}^{\mathrm{f}}$ was assumed to be equal to the matrix aspect ratio $L_{}^{\mathrm{m}}$/$W_{}^{\mathrm{m}}$ [15], and the relationship between $W_{}^{\mathrm{m}}$ and CNF volume fraction $V_{}^{\mathrm{f}}$ was given by the following equation:

$$V^{\mathrm{f}}=\frac{\pi }{4}{\left(\frac{d^{\mathrm{f}}}{W^{\mathrm{m}}}\right)}^3$$

(1)

The CNF and matrix were assumed to have isotropic elastic properties, and the constitutive equation can be written as

$$\left(\begin{array}{c}{\epsilon }_{xx}^{\delta }\\ {\epsilon }_{yy}^{\delta }\\ {\epsilon }_{zz}^{\delta }\\ 2{\epsilon }_{yz}^{\delta }\\ 2{\epsilon }_{zx}^{\delta }\\ 2{\epsilon }_{xy}^{\delta }\end{array}\right)=\left(\begin{array}{cccccc}1/E^{\delta }& -{\nu }^{\delta }/E^{\delta }& -{\nu }^{\delta }/E^{\delta }& 0& 0& 0\\ -{\nu }^{\delta }/E^{\delta }& 1/E^{\delta }& -{\nu }^{\delta }/E^{\delta }& 0& 0& 0\\ -{\nu }^{\delta }/E^{\delta }& -{\nu }^{\delta }/E^{\delta }& 1/E^{\delta }& 0& 0& 0\\ 0& 0& 0& 1/G^{\delta }& 0& 0\\ 0& 0& 0& 0& 1/G^{\delta }& 0\\ 0& 0& 0& 0& 0& 1/G^{\delta }\end{array}\right)\left(\begin{array}{c}{\sigma }_{xx}^{\delta }\\ {\sigma }_{yy}^{\delta }\\ {\sigma }_{zz}^{\delta }\\ {\sigma }_{yz}^{\delta }\\ {\sigma }_{zx}^{\delta }\\ {\sigma }_{xy}^{\delta }\end{array}\right)(\delta =\mathrm{m},\mathrm{f})$$

(2)

where ${\epsilon }_{ij}^{\delta }$ and ${\sigma }_{ij}^{\delta }$ are the strain and stress components ($i, j=x, y,z$), $E_{}^{\delta }$ is the Young’s modulus, $G_{}^{\delta }$ is the shear modulus, and ${\nu }_{}^{\delta }$ is the Poisson’s ratio.

2.2. Finite Element Analysis of Representative Volume Element

The RVE of CNF-reinforced polymer composites can be assumed to be transversely isotropic with the plane of isotropy being the $yz$-plane (see Figure 1). The solution obtained from the analysis of the RVE is equivalent to the solution of unidirectional, uniformly dispersed, CNF-reinforced polymer composites. So, for the transversely isotropic RVE, there are five independent elastic properties: longitudinal Young’s modulus $E_{\mathrm{l}}^{\mathrm{r}}$, transverse Young’s modulus $E_{\mathrm{t}}^{\mathrm{r}}$, longitudinal Poisson’s ratio ${\nu }_{\mathrm{lt}}^{\mathrm{r}}$, transverse Poisson’s ratio ${\nu }_{\mathrm{tt}}^{\mathrm{r}}$, and longitudinal shear modulus $G_{\mathrm{lt}}^{\mathrm{r}}=G_{\mathrm{tl}}^{\mathrm{r}}$. The superscript $\mathrm{r}$ denotes the RVE. In addition, transverse shear modulus $G_{\mathrm{tt}}^{\mathrm{r}}$ is calculated by $G_{\mathrm{tt}}^{\mathrm{r}}=E_{\mathrm{t}}^{\mathrm{r} }/2\left(1+{\nu }_{\mathrm{tt}}^{\mathrm{r}}\right)$. Then, the longitudinal Young’s modulus $E_{\mathrm{l}}^{\mathrm{r}}$ and Poisson’s ratio ${\nu }_{\mathrm{lt}}^{\mathrm{r}}$, transverse Young’s modulus $E_{\mathrm{t}}^{\mathrm{r}}$ and Poisson’s ratio ${\nu }_{\mathrm{tt}}^{\mathrm{r}}$, and longitudinal shear modulus $G_{\mathrm{lt}}^{\mathrm{r}}=G_{\mathrm{tl}}^{\mathrm{r}}$ can be obtained from the FEA of the RVE under longitudinal normal, transverse normal, and longitudinal shear loadings, respectively.

The appropriate constraints on the RVE under various loads have been determined by Sun and Vaidya [16]. For longitudinal and transverse normal loading, one-eighth of the RVE ($0\le x\le L_{}^{\mathrm{m}}/2,0\le y\le W_{}^{\mathrm{m}}/2,0\le z\le W_{}^{\mathrm{m}}/2$) was considered, as shown in Figure 1b. The displacement boundary conditions for the RVE under normal load in the longitudinal $x$-direction are [17]

$$\begin{array}{c}{u}_x^{\delta }\left(0,y,z\right)=0\\ u_x^{\mathrm{m}}\left(L^{\mathrm{m}}/2,y,z\right)=u_x^{*}\end{array}\}\hspace{1em}0\le y\le W^{\mathrm{m}}/2,\hspace{1em}0\le z\le W^{\mathrm{m}}/2\hspace{1em}(\delta =\mathrm{m},\mathrm{f})$$

(3)

$$\begin{array}{c}{u}_y^{\delta }\left(x,0,z\right)=0\\ u_y^{\mathrm{m}}\left(x,W^{\mathrm{m}}/2,z\right)=u_y^0\end{array}\}\begin{array}{cc}& 0\le x\le L^{\mathrm{m}}/2,\begin{array}{cc}& 0\le z\le W^{\mathrm{m}}/2\end{array}\end{array}(\delta =\mathrm{m},\mathrm{f})$$

(4)

$$\begin{array}{c}{u}_z^{\delta }\left(x,y,0\right)=0\\ u_z^{\mathrm{m}}\left(x,y,W^{\mathrm{m}}/2\right)=u_z^0\end{array}\}\begin{array}{cc}& 0\le x\le L^{\mathrm{m}}/2,\begin{array}{cc}& 0\le y\le W^{\mathrm{m}}/2\end{array}\end{array}(\delta =\mathrm{m},\mathrm{f})$$

(5)

where $u_x^{\delta }, u_y^{\delta }$, and $u_z^{\delta }$ are the displacement components in the $x$-, $y$-, and z- directions, respectively; $u_x^{*}$ is the applied uniform displacement in the $x$-direction; and $u_y^0$ and $u_z^0$ are the uniform displacements in the $y$- and $z$-directions. The uniform displacements in the $y$- and $z$-directions are determined from the condition that resultant stresses at y = W^m/2 plane and z = W^m/2 plane are zero, respectively. The longitudinal Young’s modulus $E_{\mathrm{l}}^{\mathrm{r}}=E_x^{\mathrm{r}}$ and Poisson’s ratio ${\nu }_{\mathrm{lt}}^{\mathrm{r}}={\nu }_{xy}^{\mathrm{r}}={\nu }_{xz}^{\mathrm{r}}$ are given by the following equations [16]:

$$E_x^{\mathrm{r}}=\frac{{\sigma }_{xx}^{*}}{u_x^{*}/\left(L^{\mathrm{m}}/2\right)}$$

(6)

$${\nu }_{xy}^{\mathrm{r}}=-\frac{u_y^0/(W^{\mathrm{m}}/2)}{u_x^{*}/\left(L^{\mathrm{m}}/2\right)}$$

(7)

where ${\sigma }_{xx}^{*}$ is the mechanical mean stress on the $x=L_{}^{\mathrm{m}}/2$ plane of the RVE. The mechanical mean stress ${\sigma }_{xx}^{*}$ is obtained as follows [16,17]

$${\sigma }_{xx}^{*}=\frac{{{\displaystyle \int }}_0^{W^{\mathrm{m}}/2}{{\displaystyle \int }}_0^{W^{\mathrm{m}}/2}{\sigma }_{xx}^{\mathrm{m}}\left(L^{\mathrm{m}}/2,y,z\right)\mathrm{d}y\mathrm{d}z}{{(W^{\mathrm{m}}/2)}^2}$$

(8)

Transverse normal loading can be simulated by applying the uniform displacement in the $y$- or $z$-direction. An analysis procedure like that used for the case of longitudinal normal loading can be employed for the case of transverse normal loading, and the transverse Young’s modulus $E_{\mathrm{t}}^{\mathrm{r}}=E_y^{\mathrm{r}}=E_z^{\mathrm{r}}$ and Poisson’s ratio ${\nu }_{\mathrm{tt}}^{\mathrm{r}}={\nu }_{yz}^{\mathrm{r}}$ = ${\nu }_{zy}^{\mathrm{r}}$ can be obtained.

For longitudinal shear loading, a whole RVE $(-L_{}^{\mathrm{m}}/2\le x\le L_{}^{\mathrm{m}}/2,-W_{}^{\mathrm{m}}/2\le y\le W_{}^{\mathrm{m}}/2,-W_{}^{\mathrm{m}}/2\le z\le W_{}^{\mathrm{m}}/2$) was considered. The displacement boundary conditions for the RVE under longitudinal shear load are [17]

$$\begin{array}{c}{u}_x^{\mathrm{m}}\left(x,y,-W^{\mathrm{m}}/2\right)=0\\ u_y^{\mathrm{m}}\left(x,y,-W^{\mathrm{m}}/2\right)=0\\ u_z^{\mathrm{m}}\left(x,y,-W^{\mathrm{m}}/2\right)=0\\ u_x^{\mathrm{m}}\left(x,y,W^{\mathrm{m}}/2\right)=u_x^{*}\\ u_y^{\mathrm{m}}\left(x,y,W^{\mathrm{m}}/2\right)=0\\ u_z^{\mathrm{m}}\left(x,y,W^{\mathrm{m}}/2\right)=0\end{array}\}\begin{array}{cc}& -L^{\mathrm{m}}/2\le x\le L^{\mathrm{m}}/2,\begin{array}{cc}& -W^{\mathrm{m}}/2\le y\le W^{\mathrm{m}}/2\end{array}\end{array}$$

(9)

In addition, the points on the $x=-L_{}^{\mathrm{m}}/2$ and $x=L_{}^{\mathrm{m}}/2$ planes with the same $y$- and $z$-coordinates need to be constrained to have the same displacements in all three directions as follows [16]:

$$\begin{array}{c}{u}_x^{\mathrm{m}}\left(-L^{\mathrm{m}}/2,y,z\right)=u_x^{\mathrm{m}}\left(L^{\mathrm{m}}/2,y,z\right)\\ u_y^{\mathrm{m}}\left(-L^{\mathrm{m}}/2,y,z\right)=u_y^{\mathrm{m}}\left(L^{\mathrm{m}}/2,y,z\right)\\ u_z^{\mathrm{m}}\left(-L^{\mathrm{m}}/2,y,z\right)=u_z^{\mathrm{m}}\left(L^{\mathrm{m}}/2,y,z\right)\end{array}\}\begin{array}{cc}& -W^{\mathrm{m}}/2\le y\le W^{\mathrm{m}}/2,\begin{array}{cc}& -W^{\mathrm{m}}/2\le z\le W^{\mathrm{m}}/2\end{array}\end{array}$$

(10)

The longitudinal shear modulus $G_{\mathrm{lt}}^{\mathrm{r} }=G_{\mathrm{tl}}^{\mathrm{r}}=G_{zx}^{\mathrm{r}}$ is given by [16]

$$G_{zx}^{\mathrm{r}}=\frac{{\sigma }_{zx}^{*}}{u_x^{*}/W^{\mathrm{m}}}$$

(11)

The mechanical mean stress ${\sigma }_{zx}^{*}$ is obtained by [16,17]

$${\sigma }_{zx}^{*}=\frac{{{\displaystyle \int }}_{-W^{\mathrm{m}}/2}^{W^{\mathrm{m}}/2}{{\displaystyle \int }}_{-L^{\mathrm{m}}/2}^{L^{\mathrm{m}}/2}{\sigma }_{zx}^{\mathrm{m}}\left(x,y,W^{\mathrm{m}}/2\right)\mathrm{d}x\mathrm{d}y}{L^{\mathrm{m}}{W}^{\mathrm{m}}}$$

(12)

In the analysis, the interface between the CNF and matrix was assumed to be perfectly bonding. The eight-node, three-dimensional solid elements were employed to mesh the RVE. The RVE model under longitudinal and transverse normal loadings consisted of 62,736 nodes and 1400 elements, while the model for longitudinal shear loading had 229,869 nodes and 55,680 elements.

2.3. Definition of Randomly Oriented CNF-reinforced Composite

Figure 2a shows the image of the polymer composite with randomly oriented CNF in the FEA. The global Cartesian coordinate system O-$XYZ$ and the local Cartesian coordinate systems o-$xyz$ for the composite are shown. To create a model of randomly oriented CNF-reinforced polymer composite, the mechanical properties of the RVE denoted by the superscript r were applied to each element, and all local coordinate systems were rotated randomly for the global coordinate system, as shown in Figure 2b, although the geometries of the element and unit cell are different. This random rotation of each element coordinate system was conducted based on the following equation:

$$\left(\begin{array}{c}x\\ y\\ z\end{array}\right)=\left(\begin{array}{ccc}\cos \phi & 0& \sin \phi \\ 0& 1& 0\\ -\sin \phi & 0& \cos \phi \end{array}\right)\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \omega & -\sin \omega \\ 0& \sin \omega & \cos \omega \end{array}\right)\left(\begin{array}{ccc}\cos \theta & -\sin \theta & 0\\ \sin \theta & \cos \theta & 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right)$$

(13)

where $\left(x_0, y_0, z_0\right)$ represents the initial axes before rotation. The initial orientation of the local coordinate systems was the same as that of the global coordinate system, and $\omega ,\phi ,$ and $\theta$ are the Euler angles indicating the orientation of the local coordinates ($x, y, z$) with respect to the global coordinates $\left(X, Y, Z\right)$. Furthermore, these were generated from random numbers for each element. By using the above definition, it could be possible to simulate randomly oriented CNF-reinforced polymer composites.

The polymer composite with randomly oriented transversely isotropic RVEs (randomly oriented CNFs) is isotropic material [18]. Hence, the constitutive equation for the composite can be written in the following form:

$$\left(\begin{array}{c}{\epsilon }_{XX}^{\mathrm{c}}\\ {\epsilon }_{YY}^{\mathrm{c}}\\ {\epsilon }_{ZZ}^{\mathrm{c}}\\ 2{\epsilon }_{YZ}^{\mathrm{c}}\\ 2{\epsilon }_{ZX}^{\mathrm{c}}\\ 2{\epsilon }_{XY}^{\mathrm{c}}\end{array}\right)=\left(\begin{array}{cccccc}1/E^{\mathrm{c}}& -{\nu }^{\mathrm{c}}/E^{\mathrm{c}}& -{\nu }^{\mathrm{c}}/E^{\mathrm{c}}& 0& 0& 0\\ -{\nu }^{\mathrm{c}}/E^{\mathrm{c}}& 1/E^{\mathrm{c}}& -{\nu }^{\mathrm{c}}/E^{\mathrm{c}}& 0& 0& 0\\ -{\nu }^{\mathrm{c}}/E^{\mathrm{c}}& -{\nu }^{\mathrm{c}}/E^{\mathrm{c}}& 1/E^{\mathrm{c}}& 0& 0& 0\\ 0& 0& 0& 1/G^{\mathrm{c}}& 0& 0\\ 0& 0& 0& 0& 1/G^{\mathrm{c}}& 0\\ 0& 0& 0& 0& 0& 1/G^{\mathrm{c}}\end{array}\right)\left(\begin{array}{c}{\sigma }_{XX}^{\mathrm{c}}\\ {\sigma }_{YY}^{\mathrm{c}}\\ {\sigma }_{ZZ}^{\mathrm{c}}\\ {\sigma }_{YZ}^{\mathrm{c}}\\ {\sigma }_{ZX}^{\mathrm{c}}\\ {\sigma }_{XY}^{\mathrm{c}}\end{array}\right)$$

(14)

where the superscript $\mathrm{c}$ represents the composite. The mechanical properties of the composites denoted by the superscript c were obtained using the mechanical properties of the RVE denoted by the superscript r.

2.4. Finite Element Analysis of Randomly Oriented CNF-reinforced Composite Specimen

We performed the FEA by using a one-eighth model of a JIS K 7164 tensile specimen, as shown in Figure 3. The global Cartesian coordinate system O-$XYZ$ was employed, and the local coordinate system o-$xyz$ was defined (not shown here) and rotated randomly. $L^{\mathrm{c}}$ = 215 mm represents the total length; $L_1^{\mathrm{c}}$ = 60 mm and $L_2^{\mathrm{c}}$ = 44.3 mm are the lengths of the reduced section and the grip section, respectively; $W_1^{\mathrm{c}}$ = 45 mm and $W_2^{\mathrm{c}}$ = 25 mm are the widths of the reduced section and the grip section, respectively; $T^{\mathrm{c}}$ = 2 mm is the total thickness; and $r^{\mathrm{c}}$ = 60 mm represents the curvature radius.

The boundary conditions are represented by the following equations:

$$u_X^{\mathrm{c}}\left(0, Y,Z\right)=0$$

(15)

$$u_Z^{\mathrm{c}}\left(X, Y,0\right)=0$$

(16)

$$u_X^{\mathrm{c}}\left(L^{\mathrm{c}}/2, Y,Z\right)=u_X^{**}$$

(17)

where $u_X^{\mathrm{c}}$, $u_Y^{\mathrm{c}}, u_Z^{\mathrm{c}}$ represent the displacement components in $X$-, $Y$-, and $Z$- directions, and $u_X^{**}$ denotes the applied displacement.

For the finite element model, the eight-node, three-dimensional solid elements were employed. The model consisted of 3141 nodes and 6800 elements.

From the FEA, the stress and strain of the tensile specimens were obtained by averaging these values of the nodes in the reduced section surface ($0\le X\le L_1^{\mathrm{c}}/2,0\le Y\le W_1^{\mathrm{c}}/2,Z=T^{\mathrm{c}}/2$) because the calculated solutions were different for each node due to the model of randomly oriented CNF-reinforced polymer composites. Additionally, to simulate the randomness of the solutions, an FEA was performed for each three times under the same conditions. Finally, the Young’s modulus $E_{}^{\mathrm{c}}$ obtained from the FEA was compared with the experimental results.

3. Experimental Procedure

In this work, nanocomposites consisting of an epoxy matrix filled with dry CNFs were considered. The diameters of the dry CNFs ranged from 25 to 50 nm, whereas their lengths were between 5 and 50 μm. The processing concept for CNF/epoxy composite is summarized in Figure 4. The fabricated composite samples contained CNF concentrations of 0.73, 2.2, and 3.7 vol %. Neat epoxy samples were also prepared.

Tensile tests were conducted on the CNF/epoxy composite specimen following JIS K 7164, and the stress–strain curves were measured. The stress was determined by dividing the load by the cross-sectional area of the narrow section. Test specimen dimensions are shown in Figure 3. At least five tests were performed under each condition.

4. Results and Discussion

4.1. Parameters of Representative Volume Element

The Young’s modulus of the CNF is $E^{\mathrm{f}}=$ 138 GPa [7]. The Poisson’s ratio ${\nu }_{}^{\mathrm{f}}$ was assumed to be equal to that of carbon nanotube and hence ${\nu }_{}^{\mathrm{f}}=0.2$ [17]. On the other hand, the Young’s modulus $E^{\mathrm{m}}=$ 3.23 GPa and Poisson’s ratio ${\nu }_{}^{\mathrm{m}}=0.34$ of the neat epoxy were obtained from the tensile tests.

Figure 5 gives the scanning electron microscopy (SEM) images of the surface for the 3.7 vol % and of fracture surface for the 0.73 vol % CNF/epoxy composite specimen. It was found that the agglomerated CNF clusters were almost uniformly dispersed in the epoxy matrix. The fracture of CNF clusters was observed on the fracture surface of CNF/epoxy composites without disentwining. This result indicates that the applied load was transferred to CNF clusters and that CNF clusters contributed to enhance the epoxy resin matrix during the tensile test. The CNF in the epoxy matrix tended to be agglomerated, and the apparent aspect ratio became small. From this image, the aspect ratio $L_{}^{\mathrm{f}}$/$d_{}^{\mathrm{f}}$ = $L_{}^{\mathrm{m}}$/$W_{}^{\mathrm{m}}$ = 7.5 was assumed, taking $L_{}^{\mathrm{f}}$ and $d_{}^{\mathrm{f}}$ as the length and diameter of the agglomerated CNF cluster. For an aspect ratio above 50, the predicted Young’s modulus is almost independent of the aspect ratio [17]. For simplicity here, the agglomerated CNF clusters were assumed to be isotropic and to have the same Young’s modulus and Poisson’s ratio as the CNF.

4.2. CNF/epoxy Composite Analysis

Figure 6 shows the measured and calculated stresses as a function of the measured strain for the neat epoxy and 2.2 vol % CNF/epoxy composite. Dotted lines indicate the measured data. Dashed lines are the calculated normal stress in the X-direction at a chosen point (X = Y = 0 mm and Z = T^c/2 = 1 mm here) from a linear elastic FEA. Each measured curve exhibits the linear initial region, and the calculated curves agree with the measured data very well. In the later loading stage, nonlinearity of the curve exists. This is due to the material nonlinearity. In order to determine the nonlinear stress–strain response numerically, the uniaxial stress–strain behavior, for example, based on Ref. [18] must be represented. More details can be found in our previous works on elastic-plastic finite element modeling [14,19], which is out of scope of the current work. Figure 6 also shows the calculated result for the 2.2 vol % CNF/epoxy composite from the multi-scale FEA (see solid line). It is interesting to note that the stress obtained from the multi-scale FEA was larger than that obtained from the linear FEA. It is anticipated that unexpected high stress is generated in the composite due to the inhomogeneity of the material. The predicted tensile stress ${\sigma }_{XX}^{\mathrm{c}}$ distributions obtained from the (a) linear FEA and (b) multi-scale FEA of the randomly oriented CNF/epoxy composite with the CNF volume fraction $V_{}^{\mathrm{f}}$ of 2.2 vol % and the aspect ratio $L^{\mathrm{f}}/d^{\mathrm{f}}$ = 7.5 under normal strain of 0.015 are shown in Figure 7. From Figure 7a, the uniform stress distribution was generated at the reduced section from the linear FEA. On the other hand, from Figure 7b, it was found that the non-uniform stress distribution was generated at the reduced section from the multi-scale FEA. The results indicate that if the CNFs were randomly dispersed and completely bonded to the epoxy matrix, the stress of the CNF/epoxy composite would be higher. In other words, the linear analysis may underestimate the stress.

Figure 8 shows the measured and predicted Young’s modulus for the neat epoxy, the 2.2 vol % CNF/epoxy composite, and the 3.7 vol % CNF/epoxy composite. The randomness of the solutions was confirmed for each CNF volume fraction in the composites. For the neat epoxy, the predicted Young’s modulus agreed with the measured data. On the other hand, the predicted Young’s moduli for the CNF/epoxy composites were larger than the measured data. Additionally, although the results are not shown here, when a CNF aspect ratio of 50 was used, the predicted Young’s modulus for the 3.7 vol % CNF/epoxy composite became approximately $E_{}^{\mathrm{c}}$ = 4.65 GPa. The predicted Young’s modulus $E_{}^{\mathrm{c}}$ tends to increase as the CNF volume fraction and aspect ratio increase. This implies that the dispersion of the agglomerated CNF cluster with a high spectral ratio, even in randomly oriented CNF/epoxy composites, increases the Young’s modulus.

Table 1. Comparison of the Young’s modulus of epoxy matrix nanocomposites.

Nanomaterial	Content	Modulus (GPa)	Increase Rate (%)
Dry CNF	0	3.23		This study
	0.73 vol %	3.30	2.1
	2.2 vol %	3.54	10
	3.7 vol %	3.85	19
Silica	0	2.96		[20]
	2.5 vol %	3.20	8.1
	4.9 vol %	3.42	16
	7.1 vol %	3.57	21
	9.6 vol %	3.60	22
	13.4 vol %	3.85	30
CNT	0	2.90		[21]
	0.1 wt.%	3.01	3.8
	0.2 wt.%	3.11	7.2
	0.3 wt.%	3.26	12
Aligned CNT	0	2.5		[24]
	4.5 vol %	19	660
	8.4 vol %	32	1180
	21.4 vol %	50	1900
	0	2.69		[22]
Graphene (GP)	4.0 wt.%	2.89	7.4
Surface-modified GP	4.0 wt.%	3.27	22
Clay	0	3.53		[23]
	2.0 wt.%	3.58	1.4
	5.0 wt.%	3.66	3.7
CNF slurry	0	2.20 *		[11]
	0	2.20 *	26
	0.25 vol %	3.67 *	67
	0.74 vol %	2.84 *	29

* Flexural modulus.

presents a comparison of the measured Young’s modulus of the present CNF/epoxy nanocomposite with other nanocomposites from the literature. Although dry CNFs were agglomerated in the epoxy matrix, these CNFs seem to tend to increase the Young’s modulus of the epoxy matrix as much as they do other nanomaterials [20,21,22,23], except for aligned multi-walled CNT [24]. On the other hand, the flexural modulus was significantly increased by the extremely low CNF slurry addition [11]. This increment cannot be explained by the theories for discontinuous fiber-reinforced composites. Therefore, it seems that the CNF slurry behaves not as a reinforcing fiber in the epoxy matrix but as an auxiliary agent to enhance the structure of the epoxy molecule.

Figure 9 shows the distribution of predicted tensile stress ${\sigma }_{xx}^{\delta }$ at y = 0 plane for the RVE of the composite with the CNF volume fraction V^f of 3.7 vol % and the aspect ratio $L^{\mathrm{f}}/d^{\mathrm{f}}$ = 7.5 and 50 under the applied tensile stress of 1 MPa. The stress gradients were confirmed, and the stress increased and became constant toward the center of the CNF. It was found that the stress at the center of the RVE for $L^{\mathrm{f}}/d^{\mathrm{f}}$ = 50 was larger than that of $L^{\mathrm{f}}/d^{\mathrm{f}}$ = 7.5.

There are many methods to estimate the overall mechanical and physical properties of composites [17,23,24]. In the present work, we assume the perfect bonding between the CNF and epoxy matrix because of its simplicity and use the three-dimensional finite element method. In fact, a weak boundary exists (see Figure 4). Additional modeling is required for a full interpretation of the mechanical and physical properties of CNF/epoxy composites. The main achievements of this research include the development of a finite element model that can be easily extended to complex problems (e.g., CNF-reinforced polymer composites with interfacial layers, CNF-reinforced polymer composites with slip at interfaces).

5. Conclusions

In this paper, we studied the tensile behavior of CNF/epoxy composites by performing multi-scale FEA and tensile tests. It was shown that by increasing CNF content, the Young’s modulus increases. It was also found that agglomeration of CNF decreases the aspect ratio of CNF and decreases the Young’s modulus. A correlation, obtained between the numerical and experimental results for the Young’s moduli of the CNF/epoxy composites under tension, revealed that unexpected high stress is generated in CNF/epoxy composites due to the inhomogeneity, and the inhomogeneity tends to increase the Young’s modulus. The results from this study offer some basic insights into the reinforcing mechanisms in CNF/epoxy composites under tension. These will enable the design of multifunctional nanocomposite materials based on CNF for advanced engineering applications.

All conclusions are based on the numerical predictions for CNF/epoxy composites with perfect bonding between the CNF and epoxy matrix. Some nanocomposites indicate the presence of weak boundaries. This is a challenging research area, and sooner or later, some progress will be made.