OOFEM

Material Model Library Manual
Bořek Patzák
Czech Technical University

Faculty of Civil Engineering
Department of Structural Mechanics
Thákurova 7, 166 29 Prague, Czech Republic
March 19, 2019
1
Contents
1 Material Models for Structural Analysis 3
1.1 Elastic materials . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.1 Isotropic linear elastic material - IsoLE . . . . . . . . . . 3
1.1.2 Orthotropic linear elastic material - OrthoLE . . . . . . . 3
1.1.3 General anisotropic linear elastic material - AnisoLE . . . 4
1.1.4 Hyperelastic material - HyperMat . . . . . . . . . . . . . 4
1.1.5 Hyperelastic material - Compressible Mooney-Rivlin . . . 5
1.2 Winkler-Pasternak model . . . . . . . . . . . . . . . . . . . . . . 6
1.3 Large-strain master material . . . . . . . . . . . . . . . . . . . . 7
1.4 Plasticity-based material models . . . . . . . . . . . . . . . . . . 7
1.4.1 Drucker-Prager model - DruckerPrager . . . . . . . . . . . 7
1.4.2 Drucker-Prager model with tension cut-off and isotropic
damage - DruckerPragerCut . . . . . . . . . . . . . . . . . 9
1.4.3 Mises plasticity model with isotropic damage - MisesMat 12
1.4.4 Mises plasticity model with isotropic damage, nonlocal -
MisesMatNl, MisesMatGrad . . . . . . . . . . . . . . . . . 13
1.4.5 Rankine plasticity model with isotropic damage and its
nonlocal formulations - RankMat, RankMatNl, RankMat-
Grad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.4.6 Perfectly plastic material with Mises yield condition - Steel1 16
1.4.7 Composite plasticity model for masonry - Masonry02 . . . 16
1.4.8 Nonlinear elasto-plastic material model for concrete plates
and shells - Concrete2 . . . . . . . . . . . . . . . . . . . . 23
1.5 Material models for tensile failure . . . . . . . . . . . . . . . . . . 23
1.5.1 Nonlinear elasto-plastic material model for concrete plates
and shells - Concrete2 . . . . . . . . . . . . . . . . . . . . 23
1.5.2 Smeared rotating crack model - Concrete3 . . . . . . . . . 23
1.5.3 Smeared rotating crack model with transition to scalar
damage - linear softening - RCSD . . . . . . . . . . . . . 25
1.5.4 Smeared rotating crack model with transition to scalar
damage - exponential softening - RCSDE . . . . . . . . . 25
1.5.5 Nonlocal smeared rotating crack model with transition to
scalar damage - RCSDNL . . . . . . . . . . . . . . . . . . 26
1.5.6 Isotropic damage model for tensile failure - Idm1 . . . . . 27
1.5.7 Nonlocal isotropic damage model for tensile failure - Idmnl1 34
1.5.8 Anisotropic damage model - Mdm . . . . . . . . . . . . . 39
1.5.9 Isotropic damage model for interfaces . . . . . . . . . . . 43
1.5.10 Isotropic damage model for interfaces using tabulated data
for damage . . . . . . . . . . . . . . . . . . . . . . . . . . 43
1.6 Material models specific to concrete . . . . . . . . . . . . . . . . 44
1.6.1 Mazars damage model for concrete - MazarsModel . . . . 44
1.6.2 Nonlocal Mazars damage model for concrete - MazarsMod-
elnl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
1.6.3 CebFip78 model for concrete creep with aging - CebFip78 45
1.6.4 Double-power law model for concrete creep with aging -
DoublePowerLaw . . . . . . . . . . . . . . . . . . . . . . . 45
1.6.5 Eurocode 2 model for concrete creep and shrinkage - EC2CreepMat 45
1.6.6 B3 and MPS models for concrete creep with aging . . . . 49
2
1.6.7 MPS damage model . . . . . . . . . . . . . . . . . . . . . 56
1.6.8 Microplane model M4 - Microplane M4 . . . . . . . . . . 66
1.6.9 Damage-plastic model for concrete - ConcreteDPM . . . . 66
1.6.10 CDPM2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
1.6.11 Fixed crack model for concrete - ConcreteFCM . . . . . . 75
1.6.12 Fixed crack model for fiber reinforced composites - FR-
CFCM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
1.6.13 “Nonlocal” model for SHCC . . . . . . . . . . . . . . . . . 81
1.7 Orthotropic damage model with fixed crack orientations for com-
posites - CompDamMat . . . . . . . . . . . . . . . . . . . . . . . 82
1.8 Orthotropic elastoplastic model with isotropic damage - Trab-
Bone3d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
1.8.1 Local formulation . . . . . . . . . . . . . . . . . . . . . . . 84
1.8.2 Nonlocal formulation - TrabBoneNL3d . . . . . . . . . . . 86
1.9 Material models for interfaces . . . . . . . . . . . . . . . . . . . . 91
1.9.1 Cohesive interface material - cohInt . . . . . . . . . . . . 91
1.9.2 Isotropic damage model for interfaces . . . . . . . . . . . 92
1.9.3 Simple interface material - obsolete . . . . . . . . . . . . . 92
1.9.4 Bond-slip model for reinforced concrete - BondCEB . . . 92
1.10 Material models for lattice elements . . . . . . . . . . . . . . . . 93
1.10.1 Latticedamage2d . . . . . . . . . . . . . . . . . . . . . . . 94
1.11 Material models for steel relaxation . . . . . . . . . . . . . . . . . 95
1.11.1 Model for relaxation of prestressing steel - SteelRelaxMat 95
2 Material Models for Transport Problems 98

2.1 Isotropic linear material for heat transport – IsoHeat . . . . . . . 98
2.2 Isotropic linear material for moisture transport – IsoLinMoisture 98
2.3 Isotropic material for moisture transport based on Bažant and
Najjar – BazantNajjarMoisture . . . . . . . . . . . . . . . . . . . 98
2.4 Nonlinear isotropic material for moisture transport – NlIsoMoisture100
2.5 Material for cement hydration - CemhydMat . . . . . . . . . . . 101
2.6 Material for cement hydration - HydratingConcreteMat . . . . . 105
2.7 Coupled heat and mass transport material model - HeMotk . . . 107
2.8 Coupled heat and mass transport material model - HeMoKunzel 109
2.9 Material for unsaturated flow in lattice models – LatticeTransMat 111
2.9.1 One-dimensional transport element . . . . . . . . . . . . . 111
2.9.2 Constitutive laws . . . . . . . . . . . . . . . . . . . . . . . 112
3 Material Models for Fluid Dynamic 114

3.1 Newtonian fluid - NewtonianFluid . . . . . . . . . . . . . . . . . 114
3.2 Bingham fluid - BinghamFluid . . . . . . . . . . . . . . . . . . . 114
3.3 Two-fluid material - TwoFluidMat . . . . . . . . . . . . . . . . . 115
3.4 FE2 fluid - FE2FluidMaterial . . . . . . . . . . . . . . . . . . . . 115
List of Tables
1 Linear Isotropic Material - summary. . . . . . . . . . . . . . . . . 3
2 Orthotropic, linear elastic material – summary. . . . . . . . . . . 5
3 Anisotropic, linear elastic material – summary. . . . . . . . . . . 5
3
4 Hyperelastic material - summary. . . . . . . . . . . . . . . . . . . 6
5 Compressible Mooney-Rivlin - summary. . . . . . . . . . . . . . . 7
6 Winkler Pasternak material - summary. . . . . . . . . . . . . . . 7
7 Large-strain master material material - summary. . . . . . . . . . 8
8 DP material - summary. . . . . . . . . . . . . . . . . . . . . . . . 10
9 Drucker Prager material with tension cut-off - summary. . . . . . 11
10 Mises plasticity – summary. . . . . . . . . . . . . . . . . . . . . . 13
11 Nonlocal integral Mises plasticity – summary. . . . . . . . . . . . 15
12 Gradient-enhanced Mises plasticity – summary. . . . . . . . . . . 16
13 Rankine plasticity – summary. . . . . . . . . . . . . . . . . . . . 17
14 Nonlocal integral Rankine plasticity – summary. . . . . . . . . . 18
15 Gradient-enhanced Rankine plasticity – summary. . . . . . . . . 19
16 Perfectly plastic material with Mises condition – summary. . . . 19
17 Composite model for masonry - summary. . . . . . . . . . . . . . 23
18 Nonlinear elasto-plastic material model for concrete - summary. . 24
19 Rotating crack model for concrete - summary. . . . . . . . . . . . 25
20 RC-SD model for concrete - summary. . . . . . . . . . . . . . . . 26
21 RC-SD model for concrete - summary. . . . . . . . . . . . . . . . 26
22 RC-SD-NL model for concrete - summary. . . . . . . . . . . . . . 27
23 Isotropic damage model for tensile failure – summary. . . . . . . 34
24 Nonlocal isotropic damage model for tensile failure – summary. . 37
25 Nonlocal isotropic damage model for tensile failure – continued. . 38
26 Basic equations of microplane-based anisotropic damage model . 39
27 MDM model - summary. . . . . . . . . . . . . . . . . . . . . . . . 42
28 Isotropic damage model for interface elements – summary. . . . . 44
29 Isotropic damage model for interface elements using tabulated
data for damage – summary. . . . . . . . . . . . . . . . . . . . . 45
30 Mazars damage model – summary. . . . . . . . . . . . . . . . . . 46
31 Nonlocal Mazars damage model – summary. . . . . . . . . . . . . 47
32 CebFip78 material model – summary. . . . . . . . . . . . . . . . 47
33 Double-power law model – summary. . . . . . . . . . . . . . . . . 48
34 EC2Creep material model – summary. . . . . . . . . . . . . . . . 59
35 B3 creep and shrinkage model – summary. . . . . . . . . . . . . . 60
36 B3solid creep and shrinkage model – summary. . . . . . . . . . . 62
37 MPS theory—summary. . . . . . . . . . . . . . . . . . . . . . . . 64
38 MPS damage–summary. . . . . . . . . . . . . . . . . . . . . . . . 65
39 Microplane model M4 – summary. . . . . . . . . . . . . . . . . . 66
40 Damage-plastic model for concrete – summary. . . . . . . . . . . 70
42 Fixed crack model for concrete – summary. . . . . . . . . . . . . 77
43 Fixed crack model for fiber reinforced concrete – summary. . . . 81
44 Nonlocal fixed crack model for fiber reinforced concrete – summary. 81
41 CDPM2 – summary. . . . . . . . . . . . . . . . . . . . . . . . . . 87
45 Orthotropic damage model with fixed crack orientations for com-
posites – summary. . . . . . . . . . . . . . . . . . . . . . . . . . . 88
46 Anisotropic elastoplastic model with isotropic damage - summary. 89
47 Nonlocal formulation of anisotropic elastoplastic model with isotropic
damage – summary. . . . . . . . . . . . . . . . . . . . . . . . . . 90
48 Cohesive interface material – summary. . . . . . . . . . . . . . . 91
49 Simple interface material – summary. . . . . . . . . . . . . . . . . 92
50 Bond-slip model for reinforced concrete – summary. . . . . . . . 94
4
51 Scalar damage model for 2d lattice elements – summary. . . . . . 96
52 SteelRelaxMat material model – summary. . . . . . . . . . . . . . 97
53 Linear isotropic material for heat transport - summary. . . . . . 98
54 Linear isotropic material for moisture transport - summary. . . . 99
55 Nonlinear isotropic material for moisture transport - summary. . 99
56 Nonlinear isotropic material for moisture transport - summary. . 103
57 Cemhydmat - summary. . . . . . . . . . . . . . . . . . . . . . . . 104
58 HydratingConcreteMat - summary of affinity hydration models. . 106
59 Coupled heat and mass transfer material model - summary. . . . 108
60 Parameters from Kunzel’s model. . . . . . . . . . . . . . . . . . . 109
61 Coupled heat and mass transfer material model Kunzel - summary.110
62 Material for unsaturated flow in lattice models - summary. . . . . 113
63 Newtonian Fluid material - summary. . . . . . . . . . . . . . . . 114
64 Bingham Fluid material - summary. . . . . . . . . . . . . . . . . 115
65 Two-Fluid material - summary. . . . . . . . . . . . . . . . . . . . 115
66 FE2 fluid material - summary. . . . . . . . . . . . . . . . . . . . . 116
5
1 Material Models for Structural Analysis
1.1 Elastic materials
1.1.1 Isotropic linear elastic material - IsoLE
Linear isotropic material model. The model parameters are summarized in
Tab. 1.
Description Linear isotropic elastic material

Record Format IsoLE num(in) # d(rn) # E(rn) # n(rn) # tAlpha(rn) #
Parameters - num material model number
- d material density
- E Young modulus
- n Poisson ratio
- tAlpha thermal dilatation coefficient
Supported modes 3dMat, PlaneStress, PlaneStrain, 1dMat, 2dPlateLayer,
2dBeamLayer, 3dShellLayer, 2dPlate, 2dBeam, 3dShell,
3dBeam, PlaneStressRot
Features Adaptivity support
Table 1: Linear Isotropic Material - summary.
1.1.2 Orthotropic linear elastic material - OrthoLE

Orthotropic, linear elastic material model. The model parameters are summa-
rized in Tab. 2. Local coordinate system, which determines axes of material
orthotrophy can by specified using lcs array. This array contains six numbers,
where the first three numbers represent directional vector of a local x-axis, and
next three numbers represent directional vector of a local y-axis. The local z-
axis is determined using the vector product. The right-hand coordinate system
is assumed.
Local coordinate system can also be specified using scs parameter. Then
local coordinate system is specified in so called “shell” coordinate system, which
is defined locally on each particular element and its definition is as follows:
principal z-axis is perpendicular to shell mid-section, x-axis is perpendicular to
z-axis and normal to user specified vector (so x-axis is parallel to plane, with
being normal to this plane) and y-axis is perpendicular both to x and z axes.
This definition of coordinate system can be used only with plates and shells
elements. When vector is parallel to z-axis an error occurs. The scs array
contain three numbers defining direction vector . If no local coordinate system
is specified, by default a global coordinate system is used.
For 3D case the material compliance matrix has the following form
 
1/EX −νxy /Ex −νxz /Ex 0 0 0
 −νyx /Ey 1/Ey −νyz /Ey 0 0 0 
 
 −νzx /Ez −νzy /Ez 1/E z 0 0 0 
C=   (1)
 0 0 0 1/G yz 0 0 

 0 0 0 0 1/Gxz 0 
0 0 0 0 0 1/Gxy
6
By inversion, the material stiffness matrix has the form
 
dxx dxy dxz 0 0 0

 dyy d yz 0 0 0 

 sym dzz 0 0 0 
D=  0
 (2)
 0 0 Gyz 0 0 

 0 0 0 0 Gxz 0 
0 0 0 0 0 Gxy
where ξ = 1 − (νxy ∗ νyx + νyz ∗ νzy + νzx ∗ νxz ) − (νxy ∗ νyz ∗ νzx + νyx ∗ νzy ∗ νxz )
and
dxx = EX (1 − νyz ∗ νzy )/ξ, (3)

dxy = Ey ∗ (νxy + νxz ∗ νzy )/ξ, (4)
dxz = Ez ∗ (νxz + νyz ∗ νxy )/ξ, (5)
dyy = Ey ∗ (1 − νxz ∗ νzx )/ξ, (6)
dyz = Ez ∗ (νyz + νyx ∗ νxz )/ξ, (7)
dzz = Ez ∗ (1 − νyx ∗ νxy )/ξ. (8)
Ei is Young’s modulus in the i-th direction, Gij is the shear modulus in

ij plane, νij is the major Poisson ratio, and νji is the minor Poisson ratio.
Assuming that Ex > Ey > Ez , νxy > νyx etc., then νxy is referred to as the
major Poisson ratio, while νyx is referred as the minor Poisson ratio. Note
that there are only nine independent material parameters, because of symmetry
conditions. The symmetry condition yields
νxy Ey = νyx Ex , νyz Ez = νzy Ey , νzx Ex = νxz Ez
The model description and parameters are summarized in Tab. 2.
1.1.3 General anisotropic linear elastic material - AnisoLE

Linear elastic material model with completely general material stiffness (21 in-
dependent elastic constants). The model parameters are summarized in Tab. 3.
The material stiffness matrix is a completely arbitrary symmetric 6 × 6 matrix.
The input line must contain an array with the upper triangle of this matrix.
This array has length 21 and the stiffness coefficients are listed in each row from
the diagonal to the last column.
1.1.4 Hyperelastic material - HyperMat

This material model can describe elastic behavior at large strains. A hypere-
lastic model postulates the existence of free energy potential. Existence of the
potential implies reversibility of deformations and no energy dissipation during
loading process. Here we use the free energy function introduced in [20]

1 2
K − G J 2 − 2lnJ − 1 + G (E : I − lnJ)

