- Email: [email protected]
Exact evaluation of dissipation for elastic-damage model dynamics
Computer methods in applied mechanics and engineering ELSEMER
Comput.
Methods
Appt. Mech. Engrg.
157 (1998)
1 I-18
Exact evaluation of dissipation for elastic-damage model dynamics L. Briseghella *, C. Majorana, P. Pavan Dipartimento di Costruzioni e Transporti, Facoltri di Ingegneria, Universitci degli Studi di Padova, Via F. Marzolo, 9-35131 Padova, Italy Received
25 April 1996
Abstract Some considerations about the finite strains dynamic of dissipative elasto-damage models are presented here. The dynamics of the system are treated taking care of the conservative character of the applied time-stepping algorithm in the sense of preserving the continuum balance equations. These are the linear and angular momentum laws and the expanded power of stress in respect of the Clausius-Duhem inequality. The generahsed mid-point rule for time integration and the use of a particular algorithmic form of the symmetric Piola-Kirchhoff stress tensor ensure the exact evaluation of the rate of the Helmoltz free-energy. This algorithm is known as energy-momentum method. The presented constitutive model is suited to describe polymers in finite deformations and with some accomodations appears as a way to investigate strain localization problems. Some simple examples, processed using an object-oriented f.e.m. software, show the effectiveness 0 1998 Elsevier Science S.A. of the algorithm in this last sense and the character of the constitutive model.
1. Introduction With reference to the works of many authors, we consider an elastic damage model suited to describe the behaviour of a large class of polymers. It is known that the mechanical qualities of the latter degrade themselves (presenting a softening behaviour) proportionally in some way to the maximum attained deformation. This is known also as Mullin’s effect. The damage phenomenon is evaluated using a scalar exponential function in terms of a parameter known as equivalent strain, which is a measure of the stored energy associated to the undamaged material. With this choice one obtains a very simple constitutive tensor. For the aims of the present paper we have adopted the hyperelastic constitutive model proposed by Simo and Pister [l]. Obviously, it is possible to employ different constitutive models to describe the behaviour of rubber like materials or hyperfoam materials. The typical viscous behaviour of polymers are not considered here but seem to be acommodable in the model. The use of damage criterions depending on the volumetric part of the deformation or on the deviatoric part of deformation is also possible. The aim of this paper is not to give a fear description of the polymers, but to undermark the conservative character of the algorithm and to show its effectiveness in the analysis of finite deformations dynamic problems. We also propose the present model as a possible way to describe the strain localisation phenomenon associated with softening behaviour. The adopted time integration was indicated by Simo and Tarnow [2] and is known as energy-momentum method. The method preserves, under particular imposed conditions, the exact conservation of linear momentum, angular momentum and the exact evaluation of the internal dissipation. Using a generalised mid-point rule one comes to the conclusion that the respect of the angular momentum law can be maintained only by the evaluation of external forces in the mid-point configurations. In addition, the exact estimation of internal dissipation can be obtained using specific
* Corresponding
author.
0045.7825/98/$19.00 0 1998 Elsevier PII SOO45-7825(97)00134-5
Science S.A. All rights reserved
12
L. Briseghella et al. I Comput. Methods Appl. Mech. Engrg. 157 (1998) 11-18
algorithmic expressions for the symmetric Piola-Kirchhoff stress tensor vs. strain relationship [Z]. Conservative configurations must be evaluated at each iteration and for each Gauss point. The simplicity of the model employed, however, causes no excessive additional computational costs if we compare these to the case of general hyperelastic constitutive models. The presented examples show that the conservative configurations to evaluate stresses in a conservative way do not present excessive variations. The application of the algorithm needs the derivability of the free energy function, a condition not always guaranteed in the case of dissipative constitutive models. The conservative configurations are evaluated by a systematic use of the Newton method. This method loses its applicability if the free energy function convexity vanishes. In these, rather uncommon, situations a dicotomic method is adopted. In the presented examples the damage function constants are chosen to have a marked softening behaviour. This is clearly unrealistic for the rubber-like materials but on the other hand enable us to show the potentiality of the code in the solution of softening behaviour problems.
2. The constitutive model Let B and S = q(B) be the reference and the current configurations of the body, respectively. The map p relates the material points X E B to the spatial points x = p(X) E S. The local condition J = det F > 0,being F the deformation gradient F = &rldX, ensures the regularity of the deformation. Considering the right CauchyGreen tensor C = FTF its principal invariants have the following expressions I, = tr[C]
I, = +
(tc[c-j -
I, = det[C] .
tlfc’])
(2.1)
The employed model is characterised by a monotonic function g(g s 1, g s 0) which can be used to estimate the damage degree of the material, The latter is described by a hyperelastic constitutive model if the equivalent strain is lower than its maximum value attained during time history. The free energy function is expressed by Simo [3] (2.2)
3
1c,= MC)
where w is the stored energy function
and can be expressed,
e.g. in the following
form [4]
w=,(13-1) ----(*+2p)lnI,+~p(l,-3).
(2.3)
4
Considering
the Clausius-Duhem
inequality
-$+$C>0,
(2.4)
one comes to the expressions
aw)
for the symmetric
ac
Dint =
For the parameter
g we employ
S= 2g
1
-w(c)g
Pioia-Kirchhoff
stress tensor and the internal
2 0.
a modified
dissipation
(2.5) form obtained
from the following
known expression
5 c 5,
I
with a smoothing in the neighbourhood of 5, or of .!Jy to ensure the derivability of the free energy function. This is the essential condition for the application of the mean value theorem (see Section 4.1). With this choice by hands an initial hyperelastic behaviour of the materials is possible. The parameters a, and b can take the following values
b E [O,l];
a E [O,m) .
L. Briseghella
et al.
I Comput. Methods Appl. Mech. Engrg.
The parameter b is dependent on the possible maximum attainment. The equivalent strain is defined as max =r$J;*, 5,
degree of damage and depends
13
on the rapidity
of its
(2.7)
f(lr”““, C(t)) =$&@$
a simple constitutive
tensor. The evolution
of damage
phenomenon
- C,““” S 0.
as f=O;
must
(2.8)
We can identify a damage surface by the associated loading from damage surface that is
equation.
The evolution
of damage can take place only for
$/X:C
1 g’=O constitutive
(2.9)
’
otherwise
The material
II-18
m.
This choice enables one to obtain preserve in any case the inequality
g’
1.57 (1998)
tensor C = 2 &fSlX has the following
a2w(c) 4X(5,“““)555+4
g’(5,“““) awcc> rnax 5,
ac@ac
expression
aww>as
f=O;
aflac:d>o
C=
(2.10) aZw(c) 4X( ST”“”m
otherwise
where the terms dw(C)/X, assume the form
aw
-=,(I,-
A
8’w(C)IX
aC can be easily evaluated
from the stored energy function
l)cY+$(z-c-1)
X
(2.11) and
(2.1 I)
(2.12) I being the second-order
unit tensor and (Zc.,),jLl= 1/2[(. )J. ),k + ( . ),k(. ),,I.
3. Continuum
equations
balances
We employ a Lagrangian description of the continuum. If B is the reference boundary JB the local equation of motion is written as [5] DN’)
configuration
+ p,,f = pori >
of the body with
(3.1)
where f are the body forces for unit of mass, P the first Piola-Kirchhoff stress tensor and ri the material velocity. The boundary conditions are imposed to have aJ? U 8,B = 3B and a$ fl d,B = 0. By the application of divergence theorem one obtains the following weak form of momentum balance equation
I
n
p,,I?:dB+
BP:GRAD(q)dB=
I
with the initial conditions variations of configurations. 3.1. Energy and momentum
IB
f.qdBt
t.?lddJ,
V= V, e q = q. at t = t, and being
balances
Q the functions
(3.2) of the space of admissible
equations
In the absence of external forces or in the presence of symmetries and with boundary conditions a,B = 0 the continuum is characterised by the conservation of linear and angular momentum, so we can write I;, = 0;
j,=o.
(3.3)
L. Briseghella et al. I Comput. Methods Appl. Mech. Engrg. 157 (1998) 11-18
14
In addition,
any kinetic energy variation
follows the relation
k = pint = P”“’ ,
(3.4)
where i is the time derivative of kinetic energy and Pin’, P’“’ the power of internal respectively. We also remember the validity of the Clausius-Duhem inequality
and external
forces,
4. Time discretization and discrete equations The intermediate ~V+p=PV~+, and the material
are written as
+(1-P>+%
PEKAll,
(4.1)
PE[Wl.
(4.2)
velocities
V7+fi=Pvn+, Considering
configurations
+(l
are -P)v,
these last equations
the weak form momentum
balance
The symmetry of the stress tensor enables us to establish the following variation of linear and angular momentum
becomes
conditions
for the respect of the laws of
L n+, =L,+Atj-;:a
Va~[O,ll,
(4.4)
J n+,=J,,+AtmF:U
a=1/2,
(4.5)
where f ext, mext are the resultant 4.1. Exact evaluation Using the convected
momentum
of external
forces.
of internal dissipation form of the right Cauchy-Green
c .+,=PC,+,+(l-P)C, and considering
force and the resultant
tensor
PE[O,11,
the mean value theorem
(4.6)
we obtain the following
relation
involving
the free energy function
(4.7) Therefore,
from Eq. (2.2) one comes to
(4.8) By the comparison -(*n+,
between
- #J+;s:!c,,,
we come to the conclusion symmetric Piola-Kirchhoff
this last one and the discrete form of the Clausius-Duhem -C,J=AD;;+,
inequality
,
that the last relation is satisfied only if we take the following stress tensor
(4.9) algorithmic
form for the
(4.10)
L. Briseghella et al. I Comput. Methods Appl. Mech. Engrg. 157 (1998) 11-18
In this way the algorithmic
internal
dissipation
is also specified
(4.11)
4.2. Second-order
accuracy
If one takes into account
the following
function (4.12)
which is odd about the value /? = l/2, one can also demonstrate symmetric Piola-Kirchhoff stress tensor
that the following
algorithmic
form of the
(4.13) has conservative qualities while in comparison with Eq. (4.10) ensures also a second-order accuracy. This can be viewed applying the Taylor’s series expansion to the discrete wake form of momentum balance equation. See Simo and Tamow [2] for more details. The conservative value of the parameter &, is obtained as a solution of the equation (4.14) and can be calculated in general, convex.
by the application
of a Newton method paying attention
to the fact that the function
is not,
5. Space discretization We use the Gale&in projections method for the discretization of the continuum. Let N’(X) i = 1,2, . . anode the shape functions with N’(Xi) = Sij. The admissible configurations can be expressed in terms of nodal values
&X)
= c ,=I
and the material
N’(X)u,,
velocities
(5.1)
result
“Mldl I+(X) = c N’(X)K ) I =I
(5.2)
where IC;, V, are the nodal displacements and the nodal material velocities. With the relation Vi,, =: (uk, , rcl)lAt at hand and by the application of Eq. (5.1) and Eq. (5.2) into expression (4.3) we obtain the discrete weak form
where the condition
(Y= II2 is imposed.
L. Briseghella et al. / Comput. Methods Appl. Mech. Engrg. 157 (1998) 11-18
16
The terms M” result from M” = J, N’N’p dV and are the coefficients of the mass matrix. It is demonstrated that the conditions about the conservative character of the preceding algorithmic forms of the stress tensor are extensible to the discrete space of Galerkin projections without any modifications.
6. Examples and applications The solution of the following numerical problems has showed a quadratic convergence of the algorithm. This clearly results from Table 1 in which the residual norm of the iteration of the same time step is presented. The conservative values of p are the zero of the function h(P) of Fig. 1 and are obtained by the application of a Newton scheme with /3 = 0.2 or p = 0.8 as initial values. When a Newton scheme is inapplicable a dicotomic method is used. The proposed constitutive model and the chosen constants do not describe a real material behaviour but are forced to capture some characteristic behaviours as strains localisation. For this aim the values for constants a, b of the damage law (2.6) are chosen to obtain a fast degradation of the mechanical elastic response with the amount of the equivalent strain. The damage law is shown in Fig. 2. 6.1. Elastic damage pendulum A pendulum, depicted in Fig. 3, is fixed in 0 and subjected to the gravity force. The damage-elastic model with Young modulus E = 1.8E4, Poisson ratio v = 0.4, a = 25 and b = 0.4 is adopted. For computational simplicity the pendulum is described with one only four nodes element with four Gauss points. Fig. 4 shows the deformation sequences for the first six periods of oscillation in the plane X - Z. The values of energies, depicted in Fig. 5, reveal that the total energy (as sum of gravity potential energy, internal and kinetic energy and dissipated energy) is constant. The energies diagram present an amount of dissipation and a progressive degradation of the elastic quality of the material when the pendulum mass centre passes through the vertical. After a number of cycles the pendulum has a pure elastic behaviour and the dissipation growth vanishes. In Fig. 6 is illustrated the typical relation between stress and strain.
Table 1 Residual norm for next iterations
of the same time step
Iteration no
Residual norm
Conservative
1 2 3 4 5
6.323E 1.155E 0.985E 0.423E 3.2318
7.879E 7.878E 7.878E 7.878E 7.878E
Fig. 1. Behaviour
+ + + -
006 004 001 002 006
of h(p) function. Fig. 2. Damage
-
The conservative
law for different
&
001 001 001 001 001
values of p are the intersections
values of constant
a with constant
with horizontal b = 0.
axe.
L. Briseghellaet al. / Comput.Methods Appt. hfech. E~rg.
Fig. 3. The elementary
discretization
Fig. 4. Deformed
of the pendulum.
shapes of the pendulum
The oscillations
__
4.00EiO6 2.00E+06+fi
/
__.__~_
Diss
Ekin n j&
are in the plane X - Z.
for the first six oscillations.
____ 6.00E+o6
17
157 (1998) II- 18
n
A-f?
f&f,
,
\
I\
,\
,\
n ,\
I\ ,\
$
,-._+
.-c1
1.2
1.4
1.6
-+-_m~ 2
1.6
+_---_I 2.2
2.4
PRINCIPAL STRETCH
Fig. 5. Energies
diagrams
Fig. 6. Stress-strain
6.2.Strain
localisation
for the pendulum
relation for the pendulum
of Example
6.1.
in the oscillations.
problem
A plane strain condition problem is simulated (using all the symmetries) with 42 eight nodes elements and clamping all the Z displacements as shown in Fig. 7. The damage-elastic constitutive model with Young modulus E = 2.1J3, Poisson ratio v = 0.4, a = 0.6 and b = 0.1 is here used. Nodal forces with the diagram of Fig. 7 are imposed to the upper bound of the plate. The deformed shapes of Fig. 8 show after 0.6 s the typical strain concentration bands due to the softening behaviour. The behaviour of an element of the shear bands and of a element closest to the shear bands are depicted in Figs. 9 and 10, respectively. The elements of the shear band show a strong softening behaviour while the other elements remain elastic.
Fig. 7. Discretization,
boundary
Fig. 8. Deformed
conditions
and applied forces for the plate of Example
shapes of the plate under plane strain condition.
6.2
18
L. Briseghella et al. I Comput. Methods Appl. Mech. Engrg. 157 (1998) II-18
Fig. 9. Stress-strain Fig. 10. Stress-strain
relation
for an element of the shear band,
relation for an element near the shear band
7. Conclusion We have extended the application of the so-called energy-momentum method proposed by Simo to the dynamic of a dissipative elasto-damage model. The adopted constitutive model is suited to describe, with some modifications, the behaviour of filled polymers. The time integration is thought to ensure the linear and total momentum lows respect and the exact evaluation of dissipation. This is obtained by calculating the symmetric Piola-Kirchhoff stress tensor in a particular intermediate configuration. The presented examples reveal the robustness of the algorithm and a quadratic convergence. The method seems to be a possible investigation way for dynamic problems of rubber-like materials in finite strains as well as (with obvious accommodations) for dynamic problems of brittle or quasi-brittle materials. About the latter aspect some problems related to softening behaviour are presented here.
References [I] J.C. Simo and K.S. Pister, Remarks on rate constitutive equations for finite deformation problem: Computational implications, Comput. Methods Appl. Mech. Engrg. 47 (1984) 201-215. [2] J.C. Simo and N. Tamow, The discrete energy-momentum method. Conserving algorithms for non-linear elastodynamics, 2. Angew. Math. Phys. 43 (1992) 757-792. [3] J.C. Simo, On a fully three-dimensional finite-strain visco-elastic damage model: Formulation and computational aspects, Comput. Methods Appl. Mech. Engrg. 60 (1987) 153-157. [4] J.C. Simo and T.J.R. Hughes, Elastoplasticity and viscoplasticity. Computational aspects, to appear. [5] J.E. Marsden and T.J.R. Hughes, Mathematical Foundations of Elasticity (Prentice-Hall, Englewood Cliffs, NJ, 1983).
Comput.
Methods
Appt. Mech. Engrg.
157 (1998)
1 I-18
Exact evaluation of dissipation for elastic-damage model dynamics L. Briseghella *, C. Majorana, P. Pavan Dipartimento di Costruzioni e Transporti, Facoltri di Ingegneria, Universitci degli Studi di Padova, Via F. Marzolo, 9-35131 Padova, Italy Received
25 April 1996
Abstract Some considerations about the finite strains dynamic of dissipative elasto-damage models are presented here. The dynamics of the system are treated taking care of the conservative character of the applied time-stepping algorithm in the sense of preserving the continuum balance equations. These are the linear and angular momentum laws and the expanded power of stress in respect of the Clausius-Duhem inequality. The generahsed mid-point rule for time integration and the use of a particular algorithmic form of the symmetric Piola-Kirchhoff stress tensor ensure the exact evaluation of the rate of the Helmoltz free-energy. This algorithm is known as energy-momentum method. The presented constitutive model is suited to describe polymers in finite deformations and with some accomodations appears as a way to investigate strain localization problems. Some simple examples, processed using an object-oriented f.e.m. software, show the effectiveness 0 1998 Elsevier Science S.A. of the algorithm in this last sense and the character of the constitutive model.
1. Introduction With reference to the works of many authors, we consider an elastic damage model suited to describe the behaviour of a large class of polymers. It is known that the mechanical qualities of the latter degrade themselves (presenting a softening behaviour) proportionally in some way to the maximum attained deformation. This is known also as Mullin’s effect. The damage phenomenon is evaluated using a scalar exponential function in terms of a parameter known as equivalent strain, which is a measure of the stored energy associated to the undamaged material. With this choice one obtains a very simple constitutive tensor. For the aims of the present paper we have adopted the hyperelastic constitutive model proposed by Simo and Pister [l]. Obviously, it is possible to employ different constitutive models to describe the behaviour of rubber like materials or hyperfoam materials. The typical viscous behaviour of polymers are not considered here but seem to be acommodable in the model. The use of damage criterions depending on the volumetric part of the deformation or on the deviatoric part of deformation is also possible. The aim of this paper is not to give a fear description of the polymers, but to undermark the conservative character of the algorithm and to show its effectiveness in the analysis of finite deformations dynamic problems. We also propose the present model as a possible way to describe the strain localisation phenomenon associated with softening behaviour. The adopted time integration was indicated by Simo and Tarnow [2] and is known as energy-momentum method. The method preserves, under particular imposed conditions, the exact conservation of linear momentum, angular momentum and the exact evaluation of the internal dissipation. Using a generalised mid-point rule one comes to the conclusion that the respect of the angular momentum law can be maintained only by the evaluation of external forces in the mid-point configurations. In addition, the exact estimation of internal dissipation can be obtained using specific
* Corresponding
author.
0045.7825/98/$19.00 0 1998 Elsevier PII SOO45-7825(97)00134-5
Science S.A. All rights reserved
12
L. Briseghella et al. I Comput. Methods Appl. Mech. Engrg. 157 (1998) 11-18
algorithmic expressions for the symmetric Piola-Kirchhoff stress tensor vs. strain relationship [Z]. Conservative configurations must be evaluated at each iteration and for each Gauss point. The simplicity of the model employed, however, causes no excessive additional computational costs if we compare these to the case of general hyperelastic constitutive models. The presented examples show that the conservative configurations to evaluate stresses in a conservative way do not present excessive variations. The application of the algorithm needs the derivability of the free energy function, a condition not always guaranteed in the case of dissipative constitutive models. The conservative configurations are evaluated by a systematic use of the Newton method. This method loses its applicability if the free energy function convexity vanishes. In these, rather uncommon, situations a dicotomic method is adopted. In the presented examples the damage function constants are chosen to have a marked softening behaviour. This is clearly unrealistic for the rubber-like materials but on the other hand enable us to show the potentiality of the code in the solution of softening behaviour problems.
2. The constitutive model Let B and S = q(B) be the reference and the current configurations of the body, respectively. The map p relates the material points X E B to the spatial points x = p(X) E S. The local condition J = det F > 0,being F the deformation gradient F = &rldX, ensures the regularity of the deformation. Considering the right CauchyGreen tensor C = FTF its principal invariants have the following expressions I, = tr[C]
I, = +
(tc[c-j -
I, = det[C] .
tlfc’])
(2.1)
The employed model is characterised by a monotonic function g(g s 1, g s 0) which can be used to estimate the damage degree of the material, The latter is described by a hyperelastic constitutive model if the equivalent strain is lower than its maximum value attained during time history. The free energy function is expressed by Simo [3] (2.2)
3
1c,= MC)
where w is the stored energy function
and can be expressed,
e.g. in the following
form [4]
w=,(13-1) ----(*+2p)lnI,+~p(l,-3).
(2.3)
4
Considering
the Clausius-Duhem
inequality
-$+$C>0,
(2.4)
one comes to the expressions
aw)
for the symmetric
ac
Dint =
For the parameter
g we employ
S= 2g
1
-w(c)g
Pioia-Kirchhoff
stress tensor and the internal
2 0.
a modified
dissipation
(2.5) form obtained
from the following
known expression
5 c 5,
I
with a smoothing in the neighbourhood of 5, or of .!Jy to ensure the derivability of the free energy function. This is the essential condition for the application of the mean value theorem (see Section 4.1). With this choice by hands an initial hyperelastic behaviour of the materials is possible. The parameters a, and b can take the following values
b E [O,l];
a E [O,m) .
L. Briseghella
et al.
I Comput. Methods Appl. Mech. Engrg.
The parameter b is dependent on the possible maximum attainment. The equivalent strain is defined as max =r$J;*, 5,
degree of damage and depends
13
on the rapidity
of its
(2.7)
f(lr”““, C(t)) =$&@$
a simple constitutive
tensor. The evolution
of damage
phenomenon
- C,““” S 0.
as f=O;
must
(2.8)
We can identify a damage surface by the associated loading from damage surface that is
equation.
The evolution
of damage can take place only for
$/X:C
1 g’=O constitutive
(2.9)
’
otherwise
The material
II-18
m.
This choice enables one to obtain preserve in any case the inequality
g’
1.57 (1998)
tensor C = 2 &fSlX has the following
a2w(c) 4X(5,“““)555+4
g’(5,“““) awcc> rnax 5,
ac@ac
expression
aww>as
f=O;
aflac:d>o
C=
(2.10) aZw(c) 4X( ST”“”m
otherwise
where the terms dw(C)/X, assume the form
aw
-=,(I,-
A
8’w(C)IX
aC can be easily evaluated
from the stored energy function
l)cY+$(z-c-1)
X
(2.11) and
(2.1 I)
(2.12) I being the second-order
unit tensor and (Zc.,),jLl= 1/2[(. )J. ),k + ( . ),k(. ),,I.
3. Continuum
equations
balances
We employ a Lagrangian description of the continuum. If B is the reference boundary JB the local equation of motion is written as [5] DN’)
configuration
+ p,,f = pori >
of the body with
(3.1)
where f are the body forces for unit of mass, P the first Piola-Kirchhoff stress tensor and ri the material velocity. The boundary conditions are imposed to have aJ? U 8,B = 3B and a$ fl d,B = 0. By the application of divergence theorem one obtains the following weak form of momentum balance equation
I
n
p,,I?:dB+
BP:GRAD(q)dB=
I
with the initial conditions variations of configurations. 3.1. Energy and momentum
IB
f.qdBt
t.?lddJ,
V= V, e q = q. at t = t, and being
balances
Q the functions
(3.2) of the space of admissible
equations
In the absence of external forces or in the presence of symmetries and with boundary conditions a,B = 0 the continuum is characterised by the conservation of linear and angular momentum, so we can write I;, = 0;
j,=o.
(3.3)
L. Briseghella et al. I Comput. Methods Appl. Mech. Engrg. 157 (1998) 11-18
14
In addition,
any kinetic energy variation
follows the relation
k = pint = P”“’ ,
(3.4)
where i is the time derivative of kinetic energy and Pin’, P’“’ the power of internal respectively. We also remember the validity of the Clausius-Duhem inequality
and external
forces,
4. Time discretization and discrete equations The intermediate ~V+p=PV~+, and the material
are written as
+(1-P>+%
PEKAll,
(4.1)
PE[Wl.
(4.2)
velocities
V7+fi=Pvn+, Considering
configurations
+(l
are -P)v,
these last equations
the weak form momentum
balance
The symmetry of the stress tensor enables us to establish the following variation of linear and angular momentum
becomes
conditions
for the respect of the laws of
L n+, =L,+Atj-;:a
Va~[O,ll,
(4.4)
J n+,=J,,+AtmF:U
a=1/2,
(4.5)
where f ext, mext are the resultant 4.1. Exact evaluation Using the convected
momentum
of external
forces.
of internal dissipation form of the right Cauchy-Green
c .+,=PC,+,+(l-P)C, and considering
force and the resultant
tensor
PE[O,11,
the mean value theorem
(4.6)
we obtain the following
relation
involving
the free energy function
(4.7) Therefore,
from Eq. (2.2) one comes to
(4.8) By the comparison -(*n+,
between
- #J+;s:!c,,,
we come to the conclusion symmetric Piola-Kirchhoff
this last one and the discrete form of the Clausius-Duhem -C,J=AD;;+,
inequality
,
that the last relation is satisfied only if we take the following stress tensor
(4.9) algorithmic
form for the
(4.10)
L. Briseghella et al. I Comput. Methods Appl. Mech. Engrg. 157 (1998) 11-18
In this way the algorithmic
internal
dissipation
is also specified
(4.11)
4.2. Second-order
accuracy
If one takes into account
the following
function (4.12)
which is odd about the value /? = l/2, one can also demonstrate symmetric Piola-Kirchhoff stress tensor
that the following
algorithmic
form of the
(4.13) has conservative qualities while in comparison with Eq. (4.10) ensures also a second-order accuracy. This can be viewed applying the Taylor’s series expansion to the discrete wake form of momentum balance equation. See Simo and Tamow [2] for more details. The conservative value of the parameter &, is obtained as a solution of the equation (4.14) and can be calculated in general, convex.
by the application
of a Newton method paying attention
to the fact that the function
is not,
5. Space discretization We use the Gale&in projections method for the discretization of the continuum. Let N’(X) i = 1,2, . . anode the shape functions with N’(Xi) = Sij. The admissible configurations can be expressed in terms of nodal values
&X)
= c ,=I
and the material
N’(X)u,,
velocities
(5.1)
result
“Mldl I+(X) = c N’(X)K ) I =I
(5.2)
where IC;, V, are the nodal displacements and the nodal material velocities. With the relation Vi,, =: (uk, , rcl)lAt at hand and by the application of Eq. (5.1) and Eq. (5.2) into expression (4.3) we obtain the discrete weak form
where the condition
(Y= II2 is imposed.
L. Briseghella et al. / Comput. Methods Appl. Mech. Engrg. 157 (1998) 11-18
16
The terms M” result from M” = J, N’N’p dV and are the coefficients of the mass matrix. It is demonstrated that the conditions about the conservative character of the preceding algorithmic forms of the stress tensor are extensible to the discrete space of Galerkin projections without any modifications.
6. Examples and applications The solution of the following numerical problems has showed a quadratic convergence of the algorithm. This clearly results from Table 1 in which the residual norm of the iteration of the same time step is presented. The conservative values of p are the zero of the function h(P) of Fig. 1 and are obtained by the application of a Newton scheme with /3 = 0.2 or p = 0.8 as initial values. When a Newton scheme is inapplicable a dicotomic method is used. The proposed constitutive model and the chosen constants do not describe a real material behaviour but are forced to capture some characteristic behaviours as strains localisation. For this aim the values for constants a, b of the damage law (2.6) are chosen to obtain a fast degradation of the mechanical elastic response with the amount of the equivalent strain. The damage law is shown in Fig. 2. 6.1. Elastic damage pendulum A pendulum, depicted in Fig. 3, is fixed in 0 and subjected to the gravity force. The damage-elastic model with Young modulus E = 1.8E4, Poisson ratio v = 0.4, a = 25 and b = 0.4 is adopted. For computational simplicity the pendulum is described with one only four nodes element with four Gauss points. Fig. 4 shows the deformation sequences for the first six periods of oscillation in the plane X - Z. The values of energies, depicted in Fig. 5, reveal that the total energy (as sum of gravity potential energy, internal and kinetic energy and dissipated energy) is constant. The energies diagram present an amount of dissipation and a progressive degradation of the elastic quality of the material when the pendulum mass centre passes through the vertical. After a number of cycles the pendulum has a pure elastic behaviour and the dissipation growth vanishes. In Fig. 6 is illustrated the typical relation between stress and strain.
Table 1 Residual norm for next iterations
of the same time step
Iteration no
Residual norm
Conservative
1 2 3 4 5
6.323E 1.155E 0.985E 0.423E 3.2318
7.879E 7.878E 7.878E 7.878E 7.878E
Fig. 1. Behaviour
+ + + -
006 004 001 002 006
of h(p) function. Fig. 2. Damage
-
The conservative
law for different
&
001 001 001 001 001
values of p are the intersections
values of constant
a with constant
with horizontal b = 0.
axe.
L. Briseghellaet al. / Comput.Methods Appt. hfech. E~rg.
Fig. 3. The elementary
discretization
Fig. 4. Deformed
of the pendulum.
shapes of the pendulum
The oscillations
__
4.00EiO6 2.00E+06+fi
/
__.__~_
Diss
Ekin n j&
are in the plane X - Z.
for the first six oscillations.
____ 6.00E+o6
17
157 (1998) II- 18
n
A-f?
f&f,
,
\
I\
,\
,\
n ,\
I\ ,\
$
,-._+
.-c1
1.2
1.4
1.6
-+-_m~ 2
1.6
+_---_I 2.2
2.4
PRINCIPAL STRETCH
Fig. 5. Energies
diagrams
Fig. 6. Stress-strain
6.2.Strain
localisation
for the pendulum
relation for the pendulum
of Example
6.1.
in the oscillations.
problem
A plane strain condition problem is simulated (using all the symmetries) with 42 eight nodes elements and clamping all the Z displacements as shown in Fig. 7. The damage-elastic constitutive model with Young modulus E = 2.1J3, Poisson ratio v = 0.4, a = 0.6 and b = 0.1 is here used. Nodal forces with the diagram of Fig. 7 are imposed to the upper bound of the plate. The deformed shapes of Fig. 8 show after 0.6 s the typical strain concentration bands due to the softening behaviour. The behaviour of an element of the shear bands and of a element closest to the shear bands are depicted in Figs. 9 and 10, respectively. The elements of the shear band show a strong softening behaviour while the other elements remain elastic.
Fig. 7. Discretization,
boundary
Fig. 8. Deformed
conditions
and applied forces for the plate of Example
shapes of the plate under plane strain condition.
6.2
18
L. Briseghella et al. I Comput. Methods Appl. Mech. Engrg. 157 (1998) II-18
Fig. 9. Stress-strain Fig. 10. Stress-strain
relation
for an element of the shear band,
relation for an element near the shear band
7. Conclusion We have extended the application of the so-called energy-momentum method proposed by Simo to the dynamic of a dissipative elasto-damage model. The adopted constitutive model is suited to describe, with some modifications, the behaviour of filled polymers. The time integration is thought to ensure the linear and total momentum lows respect and the exact evaluation of dissipation. This is obtained by calculating the symmetric Piola-Kirchhoff stress tensor in a particular intermediate configuration. The presented examples reveal the robustness of the algorithm and a quadratic convergence. The method seems to be a possible investigation way for dynamic problems of rubber-like materials in finite strains as well as (with obvious accommodations) for dynamic problems of brittle or quasi-brittle materials. About the latter aspect some problems related to softening behaviour are presented here.
References [I] J.C. Simo and K.S. Pister, Remarks on rate constitutive equations for finite deformation problem: Computational implications, Comput. Methods Appl. Mech. Engrg. 47 (1984) 201-215. [2] J.C. Simo and N. Tamow, The discrete energy-momentum method. Conserving algorithms for non-linear elastodynamics, 2. Angew. Math. Phys. 43 (1992) 757-792. [3] J.C. Simo, On a fully three-dimensional finite-strain visco-elastic damage model: Formulation and computational aspects, Comput. Methods Appl. Mech. Engrg. 60 (1987) 153-157. [4] J.C. Simo and T.J.R. Hughes, Elastoplasticity and viscoplasticity. Computational aspects, to appear. [5] J.E. Marsden and T.J.R. Hughes, Mathematical Foundations of Elasticity (Prentice-Hall, Englewood Cliffs, NJ, 1983).