Coupling a discrete twin model with cohesive elements to understand twin-induced fracture

The interplay between twinning and fracture in metals under deformation is an open question. The plastic strain concentration created by twin bands can induce large stresses on the grain boundaries. We present simulations in which a continuum model describing discrete twins is coupled with a crystal plasticity finite element model and a cohesive zone model for intergranular fracture. The discrete twin model can predict twin nucleation, propagation, growth and the correct twin thickness. Therefore, the plastic strain concentration in the twin band can be modelled. The cohesive zone model is based on a bilinear traction-separation law in which the damage is caused by the normal stress on the grain boundary. An algorithm is developed to generate interface elements at the grain boundaries that satisfy the traction-separation law. The model is calibrated by comparing polycrystal simulations with the experimentally observed strain to failure and maximum stress. The dynamics of twin and crack nucleation have been investigated. First, twins nucleate and propagate in a grain, then, microcracks form near the intersection between twin tips and grain boundaries. Microcracks appear at multiple locations before merging. A propagating crack can nucleate additional twins starting from the grain boundary, a few micrometres away from the original crack nucleation site. This model can be used to understand which type of texture is more resistant against crack nucleation and propagation in cast metals in which twinning is a deformation mechanism. The code is available online at https://github.com/TarletonGroup/CrystalPlasticity.


