On the Kachanov-Rabotnov continuum damage model

In this paper two partially complementary formulations of the simple phenomenological Kachanov-Rabotnov continuum damage constitutive model are presented. The models are based on a consistent thermodynamic formulation using proper expressions for the Helmholtz free energy or its complementary form of the dissipation potential. Basic features of the models are discussed and the behaviour in tensile test and creep problems is demonstrated.


Introduction
To model continuous degradation of a material Kachanov introduced in 1958 a formulation where evolution of a single internal variable continuously reduces the elastic properties [8]. Physically such variable, which he called damage index or integrity, can be interpreted as a ratio of the differential intact area element to the original area element. For the evolution of the integrity φ, he proposed the following kinetic laẇ where the superimposed dot denotes time rate and A, n are material parameters which can depend on e.g. temperature. For an undamaged material φ = 1 and during the damaging process it decreases monotonically to the value 0 in the fully damaged state. The ratio σ/φ he called the effective stress, which is the net stress acting on the undamaged area. Kachanov used his theory in predicting creep failure times, see also [9]. Rabotnov [22] generalized Kachanov's evolution equation (1) to the forṁ where k is an additional material parameter. Since then, continuum damage mechanics has developed into an important and active field of continuum mechanics exemplified by a myriad of scientific articles and numerous books, e.g. [5,12,13,16,27].

Thermodynamic formulation
In this paper two consistent thermodynamic descriptions to derive the constitutive behaviour for an elastic damaging material, closely related to the original Kachanov-Rabotnov model, are described. In both descriptions reversible behaviour is obtained from free energy functions and irreversible damage processes are captured via dissipation potential functions. In the first model the specific free energy per unit mass is assumed to be a function of the strain tensor ε and damage D as Damage is related to integrity as D = 1 − φ and it can be viewed as a ratio of damaged differential area to original differential area, see e.g. [13,Section 1.2]. In the second model stress σ and damage α are used as the independent state variables of the free energy function 2 ψ c = ψ c (σ, α).
For clarity different symbols for the damage variable are used in the models. These two free energies (3) and (4) are related to each other by the partial Legendre transformation as 3 ρψ + ρψ c = σ : ε, where ρ is the density of the material. The dissipation potential functions for the two models are denoted as which are expressed in terms of the thermodynamic forces Y and Z dual to the fluxesḊ andα, respectively. The dissipation potentials are associated with the power of dissipation γ, such that Convexity is not a necessary property of the dissipation potentials but the condition that the products (∂ϕ/∂Y )Y or (∂ϕ c /∂Z)Z are non-negative is essential for the fulfillment of the dissipation inequality. For every thermodynamically admissible process the Clausius-Duhem inequality has to be satisfied. In the absence of thermal effects it can be expressed as written in terms of the free energy (3). An alternative expression is obtained for the free energy function (4) by differentiating equation (5) by parts with respect to time and assuming constant density, which results in and substituting it into the expression of dissipation power (9) gives Considering definition (3) in the dissipation power (9) gives and using (4) in (9) results in the counterpart expression The general forms of the constitutive equations are obtained when definition (7) 1 is equated with (12) and (7) 2 with (13). In particular, defining that results in and Since these equations have to be satisfied for all thermodynamically admissible processeṡ ε, Y,σ and Z, the general constitutive equations derived from the potential functions (3), (4) and (6) have the forms and

Particular models
Two different partially complementary simple formulations to describe elastic damaging behaviour are given. The first formulation is based on the free energy (3) and the dissipation potential (6) 1 , and the particular forms studied are given as where C is the elasticity tensor. Parameters t d , r and p are associated with the damage evolution, whereas Y r is a reference value of the thermodynamic force Y that can be chosen freely. Naturally, the value chosen for Y r influences the value of characteristic damage time t d . The particular choice used for Y r is given below.
Corresponding expressions for the second formulation are based on the free energy (4) and the dissipation potential (6) 2 , and written as where parameters t c d , n and m are associated with the damage evolution. Again Z r is a reference value of the thermodynamic force Z and can be chosen freely.
Making use of eqs. (17) and choices (19) yields the following constitutive equations For the second model use of eqs. (18) and choices (20) results in the following constitutive equations Furthermore, the thermodynamic forces Y and Z can be written as A reasonable choice for the reference values Y r and Z r is where σ r = Eε r is a reference stress, and ε r is a corresponding reference strain.

