Functional role of polycrystal grain boundaries and interfaces in micromechanics of metal ceramic composites under loading

Computational Materials Science xxx (2015) xxx–xxx

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Functional role of polycrystal grain boundaries and interfaces in micromechanics of metal ceramic composites under loading V.E. Panin a,b, V.E. Egorushkin a, D.D. Moiseenko a,⇑, P.V. Maksimov a, S.N. Kulkov a,b, S.V. Panin a,b a b

Institute of Strength Physics and Materials Science SB RAS, Tomsk 634021, Russia National Research Tomsk Polytechnic University, 30 Lenin Ave, Tomsk 634050, Russia

a r t i c l e

i n f o

Article history: Received 8 June 2015 Received in revised form 28 October 2015 Accepted 30 October 2015 Available online xxxx Keywords: Polycrystals Grain boundaries Nonlinear wave flows Computer simulation Additive technologies

a b s t r a c t Surface layer and all interfaces in a solid under loading are considered as an autonomous functional subsystem where the primary plastic shears are initiated and developed. The grain boundaries nonlinear wave flows are depended on a crystal lattice curvature. Computer simulations of grain boundaries rotational wave flows for various lattices curvature were performed using modified excitable cellular automata technique. This method is offered for taking into account the grain boundary flows for the sake of computer simulation of polycrystal’s behavior under deformation and fracture as well as during processes to take place in additive manufacturing technology. Ó 2015 Elsevier B.V. All rights reserved.

1. Introduction Classical continuum and fracture mechanics are based on a linear approaches and treat plastic deformation and fracture as a single-scale processes [1,2]. In the last two decades, multiscale approaches to description of plastic deformation and fracture have been intensively developed [3,4]. In so doing, the main challenge lies in self-consistent description covering the contiguous scales: pico–nano, nano–micro, and micro–macro, and this comes to the scope of mesomechanics [5,6]. Physical mesomechanics uses a unified multiscale approach to describe plastic deformation and fracture of solids as nonlinear hierarchically organized systems [7,8]. The approach rests on the concept of structural scales of deformation and nonequilibrium thermodynamics of structural or structure–phase transformations under highly nonequilibrium conditions. This multiscale approach is followed in the paper to provide self-consistent description of plastic deformation and fracture in a loaded hierarchically organized system. It should be stressed that the nonlinear approach of physical mesomechanics is very topical question in modern materials science. According to [7,8] a deformable solid should be considered to consist of two subsystems: a 3D-translation-invariant crystal and a 2D-planar subsystems comprising surface layers and all interfaces. The 2D-planar subsystem is characterized by a cluster struc⇑ Corresponding author.

ture and absence of translation invariance [9]. Primary plastic shears in a solid under loading are developed within a planar subsystem as nonlinear waves of structural transformations [10]. Thus a polycrystalline grain boundaries and interfaces in any composites are of primary functional importance in order to justify the description of Hall–Petch equation as well as to analyze processes of mass transfer to take place in metal–ceramic composites during processes to develop in additive manufacturing technology. The functional role of polycrystal grain boundaries and interfaces in loaded metal ceramic composites will be matter of study in this paper. 2. Nonlinear translation–rotation waves of structural transformation in a 2D-planar subsystem of heterogeneous solids under loading The fundamental difference of a 3D-crystalline and a 2D-planar subsystems of a solid is in their thermodynamic state which can be expressed in terms of the extended Gibbs thermodynamic potential [11]:

F ¼ U  TS  rik eik ;

ð1Þ

where A = rikeik is a total work of deformation (here is Einstein summation convention for all repeating indices). From (1) it follows that any kind of translation–rotation structural transformations in a 2D-planar subsystem can be described by the term A = rikeik. Common used expression for the Gibbs thermodynamic potential

http://dx.doi.org/10.1016/j.commatsci.2015.10.045 0927-0256/Ó 2015 Elsevier B.V. All rights reserved.

Please cite this article in press as: V.E. Panin et al., Comput. Mater. Sci. (2015), http://dx.doi.org/10.1016/j.commatsci.2015.10.045

2

V.E. Panin et al. / Computational Materials Science xxx (2015) xxx–xxx

takes into account only hydrostatic molar volume V changes of a solid. The condition of continuity preservation in a hierarchically organized system under loading is defined by the expression N X rot Ii ¼ 0;

Let us consider a channeled defect flow in a 2D-planar subsystem along the direction L. The general coordinate system is chosen such that the Z axis is directed along L, whereas the x and y coordinates vary within the deformed layer thickness. According to [12,13], the plastic flow distribution in the local region (r < L) has the form:

