Sequential multiscale modelling of SiC / Al nanocomposites reinforced with WS 2 nanoparticles under static loading

E. I. Volkova,1,2 I. A. Jones,1,* R. Brooks,1 Y. Zhu,3 and E. Bichoutskaia2,* 1Division of Materials, Mechanics and Structures, Faculty of Engineering, University of Nottingham, Nottingham, NG7 2RD, United Kingdom 2School of Chemistry, University of Nottingham, Nottingham, NG7 2RD, United Kingdom 3College of Engineering, University of Exeter, Exeter, EX4 4QF, United Kingdom (Received 9 February 2012; revised manuscript received 21 July 2012; published 24 September 2012)


I. INTRODUCTION
Lightweight high performance ceramic/metal composites have recently attracted intense academic and industrial interest due to their high strength, ductility, and hardness, as well as the ability to withstand severe shock loadings. 1apid development of these composites, focused mainly on inclusion of ceramic nanoparticles such as B 4 C, SiC, TiB 2 , and Al 2 O 3 , offers a great potential for utilization in many critical protective applications. 2,3It has been confirmed both experimentally 4 and theoretically 5 that these composite materials can exhibit much higher strength and hardness than their parental bulk counterparts, not only under general ambient conditions but also under high shock loadings.Additional incorporation of inorganic fullerene-like (IF) nanoparticles into ceramic/metal matrices leads to a new form of multiphase nanocomposites 6 with potentially improved shock absorbing properties.Nanocomposites containing IF-WS 2 nanoparticles have already demonstrated outstanding tribological and wear properties. 7][10] Similar to IF-WS 2 nanotubes, IF-WS 2 fullerenes are hollow multilayered structures, but with approximately spherical shape.They can act as molecular absorbers in the SiC/Al nanocomposite, damping shock energy through the large interlayer separation (van der Waals gap) similar to that of the bulk interlayer distance of 0.62 nm.Thus a combination of a tough and strong SiC/Al matrix, which effectively stops fragment penetration and perforation, with energy-absorbing IF-WS 2 nanoparticles, could ultimately offer the next generation of antishock materials.
The high demand on imminent utilization of the advanced multiphased ceramic nanocomposites under severe loading conditions requires an improved understanding of the relationship between the atomic and macrostructure of these materials and its effect on the shockwave response.A detailed analysis of the influence of multiple phases on the mechanical and elastic properties could elucidate the advantages of these nanocomposites compared to the corresponding single-phase materials.The finite element modelling of hollow, faceted IF-WS 2 nanoparticles under compression was recently reported by Kalfon-Cohen et al. 11 In their work the mechanism of failure under compression was investigated using a modification of a model previously published for WS 2 nanotubes. 12he present work is motivated by a need to model the response of IF-WS 2 -containing nanocomposites to a variety of load cases, beginning with simple static loads and progressing via a range of dynamic situations, with a view to constrtucting a comprehensive explicit FE model of the material's response to shock and impact loads.It is also motivated by the need for a set of reliable and complete input data for such an analysis.In this paper, a first stage of this analysis is reported, focusing on the generation of reliable material property data for WS 2 and the prediction of elastic properties for the nanocomposite.The methodology employs a sequential multiscale modelling approach (as classified by, for example, Lu and Kaxiras 13 ) also referred to as hierarchical multiscale modeling 14 that combines atomic level investigation of the elastic properties of the WS 2 bulk material with the continuum level study of the aggregate behavior of an example SiC/Al nanocomposite impregnated with IF-WS 2 nanoparticles.This is in contrast to a concurrent multiscale approach where the models at different scales are coupled and run simultaneously, for example a coupled FE/MD simulation. 13The paper is organized in three further sections.In Sec.II the adopted methodology is described, and the results for the elastic properties of bulk WS 2 and MoS 2 (included as a test study) materials computed with density functional theory (DFT) including dispersion corrections are presented.In Sec.III, the results of the static loading finite element simulation of the multiphase nanocomposite are discussed and compared with the existing theoretical approach of Budiansky. 15In Sec.IV we analyze the obtained results and draw conclusions.

