## Abstract

In this paper, we present the development of a concurrent atomistic–continuum (CAC) methodology for simulation of the grain boundary (GB) structures and their interaction with other defects in ionic materials. Simulation results show that the CAC simulation allows a smooth passage of cracks through the atomistic–continuum interface without the need for additional constitutive rules or special numerical treatment; both the atomic-scale structures and the energies of the four different [001] tilt GBs in bi-crystal strontium titanate obtained by CAC compare well with those obtained by existing experiments and density function theory calculations. Although 98.4% of the degrees of freedom of the simulated atomistic system have been eliminated in a coarsely meshed finite-element region, the CAC results, including the stress–strain responses, the GB–crack interaction mechanisms and the effect of the interaction on the fracture strength, are comparable with that of all-atom molecular dynamics simulation results. In addition, CAC simulation results show that the GB–crack interaction has a significant effect on the fracture behaviour of bi-crystal strontium titanate; not only the misorientation angle but also the atomic-level details of the GB structure influence the effect of the GB on impeding crack propagation.

## 1. Introduction

The interaction of grain boundaries (GBs) with other defects plays a significant role in the overall properties of crystalline materials, such as ductility, strength and hardness [1,2]. The dynamics of defects usually evolve over a wide range of length scales, making the interaction between defects an inherently multi-scale phenomenon. For example, the effect of the GB–defect interaction on mechanical behaviour may depend on not only the microstructural details of the interaction sites [3,4] but also the long-range interaction of the GB with surfaces and other defects in the material. Understanding the nature of this effect thus necessitates a simulation of the material with real size and its real microstructure. Such a length scale is beyond the reach of current atomistic methods-based simulation tools. By reducing the number of degrees of freedom (DOF) in regions far away from critical regions such as the GB–defect interaction sites so as to increase the size of the computational models, concurrently coupled atomistic–continuum methods are advantageous in simulation of crystalline materials over fully atomistic methods because the former take into consideration the long-range interaction of GBs with other defects while retaining atomic details in critical GB regions.

Over past decades, there have been a number of concurrent atomistic–continuum (CAC) methods developed to investigate the interaction between GBs and other defects. Among the popular methods, the static quasi-continuum (QC) method was used to simulate the interaction of GBs with dislocations or cracks by Shenoy *et al.* [5], in which the dislocation pile-up was formed during the interaction of dislocation with the

Despite two decades of intense research efforts, however, most of the existing multi-scale methods have only been applied to the study of monatomic crystalline materials such as Cu and Al. The classical methods for coarse-graining monatomic crystals to find the elastic constants are the methods of long waves or the method of homogeneous deformations, an example of which is the Cauchy–Born approximation employed in the formulation of the QC method [9–11]. There have been interesting attempts to extend the QC method to simulate materials with multiple atoms in the primitive unit cell. However, it should be noted that an additional minimization is needed to obtain the internal position of each atom within the unit cell after obtaining the displacement of unit cells in the modified QC method [12] and QC-based method [13].

For polyatomic crystals, i.e. materials that have more than one atom in the primitive unit cell, the existence of inner displacements or internal strains increases the mathematical complexity of the formulation of the methods. As a consequence, most of the existing coupled or mixed atomistic–continuum methods that use atomistic dynamics at the nanoscale and classical elasticity at the mesoscale were formulated for monatomic crystals and hence are not applicable to polyatomic crystalline solids. Thus, these methods are not suitable for ionic crystalline materials as all ionic materials are necessarily polyatomic materials with more than one atom in the primitive unit cell. Although coarse-grained molecular dynamics (CGMD) [14] can simulate polyatomic materials in principle, extending CGMD to simulate the dynamic evolution of defects would be very difficult because the stiffness and mass matrices need to be pre-computed. Until now, there has been no successful application of other coupled atomistic–continuum methods to the simulation of defect interaction in ionic materials.

The objective of this paper is to present a new development of the CAC method for the study of ionic materials with GBs and other defects and then to demonstrate the method through simulation of the dynamic processes of crack initiation, propagation and interaction with GBs in strontium titanate (SrTiO_{3}; STO). STO is a technologically important polyatomic crystalline material that has a perovskite structure. Perovskite oxides in general have a broad application and STO in particular is one of the most commonly used substrates and a potential thermoelectric material [15]. Despite the important role of GBs and other defects in the reliability and the functionality of STO [16,17], current understanding of the structure–property relationship in this material is limited. For example, unlike the GB structures in metals, the GBs in ionic materials [18] have an open structure, in which the density of the GB is significantly less than that of the bulk. Thus, the established understanding of the effect of GBs in metals may not still be applicable to ionic materials. For STO, the atomic structures of several tilt GBs [19,20] have been obtained using experimental methods and density function theory (DFT)-based calculations. Several tilt GB structures of STO such as

