CC BY
6
DOI: 10.17516/1997-1397-2021-14-6-787-794 УДК 539.3
Molecular-dynamic Investigation of Silicon Carbide Fracture Under External Mechanical Loads
Andrey V. Utkin* Vasily M. Fomint
Khristianovich Institute of Theoretical and Applied Mechanics SB RAS
Novosibirsk, Russian Federation
Received 12.07.2021, received in revised form 22.08.2021, accepted 20.09.2021 Abstract. In this study, molecular dynamic simulations of quasi-static compression of silicon carbide nanorod, were performed. A longitudinal through defect in the form of a cylindrical channel was introduced into the central part of the nanorod. The influence of the cross sectional size of this internal channel on the strength properties was investigated.
Keywords: molecular dynamics, silicon carbide, quasi-static compression.
Citation: A.V. Utkin, V.M. Fomin, Molecular-dynamic Investigation of Silicon Carbide Fracture Under External Mechanical Loads, J. Sib. Fed. Univ. Math. Phys., 2021, 14(6), 787-794. DOI: 10.17516/1997-1397-2021-14-6-787-794.
Introduction
Silicon carbide is a promising and interesting material, owing to its perfect chemical and physical stability, great mechanical properties, low density and high thermal conductivity. It is widely used in the production of various friction mechanisms, optoelectronic devices, gas turbines and heat exchangers. Silicon carbide is also used in electronic devices designed for operation in aggressive external conditions. Composites based on silicon carbide have high specific mechanical and thermal oxidative properties, which leads to their huge potential for use in various industries, including aircraft engine building.
In recent decades, methods of producing composites based on a silicon carbide matrix reinforced with SiC fibers have been actively developed [1, 2]. Such materials combine high strength and creep resistance, which allows using them as elements of gas turbine engine structures.
Analysis and generalization of available scientific and scientific-technical information on methods of SiC composites formation showed that works are conducted in all directions of research — both development of methods of synthesis of interface coatings, and development of new and modernization of already known methods of ceramic matrix formation, research of mechanical properties of composites, definition of reliability criteria of SiC materials, ways of long-term protection of composites from oxidation, etc. [3, 4].
Despite the fact that many experiments are currently being carried out to obtain silicon carbide materials with the specified properties, there is no scientific and systematic approach to studying the properties of composite silicon carbide material. This is due to the complexity and cost of conducting experiments at the micro level, as well as the lack of a large variety of available source materials of mass production. Thus, it becomes relevant to carry out numerical experiments, in particular, molecular dynamics modeling [5, 6].
* utkin@itam.nsc.ru https://orcid.org/0000-0003-1374-5784 t fomin@itam.nsc.ru https://orcid.org/0000-0002-2811-0143 © Siberian Federal University. All rights reserved
There are a significant number of papers devoted to the study of SiC nanowires under external loads [7-10]. In the paper [8] the behavior of 3C-SiC under the different kind of quasi-static external loading (tension, compression, torsion and their different combination) was studied. It was found that Young's modulus increases with the system size and depends of the initial system temperature. Also critical stress decreases with increasing of system temperature. It was shown that short and long nanowires exhibit different modes for collapse under compression (yielding and column buckling). Critical length of nanowire was found. Below this value MD results perfectly fits by Euler column theory. Under torsion strain nanowire collapses with phase transformation from solid to amorphous structure. In [9] twinned SiC nanowires behavior was studied. Influence of segment thickness and diameter under different kind of external loading was investigated. In paper [10] nanowires of 2H-SiC and 3C-SiC structure under tensile and compression are considered. Simulation results show the difference in structural and mechanical responses of these two polytypes.
1. Methodology and physical system
The main advantage of the molecular dynamics (MD) method when modeling processes at the microlevel is that MD does not require the use of additional models, as is done in the mechanics of a deformed solid [11-14]. Thus, the proposed approach is based on Newton's law (or Hamilton's theory), while the thermodynamic and physical properties of materials are determined by the interatomic interaction potential.
Nanorods of 3C-SiC with a hexagonal cross-section were considered as a physical system in the presented numerical simulation. Nanorod circumradius was approximately 98 A and nanorod length was approximately 202 A. An artificial defect in the form of an internal through longitudinal channel of cylindrical cross-section was introduced into the center of the rod (Fig. 1). The effect of the size of this internal channel on the strength characteristics of the ceramic rod under quasi-static compression was studied. The radius of the internal channel varied in a wide range of values from 5 A to 65 A. The ratio of the length to the transverse size was chosen so that the nanorods deform by yielding. Crystal structure along [111] direction was used as a computational model.
The integration of the equation of motion was carried out using a second order velocity Verlet scheme, and the time step in numerical calculation was 10~16 s. For numerical simulation, a parallel algorithm was developed and software implemented, allowing calculations using NVIDIA general-purpose graphic processors (GPUs) supporting the CUDA (Compute Unified Device Architecture) technology [15-17].
The system was cooled using the artificial viscosity method, as a result of which the ensemble of atoms at the end of the cooling process was at a minimum of the potential energy. In this case, the initial state of the system was close to zero (10~7 K).
In this work, the Vashishta potential [18] was used to model the process of silicon carbide deformation. This potential is based on the Stillinger-Weber formalism, and is well-tested for modelling various shock-wave phenomena in ceramic substances. The potential correctly describes the fracture of various polymorphic forms of silicon carbide, and the waves of structural transformations in the substance, under the influence of intense external loading. Models based on the Vashishta potential show excellent agreement with the available experimental data.
The interatomic interaction potential consists of the sum of two- and three-body parts:
where U is the total potential energy of the system, and rij and rik are the distance between ij and ik atoms.
U = E (rij) + E Ujk (rij,Tik),
(1)
i<j
i,j<k
Fig. 1. 3C-SiC nanorod initial structure
The two-body part (1) includes the members responsible for steric size effect of the ions, Coulomb interactions, dipole charge and van der Waals interactions:
U2 (r) = Hj + exp (—r/A) - Dj exp (-r/i) - Wj, (2)
the potential constants Hj,rjij,Zi,A,Djare given in the article [18].
The three-body part (1) is responsible for structural transformations under pressure and the melting process:
Ujk (rij,rik) = R3 (rij, rik) P3 (9ijk), (3)
where
R3 (rij,rik) = Bijk exp (-Y--1--Y— ) @ (ro - rij) @ (ro - rik), (4)
\rij r0 rik r0 J
cos 0ijk - cos 0ijk)
P3 (Oijk) =--• (5)
1 + Cijk (cos Oijk - cos Oijk)
Here, Bijk defines the bond strength, Oijk is the angle formed by rij and rik vectors, Cijk, Oijk are the constants, @ (ro - r) is the Heaviside step function and ro = 2.9 A. The values of constants included in the three-body part of the interatomic interaction potential are given in the article [18].
In order to save computational resources, the two-body part of the interaction potential (2) was modified. The potential was cut off at a distance of rc = 7.35 A and shifted so that both the potential and its derivative would vanish at a distance of the cut-off radius:
J2shifted (r)=l Ц (Г) - Ц (Гс) - (г - ГС)1 f-kj , Г^ ГС (g)
ij I \ / r=rc
0, Г > Гс
To simulate the compression of the structure in the longitudinal direction, a harmonic interaction potential clamp was applied to the atoms located at the opposite faces of the nanorod [19, 20].
w (r, t) = k ((Xl - X0 + v0t)2 + y - V0)2 + (zi - z0)2) . (7)
Where x0,y0, z0 are the coordinates of the facet atoms at the initial moment but after the cooling procedure, k is the clamping stiffness, and vq is the velocity of the clamp acting on a given atoms.
2
On the left face, the clamp was at a quiescent state and its velocity was equal to zero. On the right face the clamping velocity increased linearly from zero to 5 m/s for 200000 steps of integration, and then remained constant throughout the process. Choosing a linearly increasing velocity rather than a constant one enables a more accurate simulation of the quasi-static compression process.
The following macro parameters were calculated to analyze the dynamics of the process: strain e = v0t/L0, where L0 is the initial length of the crystal, the potential part of the internal energy and the kinetic part of the internal energy of the whole nanorod and the surface atoms, the force and stresses acting on the nanorod from the clamping side.
Fig. 2. Force acting on the system versus strain. A — defect-free nanorod, B — internal channel radius is 20 A, C — internal channel radius is 35 A, B — internal channel radius is 55 A
The distribution of internal energy, temperature and stresses over the rod during compression is also monitored. Since uniaxial compression is simulated in the present task, the calculation
1
S
surface area perpendicular to the force, Fij is the longitudinal component of the force acting on atom i from atom j, assuming that the atoms are on opposite sides of the surface. The visualisation software OVITO was used to visualise the simulation results [21].
of the stresses was carried out using the following expression: axx = — Fj, where S is the
2. Results
In the compression process, the nanorod is initially deformed elastically, followed by irreversible deformation accompanied by structural changes. This can be clearly seen in Fig. 2, which shows the dependence of the force acting on the rod from the clamp side versus strain.
At the strain level of eE = 0.11, the linear growth of the force stops, and then until eP = 0.12 it behaves non-monotonically, and then starts to decrease rapidly for all four represented dependences. In the elastic region of the strain under e < eE, the temperature proportional to the internal kinetic energy is close to zero. At e > eE the temperature starts to rise sharply, indicating a displacement in the crystal lattice. Each act of irreversible deformation leads to a sharp local dissipation of energy and, consequently, to an increase in temperature (Fig. 3). Fig. 4 shows the stress-strain relations for a defect-free nanorod as well as for seven different internal channel diameters. It is easy to see that up to a certain size the presence of an internal channel has practically no effect on the strength properties of the rod (this group of almost coincident dependencies is denoted by the letter A).
Fig. 3. Kinetic temperature of the system versus strain. A — defect-free nanorod, B — internal channel radius is 20 A, C — internal channel radius is 35 A, B — internal channel radius is 55 A
The first differences in the stress-strain diagram begin to appear at a channel radius of 55 A (curve B). A further increase in the diameter of the channel leads to a sharp decrease in the compressive strength (curve C corresponds to radius 65 A). The physical system with a through channel radius of 55 A has a 37% lower weight than the solid rod, while the physical system with a channel radius of 65 A will have a 52% lower weight.
Fig. 4. Strain-stress diagrams. Group A — defect-free nanorod and nanorods with internal channel radius 5 A, 10 A, 15 A, 20 A, 40 A, B — internal channel radius is 55 A, C — internal channel radius is 65 A
It is known that the failure of ideal solids under the influence of external loads begins with the growth of micro-cracks and defects on the surface, which spread into the volume as they increase [22, 23]. An analysis of the coordination number of atoms in the rod at the initial moment of time makes it possible to determine the surface atoms — the atoms with coordination numbers less than four. The calculation of potential energy (1) only for atoms with coordination number four gives internal core energy Ucore. The surface energy will be determined as the difference
between the total energy of the system and the energy of the internal core Usurf = U — Ucore. It is also necessary to separate the surface atoms belonging to the external surface of the rod and the surface atoms in the internal channel.
Calculation of the kinetic temperature performed for the external surface atoms and for the surface atoms in the internal channel showed that the kinetic temperature in the system of external surface atoms increases earlier, indicating the appearance of irreversible structural deformations in the surface layer of atoms. However, as the size of the internal channel increases, the time interval between the beginning of deformation in the external layer and in the layer of atoms in the channel decreases. Analysis of the potential energy of the surface per particle showed that with increasing channel radius the potential energy value of the atoms is tending to the potential energy of the external surface and at radius 65 A they take on comparable values (Figs. 5, 6 and 7).
Fig. 5. Surface potential energy per particle as a function of time. A — external surface, B — internal channel surface, radius of channel is 30 A
Fig. 6. Surface potential energy per particle as a function of time. A — external surface, B — internal channel surface, radius of channel is 55 A
The ratio of surface energies at small dimensions of the internal channel explains why defects appear first on the external surface: the system needs to absorb less energy to produce changes in the structure of the substance on the external surface (Figs. 5 and 6). As the size of the internal
channel increases, the external and internal surfaces become energetically equivalent and defects can be nucleated on both the external and internal surfaces (Fig. 7). This may be the explanation for the drastic reduction in the compressive strength of the rod when the internal channel radius reaches 65 A.
-1120 -|
-1140
-1160 --
-1180 -
-1200
•io-21, J
20
40
60
80 t,Ps
Fig. 7. Surface potential energy per particle as a function of time. B — internal channel surface, radius of channel is 65 A
A — external surface,
0
Conclusion
The process of quasi-static compression of 3C-SiC silicon carbide rods with hexagonal cross-section was investigated with the method of molecular dynamics. The ratio of the length to the transverse size was chosen so that the nanorods deform by yielding.
An artificial defect in the form of an internal through longitudinal channel of cylindrical cross-section was introduced into the center of the rod. The influence of the size of this internal channel on the strength properties of a ceramic nanorod in quasi-static compression has been investigated. It was found that up to a certain size (37% weight reduction compared to a solid nanorod) the internal channel has almost no effect on the strength properties.
The present investigation was supported by Russian Science Foundation grant number 21-1900733. This research was supported in through computational resources provided by the Shared Facility Center "Data Center of FEB RAS" (Khabarovsk).
References
[1] G.S.Corman, K.L.Luthra, Handbook of Ceramic Composites, Springer US, 2005.
[2] C.P. Deck et al., Back Progress in Nuclear Energy, 57(2012), 38-45.
[3] Advances in High Temperature Ceramic Matrix Composites and Materials for Sustainable Development; Ceramic Transactions, Vol. CCLXIII Ed. Mrityunjay Singh et al. Hoboken, NJ, USA, John Wiley & Sons, Inc., 2017.
[4] 10th International Conference on High Temperature Ceramic Matrix Composites - HT-CMC 10: Book of Abstracts. Bordeaux, France, 2019.
[5] D. Frenkel, B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, Amsterdam: Elsevier Academic Press, 2001.
M.P.Allen, D.J.Tildesley, Computer Simulation of Liquids, Oxford University Press, 1987. M.A.Makeev, D.Srivastava, Phys. Rev. B, 74(2006), 165303.
Zhiguo Wang, Xiaotao Zu, Fei Gao, William J. Weber, Phys. Rev. B, 77(2008), 224113.
Zhijie Li, Shenjie Wang, Zhiguo Wang, Xiaotao Zu, Fei Gao, William J. Weber, Journal of Applied Physics, 108(2010), 013504.
H.Tsuzuki, J.P.Rino, P.S.Branicio, J. Phys. D: Appl. Phys., 44(2011), 055405.
A.V.Utkin, V.M.Fomin, AIP Conference Proceedings, 1893(2017), 030018.
A.E.Buzyurkin, E.I.Kraus, Y.L.Lukyanov, AIP Conference Proceedings, 1770(2016), 030091.
E.I.Kraus, I.I.Shabalin, Journal of Physics Conference Series, 653(2015), 012085.
Andrey E. Buzyurkin, Eugeny I. Kraus, Yaroslav L. Lukyanov, Thermal Science, 23(2019), S471-S476.
A.V.Utkin, V.M.Fomin, E.I.Golovneva, AIP Conference Proceedings, 2288(2020), 030083. A.V.Utkin, Mathematica Montisnigri, 39(2017), 101-109.
M.S.Ozhgibesov, T.S.Leu, C.H.Cheng et al., Journal Of Molecular Graphics & Modelling, 38(2012), 375-381.
P. Vashishta, R.K. Kalia, A. Nakano, J.P. Rino, Journal of Applied Physics, 101(2007), 103515.
I.F.Golovnev, E.I.Golovneva, A.V.Utkin, Procedia Structural Integrity, 13(2018), 1632-1637. DOI: 10.1016/j.prostr.2018.12.343
20] I.F.Golovnev, E.I.Golovneva, A.V.Utkin, Physical Mesomechanics, 22(2019), 420-431.
21] A.Stukowski, Modell. Simul. Mater. Sci. Eng., 18(2010), 015012.
22] A.F.Ioffe, The Physics of Crystals, New York: McGraw-Hill Book Company Inc., 1928.
23] Surface Layers and Internal Interfaces in Heterogeneous Materials, Ed. by V.E. Panin, Novosibirsk, Izd-vo SO RAN, 2006 (in Russian).
Молекулярно-динамическое моделирование разрушения карбида кремния под действием внешней механической нагрузки
Андрей В. Уткин Василий М. Фомин
Институт теоретической и прикладной механики им. С. А. Христиановича СО РАН
Новосибирск, Российская Федерация
Аннотация. Представленная работа посвящена исследованию поведения стержня карбида кремния 3C-SiC при квазистатическом сжатии. В центр стержня вносился продольный сквозной дефект в виде канала цилиндрического сечения. Было исследовано влияние поперечного размера этого внутреннего канала на прочностные характеристики стержня.
Ключевые слова: молекулярная динамика, карбид кремния, квазистатическое сжатие.
