Modelling the nucleation and propagation of cracks at twin boundaries

Fracture arising from cracks nucleating and propagating along twin boundaries is commonly observed in metals that exhibit twinning as a plastic deformation mechanism. This phenomenon affects the failure of macroscopic mechanical components, but it is not fully understood. We present simulations in which a continuum model for discrete twins and a cohesive zone model are coupled to aid the understanding of fracture at twin boundaries. The interaction between different twin systems is modelled using a local term that depends on the continuum twin variables. Simulations reveal that the resolved shear stress necessary for an incident twin to propagate through a barrier twin can be up to eight times the resolved shear stress for twin nucleation. Interface elements are used at the interfaces between all bulk elements to simulate arbitrary intragranular cracks. An algorithm to detect twin interfaces is developed and their strength has been calibrated to give good agreement with the experimentally observed fracture path. The elasto-plastic deformation induced by discrete twins is modelled using the crystal plasticity finite element method and the stress induced by twin tips is captured. The tensile stress caused by the tip of an incident twin on a barrier twin is sufficient to nucleate a crack. A typical staircase fracture path, with cracks propagating along the twin interfaces, is reproduced only if the strength of the twin interfaces is decreased to about one-third of the strength of the bulk material. This model can be used to help understand fracture caused by the activation of multiple twin systems in different materials.

by discrete twins is modelled using the crystal plasticity finite element method and the stress induced by twin tips is captured. The tensile stress caused by the tip of an incident twin on a barrier twin is sufficient to nucleate a crack. A typical staircase fracture path, with cracks propagating along the twin interfaces, is reproduced only if the strength of the twin interfaces is decreased to about one-third of the strength of the bulk material. This model can be used to help understand fracture caused by the activation of multiple twin systems in different materials.