II. CALCULATION OF THE ELASTIC PROPERTIES AT THE ATOMIC LEVEL
As the diameter of the IF nanoparticles is approximately 100 nm 16 and therefore much larger than the interlayer spacing, the curvature of the particles has been ignored, and in the DFT calculations of the structure and the elastic properties, the WS 2 and MoS 2 materials have been treated as bulk structures.Three approximations for the exchange-correlation energy have been explored, namely the local density approximation (LDA/CA-PZ 17,18 ), and the Perdew-Burke-Ernzerhof (PBE 19 ) and the Perdew-Wang (PW91 20 ) parametrizations of the generalized gradient approximation (GGA).Van der Waals interactions have been included using the semiempirical dispersion corrections to the total DFT energies.The OBS dispersion scheme 21 has been used for the LDA and PW91 functionals, and the G06 scheme 22 has been used for the PBE functional as implemented in CASTEP 5.5 quantum chemistry code. 23A unit cell consisting of six atoms has been used as shown in Fig. 1(b).An on-the-fly pseudopotential generator has been used to eliminate the core states and describe the valence electrons by nodeless pseudo-wave functions.A plane wave basis set with the cutoff energy of 440 eV has been used, and the Monkhorst-Pack grid of 12 k points has been employed to sample the Brillouin zone.This approach has been initially used on the MoS 2 bulk material for which there is ample theoretical and experimental structural data [24][25][26] (see Table I).8][29] The authors judged it essential to perform a robust test of the ab initio predictions prior to calculating the WS 2 elastic properties for which no complete set of data appears to be available in the literature.
Inclusion of dispersion interactions in DFT calculations improves the prediction of the interlayer equilibrium geometry of the MoS 2 and WS 2 bulk, yielding a better agreement with experiment.This is particularly evident for the GGA functionals (see Tables I and II).For example, the optimized GGA/PW91 + OBS values of the lattice parameters for WS 2 are calculated to be a = b = 0.318 nm, c = 1.250 nm in good agreement with experiments. 27,28The W-W distance is found to be 0.61 nm in the c direction and 0.31 nm in the ab plane, which is also consistent with the values reported in the literature. 27The GGA + OBS values of the layer thickness 2cz and the lattice constant c are closer to experimental values 27,28 than previously reported GGA evaluations 29 without dispersion correction.Similarly good agreement is seen in Table I for MoS 2 , where experimental values of the vdW gap and the Mo-S bond length are additionally available for comparison with predictions.
Having obtained the optimized structures of MoS 2 and WS 2 the full elastic constant tensors have been calculated for both structures.According to Hooke's law, 30 the relationship between stress and strain can be linearly expressed as where C ij is the stiffness matrix.
The elastic constants, in the form of the stiffness matrix, are obtained within the CASTEP package using the finite strain method described by Milman and Warren 31 and in TABLE I. Equilibrium geometry (in nm) of the MoS 2 bulk material: a and c are the lattice parameters of the unit cell, 2cz is the thickness of the layer as shown in Fig. 1(b), and d the CASTEP theory documentation. 32Prescribed strains are applied to the unit cell, the structure is optimized for each deformed state, and the stresses for each strain state are calculated.Only small values of strains were applied in order to remain within the elastic region of the compounds.The elastic stiffness coefficients C ij within a general definition of Hooke's law are then obtained by fitting the stresses as a linear function of strains.The fitting procedure is necessary because it is not possible to calculate elastic properties of the compounds from applying just one value of strain. 33In the present work, a strain of 0.007 was applied in nine incremental steps, and it was found to give a good compromise between nonlinearity and numerical errors.Since MoS 2 and WS 2 bulk materials form hexagonal crystals with a layered structure, it is assumed that these materials are transversely isotropic (in common with other materials with hexagonal structure), 34 so that they can be characterized by a set of five independent elastic constants. 35,36istinct, nonzero terms from the stiffness matrices for MoS 2 and WS 2 are presented in Tables III and IV It is straightforward to show that the single-crystal bulk modulus, K trans , of each of these transversely isotropic compounds can be calculated from the elastic constants, namely Young's moduli, E i , and Poisson's ratios, ν ij , as follows: where 37 The computed GGA values for K trans of transversely isotropic structures are very close to the experimental values 38,39 for the MoS 2 38 and WS 2 39 bulk structures obtained by fitting the pressure-volume data (derived from lattice parameter measurements under pressure loading) to the thirdorder Birch-Murnaghan equation of state (see Tables III and  IV).The LDA/CA-PZ + OBS approach, however, overestimates K trans by about 20%.The values of C 11 , C 33 , and C 44 for the MoS 2 bulk, calculated using GGA in conjunction with PBE and PW91 functionals with empirical dispersion correction, agree well (within 10% or less) with Feldman's data 40 though agreement is poorer for C 13 and irreconcilable with Feldman's surprising negative inferred value for C 12 .The LDA C 33 values are in poorer agreement with all other values including a large mismatch (up to 47%) with the experimental data of Sourisseau. 41In general, the LDA results with OBS dispersion scheme both for lattice parameters and for the elastic properties of MoS 2 are in the poorest agreement with the experimental data.The components of the elastic matrix calculated for the WS 2 bulk, using both LDA and PW91 functionals including empirical dispersion correction, are in disagreement with Sourisseau's experimental data. 41The values of C 11 are at least 60% different.PW91 with OBS dispersion scheme, however,  gives the lowest energy for both MoS 2 and WS 2 bulk materials and therefore the GGA/PW91 + OBS elastic properties are taken forward for use within the continuum finite element model in the form of the stiffness matrix components C ij .