Introduction
Intergranular fracture is a type of brittle fracture, which involves the nucleation and propagation of cracks along grain boundaries. It is a dominant failure mechanism in as-cast metals (Powell 1994), ultrafine-grained metals (Pippan and Hohenwarter 2016), during stress corrosion (Birbilis and Hinton 2011) and hydrogen embrittlement (Barrera et al. 2018;Elmukashfi et al. 2020). It takes place in materials in which the stress required to debond the grain boundary interface is lower than the stress to debond atomic planes inside the grains (Jiang et al. 2015) or to nucleate and grow microscopic voids (Tvergaard 1981;Cocks and Ashby 1982).
Plastic deformation, induced by the motion of dislocations and by twinning (Christian 2002), has an important effect on fracture. The interplay between plasticity and fracture is complex, because it depends on the load condition (Sistaninia and Niffenegger 2015;Grilli and Koslowski 2018), on the microstructure (Qian et al. 2018) and on the presence of pre-existing cracks (Duarte et al. 2018). Larger plastic deformation leads to stress relaxation and can retard the propagation of pre-existing cracks (Grilli and Koslowski 2019). In this case, a larger strain to failure is observed (Woelke et al. 2015). However, plasticity can also cause damage initiation because its anisotropic nature leads to stress variations at the microscale (Krupp 2007).
Plastic deformation at the scale of the grain-size in metals is an inhomogeneous phenomenon. It is concentrated into slip bands or twins, depending on the stacking fault energy (Frøseth et al. 2005), and can favour crack nucleation (Tanaka and Mura 1981). However, shear bands emanating from the notch of a tensile sample can delay fracture (Li et al. 2016). A plastically deforming band that impinges on a grain boundary leads to a normal stress that decreases with the band thickness (Sauzay and Moussa 2013). High stress exists in the vicinity of a twin tip, which can induce slip ahead of it when the twin tip terminates inside the grain (Sleeswyk 1962).
There is experimental evidence of the correlation between twinning and fracture (Christian and Mahajan 1995). Most metals that exhibit brittle fracture at low temperature also deform by mechanical twinning (O'Neill 1926). Microcracks in TWIP steels form at intersecting twins and propagate both along grain boundaries and along twin interfaces (Koyama et al. 2013). Microcracks have also been observed near the intersection between impinging twin bands and grain boundaries in the intermetallic TiAl (Bieler et al. 2005). The weaker boundaries have been identified as the ones in which the neighbouring grains cannot accommodate the twin strain. Cracks have been observed near the intersection between twins and grain boundaries in several BCC metals at low temperatures (Cr, W, Mo) (Marcinkowski and Lipsitt 1962;Gilbert et al. 1964). The grain size dependence of the stress to nucleate twins and cracks are similar. However, it is not clear if a twin induces crack nucleation or vice-versa. The same properties are observed in α-uranium, which is the material of interest in the present paper (Taplin 1964;Taplin and Martin 1965). α-uranium exhibits intergranular fracture in the temperature range from − 100 • C to 200 • C (Collins and Taplin 1978), while ductile fracture (Taplin and Cocks 1967) and inclusion cracking are observed at higher temperatures (Davies and Martin 1961).
Imaging techniques such as differential aperture X-ray microscopy (Balogh et al. 2013) and highresolution electron backscatter diffraction (EBSD) (Abdolvand and Wilkinson 2016) have been used to measure the elastic strain field near the twin tip and at the twin interface (Hubbell and Seltzer 1995). However, the large X-ray absorption coefficient of αuranium limits the applicability of the first technique (Hubbell and Seltzer 1995), the preparation of strainfree α-uranium surfaces for EBSD is challenging (Sutcliffe et al. 2019;Earp et al. 2018) and X-ray transmission integrates the signal over the thickness of the sample (Irastorza-Landa et al. 2016, 2017b, limiting the ability to image individual twins. In-situ techniques, which can reveal the chronological order of twin and crack nucleation and propagation, are needed but the characteristic time during which microcracks and twins nucleate is very short, as they propagate close to the speed of sound in the material (Oberson and Ankem 2005). A unique example in the literature is the measurement of the magnetic properties with microsecond time resolution, showing that twinning takes place a few microseconds earlier than crack nucleation in ironsilicon (Williams and Reid 1971).
Several open questions remain. What is the effect of strain localization due to twinning on fracture? Which combinations of grain orientations favour intergranular fracture? How does the crack propagate once nucleated? What is the chronological order of events: do propagating cracks nucleate twins or vice-versa? What is the role of dislocations? These questions are particularly important because twinning can be inhibited using prestrain at a higher temperature (Boucher and Christian 1972) or by introducing precipitates (Chun et al. 1969).
Crystal plasticity finite element (CPFE) and fracture mechanics simulations are suitable to investigate these issues. CPFE is a computational finite strain method that takes into account the plastic deformation of slip and twin systems Grilli 2016), and the reorientation of the crystal lattice due to twinning (Grilli et al. 2020d). The different orientations of neighbouring grains can be included and elastic strain incompatibility can be calculated, as well as the intergranular stress (Petkov et al. 2019).
Intergranular fracture can be described using a cohesive-zone model (Simonovski and Cizelj 2011). Elements of negligible thickness are used to describe the separation of the grain boundary interface. Traction-separation relationships for interfaces were first introduced by Barenblatt (Barenblatt 1962). Softening of the interface element takes place after a critical value of the stress is reached; this allows both crack nucleation and propagation to be simulated. Camacho and Ortiz showed that cohesive-zone modelling is compatible with finite element modelling of the bulk material and is able to reproduce arbitrary crack paths between 2D bulk elements in a mesh (Camacho and Ortiz 1996). Ortiz and Pandolfi extended the method to 3D dynamic simulations (Ortiz and Pandolfi 1999).
The recent development of continuum models to describe discrete twins ) allows us to simulate the concentration of plastic deformation in a narrow band. Twin nucleation and growth has been reproduced (Qiao et al. 2016). The coupling with phase field fracture has allowed twin nucleation induced by crack propagation from a notch to be studied (Clayton and Knap 2016;Grilli et al. 2018a). However, no models for twin-induced fracture are available.
It is fundamental that such a model can correctly predict the twin thickness and the stress for twin nucleation and growth. This is because the stress at the twin tip depends on the twin thickness (Sauzay and Moussa 2013). This has been achieved using a non-local model for the critical resolved shear stress (CRSS) for twinning, which has been validated using in-situ EBSD experiments (Grilli et al. 2020d).
In this paper, a coupled plasticity-fracture model is developed and applied to 3D polycrystal simulations. The model is calibrated using data in the literature from experiments on α-uranium . The effect of discrete twins on crack nucleation and propagation is investigated. Specific combinations of the orientations of neighbouring grains that favour or prevent intergranular fracture are investigated.
Section 2 describes the cohesive-zone model and the formulation of 3D interface elements. Section 3 describes the crystal plasticity model used for the mechanical behaviour of the grains and the coupling with a continuum model for twinning. In Sect. 4, polycrystal simulations are carried out to calibrate the parameters of the cohesive zone model. The details of crack nucleation and propagation at the grain boundaries are investigated in Sect. 5. In Sect. 6, the crack nucleation and propagation is studied in the case of noncolumnar grains. Sections 7 and 8 contain discussion and conclusions.

