Effect of interphase parameters on elastic modulus prediction for cellulose nanocrystal fiber reinforced polymer composite

Abstract In this study, the effective elastic modulus (Ec) of a cellulose nanocrystal fiber reinforced polymer composite was evaluated using the Mori-Tanaka and finite element (FE) model. The FE model was generated using a representative volume element with a periodic boundary condition. The mass fractions of the fiber in the composites (MFf) were set to 1, 2, and 3 wt.%. Elastic modulus values for interphase were input and were either uniform or exhibited a gradient. The E c for the uniform interphase region increased significantly with MFf, but was relatively low for interphase regions exhibiting a gradient. The results show that interphase parameters must be considered carefully when predicting Ec using the FE model.


Introduction
In recent years, interest in natural fibers with biodegradable and renewable characteristics has been increasing in accordance with demands for environmental preservation. In addition to their environmentally friendly characteristics, natural fibers are well suited for reinforcing polymer matrices because they possess superior stiffness relative to their weight. Among them, cellulose nanocrystals (CNCs) are natural fibers produced using the acid hydrolysis process by utilizing abundant natural resources such as wood. CNCs have a high tensile modulus, high surface area, and low density. As a result, they have attracted considerable interest. Manufacturing processes are being developed, production quantities are increasing, and numer-ous studies are being conducted [1][2][3][4][5]. This study assessed the mechanical properties of CNC composites. More specifically, the effective elastic modulus (Ec) of a CNC fiber reinforced polymer composite was evaluated using the Mori-Tanaka model and simulations using the finite element (FE) method.
Fundamentally, a CNC composite is a fiber-reinforced composite (FRC) with nano-size inclusions. A laboratory experimental study of a nano FRC is more time consuming and that of a conventional micro FRC. Furthermore, it is challenging to establish the accurate properties of an FRC because of their heterogeneous and anisotropic characteristics. Over the past several decades, analytical and numerical approaches have been developed for measuring the mechanical responses of FRCs while efficiently replacing experimental methods. The most commonly utilized analytical approaches include the Halpin-Tasi [6] model and Mori-Tanaka [7] model based on the Eshelby [8] solution. Compared to the numerical approach, a common feature of the analytical approach is that they can predict the effective elastic stiffness of a two-phase composite with inclusions and a matrix. Even though it is possible to simply assess the elastic modulus using a formalized formula, the limitation of this formula is that it cannot consider the interactions between adjacent inclusions. In other words, it is difficult to consider the stresses among the components inside a CNC composite using an analytical approach [9].
Numerical approaches can be used to overcome the above drawbacks of analytical approaches. Typical numerical approaches include molecular dynamics (MD) and FE simulation. First, the MD method is capable of assessing the potential energy and stress between inclusions and a matrix by calculating the equations of motion amongst atoms and molecules. However, as the MD method is atomistic, considerable time is required to generate a model with numerous embedded inclusions, indicating that it is ideal for the analysis of composites consisting of a single inclusion. Second, the FE method, which is based on continuum mechanics, is capable of performing simulation by freely controlling the parameters, including the as-pect ratio, volume/mass fraction, concentration, and orientation. Furthermore, by utilizing a representative volume element (RVE) with the periodic boundary condition (PBC) technique, with which an entire composite can be replaced with a unit cell, the FE method can overcome the computational limits caused by a large number of inclusion/matrix systems. Thus, the combination of the FE method and RVE can generate a numerical model that represents a realistic environment regarding a CNC composite. In other words, the model is capable of accurately analyzing the effective stiffness and stress contour plot between adjacent fibers within a short period [10][11][12]. Liu et al. [13] predicted progressive failure properties of carbon fiber reinforced epoxy composites by using the FE method and RVE. The results they obtained was compared with theories and experiments.
The interphase region that connects a matrix and a fiber in an FRC strongly affects the elastic properties of the composite [14,16]. However, it is difficult to precisely predict the properties of the interphase region. Thus, most studies have used the average values of the properties of the interphase region in simulations. Furthermore, there have been very few numerical studies of CNC composites. Therefore, it is necessary to establish a fundamental process for CNC composites and assess the resultant data produced using various interphase parameters. The establishment and assessment of such a process will not only be helpful for the investigation of CNC composites and the interphase region, but will also provide an FE method that can be applied to other composites with different types of embedded fibers.
This study was conducted to derive the effective elastic modulus (Ec) using an FE model for CNC fiber-reinforced polymer composites. The FE model was generated based on a three-dimensional (3D) RVE with the PBC technique.
The Ec was derived by controlling the parameters of the interphase region within the RVE. The final parameters were used to divide the interphase region into multiple gradients. Then, the derived Ec was compared with the values obtained from the Mori-Tanaka model and an experimental study. The stress interaction of components was investigated using the FE results. The novelty of this study is that it enables the application of a range of properties to the interphase region and describes a technique for generating a 3D RVE model based on this interphase region.