Introduction
Deformation-induced twinning is an important mechanism, common in bcc, hcp and low symmetry crystals (Christian 2002). It is observed when other plastic deformation mechanisms become difficult, e.g. at low temperature or high strain rate (Mahajan and Williams 1973;Christian and Mahajan 1995;Kapoor et al. 2015;Inouye and Schaffhauser 1969). Twin boundaries are crystallographic planes separating two regions: the parent crystal and the twinned crystal, which are mirror images of each other by reflection (Hull and Bacon 2011). Twinning can increase the ductility of metals and alloys. Twin migration under stress also contributes to the accommodation of plastic strain ).
At the same time, twin boundaries act as barriers for dislocation motion (Wang et al. 2020); therefore, they increase strength (Lu et al. 2004;Anderoglu et al. 2008).
There has been a growing interest in twinning because of the recent experimental demonstration that a high density of twin boundaries can be induced in Al alloys by deposition onto a single crystal substrate . Higher twinning activity has been revealed in additively manufactured steel (Yang et al. 2021;Grilli et al. 2021a) because of the nitrogen atmosphere used in the process (Pham et al. 2017). Moreover, the density of twin boundaries can be controlled by changing the scanning speed (Gao et al. 2020) and laser energy input (Do and Li 2016;Hocine et al. 2020) during selective laser melting.
To reach the twinned configuration, atoms rearrange and produce a homogeneous shear inside a twin band, whose characteristic thickness can range from nanometres to hundreds of micrometres (Zhang et al. 2009;Beyerlein and Tomé 2008;Ardeljan et al. 2017;Salem et al. 2003). If the twin intersects an obstacle, such as a grain boundary or another twin interface (barrier twin), a stress concentration is created (Sauzay and Moussa 2013;Koko et al. 2021). In some crystals, the shear deformation of the incident twin can transmit through the barrier twin (Sun et al. 1993). In this case, a secondary twin, deflected with respect to the incident twin, is created inside the barrier twin (Wardle et al. 1993). In other cases, the incident twin has a lenticular shape, whose edge terminates near the barrier twin, without penetrating it (Abdolvand and Wilkinson 2016). Twin penetration requires the following criteria: direction and magnitude of the plastic shear of the incident and secondary twin must be identical; their traces on the barrier twin plane must be parallel (Cahn 1953). Moreover, energy is required to reorient the crystal lattice inside the barrier twin.
Twin intersections have been investigated as a possible site for fracture due to the local stress concentration (Müllner et al. 1994). In many materials, fracture surfaces inside crystals intersect the twin planes. This has been observed in uranium (Field et al. 2009), copper (Boettner et al. 1964), Cu-Al (Qu et al. 2008), Cu-Zn (Zhang et al. 2012), ferritic TWIP steel (Koyama et al. 2013), austenitic stainless steel (Heinz and Neumann 1990). This phenomenon is also called "twin parting" by mineralogist. The intragranular fracture path in ferritic TWIP steel has a staircase shape: part of the crack surface lies on the twin plane and the remaining crack surface connects neighbouring twins (Koyama et al. 2013). The resulting macroscopic crack is approximately perpendicular to the load direction. It has been observed during both monotonic and cyclic load. The intersection between shear bands and twins can also be responsible for crack nucleation (Xiong et al. 2017). One hypothesis is that the energy release rate G c (Grilli and Koslowski 2019) of the interface between the twinned and untwinned regions is lower than in the bulk. However, no direct measurements of this fracture energy are available.
Other factors play an important role, such as dislocation slip, the shape of the twin and the stacking fault energy (SFE) (Zhang et al. 2012). The activation of slip planes parallel to the twin boundary and the presence of twin boundary steps induce crack formation at the twin boundary (Boettner et al. 1964). A lower SFE, obtained by alloying (Murr 1975), favours twin boundary fracture (Qu et al. 2008).
The following questions remain unanswered. What is the interaction strength between discrete twins that is necessary to prevent twin intersections? How strong is the twin interface in order to explain the intracrystalline fracture path observed? What is the effect of twin orientation on the fracture behaviour?
In this paper we use crystal plasticity finite element (CPFE) simulations (Grilli 2016;Irastorza-Landa et al. 2017b), a continuum model for discrete twins (Grilli et al. 2020b) and a cohesive zone model for intragranular fracture (Grilli et al. 2021b) to answer these questions. The material of interest is α-uranium, which is the stable phase of uranium at room temperature (Cahn 1951). The CPFE method includes the plastic deformation due to all slip and twin systems Irastorza-Landa et al. 2017a). Recently, continuum models have been developed to describe the nucleation and growth of discrete twins with a specific thickness (Grilli et al. 2020b). A twin variable ϕ β varies from 0 (untwinned region) to 1 (twinned region) when a twin of type β forms (Liu et al. , 2019. The model includes the reorientation of the crystal lattice due to twinning. Intragranular fracture is described by a cohesive zone model using a bilinear tractionseparation law (Barenblatt 1962;Camacho and Ortiz 1996;Alfano et al. 2015). This type of model has been previously applied to both brittle (Yang et al. 2019) and ductile fracture (Tvergaard 2001), and even compared directly with Gurson-type models (Siegmund and Brocks 1998b, a). Therefore, it is a suitable model to describe the brittle and semi-ductile crack growth which is characteristic of cast α-uranium at room temperature (Collins and Taplin 1978). In the present work, cohesive elements with zero thickness are used (Segurado and LLorca 2004;Segurado and Llorca 2002;Nguyen 2014b) at the interface between all pairs of bulk elements in the mesh, which describe the elastoplastic behaviour. The chosen value of the stiffness of the interface elements leads to good simulation convergence and it is higher than the stiffness of the bulk elements, thus reducing their artificial compliance.
The interaction between different twin families is described by a local term that depends on the twin variable ϕ β . In previous homogenised models that describe the hardening of a twin system based on the twin volume fraction of a second twin system (Roters et al. 2010;) discrete twins were not considered (Kalidindi 2001). We show that by changing the interaction coefficient, a transition from a model where twin intersections are possible to a model where incident twins are arrested at the barrier twin is observed.
Our model also includes a coupling between the maximum stress of the cohesive zone model and the twin variable. An algorithm is developed to detect the twin interface based on the value of the continuum twin variable in two neighbouring bulk elements. The bilinear traction-separation law in the cohesive element that connects the two bulk elements is changed accordingly.
We show that by changing the coupling coefficient between the twin variable and maximum stress, it is possible to reproduce fracture paths similar to the ones observed in experiments (Koyama et al. 2013). The effect of the twin orientation with respect to the load axis and the effect of the shear traction on damage (Ortiz and Pandolfi 1999) are also studied. The simulations shed light on the dynamics of crack nucleation and propagation at the twin boundary, which has not been studied using in-situ experiments.
Section 2 describes the crystal plasticity model for the mechanical behaviour of the bulk. Section 3 describes the cohesive zone damage model for the mechanical behaviour of fracture interfaces. In Sect. 4 the simulation results are shown. An extensive parametric analysis has been carried out to understand the effect of the interaction coefficients between twin families (Sect. 4.1), the effect of the normal traction (Sect. 4.2), the effect of the twin orientation (Sect. 4.3), the effect of the coupling coefficient between the twin vari-able and fracture stress (Sect. 4.4) and the effect of the damage induced by shear deformation (Sect. 4.5) on the fracture behaviour. Discussion and conclusions are in Sects. 5 and 6.
2 Crystal plasticity finite element framework and discrete twin model A finite strain, CPFE framework is used. The deformation gradient F is decomposed into elastic (F e ) and plastic (F p ) parts (Lee 1969): where I is the identity matrix, u is the displacement vector, X is the undeformed coordinates vector. The time evolution of the plastic deformation is defined in terms of the symmetric plastic deformation rate D p (Dunne and Petrinic 2006): This takes into account all the slip and twin systems (Kalidindi 1998). L p is the plastic velocity gradient (Asaro and Needleman 1985). s α and n α are the slip direction and normal. s β and n β are the twin direction and normal.γ α (σ ) is the stress-dependent plastic strain rate on slip system α.φ β (σ ) is the rate of change of the twin variable representing the twin system β. γ twin β is the constant plastic shear induced in a twinned region. As will be discussed in the following,φ β (σ ) depends on both the stress and the twin variable itself.
An implicit scheme is used to find the elasto-plastic decomposition of the deformation (Dunne et al. 2007). The approach used is hypoelastic (Kim 2016;Adzima et al. 2017), which means a Newton-Raphson algorithm solves for the increment of the Cauchy stress in one time increment t. First, the small strain increment ε is calculated as: where L is the velocity gradient (Grilli et al. 2018b). Then, the increment of the Cauchy stress σ is calculated using: σ 0 is the Cauchy stress at the previous increment and C is the stiffness tensor in the sample reference frame. W e is the elastic spin (Dunne and Petrinic 2006;Belytschko et al. 2014), which is a function of the elastic velocity gradient L e : L e is a function of the total velocity gradient, plastic velocity gradient and elastic deformation gradient: Equation (4) takes into account the stress increment in the corotational reference frame and the Jaumann stress rate due to the elastic spin tensor W e (Hill and Rice 1972;Adzima et al. 2017). The right-hand side of Eq. (4) is a non-linear function of the stress increment σ through Eq.
(2), therefore it can be solved by the Newton-Raphson algorithm. L and F are fixed quantities during the iterations of the Newton-Raphson algorithm. At every iteration, the stress value σ determines the updated plastic deformation gradient F p and L p , which in turn determine F e . Therefore, W e is determined by Eqs. (5)-(6) and it is used in Eq. (4). In the present formulation, W e does not contribute to the Jacobian of the Newton-Raphson algorithm. Once the stress increment σ is calculated, the Cauchy stress σ at the next time increment can be found. This is implemented as a UMAT routine for Abaqus (Smith 2009): at each time increment, the UMAT subroutine provides the deformation gradient as input and the Cauchy stress and Jacobian are required as output. The plastic deformation gradient is stored as an internal variable (Das et al. 2018).
The constitutive model for the plastic strain rate on each slip system is described by a power-law relationship betweenγ α and the resolved shear stress (RSS) τ α (Asaro and Needleman 1985): whereγ 0 and n are constants. τ c α is the critical resolved shear stress (CRSS) of the slip systems: it is calculated based on the Taylor hardening equation using a dislocation density based model , whose details are described in Grilli et al. (2020b). Therefore, the CRSS for slip is proportional to the square root of the dislocation density (Grilli et al. 2015). The state variables used in the slip model are the forest dislocation density on each slip system and the substructure dislocation density. The details of the physical interpretation of these two types of dislocation densities are reported in (Tomé et al. 2002). Hardening of the slip systems is included (Capolungo et al. 2009;Irastorza-Landa et al. 2016). The CRSS for slip and hardening coefficients have been calibrated using neutron diffraction (Grilli et al. 2020a;Earp et al. 2018) and digital image correlation (Grilli et al. 2020c) experiments.
The rate of change of the twin variableφ β (σ ) is modelled using two terms: The first term is stress-dependent and, similarly to the slip systems, is a power-law relationship. However, twin systems are activated only by a positive RSS (Ogata et al. 2005): γ 0 and n are constants. τ β and τ c β are the RSS and CRSS on the twin system β respectively.
The twinning process involves the shuffling of atoms from their initial positions in the parent lattice to their final positions in the twinned lattice. A potential barrier has to be overcome for this process. When atoms reach the maximum of the Gibbs free energy barrier, a driving force brings them to their final positions (Ishii et al. 2016). This mechanism is described by the second term in Eq. (8): where 1/ f is the characteristic time for complete twin formation. The CRSS for twinning includes different terms describing twin nucleation and the interaction between different twin families. Twin nucleation is described by the following term (Qiao et al. 2016): where τ 0 is a constant. It is a bilinear law which has a minimum corresponding to the maximum of the energy barrier for twinning. It contributes to stress relaxation during twin nucleation. Two interaction mechanisms are considered, as shown in Fig. 1a. The increase of the twin thickness is due to the interaction of mobile dislocations with the twin boundary (Chowdhury and Sehitoglu 2018). This first mechanism is represented by the green dislocations in Fig. 1a, which can move towards the direction indicated by the corresponding arrows. The presence of residual dislocations inside the twin can prevent these mobile dislocations reaching the twin boundary (Ojha et al. 2014), limiting the growth of the twin thickness. This is modelled by a non-local interaction term: where τ 0 twin-nl is an interaction constant. The integral in Eq. (12) is calculated over a cylindrical volume β centred at the point of interest with the axis parallel to the twin plane normal n β . The axis of the integration region β has length l 0 and the radius is r 0 , as shown in Fig. 1a. The interaction coefficient τ 0 twin-nl has been calibrated by comparing the number of twins and their thickness as a function of the applied strain to in-situ electron backscatter diffraction (EBSD) experiments on α-uranium (Grilli et al. 2020d).
The second mechanism is the interaction between two discrete twins belonging to different twin systems. Dislocations are present at the boundary between the twin tip and the original crystal (red dislocations in Fig. 1a). The motion of these dislocations, called twinning dislocations, leads to the propagation of the twin tip along the twin direction s β (Mahajan and Chin 1973). However, the motion of twinning dislocations is hindered by twin boundaries belonging to a different twin system (Zhao et al. 2018) because of the energy required for dislocations to enter the twinned region and to reorient the crystal lattice (Paramatmuni et al. 2020). This is modelled by a local interaction term: where τ 0 twin-loc is an interaction constant. The sum in Eq. (13) is taken over all twin systems other than the one considered, β. This represents the barrier for the motion of twinning dislocations caused by other twin systems. The term in Eq. (13) represents the main additional feature of this model compared to our previous discrete twin model (Grilli et al. 2020b). Finally, the total CRSS for twinning is given by: Reorientation of the crystal lattice due to twinning is also considered. A weighted average using the twin variable ) is used to find the effective stiffness C: where C mat is the stiffness matrix of the untwinned crystal, while C β twin is the stiffness matrix of the crystal after a twin of type β is formed. The matrices in Eq. (15) are expressed in the sample reference system. The effective stiffness is then used in Eq. (4).
The rotation matrix that transforms the untwinned crystal into the twinned crystal of type β, with normal n β , is (Zhang et al. 2010): The stiffness matrix of the twinned crystal is found by applying the rotation matrix in Eq. (16) to the stiffness matrix of the untwinned crystal: α-uranium material parameters are used for the crystal plasticity model in the following simulations. The slip and twin systems used are reported in Table 1. The elastic constants in the crystal reference frame are reported in Table 2. The parameters of the crystal plasticity model are reported in Table 3. The crystal plasticity UMAT is available in the following repository (Tarleton and Grilli 2020).

