Molecular Dynamics Simulation of Ceramic Matrix Composites Using BIOVIA Materials Studio, LAMMPS, and GROMACS

CHAPTER 5

Molecular Dynamics Simulation of Ceramic Matrix Composites Using BIOVIA Materials Studio, LAMMPS, and GROMACS This chapter has been divided into three parts as follows: Chapter 5.1: Molecular Dynamics Simulation of Carbon Nanotube-Reinforced Silicon Carbide Composites Using BIOVIA Materials Studio Chapter 5.2: Molecular Dynamics Simulation of Al/Al2O3 Metal-Ceramic Composite Using LAMMPS Chapter 5.3: Molecular Dynamics Simulation of Coaxial Boron Nitride/Carbon Nanotubes Using GROMACS

Chapter 5.1

Molecular Dynamics Simulation of Carbon NanotubeReinforced Silicon Carbide Composites Using BIOVIA Materials Studio Sumit Sharma, Pramod Kumar, Rakesh Chandra Department of Mechanical Engineering, Dr. B. R. Ambedkar National Institute of Technology, Jalandhar, India

5.1.1 Introduction Silicon carbide (SiC) composites are widely used in aerospace applications such as cowling mounts, reflectors, and combustion turbines [1]. This can be attributed to their excellent mechanical properties, high chemical stability, high thermal conductivity, and machinability. Molecular Dynamics Simulation of Nanocomposites using BIOVIA Materials Studio, Lammps and Gromacs https://doi.org/10.1016/B978-0-12-816954-4.00005-X # 2019 Elsevier Inc. All rights reserved.

227

228 Chapter 5 But their main disadvantage is their brittleness that has limited their usage in various fields. Thus, for improving the strength and toughness properties of SiC composites, various methods have been used such as the removal of residual Si and the use of second-phase particles (fibers or whiskers) as reinforcement. Several experimental works have recently confirmed the theoretically predicted outstanding mechanical properties of carbon nanotubes (CNTs). Consequently, CNTs emerge as potentially attractive materials as reinforcing elements in ceramic matrix composites. Peigney et al. [2] prepared novel CNTs, metal-ceramic nanocomposite powders, and dense materials. An investigation of the microstructural and mechanical properties of these materials was performed. The elaboration conditions of CNT-Fe-Al2O3 and CNT-M-MgAl2O4 (M: Fe, Co, Ni, and their alloys) nanocomposite powders were optimized with the aim to obtain both a great quantity of CNTs (up to 300,000 km of CNT bundles in a gram of powder) and a good quality of carbon, that is, a smaller average tube diameter or more carbon in tubular form. Because these two parameters generally vary in opposite directions, a compromise must be found. For reinforcement applications, with CNT-Fe-Al2O3 powders, it was inferred that a good compromise was obtained with 18% CH4 in the H2-CH4 gas mixture, 10 wt.% Fe, and a reducing thermal treatment performed at 1050°C with a short dwell time (6 min). Li et al. [3] prepared reaction-bonded SiCs (RBSCs) by uniformly dispersing chopped carbon fibers into bimodal SiC/C suspension. The main advantage of this method was the low cost and the mass production capabilities for fabricating large-scale and complex-shaped components (especially, RBSC mirrors) with improved mechanical properties. The effect of the chopped fiber fraction on microstructure and mechanical properties was evaluated. The maximum values of flexural strength and fracture toughness were 416 MPa and 5.1 MPa m1/2, respectively, superior to those of the monolithic RBSC formed via slip casting. The optimum fraction ratio of Cf/SiC was 30% when coarse and fine SiC particles are 60 and 10 mm in size, respectively. It was found that the interaction of fiber debonding, crack deflection, and crack bridging consumes crack energy during propagation and leads to the improved toughness. The chopped fibers had reacted with liquid silicon during reaction sintering, so little fiber pullout was observed. RBSC ceramic was developed first as a nuclear fuel cladding material in the 1950s. As a promising material, it then quickly found its applications in aeronautics, energy, electronics, nuclear, and transportation industries. Li et al. [4] developed SiC whisker-reinforced reactionbonded SiC composites for fabricating large-scale, complex-shaped structural components. The composites were prepared with the method of slip casting and liquid silicon infiltration at 1650°C. The distribution, morphology, and reinforcing behaviors of the SiC whisker in the composite were investigated. It was revealed that the introduction of SiC whisker increases the porosity of the green body and accordingly the bulk density of the composite. Whisker pullout

Molecular Dynamics Simulation of Ceramic Matrix Composites 229 was observed on the fractured surface, implying a moderate bonding strength between the whisker and matrix. After LSI the SiC whisker kept its initial diameter and morphology in the case of 15 wt.% whisker. The fracture toughness was enhanced by SiC whisker, reaching the peak value of 4.2 MPa m1/2 at the whisker fraction of 20 wt.%. Whisker pullout, whisker bridging, and crack deflection were considered as the main toughening mechanisms. A new nonaqueous gelcasting system of phenolic resin and furfuryl alcohol combined with a curing catalyst was developed for casting of RBSCs by Zhou et al. [5]. This gelling system could be carried out in air, and the surface exfoliation phenomenon that seems inherent to the acrylamide gelcasting system could also be eliminated. Polymerization of the premix solutions and rheological properties of the nonaqueous SiC suspensions were studied. After curing and subsequent pyrolysis of the concentrated SiC suspension, homogeneous SiC/carbon green body with a relatively high strength of about 18 MPa could be formed. It was found that dense complex-shaped SiC ceramic parts with flexure strength of 300
20 MPa and fracture toughness of 3.87
0.19 MPa m1/2 could be successfully produced after reaction sintering at 1700°C for 30 min under vacuum. Mei et al. [6] prepared vertically aligned CNT-reinforced SiC composites (VACNT/SiC) that were prepared via a combined chemical vapor deposition/chemical vapor infiltration process. Thermal oxidation tests showed the material to undergo virtually no mass loss up to 1400°C in air, hence indicating excellent oxidation resistance of the composites within the target temperature range of advanced ceramics applications. Microstructural observation demonstrated complete absence of matrix cracks in the surface of VACNT/SiC. A large number of VACNTs having undergone pullout were observed, a finding that indicated the tubes’ survival after exposure to the oxidative environment of 1400°C in air. Following thermal oxidation testing, only the surface SiC of the composites was found oxidized to silica, whereas the embedded MWCNT forest structure retained original morphology and integrity, exhibiting complete absence of oxidation due to the achievement of a crack-free SiC matrix providing effective isolation from high-temperature oxidative environments. Film formed by CNTs is usually called CNT film (CNTf). In the study conducted by Mei et al. [7], CNTf fabricated by floating catalyst method was used to prepare CNTf/SiC ceramic matrix composites by CVI. Mechanical and electric properties of the resulting CNTf/SiC composites with different CVI cycles were investigated and discussed, and the results revealed that the CNTf had a good adaptability to CVI method. Tensile test demonstrated an excellent mechanical performance of the composites with highest tensile strength of 646 MPa after two CVI cycles, and the strength has a decline after three CVI cycles for an excessively dense matrix, while the elastic modulus of the composite increased with the CVI cycles and reached 301 GPa after three CVI cycles. Tensile fracture morphologies of the composites were analyzed by scanning electron microscope to study the performance change laws with the CVI cycles. With SiC ceramic matrix infiltrated into the CNTf,

230 Chapter 5 enhanced electric conductivity of the CNTf/SiC composite compared with pure CNTf was also obtained, from 368 to 588 S/cm. Conductivity of the SiC matrix with free carbon forming in the CVI process was considered as the reason. Reaction sintering is one of the most attractive manufacturing processes of SiC, because of dense structure, low processing temperature, good shape capability, low cost, and high purity. However, mechanical properties of reaction-sintered SiC (RS-SiC) were typically much lower than normal sintered one. Particularly the bending strength was approximately 300 MPa. Suyama et al. [8] predicted the effect of the microstructure on the bending strength of RS-SiC. The bending strength of RS-SiC was recognized to be increased with decreasing the residual silicon (Si) size. The strength of over 1000 MPa was obtained to control the residual Si size under 100 nm. When the size of the residual Si was made small from 1.4 to 0.1 μm, the average bending strength increased from 580 to 1070 MPa. The high-strength RS-SiC showed high Young’s modulus and hardness that were equivalent to normal sintered one and high thermal conductivity that was equivalent to commercial RS-SiC. Miyamoto et al. [9] studied the effect of SiC coating on the properties of MWCNTs. These tubes were coated with a nanometer-sized SiC polycrystalline layer by the reaction of SiO(g) and CO(g). Nanometer-sized SiC powders were mixed with 1–5 vol% SiC-coated MWCNTs and successfully sintered at 1800°C by means of pulsed electric current sintering. These composites showed superior microhardness of 30.6 GPa and toughness of 5.4 MPa m1/2 compared with monolithic SiC ceramics due to the improvement of adhesion between MWCNTs and the SiC matrix by the SiC coating. The microhardness and toughness of the monolithic SiC ceramics were 25.5 GPa and 4.8 MPa m1/2, respectively. The SiC-coated MWCNTs/SiC composites showed an interesting behavior with the elastic recovery when undergoing an indentation test. SiC-coated MWCNTs/SiC composites showed elastic behavior due to the bridging effect of MWCNTs. Li et al. [10] investigated the microstructural evolution and oxidation resistance of MWCNTs by directly heating silicon powder and MWCNTs in a coke bed from 1000 to 1500°C, with the aid of X-ray diffraction (XRD), scanning electron microscopy (SEM), high-resolution transmission electron microscopy (HRTEM), and thermogravimetry-differential scanning calorimetry (TG-DSC). The results showed that the morphology and microstructure of MWCNTs did not change much after being treated from 1000 to 1200°C. An obvious SiC coating was formed on the surface of MWCNTs from 1300 to 1400°C. Up to 1500°C, nearly all the MWCNTs transformed into SiC nanowires. The oxidation resistance of the treated MWCNTs was improved compared with as-received ones. Nonisothermal kinetics showed that the oxidation activation energy of the treated MWCNTs reached 208 kJ/mol, much higher than 164 kJ/mol of as-received ones. CNT-dispersed Si3N4 ceramics with electric conductivity were developed by Tatami et al. [11] based on the lower temperature densification technique in which the key point was the