RVE generation
The RVE consists of three elements: a CNC fiber, an interphase region, and a polymer matrix ( Figure 1). The Digimat FE (version 2017, MSC) software was used because the PBC, which will be discussed later, was applied to the surface of the RVE. A periodic model is a unit cell model in which the components cut by the RVE's surface are symmetrical to the components on each of the two parallel planes (Figure 2). The mass fractions of the fiber in the composite (MF f ) were set to 1, 2, and 3 wt.%; as the generation of a fiber depends on its volume fraction (VF f ), VF f values of 0.7, 1.42, and 2.14 vol.% were applied, respectively. These values were determined by substituting the value of density (ρ) listed in Table 1 into the equation: where subscripts f and m represent the fiber and matrix, respectively.  The directional tensors of the fibers were arranged randomly in a 3D manner using an orientation distribution function (ODF). The function, Ψ(p), within the ODF includes angles θ and φ, which represent the fiber's vector, as shown in Figure 3. Vector p is described as follows: The second-order tensor in the fiber can be described as follows: A random arrangement is achieved when a 11 = a 22 = a 33 = 1 3 is satisfied in the above equation.

Geometry and material properties of components
Every component was assumed to be isotropic and homogeneous. The material properties of every component are listed in Table 1. The geometry and properties of the CNC fiber, including its elastic modulus (E f ), were taken from the literature [1], and an ellipsoid fiber with a width of 0.5 nm and a length of 5 nm (giving an aspect ratio of 10) was modeled. The elastic modulus of the matrix (Em) was taken from the value obtained from the tensile test of Polyamide 6 in a previous experimental study [4]. Numerous studies have been conducted on the parameters of the interphase region. Among them, two studies [4,15] measured the interval length and stiffness between a CNC fiber and the matrix by using atomic force microscopy (AFM). The results indicated that the radius of the interphase region was close to the width of the CNC fiber, while the stiffness of the interphase region increased nonlinearly as it transitioned from the matrix to the fiber. Based on the results obtained from an MD simulation, Sun et al. [16] applied the average value of the elastic modulus to the exponential function between E f and Em. In the present study, an interphase region with an offset of 0.5 nm from the fiber's surface was generated using the Solid-Works (2017, Dassault system) software. The elastic modulus of the interphase region (E i ) was categorized into two cases (subscript irepresents the interphase region). The first was the uniform interphase region, which had a constant modulus. The values of E i of the uniform interphase region were assigned as 1/2, 1/4, 1/8, and 1/16 of the value of E f -Em. The second was the interphase region exhibiting a gradient, for which the gradient of the equidistant  Figure 4 and the equation: where xa and x b are the minimum and maximum distances from the matrix to the corresponding element, respectively. As the Poisson's ratio of the interphase region was not available and had no significant influence on the elastic modulus of the composite [17], an intermediate value was assumed between the fiber and matrix.