ð2Þ

i¼1

where Ii is the flow of strain-induced defects at the i-th structuralscale level. The nonlinear theory of the strain-induced defects flows to develop in a 2D-planar subsystem is proposed and described in [12,13]. In doing so specific wave equations were derived for nondimensional linear defect flux J and density a (discontinuities of the displacement vector u). These relations include continuity equation for a defect medium from which it follows that the plastic flow source can be described as the defect rearrangement rate. Of critical importance are temporal variations in the medium density. In this case they are determined not by the divergence but by the curl rot of the defect flux, i.e., by its spatial inhomogeneity. The defect continuity condition is used for taking into account the absence of charges of the vortex component of the plastic strain field. A well-known quasielastic equilibrium equation that used in continuum mechanics incorporates plastic distortions in addition to elastic strain. In effect, this relation takes into account nucleation of straininduced defects in local hydrostatic tension regions generated by a stress concentrator. Moreover, the plastic distortion enters into these equations as defect sources. This is an evidence of a dual nature of defects as field sources. If this takes place the plastic distortion plays a crucial role in the wave propagation of plastic strain. It should be emphasized that although wave equations of plastic flow in solids were derived earlier in [14,15], they were not interpreted as plasticity waves. A conclusion on wave propagation of elastic–plastic perturbation in a medium is always reasoned from the problem of perturbation group velocity. With no group velocity dispersion, the wave is well defined. The medium inhomogeneity gives rise to wave packet dispersion and breakup. Hence, in the context of the single-scale approach the existence of plasticity waves is beyond question. However, if we consider a deformed solid as a multiscale system, taking into account the ‘‘shear plus rotation” pattern of plastic deformation and possibility of channeled local structural transformations, the conclusion on nonlinear plasticity and failure waves is justified. Moreover, stress concentrators in plastic shear propagation as local structural transformations are impossible to reproduce without using the nonlinear wave pattern. It is no accident that all known wave equations of plastic flow are similar to equations of electrodynamics and provide a qualitatively similar field behavior [13]. From the foregoing discussion it follows that the source of the strain-induced defect density is the shear strain vorticity induced by local structural transformations in the quasi-elastic stress concentrator region. The local structural transformations provide relaxation of both oppositely directed shear stresses and hydrostatic tensile stresses in the vicinity of stress concentrators. On completion of local relaxation, the new stress concentrator arises within the applied stress field. The wave process can be expressed analytically for channeled localized deformation under specified boundary conditions. If we put that the plastic distortion is equal to zero, the theory [12,13] describes a nonlinear wave of crack propagation. In this case, inelastic shear displacements in a hydrostatic tension zone at the crack tip increase the molar volume, and with the proviso that V > Vcr, structural phase decomposition of material ahead of the crack tip takes place, thus providing propagation of a nonlinear fracture wave.

J¼

  b1  b2 2L vðs; tÞbðs; tn Þ ln  1  rf ; r 4p

ð3Þ

where b is the binormal vector in the local coordinate system, v is the variation of curvature of the plastic strain region (crystal lattice curvature of the region) due to applied load, tn is the tangent, s is the running length of the region, b1 and b2 are the ‘‘Burgers vector” magnitudes of bulk translation and subsurface or rotation incompatibility, respectively, and rf is the gradient part of the flow due to external sources. The channeled spiral pattern of localized defect flow with varying curvature v is illustrated in Fig. 1. It is seen from Fig. 1(a) that with a low curvature v, the transverse strain rate v is low, and the spiral experiences slight torsion characterized by a large transverse wave length. This situation corresponds to the development of plastic shear in slightly nonequilibrium surface layers, low angle interfaces in a solid. As the plastic flow J curvature increases, the transverse wave length decreases greatly, whereas the transverse strain rate increases (Fig. 1(b)). With high curvature v, hydrostatic tension and structure–phase decomposition of material develop in these regions, giving rise to cracking. Experimental studies of surface mesobands of localized deformation propagating as spirals are reported in [5,16]. Examples of these spiral mesobands are given in Fig. 2. The spiral subgrain boundaries are readily observable in plastic deformation of highly curved crystal structure with low shear stability. Fig. 2(a) shows a spiral mesoband formed in Al foil in alternate bending after N = 4  105 cycles [16]. A plenty of surface shear bands, generated as spirals from the stress concentrates C, are seen in Fig. 2(b) [5]. 3. Nucleation of strain-induced defects in solids by laser pumping mechanism The multiscale description of a deformable solid as a nonlinear hierarchically organized system managed to solve the problem of nucleation of strain-induced defects. The distribution of normal and shear stresses at the interface between 3D-crystalline and 2D-planar subsystems has a periodic pattern [17]. As a consequence, nonlinear plastic shears within a planar subsystem make up the concentration fluctuation of ions localized in zones of tensile normal stresses (the Bose-condensate [9]), Fig. 3. Such ion clusters are screened by electron gas from neighbor interstitial space. It causes the formation of very high local curvature within the external layers of adjacent crystalline grains where the new allowed structural states appear in interstitial spaces (see below). It is very important that a new electron fluctuation band below the Fermi level appears in the electron spectrum corresponding to the new structural states in interstitial space. In these conditions highly excited ions and electrons in planar clusters are to be transferred to new low-energy states in interstitial space as a straininduced defect in a crystal. It is the general mechanism of any strain-induced defect nucleation similar to laser pumping. If there is no electron fluctuation band in a solid electron spectrum related to the curvatured interstitial space the development of local curvature causes brittle fracture. The kind of a strain-induced defect nucleated by the laser pumping mechanism depends upon the size of local crystal lattice curvature zone related to fluctuation of ion concentration in a 2Dplanar subsystem. In curvature zones of a small size, dislocations

