Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

50, 3, pp. 841-853, Warsaw 2012

50th Anniversary of JTAM

STATIC GRADIENT DAMAGE SIMULATIONS USING

STABILIZED FINITE ELEMENTS

Jerzy Pamin, Adam Wosatko

Cracow University of Technology, Faculty of Civil Engineering, Kraków, Poland

e-mail: jpamin@L5.pk.edu.pl; awosatko@L5.pk.edu.pl

Thehourglass control of two-fieldfinite elements for the coupledproblemof gradientdamage
is analyzed. For the equilibrium equations, stabilization is introduced according to the least-
square method. For the additional averaging equation, three proposals of stabilization are
considered, however only the γ operatormethod performswell. Attention is focused on the
formulation, implementationandspectralanalysisof four-nodedelements in two-dimensional
simulations.Basicbenchmarks of a tensile barwith an imperfection andabeam in four-point
bending are computed and discussed.

Key words: gradient damage, stabilized finite elements, localization of deformation

1. Introduction

It is common knowledge that if standard elements with reduced integration (RI) are used and
improper hourglass modes have influence on the results, it can lead to a singularity of the as-
sembled tangent operator for the finite element model and a stable numerical analysis becomes
impossible. On the other hand, full integration (FI) can be inefficient, especially for large si-
mulation problems. Furthermore, FI can cause a locking phenomenon in the mesh. Therefore,
beside ideas like selective integration and B̄ formulation, the concept of computations using one
point integration with hourglass control is attractive for non-linear analysis.
The hourglass control, understood asmesh stabilization, was firstly applied in nonlinear FE

analysis for the displacement field u, see e.g. Belytschko et al. (1984). For two-field problems
like themixed u−p formulation, amathematically motivated stabilization was introduced only
for pressure p (Pastor et al., 1997) or for both fields (Commend et al., 2004). The hourglass
control is also implemented in coupled problems, for example thermo-mechanical ones (Reese,
2003).
The issue addressed in this paper is the simulation of continuumdamage.This involves, espe-

cially in quasi-brittle materials, the loss of well-possedness of the governing partial differential
equations, hence regularization is necessary, cf. de Borst et al. (1993).
One of the successful regularized models of quasi-brittle failure is the gradient-enhanced

damage formulation (Peerlings et al., 1996). To solve this coupled problem of equilibrium and
nonlocal averaging, two-field finite elements are used. Quadratic interpolation of the displace-
ments and linear of the averaged strain can be introduced, hence a different number of degrees of
freedom (dofs) at the corner and midside nodes is required. Although this interpolation option
seems to be optimal, other possibilities can give stable results as shown in Pamin et al. (2003),
since the analyzed problem is coupled rather thanmixed. It ismentioned in Simone et al. (2004)
that the so-called inf-sup condition does not have to be obeyed in this case. Oscillations which
may appear for some secondary fields, for example the stress field, have a local character, i.e.
occur only in zones of strong damage variation.Moreover, a linear interpolation of both fields is
tempting due to the possibility of using one sampling point, of coursewith the hourglass control.



842 J. Pamin, A.Wosatko

Especially in three-dimensional (3D) simulations, to reduce the computational cost andavoid
locking phenomena, it would be advantageous to use eight-noded elements with linear interpola-
tion of both the displacement vector field and the averaged strain field, and one-point Gaussian
integration. Before attacking the 3D problem, this article covers the results obtained for the 2D
four-noded gradient damage element with one-point reduced integration.

Variational andmatrix equations for the gradient damagemodel are recollected in Section 2
in order to derive stabilization terms in Section 3. The analysis of the rank of stiffnessmatrix for
the two-field four-nodedFE is presented in Section 4. In next sections, we describe the results of
two localized deformation tests, namely one-dimensional tension of a bar with an imperfection
and four point bending of a concrete beam. Some final remarks are listed in Section 7.

2. Variational and matrix equations for gradient damage

We begin our considerations from a recollection of the governing equations, weak forms and
space discretization in gradient damage (Peerlings et al., 1996). The whole set of small strain
equations is as follows

L
T
σ+b=0 ǫ=Lu σ=(1−ω)Eǫ (2.1)