Cohesive zone model and interface elements
Zero thickness interface elements are assigned to grain boundaries. An algorithm has been developed to automatically generate such interface elements in an arbitrary polycrystal (Grilli 2020). First, the mesh is generated with 3D hexahedral elements and zero thickness grain boundaries. Hexahedral elements are preferred to linear tetrahedral elements in crystal plasticity finite element simulations because they do not exhibit volumetric locking and a consequent stiff response during bending (Cheng et al. 2016). Then, nodes at the grain boundaries, which are shared by two or three grains, are duplicated or triplicated respectively. An example is shown in Fig. 1a. The original nodes 1, 2, 3, 4 belong to a grain boundary between grain 1 and grain 2. During the duplication, nodes 5, 6, 7, 8 are created, which have the same undeformed coordinates as nodes 1, 2, 3, 4 respectively. Grain 1 retains the original nodes 1, 2, 3, 4 while nodes 5, 6, 7, 8 are assigned to grain 2. Therefore, nodes 1, 2, 3, 4 belong to one hexahedral element in grain 1, while nodes 5, 6, 7, 8 belong to an element in grain 2.
An interface element is created that includes all the nodes from 1 to 8. This interface element connects the two hexahedral elements in grain 1 and 2. A bilinear traction-separation law, as shown in Fig. 1b, is used for the interface element to represent the cohesion forces of the grain boundaries (Barenblatt 1962). Because of the arbitrary orientation of each interface element in space, it is necessary to define a local reference system in which the separation vector Δ is defined. This will be called the mid-plane reference system and it is defined as follows.x m is the unit vector that connects the midpoint m1 between node 1 and node 5 with the mid-point m2 between node 2 and node 6. The unit vector normal to the mid-planeẑ m is found by taking the cross product ofx m and the direction that connects the mid-points m1 and m4. The third unit vectorŷ m is calculated as: The three orthogonal unit vectors x m ,ŷ m ,ẑ m constitute the mid-plane reference system; its origin is in the middle of the quadrilateral with vertices m1, m2, m3, m4. This choice allows the traction-separation law in this coordinate system to be defined, which is independent of the orientation of the interface element. In the following, the separation vector Δ = (Δ s1 , Δ s2 , Δ n ) is defined in the mid-plane reference system and can be calculated on the x m ,ŷ m plane using shape functions, as reported in Appendix A. The cohesive force vector T = (T s1 , T s2 , T n ) is also expressed in the mid-plane reference frame. The traction/separation law is given by: where 0 < D < 1 is the damage variable, describing the degradation of the normal stiffness K n . Δ 0 is the normal separation at which damage starts and Δ f is the normal separation at which damage is complete, as shown in Fig. 1b. Given a damage D, Δ p (D) is the normal separation at which the maximum normal traction is reached. If Δ n > Δ p (D), damage increases and normal traction decreases. This formulation includes contact (Δ n < 0) during which damage does not affect the normal cohesive force. In the present model, shear separation does not induce damage, but shear cohesive force is degraded by existing damage: where Δ s = Δ 2 s1 + Δ 2 s2 is the magnitude of the shear separation. If D = 1, there is no shear traction, even in the case of contact (Δ n < 0), i.e. the contact is assumed to be frictionless. The components of the shear traction along the axesx m andŷ m are: The damage nucleation induced by shear separation is relevant in brittle materials Ravichandran 1996, 1997) but has been neglected in the description of intergranular fracture in metals (Yamakov et al. 2006;Simonovski and Cizelj 2011).
It is of the utmost importance that the nodes of the cohesive elements are arranged as in Fig. 1a. If the positions of nodes 2 and 4 are exchanged, and the same for nodes 6 and 8, the unit vectorẑ m would point in the opposite direction. In that case, according to the definition of the separation vector in Appendix A, interface opening (Δ n > 0) and contact (Δ n < 0) would be inverted and not reproduced correctly. The algorithm that generates the interface elements takes this issue into account (Grilli 2021).
The principle of virtual work provides the weak form of the equilibrium equations (Park and Paulino 2012): (6) where is the representative volume, c is the cohesive interface and is the boundary. The δ symbol represents the variation of a variable. ε is the strain tensor and σ is the Cauchy stress. u G is the displacement vector, expressed in the global reference frame, and T ext is the external force. The second term on the left-hand side of Eq. (6) is introduced in the commercial finite element code Abaqus by using the UEL subroutine (Smith 2009). This requires us to define the contribution to the internal force vector f e coh for each cohesive element. Therefore, the variation δΔ has to be expressed as a function of the variation of the displacement vector δu G . This is done by using the B matrix defined in Eq. (9. 11) of Appendix A: where δU = (δu 1 , δu 2 , δu 3 , δu 4 , δu 5 , δu 6 , δu 7 , δu 8 ) T is a [24×1] vector that contains the variation of the displacement components of the eight nodes of the cohesive element. Since δU is expressed in the mid-plane reference frame, a transformation matrix M is needed. M is a [24 × 24] block diagonal matrix. It has eight equal [3 × 3] main-diagonal blocks. These blocks contain the [3 × 3] rotation matrix that transforms vectors from the global reference frame to the mid-plane reference frame. Therefore: where δU G is the vector δU expressed in the global reference frame. Therefore, the contribution of a single cohesive element e to the second term on the left-hand side of Eq. (6) is given by: where f e coh is the internal [24 × 1] force vector of the cohesive interface. The vector δU G can be taken outside of the integral because it depends only on the displacements at the nodes. f e coh is calculated at 4 Gauss points, as described in Appendix A.
The Jacobian ∂ f e coh /∂U G which is calculated in the UEL subroutine is reported in Appendix B. The cohesive zone model parameters used in the following are reported in Table 1. The value of σ max is calibrated in Sect. 4.

