A theoretical and computational investigation of mixed mode creep crack growth along an interface

In this paper, we propose a theoretical framework for studying mixed mode (I and II) creep crack growth under steady state creep conditions. In particular, we focus on the problem of creep crack growth along an interface, whose fracture properties are weaker than the bulk material, located either side of the interface. The theoretical framework of creep crack growth under mode I, previously proposed by the authors, is extended. The bulk behaviour is described by a power-law creep, and damage zone models that account for mode mixity are proposed to model the fracture process ahead of a crack tip. The damage model is described by a traction-separation rate law that is defined in terms of effective traction and separation rate which couple the different fracture modes. Different models are introduced, namely, a simple critical displacement model, empirical Kachanov type damage models and a micromechanical based model. Using the path independence of the C∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C^{*}$$\end{document}-integral and dimensional analysis, analytical models are developed for mixed mode steady-state crack growth in a double cantilever beam specimen (DCB) subjected to combined bending moments and tangential forces. A computational framework is then implemented using the Finite Element method. The analytical models are calibrated against detailed Finite Element models and a scaling function (Ck\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{k}$$\end{document}) is determined in terms of a dimensionless quantity ϕ0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{0}$$\end{document} (which is the ratio of geometric and material length scales), mode mixity χ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi $$\end{document} and the deformation and damage coupling parameters. We demonstrate that the form of the Ck\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{k}$$\end{document}-function does not change with mode mixity; however, its value depends on the mode mixity, the deformation and damage coupling parameters and the detailed form of the damage zone. Finally, we demonstrate how parameters within the models can be obtained from creep deformation, creep rupture and crack growth experiments for mode I and II loading conditions.

against detailed Finite Element models and a scaling function (C k ) is determined in terms of a dimensionless quantity φ 0 (which is the ratio of geometric and material length scales), mode mixity χ and the deformation and damage coupling parameters. We demonstrate that the form of the C k -function does not change with mode mixity; however, its value depends on the mode mixity, the deformation and damage coupling parameters and the detailed form of the damage zone. Finally, we demonstrate how parameters within the models can be obtained from creep deformation, creep rupture and crack growth experiments for mode I and II loading conditions.
Keywords Creep · Crack · Mixed mode · C*integral · Damage zone model · Traction-separation rate law (TSRL) · Double cantilever beam (DCB) · Dimensionless analysis Nomenclature x i The Cartesian material and spatial coordinates (i = 1, 2, 3) (m) C * The rate of the J -integral (J/m 2 /h) C s The separation history function of the simple model C k The separation history function of model k σ i j The Cauchy stress tensor (MPa) ε i j The engineering strain tensor ε i j The engineering strain rate tensor (1/s) κ The curvature of a beam (1/m) κ The curvature rate of a beam (1/(m s)) (•) el The elastic component of the quantity (·) (•) cr The creep component of the quantity (·) s i j The deviatoric part of Cauchy stress tensor (MPa) σ e The von Mises equivalent stress (MPa) E Young's modulus (MPa) σ 0 The reference stress (MPa) ε 0 The strain-rate at the reference stress σ 0 (1/s) n The rate sensitivity parameter of the bulk material u i The displacement vector (i = 1, 2, 3) (m) n i The unit normal vector (i = 1, 2, 3) T i The traction vector (i = 1, 2, 3) (MPa) T e The effective traction (MPa) δ i The displacement jump vector across the damage zone (i = 1, 2, 3) (m) δ i The displacement jump rate vector across the damage zone (i = 1, 2, 3) (m/s) δ The representative separation across the damage zone (m) δ e The effective displacement jump across the damage zone (m) δ e The effective displacement jump rate across the damage zone (m/s) δ c The representative separation at which damage initiates in the damage zone at the crack tip (m) δ f The representative separation at failure in the crack tip (m) δ m i The maximum displacement jump rate vector in the crack tip (i = 1, 2, 3) (m/s) δ m e The maximum effective displacement jump rate at the crack tip (m/s) T 0 The reference traction of the damage zone (MPa) δ 0 The separation rate at the reference traction T 0 (m/s) m the rate sensitivity exponent of the damage zone ω A scalar damage parameter α t The tangential coupling parameter α n The normal coupling parameter g t The tangential coupling function of the micromechanical model g n The normal coupling function of the micromechanical model β A material parameter of the exponential damage law f The current area fraction of the pores f 0 The initial area fraction of the pores f c The coalescence area fraction of the pores 2 l The spacing between two adjacent pores (m) h The current height of a pore (m) h 0 The initial height of the pores (m) a The crack length (m) a The steady state crack velocity (m/s) a The dimensionless steady state crack velocity φ 0 The ratio of geometric to material length scales λ The characteristic geometric length scale (m)