III. ELASTIC PROPERTY CALCULATION FOR NANOCOMPOSITE AT THE NANOPARTICLE SCALE A. Finite element model of the nanocomposite
As a first stage towards understanding the behavior of the nanocomposite under a range of load cases, the elastic properties obtained from DFT calculations have been subsequently used in simulations of the nanocomposite using the finite element method (FEM). 42The FEM involves discretizing a structure into a large number of regions or elements, each associated with an appropriate material model, and the deformation of the structure is defined in terms of displacements of the nodes to which each element is linked.The relationship between nodal displacements and forces takes the form of a stiffness matrix calculated from the mechanical properties and geometry of each element.In the present (quasistatic) problem, an implicit solution method is used in which the unknown quantities (unrestrained displacements and reaction forces at restrained nodes) are found using linear algebra techniques.Here, the prediction of the nanocomposite's behavior has been performed using a unit cell model established within the ABAQUS finite element system, 43 which is a commercially available suite of software for solving engineering mechanics problems using the techniques described above.The true structure will have randomly distributed particles of WS 2 and SiC incorporated in an Al matrix with no particular directionality.As a first approximation each WS 2 particle is assumed to have a perfect hollow spherical shape having transversely isotropic properties such that the local x and y directions are tangential to the spherical surface and the local z direction is radial.In order to represent such a material via a unit cell model of a manageable size, a regular structure has been adopted as a first approximation which for convenience initially assumes an interpenetrating simple cubic configuration [Fig.1(c)] identical to caesium chloride crystal structure.The mechanical properties, particle dimensions, and volume fractions used to generate the ABAQUS models are presented in Table V.The WS 2 properties were defined in the model using a local spherical coordinate system.
The bounding planes of the unit cell are defined as follows: where x, y and z are the Cartesian coordinates of the unit cell, and h 1 , h 2 and h 3 are respectively its dimensions in the three Cartesian directions.In the present case h 1 = h 2 = h 3 = 100 nm.The planar and rotational symmetries of the unit cell allow a simplified form of periodic boundary conditions, which, for the application of direct strains, are implemented as follows [Fig.2(a)]: where u, υ, and ω are uniform displacements applied to the opposing faces perpendicular to the three Cartesian axes, and u 0 is a displacement value chosen to achieve an appropriate level of strain, in the present case 0.001 or 0.1%.In a similar manner, the following boundary conditions were used for application of shear strains [Fig.2(b)]: Figure 3 depicts the distorted shape of the unit cell, with the displacements exaggerated for clarity and with the stresses shown as color contours.Noting that the configuration of the unit cell is symmetric over the three Cartesian directions, the components of the stiffness matrix for the nanocomposite were recovered from the reactions at the unit cell boundaries using the following relations: where ε x = 2u 0 h 1 is the extensional strain in the x direction and γ xy = 2u 0 h 2 is engineering shear strain in the xy plane.RF x , RF y , and RF z are reaction forces occurring in the x, y, and z directions, respectively.A x , A y , and A z are the projected areas (taken to be equal in the present case) of the unit cell normal to the x, y, and z directions respectively so that A x = h y h z etc. Mesh convergence was demonstrated by confirming that C 11 changed by − 0.1% as the number of elements was approximately doubled from 98804 to 216215.