Please cite this article in press as: V.E. Panin et al., Comput. Mater. Sci. (2015), http://dx.doi.org/10.1016/j.commatsci.2015.10.045

V.E. Panin et al. / Computational Materials Science xxx (2015) xxx–xxx

3

Fig. 1. Localized plastic deformation configuration and velocity as a function of the curvature v of the deformed region: v1 < v2 [13].

Fig. 2. Spiral nonlinear wave subgrain boundary as local structure–phase decomposition in bended aluminum foil, alternate bending, N = 4  105 cycles, T = 293 K [16] (a); generation of nonlinear waves of surface shear bands from the angular point C in Pb + 1.9% Sn polycrystal alloy; tension, V = 1.1%/min, T = 543 K, e = 30%, 1000 [5] (b).

understanding of polycrystalline solid plasticity at very low or high temperatures of loading. At last for very large local crystal lattice curvature zone, the number of interstitial bifurcation vacancies becomes so high that they form a crack and initiate fracture. This process is classified as the structural phase collapse of a loaded solid [21].

4. Interstitial bifurcation vacancies at the interface between a 2D-planar and a 3D-crystalline subsystems Fig. 3. Laser pumping mechanism of the Bose-condensate formation in the a, b, c, d zones of the planar subsystem AB which act as the sources of dislocation generation in a 3D-crystal [18].

nucleate. If the curvature zones are of middle size, a crystal under loading generates shear bands which are translation–rotation strain-induced defects. They transform an elastic curvature zone of a solid into plastically misoriented one. The stresses in a large extended curvature zone in a crystal under loading can be relaxed by the disclination nucleation. The new kind of nonequilibrium vacancies named as interstitial bifurcation vacancies plays dominant role in this process [19]. They allow forming a disclination as a strain-induced defect of high misorientation. Well-known pile-up analysis of the Hall–Petch equation r = r0 + kd1/2 does not take into account the primary role of nonlinear wave flows within grain boundaries of a deformed polycrystal. If it will be done, the Hall–Petch equation will take the form r = r0 + k1d1/2 + k2d1 [20]. It is very important result for

It should be emphasized that in the vicinity of interface zones where positive ion clusters are formed, bifurcation structural states (BSS) not found in the initial crystal appear in interstitial space of curvatured crystal structure [21,22]. Fig. 4 illustrates the scheme of BSS formation [19], which is based on the very important result obtained for a onedimensional crystal [22]. According to [22], with a unique minimum of the pair interaction potential of atoms, an increase in local interatomic spacing gives rise to bifurcation minima of the potential of the particle system. This point is of fundamental importance for the validation of the mechanism of formation of strain-induced defects in a curvatured crystal lattice at the interface between a 2D-planar and a 3D-crystalline subsystems. As is shown in Fig. 3, flows of structural transformation at grain boundaries form clusters of positive ions in zones of tensile normal stresses. An excess positive charge in grain-boundary clusters should be screened by electron gas of adjacent 3D-crystalline layers (zone 2 in Fig. 4).

Please cite this article in press as: V.E. Panin et al., Comput. Mater. Sci. (2015), http://dx.doi.org/10.1016/j.commatsci.2015.10.045

4

V.E. Panin et al. / Computational Materials Science xxx (2015) xxx–xxx