Molecular Dynamics Simulation of Ceramic Matrix Composites 231 addition of both TiO2 and AlN and Y2O3 and Al2O3 as sintering aids. This new ceramic with a small amount of CNTs exhibited very high electric conductivity in addition to high strength and toughness. Since Si3N4 ceramics with Y2O3-Al2O3-TiO2-AlN were originally used as a wear material, electrically conductive Si3N4 ceramics are expected to be applied for high-performance static electricity-free bearings for aerospace and other high-performance components. Wilhelm et al. [12] studied the dependence of the mechanical properties of silicon-infiltrated SiC (SiC-Si) composite on the grain size of the original SiC powder. The effects of submicron SiC powders on the bending strength were analyzed. Using α-SiC powders with a mean particle size of 12.8, 6.4, 4.5, and 3 mm a clear linear enhancement of the bending strength with the decrease of SiC particle size was observed. However, a further decrease of the SiC particle size (from 3 to 0.5 mm) brought no increase of the strength and toughness, respectively. Esfehanian et al. [13] developed Cf/(Mo, Ti) Si2-SiC composite using the melt infiltration technique. Two types of carbon fiber preforms were infiltrated with an alloyed melt of Si, Ti, and MoSi2. To prevent the melt from reacting heavily with the fibers and keep the reinforcing effect of the preforms, the melting point of this alloy should be lower than 1600°C. An induction tube furnace was used to infiltrate the carbon preforms. The carbon preform was mechanically driven into the melt and kept there for 10–15 min. The infiltrated samples were studied using XRD and SEM. Silicon-free composites could be obtained using this technique. SiC was the major phase of the matrix. MoSi2 and solid solutions of Si, Mo, and Ti could be found in the matrix. RBSC composites are fully dense materials fabricated by the infiltration of compacted mixtures of SiC and carbon by molten silicon. Free carbon is usually added in the form of an organic resin that undergoes subsequent pyrolysis. The environmentally unfriendly pyrolysis process and the presence of residual silicon are serious drawbacks of this process. Frage et al. [14] described an alternative approach that minimizes the residual silicon fraction by making use of a multimodal particle size distribution, to increase the green density of the preform’s prior infiltration. The addition of boron carbide provided an alternative source of carbon, thereby eliminating the need for pyrolyzed organic compounds. The residual silicon fraction in the RBSC composites, prepared according to the novel processing route, was significantly reduced. Their mechanical properties, in particular the specific flexural strength, were about 15% higher than the value reported for RBSC composites prepared by the conventional approach. In an exhaustive review by Mei et al. [15] the processing techniques of CNT/SiC ceramic matrix composites (CMCs) and CNT/SiC hybrid structures were summarized. Current CNT/SiC CMC fabrication methods include (i) disperse mix densification, (ii) polymer-derived ceramic (PDC), (iii) CVI of CNT preforms, and (iv) LSI; usually, the disperse mix densification process is called a ceramic route, CVI is a gas-phase route, and PDC and LSI belong to liquid-phase route. During these processes, there are three key challenges we should

232 Chapter 5 face: homogeneous dispersion of CNTs in the SiC matrix, minimization of the mechanical damage of the tubes while processing, and an optimum bonding between CNTs and ceramic matrix. Generally, ultrasonication and colloidal processing can generate a homogeneous dispersion of CNTs in SiC matrix; however, ball milling and excessive CNT volume fraction easily lead to agglomeration; on the other hand, ball milling could damage CNTs and hot pressing (HP) and liquid-phase sintering (LPS), producing a degradation of CNTs that might lead to a less efficient reinforcement. Some rapid densification techniques such as spark plasma sintering (SPS) and microwave process have been widely used to fabricate CNT/SiC CMCs by a ceramic route due to their volumetric heating and faster processing. Using PDC method, mixing the preceramic polymers and CNTs is usually the first step; however, in recent years, with the development of CNT assemblies, some CNT fiber/SiC and CNT sheet/SiC composites were successfully fabricated by cyclic precursor infiltration and pyrolysis processes. CNT/SiC CMCs produced by CVI of SiC matrix into porous CNT arrays and sheets have been demonstrated to possess a unique microstructure and excellent micromechanical properties. However, these preforms are usually microscale; therefore fabrication of macroscopic CNT/ SiC composites by CVI process requires that the nanoscale fillers could form macroscopic architectures. In contrast to CVI and PDC, where the SiC matrix buildup does not influence or damage the CNTs, in LSI the molten Si is highly reactive to the C matrix and to the CNTs. Therefore a direct contact of Si melt and CNTs has to be avoided carefully; on the other hand, in the LSI process the raw materials usually contain SiC particles. CNT/SiC hybrid structures contain CNTs on SiC substrates and SiC-coated CNT composites, which were usually fabricated by in situ grown CNTs on the surfaces of various SiC substrates and carbothermic reduction synthesis SiC on CNTs; these CNT/SiC hybrid structures can have wide applications in composites as a micro-/nanoscale hybrid reinforcement and electronic devices. From the earlier studies, it could be inferred that the CNT-reinforced SiC composites have been prepared and tested widely for their mechanical properties. However, to accurately predict the properties of these composites, MD simulation could be of immense help. Through MD simulation the mechanical behavior of these composites can be studied, and the experimental studies can be analyzed. In this study the effect of CNT reinforcement on SiC matrix has been studied at the atomic level by using MD approach. For this purpose the BIOVIA Materials Studio [16] software has been used for modeling and analysis of the composites.

5.1.2 MD Methodology All the MD simulations have been performed using BIOVIA Materials Studio [16]. Fig. 5.1.1 shows the MD model of MWCNT-reinforced SiC composite. The MWCNT has inner tube as armchair (6,6) and outer tube as chiral (16,6). The volume fraction (Vf) of MWCNT in SiC ˚ 3 containing a total of 50,450 matrix was 5.6%. The cell size was 76.95  76.95  100.5 A atoms. A number of similar models with varying Vf were modeled using the “Build” module

Molecular Dynamics Simulation of Ceramic Matrix Composites 233

Fig. 5.1.1 MD model of MWCNT-reinforced SiC composite containing 50,450 atoms and having Vf ¼ 5.6%.

of BIOVIA Materials Studio [16]. Two such models having Vf of 9% and 16% have been shown in Fig. 5.1.2. A number of such models were simulated for predicting the effect of varying Vf on the properties of SiC-based composites. Before applying any strain on these models, they were subjected to thermal stabilization at 300 K and zero pressure using the constant number of atoms (N), constant pressure (P), and constant temperature ensemble (T), that is, the NPT ensemble. After completing the construction of models the simulation procedure consists of geometry optimization, dynamics, and Forcite mechanical property calculation.

5.1.2.1 Geometry Optimization The geometry optimization procedure is necessary because the initial configuration is generally in a higher state of energy. To bring it into a state of equilibrium, we have conducted the geometry optimization with the parameters as listed in Table 5.1.1. A frequent activity in MD simulation is the optimization or minimization (with respect to potential energy) of the system being examined. For instance, it is often desirable to optimize a structure after it has been sketched, since sketching often creates the molecule in a high-energy configuration and starting a simulation from such an unoptimized structure can lead to erroneous results. There are a number of optimization techniques available in BIOVIA Materials Studio, namely, steepest descent method, conjugate gradient, and Newton-Raphson method. In the steepest descent method the line search direction is defined along the direction of the local downhill gradient. Each line search produces a new direction that is perpendicular to the previous gradient; however, the directions oscillate along the way to the minimum. This inefficient behavior is characteristic of steepest descents, especially on energy surfaces having

234 Chapter 5

Fig. 5.1.2 MD models of MWCNT-reinforced SiC composite having (A) Vf ¼ 9% and (B) Vf ¼ 16%. Table 5.1.1 Geometry optimization parameters for MWCNT-SiC composites S. No. 1. 2. 3. 4. 5. 6.

Parameter

Value

Algorithm Quality convergence tolerance Energy convergence tolerance Force convergence tolerance Displacement convergence tolerance Maximum number of iterations

Steepest descent Fine 104 kcal/mol ˚ 0.005 kcal/mol/A ˚ 5  105 A 5000

narrow valleys. Convergence is slow near the minimum because the gradient approaches zero, but the method is extremely robust, even for systems that are far from being harmonic. It is the method most likely to generate the true low-energy structure, regardless of what the function is or where the process begins. Therefore the steepest descent method is often used when the

