A theoretical and computational framework for studying creep crack growth

In this study, crack growth under steady state creep conditions is analysed. A theoretical framework is introduced in which the constitutive behaviour of the bulk material is described by power-law creep. A new class of damage zone models is proposed to model the fracture process ahead of a crack tip, such that the constitutive relation is described by a traction-separation rate law. In particular, simple critical displacement, empirical Kachanov type damage and micromechanical based interface models are used. Using the path independency property of the \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}C∗-integral and dimensional analysis, analytical models are developed for pure mode-I steady-state crack growth in a double cantilever beam specimen (DCB) subjected to constant pure bending moment. A computational framework is then implemented using the Finite Element method. The analytical models are calibrated against detailed Finite Element models. The theoretical framework gives the fundamental form of the model and only a single quantity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{C}_k$$\end{document}C^k needs to be determined from the Finite Element analysis in terms of a dimensionless quantity \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}ϕ0, which is the ratio of geometric and material length scales. Further, the validity of the framework is examined by investigating the crack growth response in the limits of small and large \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}ϕ0, for which analytical expression can be obtained. We also demonstrate how parameters within the models can be obtained from creep deformation, creep rupture and crack growth experiments.

growth response in the limits of small and large φ 0 , for which analytical expression can be obtained. We also demonstrate how parameters within the models can be obtained from creep deformation, creep rupture and crack growth experiments.
Keywords Creep · Crack · C*-integral · Damage zone model · Traction-separation rate law (TSRL) · Double cantilever beam (DCB) · Dimensionless analysis Nomenclature 2 l The spacing between two adjacent pores β A material parameter of the exponential damage law δ c n The critical normal displacement jump in the damage zone at the crack tip δ f n The normal displacement jump at failure in the crack tip δ i The displacement jump vector across the damage zone (i = 1, 2, 3) δ 0 The separation rate at the reference traction T 0 δ m n The maximum normal displacement jump rate vector in the crack tiṗ δ i The displacement jump rate vector across the damage zone (i = 1, 2, 3) ε 0 The strain-rate at the reference stress σ 0 a The steady state crack velocitŷ C k The separation history function of model k λ The characteristic geometric length scale (•) cr The creep component of the quantity (•) (•) el The elastic component of the quantity The dimensionless steady state crack velocity φ 0 The ratio of geometric to material length scales σ 0 The reference stress σ e The von Mises equivalent stress σ i j The Cauchy stress tensor ε i j The engineering strain tensor a The crack length C * The rate of the J -integral C s The separation history function of the simple model E Young's modulus 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 h The current height of a pore h 0 The initial height of the pores m The rate sensitivity exponent of the damage zone n The rate sensitivity parameter of the bulk material n i The unit normal vector (i = 1, 2, 3) s i j The deviatoric part of Cauchy stress tensor T i The traction vector (i = 1, 2, 3) T 0 The reference traction of the damage zone u i The displacement vector (i = 1, 2, 3) The Cartesian material and spatial coordinates (i = 1, 2, 3)

Introduction
At elevated temperature, creep crack growth (CCG) is one of the most common failure mechanisms in many engineering applications, e.g. structural components, similar and dissimilar metal welds etc. This problem has received much attention over the last forty years due to the importance in designing structures with high integrity and safety. Hence, developing analytical models for steady-state crack growth which can be calibrated against detailed Finite Element models is of great interest. Further, an assessment of the effect of different material parameters and damage development processes on the crack growth behaviour can be provided using such models. Studying creep crack growth has a long history in the literature. A major feature of these studies is the development of a parameter that characterizes the crack tip fields as well as crack propagation. Under steady state creep conditions, 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) can be used to characterize the crack tip fields and creep crack growth. 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. Riedel and Rice (1980) studied the transition from short-time elastic to long-time creep behaviour assuming that primary creep is negligible (small-scale creep conditions). They introduced a parameter C(t) that describes the strain, strain-rate and stress fields within a creep zone that forms about the crack tip. Their analysis also provides a characteristic time for the transition to the steady state stress field (i.e. the time for C(t) to equal C * ). Later, Ehlers and Riedel (1981) proposed a relation between C(t) and C * . Saxena (1986) proposed a new parameter C t which can be measured easily in comparison with C(t). Bassani et al. (1988) compared these two parameters and concluded that C t characterizes crack growth rates much better than C(t). Further, the C(t) parameter is found to be more suitable for characterizing a stationary crack and C t is related to a rapidly propagating crack. In the primary creep regime, Riedel (1981) suggested a new parameter C * h as an analogy to the C * -integral. Further, Leung and McDowell (1990) included the primary creep effects in the estimation of the C t parameter. To this end, the C * , C(t), C t and C * h parameters are generally accepted and widely used in studying creep crack growth.
Under creep conditions, cracks in polycrystalline materials advance as a result of the growth of damage ahead of the crack tip (generally in the form of discrete voids or microcracks, which form primarily at grain boundaries). In the vicinity of a macroscopic primary crack tip, secondary micro-cracks are formed as a result of intensive void growth and coalescence and/or an accumulation and growth of micro-cracks. These secondary cracks propagate and coalesce creating the new crack surfaces, allowing the primary macroscopic crack to advance along an interface or interconnected grain-boundaries. The growth of damage can influence the constitutive properties of the material and therefore the details of the near crack tip stress and strainrate fields. Early models of creep crack growth either assumed that the stress (e.g. Riedel 1981;Tvergaard 1984) or strain-rate field (e.g. Cocks and Ashby 1981;Nikbin et al. 1984) is the same as that for the undamaged material and used either empirical or mechanistic damage growth laws to determine the crack growth rate. In the strain based models the critical damage at the crack tip is expressed in terms of a material ductility (strain to failure) which is a function of the local stress state. Extensions of this approach within a finite element framework (employing models in which the constitutive relationships for deformation are not influenced by the presence of damage) have been undertaken by Nikbin et al. (1976Nikbin et al. ( , 1984, Yatomi and Nikbin (2014). Studies of the influence of damage on the nature of the crack tip fields and crack growth process where damage influences the deformation response have been undertaken by Riedel (1987) and Bassani and Hawk (1990) for empirical Kachanov (Kachanov 1958;Rabotnov 1969) type continuum damage mechanics models. More recently, the full interaction between deformation and damage development and how this influences the crack growth process has been modelled directly using the finite element method, using both mechanistic and empirical models for the growth of damage (e.g. Onck and van der Giessen 1998;Wen and Shan-Tung 2014).
In each of the above referenced studies damage development and its influence on crack growth is modelled as a continuum process. 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, i.e. they introduce a physically meaningful length scale that is related to the dissipative mechanisms responsible for damage development. A damage zone model of this type describes the fracture process in the vicinity of the crack tip as a gradual surface separation process, such that the normal and shear tractions 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). The first use of cohesive zone models in a finite element environment was undertaken 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;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). Rate-dependent and rateindependent models as well as physically based and phenomenological models have been employed. However, to the authors' knowledge, apart from the work of Onck and van der Giessen (1998), van der Giessen and Tvergaard (1994), Thouless et al. (1983) and Yu et al. (2012) damage zone type models have not been used to study the development of creep damage and/or creep crack growth.
In this paper, a theoretical and computational framework for creep crack growth is presented in which we assume that all the damage is concentrated in a narrow zone directly ahead of the growing crack tip. The objective is to model crack propagation in materials that exhibit steady state creep behaviour outside of the damage zone and to investigate the effect of different material parameters, forms of damage zone constitutive law and damage development processes on the crack growth behaviour. A theoretical framework is initially introduced in which the constitutive behaviour of the bulk material is described by powerlaw creep. A new class of damage zone model is proposed to model the fracture process such that the constitutive relation is described by a traction-separation rate law. More specifically, three different models, i.e. a simple critical displacement model, Kachanov type empirical models and a micromechanical based interface model are investigated, which mirror the types of models employed in the continuum models of creep crack growth described above. We follow the recent approach of Wang et al. (2016), who studied cleavage failure in creeping polymers, in which we keep the material descriptions and geometric configurations as simple as possible to explore the relationship between the form of the constitutive model, material parameters and crack growth. With this in mind, we concentrate initially on the behaviour of a double cantilever beam specimen (DCB) of infinite length subjected to a constant pure bending moment, in which C * remains constant as the crack grows and the crack growth rate eventually achieves a steady-state. We concentrate on the behaviour in the steady state. By invoking the path independence of the C * -integral and choosing contours in the far field and surrounding the damage zone we demonstrate how a simple analytical expression for the crack growth rate can be obtained in terms of C * , damage zone material parameters and a dimensionless scaling parameter that is a function of the ratio of characteristic geometric and material length scales, that can be determined using the finite element method. The theoretical framework is presented in Sect. 2. 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. 3, with the crack growth results for the different interface models presented and discussed in Sect. 4.

Background
Consider a solid containing a stationary crack that is subjected to a constant load. The solid is assumed to exhibit elastic behaviour, together with primary, secondary and tertiary creep. Creep deformation evolves with increasing time and this evolution can be divided into different distinct stages. These stages have been described and evaluated by Bassani and Hawk (1990). Initially, a small-scale creep zone, i.e. small in comparison with the physical characteristic length of the body, is formed in the vicinity of the crack tip. In this stage, the material deforms by primary creep inside the creep zone and remains elastic elsewhere. Following the development of the primary creep zone, a secondary (steady-state) creep zone develops as a smaller region inside the primary creep zone. Thereafter, the primary and secondary creep zones continue to expand at the cost of the elastic and primary zones, respec-tively. During this process, damage accumulates in the crack tip region, which may lead to crack propagation if a critical condition is met. Hence, crack propagation may take place at different instants during the evolution of the near tip stress and strain-rate fields. The crack propagation scenarios are characterized by the nature of the crack tip fields at these instants: (i) the small-scale creep zone is formed surrounded by the elastic medium, (ii) the primary creep zone is large enough but remains surrounded by the elastic medium, (iii) the secondary creep zone is formed inside the primary creep zone but both zones remain surrounded by the elastic medium, (iv) the secondary creep zone is expanding inside the primary creep zone which dominates, (v) the secondary creep zone dominates.
This study concerns crack propagation in creeping materials under steady state conditions (type v). In this case, the C * -integral can be used to characterize the creep crack growth behaviour. Note also, that as a crack grows in an elastic/creeping material, in the absence of damage, a zone develops ahead of the crack tip in which the stresses are determined by the elastic and creep properties of the material (Hui and Riedel 1981) who's size is a function of the crack velocity. For steady state behaviour the size of this zone must be small compared to the size of the crack tip damage process zone. Under these conditions, the path independent property can be used to obtain a direct relationship between the far field loading and the fracture process parameters.
In the following sub-sections we describe the constitutive relationships for the bulk continuum response and introduce a number of different models to describe the response ahead of the crack tip within the damage zone. We concentrate on mode I crack growth and only present relationships for the opening mode, although a description of the shear response is also required for the computational studies presented later. We consider the crack growth process in Sect. 3, where we use the path independence of the C * -integral to relate the near crack tip damaging processes to the far field loading. In the theoretical models presented below we concentrate on steady state crack growth, where creep dominates the material response, i.e. we need not consider the elastic response. Similarly we do not need to consider any elastic/reversible contributions to the deformation within the damage zone.

Creep deformation and characterization of the remote field
Consider a body containing a crack and 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. The bulk material is assumed to exhibit steady-state creep behaviour and is defined by the constitutive laẇ where σ i j is 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 = 3 2 s i j s i j 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,Ḋ, is given bẏ where ψ is the dual rate potential andε e = 2 3ε i jεi j . Crack propagation is assumed to be determined by an interface model such that the propagation takes place along a fictitious interface surface, Γ int . As the crack advances separation occurs along the interface to create two surfaces. Hence, 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 surface and the corresponding tractions are defined by the vectors are the displacement rates either side of the interface, and T + i = σ i j n + j , T − i = σ i j n − j . In order to analyse the problem the C * -integral is determined on the inner and the 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. We consider the situation where a damage zone extends along the x 1 axis directly ahead of the crack tip, see Fig. 2. The first term of Eq. (5) then vanishes, since dx 2 = 0. The contribution of the damage region to the C * -integral along the inner path, C * in , is then evaluated as Fracture zone x 1 x 2 Fig. 1 The schematic of interface crack model in creeping solid. The figure illustrates the definition of the interface surface Γ int . inner path Γ in and the outer path Γ out The cohesive zone for pure mode I crack propagation: a schematic of the cohesive zone; and b the normal tractionseparation (T n − δ n ) and the separation rate-separation (δ n − δ n ) distribution along the cohesive zone. l int is the length of the cohesive zone, δ c n is the critical displacement, δ f n is the displacement at failure, and σ c n is the cohesive strength whereδ m n is the opening rate at the tip of the crack. In order to simplify the relationships used in subsequent analysis we omit the superscript "+" from the traction. Using the path-independence property of C * (ie, C * out = C * in ) provides a relationship between the far field loading and the behaviour within the damage zone. In order to complete the analysis we need a constitutive relationship for the damage zone that relateṡ δ i to T i . In this paper we concentrate on mode I loading and therefore only need to consider the components of traction and displacement-rate normal to the damage zone. Then In the following sub-section we present a number of different constitutive relationships for the damage zone response.

Damage zone models for creep crack growth
In this section, we present models for the damage zone. We limit our description to mode I loading. More general relationship for mode II and mixed mode loading are described elsewhere . For each of the models here we assume that the relationship betweenδ n and T n can be expressed in the form of a power-law. Further, these models are capable of predicting similar damaging processes and crack growth behaviour through the appropriate selection of material parameters, see "Appendix B".

Simple critical displacement model
The simplest form of traction-separation rate law is defined by a direct power-law relationship between the normal traction T n and the separation rateδ n . Further, damage does not influence the opening rate and failure is achieved when δ n achieves a critical separation δ f n . Thus, the constitutive model takes the following forṁ where T 0 is a reference traction (equivalent to 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). It should be noted that the traction becomes zero when the separation exceeds the critical value, i.e. T n = 0 if δ n ≥ δ f n .

Kachanov damage type empirical model
In this model damage is assumed to influence the constitutive response and a single scalar damage parameter ω is introduced to incorporate the effect of damage. The damage parameter is assumed to evolve monotonically from 0 to 1, i.e. from an undamaged to a fully damaged state. Following Kachanov (1958) and Lemaitre and Chaboche (1994) we assume that the separation rate is a function of an effective tractionT n , which is related to the normal traction T n and ω bȳ The constitutive response is then simply obtained by replacing T n byT n in Eq. (8) to givė The damage is assumed to be determined by the normal separation δ n , i.e. ω = ω (δ n ), and the damage evolution law may take different forms depending on the damage mechanism(s). In this study, we propose two different damage models, namely linear and exponential models. The different evolution laws are: -The linear damage model: -The exponential damage model: where δ c n is the separation at which damage initiates (for separations less than this value ω = 0 and the constitutive response is given by Eq. (1)), as before δ f n is the 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. Figure 3 below shows the traction-separation law for different separation rates for the case of linear and exponential damage laws.

Micromechanical based model
This model is based on the creep extension of Yalcinkaya and Cocks (2015) micromechanical damage zone model for ductile fracture described by Cocks et al. (2017), which are both derived from the creep cavitation model of Cocks and Ashby (1980). These models are based on the growth of an array of pores idealized as cylinders. The relation between the macroscopic traction and separation and the microscopic stress and strain is then obtained using classical bounding theorems. The radius and height of the pores at a given instant are denoted as r and h, respectively, and the mean spacing is 2l, see Fig. 4. Thus, the pores are characterized by their area fraction in the plane of the cavitated zone, i.e. by f = (r/l) 2 . Further, the representative volume element is assumed to be fully constrained in the radial direction and the deformation is only controlled by the normal separation, i.e. l = const. andḣ =δ n . This simplification of the void profile captures the major features of the evolving geometry (such as area fraction of pores and pore aspect ratio) while allowing simple analytical expressions for the evolution of damage to be derived. More general forms of model are discussed by Cocks et al. (2017)-this is the simplest form of model of this class and it is directly equivalent in form to classical rate dependent cohesive zone models, including those described above. The resulting expression for the 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 idealization of a pore as a cylinder with height h and diameter equal to the pore to pore spacing 2l and the macroscopic normal traction T n and separation δ n opening rate iṡ whereḡ n = g n ( f )/g 0 and and g 0 is a parameter that can be used to provide a similar rupture time to the Kachanov models under the same stress level and the same value of the material parametersδ 0 , T 0 and δ f n (see "Appendix B"). The matrix material is incompressible, therefore the total rate of change in volume is equal to the rate of change in pore volume. Hence, the pore area fraction evolves at a ratė The initial pore area fraction and height are assumed to be f 0 and h 0 , respectively. Further, the pores are assumed to coalesce and reach a complete failure when f = f c . The direct integration of Eq. (15) gives the separation at failure: The traction-separation rate relation for different separation rates is illustrated in Fig. 5.

Creep crack growth in a double cantilever beam
In this section, pure mode-I 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 2H , respectively, and the crack length is denoted by a. Each arm of the specimen is subjected to a constant moment M per unit depth. 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  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 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.
In order to determine the relationship between the loading and fracture parameters, the path independence of the C * -integral is used 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. 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 damage zone 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. We do not know the form of this function a priori, but to determine the crack growth rate we only need to determine a single quantity-the resulting integral, which can be determined from a single piece of information from a finite element analysis of the problem. The details of this process are given below.
3.1 The C * -integral in the outer path Γ out C * in the outer path can be determined in a number of different ways-for example by direct evaluation of the integral of Eq. (5) or by determining the rate of change of the rate analogue of the total potential energy, , with crack length, where We adopt the second of these approaches here. For an element of beam under bending the curvature rate is given bẏ where M 0 = 2n 2n + 1 σ 0 H 2 4 , and η = 1 for plane stress and √ 3/2 for plane strain. For the DCB specimen of The factor of 2 arises because there are 2 beams andθ is the rotation rate at the end of one of the beams. If we define a reference stress, such that and where f n (n) = 4n (2n + 1) (n + 1) and λ = H/2 is the characteristic length scale for the DCB specimen.
3.2 The C * -integral in the inner path Γ in Using the definition in Eq. (7), the C * -integral in the inner path Γ in can be written as where the function d k depends on the form of interface model adopted, with k indicating the model, i.e. s ≡ simple, kl ≡ Kachanov linear, ke ≡ Kachanov exponential and m ≡ micromechanical models. For each of these models d k is given by Apart for the simple model the integral requires a knowledge of the stress history experienced by each material point in the damage zone, which is not known a priori. For the simple model the integral of Eq. (22) can be readily determined: For all the remaining models 0 ≤ d k ≤ 1, and the resulting integral is therefore less than or equal to that given by Eq. (24). We assume that the integral for each model can be approximated by where the subscript k again identifies the model and α k falls in the range 0 ≤ α k ≤ 1.

The crack tip opening displacement rate
Equating the values of C * given by the inner and outer contours, i.e. Eqs. (21) and (25), allows the crack tip opening displacement rateδ m n to be expressed as a function of the applied loading.

The analysis of a steadily propagating crack
As noted earlier, under a constant applied moment the crack velocity will, after an initial transient, achieve a steady state, in which it reaches a constant valueȧ. 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 1981) The are two characteristic length scales in this problem, the geometric length scale λ and the material length scaleδ 0 /ε 0 . We can therefore write the separation rate in the forṁ where Λ k is a dimensionless function that depends on the interface model whose detailed form depends on φ 0 , the ratio of geometric and material length scales, as shown in Fig. 8b. For each of the models described in Sect. 2, failure of an element occurs when the separation across the damage zone reaches a critical value, δ f n . Integrating the displacement rate as an element is convected towards the crack tip gives 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ȧ where we have used Eq. (27) to substitute for dt and Eq. (28) to substitute forδ n . The dimensionless function C k (φ 0 , n, m) is only a function of φ 0 , n, m and the detailed form of model for the damage zone. Substituting forδ m n using Eq. (24) gives the steady state crack velocity.
whereδ f n = δ f n /λ. 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 n . This is only true ifĈ k is only a function of n and m or φ 0 is constant for the range of conditions of interest. Note that for the DCB specimen of Fig. 7 and m = n where B =ε 0 /σ m 0 is a material property. Then for a series of experiments in which the geometry is kept constant we would expectȧ to be proportional to C * m m+1 , but the constant of proportionality could be different for a different choice of beam height H . For more gen-eral cracked geometries σ 0 ,ε 0 and λ change as a crack grows and therefore φ 0 also changes. This needs to be taken into account in any model and description of the crack growth process. We consider this feature of the response further below.
In order to determine the crack growth rate we need to evaluate the quantityĈ k . We can determine this using the finite element method. We need not determine the distributions d k and Λ k ahead of the crack tip. We can determine directly by equating the numerically determined crack velocity with the prediction of Eq. (30).

Numerical implementation of the governing equations
The initial-boundary value problem described in Sect. 3 is numerically solved using the FE (Finite Element) code ABAQUS (Abaqus 2016). A nonlinear quasistatic 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 an elastic-rate dependent opening model for the damage zone. Elastic constitutive components have been added to the constitutive relationships of Eqs. (1), (8), (10) and (13), with the values of the elastic components chosen to have limited influence on the computed results. The geometry of the double cantilever beam (DCB) specimen shown in Fig. 6 is discretised, and a typical finite element mesh is shown in Fig. 9. Only one half of the specimen is analysed due to the symmetry of the problem. The dimensions are taken as L = 100 mm, H = 10 mm, and B = 1 mm. The initial crack is assumed to be a = 40 mm, and the crack propagation The damage zone elements are modelled with zero initial thickness such that the top and bottom face nodes coincide. The mesh has 11,634 elements, of which 11,279 are bulk elements and 355 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 convergence study on the mesh refinement was carried out for different values of φ 0 and interface parametersδ f n = 0.004, β = 1.0, f 0 = 0.01 and f c = 0.91. We found that an interface element length of l inte = 0.1 mm is necessary to obtain converged solutions for the range φ 0 ∈ 10 −5 − 10 4 . The interface stiffness K n = K t = 10 6 MPa · mm is selected such that the elastic deformation is negligible (see "Appendix A").
The numerical analysis was performed for different combinations of the dimensionless parameters defined above to confirm that the functional form of Eq. (30) is valid. (The model parameters are chosen in such a way that the dimensionless parameters are controlled.) The relative normal separation displacement, Δu 2 , between each pair of initially coincident nodes in the interface (x 2 = 0) is computed and recorded during the analysis. The crack tip position, x tip , is defined by Δu 2 = δ f n , and the crack tip velocity is determined using forward differencing as 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 crack growth
Several analyses have been performed for different combinations of the dimensionless parameters and damage zone properties. The crack tip position, transient and steady state crack propagation velocity have been obtained for all the combinations. For the simple interface model the parameters n = m = 9, φ 0 = 1.0 andδ f n = 0.004 are used to illustrate the different results. We first describe the results under plane stress conditions. Figure 10i-iv show the distribution of the effective creep strainε cr e at four different instants and crack velocities. At t = 0 the creep deformation is zero everywhere and elastic deformation prevails. As time passes x 2 [mm] x 1 [mm] (t > 0) creep deformation evolves in the bulk, initially primarily in the vicinity of the crack tip, as well as damage along the interface, leading eventually to crack growth when the critical opening is achieved at the crack tip. Figure 11a, 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 and after a short time (10 h) a lower steady state velocity is achieved. In this case a steady velocity ofȧ = 0.311 mm/h 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 and as the crack advances the creep deformation dominates where the stress relaxes leading to a slower propa-gation rate. Additionally the damage 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 propagation before reaching a steady state where a fully developed damage zone is achieved.

The C s -function
It proves instructive to concentrate initially on the response for the simple damage zone model of Sect. 2.3.1. The Finite Element analysis is used to determine the steady crack velocityȧ and then for a given set of input parameters the C s -function can be evaluated from Eq. (30): The appropriateness of the dimensionless analysis has been examined using the same set of dimensionless parameters with different model parameters, e.g. the same value of φ 0 with different combinations ofε 0 , H andδ 0 .

The physical limits and the validity of the framework
Before evaluating the computational results in detail it is instructive to examine the response in the limits of small and large φ 0 . The first extreme is when the interface is very stiff in comparison with the bulk material (the bulk material creeps faster than the interface, i.e. ε 0 δ 0 /λ and φ 0 → ∞). The other extreme occurs when the interface creeps faster than the bulk material, i.e. the interface is very compliant (ε 0 δ 0 /λ and φ 0 → 0). In this analysis we consider the simple damage zone model of Sect. 2.3.1. When an interface is very stiff in comparison with the bulk material the deformation along the interface is negligible and it does not influence the stress state in the body. The tractions seen by the damage zone are determined by the stress distribution in the bulk material, and can be expressed in terms of the C * -integral (provided the damage zone is small compared to the region in which the HRR field dominates; Hutchinson 1968; Rice and Rosengren 1968). The HRR stress field is defined as where I n is an integration constant that depends on n andσ i j is a dimensionless function of n and θ . The values of these parameters are given for the cases of plane stress and plane strain conditions by Hutchinson (1968). It follows that the normal traction along the interface is given by In this analysis, we limit ourself to the case of T 0 = σ 0 and m = n. Therefore the opening separation rate for the simple model is evaluated from Eq. (8) aṡ The critical opening separation is determined by integrating the separation rate, in a similar way to in Eq. (8), as where r = x 1 at θ = 0 and r c is the size of the fracture zone which is very small in the case of stiff interface (r c → 0). Rearrangement of Eq. (38) gives the dimensionless velocity as wherer c = r c /λ. By comparing this equation with Eq. (30), the C s -function for the case stiff interface becomes For the case of n = m = 9 and θ = 0, the integration constants are I 9 ≈ 3.025 andσ θ (9, 0) ≈σ θ (13, 0) ≈ 1.2 for the case of plane stress and I 9 ≈ 4.6 and σ θ (9, 0) ≈σ θ (13, 0) ≈ 2.6 for the case of plane strain (Hutchinson 1968). Thus, the C s -functions for the cases of plane stress and plane strain conditions are C s = 45.6 ·r 0.1 c · φ −0.9 0 and C s = 7.1 · 10 4 ·r 0.1 c · φ −0.9 0 , respectively.
The other limit is when the interface is too compliant in comparison with the bulk material which can be regarded as rigid. Hence, in the case of an infinite DCB specimen, the equilibrium between the applied moment and the traction along the infinite damage zone suggests that the traction will tend to zero and there will be no crack propagation. on the other hand, when the specimen is finite a non zero traction along the finite damage zone. The deformation along the damage zone can directly be related to the angular deflection at the end of the beam. Thus, the separation at the crack tip is obtained as where W = L − a is the length of remaining ligament during steady state propagation. The opening displacement at the tip of the propagating crack is constant and equal to the critical value δ f n = 0, thereforeδ n = 0, anḋ Similarly, the separation at the crack tip can be written in this forṁ The separation rate in the damage zone is given bẏ Now the balance by the internal and external work rates gives Introducing Eq. (43) W is computed from the finite element analysis as the remaining ligament length when a crack reaches a steady state propagation. Hence, for the case of n = m = 9 and using the computationally obtained average valueW ≈ 0.8, the C s -functions for plane stress and plane strain conditions are C s = 10.5 × 10 −15 · φ −0.9 0 and C s = 38.3 × 10 −15 · φ −0.9 0 , respectively. Another limitation comes from the time scale of the crack propagation as mentioned in Sect. 2.1. C * represents the near crack tip field when a crack propagates slowly. As the crack velocity increases elastic deformation becomes increasingly important in the vicinity of the crack tip and a zone in which both elastic and creep deformation determines the response becomes increasingly significant. If this zone becomes comparable in size to the damage zone, then C * can no longer be used as a parameter for characterization of the near tip filed and damage growth process. Cocks and Julian (1991) studied this limit and proposed conditions for the dominance of C * . They demonstrate that C * controls crack growth provided the following condition is satisfied where 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. (50) with Eq. (30) as This expression implies that for particular values ofδ f n and σ 0 /E there is a maximum velocity for which C * is a valid measure. Thus, for the case of n = m = 9, the valid C s -function for plane stress and plane strain conditions are C s ≤ 4.28 ·r respectively. Elasticity is only relevant in the computational models and this relationship can be used to assess whether the conditions employed in the FE models are consistent with the assumptions of the analytical model presented in Sect. 2. We need to be careful, however, when using this expression. It is derived from analyses in which damage development is assumed to not influ-ence the near tip fields. As illustrated above, the size of the damage zone increases with decreasing φ 0 and for small φ 0 the near tip fields given by the classical continuum analysis are no longer valid. The relationship of Eq. (49) is therefore only valid in the limit of large φ 0 where the development of damage has limited effect on the crack tip fields. It is also important to emphasise here that although, the HRR field is no longer valid for small φ 0 , C * is still a valid parameter for characterizing crack growth.
In order to evaluate the proposed framework, C s has been determined from (34) for φ 0 in the range [10 −10 , ×10 5 ] and compared with the limiting results presented above. The rate sensitivity parameters are taken to be n = m = 9. Figure 12a, b show the relationship between C s and φ 0 for plane stress and strain conditions, respectively.
Over the range of the data, the results can be fit using two separate power-law relations. Under plane stress conditions this relation is C s = 0.45 φ −0.06 0 over the range of values φ 0 ∈ [10 −10 , 8 × 10 −2 ] and C s = 0.09 φ −0.67 0 for the range φ 0 ∈ [8 × 10 −2 , 10 3 ], see the dashed lines in Fig. 12a. The transition between the power law relations occurs over the range 10 −4 ≤ φ 0 ≤ 10 0 . For a given value of φ 0 , C s lies between the two limiting values. The power-law fit for high values of φ 0 is slightly shallower than that for the stiff limit described above, indicating that response tends to this limit for values of φ 0 in excess of 10 6 . In this limit the rate of deformation in the damage zone becomes very small compared to that in the surrounding matrix, which determines the stress distribution ahead of the crack tip and therefore the rate of growth of damage. There is no evidence of the data merging to the limiting result for low values of φ 0 , but the values of φ 0 required to reach this limit are much lower than values we would expect from physical arguments (in this limit the material length scale is significantly greater that the geometric length scalein practice we would expect any characteristic material length scale to be less than the geometric length scale for the cracked body, i.e. we would expect φ 0 to be greater than 1). The power-law range for φ 0 greater than 8×10 −2 is therefore more representative of the physical behaviour of engineering components, so we concentrate on the relation for this regime here. Substituting this relationship into Eq. (30) gives the dimensionless velocity  Fig. 12 The relation between C s -function and φ 0 parameter and the physical limits in the case of n = m = 9: 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 n , and the dashed lines show the power law fit a = 1.18 × 10 −2 · φ −0.77 where we have substituted for C * using Eq. (21) to provide a relationship in terms of the reference stress σ 0 . Figure 12a also shows a series of lines below which elastic effects can be ignored for σ 0 /E = 8 × 10 −6 (as used in the computations) and different values of critical crack tip opening displacement, i.e. below which inequality (51) is satisfied. As noted earlier this relationship is only valid for large values of φ 0 (say greater than 8×10 −2 ). In this regime the computational results lie below this series of lines, indicating that the theoretical structure presented in Sect. 4 provides a valid framework for modelling the crack growth behavior. We can repeat the analysis for plane strain conditions, see Fig. 12b. In the case of plane strain we again find that the results can be fit using two power-law relationships: C s = 0.55 φ −0.05 0 over the range of values φ 0 ∈ [10 −10 , 10 −2 ] and C s = 0.19 φ −0.29 0 for the range φ 0 ∈ [10 −2 , 10 4 ], see the dashed lines in Fig. 12b. Further, the transition between the power law relations takes place in the range 10 −4 ≤ φ 0 ≤ 10 0 . The latter relation gives a dimensionless crack growth rate of The comparison between the FE results and the physical limits is also shown in Fig. 12b, which are again bounded by the physical limits of Eqs. (40) and (49), with the results asymptoting to Eq. (40) at large values of φ 0 . The large limit φ 0 gives a faster crack growth rate in plane strain than plane stress (i.e. C s is larger for a given value of φ 0 ) due to the higher stress levels ahead of a plane strain crack. The difference in slope between this limit and the computational results is greater than that observed for plane stress, but the value of φ 0 where the two curves meet is about two orders of magnitude higher. As for plane stress, the results lie in a regime where elastic effects can be ignored.