In this work, we focus our study on bi-crystal models for two reasons. First, experimental measurements and DFT calculation results of bi-crystals are available and can be used to test the applicability of the CAC method in reproducing the atomic-scale details of the GB structure. Second, the bi-crystal configuration allows for the study of specific GB structures and the GB structure–property relationship, and hence can be used to identify the mechanism of GB–defect interaction and the effect on the dynamic fracture behaviour. This paper is organized as follows: following the introduction, the CAC methodology will be introduced in §2; the computational models and simulation results of bi-crystal models containing cracks and four different tilt GBs will be presented in §3; comparisons between the CAC and DFT results of the atomic structures of different GBs and those between the CAC and molecular dynamics (MD) simulation results of the GB–crack interactions will also be presented. This paper will end with a summary and discussion in §4.

## 2. The concurrent atomistic–continuum methodology

The CAC methodology combines (i) a continuum field representation of microscopic balance equations for general crystals with fully atomistic information built into the formulation and (ii) a modified finite-element (FE) method employing primitive unit cell-shaped three-dimensional solid elements. The theoretical formulation follows the formalism developed by Kirkwood [22], Irving & Kirkwood [23] and Hardy [24], and is an extension of Irving and Kirkwood's statistical mechanics procedure [23] to a two-level structural description of crystalline materials [25–27]. This two-level description results in a link between the continuously distributed local density function ** a**(

**,**

*x*

*y*^{α},

*t*) in the

*physical space*and the

*phase space*dynamical function

**(**

*A*

*R*^{kξ},

*V*^{kξ}),

*R*^{k}are the positions of the mass centre of the

*k*th unit cell, Δ

*r*^{kξ}are the internal positions of the

*ξ*th atom in the unit cell

*k*,

*R*^{kξ}=

*R*^{k}+Δ

*r*^{kξ},

**is the physical space coordinate and**

*x*

*y*^{α}(

*α*=1,2,…

*ν*, where

*ν*is the number of the atoms in a unit cell) is the internal variable describing the internal position of atom

*α*embedded within the lattice cell at

**; the first localization function,**

*x**δ*(

*R*^{k}−

**), implies that over the entire physical space all the unit cells can be found; then, for each unit cell**

*x**k*, the second localization function,

*r*^{kξ}to be

*y*^{α}. The new formulation describes the structure of a crystalline material in terms of continuously distributed lattice cells, but with a group of discrete atoms situated in each lattice cell, in exactly the same way as the solid-state physics approach in describing the structure of all crystals.

As exact consequences of the new definition of local density functions and Newton's laws, Chen [27] obtained the field representation of microscopic balance laws
*ρ*^{α}, *ρ*^{α}(** v**+Δ

*v*^{α}) and

*ρ*

^{α}

*e*

^{α}are the local density of mass, the linear momentum and the total energy, respectively;

**+Δ**

*v*

*v*^{α}=

*ρ*

^{α}(

**+Δ**

*v*

*v*^{α})/

*ρ*

^{α}is the atomic-level velocity and

**is the velocity field;**

*v*

*f*^{α}is the external force field;

*t*^{α}and

*q*^{α}are, respectively, the momentum flux and heat flux due to the homogeneous lattice deformation and are consistent with the Cauchy stress and heat flux concepts in homogenized continua;

*τ*^{α}and

*j*^{α}are, respectively, the momentum flux and heat flux due to the reorganizations of atoms within the lattice cells and represent the components of stress and heat flux due to atomic-scale inhomogeneities.

Following the classical definitions of kinetic temperature and the internal force density [25,26], the instantaneous balance equation of linear momentum has been rewritten as [28]
*m*^{α} being the mass of the *α*th atom, Δ*V* (**x**) the volume of the unit cell located at ** x**,

**u**

^{α}(

**x**) the displacement of the

*α*th atom embedded within the unit cell located at

**,**

*x**γ*

^{α}=

*m*

^{α}/

*MΔV*,

*k*

_{B}Boltzmann constant. Generally, the internal force density

