Quantifying the High-Temperature Separation Behavior of Lamellar Interfaces in γ-Titanium Aluminide Under Tensile Loading by Molecular Dynamics

γ-titanium aluminide (TiAl) alloys with fully lamellar microstructure possess excellent properties for high-temperature applications. Such fully lamellar microstructure has interfaces at different length scales. The separation behavior of the lamellae at these interfaces is crucial for the mechanical properties of the whole material. Unfortunately, quantifying it by experiments is difficult. Therefore, we use molecular dynamics (MD) simulations to this end. Specifically, we study the high-temperature separation behavior under tensile loading of the four different kinds of lamellar interfaces appearing in TiAl, namely, the $\gamma /{\alpha }_2$, $\gamma /{\gamma }_{\text{PT}}$, $\gamma /{\gamma }_{\mathrm{TT}}$, and $\gamma /{\gamma }_{\mathrm{RB}}$ interfaces. In our simulations, we use two different atomistic interface models, a defect-free (Type-1) model and a model with preexisting voids (Type-2). Clearly, the latter is more physical but studying the former also helps to understand the role of defects. Our simulation results show that among the four interfaces studied, the $\gamma /{\alpha }_2$ interface possesses the highest yield strength, followed by the $\gamma /{\gamma }_{\mathrm{PT}}$, $\gamma /{\gamma }_{\mathrm{TT}}$, and $\gamma /{\gamma }_{\mathrm{RB}}$ interfaces. For Type-1 models, our simulations reveal failure at the interface for all γ/γ interfaces but not for the $\gamma /{\alpha }_2$ interface. By contrast, for Type-2 models, we observe for all the four interfaces failure at the interface. Our atomistic simulations provide important data to define the parameters of traction–separation laws and cohesive zone models, which can be used in the framework of continuum mechanical modeling of TiAl. Temperature-dependent model parameters were identified, and the complete traction–separation behavior was established, in which interface elasticity, interface plasticity, and interface damage could be distinguished. By carefully eliminating the contribution of bulk deformation from the interface behavior, we were able to quantify the contribution of interface plasticity and interface damage, which can also be related to the dislocation evolution and void nucleation in the atomistic simulations.

1 Introduction

Rising concern for environmental protection demands the quest for lightweight high-temperature materials to improve fuel efficiency in civil aviation. In recent times, titanium aluminide (TiAl) alloys (Dimiduk, 1999) are evolving as a potential material for rotating components operating at elevated temperatures, especially turbine blades. Compared to other high-temperature materials, alloys based on γ-TiAl (Appel et al., 2000) possess a high specific modulus and melting point, a low density, and excellent corrosion resistance. In particular, the specific strength of γ-TiAl is higher than that of superalloys at elevated temperatures, making it an attractive engineering material (Clemens and Mayer, 2013). γ-TiAl can exhibit a duplex, equiaxed, nearly lamellar (NL), or fully lamellar (FL) microstructure. Among them, the FL microstructure exhibits superior creep resistance and fracture toughness (Bewlay et al., 2016). It is characterized by thin lamellae, each of which belongs to one of the two constituent intermetallic phases, the γ phase (ordered face-centered tetragonal TiAl) and the ${\alpha }_2$ phase (ordered hexagonal close-packed Ti₃Al). In the FL microstructure, a lamella generally consists of multiple domains within which an ordered crystallographic structure is preserved. On the other hand, multiple lamellae form the so-called (grain-shaped) colonies, leading to interfaces at different length scales, namely, colony level, domain level, and lamellar level. Two types of interfaces in FL TiAl substantially affect its high-temperature deformation properties: the colony boundaries on the larger scale and the lamellar interfaces on the lowest scale, which are either of the type ($\gamma /\gamma$ or $\gamma /{\alpha }_2$).

Continuum mechanical modeling of the micromechanics of materials has attracted considerable attention over the last decades. Appel et al. (2016) recently reviewed the different modeling approaches to TiAl proposed so far. Crystal plasticity models can help to account for specific deformation mechanisms (slip and twinning) and for describing the anisotropic material behavior in TiAl. Recently, such models (Schnabel and Bargmann, 2017; Ji et al., 2018; Schnabel and Scheider, 2020) were used to predict the colony boundary strengthening coefficient as a function of lamella thickness and to map the deformation in lamellar TiAl. Microstructure-informed multi-scale models bear promise to accelerate the alloy development in γ-TiAl. To increase their accuracy, they should also include the role of interfaces and alloying elements. Obtaining the parameters required for continuum mechanical interface models from experiments is difficult so that molecular dynamics (MD) simulations appear to be the most promising tool to this end.

In MD simulations of TiAl, the different intermetallic phases ($\gamma$ and ${\alpha }_2$) are modeled with their respective crystallographic systems (i.e., γ: face-centered tetragonal and ${\alpha }_2$: hexagonal close-packed) and appropriate orientations as well as with suitable interatomic interaction potentials (Zope and Mishin, 2003; Kim et al., 2016). In recent time, several MD studies have been devoted to deepen our understanding of the deformation behavior of γ-TiAl (Zhou et al., 2004; Xie et al., 2015; Kanani et al., 2016; Wu et al., 2016; Li et al., 2017; Feng et al., 2018; Hui et al., 2018; Cao et al., 2019; Ding et al., 2019; Feng et al., 2019; Li et al., 2019; Li et al., 2020) in combination with experiments. Kanani et al. (2016) performed MD shear simulations of a distinct $\gamma /\gamma$ interface, observing different deformation mechanisms and strong in-plane anisotropy of shear strength. Wu et al. (2016) studied the influence of surface defects on the mechanical behavior of the γ phase under tensile deformation. They (Wu et al., 2016) observed for the defect-free bulk sample brittle fracture, but plastic deformation before fracture for samples with surface defects normal to the loading direction. In an extensive study using first principles, Ji et al. (2018) obtained a deformation map for plastic deformation modes using the orientation-dependent effective energy barrier. Li et al. (2019) studied the homogeneous tensile deformation behavior of the $\gamma /{\alpha }_2$ interface at 0.01 K and found an influence of boundary condition on the observed failure mode with plastic deformation for free boundaries and fracture for constrained boundaries. These observations can also be interpreted by the Schmid factor analysis of the corresponding slip system. A recent MD study (Ding et al., 2019) has been devoted to the effect of grain size and temperature, however, for polycrystalline TiAl. A decreasing Young’s modulus and average flow stress was observed with increasing temperature. Together, these studies can help to understand experimental findings about room temperature ductility. However, more atomistic insights are still necessary to understand the thermomechanical deformation behavior of a single lamellar interface in γ-TiAl at elevated service temperatures. Our understanding of the combined microscopic effects and mechanisms occurring in plasticity and fracture of TiAl still remains limited (Appel et al., 2016), even though some small-scale testing has been conducted revealing the failure mechanisms in the microstructure (Werwer et al., 2007; Barbi et al., 2012).

The finite element method is frequently used in mesoscale and macroscale models of TiAl. Such models should include the deformation behavior of bulk constituent phases as well as different interfaces ($\gamma /{\alpha }_2$ and $\gamma /\gamma$). In γ-TiAl, the anisotropy arising from the complex microstructure can be captured by crystal plasticity models (Schnabel and Bargmann, 2017; Schnabel and Scheider, 2020). For modeling the interfaces between colonies, domains, and lamellae, cohesive zone models are promising tools. These phenomenological models initially developed to handle discontinuities due to crack propagation and fracture (Hillerborg et al., 1976; Li and Chandra, 2003; Scheider, 2018) can be embedded in finite element models by the so-called interface elements to also describe the material separation behavior at both macroscale and micro-scale interfaces (Scheider, 2009a; Simonovski and Cizelj, 2015; Scheider, 2018). For lamellar γ-TiAl, micromechanical lamellae separation has been studied in Werwer and Cornec (2000), Werwer et al. (2007), and Wei et al. (2009). Cohesive zone models rely on a traction–separation law (TS law). On the macroscale, its parameters can be identified by appropriate mechanical tests (Cornec et al., 2003). However, debonding processes are difficult to observe with in-situ experiments on the nanoscale of atomistic lattice boundaries. Classical MD simulations are typically more suitable on that scale. The identification of TS laws from atomistic simulations instead of experimental tests has started in the early 2000s and has remained an active field of research since then (Yamakov et al., 2006; Zhou et al., 2008; Zhou et al., 2009; Krull and Yuan, 2011; Möller et al., 2018). For high-temperature alloys, such as γ-TiAl, in particular, temperature dependency has to be included in such studies. An important aspect when the TS law is derived from micromechanical or atomistic simulations is that usually the damage in a finite volume (called the damage zone) is taken as a basis for the separation behavior. The elimination of the “bulk” response from the damage zone (i.e., solely considering the interface damage) has been addressed here not only based on a method already used by the authors (see (Scheider, 2009b, Scheider, 2009a)) but also in a similar fashion employed by others (see, e.g., (Paggi and Wriggers, 2011)).

In this article, we perform MD simulations of four single lamellar interfaces (one $\gamma /{\alpha }_2$ and three $\gamma /\gamma$ in different orientations) in γ-TiAl under tensile load from room temperature to elevated service temperatures. This way, we first obtain general atomistic insights into the separation behavior of single lamellar interfaces. Second, we obtain atomistically informed traction–separation (TS) curves for each material interface, depending on temperature, strain rate, and preexisting defects. These are used to formulate an atomistically informed cohesive zone model, which can be used in future finite element mesoscale simulations of TiAl.