where (2.1)1 are equilibriumequations, (2.1)2 are kinematic equations and (2.1)3 are constitutive
relations for the scalar damage theory. In the above equations L is the matrix of differential
operators, σ is the stress tensor in vector form, b is the body force vector, ǫ is the strain tensor
in vector form, u is thedisplacement vector, ω is the scalar damageparameterwhich grows from
0to 1and E is the elastic stiffness operator.Thepostulate of strain equivalence (Lemaitre, 1971)
is adopted, so that the effective stress tensor σ̂ = Eǫ as a fictitious undamaged counterpart
of stresses is distinguished.We assume adequate boundary conditions. Thewell-known problem
with mesh-dependent results is overcome by introducing the following differential equation for
an averaged strain measure ǭ

ǭ− c∇2ǭ= ǫ̃ (2.2)

It is assumed that the homogeneous natural boundary condition nT∇ǭ = 0 holds. A damage
loading function fd = ǭ− κd = 0 in strain space is used and standard loading-unloading
conditions are applied. The damage ω grows with the history parameter κd. The gradient
enhancement guarantees that the damage theory is nonlocal and the results become mesh-
independent. The parameter c > 0 has a unit of length squared and it is connected with an
internal length scale l of the material. The relation c= l2/2 is derived for instance in Askes et
al. (2000).

To solve the problem, apart from the approximation of displacements u, we additionally
discretize the averaged strain measure ǭ. The weak forms of equilibrium equations (2.1)1 and
averaging equation (2.2) are as follows

∫

B

δǫTσ dV =

∫

B

δuTb dV +

∫

∂B

δuTt dS

∫

B

δǭǭ dV +

∫

B

(∇δǭ)Tc∇ǭ dV =

∫

B

δǭǫ̃ dV
(2.3)

Here t denotes tractions. Next we can apply the spatial interpolation of displacements u=Na
and the averaged strain measure ǭ=hTe where N and h contain respective shape functions.
Further, the strain field is approximated by the relation ǫ=Ba and the gradient of the averaged



Static gradient damage simulations... 843

strain by ∇ǭ=qTe, where B and q contain derivatives of the respective shape functions. Now
we obtain two variational equations

R1(δa,a,e) = δa
T

∫

B

BTσ(a,e) dV − δaT
∫

B

NTb dV − δaT
∫

∂B

NTt dS =0

R2(a,δe,e) = δe
T

∫

B

(hhT+ cqqT)e dV − δeT
∫

B

hǫ̃(a) dV =0
(2.4)

After linearization of this gradient damage formulation, we introduce the following submatrices
and vectors as in Peerlings et al. (1996)

Kaa =

∫

B

B
T(1−ωi)EB dV Kae =−

∫

B

G
i
B
T
σ̂
i
h
T dV

Kea =−

∫

B

h[sT]iB dV Kee =

∫

B

(hhT+ cqqT) dV

f
i+1
ext =

∫

B

N
T
b
i+1 dV +

∫

∂B

N
T
t
i+1 dS fiint =

∫

B

B
T
σ
i dV

f
i
ǫ =

∫

B

hǫ̃i dV fie =Keee
i

(2.5)

where

G =
[ ∂ω
∂κd

][∂κd

∂ǭ

]
sT =

dǫ̃

dǫ

Based on iteration i, the increments of the primary fields da and de are computed for iteration
i+1
[
Kaa Kae
Kea Kee

][
da
de

]
=

[
f
i+1
ext − f

i
int

fiǫ − f
i
e

]
(2.6)

The tangent operator in the above matrix equation is non-symmetric. The gradient damage
model has been applied with success in the simulations of fracture in various materials, e.g.
composites (Geers, 1997) or concrete (Geers et al., 2000; Peerlings et al., 1996). Its version
refined by coupling to hardening plasticity (de Borst et al., 1999) enables the incorporation of
the physically observed irreversible strains.

3. Stabilization terms in variational equations

3.1. Stabilization of equilibrium equations

Themixed u−p formulation inCommendet al. (2004) is the startingpoint for thederivations
below. It is possible to apply one integration point and control hourglass modes in the solution
of equilibrium equations by the Galerkin least-square (GLS) method as in Zienkiewicz et al.
(2005). A stabilization term is added to Eq. (2.4)1

R1(δa,a,e)+R
stab
1 (δa,a,δe,e) = 0 (3.1)