Molecular Dynamics Simulation of Ceramic Matrix Composites 235 gradients are large and the configurations are far from the minimum. This is commonly the case for initial relaxation of poorly refined crystallographic data or for graphically built models. The reason that the steepest descent method converges slowly near the minimum is that each segment of the path tends to reverse progress made in an earlier iteration. It would be preferable to prevent the next direction vector from undoing earlier progress. This means using an algorithm that produces a complete basis set of mutually conjugate directions such that each successive step continually refines the direction toward the minimum. If these conjugate directions truly span the space of the energy surface, then minimization along each direction in turn must, by definition, end in arrival at a minimum. The conjugate gradient algorithm constructs and follows such a set of directions. As a rule, N2 independent data points are required to solve a harmonic function with N variables numerically. Since a gradient is a vector with N variables, the best we can hope for in a gradient-based minimizer is to converge in N steps. However, if we can exploit second-derivative information, an optimization could converge in one step, because each second derivative is an N  N matrix. This is the principle behind the variable metric optimization algorithms of which NewtonRaphson is perhaps the most commonly used. Another way of looking at Newton-Raphson is that, in addition to using the gradient to identify a search direction, the curvature of the function (the second derivative) is also used to predict where the function passes through a minimum along that direction. Since the complete second-derivative matrix defines the curvature in each gradient direction, the inverse of the second-derivative matrix can be multiplied by the gradient to obtain a vector that translates directly to the nearest minimum. In this chapter, we have used the steepest descent algorithm. The parameters used in the optimization of the nanostructures have been shown in Table 5.1.1.

5.1.2.2 Dynamics Once an energy expression and, if necessary, an optimized structure have been defined for the system of interest, a dynamics simulation can be run. The basis of this simulation is the classical equations of motion that are modified, where appropriate, to deal with the effects of temperature and pressure on the system. The main product of a dynamics run is a trajectory file that records the atomic configuration, atomic velocities, and other information at a sequence of time steps that can be analyzed subsequently. In this study, we have used the NVT ensemble, that is, constant number of atoms, volume, and temperature. The parameters used in simulation have been tabulated in Table 5.1.2. Through this dynamic run the system is allowed to evolve with time to study the interactions between the atoms. For all the simulations, periodic boundary conditions (PBCs) were applied to the simulation box so as to make it as a representative volume element imitating the bulk phase of the nanocomposite. The equations of motion of polymer and nanotube atoms

236 Chapter 5 Table 5.1.2 Dynamics run parameters for MWCNT-SiC composites S. No 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11.

Parameter

Value

Ensemble Initial velocity Temperature Time step Total simulation time Number of steps Frame output Thermostat Collision ratio Energy deviation Repulsive cutoff

NVT Random 298 K 1 fs 100 ps 100,000 50,000 Andersen 1 5  1012 kcal/mol ˚ 6A

were integrated by the velocity-Verlet integrator. The temperature of nanocomposites during simulation was controlled using Andersen’s thermostat.

5.1.2.3 Mechanical Properties The dynamic properties of SiC/MWCNT composites have been computed using “Forcite” module. A constant strain of 0.005 has been applied to the simulation cell. Before calculating the mechanical properties, geometry optimization has been carried out to achieve the minimum energy configuration after the dynamics run. The simulation parameters for Forcite mechanical property calculation are shown in Table 5.1.3. Forcite mechanical property calculations use the “constant strain” approach. The process starts by removing symmetry from the system, followed by an optional reoptimization of the structure, where the cell parameters can be varied. Optimization at this stage is always advised as incorrect results can be obtained if the structure is far from its lowest-energy configuration. For each configuration a number of strains are applied, resulting in a strained structure. The resulting structure is then optimized, keeping the cell parameters fixed. For example, Number of steps for each strain ¼ 4 Max. strain amplitude ¼ 0.003 Strain patterns 100000, 010000 This defines a range of values {0.003, 0.001, 0.001, and 0.003} that are applied to each strain pattern: strain pattern 100000 gives e ¼ {0.003, 0,0,0,0,0}, {0.001, 0,0,0,0,0}, {0.001, 0,0,0,0,0}, and {0.003, 0,0,0,0,0} strain pattern 010000 gives e ¼ {0,0.003,0,0,0,0}, {0, 0.001,0,0,0,0}, {0, 0.001,0,0,0,0}, and {0, 0.003,0,0,0,0}

Molecular Dynamics Simulation of Ceramic Matrix Composites 237 Table 5.1.3 Mechanical property simulation parameters for MWCNT-SiC composites S. No 1. 2. 3. 4. 5. 6. 7.

Parameter

Value

Number of strains Maximum strain Preoptimize structure Algorithm Maximum number of iterations Force field Repulsive cutoff

6 0.005 Yes Steepest descent 10 Compass ˚ 6A

Each strain pattern represents the strain matrix in Voigt notation. It is converted to the strain matrix E, such that E(0,0) ¼ e(0), E(1,1) ¼ e(1), E(2,2) ¼ e(2), E(2,1) ¼ E(1,2) ¼ 0.5*e(3),… These are then used to generate the metric tensor G: G ¼ H00 ½2E + IH0 , 0

where H0 is formed from the lattice vectors, I is the identity matrix, and H0 is the transpose of H0. From G the new lattice parameters can be derived; these are then used to transform the cell parameters (fractional coordinates are held fixed). Following these steps the structure is optimized, and the stress is calculated. A stiffness matrix is built up from a linear fit between the applied strain and resulting stress patterns. In the case of a trajectory, this is averaged over all frames. The PBCs were used in all the simulations. For SiC matrix the embedded atom method (EAM)-based potential [17] has been used for describing the interatomic interactions. The Tersoff potential [18] has been used for describing the atomic interactions within SiC and between SiC and MWCNT.

5.1.3 Results and Discussion Fig. 5.1.3 shows the spread of bulk densities with MWCNT mass fractions from 0 to 20 wt.%. A small decrease in the bulk density was noted after the introduction of MWCNTs, which could be attributed to an increase in the amount of low-density free silicon. As a rule, strongly agglomerated CNTs could lead to porosity and introduce defects in the composites; it is therefore hard to obtain compact CNT-reinforced ceramic composites especially with a high concentration of CNT content [11]. After reactive boding, pores in the green body, which were caused by the addition of MWCNTs, were filled with silicon via capillary pressure, and this could be the reason for the observed decrease in bulk density. As the amount of MWCNTs increased, more amorphous carbon was formed from the decomposition of phenolic resin. This carbon reacted with liquid silicon, and SiC was generated, which led to a slight increase in density. The bulk density of composites with 15 wt.% MWCNTs was reduced to 2.95 g/cm3, which is lower than that of the RBSC. When the mass fraction of MWCNTs reached 15 wt.%, an excessive amount of phenolic resin was introduced. Upon carbonization, that decreased the

238 Chapter 5 3.06 3.04

Bulk density (g/cm3)

3.02 3.00 2.98 2.96 2.94 2.92 2.90 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 22.5 25.0 27.5 30.0 Mass fractions of MWCNTs (%)

Fig. 5.1.3 Variation of bulk density of MWCNT-SiC composite with MWCNT mass fraction.

infiltration of silicon. A similar trend was obtained experimentally for the bulk density by Song et al. [19]. The variation of bending strength and fracture toughness exhibited by the composites (MWCNT content ranges from 0 to 20 wt.%) has been shown in Figs. 5.1.4 and 5.1.5. With an increase in the mass fraction of MWCNTs from 1 to 10 wt.% the bending strength of the 400

Bending strength (MPa)

350

300

250

200

150

100 0.0

2.5

5.0

7.5 10.0 12.5 15.0 17.5 20.0 22.5 25.0 MWCNTs weight (%)

Fig. 5.1.4 Variation of bending strength of MWCNT-SiC composite with MWCNT mass fraction.

Molecular Dynamics Simulation of Ceramic Matrix Composites 239 7.50

Fracture toughness (MPa m1/2)

7.00 6.50 6.00 5.50 5.00 4.50 4.00 3.50 3.00 2.50 2.00 0.0

2.5

5.0

7.5 10.0 12.5 15.0 17.5 20.0 22.5 25.0 MWCNTs weight (%)

Fig. 5.1.5 Variation of fracture toughness of MWCNT-SiC composite with MWCNT mass fraction.

composites gradually increased to 350 MPa. However, when the MWCNT mass fraction increased to 25 wt.%, the bending strength decreased from 350 to 275 MPa, which is only a little higher than that of the matrix. Fig. 5.1.5 shows that the fracture toughness reaches a maximum value of 7 MPa m1/2 (at an MWCNT fraction of 10 wt.%). This indicated that the MWCNTs in the composites preserved their outstanding reinforcement properties. In this work, a larger increase in fracture toughness was obtained compared with that of RBSC combined with the same content of SiC whiskers [4] or chopped fiber [3]. Similar to the behavior observed for bending strength the fracture toughness declines with a mass fraction of 10 wt.%. This decrease in mechanical properties could be associated with agglomeration or damage to the MWCNT and residual carbon.

5.1.4 Conclusion MWCNTs/SiC ceramic matrix composites were successfully modeled by using BIOVIA Materials Studio [16]. This is the first study that has used the MD approach for studying the properties of MWCNT-SiC composites. The effect of increase in MWCNT volume (or mass) fraction on the properties of SiC matrix composites was predicted by varying the mass fraction from 0 to 25 wt.%. As a consequence of their nanosizes and high slenderness ratios the introduction of MWCNTs resulted in an increase in the porosity and the number of defects, which led to a reduction in bulk density of the composites. However, a large number of MWCNTs were preserved in the composites due to the protective SiC coating. For the composite containing 10 wt.% MWCNTs the bending strength and fracture toughness increased