2 Materials and Methods

2.1 Molecular Dynamics Model

At the micro- or nanoscale, materials can be modeled as many-particle systems. MD is a numerical technique to study the evolution of such systems in time by computing the trajectories of the individual particles. This is accomplished by numerical integration of Newton’s equation of motion. If the particles represent atoms, MD simulations provide an atomic resolution (Cai et al., 2012). Herein, all MD simulations were performed using the open source package LAMMPS large-scale atomic/molecular massively parallel simulator (Plimpton, 1995) distributed by the Sandia National Laboratories. For metallic materials, interatomic interactions in MD are often modeled in a (semi)empirical manner using the so-called embedded-atom method (EAM) of Daw and Baskes (1984). Following this approach, we define the total potential energy $\left(U_{Total}\right)$ of the many-particle (N atoms) system as

$U_{Total}=\underset{i=1}{{\displaystyle \sum }^N}\left[\frac{1}{2}{\displaystyle \sum }_{j\ne i}{\varphi }_{ij}\left(r_{ij}\right)+F_i\left({\displaystyle \sum }_{j\ne i}{\rho }_j\left(r_{ij}\right)\right)\right].(1)$

Here, $r_{ij}$ and ${\varphi }_{ij}$ are the distance and interaction potential between the ith and jth atom. N indicates the total count of atoms in the system. $F_i$ denotes the embedding energy of the ith atom due to the cumulated host electron density ${\displaystyle {\sum }_{j\ne i}{\rho }_j\left(r_{ij}\right)}$ provided by the surrounding atoms. The pair interaction potential and electron density are summed over the neighboring j atoms defined by the cutoff. For γ-TiAl, we used the EAM potential from Zope and Mishin (2003), which was validated by a comparison of experiments and ab initio simulation data and found to yield, for example, accurate lattice constants, cohesive energy, and elastic constants. In particular, its ability to predict the lattice thermal expansion makes this potential suitable for MD simulations at high temperatures. Therefore, it has also been used by Xie et al. (2015), Kiselev and Zhirov (2014) for studying the synthesis of intermetallic γ-TiAl from Al and Ti plates at 1400 K and the solidification of ${\alpha }_2$-Ti₃Al from the liquid phase at 2000 K. Remarkably, this potential correctly predicts the ${\alpha }_2$-Ti₃Al equilibrium structure although the relevant information was not included in the fitting database. For these reasons, we chose this EAM potential of Zope and Mishin (2003), whose developers underline its suitability for large-scale atomistic simulations of plastic deformation and fracture in γ-TiAl even at high temperatures.

2.2 Cohesive Zone Model

In a cohesive zone model, interfaces are idealized as two surfaces initially attached to each other without any material volume in between. Tractions on the interface are defined by a vector field $\mathbf{t}$ and may cause a relative displacement of the two surfaces forming the interface, which can be modeled by the so-called separation vector $〚\mathbf{u}〛={\mathbf{u}}^{+}-{\mathbf{u}}^{-}$ with ${\mathbf{u}}^{+},{\mathbf{u}}^{-}$ being the displacements of the two surfaces forming the interface, respectively. The relation between traction vector and separation vector can be formulated by the so-called TS law $\mathbf{t}=\mathbf{t}\left(〚\mathbf{u}〛\right)$, which can be understood as the mechanical constitutive laws of the interface itself. A host of different ways to define such laws has been proposed (see, e.g., the reviews in Park and Paulino (2011), Scheider (2018)). Here, a recently developed model is used, which splits the separation additively into a reversible (elastic) and an irreversible (plastic) part

$〚\mathbf{u}〛={〚\mathbf{u}〛}^{el}+{〚\mathbf{u}〛}^{pl}.(2)$

Similar to bulk plasticity models, the interface traction depends only on the elastic part of the separation ${〚\mathbf{u}〛}^{el}$. A large variety of TS laws is in principle possible. Herein, we focus on the mechanical behavior of the interlamellar interfaces in TiAl normal to the interfaces under normal loading only. Therefore, we only consider the normal components of the traction and the separation vector. The normal components of the full, elastic, and plastic separation vector are denoted by $〚u_n〛,〚u_n^{el}〛,\text{ and }〚u_n^{pl}〛$. Analogously, the normal component of the traction vector is denoted by $t_n$. The simplest relation of separation and traction is a linear one, which is also very common for cohesive zone models since the elastic separation is usually small. Since the traction is the derivative of the elastic potential (per unit area) W, we start with the following second order polynomial

$W=\left(1-d\right)\underset{W_0}{\underbrace{\frac{1}{2}{C}_n{〚u_n^{el}〛}^2}}.(3)$

Here, $W_0$ describes the elastic energy in the absence of any damage, and the damage parameter d ranges from $d=0$ for the damage-free state to $d=1$ at failure. The material parameter $C_n$ defines the interface stiffness. This relation can easily be generalized to three dimensions. The traction vector follows from the above elastic energy as

$t_n=\frac{\partial W}{\partial 〚u_n^{el}〛}=\left(1-d\right)C_n〚u_n^{el}〛.(4)$

Following the classical setting of damage and plasticity theory (Maugin, 1992), we define the so-called fictitious undamaged stress ${\tilde{t}}_n=t_n/\left(1-d\right)$. The plastic regime of our traction separation law is assumed to be defined by a yield function

$g\left(t_n,d\right)={\tilde{t}}_n-t_0+D\left(d\right)=0.(5)$

Here, $t_0$ denotes an initial yield stress and $D\left(d\right)$ is reduction due to damage. The results of our MD simulations indicated that plastic hardening of the interface is negligible, and it is thus not represented in Eq. 5. This also implies that the Helmholtz free energy (per unit area) of the interface is W from Eq. 3, which is thus used in Eq. 9 to compute the thermodynamic driving force of the single internal variable d of our model. The flow rule associated to the yield function g provides the plastic flow rate normal to the interface

$〚{\dot{u}}_n^{pl}〛=\dot{\lambda }\frac{\partial g}{\partial t_n}=\dot{\lambda }\frac{\partial {\tilde{t}}_n}{\partial t_n}=\frac{\dot{\lambda }}{1-d},(6)$

with $\dot{\lambda }$ being the plastic multiplier. The accumulated plastic separation α can be defined as

$\alpha \left(t\right):={{{\displaystyle \int }}^{\text{}}}_0^t\dot{\lambda }\left(\tau \right)/\left[1-d\left(\tau \right)\right]\text{d}\tau .(7)$

The damage-related reduction of yield stress D in Eq. 5 is defined herein heuristically as

$D\left(d\right)=\frac{1}{{[1-d]}^p}\frac{S}{2}{\left[\frac{W_0}{S}\right]}^2,(8)$

with two material parameters p and S. The evolution of the internal variable d can be defined in a thermodynamically consistent manner (see, e.g., Maugin, 1992) via its thermodynamic driving force

$Y=\frac{\partial W}{\partial d}=-W_0,(9)$

which by standard thermodynamic arguments from classical plasticity theory directly leads to the evolution equation

$\dot{d}=-\dot{\lambda }\frac{\partial g}{\partial Y}=\frac{\dot{\lambda }}{{[1-d]}^p}\left[\frac{Y}{S}\right],\text{if}\alpha \ge {\alpha }_{nucl}(10)$

Here, the initial value of d is zero, and ${\alpha }_{nucl}$ represents an additional model parameter below which no damage at all is assumed to occur. Moreover, it can be shown that

$d=1-{\left[1-\frac{\alpha -{\alpha }_{\text{nucl}}}{{\alpha }_{\text{fail}}-{\alpha }_{\text{nucl}}}\right]}^{\frac{1}{p+1}},(11)$

where ${\alpha }_{\text{fail}}$ denotes the plastic separation at failure, that is, at $d=1$. Note that with Eq. 11, ${\alpha }_{\text{fail}}$ can be used to replace effectively the model parameter S. Thus, our cohesive zone model can be written using only the following five material parameters: $C_n$ for characterizing the elastic properties of the interface and $t_0,p,{\alpha }_{\text{fail}},\text{ and }{\alpha }_{\text{nucl}}$ for characterizing plasticity and damage.

While the cohesive zone model has been implemented into the commercial finite element code ABAQUS^®, a python implementation of the above equations is sufficient for the calculations presented here. The python class for the cohesive zone model is available at Scheider (2020).

2.3 Molecular Dynamics Simulation—Preprocessing