Cohesive zone damage model
The mechanical behavior of a zero thickness interface element, also called a cohesive element, is described by  Table 2 Elastic constants (GPa) at room temperature for the orthorhombic structure of α-uranium (Fisher and McSkimin 1958;Beeler et al. 2013) in Voigt notation in the crystal reference frame (untwinned crystal lattice) a bilinear traction-separation law (Grilli et al. 2021b). In the present model we introduce the coupling between the maximum stress that the interface can withstand and the twin variable jump across the interface, which is the discontinuity of a twin variable between neighbouring bulk elements.
Zero thickness interface elements with 8 nodes are generated at the interfaces between all hexahedral bulk elements (Nguyen 2014b). Four of these nodes are shared with a face of the first neighbouring bulk element and the other four with a face of the second neighbouring bulk element. A reference system is defined on the mid-plane between the two faces of the bulk elements that are connected by an interface element. Details are reported in (Grilli et al. 2021b). The separation vector = ( s1 , s2 , n ) is defined in the mid-plane reference system. n represents the normal separation between the bulk elements; s1 and s2 represent the components of the relative displacement which lie parallel to the interface plane, also called shear separation in the following. This approach is useful to model arbitrary crack paths, which are normally obtained using continuum damage models (Sistaninia and Niffenegger 2015) or phase field fracture models (Francfort and Marigo 1998;Duarte et al. 2018;Grilli et al. 2018a).
Given an interface element and the two neighbouring bulk elements e 1 and e 2 , the twin variable jump across the interface is defined as: where ϕ β (e 1 ) and ϕ β (e 2 ) are the twin variables in the bulk elements, averaged over the integration points. The max operator in Eq. (18) is used to identify the twin interfaces of all twin systems using the unique scalar quantity ϕ. The twin variable jump determines the maximum normal stress that the cohesive element can withstand during the load history: where σ 0 max is the strength in the absence of twins and c ft is a coupling coefficient.
In the following, the bilinear traction-separation law used to calculate the normal and shear tractions at the interface between bulk elements is reported. The traction vector is used to calculate the contribution to the virtual work due to the interface elements, as described in details in (Grilli et al. 2021b;Park and Paulino 2012). This contribution is introduced in Abaqus as a UEL subroutine (Smith 2009).
The calculation of the normal and shear tractions requires the definition of a damage variable D, which is given in the following. In the present model, damage can be induced by normal and shear separation. Therefore, an effective interface separation is defined for n > 0: where: and β is a coefficient that determines the ratio between the shear and normal tractions that induce damage. β has been determined experimentally by studying the crack path during compression tests under confinement Ravichandran 1994, 1996). β 2 is usually lower than 0.5 (Ortiz and Pandolfi 1999), however no precise value is available in the literature for the material of interest. Therefore, a parametric study is carried out in Sect. 4.5. As shown in Fig. 1b, an effective traction T eff is calculated using a bilinear law. When the damage variable D is zero, at the beginning of the simulation, T eff grows linearly with eff . A critical effective separation is defined as: where K n is the normal stiffness. When eff becomes greater than 0 ( ϕ), damage develops and T eff decreases linearly with increasing eff , as shown by the part of the curve with a negative slope in Fig. 1b. f is the effective separation for complete damage: if eff ≥ f , the effective traction T eff becomes zero. Damage is considered an irreversible process. Therefore, for an arbitrary history of the variable eff , which can increase and decrease, the maximum value of T eff that can be reached is the value corresponding to eff = p , where p is the maximum between 0 ( ϕ) and the maximum value of eff reached during the load history. Also, p cannot become larger than f . This is shown by the black bilinear curve in Fig. 1b. p can be expressed as a function of the damage variable 0 < D < 1 and the twin variable jump ϕ (Grilli et al. 2021b): Notice that, if D = 0, p (0, ϕ) = 0 ( ϕ), and if D = 1, p (1, ϕ) = f . Equation (23) (D, ϕ). It is equivalent to assuming that the maximum traction that can be reached along the black bilinear curve in Fig. 1b, at eff = p , is K n (1 − D) p . Therefore, the stiffness of the cohesive interface element degrades by a factor (1 − D) (Skamniotis et al. 2019).
The time evolution of the damage and the calculation of the normal and shear tractions are explained in the following. Assume the damage variable and twin variable jump are known at time t. The damage at time t + dt is updated as: Therefore, damage increases only if p (t) < eff (t + dt) < f . Once the damage is updated, the normal traction at time (t + dt) is calculated as: The condition n (t + dt) < 0 represents contact, in which case the stiffness is not degraded. The shear traction is calculated as: where G s is the shear stiffness. The two components of the shear traction that lie parallel to the interface plane are given by: These components correspond to the shear separation given by s1 and s2 . The derivatives of the traction components with respect to the separation vector components are calculated in A. They are necessary to calculate the Jacobian matrix that ensures convergence of the finite element solver. The traction vector and Jacobian are transformed from the mid-plane reference frame to the global reference frame using rotation matrices, as detailed in (Grilli et al. 2021b). The parameters used for the cohesive zone damage model are reported in Table 4.
The choice of the value of the stiffness K n is important to avoid the artificial compliance that interface elements could introduce. Previous approaches employ irreversible cohesive elements that are activated during the simulation upon satisfaction of a failure criterion (Ortiz and Pandolfi 1999). Discontinuous Galerkin formulations of the shape functions of cohesive elements have also been used (Nguyen 2014a). In the present study, it is necessary to have all the interface elements active during the entire simulation because their maximum strength is dependent on the twin variable jump ϕ, which is a dynamically changing quantity. Therefore, crack nucleation is not only determined by the stress state but also by the twin interfaces, which are identified by the algorithm described by Eq. (18).
On the other hand, the stiffness K n cannot be too large because this causes numerical problems, such as spurious traction oscillations (Schellekens and De Borst 1993;Turon et al. 2006). The selected value of K n must be greater than E/l e , where E = 226 GPa is the Young's modulus of single crystal α-uranium (Fisher and McSkimin 1958), averaged over all directions . l e is a characteristic length of the bulk material between two interface elements: initially, l e is of the order of the element size (0.25 µm), which gives E/l e ≈ 904 GPa/µm. This is less than the value of K n used in the simulations, given in Table 4. Eventually, as will be shown in the following sections, damage localises at the twin interfaces, whose spacing is of the order of 10 µm. When this happens, l e becomes the effective length between two interface elements that accommodate a significant fraction of the crack opening, i.e. the interface elements at the twin boundaries. Therefore, after damage localisation, the selected value of K n becomes much greater than E/l e .
The implicit Abaqus solver is used for the simulations. Convergence of the fracture model without any numerical damping is possible because the condition proposed by Gao and Bower is satisfied (Gao and Bower 2004). That condition requires that the energy dissipated during the separation of an interface is larger than the bulk elastic energy released. If this condition is not satisfied a snap-back instability occurs. The condition can be expressed as (Gao and Bower 2004): where l e is used as the characteristic length because the interface elements are present between all pairs of neighbouring elements. The condition in equation (29) is largely fulfilled by the parameters used in the simulations, therefore no arc-length method is needed for convergence.
Simulations with and without cohesive elements are carried out to find the effect of the cohesive elements on the mechanical response of the representative volume. The results will be shown in Sect. 4.2.
The UEL implementation of the cohesive zone damage model, coupled with the twin variable, and input files are available in a Github repository (Tarleton and Grilli 2020). Some utilities to generate the cohesive elements are also available (Grilli 2020;Grilli et al. 2021c).