240 Chapter 5 from 240 to 350 MPa and from 3.75 to 7.0 MPa m1/2, respectively. MWCNT pullout and bridging are the major mechanisms behind the enhanced toughness of the composites. The study could be extended further by including the effect of MWCNT aspect ratio (l/d), use of functionalized MWCNTs for improving the bonding strength between the CNTs and SiC matrix, and prediction of damping properties using BIOVIA Materials Studio.

References [1] S.K. Lee, Y.C. Kim, C.H. Kim, Micro-structural development and mechanical properties of pressure-less sintered SiC with plate-like grains using Al2O3–Y2O3 additives, J. Mater. Sci. 29 (1994) 5321–5326. [2] A. Peigney, C. Laurent, E. Flahaut, A. Rousset, Carbon nanotubes in novel ceramic matrix nanocomposites, Ceram. Int. 26 (2000) 677–683. [3] S. Li, Y. Zhang, J. Han, Y. Zhou, Fabrication and characterization of random chopped fiber reinforced reaction bonded silicon carbide composite, Ceram. Int. 38 (2012) 1261–1266. [4] S. Li, Y. Zhang, J. Han, Y. Zhou, Fabrication and characterization of SiC whisker reinforced reaction bonded SiC composite, Ceram. Int. 39 (2013) 449–455. [5] Y. Zhang, Z. Yuan, Y. Zhou, Gelcasting of silicon carbide ceramics using phenolic resin and furfuryl alcohol as the gel former, Ceram. Int. 40 (2014) 7873–7878. [6] H. Mei, Q. Bai, K.G. Dassios, et al., Oxidation resistance of aligned carbon nanotube-reinforced silicon carbide composites, Ceram. Int. 41 (2015) 12495–12498. [7] H. Mei, J. Xia, D. Zhang, D. Han, S. Xiao, L. Cheng, Highly conductive and high-strength carbon nanotubes film reinforced silicon carbide composites. Ceram. Int. (2017), https://doi.org/10.1016/j.ceramint.2017.04.022. [8] S. Suyama, T. Kameda, Y. Itoh, Development of high-strength reaction-sintered silicon carbide, Diamond Relat. Mater. 12 (2003) 1201–1204. [9] Y. Miyamoto, Y. Morisada, Y. Takaura, K. Hirota, N. Tamari, Mechanical properties of SiC composites incorporating SiC-coated multi-walled carbon nanotubes, Int. J. Refract. Met. Hard Mater. 25 (2007) 322–327. [10] Y. Li, M. Luo, S. Jin, S. Sang, L. Zhao, Microstructural evolution and oxidation resistance of multi-walled carbon nanotubes in the presence of silicon powder at high temperatures, J. Mater. Sci. Technol. 28 (2012) 599–605. [11] J. Tatami, T. Katashima, K. Komeya, T. Meguro, T. Wakihara, Electrically conductive CNT-dispersed silicon nitride ceramics, J. Am. Ceram. Soc. 88 (2005) 2889–2893. [12] M. Wilhelm, M. Kornfeld, W. Wruss, Development of SiC-Si composites with fine grained sic microstructures, J. Eur. Ceram. Soc. 19 (1999) 2155–2163. [13] M. Esfehanian, J. Gunster, F. Moztarzadeh, J.G. Heinrich, Development of a high temperature Cf/XSi2–SiC (X ¼ Mo, Ti) composite via reactive melt infiltration, J. Eur. Ceram. Soc. 27 (2007) 1229–1235. [14] N. Frage, S. Aroati, M. Cafri, H. Dilman, M.P. Dariel, Preparation of reaction bonded silicon carbide (RBSC) using boron carbide as an alternative source of carbon, J. Eur. Ceram. Soc. 31 (2011) 841–845. [15] H. Mei, D. Han, S. Xiao, K.G. Dassios, L. Cheng, A review on the processing technologies of carbon nanotube/ silicon carbide composites. J. Eur. Ceram. Soc. (2012), https://doi.org/10.1016/j.jeurceramsoc.2018.04.033. [16] Dassault Syste`mes BIOVIA, Materials Studio, 7.0, Dassault Syste`mes, San Diego, 2017. [17] Y. Mishin, D. Farkas, M.J.M.A. Papaconstantopoulos, Interatomic potentials for monoatomic metals from experimental data and ab initio calculations, Phys. Rev. B 59 (1999) 3393–3407. [18] J. Tersoff, Modeling solid-state chemistry: Interatomic potentials for multicomponent systems, Phys. Rev. B: Condens. Matter 39 (1989) 5566–5568. [19] N. Song, H. Liu, J. Fang, Fabrication and mechanical properties of multi-walled carbon nanotube reinforced reaction bonded silicon carbide composites, Ceram. Int. 42 (2016) 351–356.

Molecular Dynamics Simulation of Ceramic Matrix Composites 241 Chapter 5.2

Molecular Dynamics Simulation of Al/Al2O3 Metal-Ceramic Composite Using LAMMPS Amit Bansal*, Prince Setia†, Sumit Sharma‡ *

Department of Mechanical Engineering, I.K. Gujral Punjab Technical University, Kapurthala, India, Department of Materials Science and Engineering, Indian Institute of Technology, Kanpur, India, ‡ Department of Mechanical Engineering, Dr. B. R. Ambedkar National Institute of Technology, Jalandhar, India †

Metal matrix composites are widely used in aerospace industry and in high-temperature applications such as boilers. High hardness and strength to weight ratio, good corrosion, and abrasion resistance have made these materials highly suitable for space, military, and automotive applications. The interfacial properties of these composites play an important role in determining the strength of these composites. A microstructural analysis of the interfacial properties can help in the determination of mechanical properties of these composites. Atomistic approaches such as density functional theory (DFT) analysis and MD have proved to be very useful in the prediction of properties of nanomaterials. Sharma et al. [1] used the MD method for studying the effect of CNT and graphene reinforcements in a metal matrix. BIOVIA Materials Studio 7.0 was used as a tool for finding the modulus and thermal conductivity of nanocomposites. Two different computational models, single-layer graphene/nickel (SLGS/ Ni) and carbon nanotube/nickel (CNT/Ni) composites, were examined to study the effect of nanofiller geometry on Young’s modulus and thermal conductivity of these nanocomposites. Effect of increase in temperature on Young’s modulus was also predicted using MD. The results showed that for the same volume fraction (Vf), Young’s modulus of the SLGS/Ni composite was higher than that of the CNT/Ni composite, showing that SLGS is a better reinforcement than CNT for Ni-based composites. Ward et al. [2] used MD to examine the deformation and failure of Al-Si interfaces and of nanocomposites consisting of Al and Si nanograins loaded in tension. The modified EAM was used to describe the forces and energies in the system and accounted for the directional dependence of Si bonding. Various Al/Si interfaces with a common [011] axis parallel to the interface and polycrystals of eight hexagonal grains of 5-nm diameter with [011]-oriented columnar grains with zero, one, or two Si particles were examined to understand the effect of Si additions on the deformation and failure of the material. Three main features were demonstrated. First, imperfectly matched Al/Si interfaces had relatively high tensile strengths and fracture energies that were larger than the ideal Griffith fracture energies. Second, different