In γ-TiAl processing, the target microstructure is often achieved through the controlled cooling of the high-temperature α phase, which results in the precipitates of γ-plates through the reactions $\alpha \to \alpha /\gamma \to {\alpha }_2/\gamma$ and $\alpha \to {\alpha }_2\to {\alpha }_2/\gamma$ (Kim and Dimiduk, 1991). Such nucleation of the γ phase from the parent α or ${\alpha }_2$ phase leads to four different lamellar interfaces. In addition to the $\gamma /{\alpha }_2$ interface, there exist three different γ/γ interfaces in the FL microstructure (Jin and Gray, 1997; Appel et al., 2011; Kanani et al., 2016), which represent different rotational variants. Herein, we studied all the four different interfaces. To this end, we constructed atomistic simulation models (refer Figure 1) whose two halves were a combination of intermetallic constituent phases representing two different lamallae. The interface between both phases represented the respective interlamellar interface. Based on the general setting illustrated in Figure 1, we designed the simulation domains in the following way. The face-centered tetragonal crystallographic structure with a₀ = 3.998 Å and c/a₀ = 1.047 (Zope and Mishin, 2003) of the γ-TiAl phase was represented by 96 atomic layers along the X-direction [1–1 0] and Y-direction [1 1–2], respectively, and by 48 atomic layers in the Z-direction [1 1 1]. Likewise, the hexagonal close-packed crystallographic structure with a₀ = 5.7884 Å and c/a₀ = 0.820 of the ${\alpha }_2$-Ti₃Al phase was represented by 96 atomic layers in the X-direction [1–1 0], 64 atomic layers in the Y-direction [1 1 0], and 48 atomic layers in the Z-direction [0 0 1]. Here, the respective number of atomic layers was chosen to ensure periodicity in all directions. The $\gamma /{\alpha }_2$ interface was aligned with the close-packed plane and direction, that is, ${\left\{111\right\}}_{\gamma }\Vert \left(0001\right){\mathit{\alpha }}_{\mathrm{2}}$ (Appel et al., 2011). The three rotational variants of the γ phase were achieved by a rotation around ${\langle 111\rangle }_{\gamma }$ by 60°, 120°, and 180°, respectively, resulting in a pseudo twin (PT), rotational boundary (RB), and true twin (TT) configuration, denoted as $\gamma /{\gamma }_{\mathrm{PT}}$, $\gamma /{\gamma }_{\mathrm{RB}}$, and $\gamma /{\gamma }_{\mathrm{TT}}$ interfaces, respectively.

FIGURE 1

FIGURE 1. Schematic showing the two types of single interface models used herein for atomistic simulations. (A) Type-1 models assume an initially defect-free (ideal) structure, whereas (B) Type-2 models include a preexisting penny-shaped void across the interface. Atoms belonging to the void region are removed.

In our study, emphasis was placed on understanding the separation behavior of the single interlamellar interface subject to tensile loading. Accordingly, we consider a defect-free interface model designated as Type-1 (refer Figure 1) to assess the theoretical strength and deformation behavior. Further, we considered the single interface models with a defect, which we refer to as the Type-2 model (refer Figure 1), to study the influence of the preexisting defect at the interface. Here, a line defect, that is, dislocation, was avoided on purpose as we were interested in neither the dislocation dynamics nor the interaction of existing dislocation with chemical species. Moreover, the simplified model of single dislocation poses several limitations: stability of the dislocation at elevated temperature, the bias introduced by the initial dislocation type, and unrealistic starting dislocation density. To circumvent these limitations, we chose a void type defect instead. Further, the preexisting void can be associated with the crack at the interface, which is often considered in continuum fracture studies to understand the crack propagation under different loading conditions. From the atomistic standpoint, such defect may act as dislocation sources during deformation, overcoming any modeling bias, as we discussed before. Accordingly, the Type-2 model includes a penny-shaped void at the interface with a radius of 1.6 nm and a height of 1 nm (refer Figure 1). The center of this void was chosen to coincide with the center of the simulation domain so that the void covers two atomic layers in each of the two intermetallic phases in the simulation model. The influence of the initial void shape and size on the interface separation behavior falls beyond this study’s scope. In test simulations, we found that a simulation domain of the size 13.7 nm × 8 nm × 22 nm was a good trade-off between computational cost and the requirement to reduce boundary effects. Using such domains, Type-1 models consist of 147,456 atoms in total and Type-2 models of 147,002 (after removing atoms belonging to the void region).

In Figure 2, the first row shows the constructed atomistic models of the four interlamellar interfaces in TiAl after an initial energy minimization. The second row displays in this initial configuration for each atom the so-called centrosymmetry parameter (CSP), which measures the local lattice disorder around the atom, for example, due to defects, surfaces, or interfaces (Stukowski, 2009). For an atom whose direct neighborhood forms a perfect lattice, the CSP is zero. By contrast, the CSP takes on high values if the atom constitutes the core of a dislocation or is located next to an interface or surface. The CSP was computed using OVITO (Stukowski, 2009). Here, blue atoms represent the γ phase, green atoms represent either an interface or the ${\alpha }_2$ phase, and red atoms represent the boundary of the simulation domain. Note that at the $\gamma /{\gamma }_{\mathrm{RB}}$ interface, the original stacking order of the γ-phase is preserved and thus no elevated CSP was observed.

FIGURE 2

FIGURE 2. First row: simulation system for the defect-free interfaces ((A)$\gamma /{\alpha }_2$, (B)$\gamma /{\gamma }_{\mathrm{PT}}$, (C)$\gamma /{\gamma }_{\mathrm{RB}}$, and (D)$\gamma /{\gamma }_{\mathrm{TT}}$) after energy minimization. Second row:computed centrosymmetry parameter (CSP) to categorize the bulk constituent phase and interface. Here, atoms with blue color code: γ; green color code: either the interface or ${\alpha }_2$ phase; red color code: surface atoms.

The preprocessing of the simulation samples was performed in several steps. In the first step, Type-1 (defect-free) interface models were energetically minimized using a conjugate gradient algorithm in LAMMPS large-scale atomic/molecular massively parallel simulator (Plimpton, 1995). For Type-2 models, now the atoms in the void regions were removed and energetically minimized. For both Type-1 and Type-2 models, the second step was assigning velocities to all the atoms using a Gaussian distribution to initialize the start temperature. Subsequently, the samples were heated up to the target temperature and equilibrated at constant temperature and pressure for 40 ps in case of Type-1 models and 100 ps in case of Type-2 models. Finally, in the third step, the samples were equilibrated at constant volume and target temperature for 100 ps so that the kinetic energy is well distributed within the model. The Nosé–Hoover thermostat was used for all the thermal equilibrations. A time step size of 1fs and periodic boundary conditions in all directions were applied in the preprocessing simulations.

2.4 Molecular Dynamics Simulation—Deformation

After the preprocessing as described in section 2.3, the deformation according to the applied strain rate was simulated. Three different regions were distinguished for imposing boundary conditions: a loading region, a fixed region, and a mobile region in between. These regions are depicted in Figure 3. Here, we choose a boundary-driven deformation technique to perform the tension test. Accordingly, we selected a few layers of atoms (spanning 2 nm) far away from the interface at the boundary of the simulation box in the Z-direction. The loading region’s thickness is carefully chosen to ensure the smooth transfer of displacement rate into the system without any boundary artifact. In the loading region, we prescribe similar to Zhou et al. (2008), Gupta et al. (2016) a far-field displacement $u_Z^{glo}$ which increases over time. A rate ${\dot{u}}_Z^{glo}=2\text{Å}/\text{ps}$ corresponds to a strain rate of ${\dot{\epsilon }}_{ZZ}^{glo}=1\cdot {10}^{10}{s}^{-1}$ and ${\dot{u}}_Z^{glo}=0.2\text{Å}/\mathrm{ps}$ to ${\dot{\epsilon }}_{ZZ}^{glo}=1\cdot {10}^9{s}^{-1}$ (both cases were studied in separate simulations) in the Z-direction was applied to all atoms in the loading region. It is a well-known fact that MD simulations possess severe time scale limitation and inherently suffer high strain rate in comparison to the experiment. Following an investigation conducted by Wu et al. (2016), who showed that the effect of the strain rate is reduced for rates smaller than ${\dot{\epsilon }}_{ZZ}^{glo}={10}^9{s}^{-1}$ in γ-TiAl, we chose the above mentioned strain rates in this study. In the fixed region, the same procedure was used but the displacement rate was set to zero, such that atoms in this region remained stationary. The positions and velocities of the atoms in the loading region were updated based on the assigned constant displacement rates, whereas for the atoms in the mobile region, they were gained in the usual way by numerical integration in time with a time step size of 2 fs. A constant target temperature and periodic boundary conditions in X- and Y-directions were maintained during the deformation simulations. A nonperiodic boundary condition was applied in the Z-direction. Thus, the loading can be described as isothermal uniaxial straining.

FIGURE 3

FIGURE 3. In the loading region (top/green), a constant velocity is prescribed in the Z-direction for all atoms. In the fixed region (bottom/orange), atoms are kept fixed at their position. In the mobile region (the whole rest of the simulation domain), no specific position or velocity was prescribed, but a local observation region with a height $h^{\text{loc}}$ is defined in order to evaluate the local separation behavior.

Figure 3 highlights a special zone within the mobile region, in the following denoted as the local observation region. It stretches over 3.5 nm in both directions from the interlamellar interface, that is, over in total $h^{\text{loc}}=7$ nm. The role of this region is to observe the local separation behavior; hence it does not affect the simulation but supports the postprocessing of the results.

The local atomic stress tensor ${\widehat{\mathbf{S}}}^i$ (units: stress times volume) of each atom i is calculated by

${\widehat{\mathbf{S}}}^i=-m^i{\mathbf{v}}^i\otimes {\mathbf{v}}^i-\frac{1}{2}{\displaystyle \sum }_{j\in \mathcal{N}}\left({\mathbf{r}}^i-{\mathbf{r}}^j\right)\otimes {\mathbf{F}}^{ij}(12)$