Model 1
Now the problem is solved in the uniaxial constant strain rate tensile test ε(t) =ε 0 t. The evolution equation for the integrity variable φ has to be integrated from the undamaged state, φ = 1 to a damaged state 0 < φ < 1, or even to the fully damaged state φ = 0. Using (32) yields while in the case p = −1, the result is Note thatε 0 t d is dimensionless. For the case k = 1, the stress-strain relation is Observe that a limiting strain ε frac exists for p > −1, i.e.the tensile specimen will break when ε frac .
In the following figures the behaviour of the model is investigated by varying the parameters. The ultimate tensile stress, i.e. the fracture stress σ frac can be found to occur at the strain In figure 1(a) the parameter r is varied while keeping the other parameters p and t d fixed. Increasing the r-parameter increases the ultimate tensile strength, however, it also increases the "brittleness".
In figure 1(b) the parameter p is varied while keeping the other parameters r and t d fixed. Increasing the p-parameter decreases the ultimate tensile strength, however, it also increases the "brittleness". It can be seen that if p < −1 the model shows terminal phase ductility, thus σ → 0 when ε → ∞.
If the loading rate is increased and the other parameters are constant, the behaviour is similar but the ultimate stress is increasing with increasing loading rate, see figure 2.
Dependency of the fracture stress σ frac on the parameters is shown in figures 3, 4. Note that a logarithmic scale is used for the vertical axes in figure 4.

Model 2
For the second model (20) integration of the constitutive equation results in If m = 2n + 1, then the solution for the damage α is It is easy to observe by comparing results (35) and (39) that in this particular case the two models yield identical stress-strain responses if n = r, m = p + 2n + 2 and t c d = t d .
However, it should be noticed that the damage variables are not identical, i.e. D = α.

Creep test
Next, the models are analysed in a creep test where the stress is kept constant, σ = ζσ r , where ζ is a non-dimensional constant. Therefore it is advisable to use the stress form of the thermodynamic force (30). Then for the first model Integration results in The fracture time t frac is obtained when φ = 0, resulting in The effect of the r-parameter is shown in figure 5 in a double logarithmic scale. For the second model the damage variable α grows unbounded at time Creep strain results are shown in figure 6 for ζ = 0.5 and relations (41) between the model parameters are used. Similar to the case of constant strain-rate straining, the two models result in equivalent behaviour.

On parameter estimation
Since calibration of elasticity parameters poses no problem, only the determination of parameters related to the damage evolution is discussed. There are three parameters r, p and t d to be calibrated. However, the p-parameter influences only the material's postpeak behaviour and near the region of complete failure. Thus the parameter p is chosen in advance based purely on computational convenience. The remaining two "real" parameters r and t d can be determined from two tensile/compression tests performed with different strain-rates. Denotingε 01 andε 02 the two test strain-rates and σ frac,1 , σ frac,2 the corresponding fracture stresses, from (37) it is found that The parameter t d is then obtained from either of the failure tests as , i = 1 or 2, where β = 2r + 1 2r + p + 2 2r+p+2 (p+1)(2r+1)

Dispersion analysis
Dispersion is the phenomenon that harmonic waves with different wave lengths or frequencies propagate with different velocities. The ability to transform the shape of waves is a necessary though not sufficient condition for continua to capture localisation phenomena [25]. In a classical strain-softening continua, the waves are not dispersive, which means that the continuum is not able to transform propagating waves into stationary localisation waves [24]. In dispersion analysis, a single linear harmonic wave is considered and the displacement field u for an infinitely long 1-D continuum has the form in which U is the amplitude, k is the wave number and ω is the angular frequency.
In [2] dispersion analysis of the elastic-damaging material was done for the Kachanov-Rabotnov model in the special case p = 1. Here, the analysis is generalized to arbitrary positive values of the p-parameter.
The equation of motion for a uniform bar is where ρ is the mass density of the material. For the dispersion analysis, the equations (49) and (21)-(32) are written in a non-dimensional form by defining the following nondimensional quantities: where L a length parameter, e.g. the wavelength corresponding to the lowest natural frequency and c e is the speed of an elastic wave. In addition, it is convenient to define the relative strain e = ε/ε r .
The following non-dimensional time is also used Using the non-dimensional quantities, the equation of motion (49) takes the form In the sequel, the superimposed dot represents the derivative w.r.t. the non-dimensional time τ and the prime the derivative w.r.t. the non-dimensional spatial coordinate ξ. If needed, the superimposed circle will denote the derivative w.r.t. real time, i.e.
The constitutive equations (21)-(32) take the form and the non-dimensional form of the wave is whereŪ = U/L,k = k/L andω = ωt e . In the following, the strain form of the damage evolution equation is chosen, i.e.
Linearisation of the divergence of the stress at state φ * , s * , e * has the form where φ ′ is obtained from the linearised equations of motion For the dispersion relation the rate form of the linearised equation of motion reads is needed. Divergence of the integrity rate is obtained aṡ where Taking into account definitions (59)-(65) finally results in Substituting the waveform (58) into equation (66) gives This equation can only be satisfied if the wavenumber is complex, i.e.k =k r +ᾱi, which means that the harmonic wave is attenuated exponentially when traversing through the bar asū Substituting the expression for a damped harmonic wave (68) into (66) yields i ω 3 − φ * ω (k 2 r −ᾱ 2 ) + 2(2r + p)g * krᾱ + g φω 2 + 2φ * ωkrᾱ + (2r + p)g * (k 2 r −ᾱ 2 ) = 0. (69) Since both the real and imaginary part of this expression must vanish, it will result in two equations from which the wavenumberk r and the damping coefficientᾱ can be solved: g φω 2 + 2φ * ωkrᾱ + (2r + p)g * (k 2 r −ᾱ 2 ) = 0.
Multiplying equation (70) by (2r + p)g * and (71) by φ * ω and summing the resulting equations by parts enables elimination of theᾱ 2 -term, resulting in Substituting this back into (70) gives where Clearly the last term in (73) is positive, but the sign of a 2 is positive only if This is true if r > 0 and p > 0. Thus, the solution of (73) can be written as k 2 r = 1 2 |a 2 | sign(a 2 ) + 1 + 4(a 0 /a 2 ) 2 .
The limit for the damping factorᾱ is In fig. 7 the phase velocity c and the group velocity c R are plotted as a function of the angular frequency ω with different values of the p-parameter. As can be noticed in the figures the bar shows both normal dispersion, i.e. c > c R , and anomalous dispersion behaviour c < c R [1].
In fig. 8 the decay parameter α is plotted as a function of the angular frequency ω for different values of the p-parameter.