Mesh and periodic boundary conditions
The generated RVE solid model was imported into the Hy-perMesh (version 2017, Altair) software and meshed. Assuming that the RVE was a continuum material, each component was perfectly bonded. As a preprocessing step for the PBC, the 2D mesh for the surface of the RVE was projected as an identical shape between parallel surfaces. A deviation in the mesh size, ranging from a minimum of 0.1 nm to a maximum of 1 nm, was applied to consider the surface curvature of the components ( Figure 5). The element type was a linear tetrahedron (C3D4), and the number of elements per RVE model ranged from approximately 4,000,000-5,000,000. The meshed model was imported into ABAQUS (version 2017, Dassault Systems) for FE analysis, which was the last step. The RVE can function as a periodic unit cell when the PBC is imposed upon its surface. The constraint equation was applied to every node of the surface of the RVE.  The basic concept of the equation is as follows: where X, Y, and Z represent the displacements corresponding to each axis, and + and − represent the two sides of the orthogonal to each axis ( Figure 6). Furthermore, the assigned deformation value (δ) in the above equation adjusts the strain on the X axis (ϵ 11 ). This implies that δis imposed upon every node on the surface of X + and is unified to 1 nm. It means that the surface of X + is shifted by 1 nm in the X + direction ( Figure 6). The effective elastic modulus in the X + direction (E 11 ) based on ϵ 11 and Hooke's law is determined using the equation:
The size of the RVE was determined through a validation test. Cubes with five different sizes ranging from 60 to 140 nm were designated, and nine models of each size were generated for validation. Then, the effective elastic modulus of the composite (Ec) was derived by calculating the average of the nine models. A 1 wt.% MF f with no interphase region was selected for simplicity.

Mori-Tanaka model
Among the various analytical approaches, the Halpin-Tsai [6] model has been widely used to measure the effective stiffness of composites. However, the Mori-Tanaka model produces more reasonable results when measuring composites containing inclusions with a large aspect ratio [18,19]. Therefore, this study measured the value of Ec using the Mori-Tanaka model.
The Mori-Tanaka model starts with a two-phase composite with a number of inclusions (I), all which have the same material, shape, and orientation, embedded within a matrix. The relationship between the volume fractions of the two phases is V 0 + V 1 = 1, where subscripts 0 and 1 represent the matrix and inclusion phases, respectively. The relationship between the volume averages of the strain field within the RVE is as follows: where ω without the subscript is the average value of the entire RVE. Any homogenization model can be described by the strain concentration tensor, B ϵ : The volume average of the strain over all the inclusions is related to tensor B ϵ . The volume average of the strain over the entire RVE is related to tensor A ϵ , with A ϵ being calculated as: The above model is calculated as a closed form by using the Eshelby [8] solution. The strain field inside an ellipsoid inclusion is uniform and related to the macrostrain as follows: where H ϵ is the strain concentration tensor for a single inclusion and is defined as follows: Bevensite et al. [20] suggested a useful interpretation of the Mori-Tanaka model. Each inclusion (I) in the actual RVE behaves as if it were isolated in the actual matrix. Finally, the macro effective stiffness, C, in the Mori-Tanaka model can be expressed as follows: The values of E i and VF i were set such that they were identical to the uniform interphase conditions of the FE model. Figure 7 presents the graph obtained for the RVE size validation. The Ec became stable once the RVE edge length reached 100 nm. This is because the number of embedded fibers decreases with the size of the RVE. This can result in tensors having fibers with a certain orientation instead of random orientations. As a result, for each variation in the parameters, nine RVE models with an edge length of 100 nm were generated to calculate nine values of E 11 , which were then averaged to calculate Ec. As the PBC was applied to control the nodes on the surface of the RVE, all the parallel surfaces, such as those shown in Figure 8, exhibited an identical stress and deformation shape. The obtained Ec data are presented in Table 2. First, in the no-interphase section shown in Figure 9a, the Mori-Tanaka and FE models were in good agreement in terms of Ec, with an error rate of 0.15%. However, the error rate increased considerably with E i . Based on this, we concluded that a high E i is the one of the key reasons for the difference between the Mori-Tanaka and FE models. As MF f in-    creased, the Ec of the FE model increased rapidly, while the Ec of the Mori-Tanaka model increased steadily. It can be speculated that the reason for this is high stress. As shown in Figure 10a, 10c, and 10e, which depict von Mises stress contour plots of the FE model with a uniform interphase region, a high stress zone between adjacent fibers was observed. As there is a geometric effect, an increase in MF f implies an increase in the number of adjacent fibers and a wider zone of high stress. All the models had a constant ϵ 11 , therefore, the development of high stress resulted in an increase in the Ec of the FE model. As the Mori-Tanaka model cannot reflect the above effect, Ec increased steadily. As shown in Figure 9b, the FE model with the interphase region exhibiting a gradient provided a noticeably different result compared with the FE model with the uniform interphase region and the Mori-Tanaka model. As E i changed from a linear function to a higher-degree function, the reduction rate of Ec fell with MF f . This implies that the slope of Ec/MF f would be lower than that of the uniform interphase region. This result is similar to that of the experimental study. A comparison of Figure 10b, 10d, and 10f with Figure 10a, 10c, and 10e, respectively, shows that the high stress zone decreased. As a result, it was determined that a lower slope of Ec/MF f would be obtained using the model with the interphase region exhibiting a gradient.