at each discrete time step in the simulation according to Thompson et al. (2009). Here, $\otimes$ denotes a tensor product and $\mathcal{N}$ is the set of atoms in the model. The first term on the right-hand side represents a kinetic, thermal contribution with $m^i$ being the mass of the ith atom and ${\mathbf{v}}^i$ its velocity vector. The second term refers to the pair interactions between the ith atom and its surrounding atoms with ${\mathbf{r}}^i\text{ and }{\mathbf{r}}^j$ the position vectors of the ith and jth atoms, respectively, and ${\mathbf{F}}^{ij}$ is the interaction force vector between the ith and the jth atom. For the transition to the continuum scale, the atomic stresses have to be divided by the current subdomain volume $V_{\text{box}}$, giving

$\mathbf{S}=\frac{1}{V_{\text{box}}}{\displaystyle \sum }_{i\in {\mathcal{N}}_{\text{box}}}{\widehat{\mathbf{S}}}^i,(13)$

where the index i runs over all atoms in the subdomain within the simulation model considered.

For mapping atomic stresses to the continuum level, we used two types of such averaged stresses, a global stress, ${\mathbf{S}}^{\text{glo}}$, and a local stress, ${\mathbf{S}}^{\text{loc}}$. For ${\mathbf{S}}^{\text{glo}}$, we averaged across the whole simulation domain in the deformed configuration (initially, $V_{\text{box}}=22\times 13.7\times 8{\text{nm}}^3$), which provides a quantity similar to the pressure tensor in the MD simulation model. For ${\mathbf{S}}^{\text{loc}}$, only the atomic stresses in the local observation region were averaged (initially, $V_{\text{box}}=7\times 13.7\times 8{\text{nm}}^3$). The $ZZ$ component of this stress tensor, $S_{\text{ZZ}}^{\text{loc}}$, was assumed to correspond to the component of the interface traction vector normal to the interface ($t_n$, as introduced in subsection 2.2) in the cohesive zone model (Dandekar and Shin, 2011). Therefore, whenever $S_{\text{ZZ}}^{\text{loc}}$ is considered, the term (normal) traction will be used. $S_{\text{ZZ}}^{\text{loc}}$ and $S_{\text{ZZ}}^{\text{glo}}$ are expected to be equal since both correspond to the same total loading imposed by the boundary conditions at upper and lower boundaries. Nevertheless, we use these two different quantities in the following to match the variables that we use for discussing the deformation. For the deformation, it is favorable to distinguish between the one in the local observation region and the global one because the latter may, especially in the inelastic regime, also contain some bulk effects which we would like to eliminate from our discussion of interface mechanics.

For relating the atomistic to the continuum level, not only stresses but also strains and the interface separation have to be discussed. The global MD-based strain in the Z-direction ${\epsilon }_{ZZ}^{glo}$ is the change of the total simulation box size in the Z-direction divided by its initial length in the Z-direction. The MD-based (normal) interface separation in the local observation region is the difference between the average displacement in the Z-direction of the atoms in the layer above and below the local observation region (refer Figure 3), that is,

$〚u_n^{loc}〛=u_Z^{lo{c}^{+}}-u_Z^{lo{c}^{-}}.(14)$

Both the MD-based tractions and the interface separation were computed in our simulations every 100 fs.

2.5 Identification of Damage Specific Interface Separation

The separation $〚u_n^{loc}〛$ defined by the difference between $u_Z^{lo{c}^{+}}$ and $u_Z^{lo{c}^{-}}$ according to Eq. 14 depends on the size of the local observation region. The reason is that the deformation in this region is not only governed specifically by interface effects but also governed to some extent by deformation of the bulk material, in particular its elastic strain, which occurs independently of the interface and which should not be modeled as part of the deformation resulting from the interface. To eliminate such bulk effects, we follow the ideas used already in Scheider (2009b, a). That is, we first assume that the interface-specific effects such as damage, which we want to capture by our cohesive zone model, occur only in the neighborhood of the interface, that is, in the local observation region. Outside this region, we assume the deformation to be equal to the bulk deformation that would result from a global loading ${\mathbf{S}}^{\text{glo}}$ also in the absence of the interface. This bulk deformation is thus assumed to correspond to a strain

${\epsilon }_{ZZ}^{\text{bulk}}=\frac{u_Z^{glo}-〚u_n^{loc}〛}{h_0-h^{\text{loc}}}.(15)$

Now, to isolate the interface-related part of the deformation in the local observation region, we have to subtract from the overall deformation $〚u_n^{loc}〛$ of this region the part attributed to the general bulk deformation, giving

$〚u_n〛=〚u_n^{loc}〛-{\epsilon }_{ZZ}^{\text{bulk}}{h}^{\text{loc}}.(16)$

With the interface-related separation $〚u_n〛$ and an associated traction $S_{\text{ZZ}}^{\text{loc}}$, the parameters of the cohesive zone model from section 2.2 can be identified by fitting Eqs 4–10 to the results of the MD simulations.

Moreover, the total dislocation density in the MD simulations was calculated using the dislocation analysis tool provided in the visualization software OVITO (Stukowski, 2009; Stukowski and Albe, 2010) and used as a measure of plasticity on the continuum level. The starting point of micro-plasticity in a continuum mechanical sense can then be identified by a significant increase in dislocation density.

3 Results and Discussion

In the following, we first present in section 3.1 the results of the MD simulations for different interlamellar interface models, considering the effect of temperature, strain rate, and preexisting defects. In subsequent section 3.2, the calculation of TS parameters for cohesive zone models from the results of the MD simulations is pointed out. Finally, in section 3.3, we discuss the results of this procedure in more detail.

3.1 Molecular Dynamics Simulation Results

In the following, the separation behavior of $\gamma /{\alpha }_2$ and $\gamma /\gamma$ interfaces under normal loading in MD simulations is compared for the Type-1 and Type-2 interface models for target temperatures $T\in \left\{300K,500K,700K,900K\right\}$ and for the two different strain rates ${\dot{\epsilon }}_{ZZ}^{glo}={10}^9/s$ and ${\dot{\epsilon }}_{ZZ}^{glo}={10}^{10}/s$ imposed by the velocities enforced on the atoms in the loading region of the simulation model (refer Figure 3).

Type-1 Interface

Figure 4 shows the global stress normal to the interface ${\mathbf{S}}_{ZZ}^{\text{glo}}$ for Type-1 models of the four single interlamellar interfaces according to Eq. 13 vs. the strain of the simulation box. One can observe a monotonous relation between stress and strain for all temperatures, interfaces, and strain rates up to a point, denoted as yield stress, which is characterized by significant nucleation of dislocations. In Figure 4, this point is marked by symbols (triangles for ${\dot{\epsilon }}_{ZZ}^{glo}={10}^{10}/s$ and circles for ${\dot{\epsilon }}_{ZZ}^{glo}={10}^9/s$). Briefly after this point, we observe either immediately or after a short regime of plastic hardening a drop of stress. The yield stresses for all interface types and temperatures are summarized in Table 1. Generally, yield stress decreases with increasing temperature. At any given temperature, the ${\gamma /\alpha }_2$ interface was found to exhibit the highest yield strength, followed by the $\gamma /{\gamma }_{\mathrm{PT}}$ and $\gamma /{\gamma }_{\mathrm{TT}}$ interfaces and finally by the (weakest) $\gamma /{\gamma }_{\mathrm{RB}}$ interface. The observed yield stresses ranged between 5.57 and 12.67 GPa. Of course, one has to keep in mind that these values are theoretical values for ideal, perfect interfaces.

FIGURE 4

FIGURE 4. Stress–strain curves under normal loading for Type-1 models of (A)$\gamma /{\alpha }_2$, (B)$\gamma /{\gamma }_{\mathrm{PT}}$, (C)$\gamma /{\gamma }_{\mathrm{RB}}$, and (D)$\gamma /{\gamma }_{\mathrm{TT}}$ interfaces for different temperatures and strain rates. Here, the inset shows the section of the corresponding interfaces normal to the close-packed direction. The solid circle and triangle markers represent the dislocation nucleation point for the strain rates ${\dot{\epsilon }}_{ZZ}^{glo}={10}^9/s$ (solid lines, circles) and ${\dot{\epsilon }}_{ZZ}^{glo}={10}^{10}/s$ (dashed lines, triangles).

TABLE 1

TABLE 1. Yield stress of Type-1 models for different temperatures and interface types at strain rates of ${\dot{\epsilon }}_{ZZ}^{glo}={10}^9/s$ and ${\dot{\epsilon }}_{ZZ}^{glo}={10}^{10}/s$.

Some interface models ($\gamma /{\gamma }_{\mathrm{PT}}$ and $\gamma /{\gamma }_{\mathrm{TT}}$) show a sharp drop in stress after reaching the peak value at room temperature for low strain rate (refer Figure 4), suggesting a brittle-like fracture failure characterized by interface separation. Oscillations in the stress–strain curves are observed for all the interfaces after failure. Especially, for the higher strain rate (studied here), the vibration effect is pronounced, visible by substantial oscillations of the stress–strain curve, which is understandable noting that these oscillations are dynamic effects increasing thereby the more the conditions are different from the quasi-static case.