Crystal plasticity model
The behaviour of the grains is described using a crystal plasticity model coupled with a continuum model for discrete twins (Grilli et al. 2020b). A variable ϕ is used to represent discrete twins. It can increase from 0 (untwinned region) to 1 (fully twinned region). The model is introduced in Abaqus by using a UMAT subroutine (Smith 2009). An implicit CPFE framework is used to calculate the increment of the Cauchy stress σ at each time step (Dunne et al. 2007). The rate of change of the corotational stress tensor is calculated using Hooke's law (Dunne and Petrinic 2006;Sakano et al. 2020): where C is the fourth order elasticity tensor (Fisher and McSkimin 1958). ε e , ε p and ε are the elastic, plastic and total strain. C is interpolated between the untwinned crystal lattice and the twinned crystal lattice when ϕ grows from 0 to 1. Plastic deformation is described using the 8 most active slip systems of α-uranium at room temperature (McCabe et al. 2010). The twin system [310] (130) is used in the following simulations (Zhou et al. 2016).
The plastic strain rate is calculated by summing the contributions of the slip ratesγ α (σ ) on the slip systems and the contribution of the twinning rate (Kalidindi 1998): where s α and n α are the slip direction and normal of the slip system α. s β and n β are the twin direction and normal, and γ twin β is the total shear produced by twinning. Given the total strain rateε in one time step, Eqs. (10)-(11) represent a non-linear system of equations in which the unknown is the Cauchy stress increment. They are solved using a Newton-Raphson algorithm (Dunne and Petrinic 2006).
The slip ratesγ α (σ ) are calculated using a powerlaw relationship (Grilli et al. 2015;Irastorza-Landa et al. 2017a): whereγ 0 and n are constants. τ α is the resolved shear stress (RSS) on slip system α and τ c α is the CRSS (Roters 2011; Jafari et al. 2017). This rate-dependent law allows us to find a unique decomposition of the plastic strain increment into slip increments (Zamiri and Pourboghrat 2010). τ c α determines the hardening of the slip system and it is evolved in time using a dislocation density-based model, whose details are reported in (Beyerlein and Tomé 2008;McCabe et al. 2010;Grilli and Cocks 2019;Grilli et al. 2020a). The evolution of the dislocation densities is based on multiplicationannihilation rate equations (Kocks and Mecking 2003;Grilli et al. 2018b). This hardening model has been able to reproduce the elastic lattice strain in polycrystals measured using neutron diffraction (Grilli et al. 2020a) and the strain field inside individual grains measured using digital image correlation (Grilli et al. 2020c).
The time evolution of ϕ determines the nucleation and growth of discrete twins. It has two contributions: the first term is a power-law relationship similar to the one used for slip (McCabe et al. 2010): This represents stress-induced twin nucleation, which takes place only when the RSS on the twin system is positive. The second term in Eq. (13) describes the driving force moving the atoms towards their equilibrium position in the twinned crystal lattice, after they have passed the energy barrier for twinning (Liu et al. 2019): (15) ϕ = 0.5 corresponds to the maximum of the energy barrier and 1/ f is a characteristic time interval during which the twin process completes.
The CRSS for twinning τ c β (ϕ) has two contributions: The first term τ 0 (ϕ) models the twin nucleation and stress relaxation during twin propagation, as observed experimentally (Lynch et al. 2014), using a bilinear law: where τ 0 is a constant. τ twin (ϕ) models the interaction between the twin interface and mobile dislocations (Ojha et al. 2014). When a mobile dislocation intersects a pre-existing twin boundary, it is decomposed into a sessile residual dislocation and a twinning dislocation. New twin layers are formed by the motion of these twinning dislocations, i.e. twin growth. The accumulation of residual dislocations prevents further mobile dislocations from interacting with the twin interface, leading to a higher CRSS for twin growth. The specific interaction depends on the dislocations character and on the direction from which they approach the twin interface . In this paper, we model this process with a non-local term that quantifies the total thickness of the twinned region in the neighbourhood of a point, which is also proportional to the density of residual dislocations: where τ 0 twin is a constant. In order to quantify the thickness of the twinned region, the integration region c is chosen as a cylinder with height l 0 and radius r 0 , Only integration points where the twin has overcome the energy barrier (ϕ > 1/2) are considered when evaluating the integral.
As shown in Fig. 2, given an axis X perpendicular to the twin planes, the term τ twin (ϕ) can be evaluated. If the point considered moves towards the right and the cylindrical volume intersects the second twin on the right, the term τ twin (ϕ) increases. This prevents the further growth of the twin on the left. The smooth increase of the term τ twin (ϕ) in the presence of the pre-existing twins leads to the formation of smooth twin interfaces, in which ϕ varies from 1 to 0 along a certain length scale. This prevents the formation of sharp interfaces between twinned and untwinned regions that could not be captured properly by a finite element model.
It has been shown that this formulation can reproduce the nucleation and growth of discrete twins at the correct value of stress and with thickness as observed in EBSD experiments (Grilli et al. 2020d).
α-uranium (orthorhombic) material parameters are used for the plasticity model in the following simulations (Grilli et al. 2020a). They are reported in Table  2. The crystal plasticity code is available in the following repository: https://github.com/TarletonGroup/ CrystalPlasticity (Tarleton 2020).