The term Rstab1 can be defined according to theGLSmethod and for continuous fields is written
as follows

Rstab1 (δu,u,δǭ, ǭ)=
nel∑

e=1

∫

Be

[LTσ(δu,δǭ)]Tχ1[L
T
σ(u, ǭ)+b] dV (3.2)



844 J. Pamin, A.Wosatko

The stabilization scaling matrix χ1 is assumed as

χ1 =
χh2e
2G
I (3.3)

Here χ is an arbitrary, but possibly small value, he is a characteristic dimension of the finite
element, for example its diagonal, and G is the shear modulus. A further explanation of the
definition of χ1 and the analysis of units are given in Commend et al. (2004).

The weighting part of the stabilization term for the equilibrium equations can be written in
the following manner

P1(δu,δǭ)=L
T
σ(δu,δǭ)=LT[(1−ωi(δǭ))Eǫ(δu)] (3.4)

As was done with the plastic part in Commend et al. (2004), to avoid the linearization of the
damage part in Eq. (3.4), the damage contribution is omitted and only the elastic one kept.
Further, instead of P1(δu,δǭ), its discrete counterpart P1(δa) is introduced

P1(δa) =L
T
EBδa=Geaδa (3.5)

where an additional matrix is defined: Gea =L
TEB.

The equilibrium is obtained in an iterative process and the stress is decomposed as
σi+1 =σi+dσ. The definition of Rσ,i+1 equal to L

T
σi+1(u, ǭ) is introduced to obtain

Rσ,i+1 =L
T(σi+dσ)=Rσ,i+dRσ (3.6)

Constitutive equation (2.1)3 of the damage theory can be rewritten in rate form

σ̇=(1−ω)Eǫ̇− ω̇Eǫ (3.7)

Further definitions are introduced: Eaa = (1−ω
i)E and Eae =−G

iEǫi, so that the following
expression for dRσ is derived

dRσ =L
T(EaaBda+Eaeh

Tde) (3.8)

Finally, we obtain the residuum

Rσ,i+1 =Rσ,i+G
d
ada+G

d
ede (3.9)

In the above relation, we have: Gda =L
TEaaB and G

d
e =L

TEaeh
T.

According to these derivations, the stabilization term (3.2) is equal to

Rstab1 =
nel∑

e=1

∫

Be

(Geaδa)
T
χ1(Rσ,i+G

d
ada+G

d
ede+b) dV (3.10)

The following submatrices and vectors are defined

K̃aa =

∫

Be

(Gea)
T
χ1G

d
a dV K̃ae =

∫

Be

(Gea)
T
χ1G

d
e dV

f̃ =

∫

Be

(Gea)
T
χ1(Rσ,i+b) dV

(3.11)



Static gradient damage simulations... 845

3.2. Stabilization of averaging equation

In a similar manner, a stabilization term can be added to the averaging variational equa-
tion (2.4)2

R2(a,δe,e)+R
stab
2 (δa,a,δe,e) = 0 (3.12)

In fact, in Commend et al. (2004) the sign before the term R2 is changed to preserve the
positive definiteness of the tangent operator. Here the sign remains positive. Apparently simple
to perform, an analogical GLS method seems questionable for the averaging equation. After
discretization, it turns out that for rectangular elements ∇2hT = 0, hence this method does
not remove spurious singular modes. Clearly, in this case K̃ee is defined by the product hh

T

like the original operator Kee, and the hourglass control fails. It is shown in detail in Wosatko
(2008) and confirmed by a spectral analysis. Similarly, for the idea taken from Harari et al.
(2002), where the gradient GLS method is considered, the stabilization operator K̃ee in the
matrix equation arises as a result of multiplying matrix q and its transposition. An identical
term is included in Kee, so that the hourglass control cannot work.
The third approach is a method based on the idea presented in Belytschko et al. (1984).

If RI is used in a four-noded quadrilateral, the results are stabilized properly by means of the
so-called γ operatormethod. Analogically to the analysis performed in Belytschko et al. (1984)
for the Laplace equation, in this approach Rstab2 is taken into account as

Rstab2 (δǭ, ǭ)=
nel∑

e=1

∫

Be

δg̃Tχ2g̃ dV (3.13)

where χ2 is defined as follows

χ2 =
χh2e
2c

(3.14)