B. Results of the finite element simulations
The stiffness and compliance components of the nanocomposite recovered by this method are It should be noted that the assumed particle configuration is an idealization, which acts as a first approximation to what is in reality a random structure with no preferential directions.Therefore, the model would not be expected to produce the isotropic properties of the true random material, and the predicted properties are those of a cubic material which needs three independent elastic constants to be defined.Specifically, its elastic behavior is uniquely defined by any three of the follows constants, which can be calculated from the cubic stiffness matrix: 37 The model yields the following values: K = 100.34GPa, G = 48.49GPa, E = 149.90GPa, and ν = 0.2510.The realistic material with a random structure can then be considered to consist of an isotropic aggregate with a locally cubic structure.Voigt and Reuss 44,45 approximations give the theoretical maximum and the minimum values of the average isotropic elastic moduli, respectively.The Voigt approximation assumes that the uniform strain in the compounds is equal to external strain and the Reuss approximation assumes that the uniform stress in the compounds is equal to external stress.The Reuss bulk modulus K R and the Voigt bulk modulus K V are equal for cubic materials 46 and are given by and the Reuss shear modulus G R and the Voigt shear modulus G V are given by Hill 46 demonstrated that the approximations of Voigt and Reuss give upper and lower bounds of the isotropic moduli, respectively, and recommended that the realistic approximation of the isotropic moduli of isotropic aggregates is the arithmetic mean of these limits.Hence the elastic moduli, K iso and G iso of the isotropic material can be approximated as Since these two constants now fully define the isotropic material, the values of Poisson's ratio ν iso and Young's modulus E iso are obtained by substituting the calculated values of K iso and G iso into The obtained values are K iso = 100.34GPa, G iso = 52.79GPa, E iso = 134.74GPa, and ν iso = 0.2762.The authors are not aware of an analytical solution which fully models a particulate with hollow, spherical, transversely isotropic inclusions.However, Budiansky's method, 15 in turn based upon Eshelby's inclusion technique, 47 finds the aggregate elastic constants of an inhomogeneous material with multiple isotropic phases taking the form of ellipsoidal or (as a special case) spherical inclusions embedded in a matrix phase; despite the assumption of inclusion shape, the solution reduces to a form symmetric across all phases including the matrix.Budiansky's solution is used here to obtain two estimates of the overall modulus of the nanocomposite.Both estimates treat the hollow interior of the nanoparticles as a fourth phase of zero modulus and volume fraction 2.4% and use the moduli for the SiC and the Al phases directly from Table V.They also both treat the nanoparticles as isotropic and with a volume fraction equal to the true volume fraction of their solid component (17.1%).It is also assumed that their Poisson's ratio is 0.22, equal to ν 12 for the WS 2 from Table V.However, one estimate assumes that the Young's modulus for the nanoparticles takes a value of 224 GPa numerically equal to E 1 = E 2 , while the other assumes a value of Young's modulus of 60 GPa which gives the same bulk modulus as for the transversely isotropic WS 2 in Table IV.The results of this comparison are presented in Table VI.While no attempt has been made to explore whether these estimates are indeed bounds on the true value of the nanocomposite modulus, it is reassuring that the estimate of average isotropic E iso calculated from the FE model lies close to the one of the estimates of Budiansky with E(WS 2 ) = 60 GPa, while the FEM estimate of Poisson's ratio is within 7% of the Budiansky estimates.It is seen in Fig. 3(a), and from closer examination of the FEM output, that the highest stresses in the direct loading situation occur as circumferential direct stresses in the loading direction, at the surface of the nanoparticles.It is noteworthy that the transverse and shear flexibility of the WS 2 prevents transmission of significant load to the inner layers.Similarly, the shear loads are carried as shear stresses at the surfaces of the nanoparticles [Fig.3(b)].The stresses in the ceramic reinforcement, and in the metallic matrix, are considerably lower than in the WS 2 nanoparticles, and it would appear that the load carrying mechanism of the nanoparticles acts to avoid stress concentrations in either of these two phases at the expense of causing a relatively large stress concentration in the nanoparticle (Fig. 4).While the response of WS 2 particles to Hertzian crushing loads has been explored elsewhere, 49 we are not aware of work characterizing their response to the present load cases and with present application.
The case study presented here acts as a feasibility study in which the modelling of IF particles with local spherically orthotropic properties has been achieved, and the results of the analysis demonstrated to be consistent with alternative approaches.This lays the foundations for undertaking more sophisticated analyses, involving larger numbers of particles and more specialized (explicit) solution methods, in order to model the response of the nanocomposite to more challenging load cases such as shock and impact.In order to model the dynamic and failure response of the composite under high rate loading conditions, it will also be necessary to consider the strength of the IF particles and of the interfaces between particles and matrix material.