Introduction
In many engineering applications at elevated temperature, structural components exhibit significant timedependent inelastic deformation, which might lead to nucleation and/or propagation of flaws. Studying creep crack growth (CCG) has been an active area of research over the last four decades, with the aim of designing structures with high integrity and safety. In this paper, we aim to devise analytical models for steady-state crack growth under mixed mode loading conditions which can be calibrated against detailed Finite Element simulations. Such models can be used to investigate the effect of different material parameters and damage models on the crack growth behaviour under complex loading conditions. The majority of early studies on creep crack growth focused on pure mode I loading conditions, with the aim of developing a parameter that characterises the near crack tip fields as well as crack propagation. The so called C * -integral (Landes and Begley 1976;Nikbin et al. 1976;Ohji et al. 1976), i.e. the creep J -integral (Rice 1978), was introduced to characterise the crack tip fields and creep crack growth under steady state creep conditions. In particular, it provides descriptions of the strain-rate and stress singularities at the crack tip and a correlation of experimental crack growth rate data (Taira et al. 1979;Riedel and Rice 1980). Moreover, the C * -integral is path independent for contours in which the material properties only vary in the direction perpendicular to the direction of crack growth within the family of contours considered. Under transient conditions, i.e. the transition from short-time elastic to longtime creep, other parameters such as C(t), C t and C * h are found to characterise both stationary and propagating cracks with a variable degree of suitability (Riedel and Rice 1980;Riedel 1981;Saxena 1986).
Under mixed mode loading and steady state conditions, the C * -integral remains a valid parameter (Chambers et al. 1992). Chambers et al. (1992) assumed that only a part of the C * -integral, i.e. that part related to mode I, is responsible for propagating the crack. They defined the damage at the crack tip by a function of the von Mises equivalent stress and then used the amplitude of the singular von Mises stress field to determine an effective value of C * by comparing mode I and mixed mode loading conditions. Further, under transient creep conditions, the role of mode mixity on the stress field ahead of a crack tip was investigated by Brockenbrough et al. (1991) using the C(t) parameter.
The damage process ahead of a propagating crack tip in polycrystalline materials under creep conditions is controlled by the formation of discrete voids or micro cracks at grain boundaries. As these voids and/or micro cracks grow and coalesce, they form secondary cracks which coalesce creating the new crack surfaces, allowing the primary macroscopic crack to advance along an interface or interconnected grain-boundaries (Riedel 1987). These secondary cracks are usually not perpendicular to the loading direction, which can generate mixed mode loading conditions locally at the crack tip. In some materials, an extensive damage zone can develop in the vicinity of a crack tip in which many boundaries are cavitated. The interaction between these different cavitated and microcracked regions is important and can lead to a significant reduction of stress locally and macroscopic dilation (i.e. diffused damage). In other materials, damage is confined to a small region directly ahead of the crack tip such that the dilation due to the growth of the damage is then very localised. In the presence of interfaces, whose fracture properties are weaker than the bulk material, e.g. in dissimilar metal welds (Yamazaki et al. 2008;Laha et al. 2012;Hu et al. 2019), a crack might be forced to grow along an interface, which might lead to mixed mode loading at the crack tip, regardless of the global loading conditions. Mode mixity can influence damage development.
The mode mixity effect is usually incorporated using effective measures such as the effective stress, effective creep strain, or accumulated creep strain at a critical distance ahead of the crack tip. In these models the deformation of the bulk material is assumed not to be influenced by the presence of damage. Using empirical Kachanov type continuum damage mechanics models (Kachanov 1958;Rabotnov 1969), the single scalar damage parameter is assumed to be determine by the stress state which allows the effect of mode mixity to be taken into account (Nikbin et al. 1976(Nikbin et al. , 1984Yatomi and Nikbin 2014). The development of the damage process zone at a crack tip has been modelled by Onck and van der Giessen (Onck and van der Giessen 1998;Van Der Giessen and Tvergaard 1994) using the Finite Element Method, together with both mechanistic and empirical models. The latter models fully describe the interaction between deformation and damage development and account for mode mixity. In these models, the contribution of the different fracture modes (i.e. the coupling) to the deformation and the damage in the damage zone are assumed to be the same.
Another method of modelling crack propagation is through the use of interface cohesive or damage zone models. Interface damage zone models of this type provide a coupling between the local separation rate across an interface and bulk deformation processes. They introduce a physically meaningful length scale that is related to the dissipative mechanisms responsible for damage development. Damage zone models of this type describe the fracture process in the vicinity of the crack tip as a gradual surface separation process, such that the normal and shear traction at the interface resist separation and relative sliding. The cohesive/damage zone modelling approach has its origins in the pioneering work of Dugdale (1960) and Barenblatt (1962) and was later implemented in a Finite Element environment by Hillerborg et al. (1976). Several models have been proposed in the literature, wherein a variety of materials and applications have been successfully investigated (Camacho and Ortiz 1996;Corigliano et al. 2003;Elmukashfi and Kroon 2014;Hui et al. 1992;Knauss 1993;Needleman 1987Needleman , 1990Rahul-Kumar et al. 1999;Rice and Wang 1989;Tvergaard 1990;Xu and Needleman 1993). These first attempts focused on studying the mode I fracture process, i.e. the relation between the traction and separation that are normal to the fracture surface. Later, it was extended to mode II fracture, in which the tangential traction and sliding are defined. In modelling mixed mode loading conditions, two approaches are often used to define the interaction between different fracture modes. The first approach is by fully coupling the normal and tangential components of the traction and separation through introducing an effective traction and/or separation parameters (Tvergaard and Hutchinson 1996;Pardoen et al. 2005). The damage is evaluated using a representative separation that determines the contribution of the different fracture modes and yields different fracture energies for the different modes. The second approach assumes uncoupled behaviour between the different modes, i.e. uncoupled normal and tangential deformation, but the total fracture energy is taken to be the sum of the fracture energy components of the different modes (Kafkalidis and Thouless 2002;Li et al. 2006). Experimental studies show that the fracture process depends on the properties of the different materials, and the coupling between the different modes strongly depends on the nature of the process zone.
The objective of this paper is to study a simple geometry, in which we project the damage onto a single plane directly ahead of the crack tip, and a series of loading conditions to provide fundamental information about the effect of mode mix on crack growth. Simplifying the model in this way allows us to undertake a rigorous analysis that provides important insights that will ultimately help us to explore and explain more complicated and realistic scenarios. Therefore, we developed a theoretical and computational framework for creep crack growth under mixed mode conditions. The crack growth is assumed to be determined by the damage in a concentrated narrow zone directly ahead of the crack tip (i.e. the damage in not diffused in the bulk material) in a weak interface that is parallel to the crack. A theoretical framework is initially introduced in which the constitutive behaviour of the bulk material is described by a power-law creep and the damage zone model constitutive relation is described by a traction-separation rate law. The damage models are formulated in terms of effective traction and separation parameters to account for mode mixity effects. The damage process is assumed to be determined by a representative separation which determines the damage dependence on different modes. We propose three different types of models, i.e. a simple critical displacement model, Kachanov type empirical models and a micromechanical based interface model, as extensions to the pure mode I models described in Elmukashfi and Cocks (2017). We use the same approach as described in Elmukashfi and Cocks (2017). We initially evaluate the behaviour of a double cantilever beam specimen (DCB) of infinite length subjected to combined bending moments and tangential forces. For this type of specimen, C * remains constant as the crack grows and the crack growth rate eventually achieves a steady state. By exploring the steady state response, we can establish the relationship between the steady state crack growth rate, C * and characteristic material and geometric length scales. The resulting models can then be used to evaluate more complex loading conditions and geometries. For the double cantilever beam, reasonably straightforward analytical expressions can be obtained for C * . In the steady state, we can invoke the path independence of the C * -integral and choose contours in the far field and surrounding the damage zone. This allows simple analytical expressions for the crack growth rate to be determined, which can be expressed in terms of C * , damage zone material parameters, a dimensionless scaling parameter that is a function of the ratio of characteristic geometric and material length scales, and mode mixity. We use the Finite Element Method to determine the magnitude of the scaling parameter.
The theoretical framework is presented in Sect. 2 and the damage zone models are developed in Sect. 3. The analysis of creep crack growth in the double cantilever beam specimen and the Finite Element implementation of the damage zone model are described in Sect. 4, with the crack growth results for the different interface models presented and discussed in Sect. 5. Having established the theoretical framework and structure of the crack growth laws, we evaluate the behaviour of more representative testing geometries in Sect. 5.

Theoretical analysis of mixed mode creep crack
In this section, we investigate the theoretical aspects of mixed mode creep crack propagation. In particular, we adopt the theoretical framework introduced in Elmukashfi and Cocks (2017) and consider mixed mode crack propagation in creeping materials under steady state conditions. A full description and evaluation of the nature of crack tip fields for stationary and growing cracks is provided in Elmukashfi and Cocks (2017) for mode I conditions. Here the focus is on mixed mode and mode II conditions, which has achieved much less attention in the literature. Crack propagation is assumed to take place along an interface whose fracture properties are weaker than the bulk material (i.e. the crack propagates along the interface regardless of the far field and local mode mixities). Further, the path independence property of C * will be used to obtain a direct relationship between the far field loading and the fracture process parameters. It should be noted that for crack growth in an elastic/creeping material, path independence of C * still holds provided the size of the zone in which elastic deformation is important is smaller than the size of the crack tip damage process zone (Hui and Riedel 1981).
Consider a body which contains a crack and is subjected to a constant far field loading, see Fig. 1. A common Cartesian coordinate system for the reference and deformed configurations x i , i = 1, 2, 3, is assumed. Directly ahead of the crack there is an interface along which the crack propagates. The bulk material is assumed to exhibit steady-state creep behaviour which is defined by the constitutive laẇ where σ i j is the Cauchy's stress tensor,ε i j is the strain rate tensor, s i j = σ i j − 1 3 σ kk δ i j is the stress deviator, σ e is the von Mises equivalent stress, σ 0 is a reference stress,ε 0 is the strain rate at the reference stress and n is the rate sensitivity parameter. Φ is the stress potential The energy dissipation rate per unit volume,Ḋ, is given bẏ where Ψ is the dual rate potential andε e = 2 3ε i jεi j . In this paper, we limit our consideration to situations in which the creep properties in the bulk material either side of the interface are the same. In dissimilar metal welds different bulk materials with different properties exist either side of the interface. The results presented here can be readily extended to consider this more general situation, but in order to identify the major features of the crack growth behaviour under mixed mode loading conditions we focus on the response of this simpler configuration.
Damage grows along the interface, driven by the local stress, and the crack advances when the damage at a material point ahead of the crack achieves a critical value. Hence, separation occurs along the interface to create two surfaces as the crack advances. Therefore, a material point along the interface is defined by the two normal vectors n + i and n − i , where n + i = n − i , i.e. the initially intact material point splits into two points with unit normals acting opposite to each other and into the material on either side of the interface. The displacement-rate jump across the damage zone and the corresponding traction are defined by the vectorṡ i are the displacement rates either side of the interface, and T + i = σ i j n + j and T − i = σ i j n − j , respectively. The problem is analysed by equating the values of C * determined on the inner and outer contours in Fig. 1. The C * -integral is defined as where Γ is an arbitrary contour around the tip of the crack with unit outward normal n i , T i = σ i j n j is the traction on ds andu i is the displacement rate. The C *integral in the outer path, C * out , is determined by the far field loading. For an interface that extends along the x 1 -axis, the C * -integral along the inner path, C * in , is then evaluated as whereδ m i is the displacement rate at the tip of the crack and T i = T + i = −T − i . In the case of mixed mode loading, the components of the traction and displacement rates tangential and normal to the interface are considered, i.e. (T t ,δ t ) and (T n ,δ n ), respectively, as illustrated in Fig. 2. We postulate that there is an effective traction T e and an effective separation rateδ e , that can be used to represent the constitutive relationship for the interface. (The effective traction and effective separation rate are assumed to be functions of the traction and separation rate components, respectively). Note also that the effective traction and separation rate are power conjugates. Thus, an equivalent expression to C * in in Eq. (6) is whereδ m e is the effective displacement rate at the tip of the crack. The effective separation δ e is then defined as the integration of the effective rate, i.e. δ e (t) = t 0δ e dt. Hence, a relationship between the far field loading and the behaviour of the interface can be determined using the path-independence of C * , i.e. C * out = C * in .

Damage zone models for mixed mode creep crack growth
In this section, we introduce a number of different models to describe the response within the damage zone under mixed mode conditions. The traction-separationrate law is assumed to take a direct power law relationship between the effective traction T e and the effective separation rateδ e . We can express the resulting constitutive relationships in a similar form to that for power- The damage zone along the interface for mixed mode crack propagation: a schematic of the damage zone and the traction and separation vectors, T and δ, respectively; and b the effective traction-separation rate (T e -δ e ) and the effective separation rate-separation (δ e -δ e ) distribution along the damage zone. l is the length of the interface zone, δ c e is the critical effective separation, δ f e is the effective separation at failure,δ m e is the maximum effective separation rate and σ c is the interface strengthd law creep of the bulk described in Sect. 2. Thus, the constitutive model takes the following forṁ where T 0 is a reference traction (equivalent to σ 0 of Eq. (1)),δ 0 is the separation rate at this traction and m is an exponent (which can have a different value to n in Eq. (1)). Throughout this paper, we set T 0 = σ 0 , so that the two reference stress like quantities are the same. A traction potential, φ, can now be defined as It follows that the energy dissipation rate per unit area, d, and dual separation rate potential, ψ, are, respectively, given bẏ and The constitutive relations can also be expressed in terms of the traction and separation rate components aṡ and where the index notation i denotes the normal and tangential directions, i.e. i ∈ [t, n], which will be used in the subsequent descriptions.
In what follows, a range of traction-separation-rate laws are introduced for mixed mode conditions which generalise the models described by Elmukashfi and Cocks (2017).

Simple critical displacement model
A simple form of the effective traction can be taken in the form of the resultant traction, i.e. the square root of the sum of squares of the normal and tangential traction, T t and T n , respectively. Thus, the effective traction can be written in the general form where α t and α n are parameters that define the coupling between tangential and normal deformation as illustrated in Fig. 3. Thus, the separation rates in Eq. (12) take the following forṁ  (14), the effective separation rate,δ e , can be determined aṡ The traction in Eq. (13) is then given by The coupling between the tangential and normal deformation in the interface is determined by the parameters α t and α n as illustrated in Fig. 3 and Eqs. (14)-(17). We assume that the damage does not influence the separation rates and failure is achieved when a suitable separation measure achieves a critical value. Hence, we introduce a representative separation,δ, that is a function of the tangential and normal separations, i.e. δ =δ(δ t , δ n ), which can be used to define the local failure condition asδ =δ f whereδ f is the representative separation at failure. It is worth mentioning that the traction T t , T n and T e , become zero at the position where local failure is achieved. The form of the representative separationδ determines the contribution of the different fracture modes to the damage process. For example, assuming thatδ = δ e , the parameters α t and α n determine the coupling between both the deformation and local failure. Moreover, other forms can also be used to defineδ, for example, for the cases that mode I controls the damage process,δ ≈ δ n .

Kachanov type empirical model
In this model, damage is assumed to influence the constitutive response. The damage is introduced using single or multiple scalar damage parameters, ω i , such that each parameter is assumed to evolve monotonically from 0 to 1, i.e. from an undamaged to a fully damaged state. We assume that the effective traction takes the quadratic form of Eq. (14) with the tangential, T t , and normal, T n , traction, replaced byT t andT n , respectively, where, following Kachanov (1958) and Lemaitre and Chaboche (1994), Hence, the effective traction becomes Thus, the constitutive response is then simply obtained by replacing α i by α i (1 − ω i ) in Eqs. (15) and (17), to givė and Similarly, the effective separation rate becomeṡ The coupling between the tangential and normal deformation in the interface is controlled by the terms α t (1 − ω t ) and α n (1 − ω n ). Thus, initially in the absence of damage, the coupling is defined by the parameters α t and α n . After damage onset, i.e. ω t > 0 and/or ω n > 0, the rate at which each damage parameter evolves controls the coupling. It is worth noting that in the case of isotropic damage, i.e. ω t = ω n = ω, the damage does not affect the coupling between the tangential and normal deformation which is determined by the parameters α t and α n .
The simplest damage representation is obtained by considering isotropic damage. Here, we assume that the damage is related to the representative separationδ, i.e. ω = ω(δ). Further, the damage evolution law may take different forms depending on the damage mechanism(s). General forms of the different damage models proposed by Elmukashfi and Cocks (2017), i.e. linear and exponential models, can be expressed as: An exponential damage model: whereδ c is the representative separation at which damage initiates (for separations less than this value ω = 0 and the constitutive response is given by Eqs. (15) and (17)),δ f is the representative separation at failure and β is a material parameter. It should be noted that for the exponential damage law, the traction does not necessarily decrease smoothly to zero at failure, but an abrupt response may result. As discussed in Sect. 3.1, the form of the representative separationδ defines the contributions of the different fracture modes to the damage process. It should be mentioned that, in the case of anisotropic damage, the form of the damage evolution laws may affect these contribution. In this study, we limit our considerations to the cases of isotropic damage and the representative separationsδ = δ e and δ = δ n .

Micromechanical based model
This model is a generalisation of the previously proposed pure mode I model employed in Elmukashfi and Cocks (2017). In this model, the damage is modelled as an array of pores idealised as cylinders, which grow as the surrounding material creeps, see Fig. 3.3. At a given instant, the radius and height of the pores are denoted as r and h, respectively, and the mean spacing is 2 l. Thus, the pores are characterised by their area fraction in the (a) (b) Fig. 4 The micromechanical representation of the creep damage by pore growth: a pores of radius r and spacing 2l on a grain boundary subjected to microscopic stress state σ i j and the deformation is controlled by steady-state creep and b an idealisation of a pore as a cylinder with height h and diameter equal to the pore to pore spacing 2l and the macroscopic tangential and normal traction (T t , T n ) and separation (δ t , δ n ).
plane of the cavitated zone, i.e. by f = (a/l) 2 . Further, the representative volume element is assumed to be fully constrained in the radial direction and the deformation occurs in the normal and tangential directions, i.e. l = const. In the case of pure normal loading the pores elongate in the normal direction and also extend in the plane of the interface. Under shear loading pores become more crack like and elongate in the direction of shear. The resulting expression for the effective traction given by Yalcinkaya and Cocks (2015) is whereḡ t = g t /g 0 ,ḡ n = g n /g 0 and The constitutive response obtained using Eqs. (8) and (12) is and The effective separation rate then becomes Thus, the contribution of each fracture mode to the deformation of the interface is controlled by the functions g t ( f ) and g n ( f ), as illustrated in Fig. 5. Examining these functions, we find that the ratio g t /g n increases with increasing value of f and reaches a maximum value of g t /g n ≈ 0.5 at f = 1.0.
The evolution of the pore's area fraction and height during loading are given by Yalcinkaya and Cocks (2015) aṡ ( The initial pore's area fraction and height are assumed to be f 0 and h 0 , respectively. Further, the pores are assumed to coalesce when f = f c . These equations provide a coupling between the normal and tangential response. The area fraction of the pores grows under both normal and tangential loading, while the height of the pores increases under normal loading, but decreases in shear, as the pores become more crack like.

Mixed mode creep crack growth in a double cantilever beam
In this section, mixed mode creep crack growth in the DCB specimen shown in Fig. 6 is analysed. The length and height, of the specimen are denoted by L and 2 H , respectively, and the crack length is denoted by a. The crack propagation is assumed to take place at an interface at x 2 = 0. Each arm of the specimen is subjected to constant force and moment per unit depth, i.e. N and M 1 on the upper arm and −N and M 2 on the lower arm, respectively. It is convenient to express the moments M 1 and M 2 in terms of a bending moment M and the moments required to balance the moment due to the pair of forces N and −N , such that M 1 = M + N H/2 and M 1 = M − N H/2. Hence, pure mode I is obtained when N = 0 and pure mode II is obtained when M = 0. We assume that the overall length L → ∞ and that the height of each arm H a. Under these conditions, C * remains constant as the crack grows, thus a steady state is eventually achieved in which the crack growth rate is constant. We focus on this steady state response. The objective is to obtain a mathematical description for the relationship between the far field loading and the local damage development within the damage zone and the crack growth rate under both plane stress and plane strain conditions. Fig. 6 The schematic of the double cantilever beam specimen Fig. 7 The definition of the inner path Γ in and the outer path Γ out that are used to evaluate the C * -integral in the double cantilever beam specimen The relationship between the loading and fracture parameters is obtained using the path independence of the C * -integral as discussed in Sect. 2. The C * -integral is evaluated along the outer and inner paths indicated by the dashed lines and different colours in Fig. 7. Equating the values of C * determined from these two paths provides a relationship between the crack-tip opening rate and the applied load. As discussed in Elmukashfi and Cocks (2017), there is a single characteristic geometric length scale for this problem, which we take as λ = H/2, and in the steady state the separation rate within the interface can be expressed as a function of x 1 /λ, integrating this function as an element is convected towards the crack tip as the crack grows at constant velocity, allows the crack growth rate to be determined. The form of this function is not known, but to determine the crack growth rate we only need to determine a single quantity-the resulting integral can be determined from a single piece of information from a Finite Element analysis of the problem, namely the crack velocity. The details of this process are given below.

4.1
The C * -integral for the outer path Γ out C * for the outer path can be determined from the rate of change of the rate analogue of the total potential energy, Π , with crack length, where It should be noted that C * can also be determined by direct evaluation of the integral of Eq. (5). For an element of beam under combined bending and uniaxial loading M arm and N arm , respectively, the curvature ratė κ and axial strain rateε can be determined by solving the following nonlinear equations (for the details see Appendix A) where for plane stress and √ 3/2 for plane strain. For an arm loaded with a bending moment M arm and axial force N arm , the contribution to C * is determined as whereθ andδ are the rotation and displacement rates at the end of the beam, respectively. For the DCB spec-imen of Fig. 7, the C * is the summation of the upper and lower arms contributions. Hence, C * out is where the curvature and axial strain rates for the upper arm areκ 1 andε 1 that are determined by solving the non-linear equations of (33) It should be noted that the effective separation rate at the crack tip,δ m e , depends on the form of interface model adopted.

The effective displacement rate at the crack tip
The effective separation rate at the tip,δ m e , can be expressed in terms of the applied load using the path independence of C * , i.e. by equating Eqs. (35) and (36). This gives: where φ 0 =ε 0 λ δ 0 =ε 0 H 2δ 0 is the ratio of geometric to material length scales for the problem,C * = C * ε 0 σ 0 λ is the dimensionless value of C * and g(m) = m + 1 2 m .

The Mode mixity
The mode mixity, χ , is defined as the ratio between different modes acting of loading. In linear elastic fracture mechanics, the mode mixity, χ el , is a function of the ratio of the mode I and mode II stress intensity factors (Hutchinson and Suo 1991), K I and K II , respectively. Following Chambers et al. (1992) and Shih (1974), we express the mode mixity for an elastic material as where 0 ≤ χ el ≤ 1, with 0 corresponding to pure mode II and 1 to pure mode I. In nonlinear problems, the mode mixity, χ pl , is defined by the following ratio between the shear and normal components of stress acting on a plane directly ahead of the crack tip (Shih 1974) where r and θ are the polar coordinates of a coordinate system that is centred at the crack tip and 0 ≤ χ pl ≤ 1.
Note also that the latter definition is equivalent to the definition in Eq. (38) for linear elastic materials. For the current problem, where a damage zone forms ahead of the crack, the definition of Eq. (39) can be difficult to employ in practice. Also, in the analysis presented above we do not need to determine the detailed form of the stress field ahead of the crack tip. It proves more appropriate to express the mode mixity for the DCB in terms of the applied loading. We employ the definition: whereM = M/M max ,N = N /N max , and M max and N max denote the maximum values of the normal force and bending moment in pure mode I and pure mode II, respectively, for a given value of C * . The maximum values M max and N max can be determined by solving Eq. (35) for a given value of C * in pure mode I and mode II, respectively. In this study we primarily use the definition of Eq. (40), but also compare it with Eq. (39) where appropriate.

The analysis of a steadily propagating crack
Under constant combined load, the crack velocity will, after an initial transient, achieve a steady state, during which it grows at a constant rateȧ. We consider a coordinate system that moves with the crack tip. A material element such as P in Fig. 8a then moves along the x 1 -direction at a rate (Cocks and Ashby 1980) There are two characteristic length scales in this problem, the geometric length scale λ and the material length scaleδ 0 /ε 0 . We can therefore write the effective separation rate in the forṁ where Λ k is a dimensionless function that depends on the interface model, whose detailed form also depends on the ratio of geometric and material length scales φ 0 , the degree of constraining at the crack tip and the mode mixity, as shown in Fig. 8b. For each of the models described in Sect. 3, failure of an element occurs when the effective separation across the damage zone reaches a critical value δ f e , which is not necessarily a material parameter. Integrating the effective separation rate as an element is convected towards the crack tip gives where we have used Eq. (41) to substitute for dt and Eq. (42) to substitute forδ e . The dimensionless function C k is only a function of φ 0 , n, m, χ , degree of constraint and the detailed form of model for the damage zone (k indicates the model, i.e. s ≡ simple, kl ≡ Kachanov linear, ke ≡ Kachanov exponential (a) (b) Fig. 8 The schematics of a steadily propagating crack in viscous solid: a a material point P at distance x 1 from the moving crack tip and b the definition of the Λ k and C k dimensionless functions for point P and m ≡ micromechanical models). It should be mentioned that C k depends on the fracture mode coupling parameters, i.e. α t , α n and l/ h 0 . Note also this function reduces toĈ k employed in Elmukashfi and Cocks (2017) for the case of pure mode I. Substituting forδ m e using Eq. (37) gives the steady state crack velocity whereδ f e = δ f e /λ. We can express this relationship in a number of different forms. An alternative form that can be used to provide some insight into the material response is: where A =δ 0 /σ m 0 is a material constant for the damage zone (see Eq. (8)). The form of this equation might suggest that the crack growth rate is a function of C * , for a given value of δ f e . This is only true if C k is only a function of n and m, or φ 0 and χ are constant for the range of conditions of interest. Note that for the DCB specimen of Fig. 7 and n = m where B =ε 0 /σ m 0 is a material property. Thus, for constant value of λ, the crack velocityȧ is proportional to C * m m+1 . Generally, λ can change as a crack grows, which results in a change of φ 0 , and therefore, a change in the crack growth rate for a given value of C * . This is discussed more fully in the context of mode I loading by Elmukashfi and Cocks (2017).
In order to determine the crack growth rate, we need to evaluate the quantity C k . We can determine this using the Finite Element Method, which can be obtained by equating the numerically determined crack velocity with the prediction of Eq. (45). Further, in order to evaluate the dependence of C k on the mode mixity, the mixity parameter χ pl of Eq. (39) will be determined numerically for a given loading condition in the absence of the damage zone.

Numerical implementation of the governing equations
The initial-boundary value problem described in Sect. 4 is numerically solved using the FE (Finite Element) code Abaqus (Abaqus 2016). A nonlinear quasi-static analysis is used for the initial loading, and a nonlinear visco analysis is used for the creep crack propagation analysis. In the visco analysis implicit time integration is used to solve the FE equations and mixed implicit/explicit integration is used for the integration of the creep and damage zone equations. The FE analysis requires the solution for an elastic/creep constitutive law in the bulk and an elastic rate dependent opening model for the damage zone. Elastic constitutive components have been added to the constitutive relationships of Eqs. (1), (8), (15), (20) and (27), with the values of the elastic components chosen to have limited influence on the computed results. It is worth mentioning that the validity of the proposed model in the presence of elastic deformation in the bulk is discussed in Appendix C and the elastic properties of the interface are selected such that any artificial compliance of the interface is prevented. The geometry of the double cantilever beam (DCB) specimen shown in Fig. 7 is discretised, and a typical Finite Element mesh is shown in Fig. 9. The dimensions are taken as L = 100 mm, H = 10 mm, and B = 1 mm. The initial crack is assumed to be a 0 = 40 mm, and crack propagation is studied over a length Δa = 50 mm. The 4-node reduced integration bilinear plane stress and strain elements (CPS4R and CPE4R) are used in the discretisation for plane stress and strain conditions, respectively. A 4-node twodimensional linear damage zone element was implemented in Abaqus using the user-defined subroutine UEL. The details of the Finite Element implementation are provided in Appendix B. The Finite Element model is divided into two regions, in which the bulk and interface elements are defined. The interface elements are inserted along the crack propagation path, i.e. a ≤ x 1 ≤ L and x 2 = 0. The top and bottom faces of the damage zone elements are attached to the bulk elements, see Fig. 9b. The interface elements are modelled with zero initial thickness such that the nodes on the top and bottom faces coincide. The mesh has 30215 elements, of which 29748 are bulk elements and 467 are damage zone elements. A uniform refined element region is created adjacent to the crack and its propagation for controlling the interface element length l inte . A mesh convergence study was performed for different values of φ 0 = 1 and interface parametersδ f e = 0.004, β = 1.0, f 0 = 0.01, f c = 0.2 and h 0 = l = 0.001 mm. We found that an interface element length of l inte = 0.02 mm is necessary to obtain converged solutions for the range φ 0 = 10 −10 − 10 5 . The interface stiffness parameters K t = K n = 10 6 MPa were selected such that the elastic deformation is negligible (see Appendix C).
The numerical analysis was performed for different combinations of the dimensionless parameters defined above. The damage parameter (ω or f ) is computed and recorded at each Gauss point of the interface during the analysis. The crack tip position, x tip , is defined by ω = 1 or f = f c , and the crack tip velocity is determined where indices p and p + 1 denote variable values at instants t p and t p+1 , respectively, and Δt p = t p+1 − t p is the time increment. Further, the steady crack velocity, a, is computed by taking the average velocity over the steady propagation period.

The mode mixity analysis
It proves instructive to concentrate initially on the definition of mode mixity for different loading combinations. In order to study the effect of mode mixity on creep crack growth, we analyse the DCB specimen for different combinations of mode mixity parameter χ and C * . Different combinations of the loading parameters M and N , ranging between pure mode I and pure mode II, are obtained by solving Eq. (35) at constant C * . Fig. 10 shows the relation betweenM andN defined following Eq. (40) for different values of n and a constant value of C * . These relations are similar for both plane stress and plane strain conditions. It should be noted that for the same value of C * , the maximum values of the bending moment M max and maximum shear force N max in the plane stress and strain conditions are related by M pd max = η ps /η pd M ps max and N pd max = η ps /η pd N ps max , where sub/superscripts ps and pd represent plane stress and strain conditions, respectively. In other words, for the same values of M andN , the values of C * in plane stress and strain conditions are related by C * pd = η ps /η pd n+1 C * ps . The results show that surfaces of constant C * for different values of n in load space nest inside each other, with the surface for n = 1 forming the outer surface and n → ∞ forming the inner surface. (Note that n = 1 corresponds to linear viscous material and n → ∞ corresponds to an ideal plastic material). This type of nesting behaviour is equivalent to the nesting surfaces of constant energy dissipation rate originally identified by Calladine and Drucker (1962) for structural problems and employed by Cocks (1994) in the formulation of constitutive models for porous bodies and sintering. The relationship between the mode mixity ratios χ in Eq. (40) and χ pl defined in Eq. (39) are obtained from the analytical and FE calculations for different values of n. A mesh similar to that in Fig. 10, but without interface elements, and with several rings of 120 elements around the crack tip giving paths of Gauss points of constant radius, is used to determine the stress field and C * . The C * values are taken as the steady state value of the contour integral C(t) which is achieved when the changes are negligible over time (i.e. less than 1%). The value of C * is obtained for 10 different contours separated by a radial distance l cont = 0.01 mm. The value of C * is found to vary with the distance from the crack tip and approach a constant value from the 5th contour, i.e. r = 5 l cont . We therefore take C * as the value determined at the 5th contour. The values of C * for different mode mixities are compared with the analytical solution given by Eq. (35) and they are found to be in excellent agreement. The stress field is obtained at elemental Gauss points at different polar coordinates r and θ . The stress distribution is obtained at a distance r = 5 l cont and 6 l cont from the tip which is found to be in excellent agreement with Westergaard's solutions (Westergaard 1939) in the elastic regime, and with Chamber's and Webster's (Chambers et al. 1992) and Shih's (Shih 1974) solutions in the creep regime.
The mixity ratio, χ pl , is then determined using the definition in Eq. (39) and the stresses at r = 5 l cont for different combinations of the dimensionless loading parametersM andN , i.e. different values of χ . Figure 10a and b show the relationship between the mode mixity ratios χ and χ pl for the plane stress and strain conditions, respectively, and different values of n. In the plane stress case, the values of χ pl are larger than χ for the lower values of χ and greater than χ for the higher values. In the intermediate range of χ , χ pl takes a constant value where the range and the value depend on n. Further, χ pl increases with increasing of n for small values of χ and decreases with increasing value of n for higher values of χ . In plane strain, χ pl are larger than χ for the entire range and it increases with increasing n. The plane stress results appear to be different from the results reported by Chambers et al. (1992). In plane strain, the results are in good agreement with the results reported by Shih (1974). The difference in plane stress might be attributed to the in-plane constraint that is introduced by the far-field loading. Hence, when the in-plane constraint is significant, the higher order terms of the small scale yielding stress field become large, which leads to a change in the near-field mode mixity, i.e. the constraint changes the values of the stresses in Eq. (39). The constraint effects in plane strain appear to be minimal. Several approaches have been used to characterise the in-plane constraint levels such as the J − Q approach by O'Dowd and Shih Shih 1991, 1992) and the J − A 2 approach by Yang et al. (1993), Chao et al. (1994) in elastic-plastic materials. Using the analogy between J and C * , Budden and Ainsworth (1999) and Bettinson et al. (2001) have used C * − Q to characterise the crack tip fields under creep conditions, i.e. the stress field is σ i j = σ HRR i j is the HRR stress field (Hutchinson 1968;Rice and Rosengren 1968) and Q is a geometry dependent parameter. Several studies have shown that, for a number of different specimens, the Q parameter can take a negative value and it decreases with increase of outof-plane constraint (O'dowd and Shih 1991;Cravero and Ruggieri 2003). Thus, for Q < 0, the tangential stress σ θθ decreases and therefore reduces the nearfield mixity parameter χ pl as indicated in Eq. (39). The constraint effects are not investigated further in this study.
In the subsequent sections, we use the far-field loading mode mixity χ for the given DCB specimen under combined bending and shear force loading. These results can be converted to a dependence on χ pl by employing the relationships plotted in Fig. 11. (a) (b) Fig. 11 The relation between the local mixity ratio χ pl and the far-field loading mixity ratio χ for: a plane stress conditions; and b plane strain conditions