Fig. 4. Schematic of formation of the interstitial bifurcation states (2) in a crystal lattice curvature region.

The interatomic spacing in zone 2 given in Fig. 4 increases resulting in the formation of interstitial bifurcation structural states. The number of degrees of freedom grows sharply in near-boundary regions of the polycrystal. These interstitial structural states can play an important role of athermal vacancies in any kind of mass transfer. Cluster ions in a 2D-planar subsystem can transfer into interstitial structural states of 3D-nearboundaries layers generating the cores of strain-induced defects. This phenomenon must be taken into account at realization of additive manufacturing technologies of any kind (see below).

5. Excitable cellular automata simulation The ultimate aim of the paper is computer simulation for revealing the role of dynamic curvature of flows to propagate along grain boundaries in a deformed polycrystal. Nowadays there are lots of deep and comprehensive studies available that emphasize an important role of rotational modes of deformation [23–25]. In most of them rotations of material fragments are discussed. In the present paper the role of rotations is considered in terms of energy flow curvature. The multiscale excitable cellular automata approach developed allows simulating transfer of elastic energy perturbations from one scale (entire specimen surface) to another (grain boundaries) through explicitly accounting for local moments of forces and angular velocities of flows. The hydrodynamic notion of an active element as a control volume introduced by the authors allows the simplest representation of near-boundary deformation processes in the form of flows with one or another degree of vorticity. The kernel of the computation algorithm used in the excitable cellular automata method represents the set of relations, based on ‘‘extended Tornbull theory” [26,27]. In the framework of this method the main relation of linear velocity of elastic energy transfer is written as follows:

v ¼ @  1;

ð4Þ

where @ is elastic front propagation mobility, 1 is normal stress. Fig. 5 illustrates the calculation scheme for the angular velocity of the material in an active element at vortex formation on an inhomogeneous field of flows between adjacent volumes of the element environment (for simplicity, the scheme is represented in the plane approximation). The total angular velocity of the i-th element in Fig. 5 under the action of the flux of matter through the boundary of the k-th and lth elements (each k-th element lies in the first coordination sphere of the i-th, and each l-th lies at the intersection of the first coordi-

Fig. 5. The scheme of angular velocity calculation: 2D simplest representation.

nation spheres of the i-th and a corresponding k-th elements) is defined as follows:

xi ¼

K X L X rikl  v kl k¼1 l¼1

jrikl j2

:

ð5Þ

Here, K is the number of elements in the first coordination sphere of the i-th element, and L is the number of elements at the intersection of the first coordination spheres of the i-th element and each k-th neighbor. The change of the moment of force of the i-th element for time s is calculated in the following way:

DMi ¼

Gpr 3c Dci : 2

ð6Þ

Here, G is the shear modulus of the material contained in the ith element, rc is the element radius, and Dci is the threedimensional angular displacement, which is proportional to the total angular velocity: Dci = xis. Hybrid bistable/excitable Cellular Automata (CA) method allows calculating a number of interstitial bifurcation vacancies (IBV) on active element when local material torsion occurs (see Fig. 6). As one can see at Fig. 7, high curvature crystal lattice region gives rise to IBV formation. Simulation of mass transfer by mecha-

Please cite this article in press as: V.E. Panin et al., Comput. Mater. Sci. (2015), http://dx.doi.org/10.1016/j.commatsci.2015.10.045

V.E. Panin et al. / Computational Materials Science xxx (2015) xxx–xxx

5

Fig. 6. Scheme of simulation of interstitial bifurcation vacancies.

previous time step, then the BCA element with ‘‘IBV”-state does not change its state and elements with ‘‘LS”-state switch stochastically into ‘‘IBV” one so that required IBV density value is reached. At each time step of BCA algorithm ‘‘IBV”-states is redistributed among the elements so that interstitial bifurcation vacancies could accumulate and form defect of higher scale (i.e. crack or pore). The algorithm of redistribution of such kind consists of two spatial loops. In the first one all BCA elements are searched for ‘‘IVB”state. When j-th element is found, the number of ‘‘IBV” neighbor elements (d1) is calculated. If d1 > 0 then the states of the neighbor elements do not change. If d1 = 0 then number of elements in 2-nd (d2) and 3-rd (d3) coordination spheres is calculated. When d1 = 0, d2 = d3 = 0 then probability of element switching into ‘‘LS”-state is increasing function of absolute value of rotation solid angle (c). In this case of switching one of the neighbors located in 1-st coordination sphere (k-th element) should transit into ‘‘IBV”-state. In other words, this pair exchanges the states. Probability of such exchange (Pjk) is calculated for each k-th element in first coordination sphere. This probability is higher the closer is c direction to one of the radius-vector between given element and its neighbor ujk. Moreover, orientation of initial crystal lattice (a, b, c) is taken into account as well. So the probability function could be written as follows:

Pjk ¼ pmax

jcj 2 2 sin ð/; cÞðcos4 ðu; aÞ þ cos4 ð/; bÞ þ cos4 ð/; cÞÞ : 2p

ð7Þ

Illustration of these rules is shown in Fig. 7. When d1 = 0, d2 > 0 or d3 > 0 then switching probabilities of given element into the ‘‘LS”-state and k-th element in the first coordination sphere into the ‘‘IBV”-one are defined by the rules for abovementioned case (d1 = d2 = d3 = 0) with taking into account distance rjl between given j-th element and the element in ‘‘IBV”state (l-th) contacting with k-th element and located in second or third coordination sphere. The less is rjl the higher is probability Pjkl. Moreover the larger is number of elements in 1-st coordination sphere contacting with l-th element (Kjl) the lower is Pjkl (if Kjl = 0 then Pjkl = 0): Fig. 7. The scheme of IBV redistribution in the BCA network.

nism of their propagation in planar subsystem can be realized on the basis of bistable cellular automata method. Elements of bistable cellular automaton is to have two possible states – ‘‘lattice site” (‘‘LS”) and ‘‘IBV”. Switching rules of each CA element depend on neighbor element state. Bistable Cellular Automata method (BCA) is embedded into Excitable Cellular Automata method as follows. Each element of ECA contains a network of BCA that simulates material behavior at lower scale level. At the end of each time step of ECAalgorithm lower-scale timing loop of BCA-simulation is involved. At the first BCA time step the distribution of the states of BCAelements is defined starting from the IBV density in a given ECAelement. When IBV density becomes larger in comparison with

Pjkl ¼

r max p r jl K jl max

jcj 2 2 sin ð/; cÞðcos4 ðu; aÞ þ cos4 ð/; bÞ þ cos4 ð/; cÞÞ : 2p X XX Pjk þ Pjkl : Pj ¼ 

k

k

ð8Þ ð9Þ

l

The sum value of switching probability for the j-th bistable element (Pj) is written as follows. After calculation over all switching probabilities the second spatial loop of the algorithm starts. In the framework of this loop it is determined whether some element should switch into ‘‘LS”state in terms of its own calculated probability. If the switching occurs then one of the elements in 1-st coordination sphere should transit into ‘‘IBV”-state. This process is stochastic and defined by

Please cite this article in press as: V.E. Panin et al., Comput. Mater. Sci. (2015), http://dx.doi.org/10.1016/j.commatsci.2015.10.045

6

V.E. Panin et al. / Computational Materials Science xxx (2015) xxx–xxx

Fig. 8. Schematic representation of disclination forming in triple joint (upper picture) and simulation of angular velocity distribution along grain boundaries of different types (lower picture).

the probability distribution {Pk}, which is determined in the first spatial loop of the algorithm:

Pk ¼

X XX Pjk þ Pjkl ; j

j

ð10Þ

l

where j is responsible for each ‘‘IBV” element in the 1-st coordination sphere of given k-th element, l is responsible for ‘‘IBV” element in intersection region of the 1-st coordination sphere of k-th element and the 2-nd coordination sphere of j-th element. Therefore the BCA elements tend to accumulate into conglomerates (clusters). Total number of these clusters and their spatial distribution define value of latent defect energy. This latent defect energy is responsible for formation of decelerating fields on the path of elastic flows and crack nucleation. Once the element is switched into ‘‘crack” state, it transfers its total energy to neighbor elements and could not exchange it anymore. Thus, taking into account generation of interstitial bifurcation vacancies, the value of elastic energy mobility between i-th and k-th elements is calculated as follows:

@

n1 ik



¼ ð@0 Þik  e

ðHik þQ n1 Þ ik kB T n1 ik

;

ð11Þ

where ð@0 Þik – initial (‘‘elastic”) mobility, upper indexes are responsible for time step; – temperature in the boundary between i-th and k-th Tn1 ik element; Hik – boundary energy; Qik – latent defect energy;

Q n1 ik

¼ ðg

n1 i

þg

n1 n1 k ÞkB T ik ln

! Nn1 þ Nn1 i k ; gn1 þ gn1 i k

ð12Þ