This coefficient is calculated according to the dimensional analysis in Commend et al. (2004).
The definitions of quantities χ and he are as previously, c is connected with the internal length
parameter l. The variable g̃ denotes a certain additional gradient connected with the averaged
strain field. Discretization must be introduced in such a way that it satisfies the condition
(Belytschko et al., 1984)

g̃=0 (3.15)

for any nodal values associated with a linear function ǭ. The rank of Kee, which is equal to 3
for RI (three positive eigenvalues), should be increased to 4 as is obtained for FI. Hence, after
Belytschko et al. (1984) operator γ is adopted, which should not influence the linear fields

γ
T = a[te− (t

T
exe)qx− (t

T
eye)qy] (3.16)

where in the above relation additional element quantities must be defined (Belytschko et al.,
1984)

t
T
e = [−1,1,−1,1] x

T
e = [xe1,xe2,xe3,xe4] y

T
e = [ye1,ye2,ye3,ye4] (3.17)

Theparameter a canbe equal to 1, because it is an arbitrary constant.Thevector te constitutes
in 2D a twisted form. The nodal coordinates are gathered in the vectors xe and ye, the second
index in their components refers to respectivenodes. It shouldbeemphasized that thederivatives
of the shape functions q for the averaged strain field are separated in Eq. (3.16) to form qx



846 J. Pamin, A.Wosatko

and qy. Analogically to the previous subsection, the discretization of the stabilization term
Rstab2 is performed as follows

g̃=γTe δg̃=γTδe (3.18)

where the discretized δg̃ can be treated as the weighting part P2(δe). The linearization of the
residual can easily be derived

g̃i+1 = g̃i+dg̃= g̃i+γ
Tde (3.19)

Therefore, the term Rstab2 has the form

Rstab2 =
nel∑

e=1

∫

Be

γχ2(g̃i+γ
Tde) dV (3.20)

For the γ operator method, the additional matrix and vector are introduced

K̃ee =

∫

Be

χ2γγ
T dV f̃e =

∫

Be

χ2γg̃i dV (3.21)

in order to obtain the final matrix equation to substitute Eq. (2.6)

[
Kaa+K̃aa Kae+K̃ae
Kea Kee+K̃ee

][
da
de

]
=

[
fext − fint − f̃

fǫ− fe− f̃e

]
(3.22)

The stabilization terms are needed in both variational equations, in order to ensure a proper
quality of the FE and stable numerical results. It can be shown that an analogical approach to
the stabilization of the equilibrium equations is equivalent to the one applied in the previous
subsection in the GLS method. However, this equivalence is valid only for the four-noded FE.
The GLS method seems to be more general, since the derivation is prepared for elements with
an arbitrary number of nodes.

4. Properties of four-noded element

In this section attention is focused on the description of the spectral analysis of a single FE,
whereRI (here one integration point)without orwith stabilization is considered.The properties
of elements including different types of interpolation and integration are described in detail in
Wosatko (2008). The single FE is subjected to tension in one direction and two loading phases
are considered: elasticity and damage. The computations of the eigenproblem for the tangent
operator K in these phases are performed at the end of chosen incremental steps. The following
material dataare adopted (units are irrelevant andhenceomitted):Young’smodulus E =20000,
Poisson’s ratio ν = 0.20, the damage model with threshold κo = 0.0001 and linear softening
where complete damage is for κu =0.002. The internal length parameter c equals 1.
The signs of eigenvalues are shown in Table 1 for the whole tangent operator K and for

particular submatrices Kaa and Kee. The accepted precision is equal to 1.0
−10, so that an

absolute eigenvalue less than this limit is assumed to be zero. Obviously, three zero eigenvalues
correspond to rigid motions of the element. However, whenRI without any hourglass control is
applied,more than three zero eigenvalues appear.Theyoriginatenot only from Kaa, butalso one
extra spurious zero eigenvalue comes from Kee, so altogether six such eigenvalues are present. A
single element withRI and stabilization described in the previous section is then examined. The
stabilization scaling factor χ is equal to 0.0001. The results in Table 1 show that in this case



Static gradient damage simulations... 847

Table 1. Signs of eigenvalues for FE with one integration point

(a) Without stabilization

Elasticity Damage

sign + 0 − + 0 −

Kaa 3 5 0 3 5 0

