# On the Effect of Unit-Cell Parameters in Predicting the Elastic Response of Wood-Plastic Composites

Authors: Fatemeh Alavi; Amir Hossein Behravesh; Abbas S. Milani; Davoud Karimi
Journal: Journal of Engineering
 (2013)

---

## Abstract

This paper presents a study on the effect of unit-cell geometrical parameters in predicting elastic properties of a typical wood plastic composite (WPC). The ultimate goal was obtaining the optimal values of representative volume element (RVE) parameters to accurately predict the mechanical behavior of the WPC. For each unit cell, defined by a given combination of the above geometrical parameters, finite element simulation in ABAQUS was carried out, and the corresponding stress-strain curve was obtained. A uniaxial test according to ASTM D638-02a type V was performed on the composite specimen. Modulus of elasticity was determined using hyperbolic tangent function, and the results were compared to the sets of finite element analyses. Main effects of RVE parameters and their interactions were demonstrated and discussed, specially regarding the inclusion of two adjacent wood particles within one unit cell of the material. Regression analysis was performed to mathematically model the RVE parameter effects and their interactions over the modulus of elasticity response. The model was finally employed in an optimization analysis to arrive at an optimal set of RVE parameters that minimizes the difference between the predicted and experimental moduli of elasticity.


---

## Body

## 1. Introduction

It is well known that particles arrangement in a matrix affects the local stress/strain field in the ensuing composite, and in turn influences the macrolevel behavior of the material. Fortunately, effective mechanical properties of some heterogeneous materials rely on the average response of their microstructures and properties of their individual constituents [1]; hence, basic micromechanics theories can be sufficiently used in analyzing such materials. The connection between the macro- and mesolevel studies is traditionally viewed via the concept of the representative volume element (RVE). However, today little quantitative knowledge is available about minimum RVE sizes for various engineering materials. Several attempts have been made to determine the optimum size of an RVE [2]. Numerical-statistical, analytical approaches, and experimental observations are three methods applied by researchers to determine this size.In the numerical-statistical approach, multiple realizations, finite element simulations of materials unit cells, and appropriate statistical procedures are employed. Kanit et al. [3] proposed a quantitative definition of RVE size which was associated with a given precision in estimating the desired overall properties and the number of realizations for a given volume of microstructure. Eventually, they demonstrated how a minimal volume size for the computation of effective properties can be determined depending on the chosen precision and number of realizations. Trias et al. [4] analyzed the random distribution of fibers by means of optical microscopy, and the obtained images were used to generate realizations of statistical representative volume elements (SRVEs) at microscale. They solved finite element models of real-microstructure SRVEs with arbitrary boundary conditions to obtain probability density functions of stress, strain, and dilatational energy density and then related the results to stress tensor at any point in the macroscale by means of a two-scale approach. Other numerical-statistical approaches based on setting a tolerance for the scatter in results are given by Vinogradov [5]. Monte-Carlo simulation is a related approach which was successfully applied by Ostoja-Starzewski [6] and Gusev [7] to predict the overall elastic constants of the studied periodic composite.Among analytical approaches, Drugan and Willis [8] employed an explicit nonlocal constitutive equation by considering averaged strain fields varying with the position of material points. A micromechanics model was also proposed by Jiang et al. [1] for studying the effective elastic modulus of composites containing regularly distributed sphere particles. Three typical particle arrangements in the forms of simple cubic lattice, body-centered cubic lattice, and face-centered cubic lattice were investigated. In the case of irregular-shaped particles, Li and Wongsto [9] derived new unit cells capable of dealing with problems involving reinforcing particles of irregular geometries and local imperfections such as debonding and microcracks in the matrix. Boundary conditions for their proposed unit cells were derived from appropriate considerations of the conditions resulting from translational symmetry transformations. Giraud et al. [10] investigated the arbitrarily oriented ellipsoidal inhomogeneities to determine the macroscopic poroelastic properties of transversely isotropic geomaterials or rock-like composites. They mainly dealt with separating the effect of matrix anisotropy and of inhomogeneities in fiber orientation distribution and shape. In the RVEs based on experimental observations, experimental analysis often involves the selection of particular sample geometries for mechanical testing and subjecting specimens to image analysis after testing. Graham and Yang [11] and Shan and Gokhale [12] employed this approach for an instance. In the case of wood-plastic composites, analytical macromodels have been employed by Hugot and Cazaurang [13] to predict the WPC effective properties.Based on the above review, RVE modeling has been proven as an efficient approach to represent mechanical properties of composites including WPCs, as micromechanical modeling of a whole material structure would be computationally costly and often infeasible. On the other hand, the determination of optimum RVE parameters including its dimensions, particle size, and orientation is of major concern during such an analysis. Having appropriate RVE parameters, the analysis can further proceed to predict a composite’s effective properties. Majority of the micromodels in the current literature of WPCs encompass only one particle, which precludes the interaction effect of particles at the microscale. In this investigation, for the first time to address this issue, two particles of different geometrical specifications are included in a WPC RVE. Particles are modeled using elliptical geometries, the sizes of which together with the RVE size, location, and the angles of the ellipses center lines are considered as study parameters. Main and interaction effects of the parameters on the ensuing macrolevel modulus of elasticity (MOE) are investigated using finite element analysis. The Taguchi method of DOE (design of experiments) is employed to determine effects of RVE parameters. Regression analysis is performed to mathematically model the MOE as a function of RVE parameters. Finally, sequential quadratic programming is employed to optimize the RVE parameters to best match the modeled MOE with the prepared WPC sample with experimentally determined mechanical properties.


