- Email: [email protected]
Multiscale analysis of carbon nanotube-reinforced nanofiber scaffolds
Composite Structures 93 (2011) 1008–1014
Contents lists available at ScienceDirect
Composite Structures journal homepage: www.elsevier.com/locate/compstruct
Multiscale analysis of carbon nanotube-reinforced nanofiber scaffolds V.U. Unnikrishnan, G.U. Unnikrishnan, J.N. Reddy * Advanced Computational Mechanics Laboratory, Department of Mechanical Engineering, Texas A&M University, College Station, TX 77843-3123, USA
a r t i c l e
i n f o
Article history: Available online 4 July 2010 Keywords: Carbon nanotube Nonlinear mathematical homogenization Hyperelastic material Multiscale analysis
a b s t r a c t Biopolymers play a significant role in tissue engineering, as they simulate the physiological environment required for the development of tissue cultures. Use of carbon nanotube polymeric scaffolds for tissue engineering applications has gained attention recently due to the enhanced mechanical properties of carbon nanotubes. In this paper, a hierarchical approach by studying the atomistic properties of carbon nanotube based polymers using molecular dynamics and coupling the scales through complex multiscale nonlinear hyperelastic material-based mathematical homogenization models are developed. These homogenization methods offer a systematic and rigorous treatment of up-scaling the properties from the micro or nanoscale to the macroscales. The material constitutive properties estimated using the developed methods for nanotube polymeric scaffold materials show excellent comparison with experimental studies. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction Carbon nanotube-reinforced nanofibers as scaffolds in bioreactor and other tissue engineering applications has received attention due to their enhanced material properties. The extremely high mechanical strength of carbon nanotubes is found to result in a nanofiber capable of withstanding tremendous mechanical stresses. The properties of these nanotubes can be obtained from various experimental investigations, but a thorough understanding of the properties and its interactions with the surrounding materials can only be obtained by advanced computational models. In this paper, we analyze the mechanical properties of carbon nanotube based Poly-Vinyl Alcohol (PVA) polymeric nanofibers [1–3], hierarchically from its constituents (i.e., PVA as matrix phase and nanotube as fiber phase). PVA is a water-soluble polyhydroxy polymer and has been widely used as a scaffold supporting material for tissue engineering applications, as it is bio-compatible and nontoxic and provides mechanical stability and flexibility [4–7]. The mechanical properties of atomistic scale nanotube system are estimated by molecular dynamic (MD) simulations. Bridging the spatial scales for a continuum nanotube approximation is carried out by volume averaging of the constituent properties of the atomistic system, using equivalence of the strain energy due to deformation of the continuum system and change in potential energy during an iso-thermal mechanical straining process of the atomistic system [8]. The information obtained from the volume averaging of the nanotube ensemble is then used in the evaluation * Corresponding author. Tel.: +1 979 862 2417. E-mail address: [email protected] (J.N. Reddy). 0263-8223/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruct.2010.06.019
of the effective macroscale properties of the nanofibers. Studies on the mechanical, electrical and thermal behavior of nanotubes in literature were primarily focused on atomistic simulations, however the computational adaptation to model macroscopic problems have not been completely established [9,10]. Continuum mechanics principles could be used to bridge the different scales of interest and provide solutions not only to the fine scale problems but also to the macroscale effects. This work is significant for the fact that, there have been very few studies on homogenization of nanotube based polymer scaffolds. Therefore the significance of this paper is that we intend to carry out both linear and nonlinear homogenization of nanotube based polymeric scaffolds. The overall effective properties of the nanofiber thus obtained can be used in various mathematical analysis, like simulation of nutrient transport within the nano-scaffold, etc. The mathematical representations of the behavior of a material, which are often validated against experimental results are known as constitutive relations of the material. It has been observed that the differences between the actual behavior and computational predictions are often due to the inaccurate representation of the constitutive behavior of the material. In this work, we consider a hyperelastic material model based on nonlinear homogenization to determine the constitutive properties of the nano-scaffold material. This is a significant improvement over conventional studies on nano-scaffolds as most of these studies were based on the assumption that the scaffold material is linearly elastic [11,12,8], and therefore accurate estimation of scaffold material properties is not possible. However, it has been found that the nano-scaffold materials exhibit nonlinear material behavior and therefore nonlinear homogenization is extremely important [13,1,6]. The significance of this
1009
V.U. Unnikrishnan et al. / Composite Structures 93 (2011) 1008–1014
work lies in the use of a multi-physical modeling strategy involving complex material geometry leading to an efficient multiscale–multi-physical computational model. This paper is organized as follows. Following this introduction, the analysis of nanotube-reinforced polymeric scaffolds is carried from the lower-order to the higher-order using MD simulations and scaling the properties by various mathematical homogenization methods. The second part of the paper deals with the nonlinear homogenization methods using hyperelastic material formulation, where the properties of the polymeric scaffold are obtained by dealing with the equivalence of the strain energy of the composite from the energy of the individual components. The effect of various fiber orientations on the overall effective properties of the nanotube–polymer composite system is also discussed. The paper concludes with a discussion of the result and conclusions.
Dt Dt fi ðtÞ ¼ vi t þ vi t þ Dt 2 2 mi
ð3Þ
where Dt is the time step, and fi(t) is the total force acting on particle i at time t. From the updated particle position, the interatomic forces are computed from the gradient of the potential energy field U with respect to the atomic coordinates ri as
fi ðt þ DtÞ ¼
@U @ri ðtþDtÞ
ð4Þ
The total potential energy due to straining, which can be considered as the elastic strain energy, can be expanded in Taylor’s series for small displacements, with the initial position represented by the equilibrium configuration. The elastic moduli tensor can then be written as (no sum on repeated subscripts)
3 ! 2 6 1 X 6 d U 1 dU a b c d 1 dU a c 7 7 aij aij aij aij þ dbd ¼ a a 7 6 2 2NXa j–i 4 dr ij r ij dr ij rij dr ij ij ij 5 |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} 2
2. Analysis of nanotube–polymer composite systems
C abcd
2.1. MD simulation of single-walled nanotube
Aac ¼0
Single-walled nanotubes are usually idealized as hollow cylinder formed from rolled graphene sheets, and each of the carbon atoms forming the tubules has three nearest neighboring bonds. In a fully relaxed structure, the angles between these bonds depend on the radius of the cylinder as well as on their orientation. The fundamental CNT structure can be classified into three categories according to their helicity: armchair, zigzag, and chiral. The knowledge of structure and molecular motion in nanoparticles are essential to understand properties of practical interest, and molecular dynamics (MD) has been widely used for the determination of these properties [14,15]. MD enables the computation of trajectories of atoms by integrating Newton’s second law of motion. The Hamiltonian, which represents the total energy of the system, is a function of position and momentum of the particles, and determines the trajectories of the particles [14,16]. For an isolated N-particle system, the Hamiltonian is given as the sum of the potential and kinetic energies
Hðr N ; pN Þ ¼
1 X 2 p þ UðrN Þ 2m i i
ð1Þ
where U(rN) is the potential energy from intermolecular interactions, which is represented as a function of the spatial ordinate P 2 1 r N ; 2m i pi is the kinetic energy, which represents the momentum pi of particle i with mass m; and pi is a function of the absolute temperature. The total potential energy (U) of the computational unit cell is given by the sum total of the valence bond energies ðV Bij Þ and the non-bonding interactions ðV NB ij Þ
U¼
X Xh j
V Bij þ V NB ij
i
ð5Þ
ð2Þ
j>i
In the atomistic level the molecular interactions are modeled using pair of potentials to respond to externally applied disturbance, and the reliability of MD simulation depends on the type of potential functions used. Tersoff–Brenner empirical bond order potentials are excellent for the simulation of carbon based materials and are used in the present work [15,14]. These bond order potentials are capable of describing the changes in the bonding between atoms, but does not describe long range interactions [17]. The force of attraction or repulsion experienced by each molecule is obtained from the gradient of the potential field which is used in calculating the updated position of the atoms using the velocity-Verlet time integration scheme. At the start of each time step, the updated velocities vi(t) are calculated for each particle by
Here Cabcd are the elastic moduli; U = U(rij) is the potential energy as a function of the inter-atomic distance rij; Aac is the internal stress tensor, which at equilibrium is equal to zero; Xa is the average volume of an atom; N is the number of atoms; dab is the Dirac delta function; a, b, c, d take on the values 1, 2, 3 which are the spatial dimensions and aaij is the undeformed value of rij in the a-direction. The deformed length is given by uaij ¼ raij aaij and any change in the inter-atomic distance (as a result of deformation) from the equilibrium (minimum energy position) inter-atomic distance contributes to the strain experienced by each atom and is given by ubij ¼ aaij eab , where eab are the components of the homogeneous infinitesimal strain tensor associated with atoms i and j. In a composite structure, where the nanotube is dispersed in a polymer matrix, the orientation and morphological characteristics like diameter, length and chirality also causes a corresponding variation in the macroscopic properties of the composite. In the next step, we estimate the effective properties of the polymeric nanocomposite using mathematical homogenization based on the elastic properties of the nanotube obtained from MD simulations and the polymer obtained from experimental studies. 2.2. Homogenization of nanotube – polymeric nanofibrous scaffolds There are numerous homogenization models for the analysis of structures composed of two or materials [18–20]. In mathematical homogenization models the heterogeneous composite material is replaced by an equivalent homogeneous continuum using a suitable homogenization technique. The homogeneous continuum is based on a statistically equivalent homogeneous volume element, called the representative volume element (RVE) which represents the material in the lower-order (the microscale) and collectively represents the global structure. 2.2.1. 3D Open Cell Homogenization based on structural models was developed for foams and cellular materials based on the ‘‘unit cell” model. These models are valid for materials with porosity 70% or higher and therefore can be used in the present analysis as numerical studies have shown that polymeric scaffolds have wide range of porosity [11,8]. The effective modulus of the nanofibrous scaffold by the structural model is given as [21]
C ¼ 0:6495Cð1 UÞ3
ð6Þ
where C* is the effective modulus of the porous material with elastic modulus C and porosity U. In a 3D open cell material model the
1010
V.U. Unnikrishnan et al. / Composite Structures 93 (2011) 1008–1014
effective modulus of the porous material is related to the fiber modulus by C* = C(1 U)2 [21]. This model is also based on assumptions of high porosity; and gives good predictions of moduli for materials with porosities from 10 to 90%. Nanofibrous scaffolds usually falls within this range of porosity and therefore this model can give a conservative prediction of the effective elastic modulus [22]. It should however be emphasized that the micromechanical based homogenization methods usually give rigorous bounds compared to structural models. 2.2.2. Dewey–Mackenzie estimates For non-interacting inclusions in a composite medium, exact solutions were obtained by Dewey and Mackenzie, called the Dewey–Mackenzie (DM) estimates [23]. The DM estimates are based on the exact classical solution of a spherical cavity embedded in an infinite elastic medium but the major drawback of the DM estimates is that they are dilute approximations. In the DM estimates, the elastic modulus is a function of Poisson’s ratio of the matrix and the porosity of the scaffolds given by
C DM
3ð1 mm Þð1 5mm Þ ¼ 1þ U Cm 2ð7 5mm Þ
ð9Þ
where r = rT is the Cauchy stress tensor, b is the body force per unit mass [24,25]. Along with the energy, entropy balance equations, and boundary conditions, the set of closure equations called the constitutive relationships are required for solving the above sets of equations. Hyperelastic materials can undergo large elastic strains with essentially small compressibility and the deformation characteristic is defined by the strain energy function W [26,27]. For a hyperelastic material there exists strain energy function (W), such that the stress components (rij) are derivable as rij ¼ @@eW , where eij are ij the strain components. The components of deformation gradient @xi tensor F are given by F ij ¼ @X , where Xj and xi are the components j of the material and spatial coordinates, X and x, respectively. The strain-energy density function can then be expressed as a linear combination of the strain invariants I1, I2, and I3.
W¼
n X
wk ðI1 ; I2 ; I3 Þ
ð10Þ
k¼1
ð7Þ
where CDM is the effective DM modulus, Cm is the modulus of the matrix and mm is Poisson’s ratio of the matrix and U is the porosity of the scaffold. 2.3. Nonlinear homogenization of hyperelastic materials Let us consider a composite system comprising the matrix and fibers phases and each material being homogeneous and hyperelastic. Let the fiber be oriented a degrees along the length of the composite RVE as shown in Fig. 1. The constitutive behavior of the individual phases can be characterized by the stored-energy functions which are assumed to be objective, isotropic and strongly elliptic (rank one convex) functions of the deformation gradient F. Now, considering X0 R3 to be the reference configuration of a deformable body with no discontinuities, let the current position of a point X be X 2 X0, which refers to the material position in the body. Under the action of external or internal forces, the point X moves to the position x in the current (deformed) configuration such that x = X + u, where u is the displacement vector and the motion (or deformation mapping) is defined by the mapping /(X, t). The principle of conservation of mass, assuming that there is no addition or deletion of mass, is given by
@q þ r ðqv Þ ¼ 0 @t
@ðq v Þ þ r ðqv v Þ ¼ r r þ qb @t
ð8Þ
where q is the density, v is the velocity. The principle of conservation of linear momentum gives the equations of motion
where wk are convex scalar functions and n is the number of individual phases in the hyperelastic composite system. In a composite material with two or more participating phases (e.g., matrix, fiber, interphase, etc.), the elastic response is assumed to result from the contributions of the individual phases and their interactions. e can be diTherefore, the strain energy function of the composite W vided into three parts, namely
e ¼ Wm þ Wf þ Wfm W
ð11Þ
where W represents the strain energy and the superscripts m, f, fm denotes the contributions due to the matrix phase, fiber phase and fiber–matrix interactions respectively. However, for non-interacting phases, we may neglect the contributions from the fiber–matrix interactions. The strain energy in terms of the principal invariants are then given by
W ¼ WðI1 ; I2 ; I3 ; I4 ; I5 Þ 1 I1 ¼ tr C; I2 ¼ ½ðtrCÞ2 tr C; 2 I4 ¼ a0 C a0 ¼ k2F
ð12Þ I3 ¼ detC ð13Þ
I5 ¼ a0 C2 a0 where kF is the fiber stretch and C = FTF is the right Cauchy–Green strain tensor, {k1, k2, k3} are the principal stretches, and a0 is the direction of the normal given by a0 = [cosa, sina, 0]T. For a hyperelastic material, the stress–strain relationship is obtained by taking the derivative of the strain-energy with respect to the Green–Lagrange strain tensor E ¼ 12 ðC IÞ, and is given by
S¼2
@W p C1 @C
ð14Þ
where S is the second Piola–Kirchhoff stress [24,25] and p is the pressure, which is added to account for the incompressibility constraint [28,29]. The corresponding Cauchy stress tensor ðr ¼ 1J FSFT Þ associated with this incompressible hyperelastic composite material is given by Quapp and Weiss [30] and Reddy [24]
h
i
r ¼ 2 W1 B þ W2 ðI1 B B2 Þ þ W4 a a þ W5 ða Ba þ Ba aÞ p I ð15Þ
Fig. 1. Nanotube fiber orientation in a fiber–matrix composite system.
where Wi = oW/@Ii, i = 1, 2, 4, & 5 and p is the pressure, I is the identity tensor, and a = Fa0 (see Fig. 1). Let us now consider the matrix material to be a compressible, Neo-Hookean material, and the strain energy in the reference configuration given in Eq. (16)
V.U. Unnikrishnan et al. / Composite Structures 93 (2011) 1008–1014
l
Wm ¼
ðI1 3Þ þ
2
j 2
ðJ 1Þ2
ð16Þ
where J = det(F). When the material is incompressible (J = 1) the second term in the strain energy density function vanishes. Correspondingly, the strain energy of the fiber due to the elongation may be considered as a function of I4, and assuming it to be in the form given in Eq. (17).
( F
W ¼
C 2 ðI4 1Þ2 þ C 3 ðI4 1Þ4
I4 > 1
0
I4 6 1
ð17Þ
It can be assumed that the fiber does not experience any strain during contraction of the composite [31]. Substituting I4 ¼ k2F and ¼ 12 ðk2F 1Þ [24] into Eq. (17), the strain energy of the fiber is obtained in terms of strain and can be used to fit experimental data for fiber elongation.
WF ¼ 42F C 2 þ 4C 3 2F
ð18Þ
The Cauchy stress tensor for this incompressible finite elastic composite material is then given by Qiu and Pence [32]
r ¼ 2F
@ W @I1 @ W @I4 T F p FC1 FT þ @I1 @C @I4 @C
ð19Þ
where
@W l ¼ @I1 2
and
@W ¼ 2C 2 ðI4 1Þ þ 4C 3 ðI4 1Þ3 @I4
1011
Using the above relation, the pressure which is now a function of the fiber orientation angle a,and principal stretches k1, k2 can be obtained and can be used to estimate the total stress due to deformation in the composite. 3. Results and discussion The elastic modulus of a single-walled nanotube (SWNT) is determined by a series of MD simulations, where the minimum energy configuration of the CNT is determined using minimization of the total potential energy (PE) at 0 K. The relaxed configuration of the nanotube is now subjected to displacement increments and the structure is allowed to equilibrate over a number of time steps. The straining process induces changes in the force experienced by each atom and this force is used in calculating the updated position of the atoms. The CNT is strained at a constant rate to calculate the strain energy of deformation by minimizing the total potential energy at each increment and the difference in the potential energy gives the strain energy due to deformation. From the MD simulations, the elastic modulus of a single-walled nanotube is obtained as 0.94 TPa and this value lies closely with many theoretically and experimentally obtained elastic modulus in the literature [16]. The elastic modulus of a polymeric electrospun nanofiber with diameter from 20 to 40 nm with 100% Poly-Vinyl Alcohol (PVA) is
ð20Þ
The derivatives of I1 and I4 are given by
@I1 @I4 ¼I ¼ a0 a0 @C @C
ð21Þ
Considering uniaxial elongation of a composite with fiber oriented along the direction of deformation, the deformation gradient F is given by
2
F1
0
6 F¼4 0
F2
0
0
0
3
7 05 F3
ð22Þ
where F1, F2 and F3 are correspondingly the principal stretches k1, k2, and k3, respectively. The corresponding stresses in the composite may then be written as:
rxx ¼ l 4C 2 ð1 k21 Þ 8C 3 ð1 þ 3k21 þ 3k41 k61 Þ k21 p ryy ¼ l k22 p; rzz ¼ l k23 p
ð23Þ
Fig. 2. Variation of elastic modulus of nanotube-reinforced PVA scaffold with void volume fractions.
It should be noted that due to incompressibility constraints, the principal stretches have the following properties: k1k2k3 = 1 and ki > 0, i = 1, 2, 3 [32] and traction free conditions exists in the y- and z-directions. The principal stretches in these directions for the fiber orientation angle a = 0 may be obtained as
1 k2 ¼ k3 ¼ pffiffiffiffiffi k1
ð24Þ
Due to the traction free conditions in the y and z directions, the pressure is obtained as
p ¼ lk1 1
ð25Þ
Additionally, we may also consider the fiber in the composite to be oriented at an arbitrary angle along the x direction. Now, assuming k3 = 1 (i.e., the deformation in the z direction is negligible), the principal stretches in the other two directions can be obtained as a result of incompressibility of the material and is given in Eq. (26).
k2 ¼
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ð3 k21 Þ 2
ð26Þ
Fig. 3. Stress–strain curve and its comparison between the experimental data and manufactured solution for a single-walled nanotube.
1012
V.U. Unnikrishnan et al. / Composite Structures 93 (2011) 1008–1014
Fig. 4. Variation of k2 with k1.
approximately 10 GPa. Thus, assuming a nanotube volume fraction of 0.005, the effective modulus of the nanofiber by Mori–Tanaka method is obtained as 10.90 GPa and by Voigt bound is 14.65 GPa. The details of Mori–Tanaka method and Voigt bounds may be found in standard literature and is not discussed here [33]. These values are very close to the experimentally obtained value 15 GPa by Wong et al. [3] and 16.24 GPa (see Fig. 5a of [6]). The difference in the predicted and experimentally obtained values could be due to the micro-structural formation during electrospinning, which could not be accounted in the mathematical homogenization models. The effective property of the scaffold is now estimated using the 3D-open cell model [8] and the DM estimates. Based on the homogenization studies from the previous sections, an average modulus of 10.9 GPa for the nanofiber and mm ¼ 0:4 is considered. The variation of elastic modulus of the scaffold with respect to various void volume fractions is shown in Fig. 2 and the inset shows the variation of the elastic modulus at the void limits. Now, considering the experimental investigation of elastic properties of nanotube bundles under tensile loads, the average
Fig. 5. q Nonlinear ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi stress–strain curve for tensile loading of nanotube-reinforced matrix composite for fiber various orientation angles for l = 3.33 GPa and 12.01 GPa and k2 ¼ þ
1 ð3 2
k21 Þ.
Fig. 6. Nonlinear stress–strain curve for tensile loading of nanotube-reinforced matrix composite for fiber various orientation angles for l = 3.33 GPa and 12.01 GPa and k2 = 2 k1.
V.U. Unnikrishnan et al. / Composite Structures 93 (2011) 1008–1014
1013
obtained stress–strain curve of nanotube systems and up-scaling to the macroscale properties. The effective property of the nanofiber scaffold estimated using both the micromechanical and hyperelastic material homogenization formulations were compared with experimental data and are found to be in good agreement. Acknowledgement The research reported herein was supported by the Oscar S. Wyatt Endowed Chair account at Texas A&M University. References
Fig. 7. Stress strain variation for uniaxial extension for l = 3.33 GPa and 12.01 GPa.
elastic moduli obtained by Yu et al. was 1.002 TPa [34] and the stress–strain curve reported by Qian et al. [34] is shown in Fig. 3. The experimentally obtained values are now fitted to the derived strain energy function in Eq. (18) to obtain the values of C2 and C3 and the curve-fitted stress–strain curve is denoted by the smooth line in Fig. 3. The elastic modulus of the matrix is found to vary from 10 GPa to 33 GPa [35,36] and considering the matrix to be incompressible, the shear modulus was found to range from 3.33 GPa to 12.01 GPa. Now, considering the orientation of the fiber along the length of the composite system, various stress–strain qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi profiles are obtained for k2 ¼ 12 ð3 k21 Þ (see Fig. 4) and is given in Fig. 5. However, various experimental studies have shown that k1 and k2 are often related by k2 = 2 k1 [31] (see Fig. 4) and the corresponding stress–strain variation is shown in Fig. 6, for various orientation angles of the fiber. It can be seen from the above figures that the stress–strain curve follows a nonlinear pattern and the difference in orientation angles severely affects the overall stress profile. The estimated stress strain curve is also compared with the experimental works of Xu et al. [36] and it can be seen that there is excellent correspondence with the experimentally determined data. Further, it can be seen that greater correlation is seen for higher shear modulus, which is the experimentally determined shear modulus for neat PVA fibers. It should also be noted that the nonlinearity is maximum when the fiber is oriented along the length of the composite RVE; this is because of the dominating nonlinear behavior of the fiber in the loaded direction. Stress–strain variation under perfect uniaxial extension, where the principal stretches are related as given in Eq. (24) (also, see Fig. 4) and the fiber oriented along the direction of deformation is also studied as shown in Fig. 7 for shear moduli equal to 3.33 GPa and 12.01 GPa. It can be seen that the nonlinear effects are significant at higher strain states. 4. Conclusions A multiscale modeling strategy to estimate the effective elastic modulus of carbon nanotube based polymeric nanofiber scaffolds is developed. The analysis is carried out from the atomistic scale comprising the nanotube using MD simulations, and the properties are scaled-up using various mathematical homogenization schemes. In this study, a nonlinear hyperelastic material based homogenization is also developed utilizing the experimentally
[1] Jeong JS, Moon JS, Jeon SY, Park JH, Alegaonkar PS, Yoo JB. Mechanical properties of electrospun PVA/MWNTs composite nanofibers. Thin Solid Films 2007;515(12):5136–41. [2] Wong KKH, Hutter JL, Zinke-Allmang M, Wan W. Physical properties of ion beam treated electrospun poly(vinyl alcohol) nanofibers. Eur Polym J 2009;45(5):1349–58. [3] Wong KKH, Zinke-Allmang M, Hutter JL, Hrapovic S, Luong JHT, Wan W. The effect of carbon nanotube aspect ratio and loading on the elastic modulus of electrospun poly(vinyl alcohol)–carbon nanotube hybrid fibers. Carbon 2009;47(11):2571–8. [4] Nho YC, Moon SW, Lee KH, Park CW, Suh TS, Jung YJ, et al. Evaluations of poly(vinyl alcohol) hydrogels cross-linked under Y-ray irradiation. J Ind Eng Chem 2005;11(1):159–64. [5] Mansur HS, Costa HS. Nanostructured poly(vinyl alcohol)/bioactive glass and poly(vinyl alcohol)/chitosan/bioactive glass hybrid scaffolds for biomedical applications. Chem Eng J 2008;137(1):72–83. [6] Almecija D, Blond D, Sader JE, Coleman JN, Boland JJ. Mechanical properties of individual electrospun polymer–nanotube composite nanofibers. Carbon 2009;47(9):2253–8. [7] Asran AS, Henning S, Michler GH. Polyvinyl alcohol–collagen–hydroxyapatite biocomposite nanofibrous scaffold: mimicking the key features of natural bone at the nanoscale level. Polymer 2010;51(4):868–76. [8] Unnikrishnan GU, Unnikrishnan VU, Reddy JN. Constitutive material modeling of cell: a micromechanics approach. J Biomech Eng 2007;129(3):315–23. [9] Baughman RH, Zakhidov AA, de Heer WA. Carbon nanotubes – the route toward applications. Science 2002;297:787–92. [10] Namilae S, Chandra N, Shet C. Mechanical behavior of functionalized nanotubes. Chem Phys Lett 2004;387(4–6):247–52. [11] Lee JH, Park TG, Park HS, Lee DS, Lee YK, Yoon SC, et al. Thermal and mechanical characteristics of poly(-lactic acid) nanocomposite scaffold. Biomaterials 2003;24(16):2773–8. [12] Tan EPS, Lim C. Characterization of bulk properties of polymer nanofibrous scaffolds from nanomechanical properties of single nanofibers. J Biomed Mater Res – Part A 2006;77A(3):526–33. [13] Vigolo B, Poulin P, Lucas M, Launois P, Bernier P. Improved structure and properties of single-wall carbon nanotube spun fibers. Appl Phys Lett 2002;81(7):1210–2. [14] Griebel M, Hamaekers J. Molecular dynamics simulations of the elastic moduli of polymer–carbon nanotube composites. Comput Methods Appl Mech Eng 2004;193(17–20):1773–88. [15] Mylvaganam K, Zhang L. Nanotube functionalization and polymer grafting: an ab-initio study. J Phys Chem B 2004;108:15009–12. [16] Unnikrishnan VU, Reddy JN. Characteristics of silicon doped carbon nanotube reinforced nanocomposites. Int J Multiscale Comput Eng 2005;3(4):437–50. [17] Che I, Cagn T, Goddard WAI. Studies of fullerenes and carbon nanotubes by an extended bond order potential. Nanotechnology 1999;10:263–8. [18] Lukkassen D, Persson LE, Wall P. Some engineering and mathematical aspects on the homogenization method. Compos Eng 1995;5(5):519–31. [19] Selmi A, Friebel C, Doghri I, Hassis H. Prediction of the elastic properties of single walled carbon nanotube reinforced polymers: A comparative study of several micromechanical models. Compos Sci Technol 2007;67(10):2071–84. [20] Unnikrishnan GU, Unnikrishnan VU, Reddy JN. Tissue–fluid interface analysis using biphasic finite element method. Comput Methods Biomech Biomed Eng 2008;3(1):1–8. [21] Thelen S, Barthelat F, Brinson LC. Mechanics considerations for microporous titanium as an orthopaedic implant Material. J Biomed Mater Res 2004;69A:601–10. [22] Leenslag J, Pennings AJ. High-strength fibers by hot drawing of as-polymerized poly(L)lactide. Polym Commun 1987;28(4):9294. [23] Pabst W, Gregorova E. Effective elastic properties of alumina-zirconia composite ceramics – Part 2. Micromechanical modeling. Ceramics – Silikaty 2004;48(1):14–23. [24] Reddy JN. Introduction to continuum mechanics: with applications. New York, NY: Cambridge University Press; 2008. [25] Reddy JN. Principles of continuum mechanics. New York, NY: Cambridge University Press; 2010. [26] Agoras M, Lopez-Pamies O, Ponte Castaneda P, Brun M. Homogenization estimates for fiber-reinforced elastomers with periodic microstructures. J Mech Phys Solids 2007;44(18–19):5953–79.
1014
V.U. Unnikrishnan et al. / Composite Structures 93 (2011) 1008–1014
[27] Bonet J. Nonlinear continuum mechanics for finite element analysis. 2nd ed. Cambridge: Cambridge University Press; 2008. [28] Guo Z, Caner F, Peng X, Moran B. On constitutive modelling of porous neoHookean composites. J Mech Phys Solids 2008;56(6):2338–57. [29] Holzapfel GA, Ogden RW. On planar biaxial tests for anisotropic nonlinearly elastic solids. A continuum mechanical framework. Math Mech Solids 2009;14(5):474–89. [30] Quapp KM, Weiss JA. Material characterization of human medial collateral ligament. J Biomech Eng 1998;120(6):757–63. [31] Peng XQ, Guo ZY, Moran B. An anisotropic hyperelasic constitutive model with fiber–matrix shear interaction for the human annulus fibrosus. J Appl Mech 2006;73(5):815–24.
[32] Qiu GY, Pence TJ. Remarks on the behavior of simple directionally reinforced incompressible nonlinearly elastic solids. J Elasticity 1997;49(1):1–30. [33] Mura T. Micromechanics of defects in solids. Hague, The Netherlands: Martinus Nijhoff; 1997. [34] Qian D, Wagner GJ, Liu WK, Yu M-F, Ruoff RS. Mechanics of carbon nanotubes. Appl Mech Rev 2002;55(6):495–533. [35] Zhang X, Liu T, Sreekumar TV, Kumar S, Hu X, Smith K. Gel spinning of PVA/ SWNT composite fiber. Polymer 2004;45(26):8801–7. [36] Xu X, Uddin AJ, Aoki K, Gotoh Y, Saito T, Yumura M. Fabrication of high strength PVA/SWCNT composite fibers by gel spinning. Carbon 2010;48(7):1977–84.
Contents lists available at ScienceDirect
Composite Structures journal homepage: www.elsevier.com/locate/compstruct
Multiscale analysis of carbon nanotube-reinforced nanofiber scaffolds V.U. Unnikrishnan, G.U. Unnikrishnan, J.N. Reddy * Advanced Computational Mechanics Laboratory, Department of Mechanical Engineering, Texas A&M University, College Station, TX 77843-3123, USA
a r t i c l e
i n f o
Article history: Available online 4 July 2010 Keywords: Carbon nanotube Nonlinear mathematical homogenization Hyperelastic material Multiscale analysis
a b s t r a c t Biopolymers play a significant role in tissue engineering, as they simulate the physiological environment required for the development of tissue cultures. Use of carbon nanotube polymeric scaffolds for tissue engineering applications has gained attention recently due to the enhanced mechanical properties of carbon nanotubes. In this paper, a hierarchical approach by studying the atomistic properties of carbon nanotube based polymers using molecular dynamics and coupling the scales through complex multiscale nonlinear hyperelastic material-based mathematical homogenization models are developed. These homogenization methods offer a systematic and rigorous treatment of up-scaling the properties from the micro or nanoscale to the macroscales. The material constitutive properties estimated using the developed methods for nanotube polymeric scaffold materials show excellent comparison with experimental studies. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction Carbon nanotube-reinforced nanofibers as scaffolds in bioreactor and other tissue engineering applications has received attention due to their enhanced material properties. The extremely high mechanical strength of carbon nanotubes is found to result in a nanofiber capable of withstanding tremendous mechanical stresses. The properties of these nanotubes can be obtained from various experimental investigations, but a thorough understanding of the properties and its interactions with the surrounding materials can only be obtained by advanced computational models. In this paper, we analyze the mechanical properties of carbon nanotube based Poly-Vinyl Alcohol (PVA) polymeric nanofibers [1–3], hierarchically from its constituents (i.e., PVA as matrix phase and nanotube as fiber phase). PVA is a water-soluble polyhydroxy polymer and has been widely used as a scaffold supporting material for tissue engineering applications, as it is bio-compatible and nontoxic and provides mechanical stability and flexibility [4–7]. The mechanical properties of atomistic scale nanotube system are estimated by molecular dynamic (MD) simulations. Bridging the spatial scales for a continuum nanotube approximation is carried out by volume averaging of the constituent properties of the atomistic system, using equivalence of the strain energy due to deformation of the continuum system and change in potential energy during an iso-thermal mechanical straining process of the atomistic system [8]. The information obtained from the volume averaging of the nanotube ensemble is then used in the evaluation * Corresponding author. Tel.: +1 979 862 2417. E-mail address: [email protected] (J.N. Reddy). 0263-8223/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruct.2010.06.019
of the effective macroscale properties of the nanofibers. Studies on the mechanical, electrical and thermal behavior of nanotubes in literature were primarily focused on atomistic simulations, however the computational adaptation to model macroscopic problems have not been completely established [9,10]. Continuum mechanics principles could be used to bridge the different scales of interest and provide solutions not only to the fine scale problems but also to the macroscale effects. This work is significant for the fact that, there have been very few studies on homogenization of nanotube based polymer scaffolds. Therefore the significance of this paper is that we intend to carry out both linear and nonlinear homogenization of nanotube based polymeric scaffolds. The overall effective properties of the nanofiber thus obtained can be used in various mathematical analysis, like simulation of nutrient transport within the nano-scaffold, etc. The mathematical representations of the behavior of a material, which are often validated against experimental results are known as constitutive relations of the material. It has been observed that the differences between the actual behavior and computational predictions are often due to the inaccurate representation of the constitutive behavior of the material. In this work, we consider a hyperelastic material model based on nonlinear homogenization to determine the constitutive properties of the nano-scaffold material. This is a significant improvement over conventional studies on nano-scaffolds as most of these studies were based on the assumption that the scaffold material is linearly elastic [11,12,8], and therefore accurate estimation of scaffold material properties is not possible. However, it has been found that the nano-scaffold materials exhibit nonlinear material behavior and therefore nonlinear homogenization is extremely important [13,1,6]. The significance of this
1009
V.U. Unnikrishnan et al. / Composite Structures 93 (2011) 1008–1014
work lies in the use of a multi-physical modeling strategy involving complex material geometry leading to an efficient multiscale–multi-physical computational model. This paper is organized as follows. Following this introduction, the analysis of nanotube-reinforced polymeric scaffolds is carried from the lower-order to the higher-order using MD simulations and scaling the properties by various mathematical homogenization methods. The second part of the paper deals with the nonlinear homogenization methods using hyperelastic material formulation, where the properties of the polymeric scaffold are obtained by dealing with the equivalence of the strain energy of the composite from the energy of the individual components. The effect of various fiber orientations on the overall effective properties of the nanotube–polymer composite system is also discussed. The paper concludes with a discussion of the result and conclusions.
Dt Dt fi ðtÞ ¼ vi t þ vi t þ Dt 2 2 mi
ð3Þ
where Dt is the time step, and fi(t) is the total force acting on particle i at time t. From the updated particle position, the interatomic forces are computed from the gradient of the potential energy field U with respect to the atomic coordinates ri as
fi ðt þ DtÞ ¼
@U @ri ðtþDtÞ
ð4Þ
The total potential energy due to straining, which can be considered as the elastic strain energy, can be expanded in Taylor’s series for small displacements, with the initial position represented by the equilibrium configuration. The elastic moduli tensor can then be written as (no sum on repeated subscripts)
3 ! 2 6 1 X 6 d U 1 dU a b c d 1 dU a c 7 7 aij aij aij aij þ dbd ¼ a a 7 6 2 2NXa j–i 4 dr ij r ij dr ij rij dr ij ij ij 5 |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} 2
2. Analysis of nanotube–polymer composite systems
C abcd
2.1. MD simulation of single-walled nanotube
Aac ¼0
Single-walled nanotubes are usually idealized as hollow cylinder formed from rolled graphene sheets, and each of the carbon atoms forming the tubules has three nearest neighboring bonds. In a fully relaxed structure, the angles between these bonds depend on the radius of the cylinder as well as on their orientation. The fundamental CNT structure can be classified into three categories according to their helicity: armchair, zigzag, and chiral. The knowledge of structure and molecular motion in nanoparticles are essential to understand properties of practical interest, and molecular dynamics (MD) has been widely used for the determination of these properties [14,15]. MD enables the computation of trajectories of atoms by integrating Newton’s second law of motion. The Hamiltonian, which represents the total energy of the system, is a function of position and momentum of the particles, and determines the trajectories of the particles [14,16]. For an isolated N-particle system, the Hamiltonian is given as the sum of the potential and kinetic energies
Hðr N ; pN Þ ¼
1 X 2 p þ UðrN Þ 2m i i
ð1Þ
where U(rN) is the potential energy from intermolecular interactions, which is represented as a function of the spatial ordinate P 2 1 r N ; 2m i pi is the kinetic energy, which represents the momentum pi of particle i with mass m; and pi is a function of the absolute temperature. The total potential energy (U) of the computational unit cell is given by the sum total of the valence bond energies ðV Bij Þ and the non-bonding interactions ðV NB ij Þ
U¼
X Xh j
V Bij þ V NB ij
i
ð5Þ
ð2Þ
j>i
In the atomistic level the molecular interactions are modeled using pair of potentials to respond to externally applied disturbance, and the reliability of MD simulation depends on the type of potential functions used. Tersoff–Brenner empirical bond order potentials are excellent for the simulation of carbon based materials and are used in the present work [15,14]. These bond order potentials are capable of describing the changes in the bonding between atoms, but does not describe long range interactions [17]. The force of attraction or repulsion experienced by each molecule is obtained from the gradient of the potential field which is used in calculating the updated position of the atoms using the velocity-Verlet time integration scheme. At the start of each time step, the updated velocities vi(t) are calculated for each particle by
Here Cabcd are the elastic moduli; U = U(rij) is the potential energy as a function of the inter-atomic distance rij; Aac is the internal stress tensor, which at equilibrium is equal to zero; Xa is the average volume of an atom; N is the number of atoms; dab is the Dirac delta function; a, b, c, d take on the values 1, 2, 3 which are the spatial dimensions and aaij is the undeformed value of rij in the a-direction. The deformed length is given by uaij ¼ raij aaij and any change in the inter-atomic distance (as a result of deformation) from the equilibrium (minimum energy position) inter-atomic distance contributes to the strain experienced by each atom and is given by ubij ¼ aaij eab , where eab are the components of the homogeneous infinitesimal strain tensor associated with atoms i and j. In a composite structure, where the nanotube is dispersed in a polymer matrix, the orientation and morphological characteristics like diameter, length and chirality also causes a corresponding variation in the macroscopic properties of the composite. In the next step, we estimate the effective properties of the polymeric nanocomposite using mathematical homogenization based on the elastic properties of the nanotube obtained from MD simulations and the polymer obtained from experimental studies. 2.2. Homogenization of nanotube – polymeric nanofibrous scaffolds There are numerous homogenization models for the analysis of structures composed of two or materials [18–20]. In mathematical homogenization models the heterogeneous composite material is replaced by an equivalent homogeneous continuum using a suitable homogenization technique. The homogeneous continuum is based on a statistically equivalent homogeneous volume element, called the representative volume element (RVE) which represents the material in the lower-order (the microscale) and collectively represents the global structure. 2.2.1. 3D Open Cell Homogenization based on structural models was developed for foams and cellular materials based on the ‘‘unit cell” model. These models are valid for materials with porosity 70% or higher and therefore can be used in the present analysis as numerical studies have shown that polymeric scaffolds have wide range of porosity [11,8]. The effective modulus of the nanofibrous scaffold by the structural model is given as [21]
C ¼ 0:6495Cð1 UÞ3
ð6Þ
where C* is the effective modulus of the porous material with elastic modulus C and porosity U. In a 3D open cell material model the
1010
V.U. Unnikrishnan et al. / Composite Structures 93 (2011) 1008–1014
effective modulus of the porous material is related to the fiber modulus by C* = C(1 U)2 [21]. This model is also based on assumptions of high porosity; and gives good predictions of moduli for materials with porosities from 10 to 90%. Nanofibrous scaffolds usually falls within this range of porosity and therefore this model can give a conservative prediction of the effective elastic modulus [22]. It should however be emphasized that the micromechanical based homogenization methods usually give rigorous bounds compared to structural models. 2.2.2. Dewey–Mackenzie estimates For non-interacting inclusions in a composite medium, exact solutions were obtained by Dewey and Mackenzie, called the Dewey–Mackenzie (DM) estimates [23]. The DM estimates are based on the exact classical solution of a spherical cavity embedded in an infinite elastic medium but the major drawback of the DM estimates is that they are dilute approximations. In the DM estimates, the elastic modulus is a function of Poisson’s ratio of the matrix and the porosity of the scaffolds given by
C DM
3ð1 mm Þð1 5mm Þ ¼ 1þ U Cm 2ð7 5mm Þ
ð9Þ
where r = rT is the Cauchy stress tensor, b is the body force per unit mass [24,25]. Along with the energy, entropy balance equations, and boundary conditions, the set of closure equations called the constitutive relationships are required for solving the above sets of equations. Hyperelastic materials can undergo large elastic strains with essentially small compressibility and the deformation characteristic is defined by the strain energy function W [26,27]. For a hyperelastic material there exists strain energy function (W), such that the stress components (rij) are derivable as rij ¼ @@eW , where eij are ij the strain components. The components of deformation gradient @xi tensor F are given by F ij ¼ @X , where Xj and xi are the components j of the material and spatial coordinates, X and x, respectively. The strain-energy density function can then be expressed as a linear combination of the strain invariants I1, I2, and I3.
W¼
n X
wk ðI1 ; I2 ; I3 Þ
ð10Þ
k¼1
ð7Þ
where CDM is the effective DM modulus, Cm is the modulus of the matrix and mm is Poisson’s ratio of the matrix and U is the porosity of the scaffold. 2.3. Nonlinear homogenization of hyperelastic materials Let us consider a composite system comprising the matrix and fibers phases and each material being homogeneous and hyperelastic. Let the fiber be oriented a degrees along the length of the composite RVE as shown in Fig. 1. The constitutive behavior of the individual phases can be characterized by the stored-energy functions which are assumed to be objective, isotropic and strongly elliptic (rank one convex) functions of the deformation gradient F. Now, considering X0 R3 to be the reference configuration of a deformable body with no discontinuities, let the current position of a point X be X 2 X0, which refers to the material position in the body. Under the action of external or internal forces, the point X moves to the position x in the current (deformed) configuration such that x = X + u, where u is the displacement vector and the motion (or deformation mapping) is defined by the mapping /(X, t). The principle of conservation of mass, assuming that there is no addition or deletion of mass, is given by
@q þ r ðqv Þ ¼ 0 @t
@ðq v Þ þ r ðqv v Þ ¼ r r þ qb @t
ð8Þ
where q is the density, v is the velocity. The principle of conservation of linear momentum gives the equations of motion
where wk are convex scalar functions and n is the number of individual phases in the hyperelastic composite system. In a composite material with two or more participating phases (e.g., matrix, fiber, interphase, etc.), the elastic response is assumed to result from the contributions of the individual phases and their interactions. e can be diTherefore, the strain energy function of the composite W vided into three parts, namely
e ¼ Wm þ Wf þ Wfm W
ð11Þ
where W represents the strain energy and the superscripts m, f, fm denotes the contributions due to the matrix phase, fiber phase and fiber–matrix interactions respectively. However, for non-interacting phases, we may neglect the contributions from the fiber–matrix interactions. The strain energy in terms of the principal invariants are then given by
W ¼ WðI1 ; I2 ; I3 ; I4 ; I5 Þ 1 I1 ¼ tr C; I2 ¼ ½ðtrCÞ2 tr C; 2 I4 ¼ a0 C a0 ¼ k2F
ð12Þ I3 ¼ detC ð13Þ
I5 ¼ a0 C2 a0 where kF is the fiber stretch and C = FTF is the right Cauchy–Green strain tensor, {k1, k2, k3} are the principal stretches, and a0 is the direction of the normal given by a0 = [cosa, sina, 0]T. For a hyperelastic material, the stress–strain relationship is obtained by taking the derivative of the strain-energy with respect to the Green–Lagrange strain tensor E ¼ 12 ðC IÞ, and is given by
S¼2
@W p C1 @C
ð14Þ
where S is the second Piola–Kirchhoff stress [24,25] and p is the pressure, which is added to account for the incompressibility constraint [28,29]. The corresponding Cauchy stress tensor ðr ¼ 1J FSFT Þ associated with this incompressible hyperelastic composite material is given by Quapp and Weiss [30] and Reddy [24]
h
i
r ¼ 2 W1 B þ W2 ðI1 B B2 Þ þ W4 a a þ W5 ða Ba þ Ba aÞ p I ð15Þ
Fig. 1. Nanotube fiber orientation in a fiber–matrix composite system.
where Wi = oW/@Ii, i = 1, 2, 4, & 5 and p is the pressure, I is the identity tensor, and a = Fa0 (see Fig. 1). Let us now consider the matrix material to be a compressible, Neo-Hookean material, and the strain energy in the reference configuration given in Eq. (16)
V.U. Unnikrishnan et al. / Composite Structures 93 (2011) 1008–1014
l
Wm ¼
ðI1 3Þ þ
2
j 2
ðJ 1Þ2
ð16Þ
where J = det(F). When the material is incompressible (J = 1) the second term in the strain energy density function vanishes. Correspondingly, the strain energy of the fiber due to the elongation may be considered as a function of I4, and assuming it to be in the form given in Eq. (17).
( F
W ¼
C 2 ðI4 1Þ2 þ C 3 ðI4 1Þ4
I4 > 1
0
I4 6 1
ð17Þ
It can be assumed that the fiber does not experience any strain during contraction of the composite [31]. Substituting I4 ¼ k2F and ¼ 12 ðk2F 1Þ [24] into Eq. (17), the strain energy of the fiber is obtained in terms of strain and can be used to fit experimental data for fiber elongation.
WF ¼ 42F C 2 þ 4C 3 2F
ð18Þ
The Cauchy stress tensor for this incompressible finite elastic composite material is then given by Qiu and Pence [32]
r ¼ 2F
@ W @I1 @ W @I4 T F p FC1 FT þ @I1 @C @I4 @C
ð19Þ
where
@W l ¼ @I1 2
and
@W ¼ 2C 2 ðI4 1Þ þ 4C 3 ðI4 1Þ3 @I4
1011
Using the above relation, the pressure which is now a function of the fiber orientation angle a,and principal stretches k1, k2 can be obtained and can be used to estimate the total stress due to deformation in the composite. 3. Results and discussion The elastic modulus of a single-walled nanotube (SWNT) is determined by a series of MD simulations, where the minimum energy configuration of the CNT is determined using minimization of the total potential energy (PE) at 0 K. The relaxed configuration of the nanotube is now subjected to displacement increments and the structure is allowed to equilibrate over a number of time steps. The straining process induces changes in the force experienced by each atom and this force is used in calculating the updated position of the atoms. The CNT is strained at a constant rate to calculate the strain energy of deformation by minimizing the total potential energy at each increment and the difference in the potential energy gives the strain energy due to deformation. From the MD simulations, the elastic modulus of a single-walled nanotube is obtained as 0.94 TPa and this value lies closely with many theoretically and experimentally obtained elastic modulus in the literature [16]. The elastic modulus of a polymeric electrospun nanofiber with diameter from 20 to 40 nm with 100% Poly-Vinyl Alcohol (PVA) is
ð20Þ
The derivatives of I1 and I4 are given by
@I1 @I4 ¼I ¼ a0 a0 @C @C
ð21Þ
Considering uniaxial elongation of a composite with fiber oriented along the direction of deformation, the deformation gradient F is given by
2
F1
0
6 F¼4 0
F2
0
0
0
3
7 05 F3
ð22Þ
where F1, F2 and F3 are correspondingly the principal stretches k1, k2, and k3, respectively. The corresponding stresses in the composite may then be written as:
rxx ¼ l 4C 2 ð1 k21 Þ 8C 3 ð1 þ 3k21 þ 3k41 k61 Þ k21 p ryy ¼ l k22 p; rzz ¼ l k23 p
ð23Þ
Fig. 2. Variation of elastic modulus of nanotube-reinforced PVA scaffold with void volume fractions.
It should be noted that due to incompressibility constraints, the principal stretches have the following properties: k1k2k3 = 1 and ki > 0, i = 1, 2, 3 [32] and traction free conditions exists in the y- and z-directions. The principal stretches in these directions for the fiber orientation angle a = 0 may be obtained as
1 k2 ¼ k3 ¼ pffiffiffiffiffi k1
ð24Þ
Due to the traction free conditions in the y and z directions, the pressure is obtained as
p ¼ lk1 1
ð25Þ
Additionally, we may also consider the fiber in the composite to be oriented at an arbitrary angle along the x direction. Now, assuming k3 = 1 (i.e., the deformation in the z direction is negligible), the principal stretches in the other two directions can be obtained as a result of incompressibility of the material and is given in Eq. (26).
k2 ¼
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ð3 k21 Þ 2
ð26Þ
Fig. 3. Stress–strain curve and its comparison between the experimental data and manufactured solution for a single-walled nanotube.
1012
V.U. Unnikrishnan et al. / Composite Structures 93 (2011) 1008–1014
Fig. 4. Variation of k2 with k1.
approximately 10 GPa. Thus, assuming a nanotube volume fraction of 0.005, the effective modulus of the nanofiber by Mori–Tanaka method is obtained as 10.90 GPa and by Voigt bound is 14.65 GPa. The details of Mori–Tanaka method and Voigt bounds may be found in standard literature and is not discussed here [33]. These values are very close to the experimentally obtained value 15 GPa by Wong et al. [3] and 16.24 GPa (see Fig. 5a of [6]). The difference in the predicted and experimentally obtained values could be due to the micro-structural formation during electrospinning, which could not be accounted in the mathematical homogenization models. The effective property of the scaffold is now estimated using the 3D-open cell model [8] and the DM estimates. Based on the homogenization studies from the previous sections, an average modulus of 10.9 GPa for the nanofiber and mm ¼ 0:4 is considered. The variation of elastic modulus of the scaffold with respect to various void volume fractions is shown in Fig. 2 and the inset shows the variation of the elastic modulus at the void limits. Now, considering the experimental investigation of elastic properties of nanotube bundles under tensile loads, the average
Fig. 5. q Nonlinear ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi stress–strain curve for tensile loading of nanotube-reinforced matrix composite for fiber various orientation angles for l = 3.33 GPa and 12.01 GPa and k2 ¼ þ
1 ð3 2
k21 Þ.
Fig. 6. Nonlinear stress–strain curve for tensile loading of nanotube-reinforced matrix composite for fiber various orientation angles for l = 3.33 GPa and 12.01 GPa and k2 = 2 k1.
V.U. Unnikrishnan et al. / Composite Structures 93 (2011) 1008–1014
1013
obtained stress–strain curve of nanotube systems and up-scaling to the macroscale properties. The effective property of the nanofiber scaffold estimated using both the micromechanical and hyperelastic material homogenization formulations were compared with experimental data and are found to be in good agreement. Acknowledgement The research reported herein was supported by the Oscar S. Wyatt Endowed Chair account at Texas A&M University. References
Fig. 7. Stress strain variation for uniaxial extension for l = 3.33 GPa and 12.01 GPa.
elastic moduli obtained by Yu et al. was 1.002 TPa [34] and the stress–strain curve reported by Qian et al. [34] is shown in Fig. 3. The experimentally obtained values are now fitted to the derived strain energy function in Eq. (18) to obtain the values of C2 and C3 and the curve-fitted stress–strain curve is denoted by the smooth line in Fig. 3. The elastic modulus of the matrix is found to vary from 10 GPa to 33 GPa [35,36] and considering the matrix to be incompressible, the shear modulus was found to range from 3.33 GPa to 12.01 GPa. Now, considering the orientation of the fiber along the length of the composite system, various stress–strain qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi profiles are obtained for k2 ¼ 12 ð3 k21 Þ (see Fig. 4) and is given in Fig. 5. However, various experimental studies have shown that k1 and k2 are often related by k2 = 2 k1 [31] (see Fig. 4) and the corresponding stress–strain variation is shown in Fig. 6, for various orientation angles of the fiber. It can be seen from the above figures that the stress–strain curve follows a nonlinear pattern and the difference in orientation angles severely affects the overall stress profile. The estimated stress strain curve is also compared with the experimental works of Xu et al. [36] and it can be seen that there is excellent correspondence with the experimentally determined data. Further, it can be seen that greater correlation is seen for higher shear modulus, which is the experimentally determined shear modulus for neat PVA fibers. It should also be noted that the nonlinearity is maximum when the fiber is oriented along the length of the composite RVE; this is because of the dominating nonlinear behavior of the fiber in the loaded direction. Stress–strain variation under perfect uniaxial extension, where the principal stretches are related as given in Eq. (24) (also, see Fig. 4) and the fiber oriented along the direction of deformation is also studied as shown in Fig. 7 for shear moduli equal to 3.33 GPa and 12.01 GPa. It can be seen that the nonlinear effects are significant at higher strain states. 4. Conclusions A multiscale modeling strategy to estimate the effective elastic modulus of carbon nanotube based polymeric nanofiber scaffolds is developed. The analysis is carried out from the atomistic scale comprising the nanotube using MD simulations, and the properties are scaled-up using various mathematical homogenization schemes. In this study, a nonlinear hyperelastic material based homogenization is also developed utilizing the experimentally
[1] Jeong JS, Moon JS, Jeon SY, Park JH, Alegaonkar PS, Yoo JB. Mechanical properties of electrospun PVA/MWNTs composite nanofibers. Thin Solid Films 2007;515(12):5136–41. [2] Wong KKH, Hutter JL, Zinke-Allmang M, Wan W. Physical properties of ion beam treated electrospun poly(vinyl alcohol) nanofibers. Eur Polym J 2009;45(5):1349–58. [3] Wong KKH, Zinke-Allmang M, Hutter JL, Hrapovic S, Luong JHT, Wan W. The effect of carbon nanotube aspect ratio and loading on the elastic modulus of electrospun poly(vinyl alcohol)–carbon nanotube hybrid fibers. Carbon 2009;47(11):2571–8. [4] Nho YC, Moon SW, Lee KH, Park CW, Suh TS, Jung YJ, et al. Evaluations of poly(vinyl alcohol) hydrogels cross-linked under Y-ray irradiation. J Ind Eng Chem 2005;11(1):159–64. [5] Mansur HS, Costa HS. Nanostructured poly(vinyl alcohol)/bioactive glass and poly(vinyl alcohol)/chitosan/bioactive glass hybrid scaffolds for biomedical applications. Chem Eng J 2008;137(1):72–83. [6] Almecija D, Blond D, Sader JE, Coleman JN, Boland JJ. Mechanical properties of individual electrospun polymer–nanotube composite nanofibers. Carbon 2009;47(9):2253–8. [7] Asran AS, Henning S, Michler GH. Polyvinyl alcohol–collagen–hydroxyapatite biocomposite nanofibrous scaffold: mimicking the key features of natural bone at the nanoscale level. Polymer 2010;51(4):868–76. [8] Unnikrishnan GU, Unnikrishnan VU, Reddy JN. Constitutive material modeling of cell: a micromechanics approach. J Biomech Eng 2007;129(3):315–23. [9] Baughman RH, Zakhidov AA, de Heer WA. Carbon nanotubes – the route toward applications. Science 2002;297:787–92. [10] Namilae S, Chandra N, Shet C. Mechanical behavior of functionalized nanotubes. Chem Phys Lett 2004;387(4–6):247–52. [11] Lee JH, Park TG, Park HS, Lee DS, Lee YK, Yoon SC, et al. Thermal and mechanical characteristics of poly(-lactic acid) nanocomposite scaffold. Biomaterials 2003;24(16):2773–8. [12] Tan EPS, Lim C. Characterization of bulk properties of polymer nanofibrous scaffolds from nanomechanical properties of single nanofibers. J Biomed Mater Res – Part A 2006;77A(3):526–33. [13] Vigolo B, Poulin P, Lucas M, Launois P, Bernier P. Improved structure and properties of single-wall carbon nanotube spun fibers. Appl Phys Lett 2002;81(7):1210–2. [14] Griebel M, Hamaekers J. Molecular dynamics simulations of the elastic moduli of polymer–carbon nanotube composites. Comput Methods Appl Mech Eng 2004;193(17–20):1773–88. [15] Mylvaganam K, Zhang L. Nanotube functionalization and polymer grafting: an ab-initio study. J Phys Chem B 2004;108:15009–12. [16] Unnikrishnan VU, Reddy JN. Characteristics of silicon doped carbon nanotube reinforced nanocomposites. Int J Multiscale Comput Eng 2005;3(4):437–50. [17] Che I, Cagn T, Goddard WAI. Studies of fullerenes and carbon nanotubes by an extended bond order potential. Nanotechnology 1999;10:263–8. [18] Lukkassen D, Persson LE, Wall P. Some engineering and mathematical aspects on the homogenization method. Compos Eng 1995;5(5):519–31. [19] Selmi A, Friebel C, Doghri I, Hassis H. Prediction of the elastic properties of single walled carbon nanotube reinforced polymers: A comparative study of several micromechanical models. Compos Sci Technol 2007;67(10):2071–84. [20] Unnikrishnan GU, Unnikrishnan VU, Reddy JN. Tissue–fluid interface analysis using biphasic finite element method. Comput Methods Biomech Biomed Eng 2008;3(1):1–8. [21] Thelen S, Barthelat F, Brinson LC. Mechanics considerations for microporous titanium as an orthopaedic implant Material. J Biomed Mater Res 2004;69A:601–10. [22] Leenslag J, Pennings AJ. High-strength fibers by hot drawing of as-polymerized poly(L)lactide. Polym Commun 1987;28(4):9294. [23] Pabst W, Gregorova E. Effective elastic properties of alumina-zirconia composite ceramics – Part 2. Micromechanical modeling. Ceramics – Silikaty 2004;48(1):14–23. [24] Reddy JN. Introduction to continuum mechanics: with applications. New York, NY: Cambridge University Press; 2008. [25] Reddy JN. Principles of continuum mechanics. New York, NY: Cambridge University Press; 2010. [26] Agoras M, Lopez-Pamies O, Ponte Castaneda P, Brun M. Homogenization estimates for fiber-reinforced elastomers with periodic microstructures. J Mech Phys Solids 2007;44(18–19):5953–79.
1014
V.U. Unnikrishnan et al. / Composite Structures 93 (2011) 1008–1014
[27] Bonet J. Nonlinear continuum mechanics for finite element analysis. 2nd ed. Cambridge: Cambridge University Press; 2008. [28] Guo Z, Caner F, Peng X, Moran B. On constitutive modelling of porous neoHookean composites. J Mech Phys Solids 2008;56(6):2338–57. [29] Holzapfel GA, Ogden RW. On planar biaxial tests for anisotropic nonlinearly elastic solids. A continuum mechanical framework. Math Mech Solids 2009;14(5):474–89. [30] Quapp KM, Weiss JA. Material characterization of human medial collateral ligament. J Biomech Eng 1998;120(6):757–63. [31] Peng XQ, Guo ZY, Moran B. An anisotropic hyperelasic constitutive model with fiber–matrix shear interaction for the human annulus fibrosus. J Appl Mech 2006;73(5):815–24.
[32] Qiu GY, Pence TJ. Remarks on the behavior of simple directionally reinforced incompressible nonlinearly elastic solids. J Elasticity 1997;49(1):1–30. [33] Mura T. Micromechanics of defects in solids. Hague, The Netherlands: Martinus Nijhoff; 1997. [34] Qian D, Wagner GJ, Liu WK, Yu M-F, Ruoff RS. Mechanics of carbon nanotubes. Appl Mech Rev 2002;55(6):495–533. [35] Zhang X, Liu T, Sreekumar TV, Kumar S, Hu X, Smith K. Gel spinning of PVA/ SWNT composite fiber. Polymer 2004;45(26):8801–7. [36] Xu X, Uddin AJ, Aoki K, Gotoh Y, Saito T, Yumura M. Fabrication of high strength PVA/SWCNT composite fibers by gel spinning. Carbon 2010;48(7):1977–84.