gn1 , gn1 – number of interstitial bifurcation vacancies in i-th i k element and its k-th neighbor; Nn1 , Nn1 – number of atoms in i-th element and its k-th i k neighbor.

specimen dimension was 10  60  60 lm3; the size of a cellular element was 1 lm; the total time was 100 ls; the time step was 1 ns. The upper face was assigned to have a constant compressive strain rate of 0.03 s1. In this numerical experiment the distributions of angular velocity components and rotational strain component were calculated over the interfaces of the Al polycrystal deformed by axial compression. Fragments of the polycrystal with a triple grain junction were investigated thoroughly. The result of simulation is shown in Fig. 8. In complete accordance with the misorientation of the adjacent grains, their boundaries are involved in development of rotational waves of plastic flow and associated rotational deformation. At the boundaries of adjusted grains rotational deformation modes of grain boundary flows with widely differing amplitudes evolve. Low-amplitude rotation waves develop at the low-angle boundary of grains. High-amplitude rotational waves develop at the highangle boundary of grains. The angular velocity modulation amplitude for the wave propagating lengthways the low-angle boundary is rather small. The constant angular velocity component, which characterizes the linear velocity of the rotation wave lengthways the low-angle boundary, remains high. An important result in the simulation of angular velocities of grain boundary flows is found when they approach the triple junction of misoriented grains. For both low- and high-angle boundaries, the angular velocities of flows on either side of the boundaries of adjacent grains assume the same sign (Fig. 8). However, the signs of angular velocities at boundaries 1–2 and 1–3 near the grain triple junction are found opposite. This superposition of the angular velocities governs the generation of a disclination in the triple junction region. In the simulation of grain boundary flows considered above, the only varied parameter was the degree of misorientation of adjacent grains in a polycrystal. Naturally, of great significance is accounting for crystallographic orientation of each of the adjacent grains. It is quite clear that engineering of grain boundaries might be developed on the basis of construction of multiscale models of deformed polycrystals. It will be considered in details in our next paper.

In different way mobility could be written as follows:

ð13Þ

6. Experimental and numerical simulation evidences of the role of micromechanical processes to develop at grain boundaries and interfaces

The excitable cellular automata method was used in numerical experiment on compression of an aluminum specimen with a grain triple junction and the source located on the upper face. The

The effect of grain boundary sliding on the generation of various strain-induced defects was studied under tensile tests performed at room temperature with flat polycrystalline specimens made of

@n1 ¼ ð@0 Þik  ik

gn1 þ gn1 i k Nn1 þ Nn1 i k

!gn1 þgn1 i

k

e



Hik

kB T n1 ik

:

Please cite this article in press as: V.E. Panin et al., Comput. Mater. Sci. (2015), http://dx.doi.org/10.1016/j.commatsci.2015.10.045

V.E. Panin et al. / Computational Materials Science xxx (2015) xxx–xxx

7

Fig. 9. Surface topography of a flat Pb + 0.03% Te alloy specimen, optical micrograph, 125 (a); near-boundary zone bd of material delamination and layer by layer intrusion; high magnification 1850 (b).

Fig. 10. Mass transfer of liquid Ti into ZrO2 along high-angle grain boundaries, being enhanced by athermal interstitial bifurcation vacancies.

Pb + 0.03% Te alloy. A pattern of inhomogeneous plastic flow in a surface layer of a specimen at the strain e = 15% is given in Fig. 9. Fig. 9a is illustrative of intensive grain boundary sliding along boundary ab of coarse non-equiaxed grain A, being a site of generation of dislocation shears and moment producing a counterclockwise rotation of grain A. At fixed triple junction D, this grain demonstrates an extended zone of constrained deformation resulting in formation of the disclination be and the counter stress field rf which directs the dislocation sliding a along the grain boundary ab being opposite to the grain boundary sliding. Hindering of rotation of non-equiaxed grain A in the triple junction region produces a dramatic weakening effect on grain boundary sliding along boundary bd, which is well seen in Fig. 9b where boundary bd is retained without visible shears. Rotation of grain A regarding boundary bd occurs in the near-boundary region where the material experiences delamination due to severe lattice curvature. This is a result of the formation of interstitial bifurcation vacancies represented in Fig. 4. Thus, experimental mechanisms of complex development of plastic deformation in the surface layer of a polycrystal are in good agreement with their simulation demonstrated in Fig. 8.