Crack growth
We have performed several analyses for different combinations of the dimensionless parameters, mode mixity and interface properties. In each analysis, the crack tip position during transient and steady state growth have been obtained. To illustrate the different results, the simple interface model, with the representative separation taken to be the effective separation, i.e.δ = δ e , is adopted. Further, the parameters n = m = 9, φ 0 = 1.0, α t = α n = 1.0 andδ f e = 0.004 are used, and the model is analysed assuming plane stress conditions. Figure 12a and b show the crack tip position and velocity as functions of time, respectively. The plots show that the crack starts to propagate slowly and accelerates to a high velocity before dropping after a short time (∼ 20 hr) to a lower steady state velocity. In this case a steady velocity ofȧ = 0.103 mm/hr is obtained. The result shows that the transient velocity is higher than the steady state velocity suggesting that the stress at the crack tip is initially high due to the elastic deformation. As the crack advances creep deformation eventually dominates and the stress fully relaxes leading to a slower propagation rate. Additionally, the damage zone ahead of the crack tip is fully developed during both transient and steady propagation. The other scenario is when the damage is not fully developed during the transient stage, which may lead to a slower rate of propagation before reaching a steady state.

The C k -function and mode mixity
In order to evaluate C k , the Finite Element analysis is used to determine the steady-state crack velocityȧ. Thus, for a given set of input parameters and mode mixity, C k can be evaluated from Eq. (44): The appropriateness of the dimensionless analysis has been examined using different values of material and geometric parameters, while keeping the dimensionless quantities and mode mixity constant. The results are found to be independent of the choice of these parameters. Therefore, the dimensionless groups identified within the mathematical formulation are appropriate for the class of problem analysed. Further, the physical limits, the limits of small and large φ 0 , and the validity of the above approach are discussed in Appendix C. It proves instructive to investigate initially the form of the relation between C k and φ 0 for different values of mode mixity χ . We explore this by considering the simple interface model with α t = α n = 1.0 and rate sensitivity parameters n = m = 9. Further, we limit the investigation to the case of pure mode II loading, i.e. χ = 1. Similar to pure mode I studied in Elmukashfi and Cocks (2017), the relationship between C k and φ 0 is found to follow two separate power-law relations. Thus, the relationship takes the general form where the parameters γ k and η k depend on the range of φ 0 considered, the mode mixity χ and the deformation constraint. Figure 13a and b show the relationship between C s (i.e. recall that k = s for the simple model) and φ 0 for plane stress and strain conditions, respectively. In each case, the relationship is given for modes I and II. In the case of plane stress (Fig. 13a), for small values of φ 0 ∈ [10 −10 , 8 × 10 −2 ], γ s = 0.45 and 0.23, η s = −0.06 and −0.10, for modes I and II, respectively. The results imply that γ k and η k depend on the mode mixity such that both parameters increase with increasing χ . For larger values of φ 0 ∈ [8 × 10 −2 , 10 3 ], γ s = 0.09 and 0.07, and η s = −0.67 and −0.25, for modes I and II, respectively. Hence, γ k and η k depend strongly on the mode mixity and increase with increasing value of χ . The transition between the two power-laws approximately occurs over the range 10 −4 ≤ φ 0 ≤ 10 0 for modes I and II, which implies that the value of φ 0 at which the transition takes place does not depend on the mode mixity. C s for modes I and II is bounded by the physical limits of Eqs. (C.13) and (C.14), i.e. the stiff and compliant limits, respectively, and asymptotes to Eq. (C.13) at large values of φ 0 , see Appendix C. Further, for the adopted value of σ 0 /E = 8×10 −6 in the analysis, elastic effects can be ignored-the chain lines in Fig. 13a illustrate the upper limit to the validity of C * for different values ofδ f e (i.e. Eq. (C.16)). Now, we concentrate on the regime φ 0 ≥ 1 which is more representative of the physical behaviour of engineering materials, i.e. the characteristic material length scale is expected to be greater than the geometric length scale. For pure mode II, for φ 0 ≥ 1, the variation of C s with φ 0 is slightly shallower than the trend in the stiff limit and meets this limit at φ 0 ≥ 10 2 , see Fig. 13a, which is smaller than observed for mode I, i.e. φ 0 ≥ 10 6 . Up to these limits, the crack propagation rate in Eq. (45) can be determined using the definition in Eq. (49) aṡ This relationship shows that for a given damage zone and constant values of the geometric length scale λ and C * , the crack propagation rate increases with the increasing mode mixity, i.e. γ s and φ η s 0 increase with increasing value of χ . Note that for general crack geometries, λ changes with mode mixity and crack length, which needs to be taken into account in the crack growth process. The dependence ofȧ on χ is determined by the form of γ s and the change of λ during crack growth. Further discussion is provided below.
In the case of plane strain (Fig. 13b), for small values of φ 0 ∈ [10 −10 , ×10 −2 ], γ s = 0.45 and 0.28, and η s = −0.06 and −0.09, for modes I and II, respectively. For the larger values of φ 0 ∈ [10 −2 , ×10 −4 ], γ s = 0.19 and 0.10, and η s = −0.29, for modes I and II, respectively. Therefore, γ s and η s increase with increasing value of χ for small values of φ 0 , and only γ s increases when φ 0 > 10 −2 , wherein the parameter η s is independent of the mode mixity. Hence, as determined in Sect. 5.1, the constraint appears to control η s in the case of plane stress. The transition between the two power-laws is over 10 −4 ≤ φ 0 ≤ 10 0 for both modes I and II. As in the case of plane stress, the values of C s for modes I and II are bounded by the physical limits and lie in a regime where elastic effects can be ignored, see Appendix C. For φ 0 > 1.0 the crack propagation rate in plane strain is given bẏ The above relations imply that a faster crack growth rate occurs in plane strain than plane stress for given values of C * and φ 0 , due to the out-of-plane constraint, i.e. higher stress levels ahead of a plane strain crack. The difference between the C s values for modes I and II is larger in plane strain than plane stress, i.e. γ s for mode II is 22% lower than the mode I value for plane stress, while it is 47% lower in the case of plane strain. In order to investigate the effect of the mode mixity in detail, we determined C s (i.e. γ s and η s ) for different values of mixity. Figure 14a and b shows the relationship between γ s , η s and χ for plane stress and strain conditions. In plane stress, for smaller values of φ 0 < 1, γ s increases with increasing mode mixity χ taking a minimum value of 0.23 for mode II to a maximum value of 1.17 for mode I. γ s remains almost constant over the range 0.1 ≤ χ ≤ 0.8, due to the almost constant near-(a) (b) Fig. 13 The relation between C s -function and φ 0 parameter and the physical limits in the case of n = m = 9, α t = α n = 1.0 and mode I and II: a plane stress conditions; b plane strain conditions. The red and blue lines represent the compliant and stiff limits, respectively, the green dash-dot lines represent the C * validity limit for different dimensionless separation at failureδ f e , and the dashed lines show the power law fit field mixity χ pl of Eq. (39) over this range, as shown in Fig.14a. Similar results are obtained for larger values of φ 0 ≥ 1, with minimum and maximum values of 0.065 and 0.17 for mode II and I, respectively, and constant values over the range 0.3 ≤ χ ≤ 0.8. For smaller values of φ 0 < 1, η s remains almost constant with a slight increase from −0.09 for mode II to −0.06 for mode I. For larger values of φ 0 ≥ 1, η s decreases with increasing mode mixity χ taking a maximum value of −0.25 for mode II to a minimum value of −0.67 for mode I which can be attributed to the constraint effects. Similar results are obtained in the case of plane strain. For smaller values of φ 0 < 1, γ s has a minimum value of 0.27 for pure mode II and a maximum value of 2.13 for pure mode I and it remains almost constant over the range 0.1 ≤ χ ≤ 0.8. Similarly, for larger values of φ 0 ≥ 1, minimum and maximum values of 0.027 and 0.7 for mode II and I, respectively, and constant values over the range 0.1 ≤ χ ≤ 0.8 are obtained. η s remains almost constant for smaller values of φ 0 < 1 with a slight increase from −0.09 for mode II to −0.05 for mode I which is similar to the plane stress case. For larger values of φ 0 ≥ 1, η s remains almost constant taking a value −0.28 indicating that the constraint effects are minimal.