## 2. RVE Modeling and Designing Virtual Experiments (Simulations)

A typical wood plastic composite sample is shown in Figure1, where elliptic particles are noticeable. Accordingly, a two-dimensional plane stress RVE was constructed (Figure 2) with two elliptic particles. The sizes of the particles are assumed to be no larger than 120 micrometer in order to be within actual range of measurements and also to ensure that the two adjacent particles are separate during the subsequent simulations of the WPC unit cell. The mechanical properties of the composite are given in Table 1.Table 1
Mechanical properties of the modeled WPC composite.
Composite components
Elastic modulus (MPa)
Poisson’s ratio
Wood particles
20000
0.34
HDPE
834
0.3Figure 1
A typical wood plastic composite sample.Figure 2
A general RVE of the WPC sample.In order to understand the effect of geometrical parameters of RVE (as shown in Figure2) on the elastic response of the unit cell, nine geometrical parameters along with their levels were chosen (Table 2). The parameters (also called factors) include the width of RVE (Lx), length of RVE (Ly), location of particles (x), the major and minor diameters of each ellipse (a,b), and the inclination angle of the two elliptic particles (α). The major diameter (a) is considered as an independent parameter with four possible levels as listed in Table 2. The minor diameter, however, is obtained using the concept of constant aspect ratio, defined as the ratio of the minor diameter to the major one. Study of typical WPC samples suggests the range of aspect ratio from 2 to 5 based on Figure 1. Hence, four levels of b were adopted using the aforementioned range:
(1)b=aζ,
where ζ indicates the aspect ratio. Inclination angle of each particle designated by α (deg) is also considered to be an independent parameter. The levels of this factor in Table 2 were chosen such that they can cover possible angles between −90 and 90 degrees. The location of each particle, denoted by x, is defined as the horizontal distance of the center of particle from the middle of the RVE. The width of the RVE (Lx) is determined as
(2)Lx=2×{max(xm1,xm2)+t},
where t is the shortest horizontal distance from the left vertical side of the RVE to the particle as shown in Figure 2. xm is the horizontal distance from the center of the particle to a vertical line which osculates the ellipse at the furthest point. xm was then calculated from the following geometrical relation:
(3)xm=|a2cos(α)cos(β)-b2sin(α)sin(β)|,β=tan-1(-batan(α)),
which was obtained by maximizing the horizontal coordinate of the ellipse’s parametric equation as follows:
(4)X=xc+a2cos(α)cos(τ)-b2sin(α)sin(τ),xc and τ denote horizontal coordinate of the center of ellipse and a free angle variable, respectively.Table 2
List of factors and their levels used in the DOE study.
Factors
Level 1
Level 2
Level 3
Level 4
a
1 (μm)
60
80
100
120
a
2 (μm)
60
80
100
120
ζ
1
2
3
4
5
ζ
2
2
3
4
5
α
1 (deg)
−80
−30
0
60
α
2 (deg)
−80
−30
0
60
x
1 (μm)
0
0.75
1.5
0
x
2 (μm)
0
0.75
1.5
0
t (μm)
1.7
2.2
2.7
3.2d
X
/
d
τ
=
0 yields the optimum angle β in (3). The length of the RVE (Ly) is calculated considering the fact that the volume fraction of the RVE is that of the WPC sample. Accordingly, Ly is obtained from (5),
(5)ϑLxLy=π4(a1b1+a2b2),
whereϑ is the wood volume fraction (20%). Next, to run different possible combinations of the nine factors in Table 2 in a computationally economic manner, Taguchi’sL32′ orthogonal matrix was adopted [14]. Based on the actual values of the nine four-level factors in Table 2, the ensuing set of 32 computer experiments can be summarized as in Table 3. In each experiment (FE simulation), the modulus of elasticity is to be calculated as the response variable of the DOE study. In the FE model, particles and the matrix are meshed using CPS4R elements in ABAQUS. The “equation” constraint in the interaction module of ABAQUS was used to apply a periodic boundary condition on the RVE. To this end, the four-sided surfaces are defined as sets, and the degree of freedom of nodes on each surface (set) was constrained by an equation to the degree of freedom of a reference point. This constraint enforces all the nodes on the side surface to have the same displacement in the specified direction. As shown in Figure 3, reference points 1 and 2 are constrained in the x (horizontal) and y (vertical) direction, respectively. The wood particle and matrix materials are considered to be isotropic elastic with the mechanical properties shown in Table 1. Figure 3 shows a typical RVE simulation result for run 12. The plane thickness is considered to be unity; thus dividing the force at reference point 1 by the displacement of reference point 3 yields the effective modulus of elasticity of the RVE in the horizontal direction.Table 3
Set of runs resulted from the TaguchiL32′ and actual factor levels in Table 2.
#
a
1 (μm)
a
2 (μm)
ζ
1
ζ
2
α
1 (deg)
α
2 (deg)
x
1 (μm)
x
2 (μm)
t (μm)
MOE (GPa)
1
60
60
2
2
−80
−80
0
0
1.7
3.93
2
60
80
3
3
−30
−30
0.75
0.75
2.2
2.84
3
60
100
4
4
0
0
1.5
1.5
2.7
2.33
4
60
120
5
5
60
60
0
0
3.2
2.4
5
80
60
2
3
−30
0
1.5
0
3.2
3.16
6
80
80
3
2
−80
60
0
1.5
2.7
3.2
7
80
100
4
5
60
−80
0
0.75
2.2
4.69
8
80
120
5
4
0
−30
0.75
0
1.7
2.34
9
100
60
3
4
60
−80
0.75
1.5
3.2
2.68
10
100
80
2
5
0
−30
0
0
2.7
4.05
11
100
100
5
2
−30
0
0
0
2.2
4.659
12
100
120
4
3
−80
60
1.5
0.75
1.7
3.51
13
120
60
3
5
0
0
0
0.75
1.7
3.1
14
120
80
2
4
60
60
1.5
0
2.2
2.54
15
120
100
5
3
−80
−80
0.75
0
2.7
4.99
16
120
120
4
2
−30
−30
0
1.5
3.2
3.63
17
60
60
5
2
60
−30
1.5
0.75
2.7
2.92
18
60
80
4
3
0
−80
0
0
3.2
3.55
19
60
100
3
4
−30
60
0
0
1.7
3.47
20
60
120
2
5
−80
0
0.75
1.5
2.2
2.21
21
80
60
5
3
0
60
0
1.5
2.2
2.22
22
80
80
4
2
60
0
0.75
0
1.7
4.62
23
80
100
3
5
−80
−30
1.5
0
3.2
2.25
24
80
120
2
4
−30
−80
0
0.75
2.7
4.46
25
100
60
4
4
−80
−30
0
0
2.2
4.62
26
100
80
5
5
−30
−80
1.5
1.5
1.7
2.03
27
100
100
2
2
0
60
0.75
0.75
3.2
5.25
28
100
120
3
3
60
0
0
0
2.7
3.38
29
120
60
4
5
−30
60
0.75
0
2.7
1.86
30
120
80
5
4
−80
0
0
0.75
3.2
3.27
31
120
100
2
3
60
−30
0
1.5
1.7
3.22
32
120
120
3
2
0
−80
1.5
0
2.2
4.69Figure 3
Sample FE results of the RVE for run 12.