Kee 3 1 0 3 1 0

K 6 6 0 5 6 1

(b)With stabilization, χ=0.0001

Elasticity Damage

sign + 0 − + 0 −

Kaa+K̃aa 5 3 0 5 3 0

Kee+K̃ee 4 0 0 4 0 0

K+K̃ 9 3 0 8 3 1

the stabilization according to Eq. (3.22) ensures a proper spectrum of eigenvalues. A negative
eigenvalue for K appears after the peak during the damage progress and is related to softening
in the material model.

The analysis of eigenforms for all loading phases and considered cases can be found inWosat-
ko (2008). It is emphasized that the tangent operator K is non-symmetric, so imaginary parts of
components of eigenvectors in the spectral analysis are admitted. ForRIwithout stabilization in
the damage phase, non-zero values arise in one eigenvector for both subspaces: displacement and
averaged strain fields. The coupling in the formulation is visible for thismode and it corresponds
to the negative eigenvalue. Six zero eigenvalues of operator K are separated in such a way that
five of them originate from the part connected with submatrix Kaa and the last one – from
submatrix Kee.

When the stabilization is adopted, we observe that the spurious eigenforms for the displa-
cement field vanish. Two in-plane bendingmodes are related to the double eigenvalue equal to
8.2512. It is connectedwith the stabilization of equilibriumequations.Thecouplingof variational
equations can also be noticed for two eigenforms. Three eigenformswith zero eigenvalues correc-
tly correspond to the combination of rigidmotions. The eigenformwith the non-zero eigenvalue
and twisted mode in the averaged strain space proves that the stabilization of the averaging
equation is active. Moreover, the second field is stabilized by the γ method, where operator γ
for a rectangular element is reduced to the twist vector te.

5. Bar with imperfection

The simulation of uniaxial tension for a barwith an imperfection in themiddle is themost basic
localization test in physically nonlinearmechanics.Here the length of the bar is equal to 100mm
anddiscretization introduces 20FEs. Plane stress is assumedwith both thewidth and thickness
equal to 5mm. The material data are as follows: Young’s modulus E is 20000MPa, Poisson’s
ratio equals zero. For gradient damage, the internal length parameter is adopted c = 4mm2

(l=2.83mm) and the damage threshold is κo =0.0001. The threshold is reduced by 10% in the
middle of the bar to introduce imperfection. The damage growth equivalent to linear softening
is employed and the ultimate value κu = 0.0125. The load control and the arc length method
are used. Two integration schemes are confronted: full integration (FI) and reduced integration
(RI) with stabilization according to Eq. (3.22) and the coefficient χ=0.000001.

Generally, if FI is employed for the second or both discretized fields, the relation between
elongation u(L) and the stress σ is like for a slightly stiffer bar. The diagram for FI, presented
in Table 2, ends before the stress approaches zero, but it is connected with the fact that the
damagehistory parameter reaches the ultimate value κu andanunwanted change of the stiffness
or unloading is obtained. Two types of integration are compared in Table 2. Four or five cha-
racteristic steps are chosen according to the figures placed at the top of the table. A completely
different number of active integration points (ips) is observed as a consequence of the adopted



848 J. Pamin, A.Wosatko

integration in an FE. It is confirmed that the solution exhibits quadratic convergence in every
step for both integration cases.

Table 2.Convergence study – bar with imperfection. Two types of integration

Type FI: 2×2=4 ips RI + stabilization: 1 ips

Step No. Relative energy norm Active ips Relative energy norm Active ips

2 1.000000000000000E+00 8 1.000000000000000E+00 2
5.321187441983860E+02 8 5.763488320259804E+02 2
2.534990264376641E–04 8 6.048211715158885E–04 2
1.770959355358748E–10 8 1.304413629483436E–09 2
8.617827040306284E–23 8 6.008164724480290E–21 2

22 1.000000000000000E+00 40 1.000000000000000E+00 8
–9.779478750996960E–06 40 –3.149860222058139E–04 8
1.548406070066814E–10 40 –7.542830989034233E–12 8
3.871666510054424E–11 40 –1.885723762863035E–12 8
1.260922539766909E–19 40 4.655693122775198E–24 8