In general, one can observe that a reduction of strain rate typically reduces the yield stress for the cases studied herein (see Table 1). This observation reaffirms the reported strain rate sensitivity (Spearot et al., 2009), that is, reduction in strain rate reduces the flow/nucleation stresses, typically for materials with reduced defect activation volumes in the nanoscale regime, for example, nano-twinned copper (Zhu et al., 2007). Especially when such small-volume materials are deformed to the elastic limit, the dislocation nucleates result from strain localization. It is worth emphasizing that the reduced activation volumes impact the nucleation stresses, making them susceptible to strain rate and temperature sensitivity (Zhu et al., 2008). Interestingly, MD simulations of a pure γ-phase show such a strain rate sensitivity by the authors Wu et al. (2016), who also rationalized the observation to some of the above sources (Zhu et al., 2007, 2008; Spearot et al., 2009). Similar to our observation in $\gamma /{\alpha }_2$ interfaces, the γ-phase is the softer phase carrying the most inelastic deformation (Wu et al., 2016). Thus, the above reported strain rate sensitivity also appears suitable to explain our results. For both applied strain rates, we observe an increase in hardening with a rise in temperature, signaling that the dislocation activities giving rise to hardening are thermally driven. Furthermore, the observed strength values are in good agreement with previous studies, especially for $\gamma /{\alpha }_2$ (Li et al., 2019).

Type-2 Interface

Figure 5 shows the global stress–strain curves for the Type-2 interface models, that is, with a preexisting void. The presence of the void type defect reduces the interface strength due to high local stress concentrations around the void and hence also the global load under which interatomic bonds break. Consequently, in comparison to the Type-1 models, the yield stress is reduced for all the interfaces and the high-stress regions near the void become the dominant dislocation sources. Such readily available dislocations lead to a more pronounced hardening after the yield stress has been reached compared to the Type-1 models. As in Type-1 models, the yield stress decreases also in Type-2 models with temperature (see Table 2). And again, it is highest for the $\gamma /{\alpha }_2$ interface, followed by $\gamma /{\gamma }_{\mathrm{PT}}$, $\gamma /{\gamma }_{\mathrm{TT}}$, and finally the (weakest) $\gamma /{\gamma }_{\mathrm{RB}}$ interface at all temperatures studied. Generally, for low strain rate simulations, the void in Type-2 models reduces yield stress for all the four interfaces by around 35% for higher temperatures and more than 40% for 300 K. In Figure 5, all the Type-2 interface models show similar trends at the strain rate ${\dot{\epsilon }}_{ZZ}^{glo}={10}^{10}/s$. That is, the tensile stress increases further after the onset of dislocation nucleation (triangles in Figure 5) followed by a smooth drop indicating a rather ductile fracture process. This interpretation, where naturally substantial energy dissipation during the fracture process can be expected from its ductile nature, is also supported by the absence of major vibrations in the stress–strain curve after separation (in contrast to Type-1 models, refer Figure 4). For the strain rate ${\dot{\epsilon }}_{ZZ}^{glo}={10}^9/s$, after reaching the yield stress, typically a short dip in stress occurs followed by hardening. Whether the drop is more pronounced or the hardening (or both of them even level each other out) depends on the dislocation and void evolution. The void growth reduces the stress (softening), while the nucleation of dislocations leads to hardening. The Type-2 models illustrated in Figure 5 differ in two important aspects from the previously shown Type-1 models. First, the strain values where dislocation nucleation starts (yield strain, indicated by triangles and circles) are nearly independent of temperature for most of the interfaces studied. The reason may be the dominant role of the preexisting void for dislocation nucleation, which reduces the effect of temperature. Second, the hardening from dislocation-driven plasticity decreases with increasing temperature due to the pronounced dislocation mobility. More information on the dislocation evolution is shown in Figures 6 and will be discussed in the subsequent subsections. The observed high-temperature deformation behavior, along with the underlying dislocation-mediated mechanisms, is discussed in more detail in this article’s Supplementary Material.

FIGURE 5

FIGURE 5. Stress–strain curves under normal loading for Type-2 models (with preexisting void) of (A)$\gamma /{\alpha }_2$, (B)$\gamma /{\gamma }_{\mathrm{PT}}$, (C)$\gamma /{\gamma }_{\mathrm{RB}}$, and (D)$\gamma /{\gamma }_{\mathrm{TT}}$ interfaces for different temperatures and strain rates. Here, the inset shows the section of the corresponding interfaces normal to the close-packed direction. The solid circle and triangle markers represent the dislocation nucleation point for the strain rates ${\dot{\epsilon }}_{ZZ}^{glo}={10}^9/s$ (solid lines, circles) and ${\dot{\epsilon }}_{ZZ}^{glo}={10}^{10}/s$ (dashed lines, triangles).

TABLE 2

TABLE 2. Yield stress of Type-2 models for different temperatures and interface types at strain rates of ${\dot{\epsilon }}_{ZZ}^{glo}={10}^9/s$ and ${\dot{\epsilon }}_{ZZ}^{glo}={10}^{10}/s$.

FIGURE 6

FIGURE 6. Traction–separation curves (solid lines) from MD simulations with strain rate ${\dot{\epsilon }}_{ZZ}^{glo}={10}^9/s$ for different temperatures and Type-1 models of three $\gamma /\gamma$ interfaces (i.e., (A)$\gamma /{\gamma }_{\mathrm{PT}}$, (B)$\gamma /{\gamma }_{\mathrm{RB}}$, and (C)$\gamma /{\gamma }_{\mathrm{TT}}$). The evolution of the dislocation density is illustrated by the dashed lines. The inset figures show the traction–separation curves in the yielding regime in greater detail.

3.2 Traction–Separation Behavior of the Interfaces

As discussed above, TS curves, which form the core of cohesive zone models, can be extracted from MD simulations by fitting the parameters in the equations of the model to the $〚u_n〛$-${\mathbf{S}}_{ZZ}^{\text{loc}}$ curves observed in MD simulations. Since the investigation of strain rate effects is (at least in this high strain rate regime) of minor practical relevance, this process is discussed only for the strain rate ${\dot{\epsilon }}_{ZZ}^{glo}={10}^9{s}^{-1}$ in this section.

Type-1 Interface Models

First of all, it is important to note that not all configurations fail at the interface, namely, for the $\gamma /{\alpha }_2$ interface, fracture occurs in the interior of the γ-phase at most of the temperatures. In this article’s Supplementary Material, the evolution of the simulation domain during the imposed deformation is illustrated (refer Supplementary Figures S1, S2), which clearly shows the location and the evolution of voids and complete failure for all configurations at $T=900\mathrm{K}$. Hence, in case of Type-1 models, we focused on $\gamma /\gamma$ interfaces only, which in general show interface fracture.¹ For the $\gamma /{\alpha }_2$ interface, it is in case of a Type-1 interface model neither possible nor of interest to identify a full TS curve. Rather it is reasonable to model this interface on the continuum level as an ideal connection between a γ-phase and an ${\alpha }_2$-phase and have an additional continuum damage model for the bulk γ-TiAl crystal, which is, however, outside the scope of this investigation.

The atomistically computed local TS curves are shown in Figure 6 for the different $\gamma /\gamma$ interfaces and temperatures. Together with the TS curves, the evolution of the dislocation density is also shown, indicating an onset of plastic deformation close to the peak traction. After this peak, the interface traction is decreasing, indicating the degradation of the load-bearing ability of the sample. The substantial differences between the TS curves considered in this section and the figures (refer Figures 4, 5) shown in the previous section is a consequence of the fact that in the previous section, we considered the overall load and deformation of the global simulation domain. In contrast, in this section, we focus on the local interface separation, according to Eq. 16 with effects of the bulk deformation already eliminated. Interestingly, eliminating these effects significantly reduces the size of the elastic regime. One may also notice that the separation is not monotonically increasing in the majority of TS curves. The reason for this behavior may be that the dislocation evolution is not homogeneous throughout the structure; in the first stage after the peak stress, dislocations may sometimes move more in the bulk region, sometimes more in the local region. The respective other region is then unloaded for a short time. However, this only happens in the first 10% of stress reduction and does not change the big picture substantially. The peak atomic tractions in the loading direction decrease with increasing temperature (Figure 6), similarly as in the global stress–strain curves discussed in the previous section. The regime directly after the stress drop depends on temperature as well, as discussed in the Supplementary Material for the global stress–strain behavior.

The shape of the TS law has been fitted to the MD results shown in Figure 6, assuming a nonhardening interface for simplicity due to the brittle fracture observed in Figure 6. Anyway, for the initially undamaged interface, damage starts almost together with plasticity. Eventually, the following parameters had to be identified: the elastic parameter $C_n$ from Eq. 4, the initial yield strength $t_0$ from Eq. 5, and the three damage parameters ${\alpha }_{\text{nucl}},{\alpha }_{\text{fail}},\text{ and }p$ from Eq. 11. Parameter fitting was automated using the optimization module of Scipy (2020). The results for the temperature $T=900$ K are shown in Figure 7; the results for the other temperatures are given in Supplementary Material.

FIGURE 7

FIGURE 7. Defect-free case (Type-1 interface model): traction–separation curves as observed in MD simulations (solid black lines) compared to the ones of cohesive zone models with parameters fitted to the results of the MD simulations (solid blue lines) at a temperature $T=900$ K and strain rate ${\dot{\epsilon }}_{ZZ}^{glo}={10}^9/s$ for the three different γ/γ interfaces (i.e., (A)$\gamma /{\gamma }_{\mathrm{PT}}$, (B)$\gamma /{\gamma }_{\mathrm{RB}}$, and (C)$\gamma /{\gamma }_{\mathrm{TT}}$). The dashed lines denote the normalized dislocation density ρ observed in MD simulations (black) and the normalized damage-corrected interface plasticity $\alpha /\left[1\mathrm{-}d\right]$ in the cohesive zone model (blue).