The role of coupling on deformation and damage development
There are numerous ways to define the contributions of normal and tangential loading to deformation and damage development at the interface, as discussed earlier in Sect. 3. C * defines the stress field in the vicinity of the crack tip, which determines the deformation in the interface, although the form of the field is different for different modes of loading. For a given mode, damage development and deformation can be directly related to C * . In mixed mode loading, the remote C * and the mode mixity χ define the singular stress field ahead of the crack tip and the near field mode mixity. Thus, the deformation in the interface is determined by the coupling of the tangential and normal deformations as well as the near field mixity. Damage development due to this coupling can be expressed in terms of a damage controlling parameter, such as the representative separationδ. Hence, the contributions of the fracture modes to deformation and damage development in the interface are determined by C * , the mode mixity χ and coupling of the deformation and damage controlling parameterδ. For the phenomenological simple and Kachanov type models of Sect. 3, the contributions to the deformation of the interface from the different modes of loading are determined by the parameters α t and α n whereas their contributions to the damage development is controlled by the representative separation δ. In the micromechanical model, the contributions to deformation and C * are determined by the functions g t and g n that only depend on the pore area fraction f , and their contributions to damage development is controlled by the ratio l/ h 0 . In this section, we investigate the relationship between the remote C * , the deformation of the interface and the choice of the damage controlling parameter at different mode mixities χ . We consider the simple interface model with rate sensitivity parameters n = m = 9 and assume plane stress and strain condi-tions. The aim is to study the relationship between C s of Eq. (48) and the mode mixity χ for different combinations of the parameters α t and α n , and forms ofδ. The contributions of the normal and tangential modes of loading to the deformation of the interface is explored by studying α t /α n in the range α t /α n ∈ [0.5 − 1.5]. The interface response in Eqs. (15), (20) and (27) suggests that for given interface traction, the tangential deformation decreases with increase of the ratio α t /α n . The measure of separation that controls damage development is considered to beδ = δ e orδ = δ n . For δ = δ e , the effect of the ratio α t /α n on interface damage is similar to its effect on deformation. In the case of δ = δ n , where only normal opening controls the damaging process, the interface response remains dependent on the parameters α t and α n , which determine the deformation in the normal direction and therefore the damage development. It is worth mentioning that, when local failure is determined byδ = δ n , the crack growth law of Eqs. (44) and (45) requires δ f e to be determined. This is determined numerically using the FE analysis, i.e. δ f e = t f 0δ e dt, where t f is the time to failure for a material point along the interface under steady crack propagation conditions, i.e. whenδ = δ f n . The quantity γ s has been obtained for different values of the ratio α t /α n and mode mixity, and for the different forms of representative separationδ. It should be noted that η s does not depend on the ratio α t /α n and its values for different mode mixity are reported in Sect. 5.3. Figure 15a and b show the relationship between γ s and χ for the two measures ofδ under plane stress conditions, respectively. Forδ = δ e , Fig. 15a, C s takes two different forms for α t /α n < 1 and α t /α n ≥ 1. The result shows that increasing the ratio α t /α n above 1 reduces C s , with the magnitude of the reduction decreasing with increasing mode mixity χ . Further, reducing the ratio below 1 increases C s , with the magnitude of the increase decreasing with increasing mode mixity. The reduction in C s in the case of α t /α n ≥ 1 can be attributed to the decline of the role of the tangential deformation. For α t /α n < 1, C s takes a large value at χ = 0 and decreases with increasing mode mixity, reaching a plateau as χ increases further. It experiences a slight decrease as mode I (χ = 1) is approached. It is worth mentioning that in the cases of α t /α n < 0.75, the maximum value of C s occurs for mode II (χ = 0). Hence, the increase of C s is due to the increase of the role of tangential deformation. Moreover, the maxi-(a) (b) Fig. 15 The relationship between C s and the mode mixity χ for the case of φ 0 = 1, n = m = 9, plane stress conditions and for different representative separation: aδ = δ e ; bδ = δ n . Here the parameter γ s = C s as φ 0 = 1 mum and minimum values of C s and the range over which it is constant depend on the ratio α t /α n . Now we consider the case ofδ = δ n , shown in Fig. 15b. The results show that γ s is zero for pure mode II, as expected since in this limitδ = δ n = 0, and increases with increasing value of mode mixity χ . Similar forms of the C s -function to the caseδ = δ e are obtained for χ ∈ [0 − 1] and α t /α n ≥ 1. Comparing the distributions of C s forδ = δ e andδ = δ n in Fig. 15a and b, respectively, demonstrates that C s is always lower forδ = δ n , which can be explained by the fact that, in addition to the dependency of the deformation on the ratio α t /α n , only normal deformation in the interface contributes to the damage development, which means that only part of C * contributes to damage development and crack growth. Recall the crack growth rate in Eq. (47) and consider the different forms of γ s = C s in Fig. 15a and b for the different values of α t /α n andδ. For constant values of the geometric length scale λ, δ f n and C * , and for α t /α n > 1, the crack growth rate decreases with increasing α t /α n as the part of C * that contributes to crack growth gradually reduces. For α t /α n < 1, a faster growth rate is obtained, which is due to a larger proportion of C * contributing to crack growth.
Similar results are obtained in the case of plane strain, as illustrated in Fig. 16a and b, with two different relations obtained for γ s for α t /α n < 1 and α t /α n ≥ 1. Increasing the ratio α t /α n beyond 1 reduces γ s , while decreasing it below 1 increases γ s . The magnitude of the change (i.e. decrease or increase) decreases with increasing mode mixity χ . Comparing γ s for the caseŝ δ = δ e andδ = δ n , indicates that it is zero for pure mode II loading whenδ = δ n , and then increases with increasing mode mixity χ , approaching a similar form toδ = δ e for large values of χ . As in the case of plane stress, the dependency of γ s on the ratio α t /α n can be explained by the change in the contribution of the tangential deformation in the case ofδ = δ e , whereas the additional reduction forδ = δ n is due to limiting the damage development to the normal component of displacement. It is worth mentioning, that similar conclusions about the crack growth rate in plane strain can be drawn using the crack growth rate in Eq. (51) and the different forms of γ s in Fig. 16a and b. The effect of the coupling between deformation and damage development is similar in plane stress and plane strain.
By selecting α n = 1 for these simulations and varying α t results in the plots converging on to a single point in Figs. 15 and 16 at χ = 1. We could have alternatively set α t = 1 and varied α n to cover the same range of ratios as in these figures. Then all the curves would converge at the point χ = 0. By exploring this limit it can be shown that the points of intersection of the curves with the γ s axis in Figs. 15a and 16a at χ = 0 can be determined directly from the mode II results presented in Fig. 13, such that γ s (0) = C s α m+1 (a) (b) Fig. 16 The relationship between C s and the mode mixity χ for the case of φ 0 = 1, n = m = 9, plane stress conditions and for different representative separation: aδ = δ e ; bδ = δ n . Here the parameter γ s = C s as φ 0 = 1 C s corresponding to φ 0 = α m+1 t . This result applies in both plane stress and strain.

The effect of damage model
In this section, we investigate the detailed form of the damage zone constitutive model on crack growth for different values of mode mixity χ , and for different extents of coupling between tangential and nor-mal deformation, α t /α n . Previously in Elmukashfi and Cocks (2017), we examined the different damage models described in Sect. 3 under prescribed uniaxial stress and for the case of pure mode I crack growth. The results show that there is a point-wise linear relationship between C k for the Kachanov and micromechanical models and its value for the simple model, i.e. C k = μ k C s with μ k ≥ 1. In the stiff limit, the detailed form of the damage zone is found to be important and the largest values of C k are obtained for the micromechanical model. Further, in the case of the Kachanov models, the exponential model yields larger values than the linear model.
We consider the damage models presented in Sect. 3 and determine C k in Eq. (48) (i.e. γ s and η s ) for different mode mixities. We limit our consideration to the response close to the stiff limit, which is of more practical significance. Further, the representative separation is taken to beδ = δ e in the case of Kachanov type models and plane stress conditions are assumed to prevail. We investigate the role of the coupling between the tangential and normal deformation in the interface by investigating the effect of the ratio α t /α n in the case of the Kachanov type models and the ratio between the pores' spacing and their initial height l/ h 0 for the micromechanical model. The parameters employed for these models are n = m = 9, δ f e = 0.02 mm, β = 1, h 0 = 0.02 mm, f 0 = 0.01, f c = 0.5 and g 0 = 2.23. Further, the ratio α t /α n is taken to be in the range α t /α n ∈ [0.5 − 1.5] with α n = 1 and the ratio l/ h 0 taken to be in the range l/ h 0 ∈ [0.2 − 10] where l is changed to give the specified ratio. It is worth mentioning that the models yield the same rapture time under a prescribed uniaxial constant stress for the given parameters.
It proves instructive to express γ k as a function of the ratios α t /α n or l/ h 0 for different mode mixity χ . Figs. 17, 18 and 19 show the relationship between γ k and α t /α n for the linear and exponential Kachanov models and the micromechanical model, respectively, for different values of χ . For the linear Kachanov model in Fig. 17, for α t /α n < 1, C kl decreases with the increasing mode mixity χ for χ ≤ 0.8, and then increases as the mixity approaches mode I (χ = 1), with the magnitude of this increase depending on the ratio α t /α n , i.e. the mode I value of C kl is larger than that for mode II for α t /α n > 0.75. For α t /α n ≥ 1, C kl increases with mode mixity χ and takes almost a constant value in the range χ ∈ [0.4 − 0.6] with the Fig. 17 The relationship between C kl and the mode mixity χ for the case of φ 0 = 1, n = m = 9, plane stress conditions. Here the parameter γ kl = C kl as φ 0 = 1 Fig. 18 The relationship between C ke and the mode mixity χ for the case of φ 0 = 1, n = m = 9, plane stress conditions. Here the parameter γ ke = C ke as φ 0 = 1 maximum value occurring for mode I. For α t /α n ≥ 1, γ kl increases with increasing mode mixity χ such that at constant α t /α n the minimum value of γ kl occurs at χ = 0 and the maximum is at χ = 1. Hence, the results imply that crack growth is controlled by the mode mixity and the coupling between the normal and tangential deformation, as in the case of the simple model. Fig. 19 The relationship between C m and the mode mixity χ for the case of φ 0 = 1, n = m = 9, plane stress conditions. Here the parameter γ m = C m as φ 0 = 1 For the exponential Kachanov model, Fig. 18, for α t /α n < 1, γ ke shows similar behaviour to the linear model. γ ke is larger in mode I than mode II for α t /α n > 0.62. For α t /α n ≥ 1, γ ke decreases with increasing mode mixity and experiences a small increase at χ = 0.4 that is followed by a decrease to its lowest value at χ = 0.6. It increases again for χ > 0.6 and reaches its maximum value for pure mode I. For the micromechanical model, Fig. 19, for l/ h 0 < 4.9, C m increases with increasing mode mixity χ for χ < 0.4 and χ > 0.8, and decreases with increasing χ for 0.4 ≤ χ ≤ 0.8. For l/ h 0 ≥ 4.9, γ m initially increases for χ ≤ 0.2, followed by a decrease with increasing mode mixity χ for χ ≤ 0.8. An additional increase takes place for χ > 0.8. The maximum and minimum values of γ m occur at χ = 0.2 for l/ h 0 = 0.2 and 10, respectively. Thus, the results imply that, for l/ h 0 < 4.9, the tangential deformation dominates the damage development process and the most severe combination of the normal and the tangential deformation is achieved in the case of χ = 0.2, as discussed earlier in Sect. 3, which results in a reduction of γ m with an increase of mode mixity χ . For l/ h 0 ≥ 4.9, the role of the normal deformation increases, and therefore C m decreases significantly in the mode II limit. Further, for higher values of l/ h 0 , γ m increases with mode mixity χ . It is worth mentioning that the value We now consider the results of simulations for other values of φ 0 , which allows us to determine the parameter η k in Eq. (49). We consider the situations where α t = α n = 1 for the Kachanov damage models and h 0 = l = 0.02 mm for the micromechanical model. Fig. 20 shows the relationship between η k and χ for the different damage models for plane stress. For the linear Kachanov model, η kl decreases with increasing mode mixity χ for χ ≤ 0.8, reaching a minimum value of −0.62 at χ ≈ 0.8, followed by a slight increase, reaching a value of −0.43 for pure mode I. η kl takes a maximum value of −0.22 for pure mode II. Similar behaviour is obtained for the exponential Kachanov model with a minimum value of −0.66 at χ ≈ 0.8, and values of −0.52 and −0.22 for modes I and II, respectively. For the micromechanical model, η m remains almost constant over the range χ ≤ 0.7 taking the value −0.13. The plateau region is followed by a steep decrease, reaching the minimum value of −0.37 for pure mode I.