77 1.000000000000000E+00 12 1.000000000000000E+00 4
–6.737739431886433E–02 12 –8.738447627994261E–02 4
–1.936878713736846E–19 12 1.595955486024598E–16 4

–3.543616722096119E–27 4

122 1.000000000000000E+00 12 1.000000000000000E+00 4
–8.424036872173256E–02 12 –5.664603859365940E–01 4
–3.265570463637080E–20 12 6.614336369593637E–19 4

143 1.000000000000000E+00 4
–6.109566557762603E–01 4
3.162731767832675E–19 4

The averaged strain and damage distributions just after the peak and at the final state are
shown in Figs. 1 and 2, respectively. Both the variables have similar distributions just after the
peak load and also in the final state.

Fig. 1. Results for bar with imperfection – elements with FI



Static gradient damage simulations... 849

Fig. 2. Results for bar with imperfection – elements with RI and stabilization

6. Four point bending in concrete beam

The four-point bendingof a concrete beam is simulated.This benchmark is based on a reinforced
concrete beam, which was subjected to dynamic loading in experiments carried out by Eibl
et al.(1994) and then computed e.g. in Sluys (1994) and Dubé et al. (1996). The results of
simulations for the RC beam using the gradient damage model were presented in Wosatko
et al. (2006), and now we limit our analysis to the verification of two different combinations
of integration schemes as in the previous section. The geometry data are depicted in Fig. 3.
Two supports are introduced at the bottom of the beam, while the loading is imposed at two
points at the top. Plane stress conditions are assumed. Thematerial data are: Young’s modulus
E = 32940MPa, Poisson’s ratio νc = 0.2 and damage threshold κo = 95.6× 10

−6. In the
experiment, concrete was additionally reinforced byDramix fibers (1.2% volume). Themodified
vonMises definition is employed as the damage loading function (see de Vree et al., 1995). The
damage growth function asymptotically increases but never reaches 1 according to exponential
softening (Mazars andPijaudier-Cabot, 1989), thematerial ductility parameter η is equal to 350
and α = 0.96, which is responsible for residual stresses. Based on the results presented in
Wosatko et al. (2006), the internal length l = 16mm is adopted. Computations are performed
for three FEmeshes, namely: coarse – 56×8, medium – 112×16 and fine – 168×24 elements.
Displacement control is used. If one-point integration and stabilization is applied, the scaling
factor χ is equal to 0.0001.

Fig. 3. Four-point bending in concrete beam

From the comparison of diagrams in Fig. 4, the general tendency is observed for RI and
stabilizedFEsthata slightly smaller loadcarryingcapacity is obtained than for thecomputations
with FI. As expected for RI, the simulated behaviour gives a response as for a slightly less stiff
beam. The general character of diagrams does not differ significantly. The diagrams for the fine
mesh are the nearest to each other.

FromFig. 5 with the deformations of themeshes and Fig. 6 with the corresponding contour
plots of the averaged strainmeasure, it is noticed that two zones of localization are simulated.An
exception occurs for the case with coarse mesh and FI. Of course these zones can be associated
with cracking of the beam. It is also confirmed in Figs. 5a and 6a that, among the calculated
cases, the solution for the coarse mesh and full integration leads to the stiffest model. However,



850 J. Pamin, A.Wosatko

Fig. 4. Load-deflection diagrams – different meshes and integration schemes

Fig. 5. Deformation – influence of meshes and integration schemes

Fig. 6. Averaged strain distribution – different meshes and integration schemes

according to the plastic hinges theory in bars, both forms of deformation, both crack patterns
and hence one or two localization zones are admissible. Moreover, mesh-independent results are
obtained for both integration schemes.

The last aspect of our analysis is the verification of the impact of the parameter χ on the
obtained results. A proper value of the scaling factor is necessary to ensure a numerically stable
response as mentioned in Section 3. On the other hand, the parameter should have a possibly
small value. In previous computations in this example, χ equals 0.0001. Nowwe vary this value
from the smallest 0.00000001 up to the largest 1.0 with multiplier 100 for each case.

Thediagram inFig. 7 and the deformation of themesh inFig. 8a show that a too small value
of the scaling factor does not guarantee proper stabilization for one-point integration in FE and
even hourglass patterns can appear. For the case with a large value, for example χ= 1.0, the
response is considerably stiffer, although two localization zones are still obtained.