Polycrystal simulation and model calibration
In this section polycrystal simulations are carried out; the maximum stress and strain to failure are compared with experiments in the literature. A calibration of the cohesive zone model parameters is pro-  posed. The columnar grain structure is shown in Fig. 3a and grains are numbered from 1 to 6. The representative volume is a parallelepiped with size 60 µm × 60 µm × 0.5 µm. The mesh, shown in Fig. 3b, is constituted of hexahedral elements with average size 0.5 µm. Cohesive elements are placed along the grain boundaries, including triple junctions. Boundary conditions, modelling pure tension, are shown in Fig. 3a. u z = 0 is imposed on the surface z = 0, u y = 0 is imposed on the surface y = 0, u x = 0 is imposed on the surface x = 0. A displacement u y is imposed on the surface y = 60 µm. This displacement increases linearly with time, from 0 to 18 µm, inducing an average strain up to 30%. The applied macroscopic strain rate used in all the following simulations is 0.001 s −1 .
The orientation of the crystal lattice of the i-th grain in Fig. 3a is determined by rotation matrices R i . These rotation matrices transform a vector from the crystal lattice reference frame into the global reference frame, shown in Fig. 3a. The rotation matrix for grain number 1 is: This matrix aligns the direction of the twin system [310] (130) at 45 • in the (x,y) plane and so the Schmid factor of this twin system in grain number 1 becomes 0.5. The rotation matrix for the i-th grain is found by multiplying R 1 by a rotation around the z axis: The angle θ i for each grain is reported in Table 3. The angle θ quantifies the magnitude of the misorientation between the six grains. The value of the separation to failure Δ f is fixed at 1 µm. Simulations are carried out with different values of σ max , the damage initiation stress shown in Fig.  1b. The stress-strain curves, obtained by averaging the stress component σ yy on the load surface y = 60 µm in Fig. 3a, are shown in Fig. 4.
The maximum stress and strain to failure are strongly dependent on σ max , as shown in Fig. 4a. Tensile stresses up to 310 MPa have been measured in tensile bar experiments on coarse-grained α-uranium without reaching failure (Grilli et al. 2020a, c, d). This represents a lower limit for the calibration. As shown in Fig. 4a, this lower bound is satisfied if σ max ≥ 500 MPa. Additionally, Huddart and Harding measured a strain to failure around ε f = 20% for cast α-uranium at low strain rate (Huddart et al. 1980). This condition is approximately satisfied in the present simulations if σ max = 650 MPa. Values of σ max larger than 700 MPa leads to an unrealistically high strain to failure. Therefore, the value σ max = 650 MPa will be used in all the following simulations, which leads to a fracture energy of 1 2 σ max Δ f = 325 J/m 2 and opening of Δ 0 = 2.87 nm at damage initiation. Figure 4b shows the stress-strain curves for different values of the misorientation angle θ . The lowest misorientation θ = 5 • shows the maximum stress and the maximum strain to failure. However, the relationship between strain to failure and misorientation angle θ does not monotonically decrease. The strain to failure for θ = 45 • is slightly larger than for θ = 20 • . The dependence on the misorientation angle confirms that the previous calibration of σ max is not strongly dependent on the choice of the texture and, therefore, it is reliable. Fig. 5 shows the twin variable ϕ and fracture path at 23% average strain for three different misorientation angles θ . The discrete twins (ϕ = 1) and their misorientations between neighbouring grains can be observed. The fracture path is similar for the three misorienta- tions. The boundaries between grains 1-4 and between grains 2-3 are the first to completely damage. These are the boundaries with an orientation almost perpendicular to the load direction. They also have several discrete twins impinging on them from grains 1, 3, 4. The only difference in the fracture path is visible in Fig. 5a: a crack initiates at the grain boundary between grains 1-6 at the triple junction between grains 1, 5, 6. A discrete twin impinges on the triple junction from grain 5.