The agreement between MD simulations and the fitted cohesive zone models depicted in Figure 8 is in each case very good. Remarkably, this holds not only for the TS curves themselves but also for the normalized dislocation density in the MD simulations and the corresponding damage-corrected plastic separation $\left(1-d\right)\alpha$. This indicates that our model is not purely phenomenological but appears to capture at least to some extent the physics. The fitted parameters for the cohesive zone model are given in Table 3 for a temperature of $900$ K. The values for the other temperatures are listed in Supplementary Table S1 (refer Supplementary Material). It is worth mentioning that the parameters depend almost linearly on temperature enabling an easy transfer of our results to other temperatures in the same range.

FIGURE 8

FIGURE 8. Traction‐separation curves (solid lines) from MD simulations with strain rate $T=900$ K and strain rate ${\dot{\epsilon }}_{ZZ}^{glo}={10}^9/s$ for different temperatures and Type-2 models of different interfaces (i.e., (A)$\gamma /{\alpha }_2$, (B)$\gamma /{\gamma }_{\mathrm{PT}}$, and (C)$\gamma /{\gamma }_{\mathrm{RB}}$, (D)$\gamma /{\gamma }_{\mathrm{TT}}$). The evolution of the dislocation density is illustrated by the dashed lines. The inset figures show the traction‐separation curves in the yielding regime in greater detail.

TABLE 3

TABLE 3. Parameters of cohesive zone model for Type-1 interface models determined from fitting the results of MD simulations at $T=900$ K

Type-2 Interface Models

Figure 7 shows the TS curves for the Type-2 interface models with the preexisting void for different temperatures as observed in MD simulations. These curves generally exhibit three regimes: in the first one, traction increases up to a certain peak value. In the second one, the traction remains on a high level with some fluctuation (which is not unusual for MD simulations), followed by a third regime where it finally decreases. While the elastic behavior and the general temperature dependence are similar to the Type-1 models, also characteristic differences are observed.

First of all, in contrast to Type-1 interface models, failure is always observed at the interface itself, also for the $\gamma /{\alpha }_2$ interface, so that this interface can also be addressed by a cohesive zone model. Second, the peak traction in Type-2 interface models is generally by around 35–50% lower than for the defect-free Type-1 models. Third, in Type-2 models, dislocation density increases more gradually after the onset of plasticity, and instead of the drop of the traction directly after the onset of yielding, it remains high for some separation. This indicates a significant regime of interface plasticity without damage evolution.

As for the Type-1 models, the parameters of a cohesive zone model were fitted to the results of MD simulations for the Type-2 models (see Figure 9 for $T=900$ K). Again a very good agreement between the results of the MD simulations and the fitted cohesive zone model is observed. In particular, the damage-free plasticity region with a constant traction corresponds extraordinarily well with the increasing branch of the dislocation density (black dash-dotted line), which is also reflected by the damage-corrected plastic separation $\left(1-d\right)\alpha$ (blue dash-dotted curve). The parameter values fitted at $T=900\mathrm{K}$ are given in Table 4; those for the other temperatures are shown in Supplementary Table S2 (refer Supplementary Material). Curve fitting for the other temperatures worked similarly well as shown in Supplementary Material (refer Supplementary Figures S6–S8).

FIGURE 9

FIGURE 9. Interfaces with defects (Type-2 interface model): traction–separation curves as observed in MD simulations (solid black lines) compared to the ones of cohesive zone models with parameters fitted to the results of the MD simulations (solid blue lines) at a temperature $T=900$ K and strain rate ${\dot{\epsilon }}_{ZZ}^{glo}={10}^9/s$ for (A)$\gamma /{\alpha }_2$, (B)$\gamma /{\gamma }_{\mathrm{PT}}$, (C)$\gamma /{\gamma }_{RB}$, and (D)$\gamma /{\gamma }_{\mathrm{TT}}$ interfaces. The dashed lines denote the normalized dislocation density ρ observed in MD simulations (black) and the normalized damage-corrected interface plasticity $\alpha /\left[1\mathrm{-}d\right]$ in the cohesive zone model (blue).

TABLE 4

TABLE 4. Parameters of cohesive zone model for Type-2 interface models determined from fitting the results of MD simulations $T=900$ K.

3.3 Discussion

From the parameters identified for the different temperatures and interface types (Tables 3, 4; Supplementary Tables S1, S2 in Supplementary Material), it is possible to calculate the separation energies ${\text{Γ}}_0$ based on

${\text{Γ}}_0={{{\displaystyle \int }}^{\text{}}}_0^{{〚u_n〛}^{\text{fail}}}{t}_n\text{d}〚u_n〛,(17)$

with ${〚u_n〛}^{\text{fail}}$ as the separation at which failure of the interface is observed and $t_n=S_{ZZ}^{\text{loc}}$. The resulting values of ${\text{Γ}}_0$ for all temperatures are given in Table 5. It turns out that the energies vary slightly for the different interfaces, and in general, they have a decreasing trend with temperature. Remarkably, the prexisiting void in the Type-2 interface models seems to decrease the interface strength by around 30% (see above), but for the interface separation energy, no such clear trend is found in our simulations.

TABLE 5A

TABLE 5A. Separation energies ${\text{Γ}}_0$ [J/m²] calculated by Eq. 17 for different temperatures and (a) Type-1 interface model and (b) Type-2 interface model.

TABLE 5B

TABLE 5B.

A quantitative comparison of the cohesive zone model parameters with values from the literature is difficult due to the sparse data available. However, there are a few articles dealing with the strength of lamellar interfaces in γ-TiAl. A thorough combined experimental and numerical study was proposed by Werwer et al. (2007), where precracked PST (polysynthetically twinned) crystal specimens made of lamellar γ-TiAl were tested in different orientations at room temperature. Cohesive model parameters (i.e., cohesive strength and separation energy) for a TS law were determined by fitting fracture simulations to the experimental data. A very low cohesive strength of about 100 MPa was found, but the identified separation energies $\left({\text{Γ}}_{\text{0}}=3\dots 6\text{J}/{\text{m}}^2\right)$ are comparable with our findings (ref. Table 5A,5B).

An explanation why the strength observed in Werwer et al. (2007) is so much (roughly by a factor of 50) lower than the one determined herein ($t_0$ in Table 4) may be given on the basis of a recent study of Möller et al. (2018), who performed MD simulations of a crack tip in a brittle tungsten metal to calculate TS law parameters. They concluded that for a simulation box size comparable to the one considered in our study, the stresses are still within the K-dominated field. That is, the opening stress ahead of the crack tip is proportional to the square root of the distance, and a TS law with constant values may not be suitable to describe the interface strength properly. By contrast, they concluded that the separation energy may be largely independent of the size of the simulated domain already for a domain size comparable to the one considered herein.

In another investigation by Krull and Yuan (2011) on a similar but more ductile material, aluminum, the authors did not consider the size effect and also found unrealistically high values for the cohesive strength of more than 10 GPa.

For the computational results presented in this article, this means that the separation energies of the interface are most likely suitable for utilization in a continuum level finite element simulation where the interface may be discretized by interface elements. By contrast, for the cohesive strength, one should probably use much lower values than reported herein if the length of the cohesive finite elements used on the continuum scale is larger than the size of the K-dominated zone.

4 Summary and Conclusion

To summarize, the high-temperature deformation behavior of the Type-1 (defect-free) and Type-2 (preexisting void) single lamellar interface models (one $\gamma /{\alpha }_2$ and three $\gamma /\gamma$) has been studied here under tensile load using MD simulations. The results of these simulations were used to fit the parameters of the traction separation laws of cohesive zone models, which can then be used for continuum level models, for example, using the finite element methods.

From this study, we can draw a number of conclusions on the MD simulations alone and some further conclusions on the derivation of TS curves from the atomistic simulations. In our simulations, the $\gamma /{\alpha }_2$ interface was found to be the strongest, followed by the $\gamma /{\gamma }_{\mathrm{PT}}$, $\gamma /{\gamma }_{\mathrm{TT}}$, and ${\gamma /\gamma }_{\mathrm{RB}}$ interfaces. Particularly, the $\gamma /{\alpha }_2$ interface for the ideal (defect-free) structure was found to be stronger than the adjacent bulk material, and failure never occurred at the interface itself, but first in the bulk material. Our simulations allowed for the first time a quantification of the reduction of the yield strength as temperature increases for such single interlamellar interfaces in TiAl.

In general, we observed in our simulations a quasi-brittle behavior with a sudden failure for the defect-free interface models, while a pronounced plasticity was observed for the more realistic interface models with a preexisting void. Beside the failure behavior, the evolution of the dislocation density was quantified by our simulations and it turned out that a dislocation burst at a high stress occurs for the defect-free case, while the dislocation density evolves gradually at lower stresses for interfaces with preexisting voids.