For systems with a given homogeneous (constant) temperature field, the temperature would behave like a surface traction on the boundaries. At the coarse-grained level, with a homogeneous temperature assumption, equations (2.5) or (2.6) can only account for the effect of thermal expansion or contraction. However, as the value of the Boltzmann constant *k*_{B} (1.3806488×10^{−23} m^{2} kgs^{−2} K^{−1}) is extremely small, such a temperature effect is essentially negligible. In this work, the effect of a homogeneous temperature effect is thus ignored.

To numerically solve the governing equation, equation (2.6), we discretize a material body with finite elements consisting of a collection of unit cells and approximate the displacement field within an element by
*Φ*_{ξ}(**x**) is the shape function, and *α*th atom within the *ξ*th node. The weak form of equation (2.6) can then be obtained as

The reformulated balance laws have a continuum framework at the unit cell level but retain the atomistic information, including the arrangement and interaction of atoms, at the sub-structural level. Solutions of the balance equations include both the lattice deformation and the rearrangement of atoms relative to the lattice; hence, the Cauchy–Born assumption of affine deformation is not invoked. Note that, in most of the existing multi-scale methods [29], the description of the material is a continuous collection of material points, i.e. unit cells or atoms. However, there is no net charge for each unit cell in ionic materials; hence it is impossible for such kinds of models to fully consider the coulombic interactions if the ions embedded within each material point are not modelled explicitly. This is the reason why most existing concurrent multi-scale methods are not applicable to simulation of ionic materials and why it is necessary for each continuum point to further have internal DOF to describe the motion of ions of different charges in order to describe the ionic interaction in the material. Thus, the CAC method is advantageous over the existing multi-scale methods in the simulation of ionic materials.

In distinction from many existing multi-scale methods [29], the CAC method is based on a field representation of a complete set of balance laws, which makes it applicable to dynamic and non-equilibrium processes with the force field or interatomic potential being the only empirical input. MD of the underlying atomistic system can be achieved by pursuing discretization of balance equations in the limit to primitive unit cells in the critical region nearby, for example GBs. Assuming a collective motion of unit cells, the majority of the DOF of an atomistic system can be eliminated in the far region away from the GBs. This naturally leads to a CAC methodology [28,30], in which both the atomistic and continuum domains are governed by the same set of balance equations. One can adaptively remesh, if desired, to ‘zoom in’ or ‘zoom out’ in terms of the resolution of the underlying atomistic system; in the meantime, the CAC method allows the defects to pass from the atomistic to the continuum domain and vice versa without additional passing techniques or adaptive mesh refining to the atomic scale [31,32].

In this work, the crystalline material body is discretized into finite elements, with each element containing a collection of primitive unit cells. In the simulations conducted in this work, standard Gauss integration [32] is employed, along with linear interpolation within each element to reduce the computational cost. The central-difference method is used to solve the discretized dynamic equations. As a consequence of the analytical equivalence between the divergence of the stress tensor and the internal force density, the spatial derivatives in the formulation are eliminated. To compare the accuracy and efficiency of CAC with all-atom MD simulations, the general-purpose parallel MD code LAMMPS [33] is employed for the latter. The numerical implementation code of CAC is parallelized and run on multiple processors.

## 3. Computational models and simulation results

### (a) Interatomic potential and simulation set-up

There have been a number of interatomic potential models developed for STO. According to the studies by Thomas *et al*. [34], Benedeck *et al*. [35] and Hirel *et al*. [36], the rigid-ion potential developed by Thomas *et al*. [37] for STO behaves even better in comparing with experimental and *ab initio* data of elastic constants, lattice parameter, formation energy of defects and stacking faults energy than some shell-model potentials. Hirel *et al*. [38] successfully used this potential in an atomistic simulation to study the mechanical behaviour of edge and screw dislocations in STO. In this work, the Thomas potential [34] is employed to derive the internal force density in the CAC model, which serves as a role for the constitutive relation in continuum mechanics. Both the short-range pair interactions represented by the Born–Mayer model and the long-range coulombic interactions are included in this potential.