Crack initiation and propagation
In order to investigate the effect of twin tips impinging on a grain boundary, a three grains simulation is carried out, as shown in Fig. 6. The two grain boundaries have the same orientation with respect to the load direction; therefore, the effect of the twin orientation can be decoupled from the effect of the grain boundary orientation with respect to the load direction.
The representative volume is a parallelepiped with size 50 µm × 80 µm × 0.5 µm. The mesh, shown in Fig. 6b, is constituted of hexahedral elements with average size 0.5 µm. Three crystals are present, as shown in Fig. 6a, and cohesive elements are placed at the two grain boundaries.
Boundary conditions are shown in Fig. 6a. u y = 0 is imposed on the surface y = 0 and a displacement u y is imposed on the surface y = 80 µm. This displacement increases linearly with time. u x = 0 is imposed at the point (25, 0, 0) µm and u z = 0 is imposed at the points (25, 0, 0) µm, (0, 0, 0) µm, (50, 0, 0) µm. These boundary conditions allow the grains to rotate during tensile loading. The effect of this rotation can be studied in these simulations.
The orientation of the crystal lattice of the grains in Fig. 6a is determined by the following rotation matrices for texture 1 (T1): The twin variable ϕ and fracture path are shown in Fig. 7 for the two different textures at 12.5% average strain. The largest damage is induced at the grain boundary between a grain with twins and a grain without twins. The propagation of the crack at the grain boundary between the two grains with twins is slower.
In order to better investigate the dynamics of the crack propagation, the damage variable D and the twin variable ϕ are reported along the upper grain boundary (between grain 1 and 2) for texture 1. These two quantities are taken along a node path, as shown in Fig. 6b, and reported in the same plot at different values of the strain, as shown in Fig. 8. The discrete nature of the twinning process results in a non-uniform distribution of stress along the grain boundaries, with the stress normal to a boundary being greatest adjacent to  Fig. 6b for texture 1 a twin. Crack nucleation (Fig. 8a) takes place on the right side of a twin tip due to the tensile stress caused by twinning. Cracks nucleate at two points adjacent to two neighbouring twins independently. As the macroscopic strain is increased, and twinning intensifies, intergranular damage develops to the left of one of the twins, promoting the development of damage at the twin, leading to the formation of a crack that propagates towards the next twin (labelled left crack tip in Fig. 8b). Other cracks nucleate ahead of this crack tip as other twins reach the grain boundary (Fig. 8c), which merge with the propagating crack. The motion of the left crack tip slows down when it meets a nucleating twin (Fig. 8d). The interaction with a twin nucleus also speeds up the nucleation and growth of the twin because of the stress concentration induced.  Figure 9a reports also the crack tip positions for texture 2.
As shown in Fig. 9a, the crack tip position does not increase uniformly with time. As discussed previously, the crack at the upper grain boundary in texture 1 tends to propagate towards the left. Therefore, Fig. 9a shows a larger motion of that left crack tip. The right crack tip in texture 1 advances at the beginning, but then its position remains almost stationary for the rest of the simulation. This is due to plastic relaxation due to the presence of twins on both sides of the grain boundaries and to the rotation of grains 1 and 2. This rotation accommodates the tensile load on the top surface without further propagation of the crack towards the right.
The right crack tip in texture 2 propagates more uniformly, as shown in Fig. 9a. This is because this crack tip does not interact with any twin during its motion, as shown in Fig. 7b. Figure 9b shows the left crack tip velocity in texture 1, i.e. the time derivative of the left crack tip position. The arrows (a), (b), (c), (d) correspond to the values of the average strain in Fig. 8. Velocity peaks are present when the crack tip propagates between two twin tips, while the velocity has minima whenever the crack tip meets a twin tip or a nucleating twin.

