Skip to main content
  • Research Article
  • Open access
  • Published:

Coupling local and non-local damage evolutions with the Thick Level Set model

Abstract

Background

The Thick Level Set model (TLS) is a recent method to delocalize local constitutive models suffering spurious localization. It has two major advantages compared to other delocalization methods. The first one is that the transition from localization to fracture is taken into account in the model. The second one is that the delocalization only acts when and where needed. In other words, the TLS has no effect when the local model is stable. The former advantage was already detailed in several papers (IJNME 86:358-380, 2011, CMAME 233:11-27, 2012, IJF 174:49-60, 2012). This paper concentrates on the latter advantage.

Methods

The TLS delocalization approach is formulated as a bound on the damage gradient. The non-local zone is defined as the zone where the bound is met whereas the local zone is defined as the zone where it is not met. The boundary (localization front) between the local and non-local zone is the main unknown in the problem.

Results

Based on the new model, a 1D pull-out test is solved both analytically and numerically. Different regimes are observed in the solution as the loading progresses: fully elastic, local damage, coupled local/non-local damage and, finally, purely non-local damage.

Conclusions

The new model introduces delocalization as an inequality allowing local damage to develop in zones whereas non-local damage may develop in other zones. This reduces dramatically the cost of implementation of such models compared to fully non-local models.

Background

Although the scope of TLS application is much wider, we consider in this paper the fracture of quasi-brittle structures under quasi-static loading and under small deformation assumption. The loading is proportional to a scalar parameter. The material is modelled by a time-independent elasto-damage constitutive model with scalar damage. Due to quasi-static analysis, the loading parameter must be controlled especially when bifurcation occurs.

The TLS model was introduced in several papers [1],[2] and [3],[4]. It lies between continuum damage mechanics and fracture mechanics. Indeed, crack opening is allowed across fully damaged zones (see [2] for instance). The fully damaged zone is located by a level set. Let us note that the description above is different from a diffuse vision of the crack in which crack opening is not explicitly modeled as in the phase-field approach [5]-[7] or the variational approach to fracture [8],[9]. We are rather in the vein of transition from damage to fracture as in [10]. However, the TLS will not be in need of a cohesive zone to perform the transition. The model can be considered as a continuous transition from damage to fracture.

The main idea of the TLS for quasi-brittle fracture is to bound the spatial gradient of the damage variable d, thus avoiding spurious localization. One imposes that the spatial damage distribution satisfies at all time

df(d)onΩ
(1)

where Ω is the domain of interest. The choice of the function f(d) will be discussed in what follows. As damage evolves, one eventually wants to locate the crack, i.e. the zone for which d=1. However, finding the iso-contour d=1 for a quantity d than cannot go beyond 1 is a tedious operation. This is where the level set ingredient comes into play. Variable d is expressed in terms of a level set ϕ as depicted in Figure 1. This relation introduces a length scale lc. Finding the zone d=1, is now well-posed since the level set ϕ is not strictly limited to lc but may go beyond. With the use of the surrogate variable ϕ, condition (1) may be rewritten as

∇ϕ 1 d = d ( ϕ )
(2)
Figure 1
figure 1

An example of damage shape function.

where f(d) in (1) is related to d(ϕ) by f(d)=d(ϕ(d)) (the prime indicating the derivative of d with respect to ϕ). The function d(ϕ) is called the damage shape function and is the main ingredient of the TLS. Equation (2) above indicates that ϕ is a distance function in the zone where the constraint is active (we name this zone the localization zone). The evolution of a distance function has been analyzed and updating algorithm proposed in [11]. In the localization zone, the evolution of ϕ is non-local, indeed

ϕ=1 ϕ . ·ϕ=0
(3)

The rate of change of ϕ is thus uniform on any segment aligned with ϕ and the rate of d is given by d . = d (ϕ) ϕ . . Such segments over which ϕ . is uniform are depicted in Figure 2. In the local zone, the evolution of ϕ stems from the evolution of d and the relation d=d(ϕ).

Figure 2
figure 2

Local ( Ω) and localization ( Ω+) domains as well as cracks faces ( Γ c ). Inside Ω+, we have ║ϕ║=1. It can be noted that ϕ (and thus d) is discontinuous along the dashed line which is the so-called skeleton of the distance function.

The delocalization (1) used in the TLS is different from existing delocalization techniques. Indeed, it directly uses the norm of the damage gradient. It is thus a Hamilton-Jacobi type equation. On the contrary, damage gradient models [12]-[14] yield Laplacian damage type equation rising the question of proper boundary conditions.

The TLS shares some similarities with the so-called non-local integral approach [15],[16] in which weighted averages are performed over segments (1D), disks (2D) and spheres (3D) of fixed size. In the TLS approach, however, weighted averages are always performed on segments (Figure 2) whatever the dimension of the body and over a length which is not fixed in time but evolves from zero to a maximum length lc. Finally, note that as lc is the minimal distance between a point where d=0 and a fully damaged point, d=1, it plays the role of the fracture process zone size.

After this quick introduction of the TLS, we get to the objective of the paper. In previous TLS paper, the delocalization condition (2) was considered as an equality on the whole domain. It meant that d was zero on the domain except in zone where the gradient norm was fixed. The short-coming of this view was that uniform or smooth damage field (because of damage hardening for instance) could not be modeled prior to localization. The inequality analyzed in this paper allows a combination of local and non-local evolutions. In the literature, the possibility to combine both local and non-local approach is seldom discussed with the exception of the so-called morphing numerical technique [17],[18].

The paper is organized as follows. The TLS concept with the inequality constraint discussed above are detailed in the first section. Next, the TLS boundary value problem is set up and a dissipation analysis is carried out. A 1D pull-out is solved semi-analytically to show the main feature of the TLS solution. This 1D test is then solved numerically with the TLS to observe the influence of the parameters choice in the model. A conclusion and perspectives end the paper.