The evaluation of the damage behavior of the interface itself (with the effect of the bulk deformation having been eliminated by Eq. 16) revealed a TS behavior of the interface without significant elastic separation. Rather, the main contribution of the overall interface separation until failure came from interface plasticity and interface damage, which could be clearly distinguished in the TS law. Without elimination of the bulk deformation in our evaluation, a very different picture would have resulted where the TS curves would have exhibited a significant elastic regime. It is important to underline this because to the authors’ best knowledge, comparable studies of other groups where traction–separation laws were deduced from MD simulations did not eliminate the effect of bulk deformation from their calculations as we did above. Our study indeed suggests that doing so makes a substantial difference and would therefore be worthwhile to be taken into consideration also by other groups in future work to characterize the physical properties of interfaces more accurately. Matching cohesive zone models with the results of MD simulations, plasticity could be correlated with an increase in dislocation density in the MD simulations, which underlines that the differentiation between elastic and plastic separation in our cohesive zone models is not simply a mathematical trick to cover a wider range of functional relations but that it captures indeed a specific form of physics that can be clearly observed on the atomic scale. Particularly, the damage-corrected plastic separation $\left(1-d\right)\alpha$ parameter compares qualitatively well with the dislocation densities from the MD simulations. In fact, this study is to the authors’ best knowledge the first one that deduces not only TS curves from MD simulations but also captures the inelastic deformation of the interface in a physically apparently meaningful way.

We also studied the effect of initial defects at the interface. They were found to decrease the interface strength as expected. However, for the interface separation energy, now such clear trend is observed in our MD simulations.

A comparison with results of other groups reported in the literature suggests that our MD simulations are suitable for identifying the interface separation energy correctly. By contrast, our results for the interface strength as observed in MD simulations should be taken with caution. They seem to overestimate the interface strength compared to experimental observations, possibly due to too small a simulation domain typical for nanoscale regime often covered in MD. This problem has also been encountered by others who tried to identify interface strengths by MD simulations, for example, for aluminum. Apparently, further research is required to overcome this key problem in the multi-scale simulation of materials.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

HG is responsible for all MD simulations; IS is responsible for all cohesive zone modeling issues; CC supervised the modeling work and contributed to its interpretation as well as the preparation of the manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The financial support by the IDEA framework within the Helmholtz-Zentrum Geesthacht is gratefully acknowledged. The authors thank Rebecca Janisch for providing some of the source construction scripts and for the fruitful discussions.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmats.2020.602567/full#supplementary-material.

Abbreviations

APB, Antiphase Boundary; CSF, complex stacking fault; CSP, centrosymmetry parameter; EAM, embedded-atom method; FL, fully lamellar; MD, molecular dynamics; PT, pseudo twin; RB, rotational boundary; SISF, super lattice intrinsic stacking fault; TS, traction separation; TT, true twin.

Footnotes

¹except for the $\gamma /{\gamma }_{\mathrm{RB}}$ interface at 500 K (green curve in Figure 6B) where the strength of the interface also appears to be close to the strength of the bulk material.

References

Appel, F., Brossmann, U., Christoph, U., Eggert, S., Janschek, P., Lorenz, U., et al. (2000). Recent progress in the development of gamma titanium aluminide alloys. Adv. Eng. Mater. 2, 699–720. doi:10.1002/1527-2648(200011)2:11<699::aid-adem699>3.0.co;2-j

CrossRef Full Text | Google Scholar

Appel, F., Clemens, H., and Fischer, F. D. (2016). Modeling concepts for intermetallic titanium aluminides. Prog. Mater. Sci. 81, 55–124. doi:10.1016/j.pmatsci.2016.01.001

CrossRef Full Text | Google Scholar

Appel, F., Paul, J. D. H., and Oehring, M. (2011). Gamma titanium aluminide alloys: science and technology. Hoboken, New Jersey: John Wiley & Sons.

CrossRef Full Text | Google Scholar

Barbi, N., Diologent, F., Goodall, R., and Mortensen, A. (2012). Fracture of convoluted and lamellar α2+ γ TiAl alloys. Intermetallics 22, 176–188. doi:10.1016/j.intermet.2011.11.001

CrossRef Full Text | Google Scholar

Bewlay, B. P., Nag, S., Suzuki, A., and Weimer, M. J. (2016). TiAl alloys in commercial aircraft engines. Mater. at High. Temp. 33, 549–559. doi:10.1080/09603409.2016.1183068

CrossRef Full Text | Google Scholar

Cai, W., Li, J., and Yip, S. (2012). “Molecular dynamics,” in Comprehensive nuclear materials. Editors J. M. K. Rudy (Oxford, UK: Elsevier), 249–265

CrossRef Full Text | Google Scholar

Cao, H., Rui, Z., Chen, W., Feng, R., and Yan, C. (2019). Crack propagation mechanism of γ-TiAl alloy with pre-existing twin boundary. Sci. China Technol. Sci. 62, 1605–1615. doi:10.1007/s11431-017-9307-5

CrossRef Full Text | Google Scholar

Clemens, H., and Mayer, S. (2013). Design, processing, microstructure, properties, and applications of advances intermetallic TiAl alloys. Adv. Eng. Mater. 15, 191–215. doi:10.1002/adem.201200231

CrossRef Full Text | Google Scholar

Cornec, A., Scheider, I., and Schwalbe, K.-H. (2003). On the practical application of the cohesive model. Eng. Fract. Mech. 70, 1963–1987. doi:10.1016/s0013-7944(03)00134-6

CrossRef Full Text | Google Scholar

Dandekar, C. R., and Shin, Y. C. (2011). Molecular dynamics based cohesive zone law for describing al–sic interface mechanics. Compos. Appl. Sci. Manuf. 42, 355–363. doi:10.1016/j.compositesa.2010.12.005

CrossRef Full Text | Google Scholar

Daw, M. S., and Baskes, M. I. (1984). Embedded-atom method: derivation and application to impurities, surfaces, and other defects in metals. Phys. Rev. B 29, 6443. doi:10.1103/physrevb.29.6443

CrossRef Full Text | Google Scholar

Dimiduk, D. M. (1999). Gamma titanium aluminide alloys – an assessment within the competition of aerospace structural materials. Mater. Sci. Eng.: A 263, 281–288. doi:10.1016/s0921-5093(98)01158-7

CrossRef Full Text | Google Scholar

Ding, J., Tian, Y., Wang, L.-s., Huang, X., Zheng, H.-r., Song, K., et al. (2019). Micro-mechanism of the effect of grain size and temperature on the mechanical properties of polycrystalline TiAl. Comput. Mater. Sci. 158, 76–87. doi:10.1016/j.commatsci.2018.10.019

CrossRef Full Text | Google Scholar

Feng, R., Song, W., Li, H., Qi, Y., Qiao, H., and Li, L. (2018). Effects of annealing on the residual stress in γ-TiAl alloy by molecular dynamics simulation. Materials. 11, 1025. doi:10.3390/ma11061025

CrossRef Full Text | Google Scholar

Feng, R., Wang, M., Li, H., Qi, Y., Wang, Q., and Rui, Z. (2019). Micromechanism of cold deformation of two-phase polycrystalline Ti-Al alloy with void. Materials 12, 184. doi:10.3390/ma12010184

CrossRef Full Text | Google Scholar

Gupta, P., Pal, S., and Yedla, N. (2016). Molecular dynamics based cohesive zone modeling of al (metal)–cu50zr50 (metallic glass) interfacial mechanical behavior and investigation of dissipative mechanisms. Mater. Des. 105, 41–50. doi:10.1016/j.matdes.2016.05.054

CrossRef Full Text | Google Scholar

Hillerborg, A., Modéer, M., and Petersson, P.-E. (1976). Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement Concr. Res. 6, 773–781. doi:10.1016/0008-8846(76)90007-7

CrossRef Full Text | Google Scholar

Hui, C., Zhiyuan, R., Wenke, C., Ruicheng, F., and Changfeng, Y. (2018). Deformation mechanisms in nanotwinned γ-TiAl by molecular dynamics simulation. Mol. Simulat. 44, 1489–1500. doi:10.1080/08927022.2018.1496248

CrossRef Full Text | Google Scholar

Ji, Z.-W., Lu, S., Hu, Q.-m., Kim, D., Yang, R., and Vitos, L. (2018). Mapping deformation mechanisms in lamellar titanium aluminide. Acta Mater. 144, 835–843. doi:10.1016/j.actamat.2017.11.028

CrossRef Full Text | Google Scholar

Jin, Z., and Gray, G. T. (1997). Experimental determination of domain orientations and domain orientation relationships across lamellar interfaces in polysynthetically twinned TiAl crystals. Mate. Sci. Eng.: A 231, 62–71. doi:10.1016/s0921-5093(97)00036-1

CrossRef Full Text | Google Scholar

Kanani, M., Hartmaier, A., and Janisch, R. (2016). Stacking fault based analysis of shear mechanisms at interfaces in lamellar TiAl alloys. Acta Mater. 106, 208–218. doi:10.1016/j.actamat.2015.11.047

CrossRef Full Text | Google Scholar

Kim, Y.-K., Kim, H.-K., Jung, W.-S., and Lee, B.-J. (2016). Atomistic modeling of the Ti-Al binary system. Comput. Mater. Sci. 119, 1–8. doi:10.1016/j.commatsci.2016.03.038

CrossRef Full Text | Google Scholar

Kim, Y.-W., and Dimiduk, D. M. (1991). Progress in the understanding of gamma titanium aluminides. JOM 43, 40–47. doi:10.1007/bf03221103

CrossRef Full Text | Google Scholar

Kiselev, S. P., and Zhirov, E. V. (2014). Molecular-dynamics simulation of the synthesis of intermetallic Ti-Al. Intermetallics 49, 106–114. doi:10.1016/j.intermet.2014.01.008

CrossRef Full Text | Google Scholar