Note that, because of the involvement of the long-range coulombic interaction, simulation of ionic materials is always more expensive than that of metals. In order to provide a computationally efficient method for calculating coulombic interactions in non-periodic models, Fennell & Gezelter [39] developed a shifted-force method to emulate the long-range ordering by including the distance-dependent damping function and thus mirror the effective charge screening, based on a damped and cut-off neutralized coulombic sum proposed by Wolf *et al*. [40]. This method can drastically reduce the range of direct electrostatic interactions without losing much accuracy. In this work, the shifted-force method is implemented to calculate the coulombic interaction and a small cut-off of 15 Å is adopted in the simulations using this method. The same cut-off of 15 Å is also used for the calculation of other short-range non-coulombic interactions.

STO is a highly ionic crystal. It has a perovskite crystal structure with five atoms (one chemical formula) in the primitive unit cell and a lattice constant of 3.905 Å. Within each primitive unit cell of STO, there are three different kinds of ions, i.e. Sr, Ti and O ions, with each having different charges, as shown in figure 1*a*. In figure 1*b*, we present the cubic eight-node linear solid element employed in this work to mimic the shape of the cubic primitive unit cell of the STO perovskite structure. Embedded within each node of the element, there is a primitive unit cell of STO that has five atoms and thus 15 DOF. Note that in the CAC method connectivity between finite elements is not required. This then allows us to simulate spontaneous and arbitrary crack initiation and propagation through separation of or sliding between the finite elements without the need for special numerical treatment.

The framework with the plane strain condition used by Tadmor *et al*. [41] and Shenoy *et al*. [42] is assumed in all the CAC and MD simulations, i.e. the displacements in the out-of-plane direction are constrained even though all the atomistic calculations are conducted in three dimensions.

### (b) The equilibrium structure of the grain boundaries

GBs are planar defects that connect adjoining regions with various crystallographic orientations. The ratio between the number of lattice points in the unit cell of a *coincident site lattice* (CSL) and the number of lattice points in a unit cell of the generating lattice is called Sigma (

To simulate the GB–defect interaction by CAC, it is necessary for the method to reproduce the equilibrium or stable GB structure as a result of the interaction of atoms in the material. To obtain a stable GB structure, we first construct an initial CAC bi-crystal model of STO. The model shown in figure 2*a* contains the *grain A* by 26.6° counterclockwise and *grain B* by 26.6° clockwise around the [0 0 1] axis. The GB region is meshed with atomic resolution, while regions away from the GB are meshed with low resolution with each finite element containing 512 unit cells. The equilibrium structures of the bi-crystals are then obtained through CAC simulation of the models under various translation states but without an applied external force field. Thereafter, the GB energy is calculated for each specified microscopic translation state using the following formula [44]:
*E*^{CSL} and *E*^{Bulk} are the total energies of the calculation region with and without CSL, respectively; *A* is the GB area; *n*_{Sr}, *n*_{Ti} and *n*_{O} are the numbers of Sr, Ti and O atoms contained in the calculation region; μ_{Sr}, μ_{Ti} and μ_{O} are the energies of each Sr, Ti and O atom in a perfect crystal.

The variation of *a*, in which the microscopic translation in the *y*-direction, i.e. the direction perpendicular to the GB plane, is 0.5 Å. The *x*-axis in figure 3*a* represents the translation in the *x*-direction, i.e. the direction along the GB plane, and different symbols denote different translation in the *z*-direction, i.e. the out-of-plane direction. For the purpose of comparison, the GB energy variation with different translation states obtained using DFT calculations by Lee *et al*. [20] is plotted in figure 3*b*. It is seen that the DFT result compares well with the CAC result qualitatively, but the overall magnitude of GB energy from CAC simulation is approximately 20% larger than that from DFT calculation. Note that Benedeck *et al*. [35] have calculated the energy of the

Figure 4*a* plots the calculated GB energy as a function of the rigid-body translation for the *y*-direction being 0.5 Å , while figure 4*b* plots the *y*-direction being 1.0 Å. Again, we see that the variation trends obtained from the CAC for the cases of *et al*. [19] and Lee *et al*. [20], respectively.

Figure 5*a* presents the stable GB structure of *dx* *dy* *dz*) being (0.0 0.5 1.5), in which the unit structures around the GB plane are outlined with solid lines. It is noted that the initial translation state is not identical to the final state, indicating that a change in the state after equilibrium has occurred. The equilibrium GB structures in the STO tilt *et al*. [20], in which the main objective was to find the stable GB atomic structure. Figure 5*b* shows the DFT calculation result and figure 5*c* shows the HAADF-STEM image of the GB. A comparison between CAC and DFT results of the outlined unit structure around the GB plane shows that these results from CAC are in excellent agreement with those obtained by DFT, including both the shape of the GB structural units and the arrangement of ions within the units. Moreover, CAC simulation has also reproduced the concomitant misalignment of the ions within the (Ti, O) columns, i.e. the O ions are further inwards (towards the boundary plane) than the Ti ions after relaxation, in exact correspondence with the DFT calculation [20]. The GB unit structure composed of Sr and Ti ions obtained by CAC also compares well with that of the HAADF-STEM experimental image obtained by Lee *et al*. [20], as shown in figure 5*c*.