Results and discussion
From the above discussion, the following two facts are apparent: -when E i increases from the interphase region adjacent to the matrix to the interphase region adjacent to the fiber, the increase in Ec with MF f is smaller than in the case with the uniform E i ; -when E i increases nonlinearly from the interphase region adjacent to the matrix to the interphase region adjacent to the fiber, the increase in Ec with MF f is smaller than in the case with a linearly increasing E i .
The AFM is able to scan and measure the stiffness of the nanoscale surface in a mechanical way by using the cantilever. Referring to the AFM profile from the matrix to the fiber in the literature [4,15], the elastic modulus of interphase gradually increased (like the nonlinear curves in Figure 4), not linearly. Since the mesh in the nonlinear interphase model was divided into several layers, it was possible to apply the gradient value to the each mesh layers in the interphase region. Therefore, it is predicted for the nonlinear interphase in the FE model to produce the most consistent result. Among all the models, the interphase model with the second-degree function produced the smallest difference from the experimental study ( Figure 9). If the degree of the function is expanded to decimal places, the difference will decrease even further.
In this study, the properties of the FRC in mesoscale was analyzed by using the RVE and FE method. To predict the stress behavior of a local interphase region in detail, the MD simulation should be conducted further. As experimental methods for analyzing the local region, AFM and nanoindentation tests are recommended. In future research, the following issues should be resolved to derive a more reliable result: -as the actual experimental specimens were produced using an injection-molding machine, the fibers were likely to be oriented in a specific direction instead of randomly. An accurate tensor must be identified and applied to the fiber orientation; -the potential occurrence of voids and the agglomeration of fibers was neglected; -as CNC is an anisotropic material, its transverse modulus is different from its longitudinal modulus. This study only considered the longitudinal modulus.

Conclusions
This study evaluated changes in the Ec of a CNC fiber reinforced polymer by adjusting the parameters of the interphase region through analytical and numerical approaches. The results showed that, among all the models, the increase in Ec with MF f was a minimum for the FE model with the interphase region exhibiting a gradient. Furthermore, the Ec of the second-degree interphase model was the closest value to the results obtained experimentally. The results confirm that the properties of the interphase region need to be decided as similar as possible to experimental value when using the FE model. In order to obtain more accurate results in the future, performing the MD simulation and applying variables such as the fiber orientation will have to be considered with the FE calculation.