## 3. Results and Discussion

Modulus of elasticity (MOE) values obtained for FE runs are presented in Table3. Subsequently, main effects plots can be illustrated as in Figure 4(a), where the mean value of the MOE over levels of each factor from Table 3 has been calculated. It is observed from Figure 4(a) that the major diameter of each wood particle at 100 mm results in a peak MOE. The aspect ratio of the second particle is inversely related to the RVE elastic modulus: as the value of ζ2 increases the modulus of elasticity drops significantly. The inclination angle of the second particle also has an inverse effect on MOE. It has produced the largest value of MOE at −80 degrees. Increasing the value of x1 decreases MOE, while increasing x2 increases the response from 3.53 to 3.75 (GPa) and then decreases to 2.67 (GPa). MOE has a peak value at x2 = 0.75 (mm). It can also be inferred from Figure 4(a) that the aspect ratio and the inclination angle of the first particle together with Lx, indicate small effects on the elastic response of the RVE. This is despite the fact that the inclination angle of the second particle shows higher effect, meaning that the effect of interactions between the parameters of particles may affect their individual role on the WPC response. To arrive at a prediction model and identify the significant effects more quantitatively, a multiple regression analysis was performed using the data in Table 3.(a) Main effect plots for the mean MOE response; (b) interaction plot betweena1 and a2; (c) interaction between α1 and ζ1; (d) interaction between x1 and x2.
(a)
(b)
(c)
(d)