Figure 6*a* plots CAC simulation results of the *b* the *et al*. [19] and Lee *et al*. [20]. Thus, CAC has reproduced the atomic-level details of all four GB structures, in good agreement with DFT calculations and experimental observations.

Figure 7 compares the GB energies obtained from CAC simulations with those from MD simulations. We have built 10 more CAC models with different sizes of the atomically resolved region by only changing the smallest distance between the atomic–continuum interface and the GB, which is denoted as *d* in figure 2*b* and figure 7*b*. The corresponding MD models being compared with these CAC models are the same. Thus, the GB energies from the corresponding MD results are the same, as shown by the dashed line in figure 7*a*. It is noted that the GB energies from CAC simulations are larger than those from MD simulations when the distance *d* is smaller than 12 Å. However, the CAC results converge to the MD results when *d* is larger than 12 Å. We have also compared the entire displacement field obtained by each CAC model with the atomic displacements obtained by MD. The discrepancy between the CAC and MD results is plotted in figure 8. In figure 8*d*, the distance *d* is 12.22 Å and the discrepancy in the magnitude of the displacement field for the majority of the CAC models is less than 0.1 Å. Based on the results of displacement discrepancies in figure 8 and the GB energy curve in figure 7, we conclude that CAC can reproduce MD results when the distance *d* is larger than 12 Å.

### (c) The effect of the crack–grain boundary interaction

#### (i) A comparison of concurrent atomistic–continuum and molecular dynamics simulation results

STO is an extremely brittle crystalline material with a fracture toughness of 1.0 MPa m^{1/2}. Cleavage is the common form of brittle fracture in crystalline materials and is defined as a split along definite crystallographic structural planes without plastic flow. There are various criteria for the cleavage planes in crystals, such as crystal growth planes, most closely packed crystal planes, and the plane with the smallest bond density. For ionic materials, Schultz *et al*. [45] proposed the crystal plane neutrality condition as a criterion for cleavage plane. It suggests that only the neutral planes can experience cleavage in highly ionic crystal structures. According to this neutrality condition criterion, the {1 0 0} planes that bound the primitive unit cells in STO are the cleavage planes because either of the possible terminating planes SrO or TiO_{2} is electrically neutral.

To allow cracks to initiate and propagate at inter-element boundaries, the eight-node linear cubic solid element is employed in this work to mimic the shape of the cubic primitive unit cell of the STO perovskite structure with all the element surfaces being the cleavage planes, i.e. {1 0 0} planes. A notched bi-crystal with two neighbouring grains that have different crystal orientations is schematically illustrated in figure 9*a*. The [1 0 0], [0 1 0] and [0 0 1] directions of *grain A* are parallel with the *X*, *Y* and *Z* directions in the global coordinate system, respectively. The bi-crystal model contains a *b*, in which the finest mesh, i.e. atomic resolution, is used in regions near the predefined crack and the GB, while a uniform coarse mesh with each element containing 512 unit cells is used in the regions away from the crack and the GB. The number of DOF in the coarse-scale FE region is 1.6% of that of the corresponding atomistically resolved model.