Methods

We consider a solid body occupying a domain Ω. The external surface Ω is composed of two parts ∂Ω u u and ∂Ω u T on which the displacements λ ud and the loading λ Td are prescribed, respectively. The parameter λ is a loading parameter.

Small strains and displacements are assumed as well as quasi-static evolution. The current state is characterized by the displacement field u, from which the strain field 𝛜= 1 2 ( u+ u T ) is derived. The current state is also characterized by an internal scalar variable, the damage denoted d. In this paper, we will not consider other internal variables.

Regarding the material model, we consider a free energy ψ(𝛜,d) from which the stress tensor σ and local energy release rate Y may be derived

σ= ∂ψ ( 𝛜 , d ) 𝛜 ,Y= ∂ψ ( 𝛜 , d ) ∂d
(4)

The potential ψ is assumed for now at least convex with respect to 𝛜. The need for other properties will be discussed later.

Regarding the time-independent damage evolution, we consider a function y depending on damage and strain history (through e) such that

d . 0,y(e,d) Y c 0,(y(e,d) Y c ) d . =0
(5)

where

e=e(𝛜(τ),τt)
(6)

and Yc is some threshold. We believe the above formalism encompasses most of the damage models in the literature. To be even more general, one may consider two relations of the kind (5): one for damage in tension and a second one for damage in compression and then combine these damages into d. One has a so-called associated damage model when the y variable is Y. In this case, damage evolution is expressed in terms of the dissipation potential φ*(Y) which is the indicator function of YYc≤0:

d . φ * ( Y ) ∂Y
(7)

Such model was already considered in [1] for dissymmetric tension-compression evolution. We emphasize the fact that the TLS description is not restricted to associated damage models.

What we have described so far is a purely local damage model. This type of model is known to suffer spurious localizations meaning that the damage gradient may become infinite. The main idea of the TLS approach is to bound damage gradient as expressed in (1). In the TLS model, damage is allowed to go to 1 (but not beyond of course). The location of a crack (or fully degraded zones like in comminution problems) is defined by the set of points for which d=1. Numerically speaking, finding the set of points for which d=1 knowing that d may not go beyond 1 is not very practical. This is why the TLS expresses damage in terms of a surrogate variable ϕ whose values are not limited as depicted in Figure 1. We assume the following regularity on d(ϕ)

d ( ϕ ) C 0 ( ] , + [ ) and monotonically increasing d ( ϕ ) = 0 if ϕ 0 d ( ϕ ) = 1 if ϕ l c d ( ϕ ) C 1 ( ] 0 , l c [ )
(8)

Finding the subdomain where d=1 is equivalent to find the subdomain whose boundary is the iso-contour ϕ=lc.

In terms of the surrogate variable, ϕ, condition (1) reads

ϕ1
(9)

provided f(d) is given by

f(d)= d (ϕ(d))
(10)

For instance, if d is linear with respect to ϕ, the gradient of damage will be bounded by a constant

d=ϕ/ l c ,ϕ[0, l c ]d 1 l c
(11)

whereas for more complex function d(ϕ), the bound depends on the level of damage. For instance, for the profile shown in Figure 1, we have

d=1 ( 1 ( ϕ / l c ) ) 2 d 2 l c 1 d
(12)

For a general power law with n≥1, we obtain

d=1 ( 1 ( ϕ / l c ) ) n d n l c ( 1 d ) 1 1 / n
(13)

Whether local or non-local constitutive model should be used at a point x is based on condition (9).

ϕ ( x ) < 1 Local constitutive model at x
(14)
ϕ ( x ) = 1 Non-Local constitutive model at x
(15)
ϕ ( x ) > 1 forbidden
(16)

The first condition is the major novelty of this paper, compared to previous paper on the TLS. At any time t, the domain may thus be decomposed into three non-overlapping zones : a local zone Ω, a non-local zone Ω+ and a fully damaged zone Ωc

Ω = Ω c Ω + Ω
(17)
Ω = { x : Ω ∇ϕ ( x ) < 1 , ϕ ( x ) < l c }
(18)
Ω + = { x Ω : ϕ ( x ) = 1 , ϕ ( x ) < l c }
(19)
Ω c = { x Ω : ϕ ( x ) l c }
(20)

We define also the boundary Γc of the fully damaged zone and the interface Γ between the local and non-local zones.

Γ c = Ω c ,Γ= Ω + Ω = Ω + Ω
(21)

The boundary Γc defines the crack faces. Figure 2 shows a typical scenario of a crack appearing inside the localization zone.

Note that the volume measure of Ωc may be zero. This information is part of the solution process. We expect different shapes of Ωc in comminution and brittle crack propagation.

Eikonal equation

Condition, ║ϕ(x)║=1 is a non-linear first-order partial differential equation. It is called an eikonal equation and belongs to the Hamilton-Jacobi equation family. Among the possible solution satisfying ║ϕ(x)║=1, we will pick the one corresponding to the vanishing viscosity solution [19]. It is characterized by

ϕ(x)= min y Γ (ϕ(y)+d(x,y)),x Ω +
(22)

where d(x,y) is the length of the shortest path connecting x and y inside Ω+. The value of ϕ at xΩ+ can be thought as the minimal fare to go from Γ to x. The fare being the sum of the initial fare ϕ(y) plus the mileage from y to x. Damage on Ω+ is thus fully determined from values on Γ. A 1D example of ϕ satisfying the eikonal on a segment [ b,d] is given in Figure 3.

Figure 3
figure 3

Distribution of ϕ on a 1D domain: Ω = [ a , e ], Ω+= [ b , d ], Ω= [ a , b ] [ d , e ], Γ ={ b , d }. Slopes at 45 degrees on Ω+ indicate that ϕ behaves as a distance function (║ϕ║=1) whereas ║ϕ║<1 on Ω. Point c is the skeleton of the distance function.

The fact that damage is related to a variable satisfying the eikonal equation, the cornerstone of the level set technology [11], explains why the damage model is coined Thick Level Set. In the non-local zone, damage is modeled over a thick layer in terms of level sets.

Damage evolution

In the local zone, Ω, damage evolution is local and given by (5). In the non-local zone, Ω+, damage rate is related to ϕ . by

d . = d (ϕ) ϕ .
(23)

where ϕ . is uniform on segments aligned with ϕ, see Equation (3). We denote this space as A:

ϕ . A= a ( x ) L 2 ( Ω + ) : a . ϕ = 0
(24)

Non-local damage evolution boils down to decomposing Ω+ into a set of independent segments and finding a value ϕ . over each of them.

As in [2], we suggest to introduce averaged quantities, y , d . over each segments. This may be expressed by a projection operation.

y A : Ω + y d a d ω = Ω + y d a d ω , a A
(25)
d . A : Ω + d . a d ω = Ω + d . a d ω , a A
(26)

We note that the averages satisfy the following property

Ω + y d . dω= Ω + y d . dω
(27)

The above indicates that duality is preserved through the averaging technique. This is not often the case in delocalization techniques as discussed in [20].

The local constitutive model, (5), is then expressed in terms of the non-local quantities

d . 0, y Y c 0, y Y c d . =0
(28)

where we have assumed Yc uniform (if not it needs to be averaged by formula (25)). Finally, we write the relation giving ϕ . in terms of d . :

d . = d ' ϕ . , d A: Ω + d adω= Ω + d adω,aA
(29)

To end this section we illustrate the average formula on the 1D example depicted in Figure 3. Averages are given by

On [ b , c ] : y = b c y d ( ϕ ) d x b c d ( ϕ ) d x , d . ( x ) = b c d . d x b c d x
(30)
On [ c , d ] : y = c d y d ( ϕ ) d x c d d ( ϕ ) d x , d . ( x ) = c d d . d x c d d x
(31)

TLS boundary value problem

We are now able to define the boundary value problem. The set of admissible displacements is given by

U= u C 0 ( Ω \ Ω c ) , Ω \ Ω c ψ ( ϵ ( u ) , d ( ϕ ) ) d ω < + , u = λ u d on Ω u
(32)

The fact that the fully damage zones are removed from the domain is important. It allows the displacement to be discontinuous across Γc. Regarding the regularity of the displacement, we request that the energy, i.e. integral of ψ over Ω\Ωc, is finite. This space is not simply H1 as in elasticity since the stiffness is possibly vanishing on Γc boundary [21].

Regarding the ϕ variable, it is required to be continuous over Ω and belong to the set . The admissible set for ϕ is denoted K.

K= ϕ C 0 ( Ω ) : ϕ ( x ) = 1 , x ( Ω + Ω c ) , ϕ ( x ) < 1 , x Ω
(33)

The continuity requirement on ϕ leads to a Hadamard compatibility condition on the moving boundary Γ. Let us define the jump of a quantity f across Γ by

[ f ] Γ (x,t)= f + (x,t) f (x,t)
(34)

The exponent −/+ placed on some quantities f defined at x on Γ has the following meaning

f ± (x,t)= lim h 0 + f(x±hn(x,t),t)
(35)

where n is the outward normal to Ω+. With these notations we have

ϕ . Γ + [ ϕ ] Γ ·n v n =0
(36)

where vn is the normal velocity of Γ counted positively along n. This gives the respective evolution of domains Ω+ and Ω.

Potential energy of the domain is given by :

u U , ϕ K E pot ( u , ϕ , λ ) = Ω ψ ( ϵ ( u ) , d ( ϕ ) ) d ω + ω + ψ ( ϵ ( u ) , d ( ϕ ) ) d ω Ω T λ T d · u d a
(37)

Note that the same free energy expression, ψ, is used over Ω+ and Ω. In what follows, n is the outward normal vector to Ω on Ω and to Ω+ on Ω+. The set of admissible displacements variations is denoted as U 0 . It has the same definition as except that u is set to zero on Ω u .

Assuming, that at time t, the spatial distribution of ϕ of the two volumes Ω+ and Ω is known, the displacement field u is the field that solves the stationarity of the potential energy:

E , u pot δu=0,δu U 0
(38)

This means;

Ω σ:ϵ(δu)dω+ Ω + σ:ϵ(δu)dω Ω T λ T d ·δuda=0,δu U 0
(39)

For simplicity, we assume that the boundary Γc is traction free (no contact on crack faces). The equilibrium (39) yields the following local equations

div σ = 0 over Ω \ Ω c
(40)
σ Γ · n = 0 on Γ
(41)
σ · n = λ T d on Ω T
(42)
σ · n = 0 on Γ c
(43)

We stress the fact that Γc denotes the boundary of the fully damage zone and thus in case of a crack Γc indicates both crack lips. To complete the set of equations to be solved for a known damage distribution, we add the stress definition and kinematic relations

σ = ∂ψ ( ϵ , d ) ϵ over Ω \ Ω c
(44)
ϵ = 1 2 u + u T over Ω \ Ω c
(45)
[ u ] Γ = 0 on Γ
(46)
u = λ u d on Ω u
(47)

Finally, we need to add damage evolution equations in the local zone (5) and non-local zone (28).

Dissipation analysis and fields regularity

The goal of this section is to analyze the expression of the dissipation as well as looking at the fields regularity across the boundary Γ.

Taking into account the conservation law for the total energy during the evolution of the system, the total dissipation associated with the loading rate λ . is:

D = Ω T λ T d · u . d a + Ω u λ . u d · σ · n d a
(48)
d d t Ω ψ ( ϵ , d ( ϕ ) ) d ω + Ω + ψ ( ϵ , d ( ϕ ) ) d ω
(49)

In the above, we did not consider the energy inside Ωc because it is assumed to be zero. Indeed, no compression is considered in this zone (see (43)).

Using Leibniz formula for the time derivative of moving domains as well as the relation:

ψ . =σ:ϵ( u . )Y d .
(50)

we obtain

D = Ω T λ T d · u . d a + Ω T λ . u d · σ · n d a Ω σ : ϵ ( u . ) d ω Ω + σ : ϵ ( u . ) d ϵ + Γ [ ψ ] Γ v n d a Γ c ψ v n d a + Ω Y d . d ω + Ω Y d . d ω

Integrating the domain integral by parts in the second line above and using the equations characterizing the equilibrium state, we get

D= Γ [ ψ ] Γ v n +n·σ· u . Γ da Γ c ψ v n da+ Ω Y d . dω+ Ω + Y d . dω
(51)

During the propagation of the interface, perfect contact is assumed on Γ, that is the displacement jump across Γ must be zero at all time. As a consequence, the derivative along the moving interface of the displacement jump must be zero, [22], yielding the so called first Hadamard compatibility condition between the front velocity vn and the jump in material velocities u . Γ :

u . Γ + u Γ ·n v n =0
(52)

Equation (51), now becomes

D= Γ n· [ P ] Γ ·n v n da Γ c ψ v n da+ Ω Y d . dω+ Ω + Y d . dω
(53)

where P is the Eshelby tensor

P=ψIσ· u
(54)

The first term is the dissipation created by the interface propagation. We show now that due to damage continuity on Γ this term is zero.

Since normal stress and displacement are continuous across Γ, the product of the jump in stress and strain across Γ is zero, [23],[24]:

σ Γ : ϵ Γ =0
(55)

Let ψd(ϵ) be the density of free energy for a given value of damage and let ψ d * (σ) be its dual by the Legendre-Fenchel transform. Since the couples (ϵ+,σ+) and (ϵ,ϵ) do satisfy the constitutive model (4), we have

ψ d ( ϵ + )+ ψ d * ( σ + ) σ + : ϵ + =0, ψ d ( ϵ )+ ψ d * ( σ ) σ : ϵ =0
(56)

Summing the two relations above and using (55), we have

ψ d ( ϵ + ) + ψ d * ( σ ) σ : ϵ + + ψ d ( ϵ ) + ψ d * ( σ + ) σ + : ϵ =0
(57)

Since both terms above are greater or equal to zero (classical property of convex analysis, see [25]), we have

ψ d ( ϵ + )+ ψ d * ( σ ) σ : ϵ + =0, ψ d ( ϵ )+ ψ d * ( σ + ) σ + : ϵ =0
(58)

This implies that the couples (ϵ+,σ) and (ϵ,σ+) do satisfy the constitutive model. Assuming that the convex potential ψd(ϵ) is such that the stress associated to any strain is unique, we have

σ Γ = ϵ Γ =0
(59)

The continuity of the strain and displacement across Γ leads to the continuity of the displacement gradient

u Γ =0
(60)

leading finally to the continuity of the Eshelby tensor.

P Γ =0
(61)

Dissipation is thus reduced to

D= Ω Y d . dω+ Ω + Y d . dω Γ c ψ v n da
(62)

The dissipation must be positive. For classical models in which Y is positive, this implies that damage may only grow. Damage growth will create a growth of the fully damaged zone Ωc (and thus a negative velocity vn). The last term in (62) is thus automatically positive. Whether this term is zero or not depends on the regularity of ψ on the boundary Γc. This regularity must be assessed from the non-local constitutive model condition: y Y c 0.

Note that dissipation may also be written

D= Ω Y d . dω+ Ω + Y d . dω Γ c ψ v n da
(63)

where Y is defined by (25) (y replaced by Y). The above expression exhibits the duality between Y and d . in the localization zone.

Results

We consider a 1D axisymmetric fiber pull-out depicted in Figure 4. The fiber of radius r i is considered rigid and infinitely long. It is pulled out of a clamped circular domain of radius r e =r i +L. The only non-zero stress component is the shear stress τ satisfying the following equilibrium conditions

( τr ) , r =0τ(r)= τ ( r i ) r i r
(64)
Figure 4
figure 4

Pull-out of an infinite fiber of radius r i from a tube of radius r e . Radius r l indicates the (evolving) extent of the non-local damage zone.

The only non-zero strain is the shear strain, derivative of the displacement along the fiber direction

γ= u , r
(65)

We consider the following free energy density involving some hardening function h(d), satisfying h(1)=0. The shear stiffness is denoted μ and Yc is also a material parameter.

ψ(γ,d)= 1 2 μ(1d) γ 2 + Y c h(d)
(66)

So, state laws read

τ=μ(1d)γ,Y= 1 2 μ γ 2 Y c h
(67)

The local evolution model is given by

d . 0,Y Y c 0,(Y Y c ) d . =0
(68)

The condition Y=Yc reduces to

τ τ c = ( 1 d ) h + 1 g ( d ) , τ c = 2 μ Y c
(69)

Let us now be more precise on the type of function g(d) we will be considering. Basically, we are interested by C1 positive concave functions with a maximum value at some damage dc<1 :

g(d) C 1 ([0,1]): g <0,g(0)=1,g(1)=0, g ( d c )=0
(70)

We shall use in what follows

g(d)=(1d)exp d 1 d c
(71)

The corresponding stress strain curve is given in Figure 5. We will now search for the complete solution linking the (non-dimensional) shear stress T needed to move by a (non-dimensional) displacement U the fiber:

T= τ ( r i ) τ c ,U= u ( r i ) μ r i τ c
(72)
Figure 5
figure 5

Local constitutive model: stress versus rising strain (case d c =0 . 5).

Four regimes will be observed. They are depicted in Figure 6: elastic, local damage, local and non-local damage and finally purely non-local damage. The first two regimes may be solved analytically whereas the two last one may not. We however pursue as much as possible the analytical path. Next section is devoted to a numerical solver.

Figure 6
figure 6

Force-displacement curve in the case n = 2. The elastic part of the solution is on part [OA] of the curve. The part [AB] corresponds to the development of local damage whereas part [BC] corresponds to both local and non-local damage development. Finally, part [CD] is governed only by non-local damage. Letters A, B and C are located at T A =1,T B =1.351 and T C =1.368, respectively. The latter being the limit load.

Pure elastic regime: T [ 0,TA=1]

The displacement solution is given by

u(r)=T τ c μ r r e r i r dr=T τ c μ r i log r e r
(73)

We thus have a linear relationship between the stress and displacement

T=U log r e r i 1
(74)

Local damage regime: T [ 1,TB]

When T reaches 1 local damage starts around the fiber. Its distribution is obtained by combining (64) and (69)

T=g(d) r r i
(75)

This distribution of local damage is acceptable provided the condition below holds true

ϕ1i.e.d d (ϕ(d))
(76)

The norm of the damage gradient is maximum at r=r i and of value

d d d r = 1 r i g ( d i ) g ( d i )
(77)

where d i =d(r i ). Considering a general power law damage profile (13), the condition (76) is

g ( d i ) g ( d i ) r i l c n ( 1 d i ) 1 1 / n
(78)

Let us denote by d i B the smallest value of damage for which the condition above is violated and T B the corresponding loading. Due to the fact that g(dc)=0, it is clear that d i B will be slightly lower than dc. We note that as the material length gets bigger with respect to r i , non-locality (violation of (78)) will step in for smaller and smaller damage d i . Considering the choice (71), we get the condition

( 1 d i ) ( 1 d c ) d c d i r i l c n ( 1 d i ) 1 1 / n
(79)

For instance if n=1, we get

d i B = d c ( 1 d c ) ( l c / r i ) 1 ( 1 d c ) ( l c / r i )
(80)

As a numerical application, with dc=0.5 and lc/r i =0.1, we get d i B =0.47.

Combined local and non-local damage regime: T [ TB,TC]

For a loading higher than T B , non-local damage will develop close to the fiber. Let [ r i ,r l ] be the current extension of the non-local damage zone in which damage ranges from d i to d l following:

r(d, d i )= r i +ϕ( d i )ϕ(d)
(81)

The condition for the non-local zone to grow is Y = Y c , i.e.

r i r l (Y Y c ) d d d ϕ rdr=0
(82)

Using, (81), we may rewrite it as

d l d i (Y Y c )r(d, d i )dd=0
(83)

Now, using (69), we get

T= d l d i ( 1 d ) 2 g 2 ( d ) ( r ( d , d i ) / r i ) d d d l d i ( 1 d ) 2 ( r ( d , d i ) / r i ) 1 d d 1 / 2
(84)

Since loading is rising, so does local d l damage at r=r l , following

T= r l r i g( d l )
(85)

Given T, system (84)-(85) returns unknowns d i and d l as well as the extent of the non-local zone r l =r(d l ,d i ). We note that for T=T B , we have d l =d i and r l =r i . Let T C be the loading above which the system has no solution.

Non-local damage regime: T decreases from TCto 0

There is no solution of the problem for a loading higher than T C . When the loading decreases below T C , there is of course a possible elastic solution. Another possible solution is the further development of the non-local damage zone (while damage in the local zone no longer evolves since loading is decreasing). The system of equations to solve still involves (84)

T= d l d i ( 1 d ) 2 g 2 ( d ) ( r ( d , d i ) / r i ) d d d l d i ( 1 d ) 2 ( r ( d , d i ) / r i ) 1 d d 1 / 2
(86)

Equation (85) is now different and reads

T C = r l r i f( d l )
(87)

Indeed the damage at r l did not change from its value at load T C because the load has been decreasing afterwards.

Analysis of the displacement of the fiber

The displacement of the fiber is given by

U=T r i r e 1 ( 1 d ) r dr
(88)

As the damage around the fiber goes to 1, the integrand goes to infinity. But, at the same time the loading goes to zero. Let us study the limit of the fiber displacement for the loading going to zero. The loading is given by (84) recalled below

T= d l d i ( 1 d ) 2 g 2 ( d ) ( r ( d , d i ) / r i ) d d d l d i ( 1 d ) 2 ( r ( d , d i ) / r i ) 1 d d 1 / 2 = d l d i N ( d , d i ) d d d l d i D ( d , d i ) d d 1 / 2
(89)

Due to the property of g(d), (70), we have

g(d)=O(1d)asd1
(90)

So

0<N(d, d i )<+,d, d i [0,1]
(91)

Finally, we have

T=O 1 d i as d i 1
(92)

Note that this property does not depend on the choice of d(ϕ). Going back to the displacement expression, (88), we have

U=T d l d i ( 1 d ) 1 ( r ( d , d i ) ) 1 d r ( d , d i ) d d dd+CT
(93)

where C is a finite constant and

r(d, d i )= r i +ϕ( d i )ϕ(d)= r i + l c ( 1 d ) 1 / n ( 1 d i ) 1 / n
(94)

assuming a power law asymptotic behavior of ϕ(d) as d goes to 1:

ϕ(d)= l c (1 ( 1 d ) 1 / n )
(95)

Finally, we get

U=TO ( 1 d i ) 1 / n 1 + C =O ( 1 d i ) 1 / n 1 / 2 as d i 1
(96)

We conclude that there exists three regimes of delocalization.

lim d i 1 U = 0 if n < 2
(97)
0 < lim d i 1 U < + if n = 2
(98)
lim d i 1 U = + if n > 2
(99)

When n<2, the fiber displacement must be zero for total failure. When n=2, there exists a limit value of fiber displacement before total failure and when n>2, it takes an infinite displacement before total failure. It is interesting to note that these three regimes also exist in gradient damage models [20].

Numerical solve

Last section gave some insight on the different regimes in the solution. In order to plot the solution, we detail here a 1D numerical solver. This code is rather ad hoc for 1D problem, since we force the advance of the Γ boundary and find the corresponding loading and fields. General 2D and 3D solvers will be detailed in a forthcoming paper. We search for the solution at a set of discrete times. Consider the solution known at time t n , the solution at time tn+1 must satisfy the following equations.

Kinematics and equilibrium on ] r i , r e [.

u n + 1 U = u H 1 ( ] r i , r e [ ) : u ( r e ) = 0
(100)
γ u n + 1 = u , r n + 1
(101)
r i r e τ n + 1 γ ( u * ) r d r = τ c T n + 1 u * ( r i ) r i , u * U
(102)

State laws and d ( ϕ ) relation on ] r i , r e [.

τ n + 1 = τ γ u n + 1 , d n + 1 = 1 d n + 1 μγ u n + 1
(103)
Y n + 1 = Y γ u n + 1 , d n + 1 = 1 2 μγ u n + 1 2 Y c h d n + 1
(104)
d n + 1 = d ϕ n + 1
(105)

Non-local evolution law on ] r i , r l n + 1 [ .

a = r i r l n + 1 Y n + 1 Y c d ϕ n + 1 r d r 0 ,
(106)
b = ϕ n + 1 ( r i ) ϕ n ( r i ) 0 , ab = 0
(107)
ϕ n + 1 = ϕ n + 1 r l n + 1 + r l n + 1 r , on ] r i , r l n + 1 [
(108)

Local evolution law on ] r l n + 1 , r e [ .

Y n + 1 Y c 0, ϕ n + 1 ϕ n 0, Y n + 1 Y c ϕ n + 1 ϕ n =0
(109)

Regarding space discretization, the segment ]r i ,r e [ is discretized with a set of finite elements. Initially, the non-local zone is empty and we proceed with a classical Newton-Raphson scheme depicted in the solver flowchart without non-local zone.

Solver flowchart without non-local zone

  1. 1.

    initialization: u 0=d 0=T 0=0

  2. 2.

    elastic step: find the load step for which damage starts

  3. 3.

    load step T n+1=T n+Δ T

  4. 4.

    iterations initialization k=0

  5. 5.

    solve linear system (110) to find Δ u

  6. 6.

    after the first iteration adapt the load step so that the maximum damage increment is d inc

  7. 7.

    update using (111)-(118)

  8. 8.

    if residual ≤ tol, go to 9 else go to 5

  9. 9.

    if ║ϕ n+1║≤1, go to 3, else go to solver flowchart with non-local zone

The linear problem at each iteration reads: find ΔuU such that:

r i r e H k γ(Δu)γ( u * )rdr= τ c T n + 1 u * ( r i ) r i r i r e τ k γ( u * )rdr, u * U
(110)

where the right hand side is the residual at iteration k. Once the displacement correction is obtained, the local update of the fields is computed from

u k + 1 = u k + Δu
(111)
τ k + 1 = τ γ u k + 1 , d k + 1 , Y k + 1 = Y γ u k + 1 , d k + 1
(112)
Y k + 1 Y c 0 , d k + 1 d n 0 , Y k Y c d k d n = 0
(113)
ϕ k + 1 = ϕ d k + 1
(114)

whereas tangent operators are obtained by

H k + 1 = H γγ k + 1 η k H γd k + 1 H dd k + 1 1 H k + 1
(115)
η k + 1 = 1 , if d k + 1 d k > 0 and 0 otherwise
(116)
H γγ k + 1 = ∂τ ∂γ | k + 1 , H γd k + 1 = ∂τ ∂d | k + 1 ,
(117)
H k + 1 = ∂Y ∂γ | k + 1 , H dd k + 1 = ∂Y ∂d | k + 1
(118)

At the end of each load step, the gradient of the level set is computed. If it is below 1 everywhere the next load step is applied. If not, a non-local zone is placed and the solver flowchart with non-local zone is used.

Solver flowchart with non-local zone

  1. 1.

    initialization: r l 0 = r i

  2. 2.

    increase non-local zone: r l n + 1 = r l n +Δ r l

  3. 3.

    iterations initialization: k=0

  4. 4.

    linear solve: solve (119) to find Δ u,Δ T,Δ ϕ

  5. 5.

    load update: T k+1=T k+Δ T

  6. 6.

    update in local zone (111)-(118), and non-local zone (120)-(124)

  7. 7.

    if residual ≤ tol, go to 8, else go to 4

  8. 8.

    if domain not fully broken (d(r i )<1), go to 2, else go to 9

  9. 9.

    end

The extent of the non-local zone is imposed and one tries to find a continuous displacement and damage field satisfying the problem.

The linear symmetric problem to be solved at each iteration when the non-local zone is not empty is to find ∈uΔU,ΔϕA,ΔTR such that

r i r l n + 1 H γγ k γ ( Δu ) + H γd k d ϕ k Δϕ γ ( u * ) r d r + r l n + 1 r e H k γ ( Δu ) γ ( u * ) r d r τ c ΔT u * ( r i ) = τ c T k u * ( r i ) r i r i r e τ k γ ( u * ) r d r , u * U r i r l n + 1 H k d ϕ k γ ( Δu ) + H dd k d 2 ϕ k + Y c Y k d ϕ k Δϕ ϕ * r d r = r i r l n + 1 Y k Y c d ϕ k ϕ * r d r , ϕ * A Δϕ + η k H dd k 1 H k γ ( Δu ) d ϕ k | r l n + 1 , + = ϕ k ϕ d k | r l n + 1 , +
(119)

where ηk is evaluated following (116). The update in the local zone follows (111)-(118) whereas in the non-local zone we have

u k + 1 = u k + Δu , ϕ k + 1 = ϕ k + Δϕ ,
(120)
d k + 1 = d ϕ k + 1
(121)
τ k + 1 = τ γ u k + 1 , d k + 1 , Y k + 1 = Y γ u k + 1 , d k + 1
(122)
H γγ k + 1 = ∂τ ∂γ | k + 1 , H γd k + 1 = ∂τ ∂d | k + 1 ,
(123)
H k + 1 = ∂Y ∂γ | k + 1 , H dd k + 1 = ∂Y ∂d | k + 1
(124)

Is is interesting to note the difference between the two solver flowcharts. When the non-local zone is empty, the linear solve deals only with displacement increments and the local update deals with the damage variable. On the contrary, when the non-local zone is not empty, the linear solve involves both displacement and damage (or more precisely the surrogate ϕ variable) increment in the non-local zone (local zone being treated as before).

The mesh is built so that it is much finer in the localization zone. Node j is located at a position x(j) given by

x(j)= ( r i + ( ( j 1 ) / N ) ( r e r i ) ) 2 r e r i + r i ,j=1,,N+1
(125)

where N is the number of elements considered. Results will be shown for the following mechanical parameters:

r i =0.1m, r e =0.2m, l c =0.02m, τ c μ =1 0 4 , d c =0.5
(126)

and numerical parameters

N=200, d inc =0.02
(127)

Regarding parameter Δ r l , the non-local zone is advanced by one element at a time or smaller when damage gets close to 1 at r i . This is done in order to capture the full load-displacement curve. The formula used in the simulation is

Δ r l =max min h * , l c ϕ n ( r i ) / 2 , 1 . e 12 * l c
(128)

where h* is the size of the element adjacent to the non-local zone at time step n. The initial (n=0) non-local zone needs to be more than one-element for convergence. Between 5 and 10 elements are used.

As a final remark on the solver flowchart with the non-local zone, we noticed that in non-local zone update, it was more efficient (reduced number of iterations) to take Δ ϕ as the one ensuring damage continuity rather that picking the one coming from step 3.

Discussion

In Figure 6, the force-displacement curve in the case n = 2 is shown. The figure indicates the different regime of the solution (pure elastic, local damage, coupled and pure non-local damage). Note that snap-back is taken into account automatically since the loading is not imposed but an unknown in the numerical scheme. Profiles of ϕ and damage along the radius at different loads are depicted in Figure 7.

Figure 7
figure 7

Distributions of ϕ / l c (top figure) and damage (bottom figure) for different loadings. Details on the curves in each plot from bottom to top: Bottom curve corresponds to T=1.232, damage evolution is purely local. The next curve is for T=1.351. It corresponds to the load at which ║ϕ║=1 at r=r i and non-locality steps in. Next curve is for T=1.368, it is the limit load. The load then decreases and damage evolution is purely non-local. Last two curves are for T=1.012 and T=0.02, respectively. The latter case depicts the profile at complete decohesion of the fiber.

Figure 8 shows the influence of the delocalization parameter n. Plots confirm the analytical limit results (97). As long as damage is purely local, all curves are superposed. As non-locality steps in, the delocalization parameters n plays a role.

Figure 8
figure 8

Force-displacement curves for n = 1 (small dots), n = 2 (solid line), n = 3 (big dots).

Finally, in order to show the insensitivity of the model with respect to the discretization parameters N, dinc, we show Figure 9 the influence of the choice of the N parameter (for the case n=2 and dinc=0.02). In Figure 10, we show the influence of dinc (for the case n=2 and N=50). Note that as expected, parameter dinc has only an influence when damage is purely local (rising part of the curve). For both figures, a zoom was used. Otherwise, curves cannot be distinguished.

Figure 9
figure 9

Force-displacement curves (zoom) for n = 2 and d inc =0 . 02 with different mesh sizes: N = 200 (solid line), N = 20 stars, N = 30 circles, N = 40 plus sign.

Figure 10
figure 10

Force-displacement curves (zoom) for n = 2 and N = 50 with different values of d inc :d inc =0 . 02 (solid line), d inc =0 . 1 stars, d inc =0 . 05 circles, d inc =0 . 01 plus sign.

Conclusions

The Thick Level Set damage model allows coupling local damage evolution in some part of the domain to a non-local damage evolution in the localization zone. Damage gradient is bounded. The bound is reached in the non-local zone (localization zone) and not reached in the local one. The localization zone boundary is the main unknown in the model. It evolves ensuring damage continuity. A semi-analytical 1D solution has been developed showing different regimes in the solution (elastic, local damage, coupled local and non-local damage and finally pure non-local damage). The solution was plotted using a numerical scheme. This numerical scheme is ad hoc for the 1D problem considered. The corresponding numerical implementation for 2D and 3D cases will be the subject of a forthcoming publication.

Authors’ contributions

All authors contributed to the main ideas in the paper: the way to couple local and non-local evolutions of damage. NM came up with the analytical solution. NM and NC designed the 1D code to plot the results. All authors read and approved the final manuscript.

References

  1. Moës N, Stolz C, Bernard P-E, Chevaugeon N: A level set based model for damage growth: the thick level set approach. Int J Numer Meth Eng 2011, 86: 358–380. 10.1002/nme.3069

    Article  Google Scholar 

  2. Bernard P-E, Moës N, Chevaugeon N: Damage growth modeling using the Thick Level Set (TLS) approach: efficient discretization for quasi-static loadings. Comput Meth Appl Mech Eng 2012, 233–236: 11–27. 10.1016/j.cma.2012.02.020

    Article  Google Scholar 

  3. Stolz C, Moës N: A new model of damage: a moving thick layer approach. Int J Fract 2012, 174: 49–60. 10.1007/s10704-012-9693-3

    Article  Google Scholar 

  4. Stolz C, Moës N: On the rate boundary value problem for damage modelization by Thick Level Set. ACOME 2012 Proceeding Ho-Chi-Minh, Viet Nam; 2012. [http://hal.archives-ouvertes.fr/hal-00725635] http://hal.archives-ouvertes.fr/hal-00725635 http://hal.archives-ouvertes.fr/hal-00725635

    Google Scholar 

  5. Karma A, Kessler D, Levine H: Phase-field model of mode III dynamic fracture. Phys Rev Lett 2001,87(4):045501. 10.1103/PhysRevLett.87.045501

    Article  Google Scholar 

  6. Miehe C, Welschinger F, Hofacker M: Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations. Int J Numer Meth Eng 2010,83(10):1273–1311. 10.1002/nme.2861

    Article  MathSciNet  Google Scholar 

  7. Spatschek R, Brener E, Karma A: Phase field modeling of crack propagation. Phil Mag 2011,91(1):75–95. 10.1080/14786431003773015

    Article  Google Scholar 

  8. Francfort GA, Marigo J-J: Revisiting brittle fracture as an energy minimization problem. J Mech Phys Solid 1998, 46: 1319–1412. 10.1016/S0022-5096(98)00034-9

    Article  MathSciNet  Google Scholar 

  9. Bourdin B, Francfort GA, Marigo J-J (2008) The Variational Approach to Fracture, 5–148, , [http://link.springer.com/10.1007/s10659-007-9107-3] Bourdin B, Francfort GA, Marigo J-J (2008) The Variational Approach to Fracture, 5–148,

  10. Comi C, Mariani S, Perego U: An extended FE strategy for transition from continuum damage to mode I cohesive crack propagation. Int. J. Numer. Anal. Meth. Geomech 2007, 31: 213–238. 10.1002/nag.537

    Article  Google Scholar 

  11. Sethian JA: Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision and material science. Cambridge University Press, UK; 1999.

    Google Scholar 

  12. Maugin GA: Internal variables and dissipative structures. J Non-Equilibrium Therm 1990, 15: 173–192.

    Google Scholar 

  13. Frémond M, Nedjar B: Damage, gradient of damage and principle of virtual power. Int J Solid Struct 1996,33(8):1083–1103. 10.1016/0020-7683(95)00074-7

    Article  Google Scholar 

  14. Comi C: Computational modelling of gradient-enhanced damage in quasi-brittle materials. Mech Cohesive-Frictional Mater 1999,36(April 1997):17–36. 10.1002/(SICI)1099-1484(199901)4:1<17::AID-CFM55>3.0.CO;2-6

    Article  Google Scholar 

  15. Pijaudier-Cabot G, Bazant ZP: Nonlocal dalmage theory. J Eng Mech ASCE 1987, 113: 1512–1533. 10.1061/(ASCE)0733-9399(1987)113:10(1512)

    Article  Google Scholar 

  16. Bazant ZP, Jirasek M: Nonlocal integral formulations of plasticity and damage: survey of progress. J Eng Mech 2002,128(November):1119–1149. 10.1061/(ASCE)0733-9399(2002)128:11(1119)

    Article  Google Scholar 

  17. Lubineau G, Azdoud Y, Han F, Rey C, Askari A: A morphing strategy to couple non-local to local continuum mechanics. J Mech Phys Solid 2012,60(6):1088–1102. 10.1016/j.jmps.2012.02.009

    Article  MathSciNet  Google Scholar 

  18. Azdoud Y, Han F, Lubineau G: A Morphing framework to couple non-local and local anisotropic continua. Int J Solid Struct 2013,50(9):1332–1341. 10.1016/j.ijsolstr.2013.01.016

    Article  Google Scholar 

  19. Lions P-L: Generalized solutions of Hamilton-Jacobi equations. Pitman Advanced Publishing Program, Boston; 1982.

    Google Scholar 

  20. Lorentz E, Godard V: Gradient damage models: toward full-scale computations. Comput Meth Appl Mech Eng 2011,200(21–22):1927–1944. 10.1016/j.cma.2010.06.025

    Article  MathSciNet  Google Scholar 

  21. Chung-Min L, Rubinstein J: Elliptic equations with diffusion coefficient vanishing at the boundary: theoretical and computational aspects. Quaterly Appl Math 2006, 64: 725–747.

    Google Scholar 

  22. Pradeilles-Duval RM, Stolz C: Mechanical transformations and discontinuities along a moving surface. J Mech Phys Solid 1995,43(1):91–121. 10.1016/0022-5096(94)00061-9

    Article  MathSciNet  Google Scholar 

  23. Hill R: Energy-momentum tensors in elastostatics: some reflections on the general theory. J. Mech. Phys. Solids 1986,34(3):305–317. 10.1016/0022-5096(86)90022-0

    Article  MathSciNet  Google Scholar 

  24. Stolz C: On micro-macro transition in non-linear mechanics. Materials 2010,3(1):296–317. 10.3390/ma3010296

    Article  Google Scholar 

  25. Rockafellar RT: Convex analysis. Princeton University Press, USA; 1970.

    Book  Google Scholar 

Download references

Acknowledgements

The support of the ERC Advanced Grant XLS no 291102 is greatfully acknowledged. Professor Antonio Huerta is also acknowledged for his advice.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Nicolas Moës.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ original submitted files for images

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (https://creativecommons.org/licenses/by/4.0), which permits use, duplication, 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 license, and indicate if changes were made.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Moës, N., Stolz, C. & Chevaugeon, N. Coupling local and non-local damage evolutions with the Thick Level Set model. Adv. Model. and Simul. in Eng. Sci. 1, 16 (2014). https://doi.org/10.1186/s40323-014-0016-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40323-014-0016-2

Keywords