Grain boundary mass transfer of liquid Ti into ZrO2 stimulated by athermal interstitial vacancies was investigated experimentally under process to develop at electron-beam assisted powder interaction to take place in additive manufacturing technology as well as by means of multiscale computer simulation using the hybrid discrete–continuous method of excitable cellular automata. Fig. 10 exhibits a cross section of the layered metal–ceramic composite ZrO2–Ti formed with the use of electron-beam additive manufacturing technology. In doing so, liquid Ti is deposited onto a flat ZrO2 specimen with a pulse scanning electron beam in a vacuum chamber. It is seen, that titanium is penetrated along highangle boundaries between ZrO2 grains to form a multilayer phase Ti2ZrO located behind the front of its propagation, and crystalizes in a form of a thin titanium film at the composite’s surface. Such three-layer metal–ceramics possesses high mechanical characteristics, which is due to originally structured interfaces between ceramics and metal. The interface structure can be simulated using the hybrid discrete–continuous method of excitable cellular automata, Fig. 11. In numerical experiments a specimen is simulated by cellular automaton with 100 nm active element size. Size of the specimen was equal to 4 lm  1 lm  6 lm. Upper melt layer was taken as 1 lm in width; width of ‘‘polycrystalline” Ti layer was equal to 4 lm and width of ZrO2 substrate was equal to 1 lm. Using special CA algorithm (see part 5) grain structure with size of 1 lm was introduced into the intermediate Ti layer (Fig. 11a). Grain size of ZrO2 substrate was taken larger than 6 lm; in doing so the lower layer represented one single grain. Initial temperature of whole specimen was 300 K (at the zero time step). At the first time step of the algorithm the temperature in upper melt layer was set to 2000 K. In the bottom of the specimen (lower face) the temperature was set to 300 K at all-time steps, in doing so the cooling was simulated there. Time step was 0.01 ns, full time of thermal load was equal to 1 ls. The numerical experiment revealed rotational wave flows generation along intergranular regions in Ti polycrystal (Fig. 11a and b). According to the theory [12,13] (see part 2 of this paper), it is obvious that in this case mass transfer of Zr along with oxygen from ZrO2 substrate at grain boundaries will occur. This will give rise to dissolution of Ti grains with some Ti–Zr–O phase nucleation. In doing so, one of the most substantial reasons of intensive structural transformation during processes to develop during additive (sintering) technology is rotational wave character of mass transfer.

Please cite this article in press as: V.E. Panin et al., Comput. Mater. Sci. (2015), http://dx.doi.org/10.1016/j.commatsci.2015.10.045

8

V.E. Panin et al. / Computational Materials Science xxx (2015) xxx–xxx

Fig. 11. ECA-simulation of pulsed heat loading of ‘‘liquid Ti”–‘‘polycrystal Ti”–ZrO2: structure of specimen simulated (a), distribution of X-component of local moment of force (b), distribution of Z-component of local moment of force (c).

7. Summary In this paper a loaded solid is considered to consist of two subsystems, namely, a 3D translation-invariant crystal and a 2D planar subsystem comprising surface layers and all interfaces. The 2D planar subsystem is characterized by a cluster structure and absence of translation invariance. Primary plastic shears in a solid under loading are developed within the 2D planar subsystem as nonlinear waves. A very important parameter of nonlinear wave flow is the curvature of crystal structure within the 2D planar subsystem. A phenomenological model is discussed that suggests that clusters of positive ions are formed in the 2D planar subsystem to generate local zones of lattice curvature in external layers of adjacent 3D grains. These zones are characterized by allowed structural states in interstices to which cluster ions pass from the 2D planar subsystem to the 3D crystal, forming cores of straininduced defects. A particular type of strain-induced defects to occur there (athermal vacancies, dislocations, shear bands, disclinations, etc.) is determined by the cluster size in the 2D planar subsystem and the electron structure of the 3D crystal. A modified method of excitable cellular automata is developed with consideration done for the local moments of forces. The method enables to simulate propagation of nonlinear waves of structural transformation at grain boundaries with various lattice curvature within gradient fields of internal stresses. Propagation of nonlinear waves of structural transformations at grain boundaries of a polycrystal in the triple junction region was simulated. A possibility of spiral-like propagation of nonlinear waves at high-angle grain boundaries and generation of disclinations within a 3D grain in the triple junction region was shown by means of computer simulation. The simulation data agrees well with experimental evidences on nonlinear waves propagation at grain boundaries of the Pb + 0.03%Te polycrystal in the triple junction region. Another experimental result to prove the idea of nonlinear wave flows to develop at interfaces was shown for the case of EB-stimulated formation of metal–ceramic composite (ZrO2–Ti). In doing so, the revealed data are considered as manifestation of grain boundary and interphase flows to play an important role both at structure formation and under deformation and fracture of polycrystals.