Simulation details and results
In this section single crystal simulations are carried out. The representative volume is a parallelepiped with size l m × l m × 0.25 µm, as shown in Fig. 2a. Two different values for l m are used in this section: l m = 20 µm and l m = 30 µm. The mesh is constituted of hexahedral elements and the average element size in the x-y plane is 0.2 µm for the smaller representative volume and 0.25 µm for the larger one. Such a small mesh size is needed to capture the twin thickness, which is about 2 µm, and the twin interface, which spreads over one or two elements. The effect of the ratio between l 0 and mesh size on the twin variable has been investigated in a previous work (Grilli et al. 2020b). Moreover, the element size cannot be much smaller than the separation for full damage f , otherwise, the crack opening displacement would spread over multiple neighbouring cohesive elements and the crack path would not be well defined. Therefore, 0.25 µm is a sensible choice.
Zero thickness cohesive elements are introduced at the interfaces between all hexahedral bulk elements. The elements have random shapes. A structured mesh, in which all elements are cubes, is not used because it leads to artifacts in the twin variable ϕ β , and the fracture path tends to follow the x and y directions. Therefore, elements are generated with a random arrangement and orientations thanks to the transfinite algorithm in the Gmsh software (Geuzaine and Remacle 2009).
Tension boundary conditions are used in all the following simulations, as shown in Fig. 2b. 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 x is imposed on the surface x = l m . Deformation is increased until an average strain of about 15% is reached. Therefore the maximum displacement applied on the x = l m surface is 3 µm for l m = 20 µm and 4.5 µm for l m = 30 µm. The strain rate used in all the following simulations is 10 −3 s −1 . For this order of magnitude, the simulation results are approximately strain rate independent (Tabourot et al. 1997). This is because of the small strain rate sensitivity exponent 1/n in Eq. (9). The time step used in the simulations is 0.02 s, therefore about 5000 time steps are needed to reach 10% strain.
The single crystal orientation is determined by a rotation matrix R, which transforms the axes of the crystal lattice reference frame into the axes of the sample reference frame. The reference rotation matrix is: and an additional rotation around the Z axis is added in Sect. 4.3 to investigate the dependence of the fracture behaviour on the twin orientation. The rotation matrix R reorients one twin system (β = 2) at 45 • with respect to the load direction. This will be called the primary twin system. Both {130} twin systems (β = 1) and (β = 2) are active for this crystal orientation (Grilli et al. 2020b) and tensile load along the X axis. The twin system β = 1 will be called the secondary twin system in the following.