A constant velocity of 5 m s^{−1} is applied to the CAC model. For the purpose of comparison, an all-atom MD simulation of the same specimen with same GB structure is also performed. The CAC and MD simulation results are presented and compared in figure 10. Note that, in the CAC result, the solid black lines in figure 10*a*,*c*,*e* represent the edges of the coarse elements. The positions of atoms within the coarsely meshed elements are calculated through interpolation of obtained FE nodal displacements. It can be seen from figure 10*a* that, when the applied strain in the CAC simulation reaches 2.87%, the crack has propagated into the coarse-meshed FE region, demonstrating that the crack can propagate from the atomistically resolved region to the coarse-meshed FE region without additional constitutive roles or numerical technique. From the crack surfaces, we see that the crack propagates in an unzipping pattern, leaving square-shaped steps composed of the {1 0 0} cleavage planes. A similar phenomenon has been reported by Matsunaga *et al*. [46], who observed in both *in situ* TEM experiments and MD simulations that the atomic-level crack surfaces of MgO are not flat but contain a number of square-shaped steps with a few {0 1 0} atomic layers in height. It is noted that this unzipping pattern disappears in the coarse-meshed FE region and this is due to the effect of the linear shape function used within the coarse elements. Figure 10*b* presents the MD result at the same applied strain, and the same unzipping pattern is also observed. When the applied strain reaches 4.14%, we see from figure 10*c* that the crack has been transmitted from the coarse-meshed FE region into the atomistically resolved region and the unzipping pattern of crack propagation appears again. When the applied strain increases to 4.52%, the interaction between the crack and the GB has resulted in the formation of voids between the two neighbouring grains, as shown in figure 10*e*; thereafter, the voids along the GB make the two grains separate from each other. Figure 10*f* confirms that the same process and interaction mechanism are reproduced in the MD simulation.

The critical stress intensity factor *K*_{IC} was calculated in both CAC and MD simulations.

*K*_{IC} in the CAC simulation is *K*_{IC} in the the MD simulation is *K*_{IC} is 6.846%. The calculated averaged crack speed in the left grain is 447.17 m s^{−1} in the CAC simulation, and 478.80 m s^{−1} in the MD simulation. The discrepancy in the averaged crack speed is 6.606%. These discrepancies are caused by the use of coarse-graining in continuum FE regions.

Figure 11 compares the CAC and MD results of the averaged stress in the tensile direction. It can be seen that the two stress–strain curves agree well qualitatively and quantitatively, and the maximum stresses predicted by CAC and the MD are very close.

The detailed comparisons of the CAC and the corresponding MD results presented in figures 10 and 11 show that, with a large amount of DOF being reduced in the coarse-scale FE region, CAC has reproduced not only the local atomic arrangements in the crack surfaces and the GB during the deformation processes but also the overall stress–strain response of the specimen, in good agreement with MD simulations. Moreover, the MD results of the dynamic processes of the crack initiation and propagation, the interaction between the crack and the GB and the essential fracture behaviour of the bi-crystal STO have been reproduced by CAC.

The computer model shown in figure 9 is a special case in which the loading direction and the crack direction are parallel to the 〈1 0 0〉 crystallographic directions. In order to test the feasibility of the CAC method in simulating the crack propagation in general situations, the computer model shown in figure 12 is constructed such that the loading direction and the crack direction are not parallel to the [010] and the [100] crystallographic directions, respectively, while the same GB and loading conditions as those in figure 9 are used. In figure 12*a*,*b*, we present the CAC and MD models before applying tensile loading, respectively. The deformed CAC and MD configurations are plotted in figure 12*c*,*d*, respectively. It can be seen from figure 12 that the CAC method is still capable of predicting the crack path that is comparable with the MD result in situations where the initial crack and loading directions are not aligned with the crystallographic directions.

#### (ii) Interaction of crack with different tilt grain boundaries

To identify the mechanisms and the effect of the interaction between cracks and GBs, four different submicrometre-sized bi-crystal models are constructed for CAC simulation. The different orientation angles in the models are obtained by rotating grain B in figure 9*a* while keeping grain A fixed. Specific translation states are then applied between the two grains in the models. Details of those models are listed in table 1. The first two models contain the

Figure 13 plots the averaged stress of the four bi-crystal models under tensile loading obtained by CAC simulations. It can be seen that the stress–strain curves for the four different models overlap at the initial loading stage until the tensile stress is about to reach the maximum value. When the strain comes to be about 0.042, we can clearly observe a change in the stress–strain curve's slope for each model. This change of the curve slope is caused by the interaction between the propagating crack and the GB. Additionally, the slope change is different among the four models due to different interaction mechanisms. Thereafter, when the applied strain reaches 0.0525, the crack almost propagates across the whole specimen. We observe that the corresponding stresses, i.e. the fracture strength, decrease in the order of

Figure 14 presents snapshots of the atomic arrangements obtained by CAC simulations of the four models at the loading strain of 4.52%. It is seen from figure 14*a* that the propagating crack in the *b*, we observe that the propagating crack in the

From figure 14*c*, an intra-granular fracture in the adjacent grain can be observed in the *ad hoc* interfaces or additional criteria. Figure 14*d* shows the similar intra-granular fracture pattern in the