### 3.1. Mathematical Model of MOE

To obtain a mathematical model of MOE, a polynomial of order two was employed. There are 35 double-interaction terms in this analysis. However, there was no need to consider all the model terms. Some terms have significant interactions, and some have not. For instance, the larger ellipse diameter and the angle of particle have high interaction. For small angles, an increase in the larger diameter significantly affects the MOE. However, for large values of particle angles the increase in the larger diameter will have an insignificant effect on MOE, as MOE is calculated in the horizontal direction. In the same fashion, it can be assumed that there is a significant interaction between the larger diameter and the rest of the parameters. It is noteworthy that the location of the particles with respect to the center of the RVE does not possess a significant interaction with the rest of parameters. However, the locations of two particles interact with each other significantly. After removing terms with negligible effects, the prediction model reads the following:(6)MOE=34.2+0.171a1-0.388a2-2.55ζ1-1.35ζ2+0.130α1+0.0294α2-1.21x1-2.01x2-7.06t-0.000168a1a2-0.000654a1α1+0.000159a1α2-0.000514a2α1-0.00003a1t-0.000661a2α2+0.0750a2t+0.202ζ1ζ2+0.0156α1t-0.000394α1α2+0.00646α2t-0.00551a1ζ1-0.0362a1ζ2+0.0161a2ζ1+0.0363a2ζ2-0.0146ζ1α1-0.000759ζ1α2-0.00575ζ2α1-0.00207ζ2α2+1.24x2x1.The goodness of fit (R2 value) of 99.7% and the adjusted R2 of 95.3%, were obtained which demonstrate that the regression equation fits the data reasonably well. The simulated and fitted values are illustrated in Figure 5. Equation (6) was next used in an optimization analysis in MATLAB to find the optimal set of RVE parameters that can match the modulus of elasticity obtained through the prediction model and that of a tension test on the WPC sample with a 20 weight percentage of wood. Tensile property of the specimen (MOE) was determined according to ASTM D638-02a type V. The obtained stress-strain data from the tensile test are shown in Figure 6, from which the MOE in the elastic region can be estimated as follows. The curve fitting toolbox of MATLAB was used to fit the hyperbolic material model by Zawlocki [15] to the data points in Figure 6. This material model is described as follows:
(7)σ=m×tanh(nε),
or
(8)∂σ∂ε=mncosh2(nε),
in which m and n are coefficients that correlate stress (σ) and strain (ε) based on hyperbolic formulation. Asε approaches zero, it can be stated mathematically that
(9)∂σ∂ε→m×n=MOE(modulusofelasticity).m and n were found from the above mentioned curve fitting such thatm×n=MOE=4.8Gpa. To find the optimal set of RVE parameters that minimizes the difference between the FE simulation (predicted value) and the experimental value (4.8Gpa), the sequential quadratic programming (SQP) method [16] in MATLAB optimization toolbox was employed.Figure 5
Simulated versus fitted MOE at different RVEs.Figure 6
Stress-strain data from the tensile test on the WPC specimen (dots) and the prediction (solid line) by a hyperbolic material model.