Twin-twin interaction coefficient
Simulations with different values of the local interaction coefficient between different twin systems τ 0 twin-loc in Eq. (13) are carried out. This is useful to find the critical value of τ 0 twin-loc for which the transmission of an incident twin through a barrier twin is prevented.
Cohesive elements are not present in these simulations, therefore fracture cannot take place. Figure 3 shows the sum of the twin variables ϕ 1 +ϕ 2 at 5% strain deformation. The violet colour indicates a value, ϕ 1 + ϕ 2 = 2, which is a twin intersection. This is visible inside the white dashed square in Fig.  3a. The area occupied by twin intersections decreases if τ 0 twin-loc is increased. A value τ 0 twin-loc = 200 MPa is sufficient to limit considerably the twin interpenetration and it will be used in all the following simulations. This is because twin intersections are not commonly observed in αuranium (Zhou et al. 2016).

Twin induced fracture and normal traction
A simulation with l m = 30 µm and with cohesive elements is carried out to understand the correlation between discrete twins and fracture. The discrete twins and fracture path at 8.5% strain are shown in Fig. 4. Fracture appears mainly along the interface of the primary twin system, which is oriented at 45 • with respect to the X axis. Some other crack surfaces follow the interface of the secondary twin system or are normal to the load direction.
It turns out that the main crack nucleates at the intersection between two intersecting discrete twins, as shown in the magnified image in Fig. 4. The secondary twin develops later in the simulation; therefore, it is considered an incident twin. The primary twin acts as a barrier twin. The crack nucleates in a position where the twin shear induces a tensile stress on the interface of the barrier twin because of the lower σ max in Eq. (19).
For the same reason, after nucleation, the crack propagates along the interface of the barrier twin towards the bottom of the representative volume. As shown in Fig.  4, above the nucleation site, the crack initially tends to propagate along the interface of the incident twin because of the lower σ max , but later deviates along a path that is approximately perpendicular to the loading direction.
To understand the dynamics of crack nucleation and propagation, it is useful to correlate the twin variable ϕ 1 of the incident twin with the normal traction T n along the fracture surface indicated by the arrow in Fig. 4. This normal traction can be found from the values in the cohesive elements along that surface and is calculated using Eq. (25). A node path is selected, which includes nodes in the front surface (the surface with the highest z coordinate) in Fig. 4; therefore, the selected interface elements are determined unambiguously as 28 N. Grilli et al. The normal traction and twin variable along the developing fracture surface are shown in Fig. 5. Different strain values are shown: at 6.5% strain, the twin variable ϕ 1 has increased to almost 0.5, which is the threshold for twin nucleation; at 7.0% strain the incident twin is almost fully nucleated on the right side of Fig. 5b; the stress decreases in the centre of the twin but remains high on the left side of it. This causes a crack to nucleate on the left side of the twin, where a tensile stress is present, as shown in Fig. 4. During crack nucleation, at 7.5% and 8.0% strain, the traction decreases because of the bilinear law in Eq. (25).
The same simulations are made with and without cohesive elements to understand the effect of the compliance of the cohesive elements on the simulation results. Stress strain curves are shown in Fig. 6; the stress is averaged over the right surface in Fig. 4, where the load is applied. The stress strain curves with and without the cohesive elements are very similar up to the strain level at which cracks start to nucleate. The initial stress drop is due to the nucleation and propagation of the twins on the primary twin system (ϕ 2 ). Meanwhile, the stress with cohesive elements decreases up to complete fracture at about 10% strain.
It is noteworthy that the stress does not drop suddenly but it decreases slowly with increasing strain. This is caused by two factors: the first factor is the plastic deformation that can accommodate a large fraction of the applied displacement; the second factor is caused by the small representative volume, compared with a macroscopic sample. In a larger representative volume, the displacement applied after crack nucleation would be mostly accommodated by the crack propagation and it would cause only a minor increase of the average strain along the load direction up to complete failure. Fracture starts when the twin growth of the primary twin system ϕ 2 is almost complete, but plastic slip continues at all strain levels. Therefore, a large fraction of the energy dissipated is due to plasticity, while a smaller fraction is due to fracture. The dissipated energy and fracture energy will be discussed in Sect. 5.
Although the interface elements introduce some compliance in the elastic regime, the value of K n chosen does not have an influence on the evolution of plastic deformation, twinning and damage development in the material. Moreover, the twinning process induces inhomogeneous deformation and local stress concentration, which are the main factors affecting the location of damage nucleation, as shown in Fig. 4.  Fig. 4 at a 6.5%, b 7.0%, c 7.5%, d 8.0% strain