Bicrystal simulation with non-columnar grains
All the previous simulations were carried out using columnar grains. A bicrystal representative volume with two non-columnar grains is investigated to understand the effect of discrete twins on intergranular fracture for the case in which the grain boundary surface normal is not coplanar with both the twin direction and twin plane normal.
The representative volume and mesh are shown in Fig. 10a. Figure 10b shows only the edges of the two grains and clarifies the orientation of the grain boundary. The size of the representative volume is 35 µm Fig. 11 Bicrystal simulation with non-columnar grains: twin variable ϕ and damaged plane (blue). Only larger twins with elements in which ϕ > 0.4 are shown. The crack is constituted of elements with D > 0.9. a 18%, b 23%, c 30%, d 38% strain × 25 µm × 10 µm. The mesh consists of hexahedral elements with average size 0.6 µm.
Boundary conditions for pure tension are shown in Fig. 10b. u z = 0 is imposed on the surface z = 0, u y = 0 is imposed on the surface y = 0, u x = 0 is imposed on the surface x = 0. A displacement u y is imposed on the surface y = 25 µm. This displacement increases linearly with time.
The orientation of the crystal lattice of grains 1 and 2 in Fig. 10a is: Since the grain boundary is almost perpendicular to the load direction, the maximum stress of the interface elements has been increased to σ max = 800 MPa in order to observe the complete development of discrete twins before failure.
The time evolution of the twin variable ϕ is shown in Fig. 11. The main discrete twins, consisting of at least 50 elements with ϕ > 0.4, are shown and numbered. The blue region represents the crack, i.e. the surface of the grain boundary on which D > 0.9.
Twins 1 and 2 are the largest discrete twins, which nucleate in the central part of the geometry on different planes. They therefore impinge on the grain boundary at two different locations. The tensile stress induced by twins 1 and 2 on the grain boundary is maximum between the two twin tips, and a crack first nucleates here (Fig. 11a). The crack nucleates near the surface z = 0. This is due to the orientation of the grain bound-ary: the tip of twin 1, which is the most important in determining crack nucleation, is closer to the grain boundary on the surface z = 0 than on the surface z = 10 µm. The fracture plane propagates first along the direction parallel to the twin plane, then along the grain boundary (Fig. 11b). The damage surface propagates much faster towards the right (Fig. 11c). This is because twin 2 is already completely formed when the crack propagates near its tip, while the nucleation of twin 5 at the crack tip slows down its propagation (Fig. 11c). The stress is not high enough to propagate twin 4, because the grain boundary above it is completely damaged. The stress concentration induced by the left crack tip leads to the nucleation and propagation of twin 5 before the grain boundary is completely damaged (Fig. 11d)).

Discussion
The model calibration carried out in Sect. 4 shows that the knowledge of a lower bound on the maximum stress and the strain to failure is sufficient to precisely calibrate the parameter σ max of the cohesive zone model. This is due to the strong dependence of the strain to failure on σ max in the coupled discrete twin-fracture model. This type of calibration has been used for mode I tearing of plates (Woelke et al. 2015). Models including void nucleation and growth are also calibrated using the experimental strain to failure (Perzyna 1984).
Since intergranular fracture depends mostly on the stress normal to the grain boundary interface (Simonovski and Cizelj 2011), the grain boundaries that are predicted to fail earlier are the ones perpendicular to the load direction.
The texture plays an important role on crack nucleation and propagation, as shown in Sect. 5. If two neighbouring grains have similar orientations, the stress concentration at the grain boundary induced by twin tips in one grain is relaxed by twins in the neighbouring grains. If a grain boundary has one neighbouring grain that twins and another neighbouring grain that does not twin, stress relaxation cannot take place. In this case, failure of the grain boundary is faster. This behaviour is similar to brittle materials with highly oriented microstructures, in which the fracture toughness approaches that of a single crystal (Schultz et al. 1994). It is also consistent with the experimental results on TiAl (Bieler et al. 2005), in which the weaker grain boundaries are the ones in which the neighbouring grain cannot accommodate the twin strain.
The present simulations reveal an additional effect of the texture. As shown in Figs. 7b and 11, a propagating crack can prevent the growth of twins in grains that are favourably oriented for twinning. This happens if the deformation can be accommodated by other mechanisms, such as grain rotation in Fig. 7b, or if the grain boundary is almost perpendicular to the load direction (twin 4 in Fig. 11).
The present simulations clarify the chronological order of twin and crack nucleation (Mahajan and Williams 1973). Twins nucleate first. Once a twin impinges on a grain boundary, the tensile stress that is present on one side of the twin tip is sufficient to nucleate a crack. The crack first propagates along a direction parallel to the twin plane until the intersection between the twin tip and the twin plane is completely damaged. Then, the crack propagates along the grain boundary, as shown in Sect. 6. If multiple twins impinge on the same grain boundary, multiple microcracks can nucleate in correspondence to the different twin tips and then merge, as shown in Sect. 5. The propagating crack can induce stress concentrations that are sufficient to nucleate a twin starting from the grain boundary. Patterns of twins extending from the grain boundaries towards the centre of the grains are commonly observed (Gussev et al. 2018). The stress relaxation induced by the nucleation of these twins slows down the propagation of the crack tip.
In the crystal plasticity model used, dislocation slip is present in both the twinned and untwinned regions. The simulations show crack nucleation near twin tips and there is no strong correlation with slip, which is a more uniform plastic deformation mode in the simulated grains. The stored dislocation energy density has been identified as an important factor to explain the locations of the observed twins (Paramatmuni et al. 2020), but in the present simulations the dislocation density does not affect the critical stress for twin nucleation. Including slip bands in the present model would be necessary to describe metals in which twin-induced microcracks are observed, while the actual failure is slip-induced (Reid 1981).
Twin boundary cracking (twin parting) is another observed fracture mechanism that we have not specifically investigated in this paper (Cahn 1953;McCabe et al. 2008).

Conclusions
In this paper, we have coupled a discrete twin model and a cohesive zone model to understand the effect of the stress induced by twins on intergranular fracture. The dynamics of twin-induced intergranular fracture is clarified. Microcracks nucleate near the intersection between a twin tip and a grain boundary, in the area where the twin tip induces a high tensile stress. The subsequent crack propagation can nucleate more twins, which then propagate away from the grain boundary. This mechanisms slows down crack propagation. Therefore, a grain boundary, surrounded by two grains that are favourably oriented for twinning, is more resistant against fracture.
An algorithm that can generate the interface elements along grain boundaries of arbitrary polycrystals has been developed. This model can be used to predict the strain to failure and stress-strain curve of cast metals, given the specific orientation of the grains.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix A: Shape functions of the cohesive element
The separation vector Δ = (Δ s1 , Δ s2 , Δ n ), expressed in the mid-plane reference frame, is a function of the coordinates on the x m ,ŷ m plane. Given the deformed coordinates of the nodes x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 , x 8 of the cohesive element in Fig. 1, the separation vectors at the mid-points m1, m2, m3, m4 are given by: These separation vectors can be expressed also as a function of the displacement vectors of the nodes u 1 , u 2 , u 3 , u 4 , u 5 , u 6 , u 7 , u 8 because nodes in the pairs (1, 5), (2, 6), (3, 7), (4, 8) are coincident in the undeformed configuration.
Four shape functions are used to interpolate the separation vector on the x m ,ŷ m plane (Cook 1995): where −1 < ξ < 1, −1 < η < 1 are the coordinate of a reference square element. Therefore, the separation vector can be written as: (9. 9) Eqs. (9. 1)-(9. 4) and (9. 9) are linear relationships between the separation vector and the displacement vectors of the nodes. Therefore, it is possible to build a [24 × 1] vector: U = (u 1 , u 2 , u 3 , u 4 , u 5 , u 6 , u 7 , u 8 ) T (9. 10) which contains the three displacement components of the eight nodes, and to define a [3 × 24] matrix B such that (Alfano et al. 2015): The B matrix is a function of the coordinates (ξ, η). Four Gauss points are used for the calculation of the internal force vector in Eq. (9). Their coordinates in the reference element are reported in Table 4.

Appendix B: Jacobian calculation
The Jacobian contribution of one cohesive element e is found by calculating the derivative of f e coh in Eq. (9) with respect to the nodal displacement vector U G : The derivative of the rotation matrices in M are neglected to increase the computational efficiency. This Jacobian leads to good convergence if no rigid body rotations are applied to the representative volume.
Using the relationship between the separation vector Δ and the nodal displacement vector U G : which is a non-symmetric [24 × 24] matrix due to the lack of symmetry in ∂ T /∂Δ. Therefore the unsymmetric solver option is specified in the Abaqus input file (Smith 2009). The integral is evaluated using the four Gauss points of the cohesive element. The Jacobian in the mid-plane reference frame ∂ T /∂Δ is a [3 × 3] matrix that is obtained from the derivatives of the traction vector components in Eqs.
(2)-(5): (10. 5) The derivatives along the two shear components can be found as: (10. 9) The off-diagonal derivatives, i.e. derivatives of the normal traction with respect to shear separation and viceversa, are not all zero because an increase of normal displacement leads to an increase in the damage and consequent decrease of the shear traction: If Δ n is increasing and Δ n > Δ p , then damage increases, otherwise the damage variable is constant: (10. 14) Therefore, the derivative of the shear traction with respect to the normal displacement are found from Eq. : (10. 15)