IV. SUMMARY
In conclusion, the present work has provided a rigorous theoretical prediction of MoS 2 and WS 2 elastic properties.For MoS 2 where more complete existing data are available, excellent agreement has been observed, giving confidence in the predictions for WS 2 .The structural geometric data for both compounds agree well with the existing published data.The present work also incorporates the predicted elastic properties for WS 2 into the first mechanical analysis of WS 2 /ceramic nanocomposites.This work both explores the load carrying mechanisms of the IF particles, and serves as a feasibility study for testing the mesh generation, boundary condition, and property definition procedures required in more complex analyses.Future work will concentrate upon extending the mechanical analysis to consider more practical particle architectures and load cases, and to incorporate failure predictions via interface properties derived from MD calculations.

FIG. 1 .
FIG. 1. (Color online) Sequential multiscale approach: (a) TEM lattice image of IF particle showing layered structure as adopted from; 16 (b) a unit cell of the bulk material used in DFT: tungsten (molybdenum) atoms are denoted in green and sulfur atoms are denoted in yellow; (c) interpenetrating simple cubic unit cell of the nanocomposite used in the mechanical FE model: a SiC nanoparticle is in the center of the unit cell; IF-WS 2 nanoparticles are in the corners of the Al matrix.
, respectively.Our assertion that the materials are transversely isotropic is supported by noting that C 11 = C 22, C 13 = C 23, C 44 = C 55 , and C 66 = C 11 −C 12 2 to an accuracy of at least 1%.Both materials show a high degree of elastic anisotropy, with the highest stiffness constants being C 11 = C 22 along the a axis and b axis, respectively, where deformation involves bond bending and bond stretching.The stiffness constant C 33 is significantly lower (approximately by a factor of 5) because it involves weak interlayer van der Waals forces.

FIG. 2 .
FIG. 2. Periodic boundary conditions for the FE model.The uniform displacements are (a) applied to the left and right surfaces for the FE simulations under direct stresses and (b) to the top and the bottom surfaces for obtaining shear components.By contrast, the remaining surfaces are fixed against any movements.

5 (C 11 −
FIG. 3. (Color online) Three-dimensional distorted view of the FE unit cell under (a) direct stresses and (b) shear stresses.Deformation scale factor is 100.A uniform strain of 0.1% has been applied in both cases.

TABLE II .
(Mo-S) is the length of molybdenum-sulfur bond.Equilibrium geometry (in nm) of the WS 2 bulk material: a, c, 2cz and d are defined in TableI.

TABLE III .
Stiffness matrix, C ij (in GPa), and bulk modulus, K trans (in GPa), of the MoS 2 bulk obtained using Eq.(2).Duplicate and zero elements of C ij are omitted.

TABLE IV .
Stiffness matrix, C ij (in GPa), and bulk modulus, K trans (in GPa), of the WS 2 bulk obtained using Eq.(2).Duplicate and zero elements of C ij are omitted.

TABLE V .
Mechanical properties of materials used for the FEM analysis.The elastic properties of WS 2 are given here for comparison, and are an alternative representation of the stiffness matrix components presented in Table IV for PW91 + OBS, which were used directly as the FE input data.

TABLE VI .
15mparison of elastic properties of the isotropic composite from the DFT/FEM approach and Budiansky analytical model.15