The simulation results of the GBs with four different misorientation angles listed in table 1 show that the maximum averaged stress decreases with the decreasing misorientation angle, as indicated in the stress–strain curves in figure 13. This phenomenon can be explained through the crack–GB interactions and the deformation mechanisms. For the two

## 4. Summary and discussions

This paper has presented a new development in the CAC methodology for the study of interactions between defects in ionic crystalline materials and has demonstrated the method through simulations of the structures of four [0 0 1] tilt GBs and the dynamic processes of crack–GB interactions in STO bi-crystals. The CAC bi-crystal models are meshed with atomic resolution in the GB region while with coarse-scale finite elements in regions away from the GB. It has been demonstrated that the CAC allows a smooth passage of cracks through the atomistic–continuum interface without the need for special numerical treatment or additional constitutive rules for dictating crack initiation or propagation. Although 98.4% of the DOF of the simulated atomistic system have been eliminated by coarse-graining in the far regions, the CAC method has reproduced both structural and dynamic properties of the GBs as well as the overall material behaviour of the bi-crystal models.

The atomic-scale GB structures and the GB energies of four STO [0 0 1] tilt GBs have been reproduced for the first time by a CAC method, in good agreement with those obtained from existing experiments and DFT calculations. The ability of the CAC to reproduce the accurate GB structure, i.e. the relative positions of atoms in the GB, of the ionic material STO originates from the unified treatment of the atomistic and continuum regions in the CAC methodology. We would like to emphasize that it is necessary for a method to first obtain the accurate local atomic structure of GBs for the study of the GB–defect interaction.

In the simulation of the dynamic processes of the crack–GB interaction, the CAC results, including stress–strain responses, the dynamic processes of crack initiation and interactions with GB and the mechanisms of the fracture, are all comparable with the corresponding MD results. In addition, CAC simulation results demonstrate that the criteria for cleavage planes in ionic material are different from those in metals. The unzipping and zigzag fracture patterns with the crack surface being {1 0 0} neutral planes found in the simulations confirm that the neutrality condition criterion for cleavage planes holds true for STO crystalline materials and the established rules and deformation mechanisms in metals no longer apply to ionic materials because of the coulombic interactions between ions.

The CAC simulations have predicted that the fracture strength of different GBs decreases in the order of

Note that the CAC method differs from most of the existing multi-scale methods in that it also solves the internal deformation embedded within each material point in the coarsely meshed FE domain. We have previously explained the necessity of including the internal DOF of atoms from the three sublattices in STO in the method in order to describe the proper coulombic interaction between ions of different charges. The simulation results presented in this work have demonstrated a significant effect of the atomic-scale arrangement in the GB region on the crack–GB interaction and consequently on the deformation mechanisms and fracture strength of the models. This further confirms the need to include the internal DOF of atoms in the method for simulation of STO. In addition, different from that in metals, the GB structures in STO are open structures, and the density of the GB region is hence significantly lower than that in bulk regions. This kind of open structure is common in ionic materials such as MgO [18] and *Mg*_{2}*SiO*_{4} [44] due to the existence of electrostatic interaction between ions with different charges. The CAC simulations have reproduced these open GB structures and revealed their effect on the mechanical behaviour. The ability to capture such atomic-scale behaviour by CAC is due to the atomistic information (the arrangement and the interaction of atoms) built in the formulation. It is thus reasonable to believe that the CAC method can be used for simulation and prediction of complex material behaviour of polyatomic ionic crystals with GBs and other defects.

## Authors contributions

S.Y. designed and carried out the simulations, conducted the data analysis and wrote the manuscript with the assistance and supervision of Y.C.

## Funding statement

This material is based upon the work supported by the Department of Energy under award no. DOE/DE-SC0006539. The CAC computer code in its present form is a culmination of developments supported in part by the National Science Foundation under award nos. CMMI-1233113 and CMMI-1129976. The MD simulation work was also supported in part by the National Science Foundation through TeraGrid resources provided by Kraken under grant no. DMR 110065. The CAC simulations were performed at the High Performance Computing Center at the University of Florida.

## Conflict of interests

We have no competing interests.

## Acknowledgements

The authors would like to thank three anonymous reviewers for their constructive comments. S.Y. also would like to thank Prof. Liming Xiong and Ms Xiang Chen for many discussions on this work.

- Received October 5, 2014.
- Accepted January 9, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.