### 3.2. Optimization of RVE Parameters

A general optimization formulation considered is as follows:(10)minimizef(x)subjectto:Ci(x)=0i=1,…,M,Ci(x)≤0i=M+1,…,N,
where x is the vector of design parameters, f(x) is the objective function which is the difference between the values obtained from (6) and the experimentally obtained modulus of elasticity, and the vector function Ci(x) is a vector containing the values of the equality and inequality constraints evaluated at x. The only constraints considered here is the lower and upper bounds of the parameters as follows:
(11)lowerbound=[60,60,2,2,-80,-80,0,0,1.7]T,upperbound=[120,120,5,5,60,60,1.5,1.5,3.2]T.The general solution to the above constrained optimization problem is to transform the problem into an easier subproblem that can be solved and used as the basis of an iterative process. Active set strategy (also known as projection method) is employed to perform the iterations. The solution procedure includes two phases. The first phase involves the calculation of a feasible point, if any. The second phase involves the generation of an iterative sequence of feasible points that converge to the solution (more details can be found in [16–18]).The resulting optimum sets of RVE parameters are presented in Table4. Equation (6) gives MOE = 4.79 GPa, and FE simulation yields MOE = 4.83 GPa for the optimal values presented in Table 4. It is clear from Table 4 that in the obtained optimum set, the parameters for the two wood particles are not necessarily identical. This goes back to the present effect of interactions between geometrical parameters of one particle with the geometrical parameters of the other particle on the overall response of the RVE, which in fact was the main motivation of the work to illustrate through the performed case study.Table 4
The obtained optimal set of RVE parameters.
RVE parameters
Optimal values
a
1 (μm)
93
a
2 (μm)
120
ζ
1
4
ζ
2
5
α
1 (degree)
60
α
2 (degree)
−80
x
1 (μm)
1
x
2 (μm)
0.95
t (μm)
2.7


## 3.1. Mathematical Model of MOE

To obtain a mathematical model of MOE, a polynomial of order two was employed. There are 35 double-interaction terms in this analysis. However, there was no need to consider all the model terms. Some terms have significant interactions, and some have not. For instance, the larger ellipse diameter and the angle of particle have high interaction. For small angles, an increase in the larger diameter significantly affects the MOE. However, for large values of particle angles the increase in the larger diameter will have an insignificant effect on MOE, as MOE is calculated in the horizontal direction. In the same fashion, it can be assumed that there is a significant interaction between the larger diameter and the rest of the parameters. It is noteworthy that the location of the particles with respect to the center of the RVE does not possess a significant interaction with the rest of parameters. However, the locations of two particles interact with each other significantly. After removing terms with negligible effects, the prediction model reads the following:(6)MOE=34.2+0.171a1-0.388a2-2.55ζ1-1.35ζ2+0.130α1+0.0294α2-1.21x1-2.01x2-7.06t-0.000168a1a2-0.000654a1α1+0.000159a1α2-0.000514a2α1-0.00003a1t-0.000661a2α2+0.0750a2t+0.202ζ1ζ2+0.0156α1t-0.000394α1α2+0.00646α2t-0.00551a1ζ1-0.0362a1ζ2+0.0161a2ζ1+0.0363a2ζ2-0.0146ζ1α1-0.000759ζ1α2-0.00575ζ2α1-0.00207ζ2α2+1.24x2x1.The goodness of fit (R2 value) of 99.7% and the adjusted R2 of 95.3%, were obtained which demonstrate that the regression equation fits the data reasonably well. The simulated and fitted values are illustrated in Figure 5. Equation (6) was next used in an optimization analysis in MATLAB to find the optimal set of RVE parameters that can match the modulus of elasticity obtained through the prediction model and that of a tension test on the WPC sample with a 20 weight percentage of wood. Tensile property of the specimen (MOE) was determined according to ASTM D638-02a type V. The obtained stress-strain data from the tensile test are shown in Figure 6, from which the MOE in the elastic region can be estimated as follows. The curve fitting toolbox of MATLAB was used to fit the hyperbolic material model by Zawlocki [15] to the data points in Figure 6. This material model is described as follows:
(7)σ=m×tanh(nε),
or
(8)∂σ∂ε=mncosh2(nε),
in which m and n are coefficients that correlate stress (σ) and strain (ε) based on hyperbolic formulation. Asε approaches zero, it can be stated mathematically that
(9)∂σ∂ε→m×n=MOE(modulusofelasticity).m and n were found from the above mentioned curve fitting such thatm×n=MOE=4.8Gpa. To find the optimal set of RVE parameters that minimizes the difference between the FE simulation (predicted value) and the experimental value (4.8Gpa), the sequential quadratic programming (SQP) method [16] in MATLAB optimization toolbox was employed.Figure 5
Simulated versus fitted MOE at different RVEs.Figure 6
Stress-strain data from the tensile test on the WPC specimen (dots) and the prediction (solid line) by a hyperbolic material model.