Dependence on twin orientation
Simulations with l m = 20 µm, β 2 = 0.2 and with different orientations of the twins with respect to the load direction are carried out to investigate the difference between the fracture patterns. The results are shown in Fig. 7.
In Fig. 7a the primary twin system is oriented at 45• with respect to the load direction. This is the same orientation analysed in Fig. 4. Therefore, the crack pattern follows a similar trend: cracks propagate along the interface of the primary twins, but also deviate along a path that is approximately perpendicular to the load direction. In Fig. 7b the primary twin plane is oriented at 22.5 • with respect to the load direction. In this case, the normal traction on the surface of the primary twin is lower, therefore the crack propagates almost perpendicular to the load direction. Crack propagation along the secondary twin system is also observed. In Fig. 7c the primary twin plane is oriented at 67.5 • with respect to the load direction. This orientation of the primary twin plane is sufficient for the fracture surface to appear exclusively on that plane.

Twin variable coupling coefficient
The parameter that affects most the crack pattern is the fracture-twin coupling coefficient c ft . Simulations with l m = 20 µm, β 2 = 0.2 and three different values of c ft have been carried out. Only the primary twin system is used to decouple the effect of c ft from the effect of twin intersection. The results are shown in Fig. 8.  Figure 8a shows the case in which no coupling is present (c ft = 0). As expected the crack pattern is not strongly correlated with the twin variable and fracture surfaces form that are perpendicular to the load direction. For c ft = 200 MPa the fracture surface is mostly perpendicular to the load direction; it passes through the central twin plane without being deviated and then it follows the twin interface of the lower twin. For c ft = 400 MPa about half of the crack length is parallel to the twin plane while the rest is perpendicular to the load direction. This is consistent with the simulations in the previous sections, in which the value c ft = 400 MPa has been used.
The strong sensitivity of the fracture path with respect to the c ft parameter allows an approximate value of c ft to be determined by comparing with experiments in the literature, as discussed in Sect. 5.
The same value of f is used for the interface elements in the bulk and at the twin interface. Physically, this is because the elongation at failure is determined by the plastic deformation. In α-uranium, plastic deformation can take place both in the untwinned and in the twinned crystal (Cahn 1953). Therefore, no large difference in the elongation to failure near a twin boundary or in the bulk material is expected. By contrast, the strength of the twin interface is much lower because of the different atomic arrangement compared with the bulk. Breaking atomic bonds at the twin interface requires lower stress.  The β 2 parameter, controlling the shear contribution to damage, is not known exactly for most metals. It is usually in the interval between 0 and 0.5 (Ortiz and Pandolfi 1999). To understand how the β 2 parameter affects the fracture path, simulations with l m = 20 µm and different values of β 2 have been carried out. The simulation results are shown in Fig. 9.
The features of the crack pattern are not strongly affected by the β 2 parameter. However, it can be noticed that higher values of β 2 lead to fracture surfaces that tend to propagate more along the twin interfaces. Therefore, increasing β 2 and increasing the fracturetwin coupling coefficient c ft has a similar effect on the fracture pattern.
Despite the different values of β 2 in Fig. 9a-c, the strain necessary for the propagation of the crack in most of the representative volume (7.8%, 7.4% and 8.4% in Fig. 9a-c respectively) has similar values.
These results show that the normal traction T n remains the most important factor for crack nucleation and propagation.