The effect of damage model
In order to investigate the effect of the detailed form of the damage zone model on the crack growth response, C k has been determined for each of the different damage zone models described in Sect. 2. The parameters employed for these models are m = 9, δ f n = 0.02 mm, β = 1.0, h 0 = 0.02 mm, f 0 = 0.01 and f c = 0.5. It should be noted that δ n is kept constant for all models.
Further, we choosed g 0 = 2.23 and β = 1.0 such that all models yield the same rupture time t f under a prescribed constant stress. (In "Appendix B" we demonstrate that under a given stress and for prescribed values ofδ 0 and δ f n the time to failure is proportional to a dimensionless quantityĨ k , see Eq. (B.6). We choose the values β and g 0 in the exponential Kachanov and micromechanical models such thatĨ k , and therefore the time to failure, is the same for all the models). In the studies of crack growth, the physical length scale of the cracked body λ = H/2 is used and the matrix rate sensitivity parameter is taken to be n = 9, as before. Here, we limit our consideration to plane stress conditions, but similar results can be obtained under plane strain. C k is determined using the same procedure as described above, by comparing the computed steady state crack growth rate with Eq. (34). As before, we can determine analytical relationships for the response in the limits of small and large φ 0 . In Sect. 4.3 we found that the analytical model in the limit as φ 0 → 0 does not provide a meaningful bound to the results and we do therefore do not present results in this limit for the remaining damage zone models described in Sect. 2. The analytical results for these models in the limit φ 0 → ∞ are presented in "Appendix C". Rather than expressĈ k as a function of φ 0 it proves instructive to express it as a function ofĈ k (which is a function of φ 0 ) for each of the models. Figure 13 shows the relationship betweenĈ k andĈ s . The results show that these relations are nonlinear. However, for each model there is point-wise a linear relationship between these two functions witĥ where μ k is a parameter that depends on the damage model used. The effect of damage is to soften the constitutive response, effectively increasing the effective separation rate across the damage zone for a given traction T n , which results in a lower effective value of φ 0 and therefore an increase of the crack growth rate. Therefore μ k is larger than 1.0 for a given value of δ f n . Values μ k are given in "Appendix C" for the three models: μ kl = μ ke = μ m = 10. Substituting these values of μ k into Eq. (55) give the straight lines plotted in Fig. 13, which also shows the computational results.
The computational results approach the analytical results for small values ofĈ s , i.e. large values of φ 0 , see Fig. 13, which corresponds to the limit where we would  Fig. 13 The relation betweenĈ k and C s (which is a function of φ 0 parameter-see Fig. 12a for the different models under plane stress conditions. The solid line represent the response given by the model presented in "Appendix C". The dashed line represent the apparent asymptotic behaviour. The chained line corresponding to μ k = 1 represents a strict lower bound to the data expect the analytical result to apply. As is increased gradually reduces for all the models and then asymptotes to a value μ kl = μ ke = μ m = 1.8. AsĈ s increases deformation in the damage zone becomes constrained and the stress relaxes more for a damaged material for a given value of φ 0 compared to a material that does not damage. Therefore, the elevation of the crack growth rate is less. In the limit that deformation is completely constrained [corresponding to the situation considered by Cocks and Ashby (1981)] the crack velocity is independent of the details of the model and only depends on the critical opening displacement δ f n ; thus in this limit we would expect μ k to equal 1 for all the models considered here. This represents a physical lower bound to μ k , which is not reached for any of the models with the results appearing to asymptote to a higher limiting value over the range of conditions considered in the computations.

Comparison with experimental data
The main objective of this paper is to identify simple constitutive models for the damaging process ahead of a crack tip in a creeping material and to identify a simple structural configuration which can be analysed rigorously to provide new insights into the relationship between damage models and the crack growth process, including the role of different characteristic material and geometric length scales. As a result, the simple geometry and loading conditions considered are not representative of laboratory test components. Nonethe-less, it proves instructive to explore how the models presented here can be calibrated against available experimental data, to determine the characteristic material and geometric length scales in these experiments and explore where this data lies with respect to the general trends identified in Figs. 14 and 15.
In this section, we consider the low alloy steel (2.25 Cr Mo steel at 538 • C) investigated by Nikbin et al. (1983). They provide data for creep deformation, creep rupture as well as creep crack growth generated using compact tension (CT) specimens, see Figs. 14 and 15. Consider the creep crack law of Eq. (31), together with the definition of φ 0 in Eq. (32) or following Eq. (26). In order to determine the crack growth rate we need to determine the characteristic length λ for the cracked geometry and the material parameters n, m, δ f n andε 0 , δ 0 (at a reference stress σ 0 ) or equivalently the material parameters A and B. The steady creep response under a constant uniaxial stress σ is given by Nikbin et al. (1983) at 538 • C: where n = 9 and B = 10 −23 MPa −9 h −1 .
In the remainder of the fitting process described in detail here we limit our consideration to the linear Kachanov model. Parallel procedures can be undertaken for the other damage models described in this paper. In determining the material parameters we assume that the damage zone model can also be used to describe damage development on grain boundaries in a uniaxial test. We further assume that damage grows primarily on boundaries normal to the direction of the applied stress. Integrating the damage growth rate equation between the limits ω = 0 at time t = 0 and ω = 1 at failure, i.e. when t = t f , then gives (see "Appendix B") Creep rupture data given by Nikbin et al. (1983) is plotted in Fig. 14, which gives m = 9, D = 2.7 × 10 22 MPa 9 · h, thus providing a relationship between two of the material parameters where A is measured in units of mm/(MPa 9 · h). In the second step, we determine another relationship between δ f n and A from fitting the creep crack growth data (ȧ vs C * ) to the model in Eq. (31). This combined with Eq. (58) provides two equations in terms of the two unknowns δ f n and A. To do this, we must represent Eq. (31) in terms of the fitting parameters δ f n and A. To do this we also need to determine the geometric length scale for the compact tension specimen employed in the crack growth studies. This requires the identification of an expression for C * . Here we employ an expression employed in the UK R5 assessment procedure (Ainsworth et al. 1987), which is equivalent in form to the relationship derived for the double cantilever beam (Eq. 21), i.e.
whereε 0 is the uniaxial strain rate at a reference stress σ 0 (Ainsworth et al. 1987). 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 component is defined by λ = K 2 I /σ 2 0 where K I is the stress intensity factor for the specimen at the applied load P. For a compact tension specimen the limit load for the case of plane stress (Miller and Ainsworth 1989) is given by where Υ L is a shape function, γ = 1.155, a is the crack length and W is the width of the specimen. The mode I stress intensity factor is given by where the shape function Υ K is given by Hence, the characteristic length scale is given by Miller and Ainsworth (1989) compared the predictions of Eq. (59) with detailed finite element calculations and suggested a modification to this expression to provide a better agreement with the computational results: where F p is a dimensionless parameter which is in the range 0.92 to 0.96 for n = 9 and a/W in the range 0.25 to 0.5. Here we use an average value of 0.94. Equation (65) effectively reduces the reference stress by a factor F p , but does not change the expression for the characteristic length. For the CT specimen tested by Nikbin et al. (1983) W = 50 mm and the initial crack length was a 0 = 12.5 mm. The crack growth rate is plotted as a function of C * in Fig. 15, which covers a 2 orders of magnitude increase in C * over the period of stable crack growth. From Eq. (64) we find that a two orders of magnitude increase in C * corresponds to an increase of crack length from a/W to 0.25 to 0.4. From Eq. (64) we find that this corresponds to a change of characteristic length from λ/W = 1.3 to 1.0. In our evaluation of the data we take an average of these values for the characteristic length, i.e. λ/W = 1.15. In order to proceed we need to determine which expressions to use forĈ k in Eq. (31), i.e. which regions of Figs. 12a and 13 the data lies in. From the creep deformation and creep rupture data presented earlier we findε 0 = 1.24×10 −24 σ 9 0 1/h and φ 0 = 2.17×10 2 /δ f n [defined after Eq. (26)], which suggests that for typical expected critical opening displacements in metals (δ f n ∈ [10 −10 − 10 −5 ] m, i.e. of the order of the mean cavity spacing) φ 0 > 1.0. Therefore, we use the power-law relation for φ 0 > 1.0 which isĈ kl = 0.40 · φ −0.494 0 . In this regime, Eq. (31) can be written in the form  (13) for the micromechanical model ‡ 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. (60), multiplied by F p = 0.94, to take into account the correction of Miller and Ainsworth (1989) whereȧ is in mm/h, δ f n is in mm and C * is in MJ/mm2/h. Fitting this expression to the data of Fig. 15 gives the critical separation δ f n . The rate parameter B can then be determined from Eq. (58). We can employ the same fitting procedure for the exponential Kachanov and micromechanical models. In these models additional parameters are required, i.e. β for the exponential Kachanov model and h 0 , f 0 and f c for the micromechanical model. For φ 0 > 1.0,Ĉ ke = 0.30 · φ −0.515 0 andĈ m = 0.49 · φ −0.373 0 for the exponential Kachanov and micromechanical models. Using a nonlinear least squares method, the crack growth rate in Eq. (66) is then fitted to the creep crack growth data. Figure 15 shows the comparison between experimental creep crack growth data and the framework predictions for different damage models (they all lie on the same straight line and yield comparable goodness of fit). The fitting parameters are shown in Table 1 for the different models. The parameters indicate that the material data falls in regime close to the stiff limit, i.e. φ 0 ∈ [10 3 − 10 4 ]. Further, the failure separation falls within the physical regime, i.e.
In the above analysis we noted that the characteristic length scale only changes by a small amount during the course of an experiment (i.e. by less than 13% from the mean value), while C * changes by over 2 orders of magnitude, thus the effect of C * swamps that due to λ and we could assume a constant value of λ when fitting the data. We need to be careful, however when using the developed model to assess the growth of defects in high temperature components. In general, any defects of interest will be much smaller than that employed in laboratory experiments, such as the CT specimens considered here. Conventional models of creep crack growth determined from experimental data do not include the effect of λ on the crack growth rate, other than how it affects C * (i.e. experimental laws generally assume thatȧ ∝ C * q , where q < 1). This has two major consequences. Firstly, the models developed from the data give a crack growth rate that is proportional to λ p , where p is of the order of 0.5. Thus as the crack size is reduced the crack growth rate reduces compared with that determined from the laboratory data for the same value of C * . Thus use of laboratory data overestimates the crack growth rate. Secondly, as λ is reduced φ 0 also reduces and damage growth ahead of a crack tip becomes more constrained, with the loading condition moving to the left on Fig. 12 and towards the top right corner on Fig. 13. This can result in a transition to the low φ 0 regime where now p ≈ 1.0. The crack growth rate is now more sensitive to the size of the defect and the laboratory data provides an even more significant overestimation of the crack growth rate.

Concluding remarks
In this study, a combined theoretical and computational study for crack growth under steady state conditions is presented. A theoretical framework is introduced in which the constitutive behaviour of the bulk material is described by power-law creep. A new class of damage zone models is proposed to model the fracture process such that the constitutive relation is described by a traction-separation rate law. Further, three different damage zone models are investigated: a simple critical displacement model; empirical Kachanov damage type models; and a micromechanical model. The path independence of the C * -integral is used to relate the far field loading to the damage growth process along a narrow zone directly ahead of the growing crack tip. A double cantilever beam specimen (DCB) subjected to constant pure bending moment is studied and analytical models are developed for pure mode-I steady-state crack growth. Further, using dimensional analysis the dimensionless functions C s andĈ k are derived which account for the detailed form of the model on the crack growth process. A computational framework is then implemented using the Finite Element method and the analytical models are calibrated against detailed Finite Element simulations, allowing the C s andĈ k functions to be determined for different combinations of dimensionless parameters. TheĈ k function depends on the ratio of geometric to material length scales for the problem, φ 0 , and the rate sensitivity exponents n and m. The C s -function is found to take two different power law forms over a wide range of values of φ 0 . Two simple models are presented by consider the response as φ 0 → 0 and φ 0 → ∞, which bound the computational results. The crack velocity asymptotes to the φ 0 → ∞ limit as φ 0 increases. Incorporation of the effect of damage into the constitutive response for the damage zone, results in a more substantial relaxation of stress ahead of the crack tip and a faster crack growth rate for a given value of φ 0 and a delay in the asymptote to the φ 0 → ∞ limit. Comparison with experimental data for a 2.25 Cr Mo steel at 538 • C shows that the different damage models are capable of fitting creep crack growth data using physically reasonable parameters. In this paper we have used simple material models for the matrix and damage zone response to analyses crack growth in a simple geometry. Nonetheless, the calculations provide a useful framework for evaluating the effect of different damaging processes on crack growth and evaluating data for more complex laboratory specimens.
x 1 whereū a is the displacement jumps of the middle surface which is related to the nodal displacement bȳ The finite element formulation is developed from principle of virtual work where the nodal forces and the stiffness matrix are computed. Hence, the principle of virtual work for the interface zone is defined as where Γ int is the interface surface, T is the traction vector across the interface and δδ is the virtual separation that are defined in the local coordinate system. The internal nodal force across the interface surface is determined from the tractions and by substitution of the finite element interpolation into Eq. (A.10) as The consistent tangent stiffness matrix is determined from the linearization of Eq. (A.11), with the result It should be noted that the nodal forces and the consistent tangent matrix are defined in the local coordinate system (t, n) and a transformation to the global system (1, 2) is needed before assembling to the finite element equations. This transformation includes mapping from the local to the global coordinates and relating the local to the global nodal numbering. where Δt is the time increment, T n i is the normal traction at the start of the increment and Θ is integration parameter (1 ≥ Θ ≥ 0).
In order to solve these equations, we use encapsulate the equations in a vector H = [H t H n ] T and define the unknown vector ΔT = [ΔT t ΔT n ] T . Using Newton-Raphson algorithm, the following linearized form of these equations is iteratively solved for the traction increment until a convergence condition is achieved where 0 is a zeros vector (2 × 1) and the Jacobian J is defined as (A.20) A.2 The cohesive zone Jacobian The definition of the cohesive zone Jacobian in Eq. (A.12) yields interface tangent modulus as