Static gradient damage simulations... 851

Fig. 7. Load-deflection diagrams – influence of scaling factor χ for medium mesh

Fig. 8. Deformation for medium mesh – influence of scaling factor χ

7. Conclusions

When one-point reduced integration is applied in the four-noded FE then, for effective compu-
tations, a stabilization must be introduced in the formulation.

Moreover, as shown inSection 3, for gradient damagefinite elementswith linear interpolation
and one integration point, hourglass controlmust be incorporated for both primary fields. In the
spectral analysiswithRIadditional spurious eigenforms,which appear for thedisplacements and
for the averaged strain measure, confirm that hourglass control (mesh-stabilization) is needed.
The stabilization terms are derived using the GLS method for the equilibrium equations and
the γ operator method for the averaging equation. This is because, as discussed in Section 3.2,
neither the GLS nor the gradient GLS methods provide effective stabilization for the latter
equation.

The analysis of a bar with an imperfection shows that the solution guarantees quadratic
convergence when RI with stabilization is used. The reduced integration results in a slightly
less stiff response of the model, which is consistent with what is known about the influence
of integration quadrature on FE results. This is also observed for the simulation of a concrete
beam in four-point bending. The influence of the stabilization scaling factor χ on the results is
examined.

The research is planned to be continued towards a 3D formulation. The elaboration of the
formulationwith an effective stabilization for the 3Dbrick element requires further research, but
the derivation in Section 3 can constitute its initial point. It turns out that for 2D simulations
the computation cost reduction in comparison to full integration is smaller than expected, but
the implementation in FEAP package (Taylor, 2001) has not been optimized for the execution
time. It is expected that the gain in 3D simulations would be much larger, even though the
solver efficiency is then the main computation cost factor.

Acknowledgment

The researchwas supportedby theEuropeanUnion through theEuropeanSocialFundwithin project

Cracow University of Technology development program – top quality teaching for the prospective Polish

engineers University of the 21st century (contract no.UDA-POKL.04.01.01-00-029/10-00).



852 J. Pamin, A.Wosatko

References

1. AskesH., Pamin J., deBorstR., 2000,Dispersion analysis and element-freeGalerkin solutions
of second- and fourth-order gradient-enhanced damage models, Int. J. Numer. Meth. Engng, 49,
811-832

2. Belytschko T., Ong J.S.-J., LiuW.K., Kennedy J.M., 1984,Hourglass control in linear and
nonlinear problems,Comput. Methods Appl. Mech. Engrg., 43, 251-276

3. Commend S., Truty A., Zimmermann T., 2004, Stabilized finite elements applied to elasto-
plasticity: I. mixed displacement-pressure fromulation,Comput. Methods Appl. Mech. Engrg., 193,
3559-3586

4. de Borst R., Pamin J., Geers M.G.D., 1999, On coupled gradient-dependent plasticity and
damage theories with a view to localization analysis,Eur. J. Mech. A/Solids, 18, 6, 939-962

5. de Borst R., Sluys L.J., Mühlhaus H.-B., Pamin J., 1993, Fundamental issues in finite
element analyses of localization of deformation,Eng. Comput., 10, 99-121

6. de Vree J.H.P., Brekelmans W.A.M., van Gils M.A.J., 1995, Comparison of nonlocal ap-
proaches in continuum damagemechanics,Comput. and Struct., 55, 4, 581-588

7. Dubé J.-F., Pijaudier-CabotG., LaBorderie Ch., 1996,Rate dependent damagemodel for
concrete in dynamics, J. of Engng. Mech., 122, 10, 939-947

8. Eibl J., Bischoff P.H., Lohrmann G., 1994, Failure mechanics of fibre-reinforced concre-
te and pre-damaged structures: dynamic loading conditions, Technical Report BRITE/EURAM
P-89-3275, Univ. of Karlsruhe, Karlsruhe, Germany

9. Geers M.G.D., 1997,Experimental Analysis and Computational Modelling of Damage and Frac-
ture, Ph.D. dissertation, EindhovenUniversity of Technology, Eindhoven

10. Geers M.G.D., de Borst R., Peerlings R.H.J., 2000, Damage and crackmodeling in single-
edge and double-edge notched concrete beams,Eng. Fract. Mech., 65, 247-261