ρ0 ψ = (9)
4 3
where K is the bulk modulus, G is the shear modulus, J is the Jacobian (deter-
minant of the deformation gradient, corresponding to the ratio of the current
7
Description Orthotropic, linear elastic material
Record Format OrthoLE num(in) # d(rn) # Ex(rn) # Ey(rn) # Ez(rn) #
NYyz(rn) # NYxz(rn) # NYxy(rn) # Gyz(rn) # Gxz(rn) #
Gxy(rn) # tAlphax(rn) # tAlphay(rn) # tAlphaz(rn) #
[ lcs(ra) #] [ scs(ra) #]
- Ex, Ey, Ez Young moduli for x,y, and z directions
- NYyz, NYxz, NYxy major Poisson’s ratio coefficients
- Gyz, Gxz, Gxy shear moduli
- tAlphax, tAlphay, tAlphaz thermal dilatation coefficients
in x,y,z directions
- lcs Array defining local material x and y axes of orthotro-
phy
- scs Array defining a normal vector n. The local x axis is
parallel to plane with n being plane normal. The material
local z-axis is perpendicular to shell mid-section.
2dBeamLayer, 3dShellLayer, 2dPlate, 2dBeam, 3dShell,
3dBeam, PlaneStressRot
Table 2: Orthotropic, linear elastic material – summary.
Description Anisotropic, linear elastic material

Record Format OrthoLE num(in) # d(rn) # stiff(ra) # tAlphax(ra) #
- stiff real array of length 21 with stiffness coefficients D11 ,
D12 , D13 . . . D66
- tAlpha real array of length 0 or 3 with thermal dilatation
coefficients in x,y,z directions
Supported modes 3dMat, PlaneStress, PlaneStrain, 1dMat
Table 3: Anisotropic, linear elastic material – summary.
and initial volume) and E is the Green-Lagrange strain. Then stress-strain law
can be derived from (9) as

∂ψ 1 2
K − G J 2 − 1 C −1 + G I − C −1

S = ρ0 = (10)
∂E 2 3
where S is the second Piola-Kirchhoff stress, E is the Green-Lagrange strain
and C is the right Cauchy-Green tensor. The model description and parameters
are summarized in Tab. 4.
1.1.5 Hyperelastic material - Compressible Mooney-Rivlin

The Mooney-Rivlin strain energy function is expressed by
1
ρ0 ψ = C1 (I¯1 − 3) + C2 (I¯2 − 3) + K(lnJ)2 (11)
2
8
Description Hyperelastic material
Record Format HyperMat (in) # d(rn) # K(rn) # G(rn) #
Parameters - material number
- K bulk modulus
- G shear modulus
Supported modes 3dMat
Table 4: Hyperelastic material - summary.
where C1 and C2 are material constants, K is the bulk modulus, J is the Ja-
cobian (determinant of the deformation gradient, corresponding to the ratio of
2 2
the current and initial volume), I¯1 = J − 3 I1 , I¯2 = J − 3 I2 , where I1 and I2 are
the first and the second principal invariants of the right Cauchy-Green deforma-
tion tensor C. Compressible neo-Hookean material model is obtained by setting
C2 = 0. Then stress-strain law can be derived from (11) as
∂ψ ∂ I¯1 ∂ I¯2
P = ρ0 = C1 + C2 + KlnJF −T (12)
∂F ∂F ∂F
where P is the first Piola-Kirchhoff stress,
∂ I¯1 2 2
= 2 F − I¯1 F −T (13)
∂F J 3 3
and
∂ I¯2 4 2
= 2I¯1 F − I¯2 F −T − 4 F · C (14)
∂F 3 J3
The first elasticity tensor is derived as
∂Pij −1 −1 −1 −1
Aijkl = = C1 A1ijkl + C2 A2ijkl + K(Fji Flk − lnJFjk Fli ) (15)
∂Fkl
where

2 −2 −1 −1 −1 2 −1 −1 −1
A1ijkl = J 3 3δik δjl + I1 Fjk Fli − 2Flk Fij + I1 Fji Flk − 2Fji Fkl
3 3
(16)
and

− 34 4 −1 8 −1 −1 4 −1
A2ijkl
= 2J I1 δik δjl + 2Fij Fkl − I1 Fij Flk − I2 Fji Flk − I1 Fji Fkl
3 9 3

4 −1 2 −1 −1 4 −1
+ Fkn Cnl Fji + I2 Fli Fjk + Fkl Fim Cmj − δik Clj − Fil Fkj + Fkm Fim δjl (17)
3 3 3
1.2 Winkler-Pasternak model

Implementation of 2D winkler-pasternak model for plate (and potentially beam)
subsoil model.
9
Description Mooney-Rivlin
Record Format MooneyRivlin (in) # d(rn) # K(rn) # C1(rn) # C2(rn) #
- K bulk modulus
- C1 material constant
- C2 material constant
Supported modes 3dMat, PlaneStrain
Table 5: Compressible Mooney-Rivlin - summary.
Description Winkler Pasternak isotropic material for subsoil interaction

Record Format winklerpasternak (in) # d(rn) # c1(rn) # c2(rn) #
- c1 C1 Winkler parameter (foundation modulus)
- c2 C2 shear interaction modulus
Supported modes 2dPlateSubSoil
Table 6: Winkler Pasternak material - summary.
1.3 Large-strain master material

In this section, a large-strain master model based on generalized stress-strain
measures is described. In the first step, strain measure is computed from equa-
tion (18).
1


 (C m − I) , if m 6= 0
 2m

E (m) = (18)

 1
 ln C, if m = 0

2
where I is the second-order unit tensor and C = F T F is Cauchy-Green strain
tensor. In the special cases when m = 0 and m = 0.5 we obtain the so-called
Hencky (logarithmic) and Biot tensor, while for m = 1 we obtain the right
Green-Lagrange strain tensor. In the second step, this strain measure enters a
constitutive law of slave material and the stress measure conjugated to the strain
measure defined in step one and appropriate stiffness matrix are computed. In
the third step, the generalized stress tensor and stiffness matrix are transformed
into the second Piola-Kirchhoff stress and the appropriate stiffness tensor. The
model description and parameters are summarized in Tab. 7.
1.4 Plasticity-based material models

1.4.1 Drucker-Prager model - DruckerPrager
The Drucker-Prager plasticity model1 is an isotropic elasto-plastic model based
on a yield function
f (σ, τY ) = F (σ) − τY (19)
1 Contributed by Simon Rolshoven, LSC, FENAC, EPFL.
10
Description Large-strain master material material
Record Format LSmasterMat (in) # m(rn) # slavemat(in) #
- m parameter defining the strain measure
- slavemat number of slave material
Supported modes 3dMatF
Table 7: Large-strain master material material - summary.
with the pressure-dependent equivalent stress

p
F (σ) = αI1 + J2 (20)
As usual, σ is the stress tensor, τY is the yield stress under pure shear, and I1
and J2 are the first invariant and second deviatoric invariant of the stress tensor.
The friction coefficient α is a positive parameter that controls the influence of
the pressure on the yield limit, important for cohesive-frictional materials such
as concrete, soils or other geomaterials. Regarding Mohr-Coulomb plasticity,
relation to cohesion, c, and the angle of friction, θ, exists for the Drucker-Prager
model
2 sin θ
α= √ (21)
(3 − sin θ) 3
6c cos θ
τY = √ (22)
(3 − sin θ) 3
The flow rule is derived from the plastic potential

p
g (σ) = αψ I1 + J2 (23)
where αψ is the dilatancy coefficient. An associated model with αψ = α would

overestimate the dilatancy of concrete, so the dilatancy coefficient is usually
chosen smaller than the friction coefficient. The model is described by the
equations
σ = D : (ε − εp ) (24)
τY = h(κ) (25)

∂g s
ε̇p = λ̇ = λ̇ αψ δ + √ (26)
∂σ 2 J2
r
2
κ̇ = kε̇p k (27)
3
λ̇ ≥ 0, f (σ, τY ) ≤ 0, λ̇ f (σ, τY ) = 0 (28)

which represent the linear elastic law, hardening law, evolution laws for plastic
strain and hardening variable, and the loading-unloading conditions. In the
above, D is the elastic stiffness tensor, ε is the strain tensor, εp is the plastic
strain tensor, λ is the plastic multiplier, δ is the unit second-order tensor, s
is the deviatoric stress tensor, κ is the hardening variable, and a superior dot
11
marks the derivative with respect to time. The flow rule has the form given
in Eq. (26) at all points of the conical yield surface with the exception of its
vertex, located on the hydrostatic axis.
For the present model, the evolution of the hardening variable can be explic-
itly linked to the plastic multiplier. Substituting the flow rule (26) into Eq. (27)
and computing the norm leads to
κ̇ = k λ̇ (29)
q
with a constant parameter k = 1/3 + 2αψ 2 , so the hardening variable is pro-
portional to the plastic multiplier. For α = αψ = 0, the associated J2 -plasticity

model is recovered as a special case.
In the simplest case of linear hardening, the hardening function is a linear
function of κ, given by
h(κ) = τ0 + HEκ (30)
where τ0 is the initial yield stress, and H is the hardening modulus normalized
with the elastic modulus. Alternatively, an exponential hardening function
h(κ) = τlimit + (τ0 − τlimit ) e−κ/κc (31)
can be used for a more realistic description of hardening.

The stress-return algorithm is based on the Newton-iteration. In plasticity,
this is commonly called Closest-Point-Projection (CPP), and it generally leads
to quadratic convergence. The implemented algorithm is convergent in any
stress case, but in the vicinity of the vertex region, quadratic convergence might
be lost because of insufficient regularity of the yield function.
The algorithmic tangent stiffness matrix is implemented for both the reg-
ular case and the vertex region. Generally, the error decreases quadratically
(of course only asymptotically). Again, in the vicinity of the vertex region,
quadratic convergence might be lost due to insufficient regularity. Furthermore,
the tangent stiffness matrix does not always exist for the vertex case. In these
cases, the elastic stiffness is used instead. It is generally safer (but slower) to
use the elastic stiffness if you encounter any convergence problems, especially if
your problem is tension-dominated.
1.4.2 Drucker-Prager model with tension cut-off and isotropic dam-

age - DruckerPragerCut
The Drucker-Prager plasticity model with tension cut-off is a multisurface model,
appropriate for cohesive-frictional materials such as concrete loaded both in
compression and tension. The plasticity model is formulated for isotropic hard-
ening and enhanced by isotropic damage, which is driven by the cumulative
plastic strain. The model can be used only in the small-strain context, with
additive split of the strain tensor into the elastic and plastic parts.
The basic equations include the additive decomposition of strain into elastic
and plastic parts,
ε = εe + εp , (32)
the stress strain law
σ = (1 − ω)σ̄ = (1 − ω)D : (ε − εp ), (33)
12
Description DP material
Record Format DruckerPrager num(in) # d(rn) # tAlpha(rn) # E(rn) #
n(rn) # alpha(rn) # alphaPsi(rn) # ht(in) # iys(rn) # lys(rn) #
hm(rn) # kc(rn) # [ yieldtol(rn) #]
- E Young modulus
- n Poisson ratio
- alpha friction coefficient
- alphaPsi dilatancy coefficient
- ht hardening type, 1: linear hardening, 2: exponential
hardening
- iys initial yield stress in shear, τ0
- lys limit yield stress for exponential hardening, τlimit
- hm hardening modulus normalized with E-modulus (!)
- kc κc for the exponential softening law
- yieldtol tolerance of the error in the yield criterion, default
value 1.e-14
- newtonIter maximum number of iterations in λ search,
default value 30
Supported modes 3dMat, PlaneStrain, 3dRotContinuum
Table 8: DP material - summary.
the definition of the yield function in terms of the effective stress,

p
f (σ̄, κ) = αI1 + J2 − τY , (34)
the flow rule p

g (σ̄) = αψ I1 + J2 , (35)
the linear hardening law
τY (κ) = τ0 + Hκ, (36)
where τ0 represents the initial yield stress under pure shear, the damage law
ω(κ) = ωc (1 − e−aκ ), (37)
where ωc is critical damage and a is a positive dimensionless parameter. More

detailed descriptioin of some parameters is in Section 1.4.1.
The dilatancy coefficient αψ controls flow associativeness; if αψ = α, an
associate model is recovered, which overestimates the dilatancy of concrete.
Hence, the dilatancy coefficient is usually chosen smaller, αψ ≤ α, and the
non-associated model is formulated.
13
Description Drucker Prager material with tension cut-off
Record Format DruckerPragerCut num(in) # d(rn) # tAlpha(rn) # E(rn) #
n(rn) # tau0(rn) # alpha(rn) # [alphaPsi(rn) #] [H(rn) #]
[omega crit(rn) #] [a(rn) #] [yieldtol(rn) #] [NewtonIter(in) #]
- E Young modulus
- n Poisson ratio
- tau0 initial yield stress in shear τ0
- alpha friction coefficient
- alphaPsi dilatancy coefficient, equals to alpha by default
- H hardening modulus (can be negative in the case of plas-
tic softening), 0 by default
- omega crit critical damage in damage law (37), 0 by de-
fault
- a exponent in damage law (37), 0 by default
- yieldtol tolerance of the error in the yield criterion, default
value 1.e-14
- newtonIter maximum number of iterations in λ search,
default value 30
Supported modes 1dMat, 3dMat, PlaneStrain
Table 9: Drucker Prager material with tension cut-off - summary.
14
1.4.3 Mises plasticity model with isotropic damage - MisesMat
This model is appropriate for the description of plastic yielding in ductile mate-
rial such as metals, and it can also cover the effects of void growth. The model
uses the Mises yield condition (in terms of the second deviatoric invariant, J2 ),
associated flow rule, linear isotropic hardening driven by the cumulative plastic
strain, and isotropic damage, also driven by the cumulative plastic strain. The
model can be used in the small-strain context, with additive split of the strain
tensor into the elastic and plastic parts, or in the large-strain context, with mul-
tiplicative split of the deformation gradient and with yield condition formulated
in terms of Kirchhoff stress (which is the true Cauchy stress multiplied by the
Jacobian).
Small-strain formulation: The small-strain version of hardening Mises plas-

ticity can be combined with isotropic damage. The basic equations include the
additive decomposition of strain into elastic and plastic parts,
ε = εe + εp , (38)
σ = (1 − ω)σ̄ = (1 − ω)D : (ε − εp ), (39)
the definition of the yield function in terms of the effective stress,

r
3 p
f (s̄, κ) = s̄ : s̄ − σY (κ) = 3J2 (σ̄) − σY (κ), (40)
2
the incremental definition of cumulative plastic strain
κ̇ = kε̇p k, (41)
the linear hardening law

σY (κ) = σ0 + Hκ, (42)
the evolution law for the plastic strain
∂f
ε̇p = λ̇ , (43)
∂s̄
the loading-unloading conditions
λ̇ > 0 f (s̄, κ) ≤ 0 λ̇f (s̄, κ) = 0. (44)
and the damage law

ω(κ) = ωc (1 − e−aκ ), (45)
In the equations above, ε is the strain tensor, εe is the elastic strain tensor, εp
is the plastic strain tensor, D is the elastic stiffness tensor, σ is the nominal
stress tensor, σ̄ is the effective stress tensor, s̄ is the effective deviatoric stress
tensor, σY is the magnitude of stress at yielding under uniaxial tension (or
compression), κ is the cumulated plastic strain, H is the hardening modulus, λ
is the plastic multiplier, ω is the damage variable, ωc is critical damage and a
is a positive dimensionless parameter.
15
Large-strain formulation is based on the introduction of an intermediate
local configuration, with respect to which the elastic response is characterized.
This concept leads to a multiplicative decomposition of deformation gradient
into elastic and plastic parts:
F = F eF p. (46)
The stress-evaluation algorithm can be based on the classical radial return map-
ping; see [20] for more details. Damage is not yet implemented in the large-strain
version of the model.
Description Mises plasticity model with isotropic hardening

Record Format MisesMat (in) # d(rn) # E(rn) # n(rn) # sig0(rn) # H(rn) #
omega crit(rn) #a(rn) #
- E Young’s modulus
- n Poisson’s ratio
- sig0 initial yield stress in uniaxial tension (compression)
- H hardening modulus (can be negative in the case of plas-
tic softening)
- omega crit critical damage in damage law (45)
- a exponent in damage law (45)
Supported modes 1dMat, PlaneStrain, 3dMat, 3dMatF
Table 10: Mises plasticity – summary.

√
VTKxml output can report Mises stress, which equals to 3J2 . When no
hardening/softening exists, Mises stress reaches values up to given uniaxial yield
stress sig0 .
1.4.4 Mises plasticity model with isotropic damage, nonlocal - Mis-

esMatNl, MisesMatGrad
The small-strain version of the model is regularized by the over-nonlocal formu-
lation with damage driven by a combination of local and nonlocal cumulated
plastic strain,
κ̂ = (1 − m)κ + mκ̄, (47)
where m is a dimensionless material parameter (typically m > 1) and κ̄ is the
nonlocal cumulated plastic strain, which is evaluated either using the integral
approach, or using the implicit gradient approach.
Integral nonlocal formulation One possible regularization technique is based

on the integral definition of nonlocal cumulated plastic strain
Z
κ̄(x) = α(x, s)κ(s) ds (48)
V
16
The nonlocal weight function is usually defined as
α0 (kx − sk)
α(x, s) = R (49)
α0 (kx − tk) dt
V
where  2
r2
 1− if r < R

R2

α0 (r) = (50)


0 if r ≥ R

is a nonnegative function, for r < R monotonically decreasing with increasing

distance r = kx − sk, and V denotes the domain occupied by the investigated
material body. The key idea is that the damage evolution at a certain point
depends not only on the cumulated plastic strain at that point, but also on
points at distances smaller than the interaction radius R, considered as a new
material parameter.
Implicit gradient formulation The gradient formulation can be conceived

as the differential counterpart to the integral formulation. The nonlocal cu-
mulated plastic strain is computed from a Helmholtz-type differential equation
κ̄ − l2 ∇2 κ̄ = κ (51)
with homogeneous Neumann boundary condition
∂κ̄
= 0. (52)
∂n
In (51), l is the length scale parameter and ∇ is the Laplace operator.
The model description and parameters are summarized in Tabs. 11 and 12.
Note that the internal length parameter r has the meaning of the radius of
interaction R for the integral version (and thus has the dimension of length)
but for the gradient version it has the meaning of the coefficient l2 multiplying
the Laplacean in (51), and thus has the dimension of length squared.
1.4.5 Rankine plasticity model with isotropic damage and its non-
local formulations - RankMat, RankMatNl, RankMatGrad
This model has a very similar structure to the model described in Section 1.4.3,
but is based on the Rankine yield condition. It is available in the small-strain
version only, and so far exclusively for plane stress analysis. The basic equations
(38)–(39) and (41)–(45) remain valid, and the yield function (40) is redefined
as
f (σ̄, κ) = max σ¯I − σY (κ) (53)
I
where σ¯I are the principal values of the effective stress tensor σ̄. The hardening
law can have either the linear form (42), or the exponential form
σY (κ) = σ0 + ∆σY (1 − exp(−Hκ/∆σY )) , (54)
where H is now the initial plastic modulus and ∆σY is the value of yield stress
increment asymptotically approached as κ → ∞. In damage law (45), parame-
ter ωc is always set to 1. If the plastic hardening is linear, the user can specify
17
Description Nonlocal Mises plasticity with isotropic hardening
Record Format MisesMatNl (in) # d(rn) # E(rn) # n(rn) # sig0(rn) #
H(rn) # omega crit(rn) #a(rn) #r(rn) #m(rn) #[wft(in) #][scal-
ingType(in) #]
- H hardening modulus
- omega crit critical damage
- a exponent in damage law
- r nonlocal interaction radius R from eq. (50)
- m over-nonlocal parameter
- wft selects the type of nonlocal weight function (see Sec-
tion 1.5.7):
1 - default, quartic spline (bell-shaped function with
bounded support)
2 - Gaussian function
3 - exponential function (Green function in 1D)
4 - uniform averaging up to distance R
5 - uniform averaging over one finite element
6 - special function obtained by reducing the 2D expo-
nential function to 1D (by numerical integration)
- scalingType selects the type of scaling of the weight func-

tion (e.g. near a boundary; see Section 1.5.7):
1 - default, standard scaling with integral of weight func-
tion in the denominator
2 - no scaling (the weight function normalized in an infi-
nite body is used even near a boundary)
3 - Borino scaling (local complement)
Supported modes 1dMat, PlaneStrain, 3dMat
Table 11: Nonlocal integral Mises plasticity – summary.
either the exponent a from (45), or the dissipation per unit volume, gf , which
represents the area under the stress-strain diagram (and parameter a is then
determined automatically such that the area under the diagram has the pre-
scribed value). For exponential plastic hardening, the evaluation of a from gf
is not properly implemented and it is better to specify a directly.
The model description and parameters are summarized in Tabs. 13–15. Note
that the default value of parameter m is equal to 1 for the integral model but for
the gradient model it is equal to 2. Also note that the internal length parameter
r has the meaning of the radius of interaction R for the integral version (and thus
has the dimension of length) but for the gradient version it has the meaning of
18
Description Gradient-enhanced Mises plasticity with isotropic damage
Record Format MisesMatGrad (in) # d(rn) # E(rn) # n(rn) # sig0(rn) #
H(rn) # omega crit(rn) #a(rn) #r(rn) #m(rn) #
- H hardening modulus
- omega crit critical damage
- a exponent in damage law
- r internal length scale parameter l2 from eq. (51)
Supported modes 1dMat, PlaneStrain, 3dMat
Table 12: Gradient-enhanced Mises plasticity – summary.
the coefficient l2 multiplying the Laplacean in (51), and thus has the dimension
of length squared.
For the gradient model it is possible to specify parameter negligible damage,
which sets the minimum value of damage that is considered as nonzero. The
approximate solution of Helmholtz equation (51) can lead to very small but
nonzero nonlocal kappa at some points that are actually elastic. If such small
values are positive, they lead to a very small but nonzero damage. If this is
interpreted as ”loading”, the tangent terms are activated, but damage will not
actually grow at such points and the convergence rate is slowed down. It is
better to consider such points as elastic. By default, negligible damage is set to
0, but it is recommended to set it e.g. to 1.e-6.
1.4.6 Perfectly plastic material with Mises yield condition - Steel1

This is an older model, kept here for compatibility with previous versions. It uses
Mises plasticity condition with no hardening and under small strain only. The
model description and parameters are summarized in Tab. 16. All its features
are included in the model described in Section 1.4.3.
1.4.7 Composite plasticity model for masonry - Masonry02

Masonry is a composite material made of bricks and mortar. Nonlinear behavior
of both components should be considered to obtain a realistic model able to de-
scribe cracking, slip, and crushing of the material. The model is based on paper
by Lourenco and Rots [17]. It is formulated on the basis of softening plasticity
for tension, shear, and compression (see fig.(1)). Numerical implementation is
based on modern algorithmic concepts such as implicit integration of the rate
equations and consistent tangent stiffness matrices.
The approach used in this work is based on idea of concentrating all the dam-
age in the relatively weak joints and, if necessary, in potential tension cracks
in the bricks. The joint interface constitutive model should include all impor-
tant damage mechanisms. Here, the concept of interface elements is used. An
19
Description Rankine plasticity with isotropic hardening and damage
Record Format RankMat (in) # d(rn) # E(rn) # n(rn) # plasthardtype(in) #
sig0(rn) # H(rn) # delSigY(rn) # yieldtol(rn) # (gf(rn) # k
a(rn) #)
- plasthardtype type of plastic hardening (0=linear=default,
1=exponential)
- H initial hardening modulus (default value 0.)
- delSigY final increment of yield stress (default value 0.,
needed only if plasthardtype=1)
- yieldtol relative tolerance in the yield condition
- gf dissipation per unit volume
Supported modes PlaneStress
Table 13: Rankine plasticity – summary.
τ
Friction
Cap mode
mode
Intermediate
surface
Tension
mode
σ
Initial surface Residual surface
Figure 1: Composite yield surface model for masonry
interface element allows to incorporate discontinuities in the displacement field

and its behavior is described in terms of a relation between the tractions and
relative displacement across the interface. In the present work, these quanti-
ties will be denoted as σ, generalized stress, and ε, generalized strain. For 2D
configuration, σ = {σ, τ }T and ε = {un , us }T , where σ and τ are the normal
and shear components of the traction interface vector and n and s subscripts
distinguish between normal and shear components of displacement vector. The
elastic response is characterized in terms of elastic constitutive matrix D as
σ = Dε (55)
For a 2D configuration D = diag{kn , ks }. The terms of the elastic stiffness

matrix can be obtained from the properties of both masonry and joints as
Eb Em Gb Gm
kn = ; ks = (56)
tm (Eb − Em ) tm (Gb − Gm )
20
Description Nonlocal Rankine plasticity with isotropic hardening and
damage
Record Format RankMatNl (in) # d(rn) # E(rn) # n(rn) # plasthard-
type(in) # sig0(rn) # H(rn) # delSigY(rn) # yieldtol(rn) #
(gf(rn) # k a(rn) #)r(rn) #m(rn) #[wft(in) #][scalingType(in) #]
1=exponential)
- H initial hardening modulus (default value 0.)
- delSigY final increment of yield stress (default value 0.)
- r internal length scale parameter l2 from eq. (51)
- m over-nonlocal parameter (default value 1.)
- wft selects the type of nonlocal weight function (see Sec-
tion 1.5.7):

bounded support)

tion (e.g. near a boundary; see Section 1.5.7):
Table 14: Nonlocal integral Rankine plasticity – summary.
where Eb and Em are Young’s moduli, Gb and Gm shear moduli for brick and
mortar, and tm is the thickness of joint. One should note, that there is no contact
algorithm assumed between bricks, this means that the overlap of neighboring
units will be visible. On the other hand, the interface model includes a com-
pressive cap, where the compressive inelastic behavior of masonry is lumped.
21
Description Gradient-enhanced Rankine plasticity with isotropic hard-
ening and damage
Record Format RankMatGrad (in) # d(rn) # E(rn) # n(rn) # plasthard-
type(in) # sig0(rn) # H(rn) # delSigY(rn) # yieldtol(rn) #
(gf(rn) # k a(rn) #)r(rn) #m(rn) #negligible damage(rn) #
1=exponential)
- H hardening modulus (default value 0.)
- delSigY final increment of yield stress (default value 0.)
- r internal length scale parameter l from eq. (51)
- m over-nonlocal parameter (default value 2.)
- negligible damage optional parameter (default value 0.)
Table 15: Gradient-enhanced Rankine plasticity – summary.
Description Perfectly plastic material with Mises condition

Record Format Steel1 num(in) # d(rn) # E(rn) # n(rn) # tAlpha(rn) #
Ry(rn) #
- E Young modulus
- n Poisson ratio
- Ry uniaxial yield stress
2dBeamLayer, 3dShellLayer 3dBeam, PlaneStressRot
Table 16: Perfectly plastic material with Mises condition – summary.
Tension mode In the tension mode, the exponential softening law is assumed
(see fig.(3)). The yield function has the following form
f1 (σ, κ1 ) = σ − ft (κ1 ) (57)
where the yield value ft is defined as
!
ft0
ft = ft0 exp − I κ1 (58)
Gf
The ft0 represents tensile strength of joint or interface; and GIf is mode-I frac-
ture energy. For the tension mode, the associated flow hypothesis is assumed.
22
hm
hb
hm
Interface elements (joints)
hm+hb
Continuum elements (brick)
Figure 2: Modeling strategy for masonry
0.2
0.15
0.1
0.05
0
0 0.2 0.4
Figure 3: Tensile behavior of proposed model (ft = 0.2 MPa, GIf = 0.018
N/mm)
Shear mode For the shear mode a Coulomb friction envelope is used. The
yield function has the form
f2 (σ, κ2 ) = |τ | + σ tan φ(κ2 ) − c(κ2 ) (59)
According to [17] the variations of friction angle φ and cohesion c are assumed
as
!
c0
c = c0 exp − II κ2 (60)
Gf

c0 − c
tan φ = tan φ0 + (tan φr − tan φ0 ) (61)
c0
where c0 is initial cohesion of joint, φ0 initial friction angle, φr residual friction

angle, and GIIf fracture energy in mode II failure. A non-associated plastic
23
potential g2 is considered as
g2 = |τ | + σ tan Φ − c (62)
σ=−1.0
σ=−0.5
σ=−0.1
1.5
0.5
0
0 0.2 0.4 0.6 0.8 1
Figure 4: Shear behavior of proposed model for different confinement levels in

MPa (c0 = 0.8 MPa, tan φ0 = 1.0, tan φr = 0.75, and GIIf = 0.05 N/mm)
Coupling of tension/shear modes The tension and Coulomb friction modes

are coupled with isotropic softening. This means that the percentage of soften-
ing in the cohesion is assumed to be the same as on the tensile strength
GIf c0 GII
f ft0
κ̇1 = λ1 + II
λ2 ; κ̇2 = λ1 + λ2 (63)
f
Gf t0 GIf c0
This follows from (58) and (60). However, in the corner region, when both yield
surfaces are activated, such approach will lead to a non-acceptable penalty. For
this reason a quadratic combination is assumed
v !2 v !2
u I
u
u G c u GII f
2 f 0 f t0
κ̇1 = (λ1 ) +
t λ2 ; κ̇2 = t λ1 + (λ2 )2 (64)
GII
f ft0 GI c
f 0
Cap mode For the cap mode, an ellipsoid interface model is used. The yield
condition is assumed as
f3 (σ, κ3 ) = Cnn σ 2 + Css τ 2 + Cn σ − σ̄ 2 (κ3 ) (65)
24
where Cnn , Css , and Cn are material model parameters and σ̄ is yield value,
originally assumed in the following form of hardening/softening law [17]
s
2κ3 κ2
σ̄1 (κ3 ) = σ̄i + (σ̄p − σ̄i ) − 23 ; κ3 ∈ (0, κp )
κp κp
2
κ3 − κp
σ̄2 (κ3 ) = σ̄p + (σ̄m − σ̄p ) ; κ3 ∈ (κp , κm ) (66)
κm − κp

κ3 − κm
σ̄3 (κ3 ) = σ̄r + (σ̄m − σ̄r ) exp m ; κ3 ∈ (κm , ∞)
σ̄m − σ̄r
with m = 2(σ̄m − σ̄p )/(κm − κp ). The hardening/softening law (66) is shown in
fig.(5). Note that the curved diagram is a C 1 continuous σ − κ3 relation. The
energy under the load-displacement diagram can be related to a “compressive
fracture energy”. The original hardening law (66.1) exhibits indefinite slope for
κ3 = 0, which can cause the problems with numerical implementation. This has
been overcomed by replacing this hardening law with parabolic equation given
by
κ3 κ3
σ̄1 (κ3 ) = σ̄i − 2 ∗ (σ̄i − σ̄p ) ∗ + (σ̄i − σ̄p ) (67)
κp κp
An associated flow and strain hardening hypothesis are being considered. This
yields
p
κ̇3 = λ3 (2Cnn σ + Cn ) ∗ (2Cnn σ + Cn ) + (2Css τ ) ∗ (2Css τ ) (68)
σp
σ2
σ1
σm
σi σ3
σr
κp κm κ
Figure 5: Hardening/softening law for cap mode
The model parameters are summarized in Tab. 17. There is one algorithmic
issue, that follows from the model formulation. Since the cap mode harden-
ing/softening is not coupled to hardening/softening of shear and tension modes
the it may happen that when the cap and shear modes are activated, the re-
turn directions become parallel for both surfaces. This should be avoided by
adjusting the input parameters accordingly (one can modify dilatancy angle, for
example).
25
Description Composite plasticity model for masonry
Record Format Masonry02 num(in) # d(rn) # E(rn) # n(rn) # ft0(rn) #
gfi(rn) # gfii(rn) # kn(rn) # ks(rn) # c0(rn) # tanfi0(rn) # tan-
fir(rn) # tanpsi(rn) # si(rn) # sp(rn) # sm(rn) # sr(rn) # kp(rn) #
km(rn) # kr(rn) # cnn(rn) # css(rn) # cn(rn) #
- E Young modulus
- n Poisson ratio
- ft0 tensile strength
- gfi fracture energy for mode I
- gfii fracture energy for mode II
- kn joint elastic property
- ks joint elastic property
- c0 initial cohesion
- tanfi0 initial friction angle
- tanfir residual friction angle
- tanpsi dilatancy
- {si, sp, sm, sr} cap parameters {σ̄i , σ̄p , σ̄m , σ̄r }
- {kp, km,kr} cap parameters {κp , κm , κr }
- cnn,css,cn cap mode parametrs
Supported modes 2dInterface
Table 17: Composite model for masonry - summary.
1.4.8 Nonlinear elasto-plastic material model for concrete plates and

shells - Concrete2
Nonlinear elasto-plastic material model with hardening. Takes into account
uniaxial stress + transverse shear in concrete layers with transverse stirrups.
Can be used only for 2d plates and shells with layered cross section and together
with explicit integration method (stiffness matrix is not provided). The model
description and parameters are summarized in Tab. 18.
1.5 Material models for tensile failure

1.5.1 Nonlinear elasto-plastic material model for concrete plates and
shells - Concrete2
The description can be found is section 1.4.8.
1.5.2 Smeared rotating crack model - Concrete3

Implementation of smeared rotating crack model. Virgin material is modeled as
isotropic linear elastic material (described by Young modulus and Poisson ratio).
The onset of cracking begins, when principal stress reaches tensile strength.
Further behavior is then determined by softening law, governed by principle of
preserving of fracture energy Gf . For large elements, the tension strength can
be artificially reduced to preserve fracture energy. Multiple cracks are allowed.
The elastic unloading and reloading is assumed. In compression regime, this
26
Description Nonlinear elasto-plastic material model for concrete plates
and shells
Record Format Concrete2 num(in) # d(rn) # E(rn) # n(rn) # SCCC(rn) #
SCCT(rn) # EPP(rn) # EPU(rn) # EOPU(rn) #
EOPP(rn) # SHEARTOL(rn) # IS PLASTIC FLOW(in) #
IFAD(in) # STIRR E(rn) # STIRR Ft(rn) #
STIRR A(rn) # STIRR TOL(rn) # STIRR EREF(rn) #
STIRR LAMBDA(rn) #
- E Young modulus
- n Poisson ratio
- SCCC pressure strength
- SCCT tension strength
- EPP threshold effective plastic strain for softening in com-
pression
- EPU ultimate eff. plastic strain
- EOPP threshold volumetric plastic strain for softening in
tension
- EOPU ultimate volumetric plastic strain
- SHEARTOL threshold value of the relative shear defor-
mation (psi**2/eef) at which shear is considered in lay-
ers. For lower relative shear deformations the transverse
shear remains elastic decoupled from bending. default value
SHEARTOL = 0.01
- IS PLASTIC FLOW indicates that plastic flow (not de-
formation theory) is used in pressure
- IFAD State variables will not be updated, otherwise up-
date state variables
- STIRR E Young modulus of stirrups
- STIRR R stirrups uniaxial strength = elastic limit
- STIRR A stirrups area/unit length (beam) or /unit area
(shell)
- STIRR TOL stirrups tolerance of equilibrium in the z
direction (=0 no iteration)
- STIRR EREF stirrups reference strain rate for Peryzna’s
material
- STIRR LAMBDA coefficient for that material (stirrups)
- SHTIRR H isotropic hardening factor for stirrups
Supported modes 3dShellLayer, 2dPlateLayer
Table 18: Nonlinear elasto-plastic material model for concrete - summary.
model correspond to isotropic linear elastic material. The model description

and parameters are summarized in Tab. 19.
27
Description Rotating crack model for concrete
Record Format Concrete3 d(rn) # E(rn) # n(rn) # Gf(rn) # Ft(rn) #
exp soft(in) # tAlpha(rn) #
- E Young modulus
- n Poisson ratio
- Gf fracture energy
- Ft tension strength
- exp soft determines the type of softening (0 = exponential,
1 = linear)
2dBeamLayer, 3dShellLayer
Table 19: Rotating crack model for concrete - summary.
1.5.3 Smeared rotating crack model with transition to scalar damage

- linear softening - RCSD
Implementation of smeared rotating crack model with transition to scalar dam-
age with linear softening law. Improves the classical rotating model (see sec-
tion 1.5.2) by introducing the transition to scalar damage model in later stages
of tension softening.
Traditional smeared-crack models for concrete fracture are known to suffer
by stress locking (meaning here spurious stress transfer across widely opening
cracks), mesh-induced directional bias, and possible instability at late stages
of the loading process. The combined model keeps the anisotropic character
of the rotating crack but it does not transfer spurious stresses across widely
open cracks. The new model with transition to scalar damage (RC-SD) keeps
the anisotropic character of the RCM but it does not transfer spurious stresses
across widely open cracks.
Virgin material is modeled as isotropic linear elastic material (described
by Young modulus and Poisson ratio). The onset of cracking begins, when
principal stress reaches tensile strength. Further behavior is then determined
by linear softening law, governed by principle of preserving of fracture energy
Gf . For large elements, the tension strength can be artificially reduced to
preserve fracture energy. The transition to scalar damage model takes place,
when the softening stress reaches the specified limit. Multiple cracks are allowed.
The elastic unloading and reloading is assumed. In compression regime, this
model correspond to isotropic linear elastic material. The model description
and parameters are summarized in Tab. 20.
1.5.4 Smeared rotating crack model with transition to scalar damage

- exponential softening - RCSDE
Implementation of smeared rotating crack model with transition to scalar dam-
age with exponential softening law. The description and model summary (Tab. 21)
are the same as for the RC-SD model with linear softening law (see section 1.5.3).
28
Description Smeared rotating crack model with transition to scalar
damage - linear softening
Record Format RCSD d(rn) # E(rn) # n(rn) # Gf(rn) # Ft(rn) # sdtransi-
tioncoeff(rn) # tAlpha(rn) #
- E Young modulus
- n Poisson ratio
- sdtransitioncoeff determines the transition from RC to SD
model. Transition takes plase when ratio of current soften-
ing stress to tension strength is less than sdtransitioncoeff
value
Table 20: RC-SD model for concrete - summary.
Description Smeared rotating crack model with transition to scalar

damage - exponential softening
Record Format RCSDE d(rn) # E(rn) # n(rn) # Gf(rn) # Ft(rn) # sdtransi-
tioncoeff(rn) # tAlpha(rn) #
Table 21: RC-SD model for concrete - summary.
1.5.5 Nonlocal smeared rotating crack model with transition to scalar

damage - RCSDNL
Implementation of nonlocal version of smeared rotating crack model with tran-
sition to scalar damage. Improves the classical rotating model (see section 1.5.2)
by introducing the transition to scalar damage model in later stages of tension
softening. The improved RC-SD (see section 1.5.3) is further extended to a
nonlocal formulation, which not only acts as a powerful localization limiter but
also alleviates mesh-induced directional bias. A special type of material insta-
bility arising due to negative shear stiffness terms in the rotating crack model is
resolved by switching to SD mode. A bell shaped nonlocal averaging function
is used.
Virgin material is modeled as isotropic linear elastic material (described
by Young modulus and Poisson ratio). The onset of cracking begins, when
principal stress reaches tensile strength. Further behavior is then determined
by exponential softening law.
The transition to scalar damage model takes place, when the softening stress
reaches the specified limit or when the loss of material stability due to negative
shear stiffness terms that may arise in the standard RCM formulation, which
29
takes place when the ratio of minimal shear coefficient in stiffness to bulk ma-
terial shear modulus reaches the limit.
Multiple cracks are allowed. The elastic unloading and reloading is assumed.
In compression regime, this model correspond to isotropic linear elastic material.
Description Nonlocal smeared rotating crack model with transition to

scalar damage for concrete
Record Format RCSDNL d(rn) # E(rn) # n(rn) # Ft(rn) # sdtransitionco-
eff(rn) # sdtransitioncoeff2(rn) # r(rn) # tAlpha(rn) #
- E Young modulus
- n Poisson ratio
- ef deformation corresponding to fully open crack
- sdtransitioncoeff determines the transition from RC to SD
model. Transition takes place when ratio of current soften-
ing stress to tension strength is less than sdtransitioncoeff
value
- sdtransitioncoeff2 determines the transition from RC to
SD model. Transition takes place when ratio of current
minimal shear stiffness term to virgin shear modulus is less
than sdtransitioncoeff2 value
- r parameter specifying the width of nonlocal averaging
zone
- regionMap map indicating the regions (currently region is
characterized by cross section number) to skip for nonlo-
cal avaraging. The elements and corresponding IP are not
taken into account in nonlocal averaging process if corre-
sponding regionMap value is nonzero.
Table 22: RC-SD-NL model for concrete - summary.
1.5.6 Isotropic damage model for tensile failure - Idm1

This isotropic damage model assumes that the stiffness degradation is isotropic,
i.e., stiffness moduli corresponding to different directions decrease proportionally
and independently of the loading direction. The damaged stiffness tensor is
expressed as D = (1 − ω)D e where ω is a scalar damage variable and D e is the
elastic stiffness tensor. The damage evolution law is postulated in an explicit
form, relating the damage variable ω to the largest previously reached equivalent
strain level, κ.
The equivalent strain, ε̃, is a scalar measure derived from the strain tensor.
The choice of the specific expression for the equivalent strain affects the shape
of the elastic domain in the strain space and plays a similar role to the choice of
30
a yield condition in plasticity. The following definitions of equivalent strain
are currently supported:
• Mazars (1984) definition based on norm of positive part of strain:
v
u 3
uX
ε̃ = t hεI i2 (69)
I=1
where hεI i are positive parts of principal values of the strain tensor ε.
• Definitions derived from the Rankine criterion of maximum principal
stress:
v
u 3
1u X
t hσ̄I i2
ε̃ = (70)
E
I=1
1 3
ε̃ = max σ̄I (71)
E I=1
where σ̄I , I = 1, 2, 3, are the principal values of the effective stress tensor
σ̄ = D e : ε and hσ̄I i are their positive parts.
• Energy norm scaled by Young’s modulus to obtain a strain-like quantity:
1p
ε̃ = ε : De : ε (72)
E
• Modified Mises definition, proposed by de Vree [26]:
s
(k − 1)I1ε 1 (k − 1)2 2 12kJ2ε
ε̃ = + I + (73)
2k(1 − 2ν) 2k (1 − 2ν)2 1ε (1 + ν)2
where
3
X
I1ε = εI
I=1
is the first strain invariant (trace of the strain tensor),

3
1X 2 1 2
J2ε = εI − I1ε
2 6
I=1
is the second deviatoric strain invariant, and k is a model parameter that

corresponds to the ratio between the uniaxial compressive strength fc and
uniaxial tensile strength ft .
• Griffith definition with a solution on inclined elipsoidal inclusion. This
definition handles materials in pure tension and also in compression, where
tensile stresses usually appear on specifically oriented elipsoidal inclusion.
The derivation of Griffith’s criterion is summarized in [14]. In impemen-
tation, first check if Rankine criterion applies
1 3
ε̃ = max σ̄I (74)
E I=1
31
and if not, use Griffith’s solution with ordered principal stresses σ1 > σ3 .
The optional parameter griff n is by default 8 and represents the uniaxial
compression/tensile strength ratio.
1 −(σ1 − σ3 )2
ε̃ = · (75)
E griff n(σ1 + σ3 )
Note that all these definitions are based on the three-dimensional descrip-
tion of strain (and stress). If they are used in a reduced problem, the strain
components that are not explicitly provided by the finite element approximation
are computed from the underlying assumptions and used in the evaluation of
equivalent strain. For instance, in a plane-stress analysis, the out-of-plane com-
ponent of normal strain is calculated from the assumption of zero out-of-plane
normal stress (using standard Hooke’s law).
Since the growth of damage usually leads to softening and may induce lo-
calization of the dissipative process, attention should be paid to proper regu-
larization. The most efficient approach is based on a nonlocal formulation; see
Section 1.5.7. If the model is kept local, the damage law should be adjusted
according to the element size, in the spirit of the crack-band approach. When
done properly, this ensures a correct dissipation of energy in a localized band
of cracking elements, corresponding to the fracture energy of the material. For
various numerical studies, it may be useful to specify the parameters of the
damage law directly, independently of the element size. One should be aware
that in this case the model would exhibit pathological sensitivity to the size of
finite elements if the mesh is changed.
The following damage laws are currently implemented:
• Cohesive crack with exponential softening postulates a relation be-
tween the normal stress σ transmitted by the crack and the crack opening
w in the form
w
σ = ft exp −
wf
Here, ft is the tensile strength and wf is a parameter with the dimension
of length (crack opening), which controls the ductility of the material.
In fact, wf = Gf /ft where Gf is the mode-I fracture energy. In the
context of the crack-band approach, the crack opening w corresponds to
the inelastic (cracking) strain εc multiplied by the effective thickness h of
the crack band. The effective thickness h is estimated by projecting the
finite element onto the direction of the maximum principal strain (and
stress) at the onset of damage. The inelastic strain εc is the difference
between the total strain ε and the elastic strain σ/E. For the damage
model we obtain
σ
εc = ε − = ε − (1 − ω)ε = ωε
E
and thus w = hεc = hωε. Substituting this into the cohesive law and com-
bining with the stress-strain law for the damage model, we get a nonlinear
equation
hωε
(1 − ω)Eε = ft exp −
wf
32
For a given strain ε, the corresponding damage variable ω can be solved
from this equation by Newton iterations. It can be shown that the solution
exists and is unique for every ε ≥ ε0 provided that the element size h
does not exceed the limit size hmax = wf /ε0 . For larger elements, a local
snapback in the stress-strain diagram would occur, which is not admissible.
In terms of the material properties, hmax can be expressed as EGf /ft2 ,
which is related to Irwin’s characteristic length.
The derivation has been performed for monotonic loading and uniaxial
tension. Under general conditions, ε is replaced by the internal variable
κ, which represents the maximum previously reached level of equivalent
strain.
In the list of input variables, the tensile strength ft is not specified directly
but through the corresponding strain at peak stress, ε0 = ft /E, denoted by
keyword e0. Another input parameter is the characteristic crack opening
wf , denoted by keyword wf.
Derivative can be expressed explicitly

ωε
∂ω (ωεf − ε f ) exp εf − ωε0
=−
∂ε εf ε exp ωε − ε0 ε
εf
• Cohesive crack with linear softening is based on the same correspon-

dence between crack opening and inelastic strain, but the cohesive law is
assumed to have a simpler linear form

w
σ = ft 1 −
wf
The relation between damage and strain can then be derived from the
cohesive law and substituing w = hωε

hωε
σ = (1 − ω)Eε = ft 1 − ,
wf
which leads to explicit evaluation of the damage variable
ε0
1−
ω= ε
hε0
1−
wf
and no iteration is needed. Parameter wf , denoted again by keyword
wf, has now the meaning of crack opening at complete failure (zero co-
hesive stress) and is related to fracture energy by a modified formula
wf = 2Gf /ft . The expression for maximum element size, hmax = wf /ε0 ,
remains the same as for cohesive law with exponential softening, but in
terms of the material properties it is now translated as hmax = 2EGf /ft2 .
The derivative with respect to ε yields
∂ω ε0
=
∂ε ε2 1 − hε0
wf
33
• Cohesive crack with bilinear softening is implemented in an approx-
imate fashion and gives for different mesh sizes the same total dissipation
but different shapes of the softening diagram. Instead of properly trans-
forming the crack opening into inelastic strain, the current implementation
deals with a stress-strain diagram adjusted such that the areas marked in
the right part of Fig. 6 are equal to the fracture energies Gf and Gf t
divided by the element size. The third parameter defining the law is the
strain εk at which the softening diagram changes slope. Since this strain
is considered as fixed, the corresponding stress σk depends on the element
size and for small elements gets close to the tensile strength (the diagram
then gets close to linear softening with fracture energy Gf t ).
• Linear softening stress-strain law works directly with strain and does
not make any adjustment for the element size. The specified parameters ε0
and εf , denoted by keywords e0 and ef, have the meaning of (equivalent)
strain at peak stress and at complete failure. The linear relation between
stress and strain on the softening branch is obtained with the damage law
εf ε0
ω= 1−
εf − ε0 ε
Again, to cover general conditions, ε is replaced by κ.

∂ω ε0 εf
= 2
∂ε ε (εf − ε0 )
• Exponential softening stress-strain law also uses two parameters ε0

and εf , denoted by keywords e0 and ef, but leads to a modified dependence
of damage on strain:

ε0 ε − ε0
ω =1− exp −
ε εf − ε0

∂ω 1 1 ε0 − ε
= + 2 ε0 exp
∂ε ε(εf − ε0 ) ε εf − ε0
• Mazars stress-strain law uses three parameters, ε0 , At and Bt , denoted

by keywords e0, At and Bt, and the dependence of damage on strain is
given by
(1 − At )ε0
ω =1− − At exp (Bt (ε − ε0 ))
ε
• Smooth exponential stress-strain law uses two parameters, ε0 and
Md , denoted by keywords e0 and md, and the dependence of damage on
strain is given by
Md !
ε
ω = 1 − exp −
ε0
This leads to a stress-strain curve that immediately deviates from linearity
(has no elastic part) and smoothly changes from hardening to softening,
with tensile strength
−1/Md
ft = Eε0 (eMd )
34
• Extended smooth stress-strain law is a special formulation used by
Grassl and Jirásek [10]. The damage law has a rather complicated form:
 m
1 ε

 1 − exp − if ε ≤ ε1
m ε

p

  
ω= − ε1 (76)
 1 − ε3 exp − h ε
 n i  if ε > ε1
ε ε 1 + ε−ε1



f ε2
The primary model parameters are the uniaxial tensile strength ft , the
strain at peak stress (under uniaxial tension) εp , and additional param-
eters ε1 , ε2 and n, which control the post-peak part of the stress-strain
law. In the input record, they are denoted by keywords ft, ep, e1, e2, nd.
Other parameters that appear in (76) can be derived from the condition
of zero slope of the stress-strain curve at κ = εp and from the conditions
of stress and stiffness continuity at κ = ε1 :
1
m = (77)
ln(Eεp /ft )
ε1
εf = m (78)
(ε1 /εp ) − 1
m
1 ε1
ε3 = ε1 exp − (79)
m εp
Note that parameter damlaw determines which type of damage law should be
used, but the adjustment for element size is done only if parameter wf is specified
for damlaw=0 or damlaw=1. For other values of damlaw, or if parameter ef is
specified instead of wf, the stress-strain curve does not depend on element size
and the model would exhibit pathological sensitivity to the mesh size. These
cases are intended to be used in combination with a nonlocal formulation. An
alternative formulation uses fracture energy to determine fracturing strain.
The model parameters are summarized in Tab. 23. Figure 6 shows three
modes of a softening law with corresponding variables.
Exponential softening Linear softening Bilinear softening
Cohesive stress, σ
Cohesive stress, σ
Cohesive stress, σ
E E E
1 1 1
Equivalent strain, Equivalent strain, Equivalent strain,
Figure 6: Implemented stress-strain diagrams for isotropic damage material.

Fracturing strain εf and crack opening at zero stress wf are interrelated through
effective thickness h of the crack band. Note that exponential softening approach
based on an exponential cohesive law is not exactly equivalent to the approach
based on an exponential softening branch of the stress-strain diagram; see the
detailed discussion of the damage laws.
Description Isotropic damage model for concrete in tension
35
Record Format Idm1 (in) # d(rn) # E(rn) # n(rn) # [tAlpha(rn) #] [equivs-
traintype(in) #] [k(rn) #] [damlaw(in) #] e0(rn) # [wf(rn) #]
[ef(rn) #] [ek(rn) #] [wk(rn) #] [sk(rn) #] [wkwf(rn) #]
[skft(rn) #] [gf(rn) #] [gft(rn) #] [At(rn) #] [Bt(rn) #] [md(rn) #]
[ft(rn) #] [ep(rn) #] [e1(rn) #] [e2(rn) #] [nd(rn) #] [max-
Omega(rn) #] [checkSnapBack(rn) #]
- tAlpha thermal expansion coefficient
- equivstraintype allows to choose from different definitions
of equivalent strain:
0 - default = Mazars, eq. (69)
1 - smooth Rankine, eq. (70)
2 - scaled energy norm, eq. (72)
3 - modified Mises, eq. (73)
4 - standard Rankine, eq. (71)
5 - elastic energy based on positive stress
6 - elastic energy based on positive strain
7 - Griffith criterion eq. (75)
- k ratio between uniaxial compressive and tensile strength,

needed only if equivstraintype=3, default value 1
- damlaw allows to choose from different damage laws:
0 - exponential softening (default) with parameters e0
and wf | ef | gf
1 - linear softening with parameters e0 and wf | ef | gf
2 - bilinear softening with (e0, gf, gft, ek) | (e0, wk, sk,
wf) | (e0, wkwf, skft, wf) | (e0, gf, gft, wk)
3 - Hordijk softening (not implemented yet)
4 - Mazars damage law with parameters At and Bt
5 - smooth stress-strain curve with parameters e0 and md
6 - disable damage (dummy linear elastic material)
7 - extended smooth damage law (76) with parameters
ft, ep, e1, e2, nd
- e0 strain at peak stress (for damage laws 0,1,2,3), limit

elastic strain (for damage law 4), characteristic strain (for
damage law 5)
- wf parameter controling ductility, has the meaning of crack
opening (for damage laws 0 and 1)
- ef parameter controling ductility, has the meaning of
strain (for damage laws 0 and 1)
- ek strain at knee point in bilinear softening type (for dam-
age law 2)
36
- wk crack opening at knee point in bilinear softening type
(for damage law 2)
- sk stress at knee point in bilinear softening type (for dam-
age law 2)
- wkwf ratio of wk/wf < 0, 1 > in bilinear softening type
(for damage law 2)
- skft ratio of sk/ft < 0, 1 > in bilinear softening type (for
damage law 2)
- gf fracture energy (for damage laws 0–2)
- gft total fracture energy (for damage law 2)
- At parameter of Mazars damage law, used only by law 4
- Bt parameter of Mazars damage law, used only by law 4
- md exponent used only by damage law 5, default value 1
- ft tensile strength, used only by damage law 7
- ep strain at peak stress, used only by damage law 7
- e1 parameter used only by damage law 7
- e2 parameter used only by damage law 7
- nd exponent used only by damage law 7
- griff n uniaxial compression/tensile ratio for Griffith’s cri-
terion
- maxOmega maximum damage, used for convergence im-
provement (its value is between 0 and 0.999999 (default),
and it affects only the secant stiffness but not the stress)
- checkSnapBack parameter for snap back checking, 0 no
check, 1 check (default)
Table 23: Isotropic damage model for tensile failure – summary.
1.5.7 Nonlocal isotropic damage model for tensile failure - Idmnl1

Nonlocal version of isotropic damage model from Section 1.5.6. The nonlocal
averaging acts as a powerful localization limiter. In the standard version of
the model, damage is driven by the nonlocal equivalent strain ε̄, defined as a
weighted average of the local equivalent strain:
Z
ε̄(x) = α(x, ξ)ε̃(ξ) dξ
V
In the “undernonlocal” formulation, the damage-driving variable is a combi-

nation of local and nonlocal equivalent strain, mε̄ + (1 − m)ε̃, where m is a
parameter between 0 and 1. (If m > 1, the formulation is called “overnonlocal”;
this case is useful for nonlocal plasticity but not for nonlocal damage.)
Instead of averaging the equivalent strain, one can average the compliance
variable γ, directly related to damage according to the formula γ = ω/(1 − ω).
The weight function α contains a certain parameter with the dimension of
length, which is in general called the characteristic length. Its specific meaning
depends on the type of weight function. The following functions are currently
supported:
37
• Truncated quartic spline, also called the bell-shaped function,
2
s2

α0 (s) = 1 − 2
R
where R is the interaction radius (characteristic length) and s is the dis-
tance between the interacting points. This function is exactly zero for
s ≥ R, i.e., it has a bounded support.
• Gaussian function
s2

α0 (s) = exp − 2
R
which is theoretically nonzero for an arbitrary large s and thus has an
unbounded support. However, in the numerical implementation the value
of α0 is considered as zero for s > 2.5R.
• Exponential function s
α0 (s) = exp −
R
which also has an unbounded support, but is considered as zero for s > 6R.
This function is sometimes called the Green function, because in 1D it
corresponds to the Green function of the Helmholtz-like equation used by
implicit gradient approaches.
• Piecewise constant function

1 if s ≤ R
α0 (s) =
0 if s > R
which corresponds to uniform averaging over a segment, disc or ball of
radius R.
• Function that is constant over the finite element in which point x is lo-
cated, and is zero everywhere else. Of course, this is not a physically
objective definition of nonlocal averaging, since it depends on the dis-
cretization. However, this kind of averaging was proposed in a boundary
layer by Prof. Bažant and was implemented into OOFEM for testing pur-
poses.
• Special function
∞
√ !
s2 + t2
Z
α0 (s) = exp − dt
−∞ R
obtained by reduction of the exponential function from 2D to 1D. The

integral cannot be evaluated in closed form and is computed by OOFEM
numerically. This function can be used in one-dimensional simulations of
a two-dimensional specimen under uniaxial tension; for more details see
[11].
The above functions depend only on the distance s between the interacting
points and are not normalized. If the normalizing condition
Z
α(x, ξ) dξ = 1
V∞
38
is imposed in an infinite body V∞ , it is sufficient to scale α0 by a constant and
set
α0 (kx − ξk)
α(x, ξ) =
Vr∞
where Z
Vr∞ = α0 (kξk) dξ
V∞
Constant Vr∞ can be computed analytically depending on the specific type of

weight function and the number of spatial dimensions in which the analysis is
performed. Since the factor 1/Vr∞ can be incorporated directly in the definition
of α0 , this case is referred to as “no scaling”.
If the body of interest is finite (or even semi-infinite), the averaging integral
can be performed only over the domain filled by the body, and the volume
contributing to the nonlocal average at a point x near the boundary is reduced
as compared to points x far from the boundary or in an infinite body. To make
sure that the normalizing condition
Z
α(x, ξ) dξ = 1
V
holds for the specific domain V , different approaches can be used. The standard
approach defines the nonlocal weight function as
α0 (kx − ξk)
α(x, ξ) =
Vr (x)
where Z
Vr (x) = α0 (kx − ξk) dξ
V
According to the approach suggested by Borino, the weight function is defined
as
α0 (kx − ξk) Vr (x)
α(x, ξ) = + 1− δ(x − ξ)
Vr∞ Vr∞
where δ is the Dirac distribution. One can also say that the nonlocal variable
is evaluated as
Z
1 Vr (x)
ε̄(x) = α0 (kx − ξk)˜
(ξ) dξ + 1 − ε̃(x)
Vr∞ V Vr∞
The term on the right-hand side after the integral is a multiple of the local
variable, and so it can be referred to as the local complement.
In a recent paper [11], special techniques that modify the averaging proce-
dure based on the distance from a physical boundary of the domain or on the
stress state have been considered. The details are explained in [11]. These tech-
niques can be invoked by setting the optional parameter nlVariation to 1, 2 or
3 and specifying additional parameters β and ζ for distance-based averaging, or
β for stress-based averaging.
The model parameters are summarized in Tabs. 24 and 25.
39
Description Nonlocal isotropic damage model for concrete in tension
Record Format Idmnl1 (in) # d(rn) # E(rn) # n(rn) # [tAlpha(rn) #]
[equivstraintype(in) #] [k(rn) #] [damlaw(in) #] e0(rn) #
[ef(rn) #] [At(rn) #] [Bt(rn) #] [md(rn) #] r(rn) # [region-
Map(ia) #] [wft(in) #] [averagingType(in) #] [m(rn) #] [scal-
ingType(in) #] [averagedQuantity(in) #] [nlVariation(in) #]
[beta(rn) #] [zeta(rn) #] [maxOmega(rn) #]
- tAlpha thermal expansion coefficient
- equivstraintype allows to choose from different definitions
of equivalent strain, same as for the local model; see Tab. 23
- k ratio between uniaxial compressive and tensile strength,
needed only if equivstraintype=3, default value 1
- damlaw allows to choose from different damage laws, same
as for the local model; see Tab. 23 (note that parameter wf
cannot be used for the nonlocal model)
- e0 strain at peak stress (for damage laws 0,1,2,3), limit
elastic strain (for damage law 4), characteristic strain (for
damage law 5)
- ef strain parameter controling ductility, has the meaning
of strain (for damage laws 0 and 1), the tangent modulus
just after the peak is Et = −ft /(εf − ε0 )
- At parameter of Mazars damage law, used only by law 4
- Bt parameter of Mazars damage law, used only by law 4
- md exponent, used only by damage law 5, default value 1
- r nonlocal characteristic length R; its meaning depends
on the type of weight function (e.g., interaction radius for
the quartic spline)
- regionMap map indicating the regions (currently region is
characterized by cross section number) to skip for nonlo-
cal avaraging. The elements and corresponding IP are not
taken into account in nonlocal averaging process if corre-
sponding regionMap value is nonzero.
- wft selects the type of nonlocal weight function:
bounded support)
— continued in Tab. 25 —
Table 24: Nonlocal isotropic damage model for tensile failure – summary.
40
Description Nonlocal isotropic damage model for concrete in tension
- averagingType activates a special averaging procedure, de-
fault value 0 does not change anything, value 1 means av-
eraging over one finite element (equivalent to wft=5, but
kept here for compatibility with previous version)
- m multiplier for overnonlocal or undernonlocal formula-
tion, which use m-times the local variable plus (1−m)-times
the nonlocal variable, default value 1
tion (e.g. near a boundary):

- averagedQuantity selects the variable to be averaged, de-

fault value 1 corresponds to equivalent strain, value 2 acti-
vates averaging of compliance variable
- nlVariation activates a special averaging procedure, de-
fault value 0 does not change anything, value 1 means
distance-based averaging (the characteristic length is lin-
early reduced near a physical boundary), value 2 means
stress-based averaging (the averaging is anisotropic and the
characteristic length is affected by the stress), value 3 means
distance-based averaging (the characteristic length is expo-
nentially reduced near a physical boundary)
- beta parameter β, required only for distance-based and
stress-based averaging (i.e., for nlVariation=1, 2 or 3)
- zeta parameter ζ, required only for distance-based aver-
aging (i.e., for nlVariation=1 or 3)
- maxOmega maximum damage, used for convergence im-
provement (its value is between 0 and 0.999999 (default),
and it affects only the secant stiffness but not the stress)
Table 25: Nonlocal isotropic damage model for tensile failure – continued.
41
1.5.8 Anisotropic damage model - Mdm
Local formulation The concept of isotropic damage is appropriate for ma-
terials weakened by voids, but if the physical source of damage is the initiation
and propagation of microcracks, isotropic stiffness degradation can be consid-
ered only as a first rough approximation. More refined damage models take
into account the highly oriented nature of cracking, which is reflected by the
anisotropic character of the damaged stiffness or compliance matrices.
A number of anisotropic damage formulations have been proposed in the
literature. Here we use a model outlined by Jirásek [16], which is based on the
principle of energy equivalence and on the construction of the inverse integrity
tensor by integration of a scalar over all spatial directions. Since the model uses
certain concepts from the microplane theory, it is called the microplane-based
damage model (MDM).
The general structure of the MDM model is schematically shown in Fig. 7
and the basic equations are summarized in Tab. 26. Here, ε and σ are the (nom-
inal) second-order strain and stress tensors with components εij and σij ; e and s
are first-order strain and stress tensors with components ei and si , which char-
acterize the strain and stress on “microplanes” of different orientations given
by a unit vector n with components ni ; ψ is a dimensionless compliance pa-
rameter that is a scalar but can have different values for different directions n;
the symbol δ denotes a virtual quantity; and a sumperimposed tilde denotes
an effective quantity, which is supposed to characterize the state of the intact
material between defects such as microcracks or voids.
Table 26: Basic equations of microplane-based anisotropic damage model
ẽ = ε̃ · n sT = ψs s=σ·n
Z Z
3 3
σ̃ : δε̃ = sT · δẽ dΩ δs · e = ds · ẽ T
δσ : ε = δs · e dΩ
2π Ω 2π Ω
Z Z
3 3
σ̃ = (sT ⊗ n)sym dΩ e = ψẽ ε= (e ⊗ n)sym dΩ
2π Ω 2π Ω
εij ~
εij
elastic ~
σij σij
compliance
constraint
constraint
kinematic
static
ei ~
ei ~
si si
ψ ψ
Figure 7: Structure of microplane-based anisotropic damage model
Combining the basic equations, it is possible to show that the components
42
of the damaged material compliance tensor are given by
e
Cijkl = Mpqij Mrskl Cpqrs (80)
e
where Cpqrs are the components of the elastic material compliance tensor,
1
Mijkl = 4 (ψik δjl + ψil δjk + ψjk δil + ψjl δik ) (81)
are the components of the so-called damage effect tensor, and

Z
3
ψij = ψ ni nj dΩ (82)
2π Ω
are the components of the second-order inverse integrity tensor. The integra-
tion domain Ω is the unit hemisphere. In practice, the integral over the unit
hemisphere is evaluated by summing the contribution from a finite number of
directions, according to one of the numerical integration schemes that are used
by microplane models.
The scalar variable ψ characterizes the relative compliance in the direction
given by the vector n. If ψ is the same in all directions, the inverse integrity
tensor evaluated from (82) is equal to the unit second-order tensor (Kronecker
delta) multiplied by ψ, the damage effect tensor evaluated from (81) is equal
to the symmetric fourth-order unit tensor multiplied by ψ, and the damaged
material compliance tensor evaluated from (80) is the elastic compliance tensor
multiplied by ψ 2 . The factor multiplying the elastic compliance√tensor in the
isotropic damage model is 1/(1 − ω), and so ψ corresponds to 1/ 1 − ω. In the
initial undamaged state, ψ = 1 in all directions. The evolution of ψ is governed
by the history of the projected strain components. In the simplest case, ψ is
driven by the normal strain eN = εij ni nj . Analogy with the isotropic damage
model leads to the damage law
ψ = f (κ) (83)
and loading-unloading conditions
g(eN , κ) ≡ eN − κ ≤ 0, κ̇ ≥ 0, κ̇g(eN , κ) = 0 (84)
in which κ is a history variable that represents the maximum level of normal

strain in the given direction ever reached in the previous history of the mate-
rial. An appropriate modification of the exponential softening law leads to the
damage law 
 1r if κ ≤ e0
f (κ) = κ

κ−e0
(85)
e0 exp ef −e0 if κ > e0

where e0 is a parameter controlling the elastic limit, and ef > e0 is another

parameter controlling ductility. Note that softening in a limited number of di-
rections does not necessarily lead to softening on the macroscopic level, because
the response in the other directions remains elastic. Therefore, e0 corresponds
to the elastic limit but not to the state at peak stress.
If the MDM model is used in its basic form described above, the compres-
sive strength turns out to depend on the Poisson ratio and, in applications to
43
concrete, its value is too low compared to the tensile strength. The model is de-
signed primarily for tensile-dominated failure, so the low compressive strength
is not considered as a major drawback. Still, it is desirable to introduce a
modification that would prevent spurious compressive failure in problems where
moderate compressive stresses appear. The desired effect is achieved by redefin-
ing the projected strain eN as
εij ni nj
eN = m (86)
1− σkk
Ee0
where m is a nonnegative parameter that controls the sensitivity to the mean
stress, σkk is the trace of the stress tensor, and the normalizing factor Ee0 is
introduced in order to render the parameter m dimensionless. Under compres-
sive stress states (characterized by σkk < 0), the denominator in (86) is larger
than 1, and the projected strain is reduced, which also leads to a reduction of
damage. A typical recommended value of parameter m is 0.05.
Nonlocal formulation Nonlocal formulation of the MDM model is based

on the averaging of the inverse integrity tensor. This roughly corresponds to
the nonlocal isotropic damage model with averaging of the compliance variable
γ = ω/(1 − ω), which does not cause any spurious locking effects. In equation
(81) for the evaluation of the damage effect tensor, the inverse integrity tensor
is replaced by its weighted average with components
Z
ψ̄ij (x) = α(x, ξ)ψij (ξ)dξ (87)
V
By fitting a wide range of numerical results, it has been found that the
parameters of the nonlocal MDM model can be estimated from the measurable
material properties using the formulas
EGf
λf = (88)
Rft2
λf
λ = (89)
1.47 − 0.0014λf
ft
e0 = (90)
(1 − m)E(1.56 + 0.006λ)
ef = e0 [1 + (1 − m)λ] (91)
where E is Young’s modulus, Gf is the fracture energy, ft is the uniaxial tensile

strength, m is the compressive correction factor, typically chosen as m = 0.05,
and R is the radius of nonlocal interaction reflecting the internal length of the
material.
Input Record The model description and parameters are summarized in

Tab. 27.
44
Description MDM Anisotropic damage model
Common parameters
Record Format Mdm d(rn) # nmp(ins) # talpha(rn) # parmd(rn) # non-
loc(in) # formulation(in) # mode(in) #
Parameters -num material model number
- D material density
- nmp number of microplanes used for hemisphere integra-
tion, supported values are 21,28, and 61
- talpha thermal dillatation coeff
- parmd
- nonloc
- formulation
- mode
Nonlocal variant I
Additional params r(rn) # efp(rn) # ep(rn) #
-r nonlocal interaction radius
-efp εf p is a model parameter that controls the post-peak
slope εf p =εf − ε0 , where εf is strain at zero stress level.
-ep max effective strain at peak ε0
Nonlocal variant II
Additional params r(rn) # gf(rn) # ft(rn) #
-r nonlocal intraction radius
-gf fracture energy
-ft tensile strength
Local variant I
Additional params efp(rn) # ep(rn) #
-efp εf p is a model parameter that controls the post-peak
slope εf p =εf − ε0 , where εf is strain at zero stress level.
Local variant II
Additional params gf(rn) # ep(rn) #
-gf fracture energy
Supported modes 3dMat, PlaneStress
Table 27: MDM model - summary.
45
1.5.9 Isotropic damage model for interfaces
The model provides an interface law which can be used to describe a damageable
interface between two materials (e.g. between steel reinforcement and concrete).
The law is formulated in terms of the traction vector and the displacement jump
vector. The basic response is elastic, with stiffness kn in the normal direction
and ks in the tangential direction. Similar to other isotropic damage models,
this model assumes that the stiffness degradation is isotropic, i.e., both stiffness
moduli decrease proportionally and independently of the loading direction. The
damaged stiffnesses are kn×(1 − ω) and ks×(1 − ω) where ω is a scalar damage
variable. The damage evolution law is postulated in an explicit form, relating
the damage variable ω to the largest previously reached equivalent “strain” level,
κ.
The equivalent “strain”, ε̃, is a scalar measure of the displacement jump
vector. The choice of the specific expression for the equivalent strain affects the
shape of the elastic domain in the strain space and plays a similar role to the
choice of a yield condition in plasticity. Currently, in the present implementa-
tion, the equivalent strain is given by
p
ε̃ = hwn i2 + βws2
where hwn i is the positive part of the normal displacement jump (opening of
the interface) and ws is the norm of the tangential part of displacement jump
(sliding of the interface). Parameter β is optional and its default value is 0,
in which case damage depends on the opening only (not on the sliding). The
dependence of damage ω on maximum equivalent strain κ is described by the
following damage law which corresponds to exponential softening:

 0 for κ ≤ ε0
ω= ε0 ft (κ − ε0 )
 1− exp − for κ > ε0
κ Gf
Here, ε0 = ft /kn is the value of equivalent strain at the onset of damage. Note
that if the interface is subjected to shear traction only (with zero or negative
normal traction), the propagation of √ damage starts when the magnitude of the
sliding displacement is |ws | = ε0 / β, i.e., when the magnitude of the shear
traction is equal to
ks ε0 ks
fs = √ = ft √
β kn β
So the ratio between the √ shear strength and tensile strength of the interface,
fs /ft , is equal to ks /kn β.
The model parameters are summarized in Tab. 28.
1.5.10 Isotropic damage model for interfaces using tabulated data

for damage
The model provides an interface law which can be used to describe a damageable
interface between two materials (e.g. between steel reinforcement and concrete).
The law is formulated in terms of the traction vector and the displacement jump
vector. The basic response is elastic, with stiffness kn in the normal direction
and ks in the tangential direction. Similar to other isotropic damage models,
46
Record Format isointrfdm01 kn(rn) # ks(rn) # ft(rn) # gf(rn) # [ max-
omega(rn) #] talpha(rn) # d(rn) #
Parameters - d material density
- kn elastic stifness in normal direction
- ks elastic stifness in tangential direction
- ft tensile strength
- gf fracture energy
- [maxomega] maximum damage, used for convergence improve-
ment (its value is between 0 and 0.999999 (default), and it affects
only the secant stiffness but not the stress)
- [beta] parameter controlling the effect of sliding part of dis-
placement jump on equivalent strain, default value 0
Supported modes 2dInterface, 3dInterface
Features
Table 28: Isotropic damage model for interface elements – summary.
this model assumes that the stiffness degradation is isotropic, i.e., both stiffness
moduli decrease proportionally and independently of the loading direction. The
damaged stiffnesses are kn×(1 − ω) and ks×(1 − ω) where ω is a scalar damage
variable.
The equivalent “strain”, ε̃, is a scalar measure derived from the displacement
jump vector. The choice of the specific expression for the equivalent strain
affects the shape of the elastic domain in the strain space and plays a similar
role to the choice of a yield condition in plasticity. Currently, in the present
implementation, ε̃ is equal to the positive part of the normal displacement jump
(opening of the interface).
The damage evolution law is postulated in a separate file that should have
the following format. Each line should contain one strain, damage pair separated
by a whitespace character. The exception to this is the first line which should
contain a single integer stating how many strain, damage pairs that the file will
contain. The strains given in the file is defined as the equivalent strain minus
the limit of elastic deformation. To find the damage for arbitrary strains linear
interpolation between the tabulated values is used. If a strain larger than one
in the given table is achieved the respective damage for the largest tabulated
strain will be used. Both the strains and damages must be given in a strictly
increasing order.
1.6 Material models specific to concrete

1.6.1 Mazars damage model for concrete - MazarsModel
This isotropic damage model assumes that the stiffness degradation is isotropic,
i.e., stiffness moduli corresponding to different directions decrease proportionally
and independently of direction of loading. It introduces two damage parameters
ωt and ωc that are computed from the same equivalent strain using two different
damage functions gt and gc . The gt is identified from the uniaxial tension tests,
while gc from compressive test. The damage parameter for general stress states
47
Record Format isointrfdm02 kn(rn) # ks(rn) # ft(rn) # tablename(rn) # [ max-
omega(rn) #] talpha(rn) # d(rn) #
- kn elastic stifness in normal direction
- ks elastic stifness in tangential direction
- tablename file name of the table with the strain damage pairs
- maxomega maximum damage, used for convergence improve-
ment (its value is between 0 and 0.999999 (default), and it affects
only the secant stiffness but not the stress)
Features
Table 29: Isotropic damage model for interface elements using tabulated data
for damage – summary.
ω is obtained as a linear combination of ωt and ωc : ω = αt gt + αc gc , where the

coefficients αt and αc take into account the character of the stress state. The
damaged stiffness tensor is expressed as D = (1 − ω)D e . Damage evolution law
is postulated in an explicit form, relating damage parameter and scalar measure
of largest reached strain level in material, taking into account the principle of
preserving of fracture energy Gf . The equivalent strain, i.e., a scalar measure
of the strain level is defined as norm from positive principal strains. The model
description and parameters are summarized in Tab. 30.
1.6.2 Nonlocal Mazars damage model for concrete - MazarsModelnl

The nonlocal variant of Mazars damage model for concrete. Model based on
nonlocal averaging of equivalent strain. The nonlocal averaging acts as a pow-
erful localization limiter. The bell-shaped nonlocal averaging function is used.
1.6.3 CebFip78 model for concrete creep with aging - CebFip78

Implementation of aging viscoelastic model for concrete creep according to the
CEB-FIP Model Code. The model parameters are summarized in Tab. 32.
1.6.4 Double-power law model for concrete creep with aging - Dou-
blePowerLaw
Implementation of aging viscoelastic model for concrete creep with compliance
function given by the double-power law. The model parameters are summarized
in Tab. 33.
1.6.5 Eurocode 2 model for concrete creep and shrinkage - EC2CreepMat

Implementation of aging viscoelastic model for concrete creep according to
Eurocode 2 for concrete structures The model parameters are summarized in
Tab. 34.
48
Description Mazars damage model for concrete
Record Format MazarsModel d(rn) # E(rn) # n(rn) # e0(rn) # ac(rn) #
[bc(rn) #] [beta(rn) #] at(rn) # [ bt(rn) #] [hreft(rn) #]
[hrefc(rn) #] [version(in) #] [tAlpha(rn) #] [equivstrain-
type(in) #] [maxOmega(rn) #]
- E Young modulus
- n Poisson ratio
- e0 max effective strain at peak
- ac,bc material parameters related to the shape of uniaxial
compression curve (A sample set used by Saouridis is Ac =
1.34, Bc = 2537
- beta coefficient reducing the effect of damage under re-
sponse under shear. Default value set to 1.06
- at, [bt] material parameters related to the shape of uniaxial
tension curve. Meaning dependent on version parameter.
- hreft, hrefc reference characteristic lengths for tension and
compression. The material parameters are specified for ele-
ment with these characteristic lengths. The current element
then will have the same COD (Crack Opening Displace-
ment) as reference one.
- version Model variant. if 0 specified, the original form
gt = 1.0 − (1.0 − At ) ∗ ε0 /κ − At ∗ exp(−Bt ∗ (κ − ε0 ));
of tension damage evolution law is used, if equal 1, the
modified law used which asymptotically tends to zero gt =
1.0 − (ε0 /κ) ∗ exp((ε0 − κ)/εf )
- equivstraintype see Tab. 23
- maxOmega limit maximum damage, use for convergency
improvement
Table 30: Mazars damage model – summary.
According to EC2, the compliance function is defined using the creep coef-
ficient ϕ as
1 ϕ(t, t0 )
J(t, t0 ) = 0
+ (92)
E(t ) 1.05Ecm
where Ecm is the mean elastic modulus at the age of 28 days and E(t0 ) is the
elastic modulus at the age of loading, t0 .
Current implementation supports only linear creep which is valid only for
stresses below 0.45 of the characteristic compressive strength at the time of
loading.
The elastic modulus at age t (in days) is defined as
h p i0.3
E(t) = exp s 1 − 28/t Ecm (93)
49
Description Nonlocal Mazars damage model for concrete
Record Format MazarsModelnl r(rn) # E(rn) # n(rn) # e0(rn) # ac(rn) #
bc(rn) # beta(rn) # version(in) # at(rn) # [ bt(rn) #] r(rn) #
tAlpha(rn) #
- E Young modulus
- n Poisson ratio
- maxOmega limit maximum damage, use for convergency
improvement
- version Model variant. if 0 specified, the original form
gt = 1.0 − (1.0 − At ) ∗ ε0 /κ − At ∗ exp(−Bt ∗ (κ − ε0 ));
of tension damage evolution law is used, if equal 1, the
modified law used which asymptotically tends to zero gt =
1.0 − (ε0 /κ) ∗ exp((ε0 − κ)/At )
- ac,bc material parameters related to the shape of uniaxial
compression curve (A sample set used by Saouridis is Ac =
1.34, Bc = 2537
- at, [bt] material parameters related to the shape of uniaxial
tension curve. Meaning dependent on version parameter.
- beta coefficient reducing the effect of damage under re-
sponse under shear. Default value set to 1.06
- r parameter specifying the width of nonlocal averaging
zone
Table 31: Nonlocal Mazars damage model – summary.
Description CebFip78 model for concrete creep with aging

Record Format CebFip78 n(rn) # relMatAge(rn) # E28(rn) # fibf(rn) #
kap a(rn) # kap c(rn) # kap tt(rn) # u(rn) #
- E28 Young modulus at age of 28 days [MPa]
- n Poisson ratio
- fibf basic creep coefficient
- kap a coefficient of hydrometric conditions
- kap c coefficient of type of cement
- kap tt coeficient of temperature effects
- u surface imposed to environment [mm2 ]; temporary here;
should be in crosssection level
- relmatage relative material age
Supported modes 3dMat, PlaneStress, PlaneStrain, 1dMat, 2dPlate-
Layer,2dBeamLayer, 3dShellLayer
Table 32: CebFip78 material model – summary.
where s is a cement-type dependent constant (0.2 for class R, 0.25 for class N
50
Description Double-power law model for concrete creep with aging
Record Format DoublePowerLaw n(rn) # relMatAge(rn) # E28(rn) #
fi1(rn) # m(rn) # n(rn) # alpha(rn) #
- E28 Young modulus at age of 28 days [MPa]
- n Poisson ratio
- fibf basic creep coefficient
- m coefficient
- n coefficient
- alpha coeficient
- relmatage relative material age
Table 33: Double-power law model – summary.
and finally 0.38 for type S), and the mean secant elastic modulus at 28 days can
be estimated from the mean compressive strength
0.3
Ecm = 22 (0.1fcm,28 ) (94)
where fcm,28 is in MPa and the resulting modulus is in GPa.

The creep coefficient is given by the expression from Annex B of the standard.
0.3
t − t0

16.8 1
ϕ(t, t0 ) = ϕRH × × (95)
fcm 0.1 + t00.2 βH + t − t0T
with

1 − henv
ϕRH = 1 + α1 × α2 (96)
0.1 (h0 )1/3

18
βH = 1.5 1 + (1.2henv ) h0 + 250 α3 ≤ 1500 α3 (97)
and
α1 = α2 = α3 = 1 for fcm ≤ 35 MPa (98)
else
0.7
35
α1 = (99)
fcm
0.2
35
α2 = (100)
fcm
0.5
35
α3 = (101)
fcm
and t0T is the temperature-adusted age according to B.9 in the code. It is

computed automatically if string temperatureDependent appears in the input
record.
51
The shrinkage deformation is additively split into two parts: drying shrink-
age εsh,d and autogenous shrinkage εsh,a . Drying shrinkage strain at time t is
computed from
t − t0
εsh,d = kh εsh,d,0 (102)
t − t0 + 0.04 h3/2
where t0 is the duration of curing, kh is a thickness-dependent parameter and
εsh,d,0 = 1.3175 × 10−6 [(220 + 110αds1 ) exp (−0.1αds2 fcm )] 1 − h3env (103)

with αds1 = 3 for cement class S, 4 for class N, and 6 for class R, and αds2 = 0.13
for cement class S, 0.12 for class N, and 0.11 for class R.
Autogenous shrinkage strain can be comptued as
h √ i
εsh,a = 2.5 (fcm − 18) 1 − exp −0.2 t × 10−6 (104)
1.6.6 B3 and MPS models for concrete creep with aging

Model B3 is an aging viscoelastic model for concrete creep and shrinkage, de-
veloped by Prof. Bažant and coworkers. In OOFEM it is implemented in three
different ways.
The first version, “B3mat”, is kept in OOFEM for compatibility. It is based
on an aging Maxwell chain. The moduli of individual units in the chain are
evaluated in each step using the least-squares method.
The second, more recent version, is referred to as “B3solidmat”. Depend-
ing on the specified input it exploits either a non-aging Kelvin chain combined
with the solidification theory, or an aging Kelvin chain. It is extended to the
so-called Microprestress-Solidification theory (MPS), which in this implementa-
tion takes into account only the effects of variable humidity on creep; the effects
of temperature on creep are not considered. The underlying rheological chain
consists of four serially coupled components. The solidifying Kelvin chain repre-
sents short-term creep; it is serially coupled with a non-aging elastic spring that
reflects instantaneous deformation. Long-term creep is captured by an aging
dashpot with viscosity dependent on the microprestress, the evolution of which
is affected by changes of humidity. The last unit describes volumetric deforma-
tions (shrinkage and thermal strains). Drying creep is incorporated either by
the “averaged cross-sectional approach”, or by the “point approach”.
The latest version is denoted as “MPS” and is based on the microprestress-
solidification theory [28] [29] [30]. The rheological model consists of the same
four components as in “B3solidmat”, but now the implemented exponential al-
gorithm is designed especially for the solidifying Kelvin chain, which is a special
case of an aging Kelvin chain. This model takes into account both humidity and
temperature effects on creep. Drying creep is incorporated exclusively by the so-
called “point approach”. The model can operate in four modes, controlled by the
keyword CoupledAnalysisT ype. The first mode (CoupledAnalysisT ype = 0)
solves only the basic creep and runs as a single problem, while the remain-
ing three modes need to be run as a staggered problem with humidity and/or
temperature analysis preceding the mechanical problem; both humidity and
temperature fields are read when CoupledAnalysisT ype = 1, only the field of
52
relative humidity is taken into account when CoupledAnalysisT ype = 2 and
finally, only temperature when CoupledAnalysisT ype = 3.
The basic creep is in the microprestress-solidification theory influenced by
the same four parameters q1 - q4 as in the model B3. Values of these parameters
can be estimated from the composition of concrete mixture and its compressive
strength using the following empirical formulae:
−0.5
q1 = 126.77f¯c [10−6 /MPa] (105)
−0.9
q2 = 185.4c f¯c
0.5
[10 −6
/MPa] (106)
4 −6
q3 = 0.29 (w/c) q2 [10 /MPa] (107)
−0.7 −6
q4 = 20.3 (a/c) [10 /MPa] (108)
Here, f¯c is the average compressive cylinder strength at age of 28 days [MPa],
a, w and c is the weight of aggregates, water and cement per unit volume of
concrete [kg/m3 ].
The non-aging spring stiffness represents the asymptotic modulus of the
material; it is equal to 1/q1 . The solidifying Kelvin chain is composed of M
Kelvin units with fixed retardation times τµ , µ = 1, 2, . . . , M , which form a
geometric progression with quotient 10. The lowest retardation time τ1 is equal
to 0.3 begof timeof interest, the highest retardation time τM is bigger than
0.5 endof timeof interest. The chain also contains a spring with stiffness E0∞ (a
special case of Kelvin unit with zero retardation time). Moduli Eµ∞ of individual
Kelvin units are determined such that the chain provides a good approximation
of the non-aging micro-compliance function of the solidifying constituent, Φ(t −
n
t0 ) = q2 ln 1 + ((t − t0 ) /λ0 ) , where λ0 = 1 day and n = 0.1. The technique
based on the continuous retardation spectrum leads to the following formulae:
0.1
1 q2 τ˜0 2τ1
= q2 ln (1 + τ˜0 ) − where τ˜0 = √ (109)
E0∞ 10 (1 + τ˜0 ) 10
1 q2 τ̃µ (0.9 + τ̃µ ) 0.1
= (ln 10) 2 where τ̃µ = (2τµ ) , µ = 1, 2, .(110)
..M
Eµ∞ 10 (1 + τ̃µ )
Viscosities ηµ∞ of individual Kelvin units are obtained from simple relation ηµ∞ =
τµ /Eµ∞ . A higher accuracy is reached if all retardation times are in the end
multiplied by the factor 1.35 and the last modulus EM is divided by 1.2.
The actual viscosities ηµ and stiffnesses Eµ of the solidifying chain change
in time according to ηµ (t) = v(t)ηµ∞ and Eµ (t) = v(t)Eµ∞ , where
1
v(t) = q3 λ0 m
(111)
q2 + t
is the volume growth function, and exponent m = 0.5. In the case of vari-
able temperature or humidity, the actual age of concrete t is replaced by the
equivalent time te , which is obtained by integrating (114).
Evolution of viscosity of the aging dashpot is governed by the differential
equation
1 ḣ
p/(p−1) ψS
η̇ + T − κT kT Ṫ (µS η) = (112)

µS T0 h q4
53
where h is the relative pore humidity, T is the absolute temperature [K], T0 =
298 K is the room temperature, and parameter p = 2. Parameter kT is different
for monotonically increasing and for cyclic temperature, and is defined as

kT m if T = Tmax and Ṫ > 0
κT = (113)
kT c if T < Tmax or Ṫ ≤ 0
in which kT m [-] and kT c [-] are new parameters and Tmax is the maximum
temperature attained in the previous history of the material point.
Equation (112) differs from the one presented in the original work; it re-
places the differential equation for microprestress, which is not used here. The
evolution of viscosity can be captured directly, without the need for micro-
prestress. What matters is only the relative humidity and temperature and
their rates. Parameters c0 and k1 of the original MPS theory are replaced
by µS = c0 T0p−1 k1p−1 q4 (p − 1)p . The initial value of viscosity is defined as
η(t0 ) = t0 /q4 , where t0 is age of concrete at the onset drying or when the tem-
perature starts changing, in the present implementation it is set relM atAge
which corresponds to the material age when the material is cast.
As mentioned above, under variable humidity and temperature conditions
the physical time t in function v(t) describing evolution of the solidified vol-
ume is replaced by the equivalent time te . In a similar spirit, t is replaced by
the solidification time ts in the equation describing creep of the solidifying con-
stituent, and by the reduced time tr in equation dεf /dtr = σ/η(t) relating the
flow strain rate to the stress. Factors transforming the physical time t into te ,
tr and ts are defined as follows:
dte
= ψe (t) = βeT (T (t)) βeh (h(t)) (114)
dt
dtr
= ψr (t) = βrT (T (t)) βrh (h(t)) (115)
dt
dts
= ψs (t) = βsT (T (t)) βsh (h(t)) (116)
dt
Functions describing the influence of temperature have the form

Qe 1 1
βeT (T ) = exp − (117)
R T0 T

Qr 1 1
βrT (T ) = exp − (118)
R T0 T

Qs 1 1
βsT (T ) = exp − (119)
R T0 T
motivated by the rate process theory. R is the universal gas constant and Qe , Qr ,
Qs are activation energies for hydration, viscous processes and microprestress
relaxation, respectively. Only the ratios Qe /R, Qr /R and Qs /R have to be
specified. Functions describing the influence of humidity have the form
1
βeh (h) = 4 (120)
1 + [αe (1 − h)]
βrh (h) = αr + (1 − αr ) h2 (121)
2
βsh (h) = αs + (1 − αs ) h (122)
54
where αe , αr and αs are parameters.
At sealed conditions (or CoupledAnalysisT ype = 2) the auxiliary coeffi-
cients βeT = βrT = βsT = 1 while at room temperature (or with CoupledAnalysisT ype =
3) factors βe,h = βr,h = βs,h = 1.
It turned out that both the size effect on drying creep as well as its delay
behind drying shrinkage can be addressed through parameter p in the governing
equation (112). In the experiments, the average (cross-sectional) drying creep is
decreasing with specimen size. Unfortunately, for the standard value p = 2, the
MPS model exhibits the opposite trend. For p = ∞ the size effect disappears,
and for p < 0 it corresponds to the experiments. It should be noted that for
negative or infinite values of p the underlying theory loses its original physical
background. If the experimental data of drying creep measured on different
sizes are missing, the exponent p = ∞ can be taken as a realistic estimate.
When parameter p is changed from its recommended value p = 2 it is ad-
vantageous to rewrite the governing differential equation to the following form

k3 ḣ
p̃ ψS
η̇f + T − κT Ṫ ηf = (123)
T0 h q4
with newly introduced parameters
p̃ = p/(p − 1) (124)
1
k3 = µSp−1 (125)
With the “standard value” p = 2 (reverted size effect) the new parameter
p̃ is also 2, for p < 0 (correct size effect) p̃ < 1 and finally with p = ∞ the
first parameter p̃ = 1 and the second parameter k3 becomes dimensionless. The
other advantage of p̃ = 1 is that the governing differential equation becomes
linear and can be solved directly.
The rate of thermal strain is expressed as ε̇T = αT Ṫ and the rate of drying
shrinkage strain as ε̇sh = ksh ḣ, where both αT and ksh are assumed to be
constant in time and independent of temperature and humidity.
There are two options to simulate the autogenous shrinkage in the MPS
material model. The first one is according to the B4 model
" w/c/0.38 #−4.5
τau
εsh,au,B4 (te ) = ε∞
sh,au,B4 1+ (126)
te
and the second one is proposed in the Model Code 2010

√
εsh,au,f ib (te ) = ε∞
sh,au,f ib 1 − exp −0.2 te (127)
For the normally hardening cement, the ultimate value of the autogenous
shrinkage can be estimated from the composition using the empirical formula
of the B4 model
−0.75 −3.5
∞ −6 a/c w/c
εsh,au,B4 = −210 × 10 (128)
6 0.38
55
Similarly, in the Model Code, the ultimate shrinkage strain can be estimated
from the mean concrete strength at the age of 28 days and from the cement grade
2.5
0.1fcm
ε∞
sh,au,f ib = −αas × 10−6 (129)
6 + 0.1fcm
where αas is 600 for cement grades 42.5 R and 52.5 R and N, 700 for 32.5 R
and 42.5 N and 800 for 32.5 N.
The model description and parameters are summarized in Tab. 35 for “B3mat”,
in Tab. 36 for “B3solidmat”, and in Tab. 37 for “MPS”. Since some model pa-
rameters are determined from the composition and strength using empirical
formulae, it is necessary to use the specified units (e.g. compressive strength
always in MPa, irrespectively of the units used in the simulation for stress). For
“B3mat” and “B3solidmat” it is strictly required to use the specified units in
the material input record (stress always in MPa, time in days etc.). The “MPS”
model is almost unit-independent, except for f¯c in MPa and c in kg/m3 , which
are used in empirical formulae.
For illustration, sample input records for the material considered in Example
3.1 of the creep book by Bažant and Jirásek is presented. The concrete mix is
composed of 170 kg/m3 of water, 450 kg/m3 of type-I cement and 1800 kg/m3
of aggregates, which corresponds to ratios w/c = 0.3778 and a/c = 4. The
compressive strength is f¯c = 45.4 MPa. The concrete slab of thickness 200 mm
is cured in air with initial protection against drying until the age of 7 days.
Subsequently, the slab is exposed to an environment with relative humidity of
70%. The following input record can be used for the first version of the model
(B3mat):
B3mat 1 n 0.2 d 0. talpha 1.2e-5 relMatAge 28. fc 45.4 cc 450.
w/c 0.3778 a/c 4. t0 7. timefactor 1. alpha1 1. alpha2 1.2 ks 1.
hum 0.7 vs 0.1 shmode 1
Parameter α1 = 1 corresponds to type-I cement, parameter α2 = 1.2 to
curing in air, parameter ks = 1 to an infinite slab. The volume-to-surface ratio is
in this case equal to one half of the slab thickness and must be specified in meters,
independently of the length units that are used in the finite element analysis
(e.g., for nodal coordinates). The value of relMatAge must be specified in days.
Parameter relMatAge 28. means that time 0 of the analysis corresponds to
concrete age 28 days. If material B3mat is used, the finite element analysis
must use days as the units of time (not only for relMatAge, but in general, e.g.
for the time increments).
If only the basic creep (without shrinkage) should be computed, then the ma-
terial input record reduces to following: B3mat 1 n 0.2 d 0. talpha 1.2e-5
relMatAge 28. fc 45.4 cc 450. w/c 0.3778 a/c 4. t0 7. timefactor
1. shmode 0
Now consider the same conditions for “B3solidmat”. In all the examples
below, the input record with the material description can start by
B3solidmat 1 d 2420. n 0.2 talpha 12.e-6 begtimeofinterest 1.e-2
endtimeofinterest 3.e4 timefactor 86400. relMatAge 28.
Parameters begoftimeofinterest 1.e-2 and endoftimeofinterest 3.e4
mean that the computed response (e.g., deflection) should be accurate in the
range from 0.01 day to 30,000 days after load application. Parameter timefactor
86400. means that the time unit used in the finite element analysis is 1 second
56
(because 1 day = 86,400 seconds). Note that the values of begtimeofinterest,
endtimeofinterest and relMatAge are always specified in days, independently of
the actual time units in the analysis. Parameter EmoduliMode is not specified,
which means that the moduli of the Kelvin chain will be determined using the
default method, based on the continuous retardation spectrum.
Additional parameters depend on the specific type of analysis:
1. Computing basic creep only, shrinkage not considered, parameters qi esti-
mated from composition.
mode 0 fc 45.4 cc 450. w/c 0.3778 a/c 4.
t0 7. microprestress 0 shmode 0
2. Computing basic creep only, shrinkage not considered, parameters qi spec-
ified by the user.
mode 1 q1 18.81 q2 126.9 q3 0.7494 q4 7.692
microprestress 0 shmode 0
3. Computing basic creep only, shrinkage handled using the sectional ap-
proach, parameters estimated from composition.
mode 0 fc 45.4 cc 450. w/c 0.3778 a/c 4. t0 7.
microprestress 0 shmode 1 ks 1. alpha1 1. alpha2 1.2 hum 0.7
vs 0.1
4. Computing basic creep only, shrinkage handled using the sectional ap-
proach, parameters specified by the user.
mode 1 q1 18.81 q2 126.9 q3 0.7494 q4 7.692
microprestress 0 shmode 1 ks 1. q5 326.7 kt 28025 EpsSinf 702.4
t0 7. hum 0.7 vs 0.1
5. Computing basic creep only, shrinkage handled using the point approach
(B3), parameters specified by the user.
mode 1 q1 18.81 q2 126.9 q3 0.7494 q4 7.692
microprestress 0 shmode 2 es0 ... r ... rprime ... at ...
6. Computing drying creep, shrinkage handled using the point approach

(MPS), parameters qi estimated from composition.
mode 0 fc 45.4 w/c 0.3778 a/c 4. t0 7. microprestress 1
shmode 3 c0 1. c1 0.2 tS0 7. w h 0.0476 ncoeff 0.182 a 4.867
kSh 1.27258e-003
Finally consider the same conditions for “MPS material”. In all the examples
below, the input record with the material description can start by
mps 1 d 2420. n 0.2 talpha 12.e-6 referencetemperature 296.
Additional parameters depend on the specific type of analysis:
1. Computing basic creep only, shrinkage not considered, parameters qi es-
timated from composition and simulation time in days and stiffnesses in
MPa.
mode 0 fc 45.4 cc 450. w/c 0.3778 a/c 4. stiffnessFactor 1.e6
timefactor 1. lambda0 1. begoftimeofinterest 1.e-2
endoftimeofinterest 3.e4 relMatAge 28. CoupledAnalysisType 0.
57
2. Computing basic creep only, shrinkage not considered, parameters qi spec-
ified by user, simulation time in seconds and stiffnesses in Pa.
mode 1 q1 18.81e-12 q2 126.9e-12 q3 0.7494e-12 q4 7.6926e-12
timefactor 1. lambda0 86400. begoftimeofinterest 864.
endoftimeofinterest 2.592e9 relMatAge 2419200. CoupledAnalysisType
0.
3. Computing both basic and drying creep, parameters qi specified by user,
simulation time in seconds and stiffnesses in MPa.
mode 1 q1 18.81e-6 q2 126.9e-6 q3 0.7494e-6 q4 7.6926e-6
timefactor 1. lambda0 86400. begoftimeofinterest 864.
endoftimeofinterest 2.592e9 relMatAge 2419200. CoupledAnalysisType
1.
ksh 0.0004921875. t0 2419200. kappaT 0.005051 mus 4.0509259e-8
Final recommendations:
• to simulate basic creep without shrinkage it is possible to use all three

models
– B3mat with shmode = 0
– B3Solidmat with shmode = 0 and microprestress = 0
– MPS with CoupledAnalysisT ype = 0
• to simulate drying creep with shrinkage using “sectional approach”, only

the first material model (B3mat) is suitable can be used (with shmode =
1)
• to simulate drying creep without shrinkage using “sectional approach”,
only the first material model (B3mat) is suitable (with shmode = 1,
mode = 1 and EpsSinf = 0.0)
• to simulate drying creep with shrinkage using “point approach” according
to B3 model there are two options:
– B3mat (with shmode = 2)
– B3Solidmat (with shmode = 2)
In order to suppress shrinkage set es0 = 0.0
• to simulate drying creep with shrinkage using “point approach” according
to MPS model there are two options:
– B3Solidmat (with shmode = 3)

– MPS with CoupledAnalysisT ype = 1
In order to suppress shrinkage set kSh = 0.0
58
1.6.7 MPS damage model
This model extends the model based on the Microprestress-Solidification theory,
described in previous section and summarized in Tab. 37 for “MPS”, by tensile
cracking.
This extension uses the isotropic damage model with Rankine definition
of the equivalent deformation defined as the biggest principal effective stress
divided by the elastic modulus. The softening law cannot deal directly with
the strain because due to creep strain increases and this would lead to further
softening without additional loading.
Two different approaches are implemented. The first one, default, (#1)
reduces the stiffness only in the directions of tension (in case the tensile strength
is exceeded). A full stiffness is restored in compression and after unloading from
tension.
The other approach (#2) is the standard isotropic damage model which
reduces the stiffness equally in all directions independently of loading. This
approach leads to faster convergence because the secant stiffness can be used
instead of the incremental viscoelastic stiffness which must be used in the first
approach. The second approach becomes useful when the loading is monotonic
or when the benefit of the accelerated computation prevails over the conse-
quences of the reduced/underestimated stiffness in compression.
Proper energy dissipation is guaranteed by the crack-band approach.
The following algorithm is used to compute the stress vector (in each time
step and until the iteration criteria are met):
1. Compute effective stress
σeff,k+1 = σeff,k + Ēk D ν (∆εk − ∆ε00k − ∆εsh,k − ∆εT,k ) (130)
where σeff,k is the effective stress vector in the preceding step, Ēk is the
incremental stiffness, D ν is a unit stiffness matrix and ∆εk , ∆ε00k , ∆εsh,k ,
∆εT,k are the increments of the total strain, strain due to creep, shrink-
age strain (sum of drying and autogenous shrinkage) and thermal strain,
respectively.
2. Compute principal effective stresses σeff,1 , σeff,2 , σeff,3
3. Evaluate equivalent deformation
ε̃ = max(σeff,1 , σeff,2 , σeff,3 )/E (131)
4. If the stress exceeds the material strength, initialize and fix the fracture
parameters ft and GF (exponential softening assumed).
Gf
εc,f = (132)
ft h
where h is the characteristic length of the finite element in the direction
of the biggest principal stress.
5. Evaluate corresponding damage if ε̃ > ε0

ε0 ε̃ − ε0
ω =1− exp − (133)
ε̃ εc,f − ε0
59
6. Compute principal nominal stresses
• approach #1:
for i = 1, 2, 3
if σeff,i > 0, σi = (1 − ω)σeff,i
else σi = σeff,i
• approach #2:
for i = 1, 2, 3
σi = (1 − ω)σeff,i
7. Construct the stress vector in the original configuration
σ = T σ princ (134)
where T is the stress transformation matrix.

The tensile strength and fracture energy can be defined either as fixed, time-
independent values (parameters f t, gf ) or as variable, depending on the equiva-
lent hydration time te . The latter case is activated by keyword timeDepF racturing;
tensile strength and fracture energy are then linked to the current value of the
mean compressive strength using (slightly modified) formulae from the fib Model
Code 2010.
For 0MPa ≤ fcm (t) ≤ 20 MPa the tensile strength is computed from a linear
function
ftm (t) = 0.015 · 122/3 · fcm (t) = 0.07862fcm (t) (135)
which starts at the origin and at fcm (t) = 20 MPa connects to expression
ftm (t) = 0.3(fcm (t) − 8 MPa)2/3 (136)
which is valid for 20 MPa ≤ fcm (t) ≤ 58 MPa. Finally, for fcm (t) > 58 MPa
ftm (t) = 2.12 ln (1 + 0.1fcm (t)) (137)
Following the guidelines from the Model Code 2010 (section 5.1.5.2), the
fracture energy Gf can be estimated from the mean compressive strength at 28
days using an empirical formula
0.18
Gf,28 = 73 × (fcm,28 ) (138)
In this function, the mean compressive strength must be provided in MPa and
the resulting fracture energy is in N/m.
The growth in fracture energy is approximately proportional to the tensile
strength. The current value of fracture energy is very simply computed as
Gf (t) = Gf,28 × ftm (t)/ftm,28 (28) (139)
The development of the mean compressive strength in time is adopted en-

tirely from the Model Code 2010 (see section 5.1.9.1)
p
fcm (t) = fcm,28 exp f ib s 1 − 28/te (140)
where f ib s is a cement-type dependent coefficient (0.38 for cement 32.5 N; 0.25

for 32.5 R and 42.5 N; 0.2 for 42.5 R and 52.5 N and R), te is the equivalent
age and fcm,28 is the mean compressive strength at the age of 28 days.
60
For some concretes the prediction formulae do not provide accurate pre-
dictions, therefore it is possible to scale the evolution of tensile strength and
fracture energy by providing their 28-day values ft28 and gf28.
61
Description EC2CreepMat model for concrete creep and shrinkage
Record Format EC2CreepMat n(rn) # [ begOfTimeOfInterest(rn) #] [ end-
OfTimeOfInterest(rn) #] relMatAge(rn) # [ timeFactor(rn) #]
stiffnessFactor(rn) # [ tAlpha(rn) #] fcm28(rn) # t0(rn) # cem-
Type(in) # [ henv(rn) #] h0(rn) # shType(in) # [ spectrum ]
[ temperatureDependent ]
- begOfTimeOfInterest determines the shortest time which
is reasonably well captured by the approximated compli-
ance function (default value is 0.1); the units are the time
units of the analysis
- endOfTimeOfInterest determines the longest time which is
reasonably well captured by the approximated compliance
function (if not provided it is read from the engineering
model); the units are the time units of the analysis
- relMatAge time shift used to specify the age of material
on the begging of the analysis, the meaning is the material
age at time t = 0;
- timeFactor scaling factor transforming the actual time
into appropriate units needed by the formulae of the eu-
rocode. For analysis in days timeFactor = 1, for analysis
in seconds timeFactor = 86,400.
- stiffnessFactor scaling factor transforming predicted stiff-
ness into appropriate units of the analysis, for analysis in
MPa stiffnessFactor = 1.e6 (default), for Pa stiffnessFactor
=1
- fcm28 mean compressive strength measured on cylinders
at the age of 28 days in MPa
- t0 duration of curing [day] (this is relevant only for drying
shrinkage, not for creep)
- cemType type of cement, 1 = class R, 2 = class N, 3 =
class S
- henv ambient relative humidity expressed as decimal
- h0 effective member thickness in [mm] calculated accord-
ing to EC2 as 2×A/u where A is the cross-section are and
u is the cross-section perimeter exposed to drying
- shType shrinkage type; 0 = no shrinkage, 1 = both drying
and autogenous shrinkage, 2 = drying shrinkage only, 3 =
autogenous shrinkage only
- spectrum this string switches on evaluation of the moduli
of the aging Kelvin chain using the retardation spectrum
of the compliance function, otherwise (default option) the
least-squares method is used
- temperatureDependent turns on the influence of tempera-
ture on concrete maturity (equivalent age concept) by de-
fault this option is not activated.
Table 34: EC2Creep material model – summary.

62
Description B3 material model for concrete aging
Record Format B3mat d(rn) # n(rn) # talpha(rn) # [ begoftimeofinter-
est(rn) #] [ endoftimeofinterest(rn) #] timefactor(rn) # rel-
MatAge(rn) # [ mode(in) #] fc(rn) # cc(rn) # w/c(rn) #
a/c(rn) # t0(rn) # q1(rn) # q2(rn) # q3(rn) # q4(rn) #
shmode(in) # ks(rn) # vs(rn) # hum(rn) # [ alpha1(rn) #] [ al-
pha2(rn) #] kt(rn) # EpsSinf(rn) # q5(rn) # es0(rn) # r(rn) #
rprime(rn) # at(rn) # w h(rn) # ncoeff(rn) # a(rn) #
- n Poisson ratio
- talpha coefficient of thermal expansion
- begoftimeofinterest optional parameter; lower boundary of
time interval with good approximation of the compliance
function [day]; default 0.1 day
- endoftimeofinterest optional parameter; upper boundary
of time interval with good approximation of the compliance
function [day]
- timefactor scaling factor transforming the simulation time
units into days
- relMatAge relative material age [day]
- mode if mode = 0 (default value) creep and shrinkage
parameters are predicted from composition; for mode = 1
parameters must be user-specified.
- fc 28-day mean cylinder compression strength [MPa]
- cc cement content of concrete [kg/m3 ]
- w/c ratio (by weight) of water to cementitious material
- a/c ratio (by weight) of aggregate to cement
- t0 age when drying begins [day]
- q1-q4 parameters of B3 model for basic creep [1/TPa]
- shmode shrinkage mode; 0 = no shrinkage; 1 = average
shrinkage (the following parameters must be specified: ks,
vs, hum and additionally alpha1 alpha2 for mode = 0 and
kt EpsSinf q5 t0 for mode = 1; 2 = point shrinkage (needed:
es0, r, rprime, at, w h, ncoeff, a)
- ks cross-section shape factor [-]
- vs volume to surface ratio [m]
- hum relative humidity of the environment [-]
- alpha1 shrinkage parameter – influence of cement type [-]
- alpha2 shrinkage parameter – influence of curing type [-]
- kt shrinkage parameter [day/m2 ]
- EpsSinf shrinkage parameter [10−6 ]
- q5 drying creep parameter [1/TPa]
- es0 final shrinkage at material point
- r, rprime coefficients
- at oefficient relating stress-induced thermal strain and
shrinkage
- w h, ncoeff, a sorption isotherm parameters obtained from
experiments [Pedersen, 1990]
63
Table 35: B3 creep and shrinkage model – summary.
Description B3solid material model for concrete creep
Record Format B3solidmat d(rn) # n(rn) # talpha(rn) # mode(in) # [ Emod-
uliMode(in) #] Microprestress(in) # shm(in) # [ begoftime-
ofinterest(rn) #] [ endoftimeofinterest(rn) #] timefactor(rn) #
relMatAge(rn) # fc(rn) # cc(rn) # w/c(rn) # a/c(rn) # t0(rn) #
q1(rn) # q2(rn) # q3(rn) # q4(rn) # c0(rn) # c1(rn) # tS0(rn) #
w h(rn) # ncoeff(rn) # a(rn) # ks(rn) # [ alpha1(rn) #] [ al-
pha2(rn) #] hum(rn) # vs(rn) # q5(rn) # kt(rn) # EpsSinf(rn) #
es0(rn) # r(rn) # rprime(rn) # at(rn) # kSh(rn) # inithum(rn) #
finalhum(rn) #
- n Poisson ratio
- mode optional parameter; if mode = 0 (default), parame-
ters q1 − q4 are predicted from composition of the concrete
mixture (parameters fc, cc, w/c, a/c and t0 need to be spec-
ified). Otherwise values of parameters q1−q4 are expected.
- EmoduliMode optional parameter; analysis of retardation
spectrum (= 0, default value) or least-squares method (= 1)
is used for evaluation of Kelvin units moduli
- Microprestress 0 = basic creep; 1 = drying creep (must
be run as a staggered problem with preceding analysis of
humidity diffusion. Parameter shm must be equal to 3.
The following parameters must be specified: c0, c1, tS0,
w h, ncoeff, a)
- shmode shrinkage mode; 0 = no shrinkage; 1 = average
shrinkage (the following parameters must be specified: ks,
vs, hum and additionally alpha1 alpha2 for mode = 0 and
kt EpsSinf q5 for mode = 1; 2 = point shrinkage (needed:
es0, r, rprime, at), w h ncoeff a; 3 = point shrinkage based
on MPS theory (needed: parameter kSh or value of kSh can
be approximately determined if following parameters are
given: inithum, finalhum, alpha1 and alpha2)
- begoftimeofinterest optional parameter; lower boundary of
time interval with good approximation of the compliance
function [day]; default 0.1 day
- endoftimeofinterest optional parameter; upper boundary
of time interval with good approximation of the compliance
function [day]
- timefactor scaling factor transforming the simulation time
units into days
- relMatAge relative material age [day]
64
- fc 28-day mean cylinder compression strength [MPa]
- cc cement content of concrete mixture [kg/m3 ]
- w/c water to cement ratio (by weight)
- a/c aggregate to cement ratio (by weight)
- t0 age of concrete when drying begins [day]
- q1, q2, q3, q4 parameters (compliances) of B3 model for
basic creep [1/TPa]
- c0 MPS theory parameter [MPa−1 day−1 ]
- c1 MPS theory parameter [MPa]
- tS0 MPS theory parameter - time when drying begins
[day]
- w h, ncoeff, a sorption isotherm parameters obtained from
experiments [Pedersen, 1990]
- ks cross section shape factor [-]
- alpha1 optional shrinkage parameter - influence of cement
type (optional parameter, default value is 1.0)
- alpha2 optional shrinkage parameter - influence of curing
type (optional parameter, default value is 1.0)
- hum relative humidity of the environment [-]
- vs volume to surface ratio [m]
- q5 drying creep parameter [1/TPa]
- kt shrinkage parameter [day/m2 ]
- EpsSinf shrinkage parameter [10−6 ]
- es0 final shrinkage at material point
- at coefficient relating stress-induced thermal strain and
shrinkage
- rprime, r coefficients
- kSh influences magnitude of shrinkage in MPS theory [-]
- inithum [-], finalhum [-] if provided, approximate value of
kSh can be computed
Table 36: B3solid creep and shrinkage model – summary.
65
Description Microprestress-solidification theory material model for con-
crete creep
Record Format mps d(rn) # n(rn) # talpha(rn) # referencetempera-
ture(rn) # mode(in) # [ CoupledAnalysisType(in) #] [ be-
goftimeofinterest(rn) #] [ endoftimeofinterest(rn) #] timefac-
tor(rn) # relMatAge(rn) # lambda0(rn) # fc(rn) # cc(rn) #
w/c(rn) # a/c(rn) # stiffnessfactor(rn) # q1(rn) # q2(rn) #
q3(rn) # q4(rn) # ksh(rn) # [ mus(rn) #] k3(rn) # [ al-
phaE(rn) #] [ alphaR(rn) #] [ alphaS(rn) #] [ QEtoR(rn) #]
[ QRtoR(rn) #] [ QStoR(rn) #] kTm(rn) # [ kTc(rn) #]
[ p(rn) #] [ p tilde(rn) #] [ alpha as(rn) #] [ eps cas0(rn) #]
[ b4 eps au infty(rn) #] [ b4 tau au(rn) #] [ b4 alpha(rn) #]
[ b4 r t(rn) #] [ b4 cem type(rn) #] [ temperInCelsius ]
- n Poisson ratio
- referencetemperature reference temperature only to ther-
mal expansion of material
- mode optional parameter; if mode = 0 (default), param-
eters q1 − q4 are predicted from composition of the con-
crete mixture (parameters fc, cc, w/c, a/c and stiffnessfac-
tor need to be specified). Otherwise values of parameters
q1 − q4 are expected.
- CoupledAnalysisType 0 = basic creep; 1 = (default) dry-
ing creep, shrinkage, temperature transient creep and creep
at elevated temperature; 2 = drying creep, shrinkage; 3 =
temperature transient creep and creep at elevated tem-
perature; for choice # 1, 2, 3 the problem must be run
as a staggered problem with preceding analysis of humid-
ity and/or temperature distribution; Following parameters
must be specified: mus or k3 (according to exponent p),
kTm (compulsory for choice #3 otherwise optional)
- lambda0 scaling factor equal to 1.0 day in time units of
analysis (e.g. 86400 if the analysis runs in seconds)
- begoftimeofinterest lower boundary of time interval with
good approximation of the compliance function; default
value = 0.01 lambda0
- endoftimeofinterest upper boundary of time interval with
good approximation of the compliance function; default
value = 10000. lambda0
- timefactor scaling factor, for mps material must be equal
to 1.0
- relMatAge relative material age = age at time when the
material is cast in the structure
66
- fc 28-day standard cylinder compression strength [MPa]
- cc cement content of concrete mixture [kg m−3 ]
- w/c water to cement weight ratio
- a/c aggregate to cement weight ratio
- stiffnessfactor scaling factor converting “predicted” pa-
rameters q1 - q4 into proper units (e.g. 1.0 if stiffness is
measured in Pa, 1.e6 for MPa)
- q1, q2, q3, q4 parameters of B3 model for basic creep
- p and p tilde replaceable parameters in the governing
equation for viscosity, default value is p = 2
- mus parameter governing to the evolution of viscosity; for
exponent p = 2, µS = c0 c1 q4 [Pa−1 s−1 ]
- k3 dimensionless parameter governing to the evolution of
viscosity replacing mus in the special case when p > 100
(then p is automatically set to ∞ which is equivalent to
p tilde = 1)
- ksh parameter relating rate of shrinkage to rate of humid-
ity [-], default value is 0.0, i.e. no shrinkage
- alphaE constant, default value 10.
- alphaR constant, default value 0.1
- alphaS constant, default value 0.1
- QEtoR activation energy ratio, default value 2700. K
- QRtoR activation energy ratio, default value 5000. K
- QStoR activation energy ratio, default value 3000. K
- kTm replaces ln h in the governing equation for viscosity
- kTc controls creep at cyclic temperature
- alpha as and eps cas0 control the ultimate value of auto-
genous shrinkage according to fib Model Code 2010; this
ultimate value can be provided either directly via eps cas0
(negative value for contraction) or using alpha as, see equa-
tion (129)
- b4 eps au infty, b4 tau au, b4 alpha, b4 r t, b4 cem type
control the evolution and the ultimate value of autoge-
nous shrinkage according to model B4; all parameters are
predicted from composition if the b4 cem type is provided;
however, this prediction can be manually overridden
- temperInCelsius this string enables to run the supplemen-
tary transport problem with temperature in Celsius instead
of Kelvin
Table 37: MPS theory—summary.
67
Description MPS damage model for concrete creep with cracking
(additionaly all parameters from Table 37)
Record Format MPSDamMat [ ft(rn) #] [ gf(rn) #] [ timeDepFractur-
ing ] [ fib s(rn) #] [ ft28(rn) #] [ gf28(rn) #] [ damlaw(in) #]
[ isotropic ] [ checksnapback(in) #]
Parameters - ft tensile strength (constant in time)
- gf fracture energy (constant in time)
- timeDepFracturing string activating fib MC 2010 predic-
tion for tensile strength and fracture energy + their time
evolution
- ft28 manual override for the fib MC 2010 prediction of
hydration-dependent tensile strength
- gf28 manual override for the fib MC 2010 prediction of
hydration-dependent fracture energy
- fib s cement-type dependent coefficient (compulsory with
timeDepFracturing)
- damageLaw traction-separation law: 0 = exponential soft-
ening (default), 1 = linear softening, 6 = damage disabled
- isotropic string activating same reduction of stiffness after
the onset of cracking both in in compression and tension
(enables to use secant stiffness matrix for faster conver-
gence)
- checkSnapBack switch for snap-back checking; 1 = acti-
vated (default), 0 = deactiavated
Supported modes 3dMat, PlaneStress, PlaneStrain
Table 38: MPS damage–summary.
68
1.6.8 Microplane model M4 - Microplane M4
Model M4 covers inelastic behavior of concrete under complex triaxial stress
states. It is based on the microplane concept and can describe softening. How-
ever, objectivity with respect to element size is not ensured – the parameters
need to be manually adjusted to the element size. Since the tangent stiffness
matrix is not available, elastic stiffness is used. This can lead to a very slow
convergence when used within an implicit approach. The model parameters are
summarized in Tab. 39.
Description M4 material model

Record Format Microplane M4 nmp(in) # c3(rn) # c20(rn) # k1(rn) #
k2(rn) # k3(rn) # k4(rn) # E(rn) # n(rn) #
Parameters - nmp number of microplanes, supported values are 21, 28
and 61
- n Poisson ratio
- E Young modulus
- c3,c20, k1, k2, k3, k4 model parameters
Table 39: Microplane model M4 – summary.
1.6.9 Damage-plastic model for concrete - ConcreteDPM

This model, developed by Grassl and Jirásek for failure of concrete under general
triaxial stress, is described in detail in [8]. It belongs to the class of damage-
plastic models with yield condition formulated in terms of the effective stress
σ̄ = D e : (ε − εp ). The stress-strain law is postulated in the form
σ = (1 − ω)σ̄ = (1 − ω)D e : (ε − εp ) (141)
where D e is the elastic stiffness tensor and ω is a scalar damage parameter. The
plastic part of the model consists of a three-invariant yield condition, nonasso-
ciated flow rule and pressure-dependent hardening law. For simplicity, damage
is assumed to be isotropic. In contrast to pure damage models with damage
driven by the total strain, here the damage is linked to the evolution of plastic
strain.
The yield surface is described in terms of the cylindrical coordinates in the
principal effective stress space (Haigh-Westergaard coordinates), which are the
volumetricp effective stress σ̄V = I1 (σ̄)/3, the norm of the deviatoric effective
stress ρ̄ = 2J2 (σ̄), and the Lode angle θ defined by the relation
√
3 3 J3
cos 3θ = (142)
2 J 3/2
2
where J2 and J3 are the second and third deviatoric invariants. The yield
69
function
2 r !2
ρ̄ σ̄V 3 ρ̄
fp (σ̄V , ρ̄, θ̄; κp ) = [1 − qh (κp )] √ + ¯ + +
6f¯c fc 2 f¯c

2 ρ̄r(θ̄) σ̄V
+m0 qh (κp ) √ + ¯ − qh2 (κp ) (143)
6f¯c fc
depends on the effective stress (which enters in the form of cylindrical coordi-
nates) and on the hardening variable κp (which enters through a dimensionless
variable qh ). Parameter f¯c is the uniaxial compressive strength. Note that, un-
der uniaxial compression characterized by axial stress σ̄ < 0, we have σ̄V = σ̄/3,
ρ̄ = − 2/3 σ̄ and θ̄ = 60o . The yield function then reduces to fp = (σ̄/f¯c )2 −qh2 .
p
This means that function qh describes the evolution of the uniaxial compressive
yield stress normalized by its maximum value, f¯c .
The evolution of the yield surface during hardening is presented in Fig. 8.
The parabolic shape of the meridians (Fig. 8a) is controlled by the hardening
variable qh and the friction parameter m0 . The initial yield surface is closed,
which allows modeling of compaction under highly confined compression. The
initial and intermediate yield surfaces have two vertices on the hydrostatic axis
but the ultimate yield surface has only one vertex on the tensile part of the
hydrostatic axis and opens up along the compressive part of the hydrostatic
axis. The deviatoric sections evolve as shown in Fig. 8b, and their final shape
at full hardening is a rounded triangle at low confinement and almost circular
at high confinement. The shape of the deviatoric section is controlled by the
Willam-Warnke function
4(1 − e2 ) cos2 θ + (2e − 1)2
r(θ) = p (144)
2(1 − e2 ) cos θ + (2e − 1) 4(1 − e2 ) cos2 θ + 5e2 − 4e
The eccentricity parameter e that appears in this function, as well as the fric-
tion parameter m0 , are calibrated from the values of uniaxial and equibiaxial
compressive strengths and uniaxial tensile strength.
(a) (b)
θ=0
4
3
θ=0
2
1
ρ/fc
-1
-2
θ = 4π/3 θ = 2π/3
-3 θ=π
-4
-4 -3 -2 -1 0
σV/fc
Figure 8: Evolution of the yield surface during hardening: a) meridional section,

b) deviatoric section for a constant volumetric effective stress of σ̄V = −f¯c /3
The maximum size of the elastic domain is attained when the variable qh is
equal to one (which is its maximum value, as follows from the hardening law,
70
to be specified in (149)). The yield surface is then described by the equation
3 ρ̄2

ρ̄ σ̄V
fp σ̄V , ρ̄, θ̄; 1 ≡ ¯2 + m0 √ r(θ̄) + ¯ − 1 = 0 (145)
2 fc 6f¯c fc
The flow rule
∂gp
ε̇p = λ̇ (146)
∂ σ̄
is non-associated, which means that the yield function fp and the plastic po-
tential
2 r !2
ρ̄ σ̄V 3 ρ̄
gp (σ̄V , ρ̄; κp ) = [1 − qh (κp )] √ + ¯ + +
6f¯c fc 2 f¯c

2 m0 ρ̄ mg (σ̄V )
+qh (κp ) √ + (147)
6f¯c f¯c
do not coincide and, therefore, the direction of the plastic flow ∂gp /∂ σ̄ is not
normal to the yield surface. The ratio of the volumetric and the deviatoric
parts of the flow direction is controled by function mg , which depends on the
volumetric stress and is defined as
σ̄V − f¯t /3
mg (σ̄V ) = Ag Bg f¯c exp (148)
Bg f¯c
where Ag and Bg are model parameters that are determined from certain as-
sumptions on the plastic flow in uniaxial tension and compression.
The dimensionless variable qh that appears in the yield function (143) is a
function of the hardening variable κp . It controls the size and shape of the yield
surface and, thereby, of the elastic domain. The hardening law is given by
qh 0 + (1 − qh 0 )κp (κp 2 − 3κp + 3) if κp < 1

qh (κp ) = (149)
1 if κp ≥ 1
The initial inclination of the hardening curve (at κp = 0) is positive and finite,
and the inclination at peak (i.e., at κp = 1) is zero.
The evolution law for the hardening variable,2
kε̇p k
κ̇p = (2 cos θ̄)2 (150)
xh (σ̄V )
sets the rate of the hardening variable equal to the norm of the plastic strain
rate scaled by a hardening ductility measure

 Ah − (Ah − Bh ) exp (−Rh (σ̄V )/Ch ) if Rh (σ̄V ) ≥ 0
xh (σ̄V ) = (151)
Eh exp(Rh (σ̄V )/Fh ) + Dh if Rh (σ̄V ) < 0

The dependence of the scaling factor xh on the volumetric effective stress σ̄V
is constructed such that the model response is more ductile under compression.
The variable
σ̄V 1
Rh (σ̄V ) = − ¯ − (152)
fc 3
2 In the original paper [8], equation (150) was written with cos2 θ̄ instead of (2 cos θ̄)2 , but
all the results presented in that paper were computed with OOFEM using an implementation
based on (150).
71
is a linear function of the volumetric effective stress. Model parameters Ah , Bh , Ch
and Dh are calibrated from the values of strain at peak stress under uniaxial
tension, uniaxial compression and triaxial compression, whereas the parameters
Eh = Bh − Dh (153)
(Bh − Dh ) Ch
Fh = (154)
B h − Ah
are determined from the conditions of a smooth transition between the two parts
of equation (151) at Rh = 0.
For the present model, the evolution of damage starts after full saturation
of plastic hardening, i.e., at κp = 1. This greatly facilitates calibration of model
parameters, because the strength envelope is fully controled by the plastic part
of the model and damage affects only the softening behavior. In contrast to
pure damage models, damage is assumed to be driven by the plastic strain,
more specifically by its volumetric part, which is closely related to cracking.
To slow down the evolution of damage under compressive stress states, the
damage-driving variable κd is not set equal to the volumetric plastic strain, but
it is defined incrementally by the rate equation

0 if κp < 1
κ̇d = (155)
Tr(ε̇p )/xs (σ̄V ) if κp ≥ 1
where  2
 1 + As Rs (σ̄V ) if Rs (σ̄V ) < 1
xs (σ̄V ) = p (156)
1 − 3As + 4As if Rs (σ̄V ) ≥ 1

Rs (σ̄V )
is a softening ductility measure. Parameter As is determined from the softening
response in uniaxial compression. The dimensionless variable Rs = ε̇− pV /ε̇pV is
defined as the ratio between the “negative” volumetric plastic strain rate
3
X
ε̇−
pV = h−ε̇pI i (157)
I=1
and the total volumetric plastic strain rate ε̇pV . This ratio depends only on
the flow direction ∂gp /∂ σ̄, and thus Rs can be shown to be a unique function
of the volumetric effective stress. In (157), ε̇pI are the principal components of
the rate of plastic strains and h·i denotes the McAuley brackets (positive-part
operator). For uniaxial tension, for instance, all three principal plastic strain
rates are nonnegative, and so ε̇− pV = 0, Rs = 0 and xs = 1. This means that
under uniaxial tensile loading we have κd = κp − 1. On the other hand, under
compressive stress states the negative principal plastic strain rates lead to a
ductility measure xs greater than one and the evolution of damage is slowed
down. It should be emphasized that the flow rule for this specific model is
constructed such that the volumetric part of plastic strain rate at the ultimate
yield surface cannot be negative.
The relation between the damage variable ω and the internal variable κd
(maximum level of equivalent strain) is assumed to have the exponential form
ω = 1 − exp (−κd /εf ) (158)
72
Description Damage-plastic model for concrete
Record Format ConcreteDPM d(rn) # E(rn) # n(rn) # tAlpha(rn) # ft(rn) #
fc(rn) # wf(rn) # Gf(rn) # ecc(rn) # kinit(rn) # Ahard(rn) #
Bhard(rn) # Chard(rn) # Dhard(rn) # Asoft(rn) # helem(rn) #
href(rn) # dilation(rn) # yieldtol(rn) # newtoniter(in) #
- E Young modulus
- n Poisson ratio
- ft uniaxial tensile strength ft
- fc uniaxial compressive strength
- wf parameter wf that controls the slope of the softening
branch (serves for the evaluation of εf = wf /h to be used
in (158))
- Gf fracture energy, can be specified instead of wf, it is
converted to wf = Gf /ft
- ecc eccentricity parameter e from (144), optional, default
value 0.525
- kinit parameter qh 0 from (149), optional, default value 0.1
- Ahard parameter Ah from (151), optional, default value
0.08
- Bhard parameter Bh from (151), optional, default value
0.003
- Chard parameter Ch from (151), optional, default value 2
- Dhard parameter Dh from (151), optional, default value
10−6
- Asoft parameter As from (156), optional, default value 15
- helem element size h, optional (if not specified, the actual
element size is used)
- href reference element size href , optional (if not specified,
the standard adjustment of the damage law is used)
- dilation dilation factor (ratio between lateral and axial
plastic strain rates in the softening regime under uniaxial
compression), optional, default value -0.85
- yieldtol tolerance for the implicit stress return algorithm,
optional, default value 10−10
- newtoniter maximum number of iterations in the implicit
stress return algorithm, optional, default value 100
Table 40: Damage-plastic model for concrete – summary.
where εf is a parameter that controls the slope of the softening curve. In fact,
equation (158) is used by the nonlocal version of the damage-plastic model,
with κd replaced by its weighted spatial average (not yet available in the public
version of OOFEM). For the local model, it is necessary to adjust softening
according to the element size, otherwise the results would suffer by pathological
mesh sensitivity. It is assumed that localization takes place at the peak of
the stress-strain diagram, i.e., at the onset of damage. After that, the strain
73
is decomposed into the distributed part, which corresponds to unloading from
peak, and the localized part, which is added if the material is softening. The
localized part of strain is transformed into an equivalent crack opening, w, which
is under uniaxial tension linked to the stress by the exponential law
σ = (1 − ω)ft = ft exp (−w/wf ) (159)
Here, ft is the uniaxial tensile strength and wf is the characteristic crack open-
ing, playing a similar role to εf . Under uniaxial tension, the localized strain can
be expressed as the sum of the post-peak plastic strain (equal to variable κd )
and the unloaded part of elastic strain (equal to ωft /E). Denoting the effective
element size as h, we can write
w = h(κd + ωft /E) (160)
and substituting this into (159), we obtain a nonlinear equation

h
1 − ω = exp − (κd + ωft /E) (161)
wf
from which the damage variable ω corresponding to the given internal variable
κd can be computed by Newton iteration. The effective element size h is ob-
tained by projecting the element onto the direction of the maximum principal
strain at the onset of cracking, and afterwards it is held fixed. The evaluation
of ω from κd is no longer explicit, but the resulting load-displacement curve of
a bar under uniaxial tension is totally independent of the mesh size. A simpler
approach would be to use (158) with εf = wf /h, but then the scaling would not
be perfect and the shape of the load-displacement curve (and also the dissipated
energy) would slightly depend on the mesh size. With the present approach,
the energy per unit sectional area dissipated under uniaxial tension is exactly
Gf = wf ft . The input parameter controling the damage law can be either the
characteristic crack opening wf , or the fracture energy Gf . If both are specified,
wf is used and Gf is ignored. If only Gf is specified, wf is set to Gf /ft .
The onset of damage corresponds to the peak of the stress-strain diagram
under proportional loading, when the ratios of the stress components are fixed.
This is the case e.g. for uniaxial tension, uniaxial compression, or shear under
free expansion of the material (with zero normal stresses). However, for shear
under confinement the shear stress can rise even after the onset of damage, due
to increasing hydrostatic pressure, which increases the mobilized friction. It has
been observed that the standard approach leads to strong sensitivity of the peak
shear stress to the element size. To reduce this pathological effect, a modified
approach has been implemented. The second-order work (product of stress in-
crement and strain increment) is checked after each step and the element-size
dependent adjustment of the damage law is applied only after the second-order
work becomes negative. Up to this stage, the damage law corresponds to a fixed
reference element size, which is independent of the actual size of the element.
This size is set by the optional parameter href. If this parameter is not speci-
fied, the standard approach is used. For testing purposes, one can also specify
the actual element size, helem, as a “material property”. If this parameter
is not specified, the element size is computed for each element separately and
represents its actual size.
74
The damage-plastic model contains 15 parameters, but only 6 of them need
to be actually calibrated for different concrete types, namely Young’s modulus
E, Poisson’s ratio ν, tensile strength ft , compressive strength fc , parameter wf
(or fracture energy Gf ), and parameter As in the ductility measure (156) of the
damage model. The remaining parameters can be set to their default values
specified in [8].
The model parameters are summarized in Tab. 40. Note that it is possible
to specify the “size” of finite element, h, which (if specified) replaces the actual
element size in (161). The usual approach is to consider h as the actual element
size (evaluated automatically by OOFEM), in which case the optional parameter
h is missing (or is set to 0., which has the same effect in the code). However, for
various studies of mesh sensitivity it is useful to have the option of specifying h
as an input “material” value.
If the element is too large, it may become too brittle and local snap-back
occurs in the stress-strain diagram, which is not acceptable. In such a case, an
error message is issued and the program execution is terminated. The maximum
admissible element size
EGf Ewf
hmax = 2 = (162)
ft ft
happens to be equal to Hillerborg’s characteristic material length. For typical
concretes it is in the order of a few hundred mm. If the condition h < hmax is
violated, the mesh needs to be refined. Note that the effective element size h
is obtained by projecting the element. For instance, if the element is a cube of
edge length 100 mm, its effective size in the direction of the body diagonal can
be 173 mm.
1.6.10 CDPM2
This model is an extension of the ConcreteDPM presented in 1.6.9. CDPM2 has
been developed by Grassl, Xenos, Nyström, Rempling and Gylltoft for modelling
the failure of concrete for both static and dynamic loading. It is is described
in detail in [9]. The main differences between CDPM2 and ConcreteDPM are
that in CDPM2 the plasticity part exhibits hardening once damage is active.
Furthermore, two independent damage parameters describing tensile and com-
pressive damage are introduced. The parameters of CDPM2 are summarised in
Tab. 41
The stress for the anisotropic damage plasticity model (ISOFLAG=0) is
defined as
σ = (1 − ωt ) σ̄ t + (1 − ωc ) σ̄ c (163)
where σ̄ t and σ̄ c are the positive and negative parts of the effective stress tensor
σ̄, respectively, and ωt and ωc are two scalar damage variables, ranging from 0
(undamaged) to 1 (fully damaged).
The stress for the isotropic damage plasticity model (ISOFLAG=1) is defined
as
σ = (1 − ωt ) σ̄ (164)
The effective stress σ̄ is defined according to the damage mechanics convention
as
σ̄ = De : (ε − εp ) (165)
Plasticity:
75
The yield surface is described by the Haigh-Westergaard coordinates: the
volumetric effective stress σV , the norm of the deviatoric effective stress ρ and
the Lode angle θ. The yield surface is
( 2 r )2
ρ̄ σ̄V 3 ρ̄
fp (σ̄V , ρ̄, θ̄; κp ) = [1 − qh1 (κp )] √ + +
6fc fc 2 fc
(166)
2 ρ̄ σ̄V 2 2
+ m0 qh1 (κp )qh2 (κp ) √ r(cos θ̄) + − qh1 (κp )qh2 (κp )
6fc fc
It depends also on the hardening variable κp (which enters through the di-
mensionless variables qh1 and qh2 ). Parameter fc is the uniaxial compressive
strength. For qh2 = 1, the yield function is identical to the one of CDPM. The
shape of the deviatoric section is controlled by the Willam-Warnke function
4(1 − e2 ) cos2 θ̄ + (2e − 1)2

r(cos θ̄) = p (167)
2(1 − e2 ) cos θ̄ + (2e − 1) 4(1 − e2 ) cos2 θ̄ + 5e2 − 4e
Here, e is the eccentricity parameter. The friction parameter m0 is given by

3 fc2 − ft2 e
m0 = (168)
fc ft e+1
where ft is the tensile strength.
The flow rule (26) is split into a volumetric and a deviatoric part, i.e., the
gradient of the plastic potential is decomposed as
∂g ∂g ∂ σ̄V ∂g ∂ ρ̄
m= = + (169)
∂ σ̄ ∂ σ̄V ∂ σ̄ ∂ ρ̄ ∂ σ̄
Taking into account that ∂ σ̄V /∂ σ̄ = δ/3 and ∂ ρ̄/∂ σ̄ = s̄/ρ̄, restricting atten-
tion to the post-peak regime (in which qh1 = 1) and differentiating the plastic
potential (147), we rewrite equation (169) as

∂g ∂mg δ 3 m0 s̄
m= = + +√ (170)
∂ σ̄ ∂ σ̄V 3fc fc 6ρ̄ fc
The flow rule is non-associative which means that the direction of the plastic
flow is not normal to the yield surface. This is important for concrete since an
associative flow rule would give an overestimated maximum stress for passive
confinement.
The dimensionless variables qh1 and qh2 that appear in (143), (147) and
(148) are functions of the hardening variable κp . They control the evolution of
the size and shape of the yield surface and plastic potential. The first hardening
law qh1 is

qh0 + (1 − qh0 ) κ3p − 3κ2p + 3κp − Hp κ3p − 3κ2p + 2κp

if κp < 1
qh1 (κp ) =
1 if κp ≥ 1
(171)
The second hardening law qh2 is given by

1 if κp < 1
qh2 (κp ) = (172)
1 + Hp (κp − 1) if κp ≥ 1
76
The evolution law for the hardening variable,
kε˙p k 2 λ̇kmk 2
κ̇p = 2 cos θ̄ = 2 cos θ̄ (173)
xh (σ̄V ) xh (σ̄V )
sets the rate of the hardening variable equal to the norm of the plastic strain
rate scaled by a hardening ductility measure, which is identical to the one used
for the CDPM.
Damage:
Damage is initiated when the maximum equivalent strain in the history of
the material reaches the threshold ε0 = ft /E. This expression is determined
from the yield surface (fp = 0) by setting qh1 = 1 and qh2 = ε̃/ε0 . From this
quadratic equation for ε̃, the equivalent strain is determined as
s 2 2 2
3ε2 ρ̄2

ε0 m 0 ρ̄ σ̄V ε0 m0 ρ̄ σ̄V
ε̃ = √ r (cos θ) + + √ r (cos θ) + + 02
2 6fc fc 4 6fc fc 2fc
(174)
Tensile damage is described by a stress-inelastic displacement law. For linear
and exponential damage type the stress value ft and the displacement value wf
must be defined. For the bi-linear type two additional parameters ft1 and wf1
are required.
For the compressive damage variable, an evolution based on an exponen-
tial stress-inelastic strain law is used. The stress versus inelastic strain in the
softening regime in compression is

εi
σ = ft exp − if 0 < εi (175)
εfc
where εfc is an inelastic strain threshold which controls the initial inclination
of the softening curve. The use of different damage evolution for tension and
compression is one important improvement over CDPM.
The history variables κdt1 , κdt2 , κdc1 and κdc2 depend on a ductility measure
xs , which takes into account the influence of multiaxial stress states on the
damage evolution. This ductility measure is given by
xs = 1 + (As − 1) Rs (176)
where Rs is  √
6σ̄V
− if σ̄V ≤ 0

Rs = ρ̄ (177)
0 if σ̄V > 0

and As is a model parameter.

Strain rate:
Concrete is strongly rate dependent. If the loading rate is increased, the
tensile and compressive strength increase and are more prominent in tension
then in compression. The dependency is taken into account by an additional
variable alphar . The rate dependency is included by scaling both the equivalent
strain rate and the inelastic strain. The rate parameter is defined by
αr = (1 − X)αrt + Xαrc (178)
77
where X is the continuous compression measure (= 1 means only compression,
= 0 means only tension).
The functions αrt and αrc depend on the input parameter fc0 . A recom-
mended value for fc0 is 10 MPa.
1.6.11 Fixed crack model for concrete - ConcreteFCM

Implementation of a fixed crack model. Uncracked material is modeled as
isotropic linear elastic characterized by Young’s modulus and Poisson’s ratio.
Cracking is initiated when principal stress reaches tensile strength. Further
loading is governed by a softening law. Proper amount energy dissipation is
guaranteed by the crack-band approach. Multiple cracking is allowed; the max-
imum number of allowed cracks is controlled by ncracks parameter. Only mutu-
ally perpendicular cracks are supported. If cracking occurs in more directions,
the behavior on the crack planes is considered to be independent. The secant
stiffness is used for unloading and reloading. In a compression regime, this
model correspond to an isotropic linear elastic material. The model supports
3 different options for stiffness reduction in shear after cracking, 2 options for
the limit of the maximum shear stress on a crack plane and 6 different laws for
postpeak behavior. The model parameters are summarized in Tab. 42.
Sample syntax for a fixed crack model with volume density 24 kN/m3 , ther-
mal dilation coefficient 12 × 10−6 K−1 , Young’s modulus 20 GPa, Poisson’s
ratio 0.2, fracture energy 100 N/m, tensile strength 2 MPa, linear softening,
constant shear retention factor β = 0.05, Collins’ shear strength (with compres-
sive strength 30 MPa, aggregate size 0.01 m) and all cracks contribute to the
shear stiffness; the analysis uses [m], [MPa] and [MN]:
ConcreteFCM 1 d 24.e-3 talpha 12.e-6 E 20000. n 0.2 Gf 100e-6 ft 2.0
softType 2 shearType 1 beta 0.05 shearStrengthType 2 fc 30 ag 0.01
lengthscale 1. multipleCrackShear
Description Fixed crack model for concrete

Record Format ConcreteFCM (in) # d(rn) # tAlpha(rn) # E(rn) #
n(rn) # [ ncracks(in) #] [ multipleCrackShear ]
[ crackSpacing(rn) #] [ softType(in) #] [ shearType(in) #]
[ shearStrengthType(in) #] [ ecsm(rn) #] [ Gf(rn) #] [ ft(rn) #]
[ beta(rn) #] [ sf(rn) #] [ fc(rn) #] [ ag(rn) #] [ length-
scale(rn) #] [ soft w(ra) #] [ soft(w)(ra) #] [ soft eps(ra) #]
[ soft(eps)(ra) #] [ beta w(ra) #] [ beta(w)(ra) #] [ H(rn) #]
[ eps f(rn) #]
Parameters - material model number
- ncracks maximum allowed number of cracks
- crackSpacing specified distance between parallel cracks
- multipleCrackShear if not given, shear stiffness computed
from the dominant crack, otherwise all cracks contribute
78
- softType allows to select suitable softening law:
0 -
no softening (default)
1 -
exponential softening with parameters Gf and ft
2 -
linear softening with parameters Gf and ft
3 -
Hordijk softening with parameters Gf and ft
4 -
user-defined wrt crack opening with parameters ft,
soft w, and soft(w)
5 - linear hardening wrt strain with parameters ft, H, and
optionally eps f
6 - user-defined wrt strain with parameters ft, soft eps,
and soft(eps)
- shearType offers to choose from different approaches for

shear stiffness reduction of a cracked element
0 - no shear reduction (default)
1 - constant shear retention factor with parameter beta
2 - constant shear factor coefficient with parameter sf
3 - user-defined shear retention factor with parameters
beta w and beta(w)
- shearStrengthType allows to select a shear stress limit on

a crack plane
0 - no stress limit (default)

1 - constant strength = ft
2 - Collins interlock with parameters fc, ag, and length-
scale
- ecsm method used for evaluation of characteristic element

size L: 1 = square root of area, 2 = projection centered, 3
= Oliver, 4 = Oliver modified, 0 (default) = projection
- beta shear retention factor
- sf shear factor coefficient
- fc compressive strength in MPa
- ag aggregate size
- lengthscale factor to convert crack opening and aggregate
size in case of Collins aggregate interlock; 1 = analysis in
meters, 1000 = in millimeters, etc.
- soft w specified values of crack opening and
- soft(w) corresponding values of traction normalized to ft
- soft eps specified values of cracking strain and
- soft(eps) corresponding values of traction normalized to
ft
- beta w specified values of crack opening and
- beta(w) corresponding values of shear retention factor
- H hardening modulus (expressed wrt cracking strain)
79
- eps f threshold for cracking strain after which traction is
zero (applicable for linear hardening only)
Table 42: Fixed crack model for concrete – summary.
1.6.12 Fixed crack model for fiber reinforced composites - FRCFCM

This material model is an extension of the ConcreteFCM described in Sec-
tion 1.6.11.
It is possible to choose from three different “classes” of fibers. This choice
is controlled by keyword fiberType: 0 = continuous aligned fibers (CAF), 1 =
short aligned fibers (SAF) and 2 = short random fibers (SRF). Currently, it is
not possible to combine more classes of fibers in one material model.
All of the above-mentioned fiber types are further defined by the mate-
rial properties and geometry. Fiber quantity is captured by the dimensionless
volume fraction Vf (as decimal). All fibers are assumed to have a circular cross-
section (shape factor for shear kfib with default value 0.9) and to possess the
same geometry characterized by the diameter Df and length Lf.
The overall elastic stiffness of the fiber-reinforced composite is calculated as
a weighted average of the moduli of matrix Em and fibers Ef and the Poisson’s
ratio is considered to be equal to the Poisson’s ratio of the matrix. Similarly
to ConcreteFCM, cracking is initiated once the tensile stress σ in matrix reaches
the tensile strength ft .
The nominal bridging stress for continuous aligned fibers which are not per-
pendicular to the crack plane σb,f,θ can be very easily obtained by multiplying
the nominal bridging stress for perpendicular fibers σb,f by two terms: the first
one reflecting lower volume of inclined fibers passing through the crack plane
and the second one capturing the snubbing effect:
σb,f,θ = σb,f cos(θ) exp(θf ) (179)
where f is a snubbing coefficient.

For the CAF perpendicular to the crack, the nominal bridging stress can be
derived as s
Ef (1 + η)τ0
σb,f = 2Vf w̄ (180)
Df
where η = (Ef Vf )/[Em (1 − Vf )].
For SRF the nominal bridging stress is
s
Ef (1 + η)τ0 w̄ Vf Ef (1 + η)w̄
σb,f (w) = 2Vf − for w̄ < w∗ (181)
Df Lf
2
Vf Lf τs (w) 2w̄
σb,f (w) = 1− for w∗ ≤ w̄ < Lf /2 (182)
Df Lf
σb,f (w) = 0 for w̄ > Lf /2 (183)

where w∗ = L2f τ0 /[(1 + η)Ef Df ]; τ0 is the bond shear strength between the
fiber and matrix for small crack openings, w < w∗ .
80
Larger pull-out displacements can lead to significant physical changes in
the fiber surface which can result into changes in the bond shear stress. This
phenomenon is captured by function τs (w) relating the frictional bond to the
crack opening and is implemented in three alternative formulations. (In order to
keep τs (w) = τ0 use fssType = 0.) In conventional FRC with ordinary concrete
matrix, the frictional bond usually decreases with increasing slip. To capture
this type of behavior we adopt the function proposed by Sajdlová (activated
with fssType = 1) reads

|b0 |w
τs (w) = τ0 1 + sign(b0 ) 1 − exp − (184)
Df
where b0 is a micromechanical parameter. In composites with high-strength
matrix and coated high-strength steel fibers (HSFRC, UHPFRC) as well as in
SHCC materials with polymeric fibers, the frictional bond-slip relation often ex-
hibits hardening; this phenomenon can be well approximated by a cubic function
(activated with fssType = 2) proposed by Kabele
" 2 3 #
w w w
τs (w) = τ0 1 + b1 + b2 + b3 (185)
Df Df Df
or an alternative formulation which results in smooth changes in the bridging

stress (activated with fssType = 3)
" 2 3 #
w̃ w̃ w̃
τs (w) = τ̃0 + τ0 b1 + b2 + b3 (186)
Df Df Df
In the last two equations b1 , b2 and b3 are micromechanical parameters and

additionally in the last equation w̃ = w − w∗ , τ˜0 = τ0 (1. − w∗ /Lf )−2 for SRF
and τ˜0 = τ0 Ef (1 + η)Df /[Ef (1 + η)Df − 2Lf τ0 ] for SAF.
The bridging stress for SRF is defined as
r !
gVf Lf τ0 w̄ w̄
σb,f (w) = 2 − ∗ for w̄ < w∗ (187)
2Df w∗ w
2
gVf Lf τs (w) 2w
σb,f (w) = 1− for w∗ ≤ w̄ < Lf /2 (188)
2Df Lf
σb,f (w) = 0 for w̄ > Lf /2 (189)
where g is the snubbing factor is defined as

1 + exp(πf /2)
g=2 (190)
4 + f2
It is possible to delay the activation of the stress in fibers using prameteter
fibreActivationOpening which can be imagined as a “lag” of the fiber-related
crack opening behind the matrix-related crack opening.
During unloading the stress in fibers does not decrease linearly to origin.
Current implementation uses a power function
M
w̄
σb,f (w) = σb,f (wmax ) (191)
w̄max
81
where wmax is the maximum crack width reached in the entire previous history
and M is a positive constant, its default value is M = 4.
The influence of crack opening and sliding on the bridging shear stress only
due to fibers is expressed as
u V̄f kGf
τb,f = V̄f kGf = γcr (192)
wmax εcr,max
where V̄f is the effective volume of fibers crossing a crack plane (Vf /2 for SRF
and Vf cos(θ) for CAF and SAF), and Gf is the fiber shear modulus. This
expression is motivated by assumption that the fibers bridging the crack planes
behave as the Timoshenko beams subjected to shear. Note that the shear stiff-
ness of fibers is not recovered upon unloading.
It has been found that in some high performance fiber reinforced cement com-
posites, fibers rupture when cracks are exposed to shearing. This phenomenon
is modeled by damage parameter ω, which accounts for the ratio of ruptured
fibers and varies between the values of 0 and 1. It is assumed that ω depends
on the maximum shear strain sustained by the protruding portions of bridging
fibers throughout the loading history. This crack shear strain can be expressed
as:

|ui (t)|
γf,max = max . . . w(t) > ∆w (193)
max (wi (t))
where ui is the crack sliding displacement (CSD) and and wi is the maximum
value of the crack opening displacement of the i-th crack. This means that the
damage does not grow if the crack closes (crack opening decreases). If more
cracks exist, the maximum contribution is considered.
Two different one-parameter damage evolution laws are currently imple-
mented. For fDamType = 0 the damage is deactivated, with fDamType = 1
damage is described by

γf
ω(γf ) = min ,1 (194)
γf c
and finally with fDamType = 2

γf
ω(γf ) = 1 − exp − (195)
γf c
where γf c (gammaCrack in the input record) is a parameter.

Since damage reduces the number of crack-bridging fibers, which is propor-
tional to the fiber volume fraction, its effect can be suitably implemented by
introducing the effective volume fraction
Vf∗ = Vf (1 − ω) (196)
The material parameters are summarized in Tables 42 (matrix) and 43 (fiber

extension).
Sample syntax for a fixed crack model reinforced with fibers with volume
density 24 kN/m3 , thermal dilation coefficient 12 × 10−6 K−1 , Young’s modulus
of the matrix 20 GPa, Poisson’s ratio of matrix 0.2, fracture energy of matrix 100
82
N/m, tensile strength of matrix 2 MPa, linear tension softening, constant shear
retention factor β = 0.05, unlimited shear strength (shearStrengthType = 0),
continuous aligned fibers, fiber volume 2%, fiber diameter 0.04 mm, Young’s
modulus of fibers 20 GPa, shear modulus of fibers 1 GPa, fiber-matrix bond
strength 1 MPa, snubbing coefficient 0.7, shear correction coefficient 0.9, deac-
tivated fiber damage, fiber act if COD exceeds 10 µm (with smoothing from
w = 8to11 µm), fiber orientation at 45 degrees in x-y plane, automatic evalua-
tion of crack spacing from composition; the analysis uses [m], [MPa] and [MN]:
FRCFCM 1 d 24.e-3 talpha 12.e-6 E 20000. n 0.2 Gf 100e-6 ft 2.0
softType 2 shearType 1 beta 0.05 FiberType 0 Vf 0.02 Df 0.04e-3
Ef 20000. Gfib 1000. tau 0 1. FSStype 0 f 0.7 kfib 0.9 fDamType 0
fibreactivationopening 10.e-6 dw0 2.e-6 dw1 1.e-6 orientationVector
3 1. 1. 0. computeCrackSpacing
Description Fixed crack model for FRC

Record Format FRCFCM input record of ConcreteFCM Vf(rn) # Lf(rn) #
Df(rn) # Ef(rn) # [ nuf(rn) #] [ Gfib(rn) #] [ kfib(rn) #]
tau 0(rn) # b0(rn) # b1(rn) # b2(rn) # b3(rn) # f(rn) #
[ M(in) #] [ fibreOrientationVector(ra) #] [ fssType(in) #]
[ fDamType(in) #] [ fiberType(in) #] [ gammaCrack(rn) #]
[ computeCrackSpacing ] [ fibreActivationOpening(rn) #]
[ dw0(rn) #] [ dw1(rn) #]
Parameters - Vf fiber content expressed as decimal
- Lf fiber length
- Df fiber diameter
- Ef fiber Young’s modulus
- nuf fiber Poisson’s ratio
- Gfib fiber shear modulus (read when nuf is not provided)
- kfib fiber cross-sectional shape correction factor
- tau 0 bond shear strength at zero slip
- b0 micromechanical parameter for fiber shear according
to Sajdlov
- b1, b2, b3 micromechanical parameter for fiber shear ac-
cording to Kabele
- f snubbing friction coefficient
- M exponent related to fiber unloading
- fibreOrientationVector vector specifying orientation for
CAF and SAF fibers
- fssType type of Fiber bond Shear Strength (bond shear
strength vs. crack opening)
0 - constant shear strength
1 - bond shear strength with parameter b0
2 - bond shear strength with parameters b1, b2, b3
3 - bond shear strength with parameters b1, b2, b3 which
leads to smooth traction-separation law
83
- fDamType type of damage law for fibers
0 - no damage
1 - damage controlled by shear slip deformation of the
crack (with gammaCrack), linear law
2 - damage controlled by shear slip deformation of the
crack (with gammaCrack), exponential law
- fiberType type of reinforcing fibers

0 - CAF (continuous aligned fibers)
1 - SAF (short aligned fibers)
2 - SRF (short randomly oriented fibers)
- gammaCrack crack shear strain parameter applicable with

fDamType = 1 or 2 (here the crack shear strain is under-
stood as the crack slip u divided by the crack opening w)
- computeCrackSpacing crack spacing is evaluated automat-
ically based on provided composition
- fibreActivationOpening crack opening at which the fibers
begin transferring bridging stress
- dw0, dw1 applicable only if fibreActivationOpening 6= 0,
then it allows to smooth the traction-separation law for
fibers; lower bound is fibreActivationOpening - dw0 and the
upper bound is fibreActivationOpening - dw1
Table 43: Fixed crack model for fiber reinforced concrete – sum-
mary.
1.6.13 “Nonlocal” model for SHCC

The material parameters are summarized in Tables 42 (matrix), 43 (fiber ex-
tension), and 44 (nonlocal extension).
Description Nonlocal fixed crack model for FRC

Record Format FRCFCMNL input record of ConcreteFCM and FR-
CFCM r(rn) # wft(in) #
Parameters - r nonlocal radius (reasonable value is several millimeters
and its maximum is Lf /2 for short fibers)
- wft nonlocal averaging function, must be set to 4 (constant
function)
Table 44: Nonlocal fixed crack model for fiber reinforced concrete
– summary.
84
1.7 Orthotropic damage model with fixed crack orienta-
tions for composites - CompDamMat
The model is designed for transversely isotropic elastic material defined by five
elastic material constants. Typical example is a carbon fiber tow. Axis 1 rep-
resents the material principal direction. The orthotropic material constants are
defined as
ν12 = ν13 , ν21 = ν31 , ν23 = ν32 , E22 = E33 (197)

G12 = G13 = G21 = G31 , G23 = G32 (198)
ν12 = ν13 ν21 = ν31 ν31 = ν21 ν13 = ν12
= , = (199)
E11 E22 E33 E11
Material orientation on a finite element can be specified with lcs optional
parameter. If unspecified, material orientation is the same as the global coor-
dinate system. lcs array contains six numbers, where the first three numbers
represent a directional vector of the local x-axis, and the next three numbers
represent a directional vector of the local y-axis with the reference to the global
coordinate system. The composite material is extended to 1D and is also suit-
able for beams and trusses. In such particular case, the lcs has no effect and
the 1D element orientation is aligned with the global xx component.
The index p, p ∈ {11, 22, 33, 23, 31, 12} symbolizes six components of stress
or strain vectors. The linear softening occurs after reaching a critical stress fp,0 ,
see Fig. 9. Orientation of cracks is assumed to be orthogonal and aligned with
an orientation of material axes [3, pp.236]. The transverse isotropy is generally
lost upon fracture, material becomes orthotropic and six damage parameters dp
are introduced.
σ
fp,0
Gf,p
σ’p
(1-dp)Ep
εp,E εp εp,0 ε
Figure 9: Implemented stress-strain evolution with damage for 1D case. Tension
and compression are separated, but sharing the same damage parameter.
The compliance material matrix H, in the secant form and including damage
85
parameters, reads
1 ν21 ν31
−E −E 0 0 0
 
(1−d11 )E11 22 33
ν12 1 ν32
 −E (1−d22 )E22
−E 0 0 0 
 11 33 
ν13 ν23 1
−E −E 0 0 0
 
(1−d33 )E33
 
H= 11 22 
1

 0 0 0 (1−d23 )G23
0 0 

1
0 0 0 0 0
 
 (1−d31 )G31 
1
0 0 0 0 0 (1−d12 )G12
(200)
Damage occurs when any out of six stress tensor components exceeds a given
strength fp,0
|σp | ≥ |fp,0 | (201)
Positive and negative ultimate strengths can be generally different but share
the same damage variable. At the point of damage initiation, see Fig. 9, one
evaluates εp,E and characteristic element length lp , generally different for each
damage mode. Given the fracture energy GF,p , the maximum strain at zero
stress εp,0 is computed
2GF,p
εp,0 = εp,E + (202)
fp,0 lp
The point of damage initiation is never reached exactly, one needs to in-
terpolate between the previous equilibrated step and current step to achieve
objectivity.
The evolution of damage dp is based on the evolution of corresponding strain
εp . A maximum achieved strain is stored in the variable κp . If εp > κp the dam-
age may grow so the corresponding damage variable dp may increase. Desired
stress σp0 is evaluated from the actual strain εp
εp,0 − εp
σp0 = fp,0 (203)
εp,0 − εp,E
and the calculation of damage variables dp stems from Eq. 200, for example
0
σ11
d11 = 1− (204)
ν21 ν31
E11 ε11 + E22 σ22 + E33 σ33
0
σ12
d12 = 1− (205)
G12 ε12
Damage is always controled not to decrease. Fig. 10 shows a typical performance
for this damage model in one direction.
The damage initiation is based on a trial stress. It becomes necessary for
higher precision to skip a few first iteration, typically 5, and then to introduce
damage. A parameter afterIter is designed for this purpose and MinIter forces
a solver always to proceed certain amount of iterations.
allowSnapBack skips the checking of sufficient fracture energy for each di-
rection. If not specified, all directions are checked to prevent snap-back which
dissipates incorrect amount of energy.
86
σ22
da
ma
ng
g e
loadi
g
a din
unlo
ε22
da
ma
ge
Figure 10: Typical loading/unloading material performance for homogenized

stress and strain in the direction ‘2’. Note that one damage parameter is com-
mon for both tension and compression.
1.8 Orthotropic elastoplastic model with isotropic dam-

age - TrabBone3d
This model combines orthotropic elastoplasticity with isotropic damage. Mate-
rial orthotropy is described by the fabric tensor, i.e., a symmetric second-order
tensor with principal directions aligned with the axes of orthotropy and prin-
cipal values normalized such that their sum is 3. Elastic constants as well as
coefficients that appear in the yield condition are linked to the principal values
of the fabric tensor and to porosity. The yield condition is piecewise quadratic,
with different parameters in the regions of positive and negative volumetric
strain.
1.8.1 Local formulation

The basic equations include an additive decomposition of total strain into elastic
(reversible) part and plastic (irreversible) part
ε = εe + εp , (206)
σ = (1 − ω) σ̄ = (1 − ω) D : εe , (207)
the yield function √

f (σ̄, κ) = σ̄ : F : σ̄ − σY (κ). (208)
loading-unloading conditions
f (σ̄, κ) ≤ 0 λ̇ ≥ 0 λ̇f (σ̄, κ) = 0, (209)
evolution law for plastic strain

∂f
ε̇p = λ̇ , (210)
∂ σ̄
87
the incremental definition of cumulated plastic strain
κ̇ = kε̇p k, (211)
the law governing the evolution of the damage variable
ω(κ) = ωc (1 − e−aκ ), (212)
and the hardening law
σY (κ) = 1 + σH (1 − e−sκ ). (213)
In the equations above, σ̄ is the effective stress tensor, D is the elastic stiffness
tensor, f is the yield function, λ is the consistency parameter (plastic multiplier),
ω is the damage variable, σY is the yield stress and s, a, σH and ωc are positive
material parameters. Material anisotropy is characterized by the second-order
positive definite fabric tensor
3
X
M= mi (mi ⊗ mi ), (214)
i=1
normalized such that Tr(M ) = 3, mi are the eigenvalues and mi the eigenvec-
tors. The eigenvectors of the fabric tensor determine the directions of material
orthotropy and the components of the elastic stiffness tensor D are linked to
eigenvalues of the fabric tensor. In the coordinate system aligned with mi ,
i = 1, 2, 3, the stiffness can be presented in Voigt (engineering) notation as
−1
1
− νE121 − νE131

E1 0 0 0

 − νE212 1
E2 − νE232 0 0 0 

 − νE313 − νE323 1
E3 0 0 0 
D= , (215)
 
1
0 0 0 0 0

 G23 
1
0 0 0 0 0
 
 G13 
1
0 0 0 0 0 G12
ml
where Ei = E0 ρk m2l k l l
i , Gij = G0 ρ mi mj and νij = ν0 ml . Here, E0 , G0 and
i
j
ν0 are elastic constants characterizing the compact (poreless) material, ρ is the
volume fraction of solid phase and k and l are dimensionless exponents.
Similar relations as for the stiffness tensor are also postulated for the com-
ponents of a fourth-order tensor F that is used in the yield condition. The yield
condition is divided into tensile and compressive parts. Tensor F is different
in each part of the effective stress space. This tensor is denoted F + in tensile
part, characterized by N̂ : σ̄ ≤ 0 and F − in compressive part, characterized by
N̂ : σ̄ ≤ 0, where P3 −2q
i=1 mi
N̂ = qP (mi ⊗ mi ) (216)
3 −4q
i=1 m i
88
χ± χ±
 
1
2 − 12
2 − 13
2 0 0 0
σ1±
 ( ) ( σ1± ) ( σ1± ) 
χ± 1 χ±
 − 21
− 23
0 0 0 
 
2 2 2
 ( )σ2± ( )σ2± ( σ2± ) 
± ±
χ31 χ32
±
 1

F  −( )
=
σ3±
2 −
( )σ3±
2
(σ3± )
2 0 0 .
0  (217)
1
 

 0 0 0 τ23 0 0 

1
 0 0 0 0 τ13 0 
1
0 0 0 0 0 τ12
In the equation above σi± = σ0± ρp m2qi is uniaxial yield stress along the i-th
principal axis of orthotropy, τij = τ0 ρp mqi mqj is the shear yield stress in the
m2q
plane of orthotropy and χ± ± i
ij = χ0 m2q is the so-called interaction coefficient, p
j
and q are dimensionless exponents and parameters with subscript 0 are related
to a fictitious material with zero porosity. The yield surface is continuously
differentiable if the parameters values are constrained by the condition
χ−
0 +1 χ+
0 +1
− 2 = . (218)
(σ0 ) (σ0+ )2
1.8.2 Nonlocal formulation - TrabBoneNL3d

The model is regularized by the over-nonlocal formulation with damage driven
by a combination of local and nonlocal cumulated plastic strain
κ̂ = (1 − m)κ + mκ̄, (219)
where m is a dimensionless material parameter (typically m > 1) and

Z
κ̄(x) = α(x, s)κ(s) ds (220)
V
is the nonlocal cumulated plastic strain. The nonlocal weight function is defined
as
α0 (kx − sk)
α(x, s) = R (221)
α0 (kx − tk) dt
V
where ( 2
r2
α0 (r) = 1 − R 2 if r ≤ R (222)
0 if r > R
Parameter R is related to the internal length of the material. The model de-
scription and parameters are summarized in Tab. 47.
89
Description CDPM2
Record Format con2dpm d(rn) # E(rn) # n(rn) # tAlpha(rn) # ft(rn) #
fc(rn) # wf(rn) # [stype(in) #] [ft1(rn) #] [wf1(rn) #]
[efc(rn) #] [ecc(rn) #] [kinit(rn) #] [Ahard(rn) #] [Bhard(rn) #]
[Chard(rn) #] [Dhard(rn) #] [Asoft(rn) #] [helem(rn) #] [dila-
tion(rn) #] [hp(rn) #] [isoflag(in) #] [rateflag(in) #] fcZero(in) #
[yieldtol(rn) #] [newtoniter(in) #]
- E Young modulus
- n Poisson ratio
- ft uniaxial tensile strength
- fc uniaxial compressive strength
- wf parameter that controls the slope of the softening
branch to be used
- stype allows to choose different types of softening laws,
default value 1:
0 - linear softening
1 - bilinear softening
2 - exponential softening
- ft1 parameter for the bilinear softening law (stype = 1)

defining the ratio between intermediate stress and tensile
strength, optional, default value 0.3
- wf1 parameter for the bilinear softening law (stype = 1)
defining the intermediate crack opening, optional, default
value 0.15
- efc parameter for exponential softening law in compres-
sion, optional, default value 100 × 10−6
- ecc eccentricity parameter, optional, default value 0.525
- kinit initial value of hardening law, optional, default value
0.3
- Ahard hardening parameter, optional, default value 0.08
- Bhard hardening parameter, optional, default value 0.003
- Chard hardening parameter, optional, default value 2
- Dhard hardening parameter, optional, default value 1 ×
10−6
- Asoft softening parameter, optional, default value 15
- helem element size h, optional (if not specified, the actual
element size is used)
- dilation dilation factor, optional, default value 0.85
- hp hardening modulus, optional, default value 0.5
- isoflag flag which allows for the use of only one damage
parameter if set to 1, default value 0
- rateflag flag which allows for the consideration of the effect
of strain rates if set to 1, default value 0
- fcZero Input parameter for modelling the effect of strain
value. Recommended value 10 MPa. Only if rateflag=1,
fcZero has to be specified.
- yieldtol tolerance for the implicit stress return algorithm,
optional, default value 1 × 10−6
- newtoniter maximum90 number of iterations in the implicit
stress return algorithm, optional, default value 100
Supported modes 3dMat, PlaneStrain
Table 41: CDPM2 – summary.

Description Orthotropic damage model with fixed crack orientations for
composites
Record Format CompDamMat num(in) # d(rn) # Exx(rn) # EyyEzz(rn) #
nuxynuxz(rn) # nuyz(rn) # GxyGxz(rn) # Tension f0 Gf(ra) #
Compres f0 Gf(ra) # [afterIter(in) #] [allowSnap-
Back(() #ia)]
- Exx Young’s modulus for principal direction xx
- EyyEzz Young’s modulus in orthogonal directions to the
principal direction xx
- nuxynuxz Poisson’s ratio in xy and xz directions
- nuyz Poisson’s ratio in yz direction
- GxyGxz shear modulus in xy and xz directions
- Tension f0 Gf array with six pairs of positive numbers.
Each pair describes maximum stress in tension and fracture
energy for each direction (xx, yy, zz, yz, zx, xy)
- Compres f0 Gf array with six pairs of numbers. Each pair
describes maximum stress in compression (given as a nega-
tive number) and positive fracture energy for each direction
(xx, yy, zz, yz, zx, xy)
Supported modes 3dMat, 1dMat
- afterIter how many iterations must pass until damage is
computed from strains, zero is default. User must ensure
that the solver proceeds the minimum number of iterations.
- allowSnapBack array to skip checking for snap-back. The
array members are 1-6 for tension and 7-12 for compression
components.
Table 45: Orthotropic damage model with fixed crack orientations for compos-
ites – summary.
91
Description Anisotropic elastoplastic model with isotropic damage
Record Format TrabBone3d (in) # d(rn) # eps0(rn) # nu0(rn) # mu0(rn) #
expk(rn) # expl(rn) # m1(rn) # m2(rn) # rho(rn) #
sig0pos(rn) # sig0neg(rn) # chi0pos(rn) # chi0neg(rn) #
tau0(rn) # plashardfactor(rn) # expplashard(rn) # exp-
dam(rn) # critdam(rn) #
- eps0 Young modulus (at zero porosity)
- nu0 Poisson ratio (at zero porosity)
- mu0 shear modulus of elasticity (at zero porosity)
- m1 first eigenvalue of the fabric tensor
- m2 second eigenvalue of the fabric tensor
- rho volume fraction of solid phase
- sig0pos yield stress in tension
- sig0neg yield stress in compression
- tau0 yield stress in shear
- chi0pos interaction coefficient in tension
- plashardfactor hardening parameter
- expplashard exponent in hardening law
- expdam exponent in damage law
- critdam critical damage
- expk exponent k in the expression for elastic stiffness
- expl exponent l in the expression for elastic stiffness
- expq exponent q in the expression for tensor F
- expp exponent p in the expression for tensor F
Table 46: Anisotropic elastoplastic model with isotropic damage - summary.
92
Description Nonlocal anisotropic elastoplastic model with isotropic
damage
Record Format TrabBoneNL3d (in) # d(rn) # eps0(rn) # nu0(rn) #
mu0(rn) # expk(rn) # expl(rn) # m1(rn) # m2(rn) # rho(rn) #
sig0pos(rn) # sig0neg(rn) # chi0pos(rn) # chi0neg(rn) #
tau0(rn) # plashardfactor(rn) # expplashard(rn) # exp-
dam(rn) # critdam(rn) # m(rn) # R(rn) #
- eps0 Young modulus (at zero porosity)
- nu0 Poisson ratio (at zero porosity)
- mu0 shear modulus (at zero porosity)
- m1 first eigenvalue of the fabric tensor
- m2 second eigenvalue of the fabric tensor
- rho volume fraction of the solid phase
- sig0pos yield stress in tension
- tau0 yield stress in shear
- chi0pos interaction coefficient in tension
- chi0neg interaction coefficient in compression
- plashardfactor hardening parameter
- expplashard exponent in the hardening law
- expdam exponent in the damage law
- critdam critical damage
- expk exponent k in the expression for elastic stiffness
- expl exponent l in the expression for elastic stiffness
- expq exponent q in the expression for tensor F
- expp exponent p in the expression for tensor F
- R nonlocal interaction radius
Table 47: Nonlocal formulation of anisotropic elastoplastic model with isotropic

damage – summary.
93
1.9 Material models for interfaces
Interface elements have to be used with material models describing the consti-
tutive behavior of interfaces between two materials (e.g. between steel reinforce-
ment and concrete), or between two bodies in contact. Such interface laws are
formulated in terms of the traction vector t and the displacement jump vector
δ.
1.9.1 Cohesive interface material - cohInt

A simple interface material with generally different stiffness in compression,
tension and shear. It is intended for simulating a contact with low stiffness
in tension and high in compression. The parameter transitionOpening specifies
initial gap embedded in an element which causes transition from tension to
compression. For example, transitionOpening=0.01 m means that there exists
the embedded gap which needs displacement jump 0.01 m to close. Traction-
separation law for normal direction takes the following form:
π k π k
t n
tn = + arctan(sm δ̄) δ̄ + − arctan(sm δ̄) δ̄ (223)
2 π 2 π
sm = smoothMag (224)
δ0 = transitionOpening (225)
δ̄ = δn + δ0 (226)
kt = kn · stiffCoeffKn (227)
where smoothMag controlls smoothing magnitude. Tangential stiffness for nor-

mal direction is found by differentiating Eq. (223)
π k π k
t n
knn = + arctan(sm δ̄) + − arctan(sm δ̄) + (228)
2 π 2 π
sm · kt δ̄ sm · kn δ̄
+ − (229)
π(sm δ̄ + 1) π(s2m δ̄ 2 + 1)
2 2
Shear stiffness remains constant during all possible loadins and there is no
influence of normal direction.
Description Cohesive interface material

Record Format cohInt (in) # kn(rn) # ks(rn) # [ stiffCoeffKn(rn) #] [ smooth-
Mag(rn) #] [ transitionOpening(rn) #]
- kn (penalty) stiffness in compression
- ks stiffness in shear
- stiffCoeffKn ratio (tensile stiffness / compression stiffness)
- smoothMag smoothing parameter for transition between
tensile/compressive behavior
- transitionOpening embedded gap in the element
Supported modes 1dInterface, 2dInterface, 3dInterface
Table 48: Cohesive interface material – summary.
94
1.9.2 Isotropic damage model for interfaces
This model is described in Section 1.5.9.
1.9.3 Simple interface material - obsolete

This model provides a simple interface law with penalty-type contact and fric-
tion. In the normal direction, the response is linear elastic, but with dif-
ferent stiffnesses in tension and in compression (stiffness kn in compression,
kn*stiffcoeff in tension). By setting kn to a high value, the penetration (over-
lap) can be reduced, in the sense of the penalty approach. By setting stiffcoeff
to 0, free opening of the gap can be allowed. The shear response is elastoplastic,
with the yield limit dependent on the normal traction. Crosssection’s width,
height or area have no influence on the results. The magnitude of the shear
traction σ T must not exceed the yield limit computed according to a Coulomb-
like friction law as the product of the (negative part of) normal traction σN and
a dimensionless coefficient of friction f c:
||σ T || ≤ h−σN · f ci (230)
σN is computed multiplying a known normal strain and a known stiffness in
tension or in compression.
The shear elastic stiffness is assumed to be kn (i.e., equal to the normal
stiffness). If the normal traction is tensile or zero, no shear traction can be
transmitted by the interface. In the future, this law will be enriched by a
cohesion-like term.
If used in connection of two nodes with interface1delement, the element-
level parameter normal determines what is compression and what is tension
(from the 1st node to the 2nd node in normal direction it is tension). If nor-
malClearance is defined, the element has tensile stiffness until the gap within
the element is closed. This feature allows simulating two faces with a gap.
When the gap is smaller than normalClearance, the element stiffness changes to
compression, see Figure 11.
Description Simple interface material

Record Format simpleintermat (in) # kn(rn) # [ fc(rn) #] [ stiffcoeff(rn) #]
[ normalClearance(rn) #]
- kn (penalty) stiffness in compression
- fc friction coefficient
- stiffcoeff ratio (tensile stiffness / compression stiffness)
- normalClearance free distance within element
Supported modes 1dInterface
Table 49: Simple interface material – summary.
1.9.4 Bond-slip model for reinforced concrete - BondCEB

The complex steel/concrete interaction can be modelled with a bond-slip rela-
tion, which describes the bond stress (tangential interface traction), τ , in terms
95
SimpleInterMat Diagram
force F
Tensile
kn·stiffCoef
normalClearance 1
Relative node displacement u
kn
Figure 11: Working diagram of SimpleInterfaceMat.
of the relative reinforcement slip (tangential interface jump), s. This inter-

face material model is based on the local bond–slip relationship for reinforced
concrete under good bond conditions outlined in the the fib Model Code for
Concrete Structures 2010 [31]. In general, the model is formulated in terms of
the interface traction vector t and the interface jump vector δ. In the model,
the interface is assumed to have only elastic stiffness kn in the normal direction.
The normal interface traction is then evaluated elastically from the normal in-
terface jump, i.e. tn = kn δn . The tangential traction tt = τ is evaluated from
the user-specified function τ (s), see Figure 12:
α
s


τ max for 0 ≤ s ≤ s1 ,



 s1
τmax for s1 < s ≤ s2 ,

τ (s) = (231)
(τmax − τf ) (s − s2 )
τ − for s < s ≤ s ,

max 2 3

(s3 − s2 )





τf for s > s3 .
For an exponent α smaller than 1 (default value 0.4), the initial tangential
stiffness of the interface, ks , is undefined and cannot be used in numerical com-
putations. Therefore, this stiffness needs to be specified manually. Should the
τmax
provided stiffness be smaller than kn = , it will be automatically adjusted
s1
to this value. Note that only elastic tangent stiffness is supported. Hence, in
2D the tangent stiffness matrix takes the form
∂tn ∂tn
 

∂t  ∂δn ∂s  kn 0
=  ∂τ ∂τ  = (232)
∂δ 0 ks
∂δn ∂s
Model description and the input parameters are summarized in Tab. 50.
1.10 Material models for lattice elements

Lattice elements have to be used with material models describing the consti-
tutive behavior in the form of vector of tractions and strains determined from
96
Figure 12: Bond stress (τ ) - reinforcement slip (s) diagram for BondCEB model.
Description Bond-slip model for reinforced concrete

Record Format bondceb (in) # kn(rn) # ks(rn) # s1(rn) # s2(rn) # s3(rn) #
taumax(rn) # [ tauf(rn) #] [ alpha(rn) #]
- kn interface elastic normal stiffness
- ks interface elastic tangential (shear) stiffness
- s1 characteristic slip value
- taumax maximum bond stress
- tauf residual bond stress at reinforcement pull-out
- alpha bond-slip curve parameter (exponent)
Table 50: Bond-slip model for reinforced concrete – summary.
displacement jumps smeared out over the element length.
1.10.1 Latticedamage2d
This is a damage lattice material used together with latticedamage2d elements.
It uses a scalar damage relationship of the form
σ = (1 − ω) De ε (233)
T T
where σ = (σn , σt , σφ ) is a vector of tractions and ε = (εn , εt , εφ ) is a vector
of strains obtained from displacement jumps smeared over the element length.
Furthermore, ω is the damage variable varying from 0 (undamaged) to 1 (fully
damaged). Also, De is the elastic stiffness matrix which is based on the elastic
modulus of the lattice material E, and a parameter γ which is the ratio of the
modulus of the shear and normal direction. The strength envelope is elliptic
(Figure 13) and determined by three parameters, ft , fq and fc . The evolution
97
of the damage variable ω is controlled by normal stress-normal crack opening
law. The three possible laws are linear, bilinear and exponential (Figure 14).
Figure 13: Strength envelope of LatticeDamage2d.
(a) (b) (c)
Figure 14: Softening types of LatticeDamage2d: (a) linear softening, (b) bilinear
softening, (c) exponential softening.
The model parameters are summarised in Tab. 51.
1.11 Material models for steel relaxation

1.11.1 Model for relaxation of prestressing steel - SteelRelaxMat
This section describes the implementation of the material model for steel relax-
ation given in Eurocode 2 (the same as in Model Code 2010) and in Bažant and
Yu (J. of Eng. Mech, 2013) which reduces to the first model under constant
strain. At variable strain history the first model uses the approach employing
the so-called equivalent time approach described in Annex D in the Eurocode
2.
The current implementation takes into account only prestress losses due to
steel relaxation, other losses (e.g. slip at anchorage, thermal dilation, friction,
etc.) need to be treated separately. The same holds for the stress transfer from
prestressing reinforcement to concrete in the region called transmission length.
On the other hand, losses due to sequential prestressing, elastic deformation
and both short-time and long-time creep and shrinkage are taken into account
automatically provided that a suitable material model is chosen for concrete.
In the first approach the stress on the end of the time step is explicitly given
by the current stress, prestressing level and the cumulative value of prestress
98
Description Saclar damage model for lattice2d
Record Format latticedamage2d (in) # d(rn) # talpha(rn) # e(rn) # a1(rn) #
a2(rn) # e0(rn) # coh(rn) # ec(rn) # stype(rn) # wf(rn) #
wf1(rn) #
- talpha Thermal exansion coefficient
- e normal modulus of lattice material
- a1 ratio of shear and normal modulus
- a2 ratio of rotational and normal modulus. Optional pa-
rameter. Default is 1.
- e0 strain at tensile strength: ft /E
- coh ratio of shear and tensile strength: fq /ft
- ec ratio of compressive and tensile strength: fc /ft
- stype softening types: 1-linear, 2-bilinear and 3-
exponential
- wf displacement threshold related to fracture energy used
in all three softening types.
Supported modes 2dlattice
Table 51: Scalar damage model for 2d lattice elements – summary.
losses. On the other hand, in Bažant’s approach it is necessary to iterate on the

material point level in order to reach equilibrium.
As a simplification the stress-strain diagram is in the current implementation
assumed to be linear (no yielding), this should be sufficient for most cases.
Under a constant strain, the evolution of prestress loss is defined as
∆σ = σinit k1 ρ1000 exp(k2 µ)(t/1000)0.75(1−µ) × 10−5 (234)
where σinit is the initial value of prestress reduced for losses during prestressing,
t is time after prestressing in hours, µ = σinit /fpk , fpk is the characteristic
strength of prestressing steel in tension, and finally k1 , k2 , and ρ1000 are material
parameters determined by the relaxation properties of the reinforcement. For
wires or cables with normal relaxation (class 1) k1 = 5.39, k2 = 6.7 and
ρ1000 = 8.0, for cables or wires with reduced relaxation (class 2) k1 = 0.66,
k2 = 9.1 and ρ1000 = 2.5, and for hot-rolled and modified rods (class 3)
k1 = 1.98, k2 = 8.0 and ρ1000 = 4.0.
The prestress σinit is not specified in the input record. It is initialized
automatically at the time instant when stress differs from zero.
The material model has one internal variable which has a meaning of cumu-
lative prestress loss when equivalent time approach is employed; otherwise its
meaning is a cumulative strain caused by relaxation.
99
Description SteelRelaxMat model for relaxation of prestressing rein-
forcement
Record Format SteelRelaxMat d(rn) # E(rn) # reinfClass(in) # [ time-
Factor(rn) #] charStrength(rn) # approach(in) # [ k1(rn) #]
[ k2(rn) #] [ rho1000(rn) #] [ tolerance(rn) #] [ relRe-
laxBound(rn) #]
- d specific weight
- reinfClass class of prestressing reinforcement (1, 2, 3)
- timeFactor scaling factor transforming the actual time
into appropriate units needed by the formulae of the eu-
rocode. For analysis in days timeFactor = 1, for analysis
in seconds timeFactor = 86,400.
- charStrength characteristic strength of prestressing steel
in appropriate units (not necessarily MPa)
- approach 0 = approach according to Bažant and Yu, 1 =
equivalent time approach according to Eurocode 2 and fib
Model Code 2010
- k1 possibility to overwrite default value given by the re-
inforcement class
- k2 possibility to overwrite default value given by the re-
inforcement class
- rho1000 possibility to overwrite default value given by the
reinforcement class
- tolerance applicable only for approach = 0; tolerance spec-
ifying the residual in the stress evaluation algorithm, de-
fault value is 10−6
- relRelaxBound ratio of stress to characteristic strength un-
der which the relaxation is zero (typically 0.4–0.5); default
value is zero
Table 52: SteelRelaxMat material model – summary.
100
2 Material Models for Transport Problems
2.1 Isotropic linear material for heat transport – IsoHeat
Linear isotropic material model for heat transport problems described by the
linear diffusion equation
∂T
c = ∇ · (k∇T ) (235)
∂t
where T is the temperature, c is the specific heat capacity, and k is the conduc-
tivity. The model parameters are summarized in Tab. 53.
Description Linear isotropic elastic material

Record Format IsoHeat num(in) # d(rn) # k(rn) # c(rn) #
- k Conductivity
- c Specific heat capacity
- maturityT0 Baseline value for material maturity, default
0
Supported modes 2dHeat
Table 53: Linear isotropic material for heat transport - summary.
2.2 Isotropic linear material for moisture transport – IsoLin-

Moisture
Linear isotropic material model for moisture transport problems described by
the linear diffusion equation3
∂h
k = ∇ · (c∇h) (236)
∂t
where h is the pore relative humidity (dimensionless, between 0 and 1), k is the
moisture capacity [kg/m3 ], and c is the moisture permeability [kg/m·s]. The
model parameters are summarized in Tab. 54.
2.3 Isotropic material for moisture transport based on

Bažant and Najjar – BazantNajjarMoisture
This is a specific model for nonlinear moisture transport in isotropic cementi-
tious materials, based on [2]. The governing equation
∂h
= ∇ · (C(h)∇h) (237)
∂t
3 Note that the symbols k and c in Eq. (236) have a different meaning than in Eq. (235). The
reason is that the nonlinear model for moisture transport described in Section 2.4 traditionally
uses c for permeability and Eq. (236) should be obtained as a special case of Eq. (239). On
the other hand, the heat conduction model from Section 2.1 was implemented earlier and the
input parameters are directly called k and c, so changing this notation now could lead to
confusion for some older input files.
101
Description Linear isotropic material for moisture transport
Record Format IsoLinMoistureMat num(in) # d(rn) # perm(rn) #
capa(rn) #
- perm moisture permeability
- capa moisture capacity
Table 54: Linear isotropic material for moisture transport - summary.
is a special case of Eq. (239), valid under the assumption that the slope of the
sorption isotherm is linear, i.e. the moisture capacity is constant. In Eq. (237),
h is the relative humidity and C(h) is the humidity-dependent diffusivity ap-
proximated by  
1 − α0
C(h) = C1 α0 + n  (238)
1−h
1 + 1−h c
where C1 is the diffusivity at saturation (typical value for concrete ≈ 30 mm2 /day),
α0 is the dimensionless ratio of diffusivity at low humidity to diffusivity at satu-
ration (typical value ≈ 0.05), hc is the humidity “in the middle” of the transition
between low and high diffusivity (typical value ≈ 0.8), and n is dimensionless
exponent (high values of n, e.g. 12, lead to a rapid transition between low and
high diffusivity). Optionally, it is possible to specify the moisture capacity.
This property is not needed for solution of the diffusion equation (237), but it
is needed if the computed change of relative humidity is transformed into water
content loss (mass of lost water per unit volume).
Description Nonlinear isotropic material for moisture transport

Record Format BazantNajjarMoistureMat num(in) # d(rn) # c1(rn) #
n(rn) # alpha0(rn) # hc(rn) # [ capa(rn) #]
- c1 moisture diffusivity at full saturation [m2 s−1 ]
- n exponent [-]
- alpha0 ratio between minimum and maximum diffusivity
[-]
- hc relative humidity at which the diffusivity is exactly
between its minimum and maximum value [-]
- capa moisture capacity (default value is 1.0)
Table 55: Nonlinear isotropic material for moisture transport - summary.
102
2.4 Nonlinear isotropic material for moisture transport –
NlIsoMoisture
This is a more general model for nonlinear moisture transport in isotropic porous
materials, based on a nonlinear sorption isotherm (relation between the pore
relative humidity h and the water content w) and on a humidity-dependent
moisture permeability. The governing differential equation solves water mass
balance in a unit volume [kg/m3 /s] and reads
∂h dα
k(h) = ∇ · [c(h)∇h] + wn (239)
∂t dt
where k(h) [kg/m3 ] is the humidity-dependent moisture capacity (k(h) = ∂w ∂h
which is derivative of the moisture content w(h) [kg/m3 ] with respect to the
relative humidity), c(h) [kg/m/s] is the moisture permeability and the sink
term wn dαdt corresponds to non-evaporable water loss due to hydration. wn is
non-evaporable water content for complete hydration per m3 [kg/m3 ] and α is
degree of hydration [-]. For the majority of cements, 1 kg of cement consumes
approximately 0.23 kg of non-evaporable water at complete hydration.
So far, six different functions for the sorption isotherm have been imple-
mented (in fact, what matters for the model is not the isotherm itself but its
derivative—the moisture capacity):
1. Linear isotherm (isothermT ype = 0) is characterized only by its slope
given by parameter moistureCapacity.
2. Piecewise linear isotherm (isothermT ype = 1) is defined by two arrays
with the values of pore relative humidity iso h and the corresponding
values of moisture content iso w(h). The arrays must be of the same size.
3. Ricken isotherm [15] (isothermT ype = 2), which is widely used for sorp-
tion of porous building materials. It is expressed by the equation
ln(1 − h)
w(h) = w0 − (240)
d
where w0 [kg/m3 ] is the water content at h = 0 and d [m3 /kg] is an
approximation coefficient. In the input record, only d must be specified
(w0 is not needed). Note that for h = 1 this isotherm gives an infinite
moisture content.
4. Isotherm proposed by Kuenzel [15] (isothermT ype = 3) in the form
(b − 1)h
w(h) = wf (241)
b−h
where wf [kg/m3 ] is the moisture content at free saturation and b is a
dimensionless approximation factor greater than 1.
5. Isotherm proposed by Hansen [13] (isothermT ype = 4) in the form

−1/n
ln h
u(h) = uh 1 − (242)
A
103
characterizes the amount of adsorbed water by the moisture ratio u [kg/kg].
To obtain the moisture content w, it is necessary to multiply the mois-
ture ratio by the density of the solid phase. In (242), uh is the maximum
hygroscopically bound water by adsorption, and A and n are constants
obtained by fitting of experimental data.
6. The BSB isotherm [4] (isothermT ype = 5) is an improved version of the
famous BET isotherm. It is expressed in terms of the moisture ratio
CkVm h
u(h) = (243)
(1 − kh)(1 + (C − 1)kh)
where Vm is the monolayer capacity, and C depends on the absolute tem-

perature T and on the difference between the heat of adsorption and con-
densation. Empirical formulae for estimation of the parameters can be
found in [27]. Note that these formulae hold quite accurately for cement
paste only; a reduction of the moisture ratio is necessary if the isotherm
should be applied for concrete.
The present implementation covers three functions for moisture perme-
ability:
1. Piecewise linear permeability (permeabilityT ype = 0) is defined by two
arrays with the values of pore relative humidity perm h and the corre-
sponding values of moisture content perm c(h). The arrays must be of
the same size.
2. The Bažant-Najjar permeability function (permeabilityT ype = 1) is
given by the same formula (238) as the diffusivity in Section 2.3. All
parameters have a similar meaning as in (238) but c1 is now the moisture
permeability at full saturation [kg/m·s].
3. Permeability function proposed by Xi et al. [27] (permeabilityT ype = 2)
reads h γh (h−1)
i
c(h) = αh + βh 1 − 2−10 (244)
where αh , βh and γh are parameters that can be evaluated using empirical

mixture-based formulae presented in [27]. However, if those formulae are
used outside the range of water-cement ratios for which they were cali-
brated, the permeability may become negative. Also the physical units
are unclear.
Note that the Bajant-Najjar model from Section 2.3 can be obtained as a spe-
cial case of the present model if permeabilityT ype is set to 1 and isothermT ype
is set to 0. The ratio c1/moistureCapacity then corresponds to the diffusivity
parameter C1 from Eq. (237).
2.5 Material for cement hydration - CemhydMat

CemhydMat represents a hydrating material based on CEMHYD3D model ver-
sion 3.0, developed at NIST [5]. The model represents a digital hydrating mi-
crostructure, driven with cellular automata rules and combined with cement
104
chemistry. Ordinary Portland cement is treated without any difficulties, blended
cements are usually decomposed into hydrating Portland contribution and in-
tert secondary cementitious material. The microstructure size can be from
10 × 10 × 10 to over 200 × 200 × 200 µm. For standard computations the
size 50 × 50 × 50 suffices.
Each material instance creates an independent microstructure. It is also
possible to enforce having different microstructures in each integration point.
The hydrating model is coupled with temperature and averaging over shared
elements within one material instance occurs during the solution. Such ap-
proach allows domain partitioning to many CemhydMat instances, depending
on expected accuracy or computational speed. A more detailed description with
engineering examples was published [25]. Tab. 57 summarizes input parameters.
The input XML file specifies the details about cement and concrete com-
position. It is possible to start all simulations from the scratch, i.e. with the
reconstruction of digital microstructure. Alternatively, the digital microstruc-
ture can be provided directly in two files; one for chemical phases, the second for
particle’s IDs. The XML input file can be created with the CemPy package, ob-
tainable from http://mech.fsv.cvut.cz/∼smilauer/index.php?id=software. The
CemPy package alleviates tedious preparation of particle size distribution etc.
The linear solver (specified as NonStationaryProblem) performs well when
the time integration step is small enough (order of minutes) and heat capacity,
conductivity and density remain constant. If not so, use of nonlinear solver is
strongly suggested (specified as NlTransientTransportProblem).
105
Description Nonlinear isotropic material for moisture transport
Record Format NlIsoMoistureMat num(in) # d(rn) # isothermType(in) #
permeabilityType(in) # [ rhodry(rn) #] [ capa(rn) #]
[ iso h(ra) #] [ iso w(h)(ra) #] [ dd(rn) #] [ wf(rn) #] [ b(rn) #]
[ uh(rn) #] [ A(rn) #] [ nn(rn) #] [ c(rn) #] [ k(rn) #] [ Vm(rn) #]
[ perm h(ra) #] [ perm c(h)(ra) #] [ c1(rn) #] [ n(rn) #]
[ alpha0(rn) #] [ hc(rn) #] [ alphah(rn) #] [ betah(rn) #]
[ gammah(rn) #] [ wn(rn) #] [ alpha(expr) #]

- isothermType isotherm function as listed above (0, 1, ...5)
- permeabilityType moisture permeability function as listed
above (0, 1, 2)
- rhodry [kg/m3 ] density of dry material (for
isothermT ype = 4 and 5)
- capa [kg/m3 ] moisture capacity (for isothermT ype = 0)
- iso h [-] humidity array (for isothermT ype = 1)
- iso w(h) [kg/m3 ] moisture content array (for
isothermT ype = 1)
- dd [-] parameter (for isothermT ype = 2)
- wf [kg/m3 ] is the moisture content at free saturation (for
isothermT ype = 3)
- b [-] parameter (for isothermT ype = 3)
- uh [kg/kg] maximum hygroscopically bound water by ad-
sorption (for isothermT ype = 4)
- A [-] parameter (for isothermT ype = 4)
- n [-] parameter (for isothermT ype = 4)
- Vm (for isothermT ype = 5)
- k (for isothermT ype = 5)
- C (for isothermT ype = 5)
- perm h [-] humidity array (for permeabilityT ype = 0)
- perm c(h) [kg m−1 s−1 ] moisture permeability array (for
permeabilityT ype = 0)
- c1 [kg m−1 s−1 ] moisture permeability at full saturation
(for permeabilityT ype = 1)
- n [-] exponent (for permeabilityT ype = 1)
- alpha0 [-] ratio between minimum and maximum diffusiv-
ity (for permeabilityT ype = 1)
- hc [-] relative humidity at which the diffusivity is ex-
actly between its minimum and maximum value (for
permeabilityT ype = 1)
- alphah [kg m−1 s−1 ] (for permeabilityT ype = 2)
- betah [kg m−1 s−1 ] (for permeabilityT ype = 2)
- gammah [-] (for permeabilityT ype = 2)
- wn [kg m−3 ] nonevaporable water content per m3 of con-
crete, default 0.23 kg/kg of cement
- alpha [-] function of degree of hydration
Table 56: Nonlinear isotropic material for moisture transport - summary.

106
Description Cemhyd - hydrating material
Record Format CemhydMat num(in) # d(rn) # k(rn) # c(rn) # file(s) #
[eachGP(in) #] [densityType(in) #] [conductivityType(in) #]
[capacityType(in) #] [castingtime(rn) #] [nowarnings(ia) #]
[scaling(ra) #] [reinforcementDegree(rn) #]
- k Conductivity
- c Specific heat capacity
- file XML input file for cement microstructure and concrete
composition
- eachGP 0 (default) no separate microstructures in each
GP, 1 assign separate microstructures to each GP
- densityType 0 (default) get density from OOFEM input
file, 1 get it from XML input file
- conductivityType 0 (default) get constant conductivity
from OOFEM input file, 1 compute as λ = k(1.33 − 0.33α)
[21]
- capacityType 0 (default) get capacity, 1 according to
Bentz, 2 according to XML and CEMHYD3D routines
- castingtime optional casting time of concrete, from which
hydration takes place. Absolute time is used.
- nowarnings supresses warnings when material data are out
of standard ranges. The array of size 4 represent entries
for density, conductivity, capacity, temperature. Nonzero
values mean supression.
- scaling components in the array scale density, conductiv-
ity, capacity in this order. nowarnings are checked before
scaling.
- reinforcementDegree specifies the area fraction of rein-
forcement. Typical values is 0.015. Steel reinforcement
slightly increases concrete conductivity and slightly de-
creases its capacity. Thermal properties of steel are con-
sidered 20 W/m/K and 500 J/kg/K.
Supported modes 2dHeat, 3dHeat
Table 57: Cemhydmat - summary.
107
2.6 Material for cement hydration - HydratingConcreteMat
Simple hydration models based on chemical affinity are implemented. The mod-
els calculate degree of hydration of cement, α, which can be scaled to the level of
concrete when providing corresponding amount of cement in concrete. Blended
cements can be considered as well, either by separating supplementary cemen-
titious materials from pure Portland clinker or by providing parameters for the
evolution of hydration degree and potential heat. Released heat from cement
paste is obtained from
Q(t) = αQpot , (245)
where potential hydration heat, Qpot , is expressed in kJ/kg of cement and for
pure Portland cement is around 500 kJ/kg.
Evolution of hydration degree under isothemal curing conditions is approx-
imated by several models. Scaling from a reference temperature to arbitrary
temperature is based on Arrhenius equation, which coincides with the maturity
method approach. The equivalent time, te , is defined as time under constant
reference (isothermal) temperature
te (T0 ) = t(T )krate , (246)

Ea 1 1
krate = exp − , (247)
R T0 T
where t is real time, T is the arbitrary constant temperature of hydration, T0 is

a reference temperature, R is the universal gas constant (8.314 Jmol−1 K−1 ) and
Ea is the apparent activation energy. Due to varying history of temperature,
incremental solution is adopted. Linear and nonlinear nonstationary solvers
are supported for all hydrationmodeltype’s. Hydration models are evaluated at
intrinsic time in each time step. Usually, intrinsic time is in the middle of the
time step.
The hydrationmodeltype = 1 is based on exponential approximation of hy-
dration degree [23]. Equivalent time increment is added in each time step. Thus
all the thermal history is stored in the equivalent time
β !
τ
α(te ) = α∞ exp − (248)
te
where three parameters τ , β and α∞ are needed. Some meaningful parameters

are provided in [23], e.g. τ = 26 · 3600 = 93600 s, β = 0.75 and α∞ = 0.90.
The hydrationmodeltype = 2 is inspired by Cervera et al. [6], who proposed
an analytical form of the normalized affinity which was refined in [7]. A slightly
modified formulation is proposed here. The affinity model is formulated for a
reference temperature 25◦ C

dα B2 α
= Ã25 (α)krate = B1 + α (α∞ − α) fs exp −η̄ krate
(249)
dt α∞ α∞
α > DoH1 ⇒ fs = 1 + P1 (α − DoH1 ) else fs = 1 (250)
where B1 , B2 are coefficients to be calibrated, α∞ is the ultimate hydration

degree and η̄ represents microdiffusion of free water through formed hydrates.
108
The function fs adds additional peak which may occur in slag-rich blended
cements with two parameters DoH1 , P1 . The solution proceeds incrementally,
where α is the unknown. During one macroscopic time step, Eq. (249) needs to
be integrated in finer inner steps. This is controlled with two optional variables;
maxmodelintegrationtime specifies maximum integration time in the loop while
minmodeltimestepintegrations specifies minimum number of integration steps.
Description HydratingConcreteMat
Record Format HydratingConcreteMat num(in) # d(rn) # k(rn) # c(rn) #
hydrationmodeltype(in) # Qpot(rn) # masscement(rn) # [ac-
tivationenergy(rn) #] reinforcementdegree(rn) #] [density-
type(in) #] [conductivitytype(in) #] [capacityType(in) #]
[minModelTimeStepIntegrations(in) #] [maxmodelintegra-
tiontime(rn) #] [castingTime(rn) #]
- d material density about 2300 kg/m3
- k Conductivity about 1.7 W/m/K
- c Specific heat capacity about 870 J/kg/K
- hydrationmodeltype 1 is exponential model from Eq. (248),
2 is affinity model from Eq. (249)
- Qpot Potential heat of hydration, about 500 kJ/kg of ce-
ment
- masscement Cement mass per 1m3 of concrete, about 200-
450
- activationenergy Arrhenius activation energy, 38400. (de-
fault)
- DoHinf degree of hydration at infinite time
- B1,B2,eta,DoH1,P1 parameters from Eq. (249)-Eq. (250)
- reinforcementDegree specifies the area fraction of rein-
forcement. Typical values is 0.015. Steel reinforcement
slightly increases concrete conductivity and slightly de-
creases its capacity. Thermal properties of steel are con-
sidered 20 W/m/K and 500 J/kg/K.
- densityType 0 (default)
- conductivityType 0 (default), 1 compute as λ = k(1.33 −
0.33α) [21]
- capacityType 0 (default)
- minModelTimeStepIntegrations Minimum integrations per
time step in affinity model 30 (default)
- maxmodelintegrationtime Maximum integration time step
in affinity model 36000 s (default)
- castingtime optional casting time of concrete, from which
hydration takes place, in s
- scaling components in the array scale density, conductiv-
ity, capacity in this order. nowarnings are checked before
scaling.
Supported modes 2dHeat, 3dHeat
Table 58: HydratingConcreteMat - summary of affinity hydration models.
109
Fig. (15) shows mutual comparison of three hydration models implemented
in OOFEM. Parameters for exponential model according to Eq. (248) are τ =
26 · 3600 = 93600 s, β = 0.75, α∞ = 0.90. Parameters for affinity model
according to Eq. (249) are B1 = 3.519e − 4 s−1 , B2 = 8.0e − 7, η = 7.4,
α∞ = 0.85.
1.0
Degree of hydration [-]
0.8
0.6
0.4
Exponential model
Affinity model
0.2 CEMHYD3D model
Calorimeter
0.0
0.1 1 10 100 1000
Hydration time at isothermal 25°C [days]
Figure 15: Performance of implemented hydration models: exponential model

from Eq. (248), affinity model from Eq. (249), CEMHYD3D model from Sub-
section 2.5.
2.7 Coupled heat and mass transport material model -

HeMotk
Coupled heat and mass transfer material model. Source: T. Krejci doctoral
thesis; Bazant and Najjar, 1972; Pedersen, 1990. Assumptions: water vapor is
the only driving mechanism; relative humidity is from range 0.2 - 0.98 (I and II
regions). The model parameters are summarized in Tab. 59.
110
Description Coupled heat and mass transfer material model
Record Format HeMotk num(in) # d(rn) # a 0(rn) # nn(rn) # phi c(rn) #
delta wet(rn) # w h(rn) # n(rn) # a(rn) # latent(rn) # c(rn) #
rho(rn) # chi eff(rn) # por(rn) # rho gws(rn) #
- d, rho material density
- a 0 constant (obtained from experiments) a 0 [Bazant and
Najjar, 1972]
- nn constant-exponent (obtained from experiments) n
[Bazant and Najjar, 1972]
- phi c constant-relative humidity (obtained from experi-
ments) phi c [Bazant and Najjar, 1972]
- delta wet constant-water vapor permeability (obtained
from experiments) delta wet [Bazant and Najjar, 1972]
- w h constant water content (obtained from experiments)
w h [Pedersen, 1990]
- n constant-exponent (obtained from experiments) n [Ped-
ersen, 1990]
- a constant (obtained from experiments) A [Pedersen,
1990]
- latent latent heat of evaporation
- c thermal capacity
- chi eff effective thermal conductivity
- por porosity
- rho gws saturation volume density
Supported modes 2dHeMo
Table 59: Coupled heat and mass transfer material model - summary.
111
2.8 Coupled heat and mass transport material model -
HeMoKunzel
The presented formulation is based on the work of Kuenzel [15]. The model
is suitable for problems with dominating water diffusion and negligible water
convection. The governing equations for temperature and humidity reads
∂Q ∂Q ∂T ∂T
= = Cv = −∇qT = ∇ (λ∇T ) + hv ∇ (δp ∇(Hpsat ))(251)
∂t ∂T ∂t ∂t
∂w ∂w ∂H
= = −∇qH = ∇ (DH ∇H + δp ∇(Hpsat )) (252)
∂t ∂H ∂t
T (K) Temperature
H (-) Relative humidity 0-1
∂Q 3
∂T ≈ Cv J/K/m Heat storage capacity per volume
∂w 3
∂H kg/m Moisture storage capacity - sorption isotherm
3
Q J/m Total amount of heat in unit volume
2
qT W/m Heat flux
λ W/m/K Thermal conductivity
hv J/kg Evaporation enthalpy of water
δp kg/m/s/Pa Water vapour permeability
psat Pa Water vapour saturation pressure
w kg/m3 Moisture content
DH kg/m/s Liquid conduction coefficient
2
Dw m /s Water diffusivity
Table 60: Parameters from Kunzel’s model.
Numerical solution leads to the system of equations

KT T KT H rT CT T CT H ṙ T qT
+ = , (253)
K HT K HH rH C HT C HH ṙ H qH
where
Z Z
KT T = B T kT T BdΩ, KT H = B T kT h BdΩ, (254)
Ω Ω
Z Z
K HT = B T khT BdΩ, K HH = B T khh BdΩ, (255)
Ω Ω
Z Z
CT T = N T cT T N dΩ, CT H = N T cT h N dΩ, (256)
Ω Ω
Z Z
T
C HT = N chT N dΩ, C HH = N T chh N dΩ, (257)
Ω Ω
112
Z Z
qT = N T q T dΓ, qH = N T q h dΓ, (258)
Γ2 Γ2
where
∆psat
kT T = λ(w) + hv · δp (T ) · H · (T ), (259)
∆T
kT H = hv · δp (T ) · psat (T ), (260)
∆psat
kHT = δp (T ) · H · (T ), (261)
∆T
∆w
kHH = Dw (H) · (H) + δp (T ) · psat (T ), (262)
∆H
cT T = Cs · ρ + Cw · w, (263)
cT H = 0, (264)
cHT = 0, (265)
∆w
cHH = (H). (266)
∆H
Note, that conductivity matrix K is unsymmetric hence unsymmetric matrix
storage needs to be used (smtype).
Description Coupled heat and mass transfer material model

Record Format HeMoKunzel num(in) # d(rn) # iso type(in) # iso wh(rn) #
mu(rn) # permeability type(in) # A(rn) # lambda0(rn) #
b(rn) # cs(rn) # [ pl(rn) #] [ rhoH2O(rn) #] [ cw(rn) #] [ hv(rn) #]
Parameters -num material model number
-d bulk density of dry building material [kg/m3 ]
-iso type=0 is isotherm from Hansen needing iso n, iso a,
=1 is Kunzel which needs iso b
-iso wh maximum adsorbed water content [kg/m3 ]
-mu water vapor diffusion resistance [-]
-permeability type =0 is Multilin h needing perm h,
perm Dw(h), =1 is Multilin wV needs perm wV,
perm DwwV, =2 is Kunzelperm needs A as water ab-
sorption coefficient [kg/m/s0 .5]
-lambda0rn thermal conductivity [W/m/K]
-b thermal conductivity supplement [-]
-cs specific heat capacity of the building material [J/kg/K]
-[pl] ambient air pressure [Pa], default = 101325
-[rhoH2O] water density [kg/m3], default = 1000
-[cw] specific heat capacity of liquid water, default = 4183
-[hv] latent heat of water phase change [J/kg], default =
2.5e+6
Supported modes 2dHeMo, 3dHeMo
Table 61: Coupled heat and mass transfer material model Kunzel - summary.
113
2.9 Material for unsaturated flow in lattice models – Lat-
ticeTransMat
Lattice transport elements have to be used with this material model. A positive
sign is assumed for liquid tension, unlike the convention of soil mechanics which
assumes compression positive.
These transport elements are idealised as one-dimensional conductive pipes.

The gradient of hydraulic head, which governs flow rate along each transport
element, is determined from the capillary pressures Pc at the two nodes.
The mass balance equation describes the change in moisture inside a porous
element as a consequence of liquid flow and solid-liquid retention. It leads to
the following partial differential equation
∂Pc Pc
c + kdiv(∇( − ρz)) = 0 (267)
∂t g
where Pc is the capillary pressure, c is the mass capacity function(s2 m−2 ), k is
the Darcy hydraulic conductivity (ms−1 ), ρ is the fluid mass density, g is the
acceleration of gravity, z is the capillary height and t is the time.
The hydraulic conductivity k consists of

k = k0 + kc (268)
where k0 is the hydraulic conductivity of the intact material and kc is the ad-
ditional conductivity due to cracking.
Darcy hydraulic conductivity k0 is defined as

ρg
k0 = κ (269)
µ
where µ is the dynamic viscosity (P a.s), κ is the permeability also called intrin-
sic conductivity(m2 ), and κr is the relative conductivity. κr is a function of the
effective degree of saturation.
The cracking part is

ρg wc3
kc = ξ κr (270)
µ 12h
where ξ is a tortuosity factor taking into account the roughness of the crack
surface, wc is the equivalent crack opening of the dual mechanical lattice and h
is the length of the dual mechanical element.
2.9.1 One-dimensional transport element

The discrete form of the differential equation for mass transport for a one-
dimensional transport element is
∂Pc
αe Pc + Ce = fe (271)
∂t
where Pc is a vector containing the nodal values of the capillary pressure, αe is
the conductivity matrix, Ce is the capacity matrix and fe is the nodal flow rate
114
vector (kg.m3 ).
The capacity matrix is

Al 2 1
Ce = c (272)
12 1 2
where c is the capacity of the material, l is the length of the transport element
and A is the cross-sectional area of the transport element.
The conductivity matrix is defined as

Ak 1 −1
αe = (273)
l g −1 1
2.9.2 Constitutive laws

The mass transport equation is based on the constitutive laws for the capacity
c and the hydraulic conductivity k.
The capacity c is defined as

∂θ
c = −ρ (274)
∂Pc
where θ is the volumetric water content (θ = VVwT with Vw the volume of wa-
ter, and VT the total volume) which is calculated by a modified version of van
Genuchtens retention model. Note that the presence of a crack in an element
does not influence the capacity in the present model.
The volumetric water content is
θ = Se (θs − θr ) + θr (275)
where θr and θs are the residual and saturated water contents corresponding to
effective saturation values of Se = 0 and Se = 1, respectively.
The effective degree of saturation Se is defined as

(
θm −θr
1 −m
Pc 1−m
Se = θs −θr 1 + a if Pc ≥ Pc(aev) (276)
1 if Pc < Pc(aev)
where θm is an additional model parameter and Pc(aev) is the air-entry value

of capillary pressure which separates saturated (Pc < Pc(aev) ) from unsaturated
states (Pc ≥ Pc(aev) ). It is intuitive that the smaller the pore size of the mate-
rial, the larger the value of Pc(aev) will be.
115
The relative conductivity κr is a function of the effective degree of saturation
and is defined as
 " 1 #m 2
m
S
 1 − 1 − θm −θe
r 
p  θs −θr 
κr = Se 
 " m1 m 
#  (277)
1 − 1 − θm1−θr
 
θs −θr
−θr
If θm = θs , we have θθms −θ r
= 1 : the equation reduces to the expression of the
relative conductivity of the original van Genuchten model.
The model parameters are summarized
Description Material for fluid transport in lattice models

Record Format latticetransmat num(in) # d(rn) # k(rn) # vis(rn) # con-
type(in) # thetas(rn) # thetar(rn) # paev(rn) # m(rn) # a(rn) #
thetam(rn) # [ ctor(rn) #]
- d fluid mass density
- k permeability (m2 )
- vis dynamic viscosity (P a.s)
- contype unsaturated flow allowed when contype=1
- thetas saturated water content
- thetar residual water content
- paev air-entry value of capillary pressure
- m van Genuchten parameter
- a van Genuchten parameter
- thetam additional model parameter for the modified ver-
sion of van Genuchtens retention model
- ctor coefficient of tortuosity (ctor = τ1 ≤ 1)
Supported modes 2dMassLatticeTransport
Table 62: Material for unsaturated flow in lattice models - summary.
116
3 Material Models for Fluid Dynamic
3.1 Newtonian fluid - NewtonianFluid
Constitutive model of Newtonian fluid. The model parameters are summarized
in Tab. 63.
Description Newtonian Fluid material

Record Format NewtonianFluid num(in) # d(rn) # mu(rn) #
- mu viscosity
Supported modes 2d, 3d flow
Table 63: Newtonian Fluid material - summary.
3.2 Bingham fluid - BinghamFluid

Constitutive model of Bingham fluid. This is a constitutive model of non-
Newtonian type. The model parameters are summarized in Tab. 64.
In the Bingham model the flow is characterized by following constitutive
equation
τ = τ 0 + µγ̇ if τ̇ ≥ τ0 (278)
γ̇ = 0 if τ̇ ≤ τ0 (279)
√
where τ is the shear stress applied to material, τ̇ = τ : τ is the shear stress
measure, γ̇ is the shear rate, τ 0 is the yield stress, and µ is the plastic vis-
cosity. The parameters for the model can be in general determined using two
possibilities: (i) stress controlled rheometer, when the stress is applied to mate-
rial and shear rate is measured, and (ii) shear rate controlled rheometer, where
concrete is sheared and stress is measured. However, most of the widely used
tests are unsatisfactory in the sense, that they measure only one parameter.
These one-factor tests include slump test, penetrating rod test, and Ve-Be test.
Recently, some tests providing two parameters on output have been designed
(BTRHEOM, IBB, and BML rheometers). Also a refined version of the standard
slump test has been developed for estimating yield stress and plastic viscosity.
The test is based on measuring the time necessary for the upper surface of the
concrete cone in the slump to fall a distance 100 mm. Semi-empirical models
are then proposed for estimating yield stress and viscosity based on measured
results. The advantage is, that this test does not require any special equipment,
provided that the one for the standard version is available.
In order to avoid numerical difficulties caused by the existence of the sharp
angle in material model response at τ = τ0 , the numerical implementation uses
following smoothed relation for viscosity
τ0
µ = µ0 + (1 − e−mγ̇ ) (280)
γ̇
where m is so called stress growth parameter. The higher value of parameter m,
the closer approximation of the original constitutive equation (278) is obtained.
117
Description Bingham fluid material
Record Format BinghamFluid num(in) # d(rn) # mu0(rn) # tau0(rn) #
- mu0 viscosity
- tau0 Yield stress
Table 64: Bingham Fluid material - summary.
3.3 Two-fluid material - TwoFluidMat

Material coupling the behaviour of two particulat materials based on rule of
mixture. The weighting factor is VOF fraction. The model parameters are
summarized in Tab. 65.
Description Two-Fluid material

Record Format TwoFluidMat num(in) # mat(ia) #
- mat integer array containing two numbers representing
numbers of material models of which the receiver is com-
posed. Material with index 0 is a material, that is fully
active in a cell with VOF=0, material with index 1 is a
material fully active in a cell with VOF=1.
Supported modes 2d flow
Table 65: Two-Fluid material - summary.
3.4 FE2 fluid - FE2FluidMaterial

Constitutive model of multiscale fluids. The macroscale fluid behavior is deter-
mined by a Representative Volume Element (RVE) which is solved for in each
integration point. Some requirements are put on the RVE, such as the it must
use the MixedGradientPressure boundary condition along its outer boundary.
This boundary condition is equivalent to that of Dirichlet boundary condition
in classical homogenization, only adjusted for a mixed control; deviatoric gra-
dient + pressure, instead of the full gradient only. Worth noting is that the
macroscale (and microscale) behavior can still be compressible, in particular is
the RVE contains pores. The triangular and tetrahedral Stokes’ flow elements
both support compressible behavior, but it perhaps be applicable other types
of flow, as long as the extended mixed control option is used when calling the
material routine (the function with gradient and pressure as input). If the RVE
doesn’t contain any pores (such that the response turns out to be incompress-
ible) then material should work for all problem classes. This could for example
be the case of a oil-water mixture. The model parameters are summarized in
Tab. 66.
118
Description FE2 Fluid material
Record Format FE2FluidMaterial num(in) # d(rn) # inputfile(s) #
- d unused
- inputfile input file for RVE problem
Table 66: FE2 fluid material - summary.
119
References
[1] H.D. Baehr and K. Stephan: Heat and Mass Transfer, Springer, 2006.
[2] Z.P. Bažant, L.J. Najjar. Nonlinear water diffusion in nonsaturated con-
crete. Materials and Structures, 5:3–20, 1972.
[3] Z.P. Bažant, J. Planas: Fracture and Size Effect in Concrete and Other
Quasibrittle Materials, CRC Press, 1998.
[4] S. Brunauer, J. Skalny, E.E. Bodor: J. Colloid Interface Sci, 30, 1969.
[5] D.P. Bentz: CEMHYD3D: A Three-Dimensional Cement Hydration and
Microstructure Development Modeling Package. Version 3.0., NIST Build-
ing and Fire Research Laboratory, Gaithersburg, Maryland, Technical re-
port, 2005.
[6] M. Cervera, J. Oliver, and T. Prato: Thermo-chemo-mechanical model for
concrete. I: Hydration and aging. Journal of Engineering Mechanics ASCE,
125(9):1018–1027, 1999.
[7] D. Gawin, F. Pesavento, and B. A. Schrefler: Hygro-thermo-chemo-
mechanical modelling of concrete at early ages and beyond. Part I: Hydra-
tion and hygro-thermal phenomena. International Journal for Numerical
Methods in Engineering, 67(3):299–331, 2006.
[8] P. Grassl and M. Jirásek. Damage-plastic model for concrete failure. In-
ternational Journal of Solids and Structures, 43:7166–7196, 2006.
[9] P. Grassl and M. Jirásek. CDPM2: A damage-plasticity approach to mod-
elling the failure of concrete. International Journal of Solids and Structures,
50:3805–3816, 2013.
[10] P. Grassl and M. Jirásek. Meso-scale approach to modelling the fracture
process zone of concrete subjected to uniaxial tension. International Jour-
nal of Solids and Structures, 47: 957–968, 2010.
[11] P. Grassl, D. Xenos, M. Jirásek and M. Horák. Evaluation of nonlocal ap-
proaches for modelling fracture near nonconvex boundaries. International
Journal of Solids and Structures, 51: 3239-3251, 2014.
[12] M. Jirásek, Z.P. Bažant: Inelastic analysis of structures, John Wiley, 2001.
[13] P.F. Hansen: Coupled Moisture/Heat Transport in Cross Sections of Struc-
tures, Beton og Konstruktionsinstituttet, 1985.
[14] E. Hoek and Z.T. Bieniawski: Brittle Rock Fracture Propagation In Rock
Under Compression, International Journal of Fracture Mechanics 1(3), 137-
155, 1965.
[15] H.M. Künzel, H.M.: Simultaneous heat and moisture transport in building
components, Ph.D. thesis, IRB-Verlag, 1995.
[16] M. Jirásek: Comments on microplane theory, Mechanics of Quasi-Brittle
Materials and Structures, ed. G. Pijaudier-Cabot, Z. Bittnar, and B.
Gérard, Hermès Science Publications, Paris, 1999, pp. 57-77.
120
[17] B. Lourenco, J.G. Rots: Multisurface Interface Model for Analysis of Ma-
sonry Structures, Journal of Engng Mech, vol. 123, No. 7, 1997.
[18] M. Ortiz, E.P. Popov: Accuracy and stability of integration algorithms

for elasto-plastic constitutive relations, Int. J. Numer. Methods Engrg, 21,
1561-1576, 1985.
[19] B. Patzák: OOFEM home page, http://www.oofem.org, 2003.
[20] J.C. Simo, T.J.R. Hughes: Computational Inelasticity, Springer, 1998.
[21] J. Ruiz, A. Schindler, R. Rasmussen, P. Kim, G. Chang: Concrete tem-

perature modeling and strength prediction using maturity concepts in
the FHWA HIPERPAV software, 7th international conference on concrete
pavements, Orlando (FL), USA, 2001.
[22] J.C. Simo, J.G. Kennedy, S. Govindjee: Non-smooth multisurface plastic-

ity and viscoplasticity. Loading/unloading conditions and numerical algo-
rithms, Int. J. Numer. Methods Engrg, 26, 2161-2185, 1988.
[23] A. K. Schindler and K. J. Folliard: Heat of Hydration Models for Cemen-
titious Materials, ACI Materials Journal, 102, 24 - 33, 2005.
[24] J.C. Simo, K.S. Pister: Remarks on rate constitutive equations for finite
deformation problems: computational implications, Comp Methods in Ap-
plied Mech and Engng, 46, 201-215, 1984.
[25] V. Šmilauer and T. Krejčı́, Multiscale Model for Temperature Distribution
in Hydrating Concrete, International Journal for Multiscale Computational
Engineering, 7 (2), 135-151, 2009.
[26] J.H.P. de Vree, W.A.M. Brekelmans, and M.A.J. van Gils: Comparison
of nonlocal approaches in continuum damage mechanics. Computers and
Structures 55(4), 581588, 1995.
[27] Y. Xi, Z.P. Bažant, H.M. Jennings: Moisture Diffusion in Cementitious

Materials, Advn Cem Bas Mat, 1994.
[28] Z.P. Bažant, S. Prasannan: Solidification theory for concrete creep. I: For-
mulation. Journal of Engineering Mechanics 115(8), 16911703, 1989.
[29] Z.P. Bažant, A.B. Hauggaard, F. Ulm: Microprestress-solidification the-
ory for concrete creep. I: Aging and drying effects. Journal of Engineering
Mechanics 123(11), 11881194, 1997.
[30] M. Jirásek, P. Havlásek: Microprestress-Solidification Theory of Concrete
Creep: Reformulation and Improvement. Cement and Concrete Research
60, 5162, 2014.
[31] International Federation for Structural Concrete (fib): The fib Model Code
for Concrete Structures 2010.
121

OOFEM

Uploaded by

Copyright:

Available Formats

OOFEM

Uploaded by

Document Information

Original Description:

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

OOFEM

Uploaded by

Copyright:

Available Formats

Material Model Library Manual

Czech Technical University

2 Material Models for Transport Problems 98

3 Material Models for Fluid Dynamic 114

Description Linear isotropic elastic material

Table 1: Linear Isotropic Material - summary.

1.1.2 Orthotropic linear elastic material - OrthoLE

dxx = EX (1 − νyz ∗ νzy )/ξ, (3)

Ei is Young’s modulus in the i-th direction, Gij is the shear modulus in

νxy Ey = νyx Ex , νyz Ez = νzy Ey , νzx Ex = νxz Ez

The model description and parameters are summarized in Tab. 2.

1.1.3 General anisotropic linear elastic material - AnisoLE

1.1.4 Hyperelastic material - HyperMat

Table 2: Orthotropic, linear elastic material – summary.

Description Anisotropic, linear elastic material

Table 3: Anisotropic, linear elastic material – summary.

1.1.5 Hyperelastic material - Compressible Mooney-Rivlin

Table 4: Hyperelastic material - summary.

1.2 Winkler-Pasternak model

Table 5: Compressible Mooney-Rivlin - summary.

Description Winkler Pasternak isotropic material for subsoil interaction

Table 6: Winkler Pasternak material - summary.

1.3 Large-strain master material

1.4 Plasticity-based material models

Table 7: Large-strain master material material - summary.

with the pressure-dependent equivalent stress

The flow rule is derived from the plastic potential

where αψ is the dilatancy coefficient. An associated model with αψ = α would

λ̇ ≥ 0, f (σ, τY ) ≤ 0, λ̇ f (σ, τY ) = 0 (28)

portional to the plastic multiplier. For α = αψ = 0, the associated J2 -plasticity

h(κ) = τlimit + (τ0 − τlimit ) e−κ/κc (31)

can be used for a more realistic description of hardening.

1.4.2 Drucker-Prager model with tension cut-off and isotropic dam-

σ = (1 − ω)σ̄ = (1 − ω)D : (ε − εp ), (33)

Table 8: DP material - summary.

the definition of the yield function in terms of the effective stress,

the flow rule p

ω(κ) = ωc (1 − e−aκ ), (37)

where ωc is critical damage and a is a positive dimensionless parameter. More

Table 9: Drucker Prager material with tension cut-off - summary.

Small-strain formulation: The small-strain version of hardening Mises plas-

the stress strain law

σ = (1 − ω)σ̄ = (1 − ω)D : (ε − εp ), (39)

the definition of the yield function in terms of the effective stress,

the linear hardening law

λ̇ > 0 f (s̄, κ) ≤ 0 λ̇f (s̄, κ) = 0. (44)

and the damage law

Description Mises plasticity model with isotropic hardening

Table 10: Mises plasticity – summary.

1.4.4 Mises plasticity model with isotropic damage, nonlocal - Mis-

Integral nonlocal formulation One possible regularization technique is based

is a nonnegative function, for r < R monotonically decreasing with increasing

Implicit gradient formulation The gradient formulation can be conceived

- scalingType selects the type of scaling of the weight func-

Supported modes 1dMat, PlaneStrain, 3dMat

Table 11: Nonlocal integral Mises plasticity – summary.

Table 12: Gradient-enhanced Mises plasticity – summary.

1.4.6 Perfectly plastic material with Mises yield condition - Steel1

1.4.7 Composite plasticity model for masonry - Masonry02