## 3.2. Optimization of RVE Parameters

A general optimization formulation considered is as follows:(10)minimizef(x)subjectto:Ci(x)=0i=1,…,M,Ci(x)≤0i=M+1,…,N,
where x is the vector of design parameters, f(x) is the objective function which is the difference between the values obtained from (6) and the experimentally obtained modulus of elasticity, and the vector function Ci(x) is a vector containing the values of the equality and inequality constraints evaluated at x. The only constraints considered here is the lower and upper bounds of the parameters as follows:
(11)lowerbound=[60,60,2,2,-80,-80,0,0,1.7]T,upperbound=[120,120,5,5,60,60,1.5,1.5,3.2]T.The general solution to the above constrained optimization problem is to transform the problem into an easier subproblem that can be solved and used as the basis of an iterative process. Active set strategy (also known as projection method) is employed to perform the iterations. The solution procedure includes two phases. The first phase involves the calculation of a feasible point, if any. The second phase involves the generation of an iterative sequence of feasible points that converge to the solution (more details can be found in [16–18]).The resulting optimum sets of RVE parameters are presented in Table4. Equation (6) gives MOE = 4.79 GPa, and FE simulation yields MOE = 4.83 GPa for the optimal values presented in Table 4. It is clear from Table 4 that in the obtained optimum set, the parameters for the two wood particles are not necessarily identical. This goes back to the present effect of interactions between geometrical parameters of one particle with the geometrical parameters of the other particle on the overall response of the RVE, which in fact was the main motivation of the work to illustrate through the performed case study.Table 4
The obtained optimal set of RVE parameters.
RVE parameters
Optimal values
a
1 (μm)
93
a
2 (μm)
120
ζ
1
4
ζ
2
5
α
1 (degree)
60
α
2 (degree)
−80
x
1 (μm)
1
x
2 (μm)
0.95
t (μm)
2.7


## 4. Conclusion

In this investigation, a new Representative Volume Element (RVE) containing two particles were employed to predict the modulus of elasticity (MOE) for a wood plastic composite. The effect of RVE geometrical parameters and their interactions on the MOE response of the RVE were studied. A regression model was used within an optimization framework to obtain the optimal set of RVE parameters that can match the experimentally determined modulus of elasticity of the WPC sample. It was found that in the presence of multiple particles in the RVE model (especially for WPCs with a high particle volume fraction), the optimum sets of parameters for the particles are not necessarily identical. This observation was attributed to the presence of possible interaction effects between the particles. Including more particles in the RVE modeling of WPC materials adds to the computational cost of simulations; however, at the same time it gives a greater degree of freedom for better prediction of the material response. Future work includes a full-scale simulation of a tensile test with the developed RVE with multiple particles along with the identified hyperbolic model. Finally, it should be mentioned that although many industrial applications may demand for three-dimensional analysis of composite structures, for the sake of simplicity a two-dimensional problem was considered in this investigation. The proposed method can be expanded to three-dimensional analysis by taking into account the following additional parameters: third direction aspect ratio, inclination angle towards third direction, and particle locations in third direction.



---

*Source: 456398-2013-02-26.xml*