Krull, H., and Yuan, H. (2011). Suggestions to the cohesive traction–separation law from atomistic simulations. Eng. Fract. Mech. 78, 525–533. doi:10.1016/j.engfracmech.2009.12.014

CrossRef Full Text | Google Scholar

Li, H., and Chandra, N. (2003). Analysis of crack growth and crack-tip plasticity in ductile materials using cohesive zone models. Int. J. Plast. 19, 849–882. doi:10.1016/s0749-6419(02)00008-6

CrossRef Full Text | Google Scholar

Li, P.-t., Yang, Y.-Q., Xia, Z., Luo, X., Jin, N., Gao, Y., et al. (2017). Molecular dynamic simulation of nanocrystal formation and tensile deformation of TiAl alloy. RSC Adv. 7, 48315–48323. doi:10.1039/c7ra10010h

CrossRef Full Text | Google Scholar

Li, W., Yin, Y., Xu, Q., Zhou, J., Nan, H., Ji, X., et al. (2019). Tensile behavior of γ/α2 interface system in lamellar TiAl alloy via molecular dynamics. Comput. Mater. Sci. 159, 397–402. doi:10.1016/j.commatsci.2018.12.043

CrossRef Full Text | Google Scholar

Li, W., Yu, W., Xu, Q., Zhou, J., Nan, H., Yin, Y., et al. (2020). Effects of γ/γ interfaces in TiAl lamellae subjected to uniaxial tensile loading. Comput. Mater. Sci. 172, 109361. doi:10.1016/j.commatsci.2019.109361

CrossRef Full Text | Google Scholar

Möller, J. J., Bitzek, E., Janisch, R., ul Hassan, H., and Hartmaier, A. (2018). Fracture ab initio: a force-based scaling law for atomistically informed continuum models. J. Mater. Res. 33, 3750–3761. doi:10.1557/jmr.2018.384

CrossRef Full Text | Google Scholar

Maugin, G. A. (1992). The thermomechanics of plasticity and fracture. Cambridge, UK: Cambridge University Press, Vol. 7.

CrossRef Full Text | Google Scholar

Paggi, M., and Wriggers, P. (2011). A nonlocal cohesive zone model for finite thickness interfaces – part i: mathematical formulation and validation with molecular dynamics. Comput. Mater. Sci. 50, 1625–1633. doi:10.1016/j.commatsci.2010.12.024

CrossRef Full Text | Google Scholar

Park, K., and Paulino, G. (2011). Cohesive zone models: a critical review of traction-separation relationships across fracture surface. Appl. Mech. Rev. 64, 060802. doi:10.1115/1.4023110

CrossRef Full Text | Google Scholar

Plimpton, S. (1995). Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 117, 1–19. doi:10.1006/jcph.1995.1039

CrossRef Full Text | Google Scholar

Scheider, I. (2018). Numerical simulation of material separation using cohesive zone models. Singapore: Springer, Vol. 1, 1–42

CrossRef Full Text | Google Scholar

Scheider, I. (2020). [Dataset]. Code for Quantifying the high-temperature separation behavior of lamellar interfaces in γ-TiAl under tensile loading by molecular dynamics. doi:10.5281/zenodo.4066828

CrossRef Full Text | Google Scholar

Scheider, I. (2009a). Derivation of separation laws for cohesive models in the course of ductile fracture. Eng. Fract. Mech. 76, 1450–1459. doi:10.1016/j.engfracmech.2008.12.006

CrossRef Full Text | Google Scholar

Scheider, I. (2009b). Micromechanical based derivation of traction-separation laws for cohesive model simulations. Procedia Eng. 1 (1), 17–21. doi:10.1016/j.proeng.2009.06.006

CrossRef Full Text | Google Scholar

Schnabel, J., and Bargmann, S. (2017). Accessing colony boundary strengthening of fully lamellar TiAl alloys via micromechanical modeling. Materials 10, 896. doi:10.3390/ma10080896

CrossRef Full Text | Google Scholar

Schnabel, J., and Scheider, I. (Forthcoming 2020). Crystal plasticity modeling of creep in alloys with lamellar microstructures at the example of fully lamellar TiAl. Front. Mater. doi:10.3389/fmats.2020.581187

CrossRef Full Text | Google Scholar

Scipy (2020). Optimization and root finding (scipy.optimize) — scipy v1.4.1 reference guide. Available at: https://docs.scipy.org/doc/scipy/reference/optimize.html.

CrossRef Full Text | Google Scholar

Simonovski, I., and Cizelj, L. (2015). Cohesive zone modeling of intergranular cracking in polycrystalline aggregates. Nucl. Eng. Des. 283, 139–147. doi:10.1016/j.nucengdes.2014.09.041

CrossRef Full Text | Google Scholar

Spearot, D. E., Tschopp, M. A., and McDowell, D. L. (2009). Orientation and rate dependence of dislocation nucleation stress computed using molecular dynamics. Scripta Mater. 60, 675–678. doi:10.1016/j.scriptamat.2008.12.037

CrossRef Full Text | Google Scholar

Stukowski, A., and Albe, K. (2010). Dislocation detection algorithm for atomistic simulations. Model. Simulat. Mater. Sci. Eng. 18, 025016. doi:10.1088/0965-0393/18/2/025016

CrossRef Full Text | Google Scholar

Stukowski, A. (2009). Visualization and analysis of atomistic simulation data with ovito–the open visualization tool. Model. Simulat. Mater. Sci. Eng. 18, 015012. doi:10.1088/0965-0393/18/1/015012

CrossRef Full Text | Google Scholar

Thompson, A. P., Plimpton, S. J., and Mattson, W. (2009). General formulation of pressure and stress tensor for arbitrary many-body interaction potentials under periodic boundary conditions. J. Chem. Phys. 131, 154107. doi:10.1063/1.3245303

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, Y., Gao, H., and Bower, A. F. (2009). Numerical simulations of crack deflection at a twist-misoriented grain boundary between two ideally brittle crystals. J. Mech. Phys. Solid 57, 1865–1879. doi:10.1016/j.jmps.2009.07.007

CrossRef Full Text | Google Scholar

Werwer, M., and Cornec, A. (2000). Numerical simulation of plastic deformation and fracture in polysynthetically twinned (PST) crystals of TiAl. Comput. Mater. Sci. 19, 97–107. doi:10.1016/s0927-0256(00)00144-0

CrossRef Full Text | Google Scholar

Werwer, M., Kabir, R., Cornec, A., and Schwalbe, K.-H. (2007). Fracture in lamellar TiAl simulated with the cohesive model. Eng. Fract. Mech. 74, 2615–2638. doi:10.1016/j.engfracmech.2006.09.022

CrossRef Full Text | Google Scholar

Wu, H., Xu, D., Wang, H., and Yang, R. (2016). Molecular dynamics simulation of tensile deformation and fracture of γ-TiAl with and without surface defects. J. Mater. Sci. Technol. 32, 1033–1042. doi:10.1016/j.jmst.2015.12.001

CrossRef Full Text | Google Scholar

Xie, Z.-C., Gao, T.-H., Guo, X.-T., and Xie, Q. (2015). Molecular dynamics simulation of nanocrystal formation and deformation behavior of Ti₃Al alloy. Comput. Mater. Sci. 98, 245–251. doi:10.1016/j.commatsci.2014.10.012

CrossRef Full Text | Google Scholar

Yamakov, V., Saether, E., Phillips, D. R., and Glaessgen, E. H. (2006). Molecular-dynamics simulation-based cohesive zone representation of intergranular fracture processes in aluminum. J. Mech. Phys. Solid 54, 1899–1928. doi:10.1016/j.jmps.2006.03.004

CrossRef Full Text | Google Scholar

Zhou, L., Hsiung, L., and Huang, H. (2004). Nucleation and propagation of deformation twin in polysynthetically twinned TiAl. Comput. Model. Eng. Sci. 16, 245–252. doi:10.3970/CMES.2004.006.245

CrossRef Full Text | Google Scholar

Zhou, X. W., Moody, N. R., Jones, R. E., Zimmerman, J. A., and Reedy, E. D. (2009). Molecular-dynamics-based cohesive zone law for brittle interfacial fracture under mixed loading conditions: effects of elastic constant mismatch. Acta Mater. 57, 4671–4686. doi:10.1016/j.actamat.2009.06.023

CrossRef Full Text | Google Scholar

Zhou, X. W., Zimmerman, J. A., Reedy, E. D., and Moody, N. R. (2008). Molecular dynamics simulation based cohesive surface representation of mixed mode fracture. Mech. Mater. 40, 832–845. doi:10.1016/j.mechmat.2008.05.001

CrossRef Full Text | Google Scholar

Zhu, T., Li, J., Samanta, A., Kim, H. G., and Suresh, S. (2007). Interfacial plasticity governs strain rate sensitivity and ductility in nanostructured metals. Proc. Natl. Acad. Sci. U S A 104, 3031–3036. doi:10.1073/pnas.0611097104

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, T., Li, J., Samanta, A., Leach, A., and Gall, K. (2008). Temperature and strain-rate dependence of surface dislocation nucleation. Phys. Rev. Lett. 100, 025502. doi:10.1103/physrevlett.100.025502

PubMed Abstract | CrossRef Full Text | Google Scholar

Zope, R. R., and Mishin, Y. (2003). Interatomic potentials for atomistic simulations of the Ti-Al system. Phys. Rev. B 68, 024102. doi:10.1103/physrevb.68.024102

CrossRef Full Text | Google Scholar