Comparison with experimental data
The analysis of the idealised DCB specimen presented above is rigorous due to C * not changing as the crack grows, allowing a true steady state to be achieved. In conventional components, whether operating in an industrial facility or tested in the laboratory, C * continually changes as a crack grows. We need to make some assumptions in order to extend the results presented above to realistic geometries. Here we assume that as a crack grows a state equivalent to the steady state described above is achieved rapidly, so that the crack velocity only depends on the current value of C * , φ 0 , etc., and not on the history of loading and crack growth. This is consistent with assumptions generally employed in the interpretation of experimental data.
In this section, we calibrate the damage models against available experimental data. In particular, we extend the fitting procedure proposed in Elmukashfi and Cocks (2017) for pure mode I to II and mixed mode cases. The main focus is to explore the relationship between damage model and the crack growth process, including the role of different characteristic material and geometric length scales, coupling of deformation in the damage zone and mode mixity. In the process we explore where the available data lies with respect to the general trends identified in Figs. 17, 18, 19 and 20. We consider the high chromium steel (Jethete M 152 steel at 550 • C) investigated by Hyde and Chambers (1991). They provide data for creep deformation and creep rupture, as well as creep crack growth for different mode mixity conditions generated using compact tension (CT) and compact mixed-mode (CMM) specimens, see Figs. 21 and 22. In the case of pure mode I, the crack growth direction was found to be predominantly straight ahead, i.e. θ = 0 • to the initial crack with a significant thinning of the specimen near the crack tip indicating large creep strains. In the case of pure mode II, the crack propagates straight ahead parallel to the initial crack, i.e. similar to pure mode I θ = 0 • , such that the displacement of the surface markings near the crack tip shows a clear evidence of the shear deformation. In the mixed mode case, the crack propagated in two directions that are θ = 0 • and θ = −90 • to the initial crack, and the latter crack was considered to be dominant because it generally led to final failure. Therefore, in this study, we limit our consideration to the cases of pure mode I and II. In order to determine the crack growth rate, we need to determine the characteristic length λ for the cracked geometries under different mode mixity χ and the material parameters, i.e. n andε 0 (at a reference stress σ 0 ) or equivalently B and the interface model parameters. In this paper we limit the fitting to the linear Kachanov model, and we assume that damage development is determined byδ = δ e . We also assume isotropic damage, i.e. ω t = ω n = ω (δ e ). Parallel procedures can be undertaken for the other models described in this paper. We need to determine n, m, δ f e andδ 0 (at a reference stress σ 0 ) or equivalently the material parameter A and the mode mixity related parameters α t and α n . For simplicity, we take α n = 1 and consider α t to be a fitting parameter. (This assumption is consistent with the mode I fitting procedure in Elmukashfi and Cocks (2017).) The steady creep response under a constant uniaxial stress σ is given by Hyde and Chambers (1991) at 550 • C: where n = 9.5 and B = 2.93 × 10 −27 MPa −9.5 h −1 . Similarly, we assume that the damage zone model can also be used to describe the damage development on grain boundaries in a uniaxial test and the damage grows primarily on the boundaries normal to the direction of the applied loading. The integration of the damage growth rate equation between the limits, ω = 0 at time t = 0 and ω = 1 at failure, i.e. when t = t f , yields Creep rupture data given by Hyde and Chambers (1991) is plotted in Fig. 21, which gives m ≈ 9 and D = 1.41 × 10 25 MPa 9.5 h. Therefore, the relationship between the material parameters is where A is measured in units of mm/(MPa 9 · h). In the following steps, we aim to determine δ f e , A and α t by fitting the creep crack growth data (ȧ vs C * ) for modes I and II to the model in Eq. (45). Firstly, we determine a relationship between δ f e and A by fitting mode I data which can be used together with Eq. (54) to determine the two parameters. Secondly, knowing δ f e and A, we obtain α t by fitting mode II data. To do this, we need to obtain the geometric length scale for the compact tension (CT) and compact mixed-mode (CMM) specimens for different mode mixities. Hence, following the procedure in Elmukashfi and Cocks (2017), we adopt the definition in the UK R5 assessment procedure (Ainsworth et al. 1987) whereε R is the uniaxial strain rate at a reference stress σ R . The reference stress is defined by where P L is the limit load for a perfectly plastic material of yield strength σ y and P is the applied load. The characteristic length scale for a given specimen and mode mixity is determined from the stress intensity factor and reference stress by λ = (K /σ 0 ) 2 where K is the stress intensity factor at the applied load P, which is taken as K = K 2 I + K 2 II . The limit load and stress intensity factor for the compact tension (CT) and compact mixed-mode (CMM) specimens can generally be expressed as where a is the crack length, W is the specimen width, Υ L and Υ K are the limit load and stress intensity factor shape functions, respectively, that depend on the specimen and mode mixity. Thus, Eqs. (56) and (57) and the definition of λ, the characteristic length can be written in the general form where Υ = (Υ L Υ K ) 2 . Now consider mode I, we follow the procedure in Elmukashfi and Cocks (2017) and use Miller and Ainsworth's (Miller and Ainsworth 1989) is a dimensionless parameter (load factor); this modification gives better agreement with computational results, which effectively reduce the reference stress by a factor of F p . For the CT specimen tested by Hyde and Chambers (1991), W = 32 mm and the initial crack length is a 0 = 16 mm. Thus, for n = 9.5 and a/W in the range 0.25 to 0.5, the load parameter is in the range F p ∈ [0.92 − 0.96]. In this investigation we use the upper bound value of 0.96. The characteristic length scale for the compact tension specimen is determined by Elmukashfi and Cocks (2017), wherein the shape function in Eq. (58) is determined from the limit load and stress intensity factor shape functions for the CT specimen: where κ = 1.155. The relation between the crack growth rateȧ and C * in Fig. 22 suggests that stable crack growth is obtained over two order of magnitudes of C * , which corresponds to an increase of crack  (59) and (60) to be λ/W = 0.74 to 0.42. In the evaluation of the data, we take the average value of the characteristic length, i.e. λ/W ≈ 0.58. To determine the crack growth rate, we need an expression for C k . Thus, substitution of A from the definition φ 0 = B A λ into Eq. (54) yields φ 0 = 7.67/δ f e where δ f e is in mm. Therefore, for the physical range of the critical separation in metals (i.e. it is of the order of the mean cavity spacing δ f e ∈ 10 −10 − 10 −5 mm), the dimensionless parameter is φ 0 > 1. It follows that we use the relationship C kl = 0.4 φ −0.494 0 for φ 0 > 1 and χ = 1. In this regime, the crack growth rate is determined from Eqs. (49) and (45) as: whereȧ is in mm/h, δ f e is in mm and C * is in mJ/mm2/h. Fitting the expression in Eq. (61) to the mode I data in Fig. 22 yields the critical separation δ f e = 67.83 μm. The rate parameter A can then be determined from Eq. (54) to be A = 4.81 × 10 −28 mm/MPa 9 h and φ 0 = 1.13 × 10 2 . a δ f e is estimated from the comparison between experimental creep crack growth data for Jethete M 152 steel at 550 • C and the proposed framework. b Note, σ 0 changes during the duration of a test, causingε 0 andδ 0 to also change. At a given instant σ 0 can be determined from Eq. (56), multiplied by F p = 0.96, to take into account the correction of Miller and Ainsworth (1989) Now we are interested in obtaining the parameter α t .
We consider the CMM specimen tested by Hyde and Chambers (1991) under pure mode II loading, of width W = 33 mm and initial crack length a 0 = 16.5 mm. It should be noted that the data consists of only three points which may affect the fitting results. However, the comparison provides an important physical insight into the coupling parameters (i.e. α t and l/ h 0 ). For the mode II case, the characteristic length scale is determined using Eqs. (58) where the limit load and stress intensity factor shape functions for a compact mixedmode (CMM) specimen under pure mode II are given by Chambers (1988) and(1991) The relation between the crack growth rateȧ and C * in Fig. 22 suggests that stable crack growth is obtained over one order of magnitude change in C * , which corresponds to an increase of crack length from a/W = 0.50 to 0.61. Similarly, the change of characteristic length over the C * range is estimated from Eqs. (58), (59) and (60) to be λ/W = 0.24 to 0.18. We take the average value of the characteristic length, i.e. λ/W ≈ 0.21. Thus, the dimensionless parameter is φ 0 = 4.11 × 10 2 . Therefore, we use C kl for φ 0 > 1 and χ = 0. Using the definition in Eq. (49) and the fact that γ kl depends on χ and the parameters α t and α n , we can write the form of C kl for n = m = 9 and mode II as C kl = γ kl (χ, α t , α n ) φ −0.223 0 . Thus, the crack growth rate is determined from Eqs. (49) and (45) as: whereȧ is in mm/h and C * is in mJ/mm2/h. Fitting the expression in Eq. (64) to the mode II data in Fig. 17 yields γ kl ≈ 0.062. Thus, using the relationship between C kl and the parameter α t in Fig. 17 for α n = 1.0, gives α t ≈ 1.14.
The same fitting procedure can be employed for the exponential Kachanov and the micromechanical models. In these models additional parameters are required, i.e. β and α t for the exponential model (similarly we assume α n = 1), and h 0 , l, f 0 and f c for the micromechanical model. We limit our consideration to the case of β = 1 for the exponential model. Then, for φ 0 > 1, C ke = γ ke (χ, α t , α n ) φ η ke 0 and C m = γ m (χ, h 0 , l) φ η m 0 with γ ke ≈ 0.03, η ke ≈ −0.515, γ m ≈ 0.49 and η m ≈ −0.373 for the case of pure mode I and η ke ≈ −0.215 and η m ≈ −0.1 for the case of pure mode II. Therefore, using a nonlinear least squares method, the critical separation for both models are determined by fitting the crack growth rate of Eqs. (45) and (49) to mode I creep data. It follows that we obtain the parameters γ ke ≈ 0.031 and γ m ≈ 0.127 for the case of pure mode II. Thus, using the relationship between γ kl and α t , and γ m and l/ h 0 in Figs. 18 and 19, respectively, we obtain α t for the exponential model and h 0 , l, f 0 and f c for the micromechanical model. It should be noted that the model predictions lie on the same straight line and yield comparable goodness of fit of ≈ 93%. The fitting parameters for the different models are given in Table 1. The parameters indicate that the material data falls in the regime φ 0 ∈ 10 1 − 10 3 , i.e. close to the stiff limit and the critical separation falls within the physical regime, i.e. δ f e ∈ 10 −10 − 10 −5 m. In our previous study (Elmukashfi and Cocks 2017), we compared the proposed model with conventional models of creep crack growth, which do not include the effect of the ratio of geometric and material length scales, i.e.ȧ ∝ C * q , where q < 1. We concluded that for large defects, as in the commonly used experimental specimens (e.g. CT, DCB, CMM, etc.), the characteristic length scale λ does not change significantly during crack growth, i.e. by < 13% for CT specimens and by < 5% for CMM specimens. Therefore, the proposed model is then equivalent to conventional models in predicting a power-law dependence between crack growth rate and C * and can be fit to data in a similar way to conventional existing models. In the case of a smaller defects, which are commonly evaluated in assessments of real components, conventional models result in an overestimation of the crack growth rate. Additionally, a reduction of the characteristic length scale λ decreases φ 0 , which may lead to a transition to the low φ 0 regime, resulting in a further overestimation of the crack growth rates. Hence, the proposed framework is useful for the cases where defects are small which are more common in engineering components.
Comparing mode I and II, the experimental data show that for the same value of C * , the crack growth rate in mode II is reduced to ≈ 25% of the growth rate in mode I. Further, the characteristic length scale for the CMM specimen in mode II is ≈ 36% of that for the CT specimen in mode I. Hence, the model predictions give a crack growth rate that is proportional to λ 1+η k where 1 + η k ≈ 0.5, therefore, the change in λ yields ≈ 40% reduction in the crack growth rate in mode II, i.e. 53% of the total 75% reduction in the crack growth rate. The remaining 22% reduction comes from the decrease in the mode mixity constraint and the coupling between the tangential and normal deformation in the interface, which are captured by the parameter γ k .