Discussion
The results in Sect. 4.1 show that the interpenetration between two twin families is strongly affected by the interaction coefficient τ 0 twin-loc . A value of 200 MPa is sufficient to limit considerably the twin intersections. According to Eq. (13), this value corresponds to the increase of the CRSS required to propagate an incident twin through a barrier twin. Experiments by Zhou et al. report evidence that the propagation of incident twins can be completely stopped by barrier twins in αuranium (see for instance Fig. 7 in Zhou et al. (2016)). But this is true only for some twin combinations. For instance, Cahn observed intersections between (172) and (130) twins (Cahn 1951).
The CRSS for twin nucleation used in this study (τ 0 = 25 MPa) was calibrated using neutron diffraction experiments during tensile tests (Grilli et al. 2020a). Therefore, the present simulations suggest that the RSS to propagate an incident twin through a barrier twin in α-uranium can be almost one order of magnitude higher than the RSS for twin nucleation. This is not the case Fig. 9 Twin variable and fracture pattern depending on the shear contribution to damage as a function of β 2 a β 2 = 0.1, 7.8% strain, b β 2 = 0.3, 7.4% strain, c β 2 = 0.4 8.4% strain. The displacement of the mesh is multiplied by 2 to highlight the crack pattern. Note that the case β 2 = 0.2 is reported in Fig. 7a for other metals, such as TiAl, in which intersecting twins are commonly observed (Sun et al. 1993).
Twin pairs that do not interpenetrate lead to stress concentration at the interface between incident and barrier twin, as shown in Sect. 4.2. The tensile stress induced by the incident twin can be sufficient to initiate fracture. This has been observed at grain boundaries in TiAl (Bieler et al. 2005), but in α-uranium this phenomenon can take place also at the twin intersections (Taplin and Martin 1965). The present simulations are necessary to understand the stress level that is sufficient to nucleate and propagate cracks.
The results in Sect. 4.4 show that a decrease of 400 MPa in the strength of the twin interface, with respect to an original value of 650 MPa in the bulk, produces a fracture pattern with the following feature: about half of the crack surface lies at the twin interface and the other half is approximately perpendicular to the load direction. The value of the fracture strength (650 MPa) was calibrated by comparing the strain to fracture found using polycrystal simulations (Grilli et al. 2021b) and the one measured in fracture experiments on α-uranium (Huddart et al. 1980). A similar fracture pattern with a staircase shape has been observed in TWIP steel by Koyama et al. (Koyama et al. 2013). Therefore, the present simulations indicate that the strength of the twin interfaces, in materials exhibiting cracks at twin boundaries, can be as low as one third of the bulk value. Atomistic simulations would be required to understand the origin of the decrease in strength of the twin interface (Yamakov et al. 2006;Sakano et al. 2020).
The value of the fracture energy of the interface elements corresponds to the area below the tractionseparation curve in Fig. 1. It is 162.5 J/m 2 with the parameters used in the simulations. This value is lower than characteristic values used for cohesive zone models of steel (Ortiz and Pandolfi 1999) and Ni alloys (Citarella et al. 2018). It is typical of low ductility metals. The total energy dissipated in the tensile load simulations can be approximated by the area below the stress-strain curve in Fig. 6. It turns out to be 878 J/m 2 , therefore only 18% of the energy is dissipated by the fracture surfaces and the rest is dissipated by plastic deformation. This is the cause of the slow decrease of the stress during displacement controlled tensile load simulations.
In the experiment by Koyama et al. (2013), the primary twin is oriented at about 65 • with respect to the tensile direction (see Fig. 9 in Koyama et al. (2013)). The fracture path has a staircase shape; first the crack propagates along the twin boundary, then it deviates at about 90 • before reaching the neighbouring twin interface. This process repeats over and over, and the final fracture surface is approximately perpendicular to the tensile axis. This mechanism is well reproduced by the present simulations. In Sect. 4.3, a simulation with the primary twin oriented at 67.5 • with respect to the tensile direction is shown (Fig. 7c). The crack surface lies mostly on the primary twin plane and it deviates at about 90 • to connect the different twin planes. In this case, the role of the weak interface of the secondary twin system is also important. Indeed, secondary twins can intersect multiple primary twins and represent weak links along which the crack can propagate. Therefore, these simulations suggest that the activation of two twin systems can favour the formation of a fracture path with a staircase shape, in which the crack surface can propagate quickly along the two different twin planes.
As shown in Sects. 4.2 and 4.5, the main stress component that affects crack nucleation at the twin interface is the normal traction. If the shear contribution to damage development is increased, the fracture path remains similar (Fig. 9). Therefore, it can be concluded that twin intersections do not cause a large shear displacement parallel to the twin interfaces, even if they can do so along other planes in the bulk material. The cohesive zone model is used in this work as an empirical model to deduce the most likely mechanism of failure, and to quantify the strength and fracture energy of the twin interface. It is calibrated using mesoscale observations, but we do not attempt to use it to describe fracture at the nanoscale. The resolution of our model is limited by the element size, which is about 0.25µm and the resolution of the fracture surface is limited by this characteristic length. Therefore, the roughness of the fracture surfaces in the simulation results in this manuscript is determined by the element size and should not be regarded as the actual roughness.
As shown by the simulation results, the fracture path at the length scale of the representative volume (about 30 µm) is not determined by the element size but it depends on the load direction and on the orientation of the twin interfaces.

Conclusions
In this paper, a continuum model for discrete twins is coupled with a cohesive zone model to understand the nucleation and propagation of cracks at twin boundaries. Zero thickness interface elements are placed at the interfaces between all hexahedral bulk elements to simulate intragranular crack surfaces with arbitrary shapes.
The maximum stress of the cohesive zone model is coupled with the continuum twin variable to model the reduction of the strength of twin interfaces. Therefore, the strength of a twin interface, when a normal traction is applied, can be determined with respect to the strength of the bulk material by comparing the simulated and observed fracture surfaces. The single crystal simulations with two active twin systems show that the strength of twin interfaces can be reduced by about 400 MPa with respect to the strength of the bulk. Moreover the simulations reveal that the RSS required to propagate an incident twin through a barrier twin in αuranium can increase by about 200 MPa with respect to the RSS for twin nucleation, which is only 25 MPa.
The developed model can help understand if the activation of multiple twin systems can lead to fracture in different materials.
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: Jacobian calculation
In this section, the derivatives of the traction components in Eqs. (25), (26), (27), (28) with respect to the separation vector components are reported.
In the following the implicit assumption that all the derivatives are calculated at time t + dt is made. This is because the damage and its derivatives are expressed at the end of the time step.
Since the traction components are functions of the damage variable, it is convenient to write the derivatives of the damage variable with respect to the separation vector components: The derivatives of the normal traction depend on the crack opening or closure condition: The derivatives of the shear traction do not depend on the sign of n (t + dt): The derivatives along the two shear directions are: The derivatives of the normal traction with respect to shear separation and vice-versa are: