Complete damage in elastic and viscoelastic media and its energetics

Weierstraÿ-Institut

für Angewandte Analysis und Sto hastik im Fors hungsverbund Berlin e.V.

Preprint

ISSN 0946  8633

Complete damage in elasti and vis oelasti media and its energeti s

Alexander Mielke1,2, Tomá² Roubí£ek3,4 , Jan Zeman5 1

Weierstraÿ-Institut für Angewandte Analysis und Sto hastik, Mohrenstraÿe 39, 10117 Berlin, Germany, E-Mail: mielkewias-berlin.de

2

Institut für Mathematik, Humboldt-Universität zu Berlin, Rudower Chaussee 25, 12489 Berlin-Adlershof, Germany

3

Mathemati al Institute, Charles University, Sokolovská 83, CZ18 675 Praha 8,

4

Institute of Information Theory and Automation, A ademy of S ien es,

Cze h Republi , E-Mail: roubi ekkarlin.m. uni. z

Pod vodárenskou v¥ºí 4, CZ18208 Praha 8, Cze h Republi

5

Department of Me hani s, Fa ulty of Civil Engineering, Cze h Te hni al University in Prague, Thákurova 7 CZ-166 29 Praha 6, Cze h Republi , E-Mail: zemanj ml.fsv. vut. z

submitted: 21st De ember 2007

No. 1285 Berlin 2007

Key words and phrases. Inelasti damage,small strain, energeti formulation. PACS. 46.15.-x, 62.20.-x .

The rst author a knowledges the support from C18-subproje t of the Resear h Center Matheon (Deuts he Fors hungsgemeins haft). The se ond author a knowledges partial support of this resear h from the grants IAA 1075402 (GA AV ƒR), LC 06052 and MSM 21620839 and 1M06031 (M’MT ƒR). The work of the third author was supported by the grant No. 106/08/1379 (GA ƒR) and by the resear h proje t MSM 6840770003 (M’MT ƒR). Support is a knowledged from the European grant MRTN-CT-2004-505226 Multi-s ale modelling and hara terisation for phase transformations in advan ed materials, too.

Edited by Weierstraÿ-Institut für Angewandte Analysis und Sto hastik (WIAS) Mohrenstraÿe 39 10117 Berlin Germany

Fax:

+ 49 30 2044975

E-Mail:

preprintwias-berlin.de

World Wide Web:

http://www.wias-berlin.de/

Abstra t A model for the evolution of damage that allows for omplete disintegration is addressed. Small strains and a linear response fun tion are assumed. The ow rule for the damage parameter is rate-independent. The stored energy involves the gradient of the damage variable, whi h determines an internal length-s ale.

Quasi-stati fully rate-independent evolution is onsidered as

well as rate-dependent evolution in luding vis ous/inertial ee ts. Illustrative 2-dimensional omputer simulations are presented, too.

1

Introdu tion

Damage, as a spe ial sort of inelasti response of solid materials, results from mi rostru tural hanges under me hani al load.

Routine omputational simulations

based on various models are widely performed in engineering, although mostly without being supported by rigorous mathemati al and numeri al analysis. We will onsider damage as a rate-independent pro ess.

This is often, although

not always, an appropriate on ept and has appli ations in a variety of industrially important materials, espe ially to on rete [Fre02, FrN95, Ort87℄, lled polymers [DPO94℄, or lled rubbers [GoS91, Mie95, MiK00℄. Being rate-independent, it is ne essarily an a tivated pro ess, i.e. to trigger a damage the me hani al stress must a hieve a ertain a tivation threshold. The mathemati al di ulty is ree ted in the fa t that only lo al-in-time existen e for a simplied s alar model or for a rate-dependent 0- or 1-dimensional model has been obtained, see [BoS04, DMT01, FKNS98, FKS99℄. The 3-dimensional situation was investigated in [FrG06, MiR06, MiRo℄ with the fo us to in omplete damage.

The main fo us of this paper is on

omplete damage, i.e. the material an ompletely disintegrate and its displa ement

ompletely loses any sense on su h regions. We show how mathemati al modeling

an be used to derive well-posed models by suppressing the use of the displa ement

u

and formulating everything in terms of stresses and energies. In Se tions 2-3 we

will negle t all rate dependent pro esses like vis osity and inertia so that the damage pro ess is quasistati and fully rate-independent. Eventually, in Se tions 4, we will

ombine a rate-independent damage pro ess with vis osity and inertia whi h are, of

ourse, rate-dependent. We onsider an anisotropi material but onne ourselves to a materials with linear

elasti response and an isotropi damage using only one s alar damage parameter under small strains (as in [BBT01, BoS04, Fre02, FrN96℄). Moreover, we use gradient-

of-damage theory [DBH96, Fre02, FrN95, FrN96, LoA99, Mau92, PMG04, StH03℄ 1

expressing a ertain nonlo ality in the sense that damage of a parti ular spot is to some extent inuen ed by its surrounding, leading to possible hardening or softeninglike ee ts, and introdu ing a ertain internal length s ale eventually preventing damage mi rostru ture development. From the mathemati al viewpoint, the damage gradient has a ompa tifying hara ter and opens possibilities for the su

essful analysis of the model. Anyhow, some investigations are still possible without gradient of damage, as shown in [FrG06℄ for in omplete damage, leading to a possible mi rostru ture in the damage prole. The goal of this arti le is to survey and further develop basi mathemati al tools fo used on omplete damage.

2

Complete quasistati damage at small strains

ϕ quadrati in terms of small-strain tensor e, parameter z , and onvex in terms of a gradient of

We will onsider spe i stored energy linear in terms of s alar damage

the damage

g: 1 κ ϕ(e, z, g) = zCe : e + |g|p + δ[0,+∞) (z), 2 p

(1)

C ∈ Rd×d×d×d is a positive-denite elasti ity tensor satisfying Cijkl = Cjikl = Cklij , d ∈ N denotes the onsidered spatial dimension, and κ > 0 is a so- alled fa tor of inuen e, and δ[0,+∞) is the indi ator fun tion of the interval [0, +∞), i.e. δ[0,+∞) (z) = 0 for z ≥ 0 while δ[0,+∞) (z) = +∞ otherwise. In this se tion, we

onsider the rst-order gradient of the damage prole ζ , hen e we put ∇ζ(t, x) in d pla e of the variable g ∈ R . Another ingredient of the damage-evolution model is

where

a spe i dissipated energy

̺(z) ˙ = where

a>0



−a z˙ +∞

if

z˙ ≤ 0,

d-dimensional by de reasing z

determines the phenomenology how mu h energy (per

volume) is dissipated by a

omplishing the damage pro ess, i.e. from 1 to 0.

(2)

elsewhere .

The value +∞ ree ts that we onsider damage as a unidire tional

pro ess, i.e. damage an only develop, but the material an never heal. Note that

̺ is

degree-1 homogeneous, whi h is related with the intended rate-independent evolution of the damage pro ess. Simultaneously, the level of the inelasti  driving for e

a is also a tivation threshold determining σi := ϕ′z (e, z, g) − divϕ′g (e, z, g) (with the

physi al dimension as stress) that triggers the damage pro ess.

As we want to fo us on omplete damage where the material an ompletely disintegrate, in the quasi-stati ase we annot have loading by dead load as e.g. gravity load. Thus we will onsider a hard-devi e loading by time-varying Diri hlet boundary onditions but zero bulk for es.

2

Considering the elasti  stress tensor

σe := ϕ′e (e, z, g)

that should be in quasistati equilibrium. Altogether formally we

onsider the problem


(∇u)⊤ +∇u div σe = 0, σe := ϕ′e (e(u), ζ, ∇ζ) = ζCe(u), e(u) = , 2  ∂ζ  + σi + σr ∋ 0, σr ∈ N[0,+∞) (ζ), ∂̺ ∂t σi := ϕ′z (e(u), ζ, ∇ζ) − divϕ′g (e(u), ζ, ∇ζ) 1 = Ce(u) : e(u) − div(κ|∇ζ|p−2∇ζ), 2 where

∂̺ denotes

the subdierential of the onvex fun tion

̺.

We also denoted

(3a)

(3b)

σr

a

0 ≤ ζ , and N[0,+∞) = ∂δ[0,+∞) denotes the normal

one. In fa t, as the evolution of ζ is unidire tional (here non-in reasing in time) and ζ will be pres ribed at the beginning, see (9) below, it always holds 0 ≤ ζ(t, x) ≤ ζ0 (x). Usually ζ0 = 1 is onsidered so even ζ ∈ [0, 1] a.e. on Q := (0, T ) × Ω. rea tion for e to the onstraint

This is indeed to be understood only formally be ause in the ompletely damaged part

ζ =0

and displa ements

u

as well as strain

e(u)

lose any sense.

Therefore, we will also onsider the regularized stored energy

κ 1 ϕε (e, z, g) = (z+ε)Ce : e + |g|p + δ[0,+∞) (z), 2 p

(4)

and then the regularized problem


div σe = 0,  ∂ζ  + σi + σr ∋ 0, ∂̺ ∂t

σe = (ζ+ε)Ce(u),

(5a)

σr ∈ N[0,+∞)(ζ),

1 σi = Ce(u) : e(u) − div(κ|∇ζ|p−2∇ζ). 2

As we have the displa ement well dened if

(5b)

ε > 0, we an easily onsider the Diri hlet

boundary onditions

u|Γ(t, x) = w(t, x) where

Γ ⊂ ∂Ω

is a part of the boundary of

Ω

(6)

where the hard-devi e loading is

applied. For simpli ity, the remaining boundary onditions are onsidered as homogeneous Neumann one:

Ce(u)ν = 0 where

ν

on

∂Ω\Γ

and

κ|∇ζ|p−2

∂ζ =0 ∂ν

denotes the outer unit normal to the boundary

∂Ω

on of

∂Ω,

Ω.

(7)

Then we dene

the Gibbs' stored energy

 κ 1   (ζ+ε)Ce(u) : e(u) + |∇ζ|p 2 p Gε (t, u, ζ) :=   +∞ 3

u|Γ = w(t, ·) and if 0 ≤ ζ a.e. on Ω, elsewhere , if

(8)

We still pres ribe an initial ondition

ζ0

for the damage prole:

ζ(0) = ζ0 . By the denition of the subdierential

̺(˜ z )},

(9)

∂̺(z) ˙ = {σ∈R; ∀˜ z ∈R : ̺(z) ˙ + (˜ z −z)σ ˙ ≤

the in lusion (5b) an equivalently be written as a variational inequality

 ∂ζ  
∂ζ  ∀˜ z: ̺ ≤ ̺(˜ z ) + z˜ − σi + σr ∂t ∂t

(10)

(t, x) ∈ Q := (0, T ) × Ω, where Ω ⊂ Rd is a onsidered domain o

upied by body and T > 0 a xed time horizon. This ould be used for a denition of a

for a.a. the

weak solution. Here, however, we an use homogeneity of

̺

to formulate a more suitable on ept

of so- alled energeti solution. By (2), we have

 ∂z  σi + σr ∈ ∂̺ ⊂ ∂̺(0) = [−a, +∞). ∂t

(11)

By the denition of the subdierential ∂̺(0) and properties of ̺, this means 0 = ̺(0) ≤ ̺(˜ z ) − (σi + σr )˜ z for any z˜ ∈ R. Written for z˜ − z instead of z˜, we have 0 ≤ ̺(˜ z −z)+(σi +σr )(˜ z −z). Further, by onvexity of ϕε (e, ·, ·), we have ϕε (e, z, g) ≤ ϕε (e, z˜, g˜) − ξ1 (˜ z −z) − ξ2 · (˜ g −g) for any (ξ1 , ξ2 ) ∈ ∂(z,g) ϕε (e, z, g). In parti ular, we 1 p−2 will use it for ξ1 = Ce : e + σr and ξ2 = (κ|∇ζ| ∇ζ). Altogether, substituting 2 e = e(u), g = ∇ζ(x) and z = ζ(x) we have Z Z


˜ ∇ζ˜ − (σi +σr ) ζ−ζ(t) ˜ ϕε e(u(t)), ζ(t), ∇ζ(t) dx ≤ ϕε e(u(t)), ζ, dx Ω Ω Z

˜ ∇ζ˜ + ̺ ζ−ζ(t) ˜ ≤ ϕε e(u(t)), ζ, dx

Ω

1,p ˜ ∀0 ≤ ζ∈W (Ω).

(12)

If ζ(t) satises (12), we say that ζ(t) is partially stable at t. Summing (5b) tested by ∂(u−w) ∂ζ ∂ζ with (5a) tested by , using −(σi +σr ) ≥ ̺( ∂ζ ) for any −(σi +σr ) ∈ ∂̺( ∂ζ ), ∂t ∂t ∂t ∂t ∂t integrating it over the time interval [0, T ], and applying by-part integration in time, we obtain formally the Gibbs-type energy balan e





Gε T, u(T ), ζ(T ) + Var̺ (ζ; 0, T ) ≤ Gε 0, u(0), ζ(0) + where of

ζ

w

Z TZ 0

 ∂w  dx dt σe : e ∂t Ω

w from (6) onto the whole Ω and the variation Var̺ ̺ (i.e. total dissipation of energy within the damage pro ess) is,

means an extension of

with respe t to

(13)

in view of (2), given by a simple formula

Var̺ (ζ; t1 , t2 ) =

 Z  a ζ(t1 , x)−ζ(t2, x) dx Ω



+∞

4

ζ(·, x) is nonde reasing on [t1 , t2 ] for a.a. x ∈ Ω, otherwise.

if

Let us denote by B and BV the Bana h spa e of everywhere dened bounded measurable and bounded-variation fun tions, respe tively. Moreover, let us abbre-

I := (0, T ), I¯ := [0, T ], Q := I × Ω, and Σ := I × Γ. ϕ(e, ·, ·) is onvex, (13) together with the partial stability

viate

It is important that, as (12) allows us to derive

ba kwards (10). This authorizes us to introdu e a denition of a solution:

Denition 2.1

(Weak/energeti solution.)

¯ W 1,2(Ω; Rd )) B(I;

and

We

all

¯ W 1,p (Ω; Rd )) ∩ BV(I; ¯ L1 (Ω)) ζε ∈ B(I;

(uε , ζε )

uε

with

∈

a weak/energeti so-

lution to the original problem (5) with the initial ondition (9) and the boundary

ondition (6)-(7) if the partial stability (12) holds for all

(i)

Z




ϕε e(uε (t)), ζε (t), ∇ζε (t) dx ≤

Ω

Z

Ω

i.e.



˜ ∇ζ˜ + ̺ ζ−ζ ˜ ε (t) dx ϕε e(uε (t)), ζ,

1,p ˜ ∀0 ≤ ζ∈W (Ω).

the energy inequality (13) holds with

(ii)

t ∈ [0, T ],

(uε , ζε )

in pla e of

(14)

(u, ζ),

(iii) (5a) is satised in the weak sense, i.e.

Z

∀t ∈ [0, T ], v ∈ W 1,2 (Ω; Rd ),

(ζε (t)+ε)Ce(uε(t)) : e(v) dx dt = 0 Q

v|Σ = 0, (iv) (6) and (9) hold with

(uε , ζε )

in pla e of

(15)

(u, ζ).

As the for e equilibrium (5a) is governed by minimization of the onvex fun tional

Gε (t, ·, ζε)

whi h also governs the evolution of

ζε ,

(5a) and the partial stability (12)

is equivalent to the full stability

Z

Ω

ϕε (e(uε (t)), ζε (t), ∇ζε (t)) dx ≤

Z

Ω

˜ ∇ζ) ˜ + ̺(ζ−ζ ˜ ε (t)) dx ϕε (e(˜ u), ζ, ˜ ∈ W 1,2 (Ω; Rd ) × W 1,p (Ω), ∀(˜ u, ζ) u˜|Γ = w(t), ζ˜ ≥ 0.

The onguration

(uε (t), ζε (t))

is alled stable at

t

if it satises (16). In fa t, under

(16), the energy inequality (13) implies even energy equality at any time







Gε t, uε (t), ζε (t) + Var̺ (ζε ; 0, t) = Gε 0, uε(0), ζ0 + with also

σe = (ζε + ε)Ce(uε ). ζ0 to be stable. The

Note that (16) at

(16)

Z tZ

t = 0

0

t,

i.e.

 ∂w  dx dt σe : e ∂t Ω

(17)

in fa t qualify through (9)

onditions (16)(17) leads to a on ept introdu ed in

[MiT99, MiT04, MTL02℄ (see also [Mie05℄ for a survey)

5

Denition 2.2

(Energeti

¯ W 1,2(Ω; Rd )) B(I;

and

solution.) We all (uε , ζε ) with uε ∈ ¯ W 1,p (Ω; Rd )) ∩ BV(I; ¯ L1 (Ω)) an energeti solution ζε ∈ B(I;

to the original problem (5) with the initial ondition (9) and the boundary ondition (6)-(7) if

t ∈ [0, T ], holds with (uε , ζε )

(i)

the stability (16) holds for all

(ii)

the energy balan e (17)

in pla e of

(u, ζ)

for all

and (iii) (6) and (9) hold with

(uε , ζε )

in pla e of

t ∈ [0, T ],

(u, ζ).

As already said, Denitions 2.1 and 2.2 are equivalent to ea h other. Under the 1,1 (I; W 1/2,2 (Γ)) (and thus onsidering an extension from hard-devi e loading w ∈ W W 1,1 (I; W 1,2(Ω)) for (13) or (17)), assuming p > d and stability of ζ0 , the existen e of a (weak) energeti solution

(uε , ζε )

proof onsists in limit passage with

is guaranteed for any

τ → 0

ε > 0,

f. [MiR06℄. The

for an approximate solution obtained

by the impli it time-dis retization with a time step

τ,

whi h leads to a re ursive

minimization problem

Minimize subje t to

for

ζτkε +ε Ce(∇ukτε ) : e(∇ukτε ) − aζτkε + 2 Ω 0 ≤ ζτkε ≤ ζτk−1 , ukτε |Γ = w(kτ ), ε ukτε ∈ W 1,2 (Ω; Rd ) , ζτkε ∈ W 1,p (Ω), Z

k = 1, ..., K := T /τ

with

ζτ0ε := ζ0 .

k = 1, ..., T /τ .

(18)

(ukτε , ζτkε ) to (18), we k that uτ ε |((k−1)τ,kτ ] = uτ ε

Having (some) solutions

assemble the pie ewise onstant interpolation for

 κ k p  ∇ζτ ε dx   p   

Likewise, we dene also

(uτ ε , ζτ ε ) so ζτ ε . Also, wτ

denotes the pie ewise

onstant interpolation of w . For the right-hand side of (19) below, we assume the ζτ ε (t) = ζτ0ε = ζ0 for t < 0, and similarly wτ (t) = w(0) and uτ ε (t) = u0τ ε 0 for t < 0, with uτ ε minimizing Gε (0, ·, ζ0 ). We have the two-sided dis rete energy

prolongation estimate:

Z tZ

 ∂w  dx dθ (ζτ ε +ε)Ce(uτ ε +w−wτ ) : e ∂θ 0 Ω
≤ Gε t, uτ ε (t), ζτ ε (t) + Var̺ (ζτ ε ; 0, t) − Gε (0, uτ ε (0), ζ0) Z tZ  ∂w  R ≤ (ζτRε +ε)Ce(uR +w−w ) : e dx dθ τε τ ∂θ 0 Ω

(19)

R holds with t = kτ for any k = 1, ..., T /τ , where (·)τ denotes fun tions retarded by R τ , i.e. [uτ ε ](t) := uτ ε (t − τ ), and where w has the meaning of an extension of the boundary onditions into Ω; f. [MiR06, Lemma 3.3℄. We are now going to formulate a suitable solution to the omplete damage problem. The essential pe uliarity is that displa ement

u

e(u) ζ = 0.

and the strain

well dened on areas that are ompletely damaged, i.e. where

are no longer

At ea h time t, we have, however, estimates on the stress (ζε (t)+ε)Ce(uε (t)) in d×d L2 (Ω; Rd×d sym ) uniform with respe t to ε > 0, where Rsym is the set of symmetri 6

d×d-matri es.

s is alled a realizable stress. The set of all L2 (Ω; Rd×d sym ), f. [BMR07, Proposition 2.8℄. R 1 A realizable stress s that minimizes s 7→ s : e(w(t)) dx is alled an ee tive 2 Ω stress at a given t. Let us remark that one an also dene an ee tive strain e ∈ L2loc ({x∈Ω; ζ(t, x) > 0}; Rd×d sym ) by Ea h weak luster point

realizable stresses is weakly ompa t in

e(t, x) = C−1

 s(t, x)  ζ(t, x)

t

for all

and a.a.

x∈Ω

su h that

ζ(t, x) > 0

(20)

d×d C−1 means the inversion of the mapping C : Rd×d sym → Rsym . It is important that e(t) is a orresponding limit of e(uε (t)) for ε → 0, f. [BMR07, Se t. 2.3℄ for details. Let us dene, for a given damage prole ζ , the ee tive stored energy as ˜ ε>0 : the so- alled Γ-limit [Dal95℄ of the olle tion {Gε (t, u, ζ)} where

g (t, ζ) :=

˜

min Gε (t, u, ζ). lim inf 1,p (Ω) ˜ 0≤ζ∈W ε→0+, ζ˜ ⇀ ζ in W 1,p (Ω) u∈W 1,2 (Ω;Rd )

(21)

The so- alled re overy sequen e of damage proles that asymptoti ally rea hes the ˜ = (ζ − δ)+ when δ → 0+ su iently slowly with respe t to value g (t, ζ) involves ζ

ε → 0+.

t and ζ there s(t, ζ) (i.e., div s = 0). Hen e, we an write Z κ p 1 s(t, ζ) : e(w(t)) + ∇ζ dx. g (t, ζ) = p Ω 2

An important result of [BMR07℄ is that for ea h

is unique

ee tive equilibrium stress

(22)

Also, we have an important formula for the power of external loading:

∂gg (t, ζ) = ∂t

Z

 ∂w  dx. s(t, ζ) : e ∂t Ω

(23)

Our denition for the omplete damage is based on the energeti -solution on ept as in Denition 2.2.

Denition 2.3

(Energeti

solution for omplete damage.)

ζ : data ϕ,

The pro ess

1,p

[0, T ] → W (Ω) is alled an energeti solution to the problem given by the ̺, w , and ζ0 , if, beside (9), also (i) ζ ∈ BV([0, T ]; L1 (Ω)) ∩ B([0, T ]; W 1,p(Ω)), (ii) it is stable for all t ∈ [0, T ] in the sense that Z

˜ ˜ g t, ζ(t) ≤ g (t, ζ) + ̺ ζ−ζ(t) dx ∀0 ≤ ζ˜ ∈ W 1,p (Ω),

(24)

Ω

(iii) and, for any

0 ≤ t1 < t2 ≤ T ,

the energy equality holds:




g t2 , ζ(t2) + Var̺ (ζ; t1, t2 ) = g (t1 , ζ(t1)) + in parti ular, the fun tion

Z t2Z t1

t 7→

R

Ω

 ∂w  dx dt, s(t, ζ(t)) : e ∂t Ω

(t)) dx s(t, ζ(t)) : e( ∂w ∂t 7

belongs to

(25)

L1 (0, T ).

Existen e of an energeti solution has been proved in [BMR07℄ by onvergen e of the above introdu ed regularization for

ε → 0.

Proposition 2.4

(Existen e of energeti solutions, onvergen e of (uε , ζε ).) p > d and w ∈ C 1 ([0, T ]; W 1/2,2 (Γ; Rd )), Then, there exist a subsequen e {εn }n∈N 1,p

onverging to 0 and a pro ess ζ : [0, T ] → W (Ω) being an energeti solution a

ording to Denition 2.3 su h that the following holds for all t ∈ [0, T ]: (i) Gεn (t, uεn (t), ζεn (t)) → g (t, ζ(t)), (ii) Var̺ (ζεn ; 0, t) → Var̺ (ζ; 0, t), 1,p (Ω), (iii) ζεn (t) → ζ(t) strongly in W 2 d×d (iv) (ζεn (t) + ε)C(e(uεn (t))) ⇀ s(t, ζ(t)) weakly in L (Ω; Rsym ). Let

Remark 2.5

(Quasi-stress.)

In fa t, we have bounded in

√

¯ L2 (Ω; Rd×d )) B(I; sym

not

(ζε +ε)Ce(uε ) but even ζε +εCe(uε ), whi h thus onverges (as a ∞ 2 d×d subsequen e) weakly* in L (I; L (Ω; Rsym )) to some χ. Let us all it quasi-stress. √ √ ζχ for the orresponding ee tive stress s and, by (20), χ = ζCe We have s = with the ee tive strain on the part with ζ > 0. Contrary to the stress itself whi h 2

onverges even L -strongly to zero on the ompletely damaged part, f. [BMR07, √ Proposition 2.5℄, ζε +εCe(uε ) need not onverge to zero on this part. only the stress

Remark 2.6

(Large strains.)

Generalization for stored energies that are non-

quadrati in terms of strain seems di ult, however. For in omplete damage (or, in other words,

ε>0

xed) we refer to [MiR06℄ where su h a model was analyzed

even at large strains and a unilateral onta t.

3

Numeri al implementation, 2D omputational simulations

In order to arrive at an implementable numeri al algorithm, we perform a spatial dis retization of the time-in remental minimization problem (18). To that end, 1,2 we introdu e nite-dimensional spa es Uh ⊂ W (Ω; R2 ) and Zh ⊂ W 1,p (Ω) and

onsider the following minimization problem:

Minimize subje t to

for

 ζτkhε +ε κ k p k k k  Ce(∇uτ hε ) : e(∇uτ hε ) − aζτ hε + ∇ζτ hε dx   2 p Ω k−1 k k 0 ≤ ζτ hε ≤ ζτ hε , uτ hε |Γ = w(kτ ),    k k uτ hε ∈ Uh , ζτ hε ∈ Zh Z

k = 1, ..., K := T /τ

with

(u0τ hε , ζτ0hε ) := (u0 , ζ0 ),

(26)

i.e. the dis retized in remental

problem leads to a non- onvex, box- onstrained optimization program. Note that the onvergen e of the fully dis rete solution to the solution of the spa e-time ontinuous problem is guaranteed thanks to abstra t approximation results available in [MiRo℄. 8

In the a tual numeri al implementation, the spatial dis retization is performed using the linear onforming nite elements, e.g. [BiS96, Bra07℄. Moreover, for omputational e ien y, we restri t our attention to

d=2

p=2

and dare to hoose

(whi h

ts with the theory presented in Se tion 2 only up to epsilon as we have required

p > d). For a given regularization parameter

ε

and the time level

k,

we express the dis rete

elds in the form

ζτkhε (x) = N ζh (x)ζ kh ,

ukτhε (x) = N uh (x)ukh , where

ukh

ζ kh

and

(27)

denote ve tors of the nodal values of displa ement and damage

parameter elds, respe tively (indi es τ ε are omitted in the sequel for the sake of ζ u brevity) and N h and N h denote the operators of pie ewise linear basis fun tions. The dis rete problem (26) an now be re-written in a fully algebrai format

 1 kT u k
k 1 kT ζ k ζT k  u K h ζ h uh + ζ h K h ζ h + f h ζ h 2 h 2  k , u 0 ≤ ζ kh ≤ ζ k−1 h,D = w D (kτ ) h

Minimize subje t to with omponents of

wD

(28)

storing the nodal displa ements on the Diri hlet part of the

boundary. The individual matri es are provided by:

K uh

(ζ h ) =

f ζh B

B uT h (x)

Z

ζ B ζT h κ(x)B h (x) dx,

Ωh

K ζh

where the

Z

=

Ωh

= −

Z

Ωh



ε+

N ζh (x)ζ h





C(x) B uh (x) dx,

(29)

(30)

a(x)N ζT h (x) dx,

(31)

operators ontain derivatives of the shape fun tions and

representation of the material stiness tensor

C;

C

is the Voigt

see e.g. [BiS96℄.

The dis rete formulation (28) leads to a (usually large-s ale) non- onvex program. k Nevertheless, re ognizing that the obje tive fun tion is quadrati separately in uh k and ζ h and exploiting the formal similarity between the ε-regularized damage model and the Fran fort-Marigo variational approa h to fra ture [BFM00℄, the problem (28)

an be e iently solved employing a variant of the alternate minimization algorithm proposed re ently by Bourdin in [Bou07, Bou℄. In the urrent ontext, the in remental version of algorithm is briey summarized in Table 1. In ea h internal iteration, the minimization problem with respe t to

u

(Step 4) redu es to the solution of

a sparse system of linear equations, while the subsequent sparse box- onstrained problem is solved using a ree tive Newton method introdu ed in [CoL96℄. The onvergen e of the alternate minimization was studied by Bourdin in [Bou07℄, where it was shown that the algorithm onverges to a riti al point of the dis retized problem in a nite number of iterations.

Of ourse, there is no guaran-

tee that the riti al point is a global minimizer of the non- onvex obje tive fun tion, whi h is a ru ial assumption of the theoreti al framework. 9

This obsta le

Table 1: Con eptual implementation of the optimization algorithm for time level (0) and an initial value of interval variable ζ . 1: Set 2:

k

j=0

repeat

j =j+1

3:

Set

4:

Solve

u(j) = arg

5:

Solve

ζ (j)

1 T u (j−1)
u u Kh ζ min uD =w D (kτ ) 2 1 T (j) u
(j) 1 T ζ T u K h ζ u + ζ K h ζ + f ζh ζ = arg min k−1 2 2 0≤ζ ≤ζ h

6:

until kζ (j) − ζ (j−1) k∞ ≤ δ

7: Set

ukh = u(j) , ζ kh = ζ (j)

an be, for example, resolved by resorting to the global sto hasti optimization approa hes [HJK00, IKLK04, MLZS00℄. Su h te hniques, however, require very large number of iterations and as su h are appli able only to very inexpensive obje tive fun tions.

Fortunately, it is possible to onstru t a feasible numeri al approa h

exploiting the two-sided energeti estimates (19). To that end, onsider the dis retized version of (19)

k Z X

jτ

 ∂w  dx dθ (ζτjhε +ε)Ce(ujτ hε+w−wτ ) : e ∂θ (j−1)τ Ω h j=1
0 ≤ Gε kτ, ukτhε , ζτkhε + Var̺ (ζτ hε ; 0, kτ ) − Gε (0, u0τ hε, ζhε ) Z Z k  ∂w  jτ X j−1 R ≤ η+ (ζτj−1 +ε)Ce(u +w−w ) : e dx dθ τ hε τ hε ∂θ (j−1)τ Ω h j=1 − η+

where

η

Z

(32)

is an energy toleran e parameter introdu ed for the numeri al implemen-

tation. The previous ondition is used to dete t lo al minimizers: if the result of k the alternate minimization algorithm ζ h fails to verify the inequality (32), the alk gorithm is restarted from the previous time level with ζ h used as an initial value k−1 for the minimization algorithm instead of ζ h . This pro edure is repeated until an admissible solution is found, see Table 2 for additional details. It is worth noting that the resulting algorithm shares similar features with the ba ktra king s heme introdu ed by Bourdin [Bou07℄ in the framework of variational fra ture theories. Performan e of the proposed algorithm will be illustrated on two ben hmark problems inspired by [SAS04℄: an inhomogeneous and a pre-not hed spe imen, see Figure 1. The orresponding geometri and material data together with the algorithm parameters are gathered in Figure 1 and Table 3, respe tively. Both stru tures are assumed to be in the plane stress state and are subje t to a proportional-in-time axially symmetri hard-devi e loading. In both ases, the spatial dis retization was performed using the unstru tured mesh generator T3D [Ryp98℄ and the problem size

[0, 1] was 100 identi al time steps (a physi al dimension of time is omitted in

was redu ed using symmetries of the spe imens. The analyzed time interval de omposed into

10

Table 2: Con eptual implementation of the time stepping pro edure. −1 0 (0) 1 : Set k = 1, ζ h = 0, ζ h = 0, ζ =0 2 : k 3 : Determine ζ h using the alternate minimization algorithm (0) . for time tk and initial value ζ (0) = ζ kh 4 : Set ζ 5 : two-sided inequality (32) is satised

repeat

if

6

:

7

:

8

:

9

:

10

:

Set

k =k+1

Set

k =k−1

else

end until k ≤ K

v(t) 0.2 m u(t)

Inhomogeneity Threshold a/2

y

0.2 m

0.5 m

u(t)

2m

1m

x

Notch

4m

(a)

(b)

v(t) 1m

Figure 1: S heme of simulated experiments; (a) inhomogeneous spe imen, (b) prenot hed spe imen

the sequel be ause of rate-independen e). Finally, for the inhomogeneous spe imen, the damage lo alization is triggered by pre-existing imperfe tions introdu ed by a redu ed a tivation threshold in the shaded area on the axis of symmetry. The resulting energeti s for the inhomogeneous spe imen is displayed in Figure 2 for a representative hoi e of the

ε

and

h

parameters. Clearly, in its basi version, the

dis rete solution obtained by the alternate minimization algorithm fails to provide an appropriate energeti solution to the problem. The two-sided inequality is satised only in the initial stage, where the spe imen stays mainly elasti . At time

t ≈ 0.61,

the damage propagates simultaneously through the spe imen, as manifested by the drop of the sum of the globally dissipated and the Gibbs energy, see Figure 2(a). Even after this instant, however, this quantity in reases, whi h is the onsequen e of the non-zero value of regularization parameter

ε.

Moreover, the damage prole

still evolves in the subsequent time levels, leading to the in rease in the dissipated energy balan ed by the ontribution of the Gibbs energy. With the ba ktra king option a tive, however, the algorithm dete ts the lo al minimizer at

t ≈ 0.61

and, following the dotted line in Figure 2(a), returns to the time

level where the in remental two-sided inequality is satised. After the ba ktra k-

11

300

upper bound

Energy [Jm−1 ]

Energy [Jm−1 ]

450 400 350 300 250 200 150 100 50 0 0

lower bound

Gε + Varρ

Gε Varρ

0.4

0.2

t

0.6

0.8

250 200 150

(a)

Gε + Varρ

lower bound

Gε

100 50 0 0

1

upper bound

Varρ 0.2

0.4

t

0.8

0.6

1

(b)

Figure 2:

0.03 m);

Global energeti s of the inhomogeneous spe imen (ε

= 5 · 10−2, h =

(a) Without ba ktra king (energy balan e fails), (b) with ba ktra king (an

approximate energeti solution)

ing stage is ompleted, the alternate minimization algorithm is apable of nding an approximate energeti solution, f. Fig. 2(b). As further illustrated by Fig. 3, evolution of the damage prole for the algorithm with ba ktra king is more gradual when ompared with the basi variant. Additional numeri al tests summarized in Figures 4 demonstrate the mesh-independent behavior of the method, i.e. the fa t that the global energeti response is almost independent of the dis retization parameter

h.

The inuen e of the energy regular-

ε → 0, the algoα rithms tries to reprodu e the one-dimensional optimal damage prole ζ(x, y) ≈ |x| ,

ization parameter

ε,

however, is mu h stronger, f. Figure 4(b). As

derived in [BMR07℄.

The same set of numeri al experiments was exe uted for the pre-not hed spe imen leading to the results appearing in Figures 5, 6 and 7.

When ompared to the

inhomogeneous spe imen, the global response shows similar trends for algorithms with and without ba ktra king. It is further onrmed by Figure 8 that the numeri al results are almost independent

Table 3: Parameter of the damage model and in remental algorithm Young's modulus, Possion's ratio,

E

ν

Fa tor of inuen e,

κ

A tivation threshold,

a

(see [FrN96℄)

Maximal pres ribed displa ement for the inhomogeneous spe imen Maximal pres ribed displa ement for the pre-not hed spe imen Time step,

τ

Damage prole toleran e,

δ

Two-sided energy inequality toleran e,

η

12

27 GPa 0.2 10 Jm−2 500 Jm−3 5 · 10−4 m 2.25 · 10−4 0.01 10−6 10−3

m

Without ba ktra king

With ba ktra king 1

0.4 0.3

0.5

0.2

0.8

0.3

0.6

0.2

0.4 0.2

0.1

0.1 0

0.4

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

0

0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

t = 0.32 0.4

0.995

0.3 0.2

0.99

0.8

0.3

0.6

0.2

0.4 0.2

0.1

0.1 0

0.4

0.985 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

t = 0.42 0.4

0.8

0.4

0.8

0.3

0.6

0.3

0.6

0.2

0.4

0.2

0.4

0.1

0.2

0.1

0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

0

2

0.2 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

t = 0.64 0.4

0.8

0.4

0.8

0.3

0.6

0.3

0.6

0.2

0.4

0.2

0.4

0.1

0.2

0.1

0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

0

2

0.2 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

t = 0.80 Figure 3: Time evolution of ζ eld for the inhomogeneous spe imen h = 0.03 m, ε = 5 · 10−2 , displa ements are s aled by a fa tor 100 and only a quarter of the spe imen is shown.

of the spatial dis retization parameter

h,

whi h is onsidered to be an essential

requirement for any damage model in the engineering literature.

The extent of

damage zone depends on the value of the regularization parameter

ε

residual energy after the omplete damage). As

ε → 0,

(related to a

however, the width of the

lo alized damage zone, displayed in Figure 9, remains still nite and insensitive to spatial dis retization.

Remark 3.1

(Clapeyron prin iple.) Similarly to [KMR06℄, it an be observed that

the work of external load is approximately equally distributed to the dissipated energy

Varρ and the stored energy Gε after the damage initiation; the ee t known as

the Clapeyron prin iple for slowly loaded bodies with vis ous damping, f. [FoT03℄. The deviation from the ideal parameters

1 : 1 ratio depends mainly on the energy regularization

ε, see Figures 4 and 7, whi h makes a ertain portion of the stored energy

unavailable to the damage pro ess. In addition, due to the lo alized hara ter of damage, only a part of the work of the external load an ontribute to the dissipative pro esses (analogously to the beginning of the loading program where no damage o

urs).

13

250 200

125

Gε

100 50

Varρ

0 0

0.4

0.2

0.6

t

75 50 25 0 0

1

0.8

(a)

Gε

Gε + Varρ

100

Gε + Varρ

150

ε = 5 · 10−2 ε = 10−2 ε = 10−4

150

h = 0.02 m h = 0.03 m h = 0.05 m

Energy [Jm−1 ]

Energy [Jm−1 ]

300

Varρ 0.2

0.4

0.6

t

0.8

1

(b)

Figure 4: Convergen e of the approximate energeti solution for the inhomogeneous −2 spe imen; (a) h → 0 m, ε = 5 · 10 , (b) ε → 0, h = 0.02 m; mesh with h = 0.05 m

ontains

493

h = 0.02

m to

triangular elements,

1, 549

h = 0.03

m orresponds to

1, 193

elements and

elements.

175 upper bound

100

Energy [Jm−1 ]

Energy [Jm−1 ]

150

lower bound

125

100 75 50

Gε

Gε + Varρ

25 0 0

Varρ 0.2

0.4

t

0.6

(a)

50

lower bound

Gε + Varρ

Gε

25 0 0

1

0.8

75

upper bound

Varρ 0.2

0.4

t

0.6

0.8

1

(b)

Figure 5: Global energeti s for the pre-not hed spe imen (ε

= 10−2, h = 0.03

m);

(a) Without ba ktra king (energy balan e fails), (b) with ba ktra king (an approximate energeti solution)

4

Damage in vis oelasti media with inertia

Finally we in lude also some rate-dependent phenomena, in parti ular vis osity and inertia.

Combination with vis osity has been addressed in Maxwellian rheology

(even with plasti ity) in [FeS03℄ and in the Kelvin-Voigt rheology in [HSS01, PPS07, SHS06, CFKSV06℄. We will onsider linear vis osity in the Kelvin-Voigt rheology, i.e. the total stress

σ

is omposed from the elasti ontribution σe := ζCe(u) as before and now also the ∂u ) where C is a positive-denite elasti ity tensor vis ous ontribution σv := ζDe( ∂t as before and D is a positive-denite vis osity tensor satisfying Dijkl = Djikl = Dklij . Note that, like the elasti response, it is natural to assume that also the vis ous response depends on the damage

ζ

and vanishes in the ompletely damaged. This

substantially diers from previous studies [FeS03, HSS01, PPS07, SHS06℄ whi h on-

14

1

1

1

0.9

0.9

0.9

0.9

0.8

0.8

0.8

0.8

0.7

0.7

0.7

0.7

0.6

0.6

0.6

0.6

0.5

0.5

0.5

0.5

0.4

0.4

0.4

0.4

0.3

0.3

0.3

0.3

0.2

0.2

0.2

0.2

0.1

0.1

0.1

0.1

0.9

0.995 0.8

0.7

0.99

0.6

0.985

0.5

0.4

0.98 0.3

0.2

0.975

0.1

0.97 0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0

1

0

0.1

0.2

0.3

t = 0.34

0.4

0.5

0.6

0.7

0.8

0.9

0

1

0

0.1

0.2

0.3

t = 0.64

0.4

0.5

0.6

0.7

0.8

0.9

1

t = 0.80

Figure 6: Snapshots of the time evolution of the ζ eld for the pre-not hed spe imens; h = 0.03 m, ε = 10−2 , displa ements are s aled by a fa tor 100 and only a half of the spe imen is displayed.

150 h = 0.05 m h = 0.03 m h = 0.02 m

50 50

Gε

Gε + Varρ

25 0 0

Varρ 0.2

0.4

t

0.6

100 75 50

(a)

Gε

Gε + Varρ

25 0 0

1

0.8

ε = 5 · 10−2 ε = 10−2 ε = 10−4

125

Energy [Jm−1 ]

Energy [Jm−1 ]

100

Varρ 0.2

0.4

t

0.6

0.8

1

(b)

Figure 7: Convergen e study for the pre-not hed spe imen; (a)

h→0

m,

ε = 10−2 ,

(b) ε → 0, h = 0.02 m; mesh with h = 0.05 m ontains 377 triangular elements, h = 0.03 m orresponds to 1, 229 elements and h = 0.02 m to 1, 773 elements.

sidered vis osity un hanged even in damaged material. Like in [PPS07, SHS06℄, we also onsider inertia related to the mass density

ρ.

Naturally, ontrary to the vis o-

elasti response, the inertial ee ts are independent of damage be ause the mass is not destroyed by damaging inter-atomi links. Thus the rate-independent evolution of the damage is now oupled with rate-dependent evolution of the displa ement. Due to the inertial ee ts, we an now impose dead loading by a bulk for e

f.

For

simpli ity, we then do not onsider any hard-devi e loading, i.e. we impose only the boundary onditions (7) with


∂2u − div σ + σ = f, v e ∂t2  ∂ζ  + σi + σr ∋ 0, ∂̺ ∂t ρ

Γ = ∅.

Altogether, formally, we onsider

σv = ζDe

∂u
, ∂t

σe = ζCe(u),

(33a)

σr ∈ N[0,+∞) (ζ)

1 σi := Ce(u):e(u) − div(κ|∇k z|p−2 ∇ζ). 2

(33b)

Of ourse, now we must pres ribe also the initial ondition on the displa ement and

15

1

1

0.9

0.9

0.8

0.8

1

0.9

0.9

0.9

0.9

0.8

0.8

0.8

0.8

0.7

0.7

0.7

0.7

0.7

0.7

0.6

0.6

0.6

0.6

0.6

0.6

0.5

0.5

0.5

0.5

0.5

0.5

0.4

0.4

0.4

0.4

0.4

0.4

0.3

0.3

0.3

0.3

0.3

0.3

0.2

0.2

0.2

0.2

0.2

0.2

0.1

0.1

0.1

0.1

0.1

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

h = 0.05

0.8

0.9

0

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

h = 0.03

m

Figure 8: Examples of the

ζ

0.8

0.9

0

1

0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

h = 0.02

m

eld distribution for

0.8

0.9

1

m

t = 1 (ε = 5 · 10−2 , h → 0

m);

displa ements are s aled by a fa tor 100 and only a half of the spe imen is displayed. 1

1

0.9

0.9

1

0.9

0.9

0.9

0.9

0.8

0.8

0.8

0.8

0.8

0.8

0.7

0.7

0.7

0.7

0.7

0.7

0.6

0.6

0.6

0.6

0.6

0.6

0.5

0.5

0.5

0.5

0.5

0.5

0.4

0.4

0.4

0.4

0.4

0.4

0.3

0.3

0.3

0.3

0.3

0.3

0.2

0.2

0.2

0.2

0.2

0.2

0.1

0.1

0.1

0.1

0.1

0.1

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0

1

0

0.1

0.2

0.3

0.4

ε = 5 · 10−2 Figure 9:

0.5

0.6

0.7

0.8

0.9

1

0

0

0.1

0.2

ε = 10−2

Examples of the

ζ

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

ε = 10−4

eld distribution for

t = 1 (h = 0.02

m,

ε → 0);

displa ements are s aled by a fa tor 100 and only a half of the spe imen is displayed.

the velo ity, so altogether we have

u(0, ·) = u0 ∈ W 1,2 (Ω; Rd ),

ζ(0, ·) = ζ0 ∈ W 1,p (Ω). We assume naturally

∂u (0, ·) = u˙ 0 ∈ L2 (Ω; Rd ), ∂t (34)

0 ≤ ζ 0 ≤ 1.

Similarly as before, let us take

ε>0

and onsider the regularized problem:

  ∂u
∂2u − div (ζ+ε)De + (ζ+ε)Ce(u) = f, ∂t2 ∂t  ∂ζ  1 + Ce(u) : e(u) − div(κ|∇z|p−2 ∇ζ) + N[0,+∞) (ζ) ∋ 0. ∂̺ ∂t 2 ρ

Its weak solution, let us denote it by

(uε , ζε ),

(35a) (35b)

an be obtained by rather standard

methods. The for e equilibrium (35a) in the weak form looks as

Z

0

T

D 2  Z   ∂ uε E ∂uε
ρ 2 , v + (ζε +ε) De + Ce(uε ) : e(v) − f · v dx dt = 0 ∂t ∂t Ω 16

(36)

2 1,2 for all v ∈ L (I; W (Ω; Rd )) with h·, ·i standing for the duality between W 1,2 (Ω; Rd )∗ 1,2 d and W (Ω; R ). Like (14) and (13), we have now the partial stability

κ ζε (t)+ε Ce(uε (t)) : e(uε (t)) + |∇ζε (t)|p dx 2 p Ω Z ˜
κ ˜p ζ+ε ˜ ε (t) dx ∀0 ≤ ζ˜ ∈ W 1,p (Ω) Ce(uε (t)) : e(uε (t)) + |∇ζ| + ̺ ζ−ζ ≤ p Ω 2

Z

for any

t ∈ [0, T ] with a from (2), and an energy inequality Z 2 ρ ∂uε (T ) + Gε (T, uε (T ), ζε (T )) dx + Var̺ (ζε ; 0, T ) 2 ∂t Ω Z ∂uε
∂uε
:e dx dt + (ζε +ε)De ∂t ∂t Q Z Z ∂uε ̺ 2 |u˙ 0| + Gε (0, u0 , ζ0 ) + f · dx dt; ≤ ∂t Q Ω 2

(37)

(38)

here we used ζ0 = 1 from (34) and, for oming from (13) to (38), we relied on (36) ∂uε ∈ L2 (I; W 1,2 (Ω; Rd )). Note that e(uε (T )) is well dened be ause for all v := ∂t ∂uε ∈ L2 (I; W 1,2 (Ω; Rd )) just due to the regularization by ε > 0. ∂t Now, as no minimization of stored energy applies, we unfortunately do not have at R 1 σ : e(w) dx for the stored energy,

f. (22). To avoid our disposal the formula like 2 Ω e R 1 usage of e(u) on the fully damaged parts, the stored energy ζCe(u) : e(u) Ω 2 R 1 √ dx an −1 alternatively be written as χ :C χe dx where we have denoted χe := ζCe(u) Ω 2 e −1 d×d d×d and, as above, C means the inversion of the mapping C : Rsym → Rsym . As in Remark 2.5, let us all χe an elasti quasi-stress; its physi al dimension is again ∂u 3 Pa=J/m as a standard stress. Similarly, to avoid usage of e( ), we introdu e the ∂t √ ζDe( ∂u vis ous quasi-stress χv := ) . ∂t

Also, let us denote the orresponding quasi-stresses for (35), i.e.

χe,ε =

p ζε +ε Ce(uε )

and

χv,ε =

p

ζε +ε De

∂uε
. ∂t

(39)

Then, in terms of these quasi-stresses, (36) rewrites to

Z

0

T

 D 2 Z p
∂ uε E −1 −1 ζε +ε χv,ε :D e(v) + χe,ε :C e(v) − f · v dx dt = 0. ρ 2 ,v + ∂t Ω

(40)

Moreover, (37) and (38) an be written as

Z

Q

κ 1 χe,ε :C−1 χe,ε + |∇ζε |p dx dt ≤ 2 p
κ ˜p ˜ ε dx dt + |∇ζ| − ̺ ζ−ζ p 17

Z

Q

˜ 1 ζ+ε χe,ε :C−1 χe,ε 2 ζε +ε ∀0 ≤ ζ˜ ∈ W 1,p (Ω)

(41)

t ∈ I and ρ ∂uε κ 2 1 (T ) + χe,ε (T ):C−1 χe,ε (T ) + |∇ζε (T )|p + δ[0,+∞) (ζε (T )) dx 2 ∂t 2 p Z + Var̺ (ζε ; 0, T ) + χv,ε :D−1 χv,ε dx dt Q Z Z 1+ε κ ∂uε ̺ 2 p |u˙ 0 | + Ce(u0 ):e(u0 ) + |∇ζ0 | dx + f· dx dt. ≤ 2 p ∂t Q Ω 2

to be satised for all

Z

Ω

(42)

∂uε We derive a-priory estimates that are independent of ε > 0 by testing (35a) by . ∂t ∂ζε ≤ 0 and symmetry and positive deniteness of C to obtain It is essential to use ∂t

 1∂ ∂uε
1 ∂ζε (ζε +ε)Ce(uε ) : e(uε ) = (ζε +ε)Ce(uε ) : e + Ce(uε ) : e(uε ) 2 ∂t ∂t 2 ∂t ∂uε
. (43) ≤ (ζε +ε)Ce(uε ) : e ∂t Thus

d dt

̺ ∂uε 2 ζε +ε Ce(uε ) : e(uε ) dx + 2 Ω 2 ∂t Z Z ∂uε ∂uε ∂uε ) : e( ) dx ≤ f · dx. + (ζε +ε)De( ∂t ∂t ∂t Ω Ω Assuming

Z

f ∈ L1 (I; L2 (Ω; Rd )), by Gronwall's inequality we obtain

∂u

ε ≤ C,

∂t L∞ (I;L2 (Ω;Rd ))

p

ζε +ε Ce(uε ) ∞ 2 ≤ Ce , L (I;L (Ω;Rd×d sym ))

p ∂uε


≤ Cv .

ζε +ε De ∂t L2 (Q;Rd×d sym )



ζε ¯ 1 ≤ C, BV(I;L (Ω))∩L∞ (I;W 1,p (Ω)) C , Ce , and Cv . ≤ Ce . From this,

with some onstants

kχe,ε kL∞ (I;L2 (Ω;Rd×d sym ))

In other words,

for

0 < ε ≤ 1,

(44)

the bounds

(45a) (45b) (45 ) (45d)

kχv,ε kL2 (Q;Rd×d ≤ Cv , sym )

and

we also obtain

∂2u

ε

2 2 ∂t L (I;W 1,2 (Ω;Rd )∗ )+L1 (I;L2 (Ω;Rd )) Z p
ζε +ε χv,ε :D−1 e(v) + χe,ε :C−1 e(v) − f ·v dx dt = sup kvkY ≤1 Q Z ≤ sup 2 χv,ε :D−1 e(v) + χe,ε :C−1 e(v) − f ·v dx dt kvkY ≤1

Q

−1

≤ 2|D |Cv + 2|C−1|Ce + 2kf kL1 (I;L2 (Ω;Rd )) . where

kukY = kukL2 (I;W 1,2 (Ω;Rd )) + kukL∞ (I;L2 (Ω;Rd )) . 18

(46)

∂χe,ε is available, whi h brings ∂t troubles by dening values of χe,ε at parti ular times in the limit. In the spirit of

Unfortunately, it does not seem that any estimate for

Denitions 2.1 and 2.3 but balan ing Helmholtz stored energy (sin e the by-part integration in time of the outer loading is no longer ne essary and advantageous) and in view of the estimates (45), we an exploit the above relations (36), (39), (41), and (42) when putting

ε = 0

for a denition of a weak/energeti solution to the

omplete-damage problem in the following way:

Denition 4.1

(Weak/energeti solution.)

We all

(u, χe , χv , ζ, E)

with

u ∈ W 1,∞ (I; L2 (Ω; Rd )),

(47a)

χe ∈ L∞ (I; L2 (Ω; Rd×d sym )), 2

χv ∈ L

(Q; Rd×d sym ), 1 ¯

¯W ζ ∈ BV(I; L (Ω)) ∩ B(I; ¯ E ∈ BV(I)

(47b) (47 )

1,p

(Ω)),

(47d) (47e)

su h that


∂u
∈ L2loc {(t, x) ∈ Q; ζ(t, x) > 0}; Rd×d sym , ∂t 2 ∂ u ∈ L2 (I; W 1,2 (Ω; Rd )∗ ) + L1 (I; L2 (Ω; Rd )) 2 ∂t

e

(48a) (48b)

a weak/energeti solution to the problem (33) with the initial onditions (34) and

Γ = ∅, if  D 2 Z p
∂ u E −1 −1 ζ χv :D e(v) + χe :C e(v) − f · v dx dt = 0 ρ 2 ,v + ∂t Ω

the homogeneous Neumann boundary ondition, i.e. (7) with

Z

T 0

for all

Z

v ∈ L2 (I; W 1,2 (Ω; Rd )),

A

(49)

if the partial stability

Z κ 1 ζ˜ 1 −1 p χe :C χe + |∇ζ| dx dt ≤ χe :C−1 χe 2 p 2 ζ A κ ˜p + |∇ζ| + ̺(ζ˜ − ζ) dx dt ∀0 ≤ ζ˜ ∈ Lp (I; W 1,p(Ω)) p

(50)

and

χe =

p ζ Ce(u)

on whi h

and

χv :=

ζ(t, x) > 0,

 ∂u  p ζ De ∂t

on any open

A⊂Q (51)

and if the energy inequality holds, i.e.

E(T ) +

Z

Ω

Z 2 ρ ∂u χv :D−1 χv dx dt (T ) + δ[0,+∞)(ζ(T )) dx + Var̺ (ζ; 0, T ) + 2 ∂t Z Q Z ∂u ̺ |u˙ 0 |2 dx + f · dx dt. (52) ≤ E(0) + ∂t Q Ω 2 19

with

E(0) =

1 Ce(u0 ) Ω 2

: e(u0 ) + κp |∇ζ0 |p dx and, for all t1 ∈ I , Z t1 Z t1 Z 1 κ E(t) dt ≥ χe :C−1 χe + |∇ζ|p dx dt. p 0 0 Ω 2

R

Remark 4.2

(53)

Let us omment this denition espe ially at the point that we laim

mu h less information on the ompletely damaged part than we did in the quasistati evolution in Se tion 2, whi h is related with what we are able to prove.

onsequen e, we also annot prove full energy balan e as an equality.

As a

Anyhow,

the granted a-priory estimates (45) and (46) give ertain solid base for engineering

al ulations and Denition 4.1 then indi ates what information we an surely read for the limit when

ε

approa hes zero. In fa t, we have bounds also on some other ∂ (ζε +ε) ∂t (Ce(uε ):e(uε )) whi h equals to χe,ε :D−1 χv,ε whi h is 2 1 bounded due to (45b, ) in L (I; L (Ω)).

derived quantities, e.g.

Proposition 4.3 L2 (Ω; Rd ),

and

p > d and f ∈ L1 (I; L2 (Ω; Rd )), u0 ∈ W 1,2 (Ω; Rd ), u˙ 0 ∈ ζ0 ∈ W 1,p (Ω), 0 ≤ ζ0 ≤ 1. Then there exists a weak/energeti Let

solution in a

ord to Denition 4.1.

∗ χe,ε ⇀ χe in L∞ (I; L2 (Ω; Rd×d sym )) d×d and χv,ε ⇀ χv in L (Q; Rsym ). Though the obtained χe need not be well dened R 1 at parti ular time levels, the stored energy Eε : t 7→ χ (t):C−1 χe (t) dx itself is Ω 2 e well dened and measurable be ause its sum with the kineti energy has a bounded

Proof. By (45b, ), we an hoose a subsequen e su h that 2

variation whi h is seen from (44) and (45 ). By Helly's prin iple, we hoose a subsequen e so that also

Eε (t) → E(t)

The limit passage in (40) uses

for all

ζε → ζ

t ∈ [0, T ].

in

Lq (Q)

with any

1 ≤ q < +∞,

whi h

follows by a generalized Aubin-Lions' theorem [Rou05, Cor.7.9℄ from the estimate ζε ∗ ¯ L1 (Ω)), and also it uses χe,ε ⇀ L∞ (I; W 1,p(Ω)) ∩ BV(I; χe in L∞ (I; L2 (Ω; Rd×d sym )) 2 d×d and χv,ε ⇀ χv in L (Q; Rsym )). in

ε ) in L2 (K; Rd×d e( ∂u sym ) ∂t on any ompa t ylinder K of the form [0, t]×K0 on whi h ζ > 0. Here we use a very spe ial stru ture of the problem that K0 ⊂ Ω su h that ζ(t) > 0 on K0 implies that, for any δ > 0, there is ε0 su h that for any 0 < ε ≤ ε0 we have ζε (t) + ε ≥ δ 1,p ¯ be ause p > d. for all x ∈ K0 ; here we used that W (Ω) is embedded into C(Ω) Thus also ζε + ε ≥ δ for all (t, x) ∈ K = [0, t]×K0 be ause ζε (·, x) is nonin reasing. Then we an pass to the limit in (39) and over A in (51) by ylinders of the form K above.

The limit passage in (39) uses also the bounds of

e(uε )

and

The limit passage in the partial stability ondition (41) in the term

Z

Q

Z ˜ 1 ˜ 1 ζ+ε −1 χe,ε :C χe,ε dx dt = (ζ+ε)Ce(uε ) : e(uε ) dx dt 2 ζε +ε Q 2

is more di ult than in the usual full stability (16) in the rate-independent ase. We must do it simultaneously with the left-hand-side term

Z

Q

1 (ζε +ε)Ce(uε ) : e(uε ) dx. 2 20

0 ≤ ζ˜ ≤ ζ and, following [BMR07, Proposition 2.10℄, put ζ˜δ := (ζ˜ − δ)+ . ˜δ (t) ≤ ζε (t) if ε > 0 is small enough (depending Then, for any xed δ > 0, we have ζ 1,p ¯ ompa tly. Simultaneously (Ω) ⊂ C(Ω) on t, however); re all that p > d so that W ˜ in W 1,p (Ω) ¯ . Indeed, let us onsider an open ǫ-neighbourhood Oǫ (t) of ζ˜δ (t) → ζ(t) ˜ x) = 0}. Then, for δ > 0 small enough, ζ˜δ > 0 on ¯ ζ(t, a ompa t set N(t) := {x∈Ω; ˜ x) − δ ¯ Ω\Oǫ (t). For a.a. x ∈ Oǫ (t)\N(t), we have either ζ˜δ (x) = 0 or ζ˜δ (t, x) = ζ(t, ˜δ (t, x) = 0 or ∇ζ˜δ (t, x) = ζ(t, ˜ x). Hen e, for δ > 0 small enough, and also ∇ζ Z Z p ˜ ˜ ˜ p dx ∇ζδ (t) − ∇ζ(t) dx = ∇ζ˜δ (t) − ∇ζ(t) Ω O (t)\N (t) Z ǫ ˜ p dx. ∇ζ(t) ≤ (54) Let us take

Oǫ (t)\N (t)

˜ p ∈ L1 (Ω) is ǫ → 0 be ause |∇ζ(t)| RT R ˜ p dx dt → 0 ∇ζ˜δ (t)−∇ζ(t) absolutely ontinuous for a.a. t ∈ [0, T ]. Then also 0 Ω Yet, the last expression an be pushed to zero with

by the Lebesgue dominated- onvergen e theorem; the ommon integrable majorant p ˜ is t 7→ k∇ζ(t)k . Lp (Ω;Rd ) Then, by the partial stability for

Z

ζε ,

we have

̺(ζε − ζ˜δ ) dx dt

Q

ζε +ε ζ˜δ +ε
κ κ Ce(uε ) : e(uε ) + |∇ζε |p − |∇ζ˜δ |p dx dt − 2 2 p p Q Z  1 ζ˜δ +ε  κ κ = 1− χe,ε : C−1 χe,ε + |∇ζε |p − |∇ζ˜δ |p dx dt. ζε +ε p p Q 2 ≥

Z

(55)

(ζ˜δ +ε)/(ζε +ε) = ζ˜δ /ζ onverges strongly in any Lq (K), q < +∞, L∞ (K) on every ompa t ylinder K of the form [0, t] × K0 where

Now we use that and weakly* in

ζ > 0, Z

K

as already used above. Then, by the weak lower semi ontinuity, we obtain

̺(ζ − ζ˜δ ) dx dt ≥

Z

K

ζ˜δ  κ κ 1 1− χe : C−1 χe + |∇ζ|p − |∇ζ˜δ |p dx dt. 2 ζ p p

(56)

˜δ → ζ˜ weakly* in L∞ (Q) be ause we proved already and use ζ p 1,p strong onvergen e in L (I; W (Ω)) and bounds in L∞ (Q). When overing A Then we pass

δ→0

involved in (50) by ylinders of the form

K,

we obtain just (50).

Limit passage in (42) is then by weak lower-semi ontinuity. Here we use also that that

Eε (t) → E(t)

and the weak lower semi ontinuity, hen e we get also (53).

2

Referen es [BBT01℄ E. Benvenuti, G. Borino, A. Tralli: A thermodynami ally onsistent nonlo al formulation for damaging materials, Euro. J. Me h. A/Solids 535-553. 21

21 (2001),

[BiS96℄ Z. Bittnar and J. ’ejnoha, Numeri al methods in stru tural me hani s, ASCE Press and Thomas Telford, New York and London, 1996. [BoS04℄ E. Bonnetti, G. S himperna: Lo al existen e for Frémond model of damage in elasti materials. Cont. Me h. Thermodyn.

16 (2004), 319335.

[BMR07℄ G. Bou hitté, A. Mielke, T. Roubí£ek: A omplete-damage problem at small strains. Zeit. für angew. Math. Phys., in print. http://dx.doi.org/10.1007/s00033-007-7064-0. [Bou07℄ B. Bourdin: Numeri al implementation of the variational formulation of brittle fra ture, Interfa e Free Bound.

9 (2007), 411430.

[Bou℄ B. Bourdin: The variational formulation of brittle fra ture: numeri al implementation and extensions. In: IUTAM Symp. Dis retization Meth. for Evolving

Dis ontinuities. (T. Belyts hko, A. Combes ure, R. de Borst, eds.), Springer, 2007. (to appear) [BFM00℄ B. Bourdin, G.A. Fran fort, J.-J. Marigo: Numeri al experiments in revisited brittle fra ture. J. Me h. Phys. Solids

48 (2000), 797-826.

[Bra07℄ D. Braess: Finite elements: Theory, fast solvers, and appli ations in solid

me hani s, third ed., Cambridge University Press, Cambridge, 2007. [CFKSV06℄ M. Campoa, J.R. Fernándeza, K.L. Kuttler, M. Shillor, J.M. Viaño: Numeri al analysis and simulations of a dynami fri tionless onta t problem with damage. Comput. Meth. Appl. Me h. Eng.

196 (2006), 476488.

[CoL96℄ T.F. Coleman and Y. Li: A ree tive Newton method for minimizing a quadrati fun tion subje t to bounds on some of the variables. SIAM J. Optim.

6 (1996), 10401058.

[Dal95℄ G. Dal Maso: An Introdu tion to

Γ- onvergen e.

Birkhäuser, 1995.

[DBH96℄ R. De Borst, A. Benallal, O.M. Heeres: A gradient-enhan ed damage ap-

6 (1996) C6, 491-509.

proa h to fra ture. J. de Physique IV

[DMT01℄ A. DeSimone, J.-J. Marigo, L. Teresi: A damage me hanism approa h to stress softening and its appli ation to rubber. Eur. J. Me h. A/Solids

20

(2001), 873-892. [DPO94℄ E.A. De Souza Neto, D. Peri , D.R.J. Owen: A phenomenologi al threedimensional rate-independent ontinuum damage model for highly lled polymers: Formulation and omputational aspe ts. J. Me h. Phys. Solids

42 (1994),

1533-1550. [FeS03℄ J.R. Fernández, M. Sofonea: Variational and numeri al analysis of the Signorini's on ta t problem in vis oplasti ity with damage. J. Appl. Math.

2

(2003), 87114. [FrG06℄ G.A. Fran fort, A. Garroni: A variational view of partial brittle damage evolution. Ar h. Rat. Me h. Anal.

182 (2006), 125152.

[Fre02℄ M. Frémond: Non-smooth thermome hani s. Springer, Berlin, 2002. [FKNS98℄ M. Frémond, K.L. Kuttler, B. Nedjar, M. Shillor: One-dimensional models of damage. Adv. Math. S i. Appl.

8 (1998), 541-570.

22

[FKS99℄ M. Frémond, K.L. Kuttler, M. Shillor: Existen e and uniqueness of solutions for a dynami one-dimensional damage model. J. Math. Anal. Appl.

229

(1999), 271-294. [FrN95℄ M. Frémond, B. Nedjar: Damage in on rete: the unilateral phenomenon.

Nu l. Eng. Des.

156 (1995), 323335.

[FrN96℄ M. Frémond, B. Nedjar:

Damage, gradient of damage and prin iple of

virtual power. Int. J. Solid Stru t.

33 (1996), 10831103.

[FoT03℄ R. Fosdi k, L. Truskinovsky: About Clapeyron's theorem in linear elasti ity. J. Elast.

72 (2003), 145172.

[GoS91℄ S. Govindjee, J.C. Simo: A mi ro-me hani ally based ontinuum damage model for arbon bla k-lled rubbers in orporating Mullins' ee t. J. Me h.

Phys. Solids

39 (1991), 87112.

[HSS01℄ W. Han, M. Shillor, M. Sofonea: Variational and numeri al analysis of a quasistati vis oelasti problem with normal omplian e, fri tion and damage.

J. Comput. Appl. Math.

137 (2001), 377398.

[HJK00℄ J. Haslinger, D. Jedelský, T. Kozubek, J. Tvrdík: Geneti and random sear h methods in optimal shape design problems, J. Glob. Optim.

16 (2000),

109131. [IKLK04℄ A. Ibrahimbegovi¢, C. Knopf-Lenoir, A. Ku£erová, P. Villon: Optimal design and optimal ontrol of stru tures undergoing nite rotations and elasti deformations, Int. J. Numer. Methods Eng.

61 (2004), 24282460.

[KMR06℄ M. Ko£vara, A. Mielke, T. Roubí£ek: A rate-independent approa h to the delamination problem. Math. Me h. Solids

11 (2006), 423-447.

[LoA99℄ E. Lorentz, S. Andrieux: A variational formulation for nonlo al damage models. Int. J. of Plasti ity

15 (1999), 119-138

[MLZS00℄ K. Matou², M. Lep², J. Zeman, M. ’ejnoha:

Applying geneti algo-

rithms to sele ted topi s ommonly en ountered in engineering pra ti e, Com-

put. Meth. Appl. Me h. Eng.

190 (2000), 16291650.

[Mau92℄ G.A. Maugin: The thermodynami s of plasti ity and fra ture. Cambridge Univ. Press, Cambridge, 1992. [Mie95℄ C. Miehe: Dis ontinuous and ontinuous damage evolution in Ogden-type large-strain elasti materials. Eur. J. Me h., A [MiK00℄ C. Miehe, J. Ke k:

14 (1995), 697-720.

Superimposed nite elasti -vis oelasti -plastoelasti

stress response with damage in lled rubbery polymers. Experiments, modelling and algorithmi implementation. J. Me h. Phys. Solids

48 (2000), 323-365.

[Mie05℄ A. Mielke: Evolution of rate-independent systems. In: Handbook of Dif-

ferential Equations: Evolutionary Di. Eqs. (Eds. C.Dafermos, E.Feireisl), pp. 461559, Elsevier, Amsterdam, 2005. [MiR06℄ A. Mielke, T. Roubí£ek: Rate-independent damage pro esses in nonlinear elasti ity. Math. Models and Meth. in Appl. S i. [MiRo℄ A. Mielke, T. Roubí£ek:

16 (2006), 177-209.

Numeri al approa hes to rate-independent pro-

esses and appli ations in inelasti ity. (Preprint No.1169, WIAS, Berlin, 2006.)

Math. Modelling Numer. Anal., submitted. 23

[MiT99℄ A. Mielke, F. Theil: A mathemati al model for rate-independent phase transformations with hysteresis. In: Models of ontinuum me hani s in analysis

and engineering. (Eds.: H.-D.Alber, R.Balean, R.Farwig), Shaker Ver., Aa hen, 1999, pp.117-129. [MiT04℄ Mielke, A., Theil, F.: On rate-independent hysteresis models. Nonlin. Di.

Eq. Appl.

11 (2004), 151189.

[MTL02℄ A. Mielke, F. Theil, V.I. Levitas:

A variational formulation of rate-

independent phase transformations using an extremum prin iple. Ar hive Rat.

Me h. Anal.

162 (2002), 137177.

[Ort87℄ M. Ortiz: An analyti al study of the lo alized failure modes of on rete.

Me h. Mater.

6 (1987), 159174.

[PPS07℄ C. Padovani, G. Pasquinelli, M. ’ilhavý, Pro esses in masonry bodies and the dynami al signi an e of ollapse. Math. Me h. Solids (2007), in print. [PMG04℄ R. H. J. Peerlings, T. J. Massart, M. G. D. Geers: A thermodynami ally motivated impli it gradient damage framework and its appli ation to bri k masonry ra king. Comput. Meth. Appl. Me h. Eng. [SAS04℄ A. Simone, H. Askes, L. J. Sluys:

193 (2004), 34033417.

In orre t initiation and propagation

of failure in non-lo al and gradient-enhan ed media, Int. J. Solid Stru t.

41

(2004), 351363. [Rou05℄ T. Roubí£ek:

Nonlinear partial dierential equations with appli ations.

Birkhäuser, Basel, 2005. [Ryp98℄ D. Rypl:

Sequential and parallel generation of unstru tured 3D meshes,

Ph.D. thesis, CTU in Prague, 1998, http://me h.fsv. vut. z/ dr/t3d.html. [SHS06℄ M. Sofonea, W. Han, M. Shillor: Analysis and Approximation of Conta t

Problems and Adhesion or Damage. Chapman & Hall, New York, 2006. [StH03℄ H. Stumpf, K. Ha kl: Mi rome hani al on ept for the analysis of damage evolution in thermo-vis oelasti and quasi-bittle materials. Int. J. Solids

Stru tures

40 (2003), 15671584.

24