11. Harari I., Frey S., Franca L.P., 2002, A note on a recent study of stabilized finite element
computations for heat conduction,Computational Mechanics, 28, 63-65

12. Lemaitre J., 1971, Evaluation of dissipation and damage in metals, [In:] Proc. I.C.M., vol. 1,
Kyoto, Japan

13. Mazars J., Pijaudier-Cabot G., 1989, Continuum damage theory – application to concrete,
ASCE J. Eng. Mech., 115, 345-365

14. Pamin J., Wosatko A., Winnicki A., 2003, Two- and three-dimensional gradient damage-
plasticity simulations of cracking in concrete, [In:] N. Bićanić et al. (Edit.), Proc. EURO-C 2003
Int. Conf. Computational Modelling of Concrete Structures, Rotterdam/Brookfield,A.A. Balkema,
325-334,

15. PastorM.,QuecedoM., ZienkiewiczO.C., 1997,Amixeddisplacement-pressure formulation
for numerical analysis of plastic failure,Comput. and Struct., 62, 1, 13-23

16. Peerlings R.H.J., de Borst R., Brekelmans W.A.M., de Vree J.H.P., 1996, Gradient-
enhanced damage for quasi-brittle materials, Int. J. Numer. Meth. Engng, 39, 3391-3403

17. Peerlings R.H.J., de Borst R., Brekelmans W.A.M., Geers M.G.D., 1998, Gradient-
enhanced damagemodelling of concrete fracture,Mech. Cohes.-frict. Mater., 3, 323-342

18. Reese S., 2003,On a consistent hourglass stabilization technique to treat large inelastic deforma-
tions and thermo-mechanical coupling in plane stress problems, Int. J. Numer. Meth. Engng, 57,
1095-1127

19. Simone A., Askes H., Peerlings R.H.J., Sluys L.J., 2003, Interpolation requirements for
implicit gradient-enhanced continuum damagemodels,Communications in Numerical Methods in
Engineering, 19, 563-572 (see also Simone et al., 2004)



Static gradient damage simulations... 853

20. SimoneA.,AskesH., PeerlingsR.H.J., SluysL.J., 2004,Corrigendum’Interpolation rquire-
ments for implicit gradient-enhanced continuum damage models’, Communications in Numerical
Methods in Engineering, 20, 163-165

21. Sluys L.J., 1994, Dynamic failure in reinforced concrete structures, [In:] G.M.A. Kusters and
M.A.N. Hendriks (Edit.),DIANAComputational Mechanics ’94, 193-203,Dordrecht, KluwerAca-
damic Publishers

22. Taylor R.L., 2001, FEAP – A Finite Element Analysis Program, Version 7.4, User manual,
University of California at Berkeley, Berkeley

23. Wosatko A., 2008, Finite-Element Analysis of Cracking in Concrete Using Gradient Damage-
Plasticity, Ph.D. dissertation, CracowUniversity of Technology, Cracow

24. WosatkoA., Pamin J.,Winnicki A., 2006,Gradient damage in simulations of behaviour of RC
bars and beams under static and impact loading,Archives of Civil Engineering, 52, 1, 455-477

25. Zienkiewicz O.C., Taylor R.L., Zhu J.Z., 2005, The Finite Element Method: Its Basis and
Fundamentals, Elsevier Butterworth-Heinemann, sixth ed.

Gradientowy model do statycznych symulacji uszkodzenia przy użyciu elementów

skończonych ze stabilizacją

Streszczenie

Wartykuleprzedstawionodwupolowe elementy skończone z kontroląpasożytniczych formdeformacji,
sformułowane dla sprzężonego problemu gradientowej mechaniki uszkodzenia. Zastosowano stabilizację
równań równowagi zgodnie z metodą najmniejszych kwadratów. Rozpatrzono trzy warianty stabilizacji
dodatkowego równania uśredniającego, z których tylko metoda wykorzystująca operator γ jest skutecz-
na. Uwaga została skupiona na wyprowadzeniu, implementacji i analizie spektralnej czterowęzłowego
elementu do symulacji dwuwymiarowych. Przedyskutowano wyniki uzyskane dla podstawowych testów
rozciąganego pręta z imperfekcją i belki czteropunktowo zginanej.

Manuscript received March 1, 2012; accepted for print March 26, 2012