242 Chapter 5 mechanisms dominated the early stages of nonlinear deformation in the Al polycrystal and Al-Si nanocomposites: the Al polycrystal showed a mix of grain boundary deformation and dislocation activity, while the Al-Si polycrystals showed significant grain boundary sliding/ shearing only at the Al-Si interfaces. The grain boundary sliding/shearing in the Al-Si nanocomposite was also reduced with increasing Si content, leading to higher yield stress. Third the failure modes differ between the polycrystal and the nanocomposite: failure in the Al polycrystal initiated at a grain boundary and propagated transgranularly, whereas failure in the Al-Si composites occurred by void damage accumulation at an Al/Si interface. The composite failure stresses were lower than the strengths of the Al-Si bimaterial interfaces due to stress concentrations in the composites and nonideality of the as-fabricated composite interfaces, but the normal stresses at the particle interfaces did reach values comparable with those of the bimaterial interfaces. These results suggested that Al-Si nanocomposites could be engineered for high hardness strength. Zhou et al. [3] used MD simulation to examine mixed mode loading leading to interfacial delamination between two brittle materials. The propagation of a model crack along the interface was simulated under various combinations of far field tensile (mode I) and shear (mode II) stresses. Local stress and local opening data under the steady-state crack propagation condition were extracted from these MD simulations and then used to determine analytic stress versus opening relations for mixed-mode-loading conditions. Their simulation results enabled the development of a qualitative picture of the traction-separation behavior and functional forms and parameters for a cohesive surface constitutive model consisting of separate normal and shear traction-separation relations. Using a reactive force field (ReaxFF), Zhang et al. [4] investigated the structural, energetic, and adhesion properties of both solid and liquid Al/α-Al2O3 interfaces. The ReaxFF was developed solely with ab initio calculations on various phases of Al and Al2O3 and Al-O-H clusters. The computed lattice constants, elastic constants, surface energies, and calculated work of separation for the solid-solid interface agreed well with earlier first-principle calculations and experiments. For the liquid-solid system, Zhang et al. [4] investigated the nonwetting-wetting transition of liquid Al on α-Al2O3 (0001). The results revealed that the evaporation of Al atoms and diffusion of O atoms in α-Al2O3 led to the wetting of liquid Al on the oxide surface. The driving force for this process was a decrease in interfacial energy. The nonwetting-wetting transition was found to lie in the 1000–1100 K range, which is in good agreement with sessile drop experiments. Based on a kind of pair potential function describing the interaction across the sapphire/copper interface and EAM potential for copper, the tensile mechanical properties of two kinds of Cu (111)/α-Al2O3(0001) interfaces with two types of initial condition, that is, single-layer and double-layer Al-terminated interface with copper atoms nearly atop Al3+ and two successive (110) planes of copper film astride Al3+ and O2 along the [112 direction, were studied by

Molecular Dynamics Simulation of Ceramic Matrix Composites 243 Shen and Yu [5] employing the quasi-continuum method. The curves of loading versus tensile displacement were plotted, and the atomistic structures at different loading stages for the two kinds of interface with the two different initial conditions were also discussed. Although it seemed that in the stable area of relaxed structures for single-layer/double-layer Al-terminated interface, copper atoms locate atop Al3+/O2, most of copper atoms rested atop O2 for the two kinds of Al-terminated interface. Meanwhile the location of interfacial imperfect bonding area had a tremendous influence on the deformation process and mode of interface failure. The strength of single-layer Al-terminated interface was relatively larger than that of double-layer counterpart. The initial positions of coppers atoms had small effects on the tensile strength. The flow stress behavior of the nanometric Al2O3 particulate-reinforced Al alloy matrix composites was investigated by Yang et al. [6] in the temperature range from 360 to 480°C, with strain rates of 0.001–1 s1 on Gleeble 1500 machine. The results showed that the flow stress increased with increasing strain rate and decreasing temperature. The flow stress of the deformed composites showed a dynamic recovery character in the temperature range from 360 to 440°C and in the strain rate range of 0.001–0.01 s1, while that deformed at 480°C showed dynamic recrystallization character. The processing map at the strain of 0.4 was obtained and exhibited three safe deformation domains: 360–400°C at 0.001–0.01 s1, 400–440°C at 0.001– 1 s1, and 440–480°C at 0.001–0.01 s1. According to the processing map and microstructure observation the optimum hot-working parameters are determined to be 480°C/0.001 s1. Siegel et al. [7] performed a series of ab initio calculations to determine the atomic structure, ideal work of adhesion (Wad), and bonding character of the Al (111)/α-Al2O3(0001) interface. Six candidate interface geometries were considered, including Al and O terminations of the oxide. Minimization of the Hellmann-Feynman forces resulted in substantial changes to the atomic structure of the metal near the interface, wherein some atoms adopted positions consistent with a continuation of the oxide’s Al-sublattice crystal structure across the interface. Consequently the lowest-energy structures (i.e., having the largest Wad) were those that facilitated this “oxide extension” mechanism. By applying several methods of analysis a thorough characterization of the electronic structure was performed. It was determined that AldO bonds constitute the primary interfacial bonding interaction. These bonds were very similar to the cation-anion bonds found in the oxide bulk and were mainly ionic yet maintained a small amount of covalent character. In addition, there was evidence of metal-cation bonding at the optimal Al-terminated interface. Taking into account recent theoretical and experimental evidence suggesting an Al termination of the clean oxide surface the calculations predicted Wad ¼ 1.36 J/m2 (local density approximation [LDA]) and 1.06 J/m2 (generalized gradient approximation [GGA]) for the optimal Al-terminated structure, which were in good agreement with the experimental value of 1.13 J/m2 as scaled to 0 K. These values were approximately an order of magnitude smaller than what was found for the optimal O-terminated interface: 10.70 (LDA) and 9.73 J/m2 (GGA). Although cleavage preferentially occurred at the interface for the Al termination, strong bonding at the O-terminated interface favored cleavage within the metal.

244 Chapter 5 Rambaut et al. [8] performed MD simulation of the α-Al2O3 crystal. Structural, thermodynamic, and vibrational properties were determined by using a modified shell model. In this model the energy conservation was ensured very accurately by assigning a small mass to each shell and treating them as dynamic variables. Simulation results were compared with experimental data, particularly with inelastic neutron scattering experiments giving the density of phonon states. Calculated structural and thermodynamic properties were in very good agreement with data obtained from the literature. Inelastic neutron scattering experiments were performed; calculated phonon characteristics compared very well with these experimental results. Liu et al. [9] studied the coherent and semicoherent Al/α-Al2O3 interfaces using MD simulations with a mixed, metallic-ionic atomistic model. For the coherent interfaces, both Al-terminated and O-terminated nonstoichiometric interfaces have been studied, and their relative stability has been established. To understand the misfit accommodation at the semicoherent interface a one-dimensional (1D) misfit dislocation model and a two-dimensional (2D) dislocation network model have been studied. For the latter case the analysis revealed an interface dislocation structure with a network of three sets of parallel dislocations, each with pure-edge character, giving rise to a pattern of coherent and stacking fault-like regions at the interface. Structural relaxation at elevated temperatures led to a further change of the dislocation pattern, which could be understood in terms of a competition between the stackingfault energy and the dislocation interaction energy at the interface. These results could serve as an input for the subsequent dislocation dynamics models to understand and predict the macroscopic mechanical behavior of Al/α-Al2O3 composite heterostructures. The atomic structure and charge distribution pattern of the Al/α-Al2O3 interface were studied utilizing MD simulations by Liu et al. [10]. To accurately describe the interactions between the atoms around the interface the charge transfer ionic and EAM potential was used. Energetically preferable Al/α-Al2O3 interface systems were first determined to study the layer structures of systems. Two
the interface 
 energetically preferable Al/α-Al2O3 interface systems with a 110 ð1 1 1ÞAl ll 1 0 1 0 ð0 0 0 1ÞAl2 O3 orientation relationship were obtained, corresponding to the atop-O Al-terminated and atop-O O-terminated relaxed models, respectively. Further studies revealed the presence of an interfacial layer, which was consistent with experimental ˚ . It results. The models predicted a thickness of the interfacial layer between 12.14 and 16.82 A was composed of aluminum suboxide (with an Al to O atomic ratio between 1:1.07 and 1:1.17). In addition, both the combination between the interfacial layer and the metallic Al layer and the interfacial layer and the ceramic α-Al2O3 were perfect. To further study the atomic structure of the Al/α-Al2O3 interfacial layer, it was isolated from the system and modeled separately. An analysis of the radial distribution function revealed that the interfacial layer inherited its structure from the α-Al2O3 moiety. The study of the charge distributions in the interface systems indicated that the charge of the Al atoms in the interfacial layer was mainly in the range from +2.1e to +2.6e, while the charge of the O atoms is at the saturated state of 2e.

Molecular Dynamics Simulation of Ceramic Matrix Composites 245 In this chapter, MD simulation of Al/Al2O3 metal-ceramic interface has been performed using the ReaxFF potential [11]. From the literature, it was found that the accurate prediction of interfacial properties depend strongly on the type of force field chosen. Because of computational limitations of the empirical potentials, bond formation and bond breakage could not be predicted accurately. Hence ReaxFF potential was used here. All the simulations were performed using the LAMMPS [12].

5.2.1 MD Simulation 5.2.1.1 Interatomic Potential ReaxFF potential [11] was used for MD simulation of the Al/Al2O3 interface. According to this function the potential energy of the system is given by Eq. (5.2.1): ESystem ¼ Ebond + ECoulomb + EvdWaals + Eover

(5.2.1)

Here, Ebond ¼ bond energy derived from bond orders. ECoulomb ¼ Coulomb interactions EvdWaals ¼ van der Waalsenergy Eover ¼ overcoordinationenergy This potential function provides a ReaxFF based on bond-order formalism in conjunction with a charge equilibration scheme. In other words, within this potential the formation and dissociation of the atomic bonding is obtained during each time step. Moreover, formability of the bonds during the simulation with charge variation of atoms arising from the bond length is the major characteristic of ReaxFF potential that distinguishes it from traditional potential functions. This characteristic is particularly helpful in the modeling of atomic bonds in the Al/ Al2O3 interface, where both metallic and ionic atomic bonds are present. It is also sufficiently general to be able to describe both metallic and ionic bonds between Al and O elements [4]. Furthermore, it is well suited to be used in commercial software and open-source MD programs such as LAMMPS.

5.2.1.2 Al and Al2O3 Models The Al FCC crystal was modeled using BIOVIA Materials Studio [13]. Using the “build ˚ 3 was built. The αsupercell” command the Al crystal having dimensions of 48  97  95 A Al2O3 was also modeled with the similar dimensions. The interface was obtained by putting the FCC-structured Al and HCP-structured α-Al2O3 atoms beside each other. The two structures were separated with an initial gap near the equilibration atomic distance in the interface. Fig. 5.2.1 shows the Al/α-Al2O3 interface. This size was chosen because it has been proved [14]

246 Chapter 5 190

97

Al

Al2O3

Fig. 5.2.1 MD model of Al/α-Al2O3 interface.

that the cell of this size resulted in high accuracy and low computational cost. The periodic boundary conditions were applied along both x and y directions for the simulation box. The material was put free along the z-axis in both directions. Utilizing different thermostat schemes the simulations were carried out under the canonical ensemble (NVT). In this approach the initial temperature was kept at 298 K during a specific equilibration time. The simulation time step was adopted equal to 0.4 fs, which was small enough for the ReaxFF potential function. The atomic charge was updated at every time step. To explore the effect of equilibration time on the final results, the stress-strain curves were obtained at simulation times of 25, 50, 75, 100, 125, 150, and 175 ps. The results showed no significant difference between the behaviors after the equilibration times of 125 ps. Therefore the equilibration time of 125 ps seemed to be enough to get an acceptable result. Within the 125 ps equilibration time the Langevin thermostat was utilized for the first 25 ps of the simulation to let a reasonable equipartitioning of the kinetic energy at the beginning. Thereafter the process was followed by a Nose-Hoover thermostat for the last 100 ps. This time interval was long enough to dissipate all the fluctuations in the system. This procedure allowed the system to reach a more physical state. After 125 ps of equilibration the loading procedure was initialized. To exert the shear and ˚ -thick regions at the two ends of the simulation box were tensile loading into the atoms, two 7-A considered as the loading regions. The velocity was imposed on the atoms in these regions at a rate of 109 (s1). To avoid the atoms to be shocked by an abrupt change in the kinematic quantities a uniform stretching was imposed to the other atoms as the initial condition. To get the simulation results, atomic quantities were ensemble averaged in both time and space at every 100 steps. The evolution of the atomic structure was visualized by using BIOVIA Materials Studio [13].

Molecular Dynamics Simulation of Ceramic Matrix Composites 247

5.2.2 Results and Discussion Young’s and shear moduli of Al/α-Al2O3 were obtained using the stress-strain curves obtained using LAMMPS. Young’s modulus was obtained as 170 GPa, while the shear modulus was found to be 68 GPa. The results were found to be in good agreement with experimental results [4] and with the Hashin-Shtrikman bounds calculated from earlier works using the ReaxFF method [4]. These results validated the ReaxFF potential and the parameters used. The results have been tabulated in Table 5.2.1. After the formation of the interface along the equilibration stage the structure was stretched to obtain the mechanical properties of the interface. The stress-strain curve was utilized as the main tool to identify the elastic modulus, shear modulus, toughness, and yield strength of the structure. The stress-strain curve was obtained using the same procedure as has been explained in
 Chapter 5.1. One such stress-strain curve obtained for Al/α-Al2O3 system with ð111Þ 011 orientation of the Al plane in the aluminum has been shown in Fig. 5.2.2. The process initiated with an elastic deformation in the initial structure. After about 20 ps, some vacancies nucleated in the interface that appeared as a deviation from linear response of the curve. Thereafter the increase of the vacancies led to an interfacial separation. At this step an abrupt fall was detected in the stress-strain curve. This process continued until the point that the system lost the load-carrying capability. Then the stress value was nearly equal to zero, and the separation occurred. The atomic forces led to some absorptions and repulsions in this stage that appeared as a spring back effect.

5.2.3 Conclusion The mechanical behavior of Al/α-Al2O3 was analyzed using uniaxial tensile load by making use of MD simulations through LAMMPS. ReaxFF potential was used for accurate prediction of interfacial properties. Young’s modulus was obtained as 170 GPa, while the shear modulus was found
to be 68 GPa. The stress-strain curve was obtained for Al/α-Al2O3 system with ð111Þ 011 orientation of the Al plane in the aluminum. From the stress-strain diagram, different stages were observed such as that of elastic deformation, vacancy nucleation, interfacial separation, the loss of the load-carrying capability, and fracture. The MD simulations could be further extended by using different orientation of the Al plane in Al, and the effect of different types of termination of Al2O3 could also be studied. Table 5.2.1 Comparison of mechanical properties of Al/α-Al2O3 Hashin-Shtrikman Bounds Experimental [4] This study

E (GPa)

G (GPa)

153–200 170

59–78 68

248 Chapter 5 5 4.5 4

Stress (GPa)

3.5 3 2.5 2 1.5 1

Elastic spring-back

0.5 0

0

20,000

40,000

60,000

Time (fs)

Fig. 5.2.2
 Stress-strain diagram for Al/α-Al2O3 system with ð111Þ 011 orientation of the Al plane in the aluminum.

References [1] S. Sharma, P. Kumar, N. Kumar, R. Chandra, Graphene/carbon nanotube reinforced nickel composites: a molecular dynamics study. Compos.: Mech. Comput. Appl.: Int. J. (2018), https://doi.org/10.1615/ CompMechComputAppl IntJ.2018 021097. [2] D.K. Ward, W.A. Curtin, Y. Qi, Aluminum–silicon interfaces and nanocomposites: A molecular dynamics study, Compos. Sci. Technol. 66 (2006) 1151–1161. [3] X.W. Zhou, J.A. Zimmerman, E.D. Reedy Jr., N.R. Moody, Molecular dynamics simulation based cohesive surface representation of mixed mode fracture, Mech. Mater. 40 (2008) 832–845. [4] Q. Zhang, T. Cagin, A. van Duin, et al., Adhesion and nonwetting-wetting transition in the Al/α-Al2O3 interface, Phys. Rev. B 69 (045423) (2004) 1–11. [5] S. Shen, W. Yu, Multiscale study on the tensile fracture of Al-terminated Cu (1 1 1)/α-Al2O3(0 0 0 1) interfaces, Comput. Mater. Sci. 48 (2010) 228–240. [6] Y. Yang, Z. Zhang, X. Zhang, Processing map of Al2O3 particulate reinforced Al alloy matrix composites, Mater. Sci. Eng. A 558 (2012) 112–118. [7] D.J. Siegel, L.G. Hector Jr., J.B. Adams, Adhesion, atomic structure, and bonding at the Al (1 1 1)/α-Al2O3 (0 0 0 1) interface: a first principles study, Phys. Rev. B 65 (085415) (2002) 1–19. [8] C. Rambaut, H. Jobic, H. Jaffrezic, J. Kohanoff, S. Fayeulle, Molecular dynamics simulation of the α-Al2O3 lattice: dynamic properties, J. Phys.: Condens. Matter 10 (1998) 4221–4229. [9] X.Y. Liu, G. Pilania, B.J. Thijsse, R.G. Hoagland, I. Lazic, S.M. Valone, Revisiting the Al/Al2O3 interface: coherent interfaces and misfit accommodation. Sci. Rep. 4 (2014), https://doi.org/10.1038/srep04485. [10] H. Mei, Q. Liu, L. Liu, X. Lai, W. She, P. Zhai, Molecular dynamics simulations of the microstructure of the aluminum/alumina interfacial layer, Appl. Surf. Sci. 324 (2015) 538–546.

Molecular Dynamics Simulation of Ceramic Matrix Composites 249 [11] A.C.T. Van Duin, S. Dasgupta, F. Lorant, W.A. Goddard, ReaxFF: a reactive force field for hydrocarbons, J. Phys. Chem. A 105 (2001) 9396–9409. [12] Official Web of LAMMPS, https://lammps.sandia.gov/doc/Manual.html. [13] Dassault Syste`mes BIOVIA, Materials Studio, 7.0, Dassault Syste`mes, San Diego, 2017. [14] A. Sazgar, M.R. Movahhedy, M. Mahnama, S. Sohrabpour, A molecular dynamics study of bond strength and interface conditions in the Al/Al2O3 metal-ceramic composites, Comput. Mater. Sci. 109 (2015) 200–208.

Chapter 5.3

Molecular Dynamics Simulation of Coaxial Boron Nitride/Carbon Nanotubes Using GROMACS Sumit Sharma, Pramod Rakt Patel Department of Mechanical Engineering, Dr. B. R. Ambedkar National Institute of Technology, Jalandhar, India

CNTs were first reported by Iijima [1]. Since then, CNTs have been attracting much attention to explore their exceptional electronic and material properties. Due to their large aspect ratios and small diameters, CNTs have emerged as potentially attractive materials as reinforcing elements in lightweight and high-strength structural composites. As a one-dimensional structure, CNTs can be thought of as one sheet or multiple sheets of graphene rolled into a cylinder. There are single or multiple layers of carbon atoms in the tube thickness direction, called single-walled CNTs (SWCNTs) or MWCNTs, respectively. According to different chiral angles, SWCNTs can be classified into zigzag (θ ¼ 0 degrees), armchair (θ ¼ 30 degrees), and chiral tubule (0 degrees < θ < 30 degrees). CNTs also have excellent thermal conductivity of about 3000 W/m/K [2]. Their main disadvantage is that CNTs oxidize at temperatures approaching above 400°C when exposed to air. This has limited their widespread usage in electronics industry. However, BNNTs could cure this problem of CNTs. BNNTs were discovered after CNTs [3]. With the rapid development in electronic science and technology, electronic systems are becoming smaller day by day and resulting in rapidly increasing power density, which has aggravated the problem of heat dissipation correspondingly. BNNTs are also hollow tube with hexagonal hole on the surface, making them possess similar fundamental properties as that of a CNT. BNNTs show better mechanical properties [4], thermal stability [5], and oxidation resistance [6]. These properties have made these tubes appropriate for hightemperature applications. The thermal conductivity of BNNTs is also slightly lower than that of CNTs, it is much higher than other metals and nonmetallic materials, so BNNTs are also excellent candidates for thermal interface material (TIM) [7].

250 Chapter 5 Recently a novel structure has been proposed containing a coaxial nanocable model consisting of a CNT core covered with BNNT [8]. The bandgap of CNT is sensitive to diameter and chirality, which makes its conductive property similar to metal or semiconductor, while BNNTs are all insulators regardless of diameter and chirality. Considering the conductivity of this novel structure, it is believed that there will be broad application prospects in nanoelectromechanical systems (NEMSs) such as nanoprobe [9] or nanocable [10]. BNNTs have the same nanostructure as CNTs but are found to exhibit significant resistance to oxidation at high temperatures. The study by Chen et al. [11] revealed that BNNTs were stable at 700°C in air and that some thin NTs having diameter less than 20 nm with perfect multiwalled cylindrical structure could survive up to 900°C. Thermogravimetric analysis revealed an onset temperature for oxidation of BNNTs of 800°C compared with only 400°C for CNTs under the same conditions. This more pronounced resistance of BNNTs to oxidation was inherited from the hexagonal BN and also depends on the nanocrystalline structure. This high level of resistance to oxidation allows promising BNNT applications at high temperatures. Despite the significant amount of research on CNTs the thermal conductivity of individual SWCNTs has not been well established. To date, only a few groups have reported experimental data for these molecules. Existing MD simulation results range from several hundred to 6600 W/m/K, and existing theoretical predictions range from several dozens to 9500 W/m/K. To clarify the several order-of-magnitude discrepancy in the literature, Lukes and Zhong [2] utilized MD simulation to systematically examine the thermal conductivity of several individual (10,10) SWCNTs as a function of length, temperature, boundary conditions, and MD simulation methodology. Nanotube lengths ranging from 5 to 40 nm were investigated. The results indicated that thermal conductivity increases with nanotube length, varying from about 10 to 375 W/m/K depending on the various simulation conditions. Phonon decay times on the order of hundreds of femtosecond were computed. These times increase linearly with length, indicating ballistic transport in the NTs. A simple estimate of speed of sound, which did not require involved calculation of dispersion relations, was presented based on the heat current autocorrelation decay. Agreement with the majority of theoretical/computational literature thermal conductivity data was achieved for the nanotube lengths considered. Discrepancies in thermal conductivity magnitude with experimental data were primarily attributed to length effects, although simulation methodology, stress, and intermolecular potential might also play a role. Quantum correction of the calculated results revealed thermal conductivity temperature dependence in qualitative agreement with experimental data. In the absence of specific partial atomic charge information for boron and nitrogen, Santosh et al. [12] studied the elastic properties of BNNTs by varying the partial atomic charges on boron and nitrogen. The authors computed Young’s (Y) and shear modulus (G) of BNNT as a function of the tube radius and number of walls using molecular mechanics calculation and predicted that the Young modulus of BNNTs increases with increase in magnitude of the partial atomic charge on B and N and could be larger than Young’s modulus of CNTs of same radius.

Molecular Dynamics Simulation of Ceramic Matrix Composites 251 This was in contrast to the earlier finding that CNTs had the largest tensile strength. Shear modulus, on the other hand, depended weakly on the magnitude of partial atomic charge and was less than the shear modulus of the CNT. The values obtained for Young modulus and shear modulus were in excellent agreement with the available experimental results. Chandra et al. [13] studied the thermal buckling behavior of precompressed BNNTs using MD simulations with Tersoff interatomic potential. The authors computed the critical buckling strains at near-zero temperature and subsequently precompressed the NTs at a certain fraction of this value followed by temperature ramping. The critical buckling temperature (Tcr) was marked by a sudden decrease of the internal force. It was observed that (i) at small to moderate lengths, Tcr was higher for chiral NTs than for either armchair or zigzag NTs; (ii) Tcr decreased with increasing diameter unlike in thermal disintegration where disintegration temperatures rose with increasing diameter; and (iii) armchair NTs had an optimal length for which Tcr was maximum. The authors qualitatively explained the reasons for each of the findings. Thermomechanical buckling occurred predominantly in two ways depending on the length of the nanotube, while the shorter NTs failed by radial instability (shell-like behavior); the longer ones invariably failed due to bending-buckling (rodlike behavior). Although many achievements have been made in the thermal properties of CNTs and BNNTs, very less research has been done on the properties of CNT-BNNT composite. In this chapter, MD simulations have been performed using GROMACS package to predict the mechanical and thermal properties of individual CNT-BNNT composite. Stress-strain curves of CNT-BNNT composite with various lengths have been obtained for finding the mechanical properties.

5.3.1 MD Simulation 5.3.1.1 Interatomic Potential Tersoff potential [14] was used to model the interactions of covalent CdC and BdN bonds. As for the nonbonding interaction between CNT and BNNT, van der Waals potential of 12-6 Lennard-Jones (LJ) given by Eq. (5.3.1) was used:    σ 12 σ 6 (5.3.1)  ϕLJ ðRÞ ¼ 4με r r Here, ε ¼ well-depth. σ ¼ equilibrium separation r ¼ distancebetween two pairedatoms μ ¼ adjustmentparameterforinteractionbetween CNT andBNNT ˚, The LJ coefficients were set to be ECdB ¼ 3.294 meV, ECdN ¼ 4.068 meV, σ CdB ¼ 3.411 A ˚ [15,16]. and σ CdN ¼ 3.367 A

252 Chapter 5

5.3.1.2 CNT-BNNT Composite MD simulation of coaxial CNTs embedded in BNNTs was performed using the open-source GROMACS 5.1.1. The effect of variation of length of the NTs on the mechanical properties ˚ . Thermal conductivity of was studied by varying the length of NTs from 21.3 to 85.2 A BNNT-CNT composites was also predicted. The inner CNT tube had a configuration of (9,0), while the outer BNNT had a configuration of (18,0). Both CdC and BdN bond lengths ˚ in these simulations. Fig. 5.3.1 shows the MD model of CNT-BNNT were 1.42 A ˚ . Various other models having lengths of 42.6, 63.9, composite having a length of 21.3 A ˚ were also simulated as shown in Fig. 5.3.2. and 85.2 A In each simulation the time step was taken as 0.5 fs, and the simulation box in each dimension was periodic so that the atoms at two ends shared the same effect of surrounding atoms with others. The constraints produced during the construction of the model were first removed using an optimization algorithm. Then the simulation system temperature was set at 298 K within the NVT ensemble for 250 ps. After the optimization the axial compressive strains were gradually applied to the top end with a fixed strain value of 0.1%, and then a 100 ps equilibration process was followed to establish a new balance. Totally the compression-equilibration circle repeated 100 times until final axial compressive strain of the system was 10%. During the compressionequilibration process the stress of the system was recorded every 50 ps to obtain the stress-strain curves. Young’s modulus was further obtained by fitting the linear part of stress-strain curves.

5.3.2 Results and Discussion Fig. 5.3.3 shows the stress-strain curves of CNT (9,0)-BNNT (18,0) composite with different original lengths. In the beginning the stress increases linearly with the increase of strain. Compression Young’s modulus of CNT (9,0)-BNNT(18,0) obtained by fitting linear region of

L = 21.3Å

Fig. 5.3.1 ˚. MD model of CNT-BNNT composite having a length of 21.3 A

Molecular Dynamics Simulation of Ceramic Matrix Composites 253

(A)

(B)

(C)

L = 85.2Å

L = 63.9Å

L = 42.6Å

Fig. 5.3.2 ˚. MD models of CNT-BNNT composite having length of (A) 85.2, (B) 63.9, and (C) 42.6 A

curves remains unchanged with different nanotube lengths. The compression Young’s modulus of CNT (9,0)-BNNT(18,0) is 0.8 TPa, which was slightly higher than the reported stretching Young’s modulus of this structure as 0.73 TPa [17]. Considering that the compression Young’s modulus of the single BNNT was about 0.72 TPa and less than the compression Young’s modulus of the single CNT 1.1 TPa [7], so the compression Young’s modulus of CNT-BNNT composite structure lying in between that of individual CNT and BNNT was reasonable. To analyze the deformation of the NTs during the compression process a typical stress-strain curve from Fig. 5.3.3 was chosen, and images were taken from different stages of the ˚ -long CNT (9,0)deformation process. Fig. 5.3.4 shows the stress-strain curve of the 21.3-A BNNT (18,0) composite. In the linear region the composite structure of CNT (9,0)-BNNT (18,0) stood upright with increased stress. When the stress dropped a concave fold appeared. When the decrease slowed down, CNT (9,0)-BNNT (18,0) composite still had only one concave fold. When the second large-amplitude stress drop appeared, CNT (9,0)-BNNT (18,0) composite showed two concave folds, then the decline slowed down, and the composite

254 Chapter 5 50 45

L = 21.3

40

Stress (GPa)

35 L = 42.6

30 25 20

L = 63.9

15 10 L = 85.2

5 0 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 Strain (%)

Fig. 5.3.3 Stress-strain diagrams obtained for CNT (9,0)-BNNT (18,0) composite for different lengths of the nanotubes.

50 45

2

40

3

Stress (GPa)

35

4

30

5

1 25 20 15 10 5

0 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 Strain (%)

Fig. 5.3.4 Stress-strain diagram and mechanical behavior obtained for CNT (9,0)-BNNT (18,0) composite ˚. having length of nanotubes as 21.3 A

Molecular Dynamics Simulation of Ceramic Matrix Composites 255 structure CNT (9,0)-BNNT (18,0) composite continued to bend, but two concave folds still appeared in stress descending smooth area. NTs with different original lengths behaved differently. For example, Fig. 5.3.3 showed that smooth region appeared after the linear region ˚ tube instead of a sharply rapid decline. It was due to of the stress increase for the longest 85.2-A the known fact that the longer tubes have enough length to bend without folds. In this chapter, imposed flux method of Jund [18] was used to find thermal conductivity. The number of layers in which the direction of flux was divided, was fixed at 40. Increasing the number of layers could increase the accuracy of the gradients, but too many layers might lead to large fluctuations in the layer temperatures. Two types of exchange method were used in the study. First was “VARIABLE” that exchanged a variable energy between one object in the hot layer and one in the cold layer. “FIXED” exchanged a constant energy between all hot objects in the hot layer and all objects in the cold layer. The amount of energy to exchange in each step when using the “FIXED” exchange type was taken as 1 kcal/mol. The flux was determined by the ratio of exchange energy and number of steps. The number of exchanges during the equilibration stage was taken as 500. During the equilibration stage a thermostat (NVT) acts on the system. The number of exchanges during the production stage was equal to 1000. The production stage was carried out at constant energy (NVE). A time step of 1 fs was used in the simulation. The number of time steps in between two exchanges was fixed at 100. Decreasing the number of steps led to higher fluxes and increased the temperature gradient. Too small values were avoided as these introduced nonlinear effects and might impact performance. The thermal conductivity (k) of the model could be calculated using Eq. (5.3.2) [19]: k¼

J A  rT

(5.3.2)

Here, J ¼ constant heat flux. A ¼ cross-sectional area. rT ¼ temperature gradient. The temperature at each layer was computed using Eq. (5.3.3) [19] as given later: Ti ¼ Here, Ti ¼ temperature of the ith layer. N ¼ the number of atoms in the layer. kB ¼ Boltzmann’s constant. mj ¼ atomic mass of atom j. pj ¼ momentum of atom j.

2 X p2j j 2m 3NkB j

(5.3.3)

256 Chapter 5 450

Thermal conductivity (W/m/K)

400 CNT (9,0)

350 300 250 200 150

BNNT (18,0)-CNT (9,0)

100 50 0 0

100 200 300 400 500 600 700 800 900 1000 1100 1200 1300

Temperature (K)

Fig. 5.3.5 Variation of thermal conductivity with temperature of CNT (9,0)-BNNT (18,0) composite and ˚. SWCNT (9,0) having length of nanotubes as 85.2 A

The average temperatures at each slab were computed from the transient temperatures during the last 200 ps. In this way a linear temperature gradient was established in the middle of the model. By calculating the established rT, thermal conductivity k was obtained using Eq. (5.3.2). Fig. 5.3.5 shows the thermal conductivities of CNT (9,0)-BNNT (18,0) composite, ˚ . In and SWCNT (9,0) varied with temperatures, and both of their original lengths were 85.2 A the temperature range of 100–1300 K, thermal conductivity of the CNT (9,0)-BNNT (18,0) composite decreased from 325 to 75 W/m/K. In addition, from the diagram, we could see that the thermal conductivity decreased with the increasing temperature for both CNT (9,0)-BNNT (18,0) composite and SWCNT (9,0). The main reason was that the heat conduction mode of the nanotube was accomplished by the phonon. As the temperature increased, more high-frequency phonons were excited, but the high-frequency phonon contributed less to the thermal transport, while the increase of the temperature reduced the average free path of the phonons, thus increasing the probability of the phonon collisions and then influencing the heat transport of the phonons. Moreover, as the temperature was increased, effect of the collision between phonons becomes saturated, which slowed down the decline in thermal conductivity. In addition, it could be seen from Fig. 5.3.5 that the thermal conductivities of the composite structure of CNT (9,0)-BNNT (18,0) were lower than those of SWCNT (9,0) at each temperature. This was mainly due to the existence of van der Waals force between the inner and outer tube, which caused the scattering between the phonons and then affected the thermal

Molecular Dynamics Simulation of Ceramic Matrix Composites 257

Thermal conductivity (W/m/K)

250

200

150

100

BNNT (18,0)CNT (9,0) length = 85.2

50

0 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 Strain (%)

Fig. 5.3.6 Variation of thermal conductivity with strain of CNT (9,0)-BNNT (18,0) composite having length of ˚. nanotubes as 85.2 A

transport of the phonons, and the thermal conductivity of the single-walled BNNTs was slightly lower than that of the SWCNTs [2], but the overall relative difference was 18.7%, so the composite structure did not significantly reduced the thermal conductivity of the NTs. Fig. 5.3.6 shows that the thermal conductivity decreased monotonically with the increasing of compression strains until the critical one (5%). When approaching the critical strain the thermal conductivity dropped sharply. After that the thermal conductivity changed little and tends to be a stable value. When the strain was bigger than 5%, CNT (9,0)-BNNT (18,0) composite could be in the depressed double-fold states, which might be the reason for the convergence of axial thermal conductivity.

5.3.3 Conclusion MD simulation of coaxial CNT embedded in BNNTs was performed using the open-source GROMACS 5.1.1. The effect of variation of length of the NTs on the mechanical properties was ˚ . Thermal conductivity of BNNT-CNT studied by varying the length of NTs from 21.3 to 85.2 A composites was also predicted. Young’s modulus was predicted to be 0.8 TPa. The intertube coupling had little effect on Young’s modulus of the CNT-BNNT composite; however, the critical strain and the maximum stress could be increased by increasing the intertube coupling strength. In the temperature range of 100–1300 K, thermal conductivity of the CNT-BNNT ˚ decreased from 325 to 75 W/m/K and was overall about composite with the length of 85.2 A

258 Chapter 5 18.7% lower than that of the individual CNT, which was due to the stronger phonon-phonon scattering. The results show that the thermal conductivity decreased monotonically with the increasing compression strains until the critical one ( 5%). When approaching to the critical strain the thermal conductivity dropped sharply. Finally, it tends to an equilibrium situation in which the thermal conductivity changed very slowly with further increase of the compression strain.

References [1] S. Iijima, Helical microtubules of graphitic carbon, Nature 354 (56) (1991) 56–58. [2] J.R. Lukes, H. Zhong, Thermal conductivity of individual single-wall carbon nanotubes, ASME J.Heat Transfer 129 (2006) 705–716. [3] A. Rubio, J.L. Corkill, M.L. Cohen, Theory of graphitic boron nitride nanotubes, Phys. Rev. B 49 (1994) 5081–5084. [4] N.G. Chopra, A. Zettl, Measurement of the elastic modulus of a multi-wall boron nitride nanotube, Solid State Commun. 105 (1998) 297–300. [5] K.M. Liew, J. Yuan, High-temperature thermal stability and axial compressive properties of a coaxial carbon nanotube inside a boron nitride nanotube, Nanotechnology 22 (2011) 085701 (1–7). [6] S.P.S. Arya, A.D. Amico, Preparation, properties and applications of boron nitride thin films, Thin Solid Films 157 (1988) 267–282. [7] T. Li, Z. Tang, Z. Huang, J. Yu, A comparison between the mechanical and thermal properties of single-walled carbon nanotubes and boron nitride nanotubes, Physica E 85 (2017) 137–142. [8] J.H. Yuan, K.M. Liew, Structural stability of a coaxial carbon nanotube inside a boron-nitride nanotube, Carbon 49 (2011) 677–683. [9] Y.D. Kuang, S.Q. Shi, P.K.L. Chan, C.Y. Chen, Self-assembly of carbon nanotubes and boron nitride nanotubes into coaxial structures, Comput. Mater. Sci. 50 (2010) 645–650. [10] Z.H. Zhang, W.L. Guo, G.A. Tai, Coaxial nanocable: carbon nanotube core sheathed with boron nitride nanotube, Appl. Phys. Lett. 90 (2007) 133103 (1–3). [11] Y. Chen, J. Zou, S.J. Campbell, G.L. Caer, Boron nitride nanotubes: pronounced resistance to oxidation, Appl. Phys. Lett. 84 (2004) 2430–2432. [12] M. Santosh, P.K. Maiti, A.K. Sood, Elastic properties of boron nitride nanotubes and their comparison with carbon nanotubes, J. Nanosci. Nanotechnol. 9 (2009) 1–6. [13] A. Chandra, P.K. Patra, B. Bhattacharya, Thermomechanical buckling of boron nitride nanotubes using molecular dynamics, Mater. Res. Express (2016). http://iopscience.iop.org/2053-1591/3/2/025005. [14] V. Verma, V.K. Jindal, K. Dharamvir, Elastic moduli of a boron nitride nanotube, Nanotechnology 18 (2007) 435711. [15] L.F.C. Pereira, B. Mortazavi, M. Makaremi, T. Rabczuk, Anisotropic thermal conductivity and mechanical properties of phagraphene: a molecular dynamics study, RSC Adv. 6 (2016) 57773–57779. [16] B. Mortazavi, Z.Y. Fan, L.F.C. Pereira, A. Harju, T. Rabczuk, Amorphized graphene: a stiff material with low thermal conductivity, Carbon 103 (2016) 318–326. [17] A.N. Enyashin, A.L. Ivanovskii, Mechanical and electronic properties of a C/BN nanocable under tensile deformation, Nanotechnology 16 (2005) 1304. [18] P. Jund, R. Jullien, Molecular dynamics calculation of the thermal conductivity of vitreous silica, Phys. Rev. B 59 (1999) 13707–13711. [19] B. Mortazavi, T. Rabczuk, Multiscale modeling of heat conduction in graphene laminates, Carbon 85 (2015) 1–7.