B. Arman, 1,2 Q. An, 1,3 S. N. Luo, 1,a) T. G. Desai,4 D. L. Tonks,1 T. Çağın,2 and W. A. Goddard III3
1Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA
2Department of Chemical Engineering, Texas A&M University, College Station, Texas 77845, USA
3Materials and Process Simulation Center, California Institute of Technology, Pasadena, California 91125, USA
4Advanced Cooling Technologies Inc., Lancaster, Pennsylvania 17601, USA
(Received 5 October 2010; accepted 8 November 2010; published online 4 January 2011)
We investigate with nonreactive molecular dynamics simulations the dynamic response of phenolic resin and its carbon-nanotube (CNT) composites to shock wave compression. For phenolic resin, our simulations yield shock states in agreement with experiments on similar polymers except the “phase change” observed in experiments, indicating that such phase change is chemical in nature. The elastic–plastic transition is characterized by shear stress relaxation and atomic-level slip, and phenolic resin shows strong strain hardening. Shock loading of the CNT-resin composites is applied parallel or perpendicular to the CNT axis, and the composites demonstrate anisotropy in wave propagation, yield and CNT deformation. The CNTs induce stress concentrations in the composites and may increase the yield strength. Our simulations suggest that the bulk shock response of the composites depends on the volume fraction, length ratio, impact cross-section, and geometry of the CNT components; the short CNTs in current simulations have insignificant effect on the bulk response of resin polymer. © 2011 American Institute of Physics. [doi:10.1063/1.3524559]
Polymers and polymer-based composites have long been explored/exploited for a wide range of engineering applications including high strain rate loading (e.g., shock waves).1–16 Despite extensive shock experiments on these materials,1–4,9–12 the underlying deformation and “phase change” mechanisms have been elusive due to the daunting complexities inherent in polymeric materials. While the challenge remains and numerical simulations of such materials are computationally intractable and expensive, reactive and nonreactive molecular dynamics (MD), coarse-grain dynamics and first-principles-based modeling/simulations are advantageous in revealing the microscopic details.6–8,13–16
Direct MD shock simulations of polymers and polymer composites are rare; some previous MD simulations explored the shock response of molecular crystals and chemistry.17 As a first attempt on direct MD simulations of shock response of polymers and polymer composites, we choose phenolic resin and its carbon-nanotube (CNT) composite in current study [Figs. 1(a) and 2]. CNTs are highly desirable as a structural component in the composites for their superior mechanical and physical properties.6–8,18 For instance, a recent MD work explored the shock response of the CNT-SiC composites modeled with the Tersoff potential.14 Our shock simulations yield the Hugoniots of phenolic resin and its CNT composites, and reveal the mechanisms for plasticity and the anisotropy in the shock response of the composites with regularly ordered CNTs. This work is presented as follows: Sec. II addresses the methodology of MD simulations and postprocessing, followed by results and discussion in Sec. III. Summary and conclusions are presented in Sec. IV.
The forcefield or interatomic potential describing the interactions in phenolic resin, CNT and their composites, is abinitio-based polymer consistent force field (PCFF).19 PCFF includes valence terms (bond, angle, torsion angle and out-of-plane angle) and nonbond interaction terms. The nonbond interaction terms account for electrostatic and van der Waals interactions. For the convenience of discussion, we denote C atoms in benzene rings as atom type C1, C in CH2 as C2, C in CNT as C1*, H in C–H as H1, and H in O–H as H2. This nonreactive forcefield is not appropriate for chemical reactions (involving bond breaking/formation), if any, induced by shock loading. Our MD simulations are performed with the large-scale atomic/molecular massively parallel simulator (LAMMPS) package.20 Periodic boundary conditions are applied along all three directions in nonshock simulations but
only along the two directions orthogonal to the shock direction in shock simulations.
We construct a unit configuration of amorphous phenolic resin with XENOVIEW (Ref. 21) via randomly placing 64 polymeric chains into a 2.7×2.9×10 nm3 supercell. Each such polymeric chain contains eight monomers [Figs. 1(a) and 2(a)]. A similar structure was studied with the reactive forcefield, ReaxFF.13 This unit configuration (7296 atoms) is equilibrated with the constant-pressure-temperature (NPT) ensemble and a time step of 0.25 fs, and then replicated by 2x2x8. The resulting configuration (233 472 atoms, or 5.4×5.7×80.6 nm3 in edge lengths) is further equilibrated at ambient conditions with a time (t) step of 0.25 fs and reaches a density of ρ<sub>0</sub>=1.12 g cm<sup>−3</sup>, and is adopted as the projectile for shock simulations. A larger configuration (2x2x12; 5.4×5.7×120.6 nm3) is also constructed and equilibrated for shock simulations.
A CNT composite unit configuration consists of a capped single-wall CNT with (10,0) chirality embedded in phenolic resin containing 64 polymeric chains, and is 2.9×2.9×10 nm3 in edge lengths. The CNT is ~0.78 nm in diameter and 7.8 nm long. The phenolic resin chains are introduced randomly around the CNT. The van der Waals distance between the resin and CNT is 0.34 nm, similar to an earlier work for CNT-polyethylene composite modeled with a Tersoff–Brenner potential and united atom model potential.22 This van der Waals distance thus induces the excluded volume between the CNT and the resin matrix, as seen in Fig. 2(b). The composite unit is first equilibrated at 0.1 K with the constant-volume-temperature (NVT) ensemble for 40 ps, followed by thermal annealing procedure at constant volume, where the system is heated with a ramp rate of 0.02 K fs−1 to 2000 K, equilibrated at 2000 K for 20 ps and then cooled to 300 K with the same ramp rate. The annealing procedure is repeated twice. Replications of the composite unit by 2x2x8 and 2x20x1 are adopted to construct configurations for longitudinal and transverse shock loading of the composites, respectively. The corresponding edge lengths are 5.7×5.9×84.1 nm3 (256 640 atoms) and 5.7×59.2×10.5 nm3 (320 800 atoms). [An example is shown in Fig. 2(b).] The resulting configurations (projectiles) are then equilibrated for 125 ps with the NPT ensemble at ambient conditions (ρ0=1.18 g cm−3) for shock simulations. The CNTs are tilted slightly as a result of relaxing these particular configurations.
The shock simulations adopt the projectile-wall geometry and microcanonical ensemble.23–25 A desired particle velocity along the shock or x-direction, up, is added to the x-component of the thermal velocities for each atom within the projectile. (The loading direction is along the direction with the longest dimension for a given supercell.) The other two directions orthogonal to the shock direction are y- and z-directions. Periodic boundary conditions are applied only along the y- and z-axes, and thus the nonimpact side of the projectile is a free surface. The bonds among the atoms on the impact and nonimpact surfaces are removed before simulation. The cell dimensions are fixed along the y- and z-directions, so the simulations mimic one-dimensional (1D) strain loading conditions as encountered in experiments. We choose wall/lj126 in LAMMPS as the wall. Upon impact, a shock wave is induced and propagates away from the wall into the projectile. For the CNT-composites, shock loading is applied either parallel or perpendicular to the CNTs, referred to as longitudinal and transverse loading, respectively. The time step for integrating the equation of motion is 0.5 fs, and run durations are up to 40 ps.
The atomic-level deformation can be characterized with the slip vector26
Here n is the number of the nearest neighbors to atom i, ns is the number of the slipped neighbors j, and xij and Xij denote the vector (between atom i and j) difference in current and reference configurations, respectively. The reference configurations are the preshock structures. Similarly, the maximum relative displacement is defined as27
The latter definition is used in our analysis, and the scalar slip is si=|si|. Another technique for characterizing shear deformation is the local von Mises shear strain28,29 but it is less revealing than si and thus not presented here.
We obtain the shock profiles of stress (σij; i , j=1, 2, and 3, or x, y, and z), temperature (T), density, and slip via 1D binning analysis.30 Pressure P follows as (1/3)(σ11+σ22+σ33), and the von Mises stress, 2τ=σ11−(1/2)(σ22+σ33), where τ is the maximum shear stress.
III. RESULTS AND DISCUSSION
Shock simulations are performed on pure phenolic resin and the CNT-resin composites along the longitudinal and transverse directions, and yield such profiles as stress, temperature, and slip as well as structure information. The stress (σ11) evolutions, plotted in the traditional x−t diagrams, illustrate wave propagation and interaction, which result in the shocked and unshocked regions as well as the release fan originated on the free surface (Fig. 3). We examine below the shock (Hugoniot) states, deformation, and related structural changes in pure resin and then the composites.
A. Shock response of phenolic resin
For phenolic resin, well supported shocks are observed in the x−t diagrams and such profiles as σ11(x) and T(x) [e.g., Figs. 3(a) and 4(a)]. For a given up, we obtain the Hugoniot state (denoted with a subscript H) values of σii,H, PH, and TH; the shock velocity us can thus be obtained from the jump condition as us=σ11,H/ Ρ0up (us can be measured directly from the shock fronts as in experiments3 but it may not represent the shock state given the complicated shock fronts.) The Hugoniot states are summarized in Figs. 5 and 6.
We observe in Fig. 5 a well-defined linear us−up relation for phenolic resin at up >1 km s−1: us=2.37 +1.58up (km s−1). However, the data points at up ≤1 km s−1 concave upward, lying below the extrapolation of the linear relation. Epoxy and methylol phenolic are two polymers similar to phenolic resin in their monomers and densities (Fig. 1); ρ0=1.192 g cm−3, 1.385 g cm−3, and 1.12 g cm−3, respectively.4 We thus compare the experiments on epoxy and methylol resin with our simulations of phenolic resin. Agreement is found approximately in the range of 1<up<3 km s−1; and the deviation from the linear extrapolation at the low up end appears to be common for all the three polymers. On the other hand, a phase change with a noticeable density increase is indicated at up3 km s−1 by the experiments but this feature is missing in the simulations.
As discussed by Carter and Marsh,4 the experimental us−up relations for a large number of polymers show three distinct regimes (I–III, with increasing up), schematically di-
vided by the two arrows in Fig. 5. Regime I shows a strong curvature, followed by regime II with a normal, linear us−up relation. As a result, extrapolation of regime II to zero up yields a us value above the ambient bulk sound speed. The detailed experimental study on polymethyl methacrylate (PMMA) is a solid example of this “general” observation.3 The curvature in Regime I can be explained with the interatomic potentials.4 Regime III is also linear, and considerable volume reduction occurs upon the II–III transition. They argued that this “phase transition” is neither polymorphic transformation in the usual crystallographic sense nor melting/vaporization, and that the breaking of covalent bonds within chains and subsequent reformation of tetravalent bonds between chains lead to large volume changes.4 The failure of our simulations to predict the II–III transition is consistent with its chemical nature inferred, since the bond breaking and formation are not allowed by the current forcefield. Therefore, reactive forcefields such as ReaxFF (Ref. 13) are necessary. The simulation results appear to be accurate up to up=3 km s−1. The shock temperature near the transition is about 1100 K (up=3 km s−1 and PH=23.4 GPa), lower than the value (~2000 K) estimated by Carter and Marsh.4
We calculate the radial distribution functions (RDFs) of phenolic resin in shocked and unshocked regions, and Fig. 7 shows the total RDFs and the corresponding coordination numbers (CN) for up=3 km s−1. Upon shock, the sharp peaks of the unshocked resin are smeared considerably; the average CN for the first neighbors is small and remains unchanged, while CN increases for the second shell and beyond. Since calculating si requires a sufficient number of nearest neighbors, we choose a cutoff distance of 2.5 Å. As an example, Fig. 8 shows a snapshot with color-coding based on si, which reveals clearly the shocked and unshocked regimes in phenolic resin (up=2 km s−1).
The dynamics of plastic deformation in shocked resin may be manifested in that of the von Mises stress 2τ. Upon shock arrival, 2τ rises rapidly to a peak value, 2τmax [as indicated by the arrow in Fig. 4(b)], and it then reaches a steady shock state value, 2τH. If 2τmax>2τH (shear stress relaxation), the shocked region undergoes plastic deformation. Figure 4(b) shows such stress relaxation due to plastic deformation, via the microscopic slip [Figs. 4(c) and 8]; and
the relaxation dynamics is nearly identical in 2τ(x) and s(x),
as expected [Figs. 4(b) and 4(c)]. In the case of plastic deformation, the shock front widths in σ11(x) [Fig. 4(a)] is much narrower than its counterparts in shear properties [Figs. 4(b) and 4(c)], similar to a shocked metallic glass.31
The elastic precursor is not definitely identified in our simulations, similar to experiments.1,3 The shock front [Fig.4(a)] shows a rapid rise followed by a rounding up to the shock plateau, a feature well documented in experiments.3 In metallic glass simulations, this rounding is related to plastic deformation.31 However, such rounding in phenolic resin occurs even at up=0.25 km s−1 (elastic shock; see below), likely due to viscoelastic behavior.3 (Both viscoelasticity and rate-dependent plasticity play a role at higher shock strengths.) Since there is no crystalline order in phenolic resin, there are no definitive structure features related to its plasticity as dislocations to crystal plasticity. The shear stress relaxation is a best indication of the elastic–plastic transition, and is absent at up−1 (thus presumably elastic). The well-defined values of 2τH increase with increasing shock strength, indicating strain hardening of the shocked resin [Fig. 6(b)]. (2τmax also increases with increasing shock strength.) Such strain or work hardening has been observed in both experiments and simulations of polymers.12,15
Atomic-level slip leads to the plastic deformation in phenolic resin; however, different types of atoms may differ in slip. Before shock, all atoms undergo slip solely due to thermal fluctuations, and s increases in the order of C1 (C2), O, H1, and H2, varying in the range of 1–2 Å (Fig. 9). Upon shock, s escalates to about 2.8 Å, 4 Å, and 5 Å for C, O, and H, respectively; and s is the same for H1 and H2 atoms, and C1 and C2 atoms (Fig. 9). The average preshock value of s is about 1.4 Å (Fig. 4). The atomic slip resistance increases in
the order of H, O, and C, and such differences may give rise to localized shear deformation, similar to the observation on a metallic glass.31 The backbone of a polymeric chain is composed of C atoms, and the stiffness of the backbone may enhance the slip resistance of C atoms. At longer range, the orientation of a segment in a chain may also affect its deformation.15 Therefore, the slip directions do not necessarily follow the presumed maximum shear stress directions (45°), as seen in Fig. 8. These structural inhomogeneities (intrachain and interchain) prevent the formation of longrange slips (so the slip deformation is localized, Fig. 8) and frustrate the plastic deformation, which may lead to strain hardening.
Plastic deformation in shock-loaded polymers has long been a subject of controversy,1,3,4 and the lack of elastic precursors certainly contributes to this debate. For example, Schmidt and Evans1 proposed that PMMA yields in a wide range of stresses (no definite shear strength), while Barker and Hollenbach3 argued that there is a definite yield point in PMMA, and the missing elastic precursor could be due to similar elastic and plastic waves velocities. In our simulations, the absence or presence of shear stress relaxations in
B. Shock response of the CNT-resin composites
Shock loading is applied to the CNT-resin composites at the same up as to the pure resin. For the particular nanocomposites explored, incorporating CNTs in phenolic resin does induce certain observable features in the mechanical behavior, and the composites show anisotropic response to shock loading in compression and shear (Figs. 3 and 10–13).
Due to the higher shock impedance of CNTs (mostly higher elastic constants and shock velocity) as compared to the resin matrix, stress concentrations are induced upon compression. Such stress concentrations are manifested as “stripes” in the x−t diagrams [Figs. 3(b) and 3(c); a stripe is indicated by an arrow]. The right-tilting stripes are due to shock enhancement by the CNTs at the internal interfaces (between the resin matrix and CNTs), and the left-tilting stripes, by the high-impedance reflecting wall. (This phenomenon is essentially reshock or double-shock.) The numbers of such stress concentrations match those of CNTs in the unshocked composites for both longitudinal and transverse loading. This reshock yields spatial fluctuations in stress at a given time, e.g., Fig. 10(a), as well as temporal fluctuations for a given position. These fluctuations are inherent in dynamic response of structured materials, and depend on the geometry of CNTs within the matrix.
Figure 10(a) compares three wave profiles where the shocks are initiated at the same position (the wall) and recorded at the same time (10 ps). The shock front for the longitudinal loading leads slightly that for the transverse loading of the composite as well as that for the pure resin. The shock velocities for the latter two are similar. A precur-
sor is also observed at the foot of the shock front for the longitudinal loading. The higher wave speed in CNTs and the CNT geometry directly lead to the differences in the shock fronts. The length ratios (the total length of CNTs to the cell length) along the shock direction for the longitudinal and transverse loading are about 0.74 and 0.26, respectively, so the CNT effect on shock velocity is more pronounced for the longitudinal loading. We obtain σ11,H via averaging the shocked regime, and σ11,H achieved in the composites is sightly higher than the pure resin. Thus, CNTs increase the compressional “stiffness” of the resin, although this increase is not pronounced given the small CNT fraction (and length ratios). The volume ratio of CNTs in the nanocomposites is about 9%. Nonshock simulations on CNT-polyethylene composites yielded similar results.8 (A recent shock simulation of CNT-SiC composites showed more pronounced effect of CNTs.14) Increasing the length and volume ratios should have a positive effect both on the shock front characteristics and shock states.
Besides the compressional stiffness, CNTs also increase the shear resistance of resin as a structure component in the composites. Figure 10(b) shows the profiles of 2τ(x) at selected time instants for pure resin and the composites, which reveal relaxation in 2τ behind the shock front to a steady state. Depending on where the shock front traverses (resin or CNTs), 2τmax can vary substantially for the composites. For up=2 km s−1, 2τmax increases from about 3 GPa in the pure resin to 4 GPa in CNT-regions; the shock front in the longitudinal loading is broader than the transverse loading due to different CNT alignment geometry [Fig. 10(b)]. (It is easier to increase the length ratio and thus the yield strength in the longitudinal loading). However, CNT appears to have diminished effect on the steady state shear strength, 2τH. (2τH is comparable for pure resin and the composites simulated here, regardless the loading direction.) We expect that increasing the volume ratio of CNTs may further improve the shear strengths at a shock front (onset of plasticity) and steady shock state.
In the CNT-resin composites, the plastic deformation is also manifested as shear stress relaxation and accompanied by slip. For transverse loading, the slip profiles are relatively smooth and the steady state slip, sH, is about 3.6 Å at up =2 km s−1, slightly lower than the pure resin 3.9 Å, while there are large fluctuations in s(x) for the longitudinal loading, and sH is about 5 Å. Such observations can be explained with the s(x) profiles of individual atom types in the composites [Fig. 11; more mobile O, H1, and H2 atoms are omitted for clarity]. The C1* atoms (CNT) slip less than C1 and C2 atoms in transverse loading but much more in the longitudinal loading. The s(x) profile of C1* atoms in the transverse case lacks the pronounced fluctuations seen in the longitudinal case. In the latter case, the slip peaks in C1* occur concurrently with those in C, H, and O atoms in the resin matrix. Thus, the CNT geometry directly affects the slip behavior of the composites.
The anisotropic deformation/damage of CNTs under shock loading is examined in more detail in Figs. 12 and 13 (using up=2 km s−1 as an example). For the longitudinal loading, the CNTs are distorted with slip and twisting as well as compression-related diameter changes; the regions near the tube ends undergo the most amounts of slip, leading to the slip peaks in s(x) [Fig. 11(b)]. For transverse loading, the most pronounced feature is the crushing of CNTs along the shock direction, and the tube ends are bulged relative to the rest of the tube [Figs. 12(b) and 13]. Thus, compression, shear and torsion are involved to different extents in the CNT deformation for both loading cases. The difference in the deformation for these two loading cases can be correlated with the CNT cross-section normal to the shock direction (impact cross-section): it is about 1 nm2 and 8.7 nm2, respectively, for the longitudinal and transverse loading. In the longitudinal loading, the small cross-section of a stiffer CNT embedded in a soft resin matrix leads to more pronounced deformation (particularly slip); the CNT ends have less constraint than its center portion and are more susceptible to slip [Fig. 12(a)], except some CNT caps. The cap itself appears more rigid overall in the transverse loading likely due to its geometry, while some highly slipped atoms are observed in the cap region (Fig. 13). A neighboring polymeric chain shows complicated slip (Fig. 13). At higher shock strengths, a shocked CNT is severely deformed/damaged beyond recognition.
We have characterized the shock response of phenolic resin and the CNT-resin composites under longitudinal and transverse loading. The simulated shock states of phenolic resin agree with the experiments but fail to predict the phase change observed in experiments, likely because such phase change involves bond breaking and formation. The plastic deformation in phenolic resin is achieved via atomic level slip accompanied by shear stress relaxation. Phenolic resin also shows strain hardening, which could be caused by the frustrated slip related to intrachain and interchain inhomogeneities. The CNT-resin composites demonstrate anisotropy in wave propagation, yield and CNT deformation/damage. The CNTs induce stress concentrations in the composites and may increase the yield strength. Our simulations suggest that the bulk shock response of the composites depends on the volume fraction, length ratio, impact cross-section, and geometry of the CNT components; the short CNTs in current simulations have insignificant effect on the bulk response of resin polymer.
We have benefited from J. W. Lawson, C. Wei (NASA Ames Research Center), J. Li (UPenn), and C. Brandl (LANL) in various ways. This work was partly supported by the Advanced Simulation and Computation (ASC) Program at LANL. LANL is operated by Los Alamos National Security, LLC for the U.S. Department of Energy under Contract No. DE-AC52-06NA25396. T.G.D. was supported by the NASA Small Business Innovation Research (SBIR) Grant under Contract No. NNX10CC69P.
1D. N. Schmidt and M. W. Evans, Nature London 206, 1348 1965.
2D. J. Pastine, J. Chem. Phys. 49, 3012 1968.
3L. M. Barker and R. E. Hollenbach, J. Appl. Phys. 41, 4208 1970.
4W. J. Carter and S. P. Marsh, “Hugoniot equation of state of polymers,”
Los Alamos National Laboratory Report No. LA-13006-MS, 1995.
5S. P. Marsh, LASL Shock Hugoniot Data University of California Press,
6K. Liao and S. Li, Appl. Phys. Lett. 79, 4225 2001.
7Y. Hu, I. Jang, and S. B. Sinnott, Compos. Sci. Technol. 63, 1663 2003.
8S. J. V. Frankland, V. M. Harik, G. M. Odegard, D. W. Brenner, and T. S.
Gates, Compos. Sci. Technol. 63, 1655 2003.
9N. K. Bourne and A. M. Milne, J. Appl. Phys. 95, 2379 2004.
10R. E. Setchell and M. U. Anderson, J. Appl. Phys. 97, 083518 2005.
11J. C. F. Millett, N. K. Bourne, Y. J. E. Meziere, R. Vignjevic, and A.
Lukyanov, Compos. Sci. Technol. 67, 3253 2007.
12J. C. F. Millet, N. K. Bourne, E. N. Brown, and G. T. Gray III, in Shock
Compression of Condensed Matter–2007, edited by M. Elert, M. D. Furnish,
R. Chau, N. Holmes, and J. Nguyen American Institute of Physics,
Melville, NY, 2007, Vol. CP955, p. 719.
13D.-E. Jiang, A. C. T. van Duin, W. A. Goddard III, and S. Dai, J. Phys.
Chem. A 113, 6891 2009.
14M. A. Makeev, S. Sundaresh, and D. Srivastava, J. Appl. Phys. 106,
15T. Ge and M. O. Robbins, J. Polym. Sci., Part B: Polym. Lett. 48, 1473
16T. R. Mattsson, J. M. D. Lane, K. R. Cochrane, M. P. Desjarlais, A. P.
Thompson, F. Pierce, and G. S. Grest, Phys. Rev. B 81, 054103 2010.
17M. L. Elert, S. V. Zybin, and C. T. White, J. Chem. Phys. 118, 9795
18R. H. Baughman, A. A. Zakhidov, and W. A. de Heer, Science 297, 787
19H. Sun, J. Phys. Chem. B 102, 7338 1998.
22C. Wei and D. Srivastava, Nano Lett. 4, 1949 2004.
23B. L. Holian, Shock Waves 5, 149 1995.
24S. N. Luo, L.-B. Han, Y. Xie, Q. An, L. Q. Zheng, and K. Xia, J. Appl.
Phys. 103, 093530 2008.
25L. B. Han, Q. An, S. N. Luo, and W. A. Goddard III, Mater. Lett. 64, 2230
26J. A. Zimmerman, C. L. Kelchner, P. A. Klein, J. C. Hamilton, and S. M.
Foiles, Phys. Rev. Lett. 87, 165507 2001.
27C. Brandl, Ph.D. thesis, Ecole Polytechnique Fédérale de Lausanne, 2009.
28J. Li, Modell. Simul. Mater. Sci. Eng. 11, 173 2003.
29F. Shimizu, S. Ogata, and J. Li, Mater. Trans. 48, 2923 2007.
30S. N. Luo, Q. An, T. C. Germann, and L. B. Han, J. Appl. Phys. 106,
31B. Arman, S.-N. Luo, T. C. Germann, and T. Çağın, Phys. Rev. B 81,
32A. Stukowski, Modell. Simul. Mater. Sci. Eng. 18, 015012 2010.