2013–2020, RFBR project No. 14-01-00789 as well as grant of the President of the Russian Federation for support of leading scientific schools No. NSh-2817.2014.1. References [1] A. Needleman, R.J. Asaro, J. Lemonds, D. Peirce, Comput. Methods Appl. Mech. Eng. 52 (1985) 689–708. [2] G.P. Cherepanov (Ed.), Fracture: A Topical Encyclopedia of Current Knowledge, Krieger Publ. Comp., Malabar, Florida, 1998. [3] W. Skrotzki (Ed.), Abstract Book of the 15th Int. Conf. on the Strength of Materials, Techn. Univ. Dresden, Germany, 2009. [4] A. Carpinteri (Ed.), Abstract Book of the 11th Int. Conf. on Fracture, Politecnico di Torino, Italy, 2005. [5] V.E. Panin (Ed.), Physical Mesomechanics of Heterogeneous Media and Computer-Aided Design of Materials, Camb. Internat. Science Publ., Cambridge, 1998. [6] G.C. Sih, Phys. Mesomech. 13 (5–6) (2010) 233–244. [7] V.E. Panin, V.E. Egorushkin, A.V. Panin, Phys. Usp. 55 (12) (2012) 1260–1267. [8] V.E. Panin, V.E. Egorushkin, A.V. Panin, Phys. Mesomech. 15 (3–4) (2012) 113– 146. [9] V.P. Maslov, Teoret. Mater. 159 (1) (2009) 20–23. [10] V.E. Panin, V.E. Egorushkin, T.F. Elsukova, Phys. Mesomech. 16 (1) (2013) 1–8. [11] L.D. Landau, E.M. Lifshits, Theory of Elasticity, Nauka, Moscow, 1987 (in Russian). [12] V.E. Egorushkin, Russ. Phys. J. 33 (2) (1990) 135–149. [13] V.E. Egorushkin, Russ. Phys. J. 35 (4) (1992) 316–338. [14] E. Kroner, Gauge Field Theories of Defects in Solids, Max Planck Inst., Stuttgart, 1982. [15] A. Kadicˇ, D.G.B. Edelen, A gauge theory of dislocations and disclinations, Lecture Notes in Physics, vol. 174, Springer, Berlin, 1983. [16] V.E. Panin, N.S. Surikova, T.F. Elsukova, V.E. Egorushkin, Yu.I. Pochivalov, Phys. Mesomech. 13 (3–4) (2010) 103–112. [17] G.P. Cherepanov, J. Appl. Phys. 78 (11) (1995) 6826–6832. [18] V.E. Panin, V.E. Egorushkin, Fundamental role of local curvature of crystal structure in plastic deformation and fracture of solids, in: 2014 AIP Conf. Proc. Physical Mesomechanics of Multilevel Systems, vol. 1623, 2014, pp. 475–478. [19] V.E. Panin, A.V. Panin, T.F. Elsukova, Yu.F. Popkova, Phys. Mesomech. 18 (2) (2015) 89–99. [20] V.E. Panin, R.W. Armstrong, Phys. Mesomech. 18 (3) (2015) 5–10. [21] V.E. Panin, V.E. Egorushkin, Phys. Mesomech. 16 (4) (2013) 267–286. [22] M.A. Guzev, A.A. Dmitriev, Phys. Mesomech. 16 (4) (2013) 287–293. [23] W.H. Mueller, An Expedition to Continuum Theory, Springer, 2014. [24] A.Yu. Smolin, E.V. Shilko, S.V. Astafurov, I.S. Konovalenko, S.P. Buyakova, S.G. Psakhie, Defence Technol. 11 (2015) 18–34. [25] S.G. Psakhie, K.P. Zolnikov, A.I. Dmitriev, A.Yu Smolin, E.V. Shilko, Phys. Mesomech. 17 (1) (2014) 15–22. [26] D.D. Moiseenko, Yu I. Pochivalov, P.V. Maksimov, V.E. Panin, Phys. Mesomech. 16 (3) (2013) 248–258. [27] V.E. Panin, D.D. Moiseenko, S.V. Panin, P.V. Maksimov, I.G. Goryacheva, C.H. Cheng, J. Appl. Mech. Tech. Phys. 55 (2) (2014) 318–326.

Acknowledgments The study was performed under support of fundamental Research Program of the State Academies of Sciences for

Please cite this article in press as: V.E. Panin et al., Comput. Mater. Sci. (2015), http://dx.doi.org/10.1016/j.commatsci.2015.10.045