Finite element analysis
For strain-softening inviscid solids localisation takes place in a plane of zero thickness. Viscosity added to either plasticity or damage models may bring in the desired non-zero material length-scale. To investigate the regularizing behaviour of the Kachanov-Rabotnov continuum damage model, a one-dimensional finite element analysis is carried out. A bar fixed at x = 0 and subjected at x = L to either a linearly increasing stress σ(t) = ησ r t/t d or a linearly increasing displacement u(L, t) = ηǫ r Lt/t d is analysed using different uniform mesh sizes.  If figures 9 and 10 damage profiles analysed by using 10, 100, 1000 and 10000 elements are shown for the models (19) and (20), respectively. The model parameters are related as in (41) with p = 1, r = 2 and the loading rate is defined with η = 0.01. A standard central difference scheme is used to integrate the equations of motion with a constant time-step equal to the critical time step of the elastic bar.
Fracture times are tabulated in table 1. Results for both models are almost equivalent. In model (20) the damage variable α has no upper bound and thus time to fracture is declared when strain exceeds the value 2ε r . Convergence of the fracture time is fast; the result with the coarsest mesh having ten elements differs less than 2% from the result obtained with the finest mesh. It is clearly seen that the model regularises the strainsoftening problem and there exists a clear material length scale.
From the numerical results, it can be concluded that the width of the localisation zone which could also be deduced from (33)-(35). Wang studied localisation of strain softening viscoplastic solids [26]. In the case of pure shearing of strain softening von Mises viscoplastic solid, he obtained the following  approximation of the material length scale where G is the shear modulus, c g = G/ρ the elastic shear wave velocity, h the softening modulus and η * the viscosity. He did not consider the material length scale dependency on the applied strain rate. However, it is known that decreasing the strain-rate shortens the material length scale [3,11,17]. For the elastic damaging Kachanov-Rabotnov type solid the material length scale is a function of the strain rate, the elastic wave speed and the parameters r, p and t d , i.e.
A numerical localisation study is performed with varying prescribed rate η and the damage localisation width is defined as the measure of the domain Ω loc defined as where D * is the damage value at fracture stress for quasi-static constant strain-rate loading and it is independent of the applied strain-rate. The width of the localisation zone is shown as a function of the loading rate η in figure 12 for three different mesh sizes h = L/100, L/1000 and L/10000 and for the cases r = 2, p = 1, and r = 4, p = 1. As it can be seen, the width of the damage localisation zone is mesh-size independent if the loading rate satisfies η 0.75.
In figure 13 the damage localisation width is shown in a semi-logarithmic plot as a function of the parameter r. As it can be seen, the localisation width converges with respect to the mesh-size and the relationship is almost linear in the logarithmic scale.

Concluding remarks
Two different formulations of the well-known Kachanov-Rabotnov continuum damage model have been investigated using two simple uniaxial loading cases, namely the constant strain rate tensile test and creep test. Moreover, dispersion analysis has been carried out. It is shown that the two models yield identical results if certain simple algebraic relations between the model parameters are satisfied. However, the models are not identical since the damage parameters in the models have different meanings. Numerical computations have also revealed that the models have intrinsic length-scales and deformation and damage localise into a zone with finite length. However, in dynamic case the models predict wave speeds higher than the elastic one, which is physically questionable. Furthermore, the damage localisation depends on the applied rate of loading and the finite element computations show that for slow loading rates the computations are mesh-size independent. Thus, there exists a certain threshold after which faster loadings result in mesh dependent solutions.
It is well known that using an integral type non-local constitutive equation or inclusion of higher order spatial displacement gradients, as well as damage gradients, brings the desired length-scale into the material model, see e.g. [6,15,18,19,20,21,23]. However, the damage evolution equation is then a partial differential equation and extra boundary conditions are required. It also adds unknowns to the global equation system and increases the computational burden. An interesting topic for further research could be an investigation to the effect of higher order time rates to the damage evolution equation.