Concluding remarks
In this study, we have extended the theoretical and computational framework presented by Elmukashfi and Cocks (2017) to the case of mixed mode creep crack growth under steady state creep conditions. We have developed a theoretical framework wherein the bulk material is described by a power-law creep law and the interface behaviour is assumed to be described by a generalised traction-rate of separation law. The damage zone models have been extended to the mixed mode loading case such that the deformation and damage processes are expressed in terms of effective traction and rate of separation which determine the coupling between normal and tangential deformation and a representative separation, respectively. A double cantilever beam specimen (DCB) subjected to a combined force and moment is studied, allowing the path independency of the C * -integral to be used to relate the far field loading to the damage growth process along the interface during steady state crack growth. The damage is assumed to be localised in a narrow zone ahead of the crack tip. Analytical models are developed for mixed mode steady-state crack growth, wherein the crack growth rate is determined by a dimensionless function, C k . A computational framework is then implemented using the Finite Element Method. The analytical models are calibrated using the detailed Finite Element model to give C k for different combinations of the dimensionless parameters and different loading conditions. The quantity C k depends on the ratio of geometric to material length scales φ 0 , the rate sensitivity exponents n and m, mode mixity ratio χ , the form of the damage law, the coupling parameters (α t , α n , h 0 and l) and the out-of-plane constraint. Similar expressions for C k (i.e. Eq. (49)) are obtained in the case of mode II loading which are bounded by the physical limits. In the case of plane stress, the parameters γ s and η s depend on the mode mixity for the entire range of φ 0 . In the case of plane strain, at the stiff limit, i.e. φ 0 → ∞, η s does not depend on the mode mixity whereas γ s depends on the mode mixity. In the compliant limit, i.e. φ 0 → 0, both γ s and η s depend on the mode mixity. C s is found to be strongly dependent on the near-field mode mixity χ pl , given by Eq. (39). Further, the coupling between the tangential and normal directions, expressed in terms of the ratio α t /α n , affects the deformation and damage development at the interface. Damage development is also determined by the definition of the representative separationδ, which achieves a critical value at the growing crack tip. Faster crack growth rates are obtained when damage generated in the interface influences the constitutive response, such as in the Kachanov and micromechanical models. Substantial stress relaxation can then occur along the interface ahead of the crack tip, which changes the nearfield mode mixity and therefore changes C k . Further, the coupling parameters, i.e. α t and α n in the Kachanov models, and h 0 and l in the micromechanical model, affect crack growth in a similar manner. In the final part of this study, we have compared the model with experimental data for a Jethete M 152 steel at 550 • C. The results show that the different damage models are capable of fitting and explaining the creep crack growth data for different mode mixities using physically reasonable parameters.
Acknowledgements The authors are grateful to funding from the EPSRC under grant EP/K007815/1. They would also like to acknowledge helpful discussions with Dr Jianan Hu.
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: Analysis of beam under steady-state creep and combined axial load and bending moment
Consider the beam element in Fig. 23 below. The length, height, and thickness of the beam are denoted by L, H and B, respectively. The beam is subjected to constant axial force and moment per unit thickness N and M, respectively. Both, plane stress and strain conditions, can be considered, such that the x 1 x 2 -plane is the plane in which loading is applied and in which we determine the stresses and strains.
The beam is assumed to exhibit steady-state creep behaviour and is defined by the constitutive law of Fig. 23 The schematic of a beam element subjected to constant axial force and moment per unit thickness N and M, respectively. The element is of thickness B Eq. (1). The deformation in the beam is given by the normal strain in 1-directionε c 11 aṡ ε c 11 =ε −κ x 2 , (A.1) whereε is the axial strain rate andκ is the rate of curvature. The axial strain and curvature rates, at given N and M, can be determined using the variation of the rate analogue of the total potential energy where the dual potential Ψ is The variation of Π with respect to the ratesε andκ, yields N N 0 = n 2 (n + 1)κ 0 κ is the traction at the start of the increment and Θ is an integration parameter. It should be noted that T e and H i are calculated at the current increment.
In order to solve these equations, we encapsulate the equations in a vector F = [F t F n ] T and define the unknown vector ΔT = [ΔT t ΔT n ] T . Using a Newton-Raphson algorithm, the following linearised form of these equations is iteratively solved for the traction increment δΔT until a convergence condition is achieved where 0 is a zeros vector (2 × 1) and the Jacobian, J, is defined as The interface tangent modulus is defined by (B.7) The elements of the tangent modulus can be determined from the variations of the incremental forms in Eqs. (B.3). The variation of F i is given by where the variations of the separation increments δΔδ t and δΔδ n can be chosen arbitrarily yielding ∂ F i /∂Δδ t = 0 and ∂ F i /∂Δδ n = 0. Hence, the variations of Eqs. (B.3) give the following set of equations (B.10) The derivative of the damage parameters are taken as ∂d i /∂Δδ j ≈ ∂d i /∂Δδ cr j . Therefore, solving Eq. B.9 gives the components of the Jacobian. and T t = σ 0 C * ε 0 σ 0 I n r 1 n+1σ r θ n, 0, χ pl .
(C.3) Therefore, the effective traction for the different models is evaluated from Eqs. (14), (19) and (25)  (C.13) The other limit is when the interface is too compliant in comparison with the bulk material which can be regarded as rigid. We have determined that the physical domain for crack growth data is far from this limit (Elmukashfi and Cocks 2017). Hence, we limit considerations to the pure mode I and mode II cases. In the case of pure mode I (Elmukashfi and Cocks 2017), C k for the case of a compliant interface is determined as where q n = 2/ (2 n − 1), W = L − a is the length of the remaining ligament during steady state propagation andW = W/λ. In the case of pure mode II, using the assumption that the bulk material is rigid, crack growth takes place when the separation at the tip is equal to the tangential deformation of the DCB arms, i.e. δ f t = δ. Thus, at that instant, the interface reaches the critical value leading to unstable crack growth. Therefore, for an infinite crack growth rate, C k is unconditionally bounded in the compliant limit, i.e. C k ≤ ∞.
Another limitation comes from the validity of C * as a parameter for characterising the near tip field and damage growth process. In particular, at higher crack velocities, in the vicinity of the crack tip, the elastic deformation becomes increasingly important. In this case, a zone in which the response is determined by both elastic and creep deformation will be formed. Hence, when the size of this zone is comparable in size to the damage zone, C * can no longer be used. Following Cocks and De Voy (1991), C * controls crack growth provided the following condition is satisfied where f n (n) = 4n (2n + 1) (n + 1) , E is Young's modulus and Z (n) = (n − 1) I n n−1 n+1 . Using this condition we derive a condition for C s function by comparing Eq. (C.15) with Eq. (44) as (C.16) Hence, for particular values ofδ f e and σ 0 /E, there is a maximum velocity below which C * is a valid measure.