Stochastic Computation of Quasibrittle Fracture
A DISSERTATION
SUBMITTED TO THE FACULTY OF THE
UNIVERSITY OF MINNESOTA
BY
Anna Gorgogianni
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
Jia-Liang Le, Advisor
August 2021
c© Anna Gorgogianni 2021
Acknowledgements
Reflecting on these challenging and rewarding years of my doctoral studies at the University of
Minnesota, I would like to firstly give my gratitude to my advisor, Professor Jia-Liang Le, for
the exciting research project he assigned to me, for exposing me to areas of knowledge I had not
previously explored, for his guidance and commitment to the progress of my research and of myself
as a researcher. Besides, I would like to express my appreciation to Professor Jan Eliáš at Brno
University of Technology for our close collaboration throughout these years, for his efforts, inspiring
enthusiasm and constructive criticism.
I thank Professor Henryk K. Stolarski, Professor Emmanuel Detournay and Professor Ellad B.
Tadmor at the University of Minnesota, for serving in my Ph.D. committee. I indeed consider
myself extremely lucky having been a student of them; I admire Professor Stolarski’s ability to
present tough concepts in a simple and meaningful way, Professor Detournay’s physical intuition
and critical thinking, Professor Tadmor’s clarity in teaching and instructional feedback that allows
students to find the answers themselves.
I deeply appreciate the friendship and support of other members of the research group, particu-
larly Dr. Anu Tripathi and Dr. Teng Man, as well as Mr. Chen Hu and Mr. Tianhao Yan. Many
thanks to Dr. Zhifeng Xu and Mr. Joshua Vievering for our collaboration.
I also acknowledge the financial support for my doctoral studies provided by the U.S. Army
Research Office, the U.S. Department of Energy, and the Department of Civil, Environmental and
Geo-Engineering (CEGE) of the University of Minnesota. I thank the Department Head, Professor
Joseph F. Labuz, as well as all faculty members at CEGE for the welcoming academic environment
they have created.
I strongly feel the need to express my gratitude to my former advisors in Greece, Professor Sofia
Papargyri-Beskou at the Aristotle University of Thessaloniki and Professor Dimitri E. Beskos at
the University of Patras. I thank them for guiding me during my very first research experience,
for further solidifying my interest to continue with advanced studies, and for holding me up during
i
ii
personal hardship. I thank all my former educators for they have largely shaped my thinking. I am
also grateful to Dino Xykis, former Ph.D. graduate of the CEGE Department at the University of
Minnesota, for his help during my arrival in the U.S. and encouragement in my doctoral studies.
I am more than thankful to my parents, Ioannis and Evangelia, for their devotion to my well-
being and personal development which has been the driving force for all my efforts. I also thank
my sister Eleni, my aunts Chariklia and Maria, and my precious friends in Greece for their love and
support that reveals the bright side of life.
To my father, for shaping my character
and guiding my thought.
iii
“If one way be better than another, that you may
be sure is nature’s way.”
Aristotle
iv
ABSTRACT
Continuum finite element (FE) modeling of quasibrittle fracture is a challenging task. It is well
known that the continuum approaches suffer from the issue of spurious mesh sensitivity due to
localization instability. The issue of mesh dependence must be addressed to ensure the robustness of
FE simulations. Numerical techniques to mitigate the spurious mesh sensitivity include localization
limiters such as nonlocal continuum models, as well as energy regularization methods such as the
crack band model. Nonlocal integral formulations and implicit gradient models, have been shown to
be largely equivalent to each other, able to recover mesh-objective failure patterns that preserve the
finite size of localization bands as well as mesh-objective global structural response. A drawback of
these techniques is the need for sufficiently small mesh sizes to capture the nonlocality or resolve the
gradient terms, which makes them inefficient for the analysis of large-scale structures. On the other
hand, the crack band model addresses solely the mesh dependence of the global structural response
only for cases of failure with severe strain localization. It is computationally efficient for large-scale
simulations since it allows use of elements of size larger than the material-dependent localization
bandwidth.
The first part of this dissertation focuses on mitigating the mesh dependence of FE simulations of
softening continua in the deterministic setting, by proposing a general energy regularization scheme.
In contrast to the crack band model, the proposed model captures the effect of the failure mechanism
on the energy regularization of the constitutive law. Considering that each finite element represents
a material element of finite size, the constitutive law assigned to it should generally depend on both
the material properties as well as the failure pattern of the corresponding element. The goal of energy
regularization methods is to ensure correct amount of energy dissipation in simulations of failure. In
cases of a purely localized damage pattern, the post-peak part of the input constitutive law must be
adjusted according to the element size such that the fracture energy is preserved. However, in cases
of diffused damage, the input constitutive model represents purely material behavior and should thus
be independent of the element size. In order to mathematically describe the actual damage pattern,
which evolves during the loading process, the present study introduces a localization parameter which
v
vi
is calculated from the spatial distribution of the damage variable and is updated at every loading
step. The localization parameter governs the energy regularization of the constitutive law, and can
capture the transition from diffused to localized damage during the failure process. The method
is cast into an isotropic damage model, and is further extended to rate-dependent behavior. The
proposed model is applied to FE simulations of static and dynamic failure of ceramic specimens.
The results show that the model is highly effective in mitigating mesh dependence of the global
response of quasibrittle structures experiencing different failure mechanisms.
An important feature of the aforementioned mechanism-based energy regularization scheme is
that it does not require any a priori knowledge of the structural failure mechanism. In the second
part of this dissertation, the model is extended to stochastic FE simulations of quasibrittle structures
with spatially varying element tensile strengths and fracture energies. Focus is placed on the case
when the finite element mesh size is larger than the crack band width of the material, as well as the
autocorrelation length of the underlying random fields of constitutive properties. For general loading
scenarios it is not possible to know the evolving failure mechanism. Therefore, use of pre-assumed
sampling distribution functions of the random element strengths following some well-established
models such as the weakest-link or fiber-bundle model, which is a common practice in the literature,
could be inconsistent with the actual failure mechanism. It is well known that the probability
distributions of fracture properties of a material element depend on the failure mechanism. For
instance, the weakest-link model is applicable for failures with damage localization, whereas the fiber-
bundle model is suitable for failures with distributed cracking. In this study, a general probabilistic
model of fracture properties of finite elements is proposed. The probability distribution of the
fracture property of strength depends on the evolving failure mechanism and can transition from a
fiber-bundle model representation of structural failure to a weakest-link model as damage localization
emerges from diffusive damage. This is achieved by introducing a dependence of the sampling
distribution functions of apparent fracture properties on numerical parameters describing the damage
pattern and introduced in the energy regularization scheme (e.g. localization parameter). Hence,
the size effects on the statistics of fracture properties arising due to the failure mechanism of an
element are taken into account. The model is applied to stochastic FE simulations of quasibrittle
structures of different geometries, subjected to different loading configurations. The simulation
results show that the mesh dependence of the output probabilistic structural response is effectively
mitigated. Most importantly, the mechanism-based energy regularization scheme and the proposed
probabilistic model are interconnected and are always consistent with the actual failure mechanism
of the structure, resulting in a unified and general framework for addressing the mesh dependence
vii
of numerical simulations of quasibrittle failure in both deterministic and stochastic settings.
Contents
List of Figures x
1 Introduction 1
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Numerical Modeling of Quasibrittle Failure . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.1 Cohesive Crack Model and Discrete Element Models . . . . . . . . . . . . . . 4
1.2.2 Crack Band Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.2.3 Nonlocal Damage Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.2.4 Phase-Field Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.3 Probabilistic Modeling of Strength Distribution of Quasibrittle Structures . . . . . . 16
1.3.1 Chain Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.3.2 Fiber-Bundle Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.3.3 Fishnet Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
1.4 Research Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2 Energy Regularization in Deterministic Numerical Simulations of Strain-Softening
Continua 24
2.1 Model Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.2 Extension to Rate-Dependent Behavior . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.3 Simulations of Quasi-Static Failure . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.4 Simulations of the Dynamic Spalling Test . . . . . . . . . . . . . . . . . . . . . . . . 40
2.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3 Rate and Size Effects on Strength Distribution of Quasibrittle Structures 46
3.1 Stochastic Discrete Model Simulations of Dynamic Failure . . . . . . . . . . . . . . . 47
viii
CONTENTS ix
3.1.1 Model Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.1.2 Simulations of Dynamic Tensile Failure of Aluminum Nitride Specimens . . . 52
3.2 Rate-Dependent Finite Weakest Link Model . . . . . . . . . . . . . . . . . . . . . . . 56
3.2.1 Model Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
3.2.2 Numerical Verification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
3.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4 Stochastic Numerical Simulations of Strain-Softening Continua 69
4.1 Model Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.1.1 Probabilistic Modeling of Element Fracture Properties . . . . . . . . . . . . . 70
4.1.2 Description of the Local Damage Pattern . . . . . . . . . . . . . . . . . . . . 77
4.2 Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.2.1 Description of Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.2.2 Unnotched beams under three-point bending . . . . . . . . . . . . . . . . . . 83
4.2.3 Notched beam under three-point bending . . . . . . . . . . . . . . . . . . . . 87
4.2.4 Unnotched beam under four-point bending . . . . . . . . . . . . . . . . . . . 90
4.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
5 Conclusions 94
Bibliography 97
List of Figures
1.1 Probabilistic analysis of failure of a strain softening bar under uniaxial tension . . . 3
1.2 Characterization of fracture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Cohesive crack model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4 Cohesive laws . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.5 a) Crack band parallel to the element sides in a square mesh, and b) Crack band
along the element diagonals in a square mesh . . . . . . . . . . . . . . . . . . . . . . 10
1.6 a) Nonlocal weighting function, and b) treatment of the weighting function near the
structural boundaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.7 A solid body Ω with the crack set Γ: a) sharp cracks, and b) approximated diffuse
crack bands . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.8 Probabilistic models of structural failure: a) chain model, b) fiber-bundle model, and
c) fishnet model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.9 a) Hierarchical model of the RVE, b) idealization of the force transmission mechanisms
inside the FPZ, and c) stress-strain diagrams of the coupled elements . . . . . . . . . 18
2.1 Two distinct crack patterns in a material element represented by a finite element: a)
localized cracking, and b) diffused cracking . . . . . . . . . . . . . . . . . . . . . . . 27
2.2 Schematic of finite element mesh for defining the localization parameter . . . . . . . 28
2.3 Uniaxial tensile stress-strain curves for different mesh sizes for the case of fully local-
ized damage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.4 Rate-dependent uniaxial tensile stress-strain behavior of a material element for the
case of fully diffused damage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.5 Rate effect on material properties: a) apparent tensile strength, b) fracture energy,
and c) Irwin’s characteristic length . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
x
LIST OF FIGURES xi
2.6 Notched and unnotched three-point bend beams considered in the simulation (all
dimensions are in mm) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.7 Evolution of the damage pattern of the notched beam . . . . . . . . . . . . . . . . . 36
2.8 Simulation results of the notched beam: a) load-displacement curves, and b) profile
of normal stress along the ligament for the diffused damage model . . . . . . . . . . 37
2.9 Snapshots of the damage pattern of the unnotched beam during the loading process 38
2.10 Simulated load-displacement curves of the unnotched beam . . . . . . . . . . . . . . 39
2.11 a) Schematic of the modified split Hopkinson bar test, b) stress wave traveling in the
bar, c) stress pulse used in the present simulations, and d) discretization of the bar . 40
2.12 Time histories of free-end velocity predicted by the localized damage model (first
row), the diffused damage model (middle row), and the present model (bottom row) 41
2.13 Time histories of energy dissipation simulated by the localized damage model (first
row), the diffused damage model (middle row), and the present model (bottom row) 43
2.14 Snapshots of the damage pattern of the bar at different times simulated by the local-
ized damage model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
2.15 Snapshots of the damage pattern of the bar at different times simulated by the present
model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.1 Representation of the discrete element model: a) domain discretization, and b) Voronoi
body and facet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.2 A typical realization of the random field h(x) . . . . . . . . . . . . . . . . . . . . . . 51
3.3 A rectangular specimen loaded by a prescribed strain rate: a) schematic of the spec-
imen, and b) weakest link model representation . . . . . . . . . . . . . . . . . . . . . 52
3.4 Simulated average nominal stress-strain curves for different specimen sizes and strain
rates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.5 Simulated damage patterns of specimen of D = 400µm at the peak load . . . . . . . 55
3.6 Modeling of strength distribution of a RVE: a) hierarchical model of static RVE, and
b) fiber-bundle model of RVE under high strain-rate loading . . . . . . . . . . . . . . 57
3.7 Simulated mean size effect curves of nominal strength at different strain rates and the
optimum fits by the rate-dependent finite weakest-link model . . . . . . . . . . . . . 63
3.8 Effect of the strain rate on the RVE size . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.9 Rate dependence of the statistical parameters of P1(σN ) . . . . . . . . . . . . . . . . 65
LIST OF FIGURES xii
3.10 Simulated size effects on the standard deviation of nominal strength at different strain
rates and its comparison with the rate-dependent finite weakest-link model . . . . . 67
4.1 A finite element with a damage zone consisting of nb number of crack bands . . . . . 71
4.2 Schematic of finite element mesh for calculation of the localization and propagation
parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.3 Geometry of unnotched beams under three-point bending considered in the simulation
(all dimensions are in mm) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.4 Mean (left) and standard deviation (right) of peak load of the unnotched beam under
three-point bending with probabilistic models I)-III) (Gf0 = 50J/m2) . . . . . . . . . 84
4.5 Mean (left) and standard deviation (right) of peak load of the unnotched beam under
three-point bending with probabilistic models I)-III) (Gf0 = 15J/m2) . . . . . . . . . 87
4.6 Notched beam under three-point bending considered in the simulation (all dimensions
are in mm) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.7 Mean (left) and standard deviation (right) of peak load of the notched beam under
three-point bending with probabilistic models I)-III) . . . . . . . . . . . . . . . . . . 90
4.8 Unnotched beam under four-point bending considered in the simulation (all dimen-
sions are in mm) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
4.9 Mean (left) and standard deviation (right) of peak load of the unnotched beam under
four-point bending with probabilistic models I)-III) . . . . . . . . . . . . . . . . . . . 92
Chapter 1
Introduction
1.1 Background
Failure of engineering structures made of brittle heterogeneous, aka quasibrittle materials, is caused
by material fracture, which is commonly referred to as quasibrittle fracture. A salient feature of
quasibrittle fracture is the formation of a fracture process zone (FPZ) whose size is not negligible
compared to the characteristic structural dimensions. The failure behavior of quasibrittle structures
is dependent on the ratio between the structure size and the FPZ size, and transitions from quasi-
plastic to perfectly brittle with increasing values of this ratio. At the scale of normal laboratory
testing, typical examples of quasibrittle materials include concrete, composites, ceramics, rock, cold
asphalt, etc.
Quasibrittle materials exhibit a strain-softening constitutive behavior, defined as decreasing stress
at increasing strain in the one-dimensional setting, which may result in strain localization. Strain
localization can by considered as an instability of the constitutive description of the material, which
may allow the emergence of discontinuities in the deformation field, with the formation of a band
of non-uniform deformation and continuing equilibrium, whereas outside the localization band the
initial conditions of homogeneous deformation continue to be satisfied [124]. The localization in-
stability was initially studied for finite strain plasticity and its onset was analyzed by the acoustic
tensor [124]. It was found that strain localization in rock masses described by von Mises plasticity
could occur at either positive or negative values of the hardening modulus, depending on the details
of the incremental constitutive relation [124]. The analysis was limited to the behavior of a material
point in an infinite space, and the stability of the initial state of uniforms strains against localiza-
1
1.1. BACKGROUND 2
tion was not examined. The above considerations were taken into account in [20–22,82], where the
condition for the onset of localization was extended to arbitrary constitutive laws and to bodies of
finite size. The effects of the size of the localization band and of unloading of the material outside
the localization region were also taken into account. Besides, the stability of inelastic structures was
studied within the framework of thermodynamics upon introducing the concept of the tangentially
equivalent elastic structure in [12]. It was shown that [12,22], out of the infinitely many equilibrium
paths emanating from a bifurcation point in the load-displacement diagram of an inelastic structure
with homogeneous properties, the one that will be actually followed by the structure, as dictated by
the 2nd Law of Thermodynamics, is not the path corresponding to homogeneous deformation, but
rather the path with strains localizing inside the narrowest possible region, the localization band,
whereas the surrounding structure experiences unloading.
Based on the preceding conclusion, in finite element (FE) simulations of softening continua,
strain localization would occur in a single finite element, or even in a subdivision of the element
defined by its Gauss points [79]. Thus, upon mesh refinement, the overall volume of the localization
band tends to zero. This implies that the energy expended to form a damage band tends to zero.
Consequently, a smaller mesh size will lead to a more brittle failure behavior of the simulated
structure, which is known as the spurious mesh sensitivity of FE simulations. Various methods,
known as the localization limiters, have been proposed to address this problem by preventing the
strains from localizing into infinitesimally small regions. These methods will be briefly discussed in
Section 1.2.
Another issue that has received very limited attention is the spurious mesh sensitivity of stochas-
tic FE simulations of quasibrittle fracture. This can be demonstrated by considering a quasibrittle
bar under quasi-static uniaxial tension, which is sub-divided into several elements along its length
(Fig. 1.1a)). Due to the strain-softening behavior, strains will localize into one element and the max-
imum load capacity of the structure will be determined by the strength of this element. Considering
the tensile strength ft and fracture energy Gf as the constitutive parameters governing the failure
behavior and thus modeling those only as statistically independent random variables, the peak load
of the bar Pmax is governed by the minimum of the random strengths of the elements in the mesh,
i.e. Pmax = A0min (fti), where A0 is the cross sectional area and fti denotes the tensile strength
of the ith element. Without introducing any mesh-size dependence of the input distributions of the
constitutive parameters, a finer mesh will render a lower mean peak load, and the probability dis-
tribution of Pmax will depend on the number of elements. The existing localization limiters cannot
address the mesh sensitivity of stochastic simulations, since they were formulated within a purely
1.1. BACKGROUND 3
deterministic framework. To demonstrate this, the cumulative distribution function (cdf) of Pmax
was calculated using the conventional crack band model [17], with some commonly employed as-
sumptions: 1) the material has a linear softening stress-strain behavior (Fig. 1.1b)), 2) the material
strength and fracture energy of each element are statistically independent random variables, and 3)
the cdf of ft of each element follows a Gauss-Weibull grafted distribution [14, 90]. Fig. 1.1c) shows
the computed cdfs of Pmax for different element sizes. It can be seen that the resulting cdf of the
peak load strongly depends on the element size. On the other hand, the cdf of the total energy
dissipation of the bar is independent of the element size (Fig. 1.1d)) since the crack band model was
used to regularize the fracture energy.
3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0
maximum load [kN]
0.0
0.2
0.4
0.6
0.8
1.0
el. size 50 mm
el. size 100 mm
el. size 200 mm
40 60 80 100 120 140
dissipated energy [N/m]
0.0
0.2
0.4
0.6
0.8
1.0
el. size 50 mm
el. size 100 mm
el. size 200 mm
p
ro
b
ab
il
it
y
1000 mm

✏
f 0t
F
Element descritization 
Stress-strain relationship 
a) b) 
c) d) 
, P
Figure 1.1: Probabilistic analysis of failure of a strain softening bar under uniaxial tension
As has been demonstrated, strain localization due to softening can lead to the spurious mesh
sensitivity of FE simulations of quasibrittle fracture, in both the deterministic and stochastic settings.
Quasibrittle structures can exhibit a complex damage pattern, featured by both diffused and localized
damage, depending on the loading configuration and structural geometry. For instance, strain-
rate effects have been shown to significantly affect the failure mechanism of quasibrittle structures
[54, 68, 91]. The prevalent damage pattern has profound implications on the regularization of the
constitutive law of finite elements. The goal of this research is to develop a general framework, that
1.2. NUMERICAL MODELING OF QUASIBRITTLE FAILURE 4
is able to take into account the evolving state of damage in a structure and address the issue of mesh
dependence under general failure patterns, for both cases of deterministic and stochastic analyses.
1.2 Numerical Modeling of Quasibrittle Failure
1.2.1 Cohesive Crack Model and Discrete Element Models
Fracture is accompanied by the formation and propagation of a macroscopic crack. At the crack
tip, on account of the stress concentration, irreversible deformation occurs and a nonlinear zone
is formed. This nonlinear zone could consist of the fracture process zone (FPZ), which exhibits
softening damage due to microcracking (and is denoted with black color in Fig. 1.2), and a plastic
(or yielding) zone surrounding the FPZ, which can exhibit perfect or hardening plasticity due to
small scale yielding and void growth. Depending on the dimensions of the FPZ compared to the
structure size, fracture can be characterized as brittle, ductile or quasibrittle. In brittle fracture,
the entire nonlinear zone, and therefore the FPZ as well, is very small compared to the structural
dimensions and can be lumped into a single point. As a consequence, the entire volume of the
structure remains elastic and the theory of linear elastic fracture mechanics (LEFM) is applicable.
In ductile fracture, observed in elastoplastic materials such as metals, the nonlinear plastic zone is
of significant size, whereas the FPZ remains negligible. In quasibrittle fracture, which is the focus
of this dissertation, it is the size of the FPZ which is non-negligible compared to the structural
dimensions, whereas the nonlinear yielding zone is very small, almost pointwise, so that it can be
considered that the boundary of the FPZ coincides with that of the plastic zone.
In the numerical modeling of quasibrittle fracture, an explicit representation is required for the
deformation mechanisms inside the FPZ. Cohesive zone models (CZM), introduced by Barenblatt [6]
and Dugdale [57] for elastic-plastic fracture in ductile metals, and for quasi-brittle materials by
Hillerborg et al. [74] in the so-called fictitious crack model, model the FPZ as a line (Fig. 1.3). The
damage behavior of the FPZ is described by a softening cohesive law, which relates the cohesive
traction to the opening displacement w (Fig. 1.4). The cohesive law is regarded as a material
property (Eq. (1.1)):
σyy = f(w) (1.1)
Eq. (1.1) can be considered as an accurate description of the constitutive behavior of the FPZ in
the case when the normal stresses σxx and σzz along the crack plane are negligible compared to the
compressive and tensile material strengths, otherwise a tensorial constitutive law would be needed
1.2. NUMERICAL MODELING OF QUASIBRITTLE FAILURE 5
Linear elastic
Nonlinear 
hardening
Softening
Linear elastic
Nonlinear 
hardening
Softening
Linear elastic
Nonlinear 
hardening
Softening
Brittle QuasibrittleDuctile
Characterization of fracture 1
Figure 1.2: Characterization of fracture
(Sections 1.2.2, 1.2.3). Two of the most important parameters of the cohesive law described by
Eq. (1.1) are the material tensile strength f(0) = ft, which is the stress level at which the cohesive
opening is activated, and the fracture energy GF , defined as the energy required to create a unit
of new surface area along the crack line. It can be shown that GF is equal to the area under the
softening curve:
GF =
∫ wf
0
f(w)dw (1.2)
where wf is a critical value of the opening of the cohesive crack, after which the stress at the
corresponding location reduces to zero. Thus, the point at which w = wf and σyy = 0 corresponds
to the tip of the real traction-free crack.
Fig. 1.4 shows examples of some commonly used decohesion relations: a) a simple linear relation,
b) a stress-separation relation with a horizontal yield plateau applicable for ductile fracture, and
c) a bilinear curve, which has been used as a realistic cohesive law for concrete by [29, 52]. The
cohesive stress σyy decays from the finite failure strength ft to zero, removing the LEFM stress
singularity at the crack tip. For the case of brittle fracture, LEFM is applicable and can provide
a good approximation to the predictions of the cohesive zone model. In brittle fracture, the most
important property of the cohesive law is the material fracture energy GF , since it fully controls
the peak load of the structure, irrespectively of the assumed shape of the softening curve. For
quasibrittle fracture, the precise shape of the stress separation law plays a significant role in the
proximity of the predicted mechanical response of the structure to experimental results [43].
1.2. NUMERICAL MODELING OF QUASIBRITTLE FAILURE 6
4
FPZ
FPZ and cohesive crack approximation
Cohesive crack 
approximation of the 
FPZ
Figure 1.3: Cohesive crack model
Cohesive laws 2
a) b) c)
Figure 1.4: Cohesive laws
Given that the material properties of fracture toughness GF and failure strength ft introduced
in cohesive laws are independent and have different dimensions, the so-called Irwin’s internal length
lch is usually introduced to characterize the length of the fracture process zone (FPZ) ahead of the
crack tip:
lch =
EGF
f2t
(1.3)
For a vanishing internal length, lch → 0, Griffith’s brittle fracture is asymptotically recovered [66].
Therefore, Barenblatt’s CZM can be used not only for cohesive fracture, but also for brittle fracture.
Another major advantage of the cohesive crack model is that it allows for the entire volume
of the structure to be treated as elastic since all the nonlinearity of the FPZ is embedded in the
boundary conditions at the crack faces. One important limitation of cohesive zone models, though,
is that they cannot account for triaxiality of the stress field ahead of the crack tip, which can lead
to premature fracture through its effect on the fracture toughness of the material. Besides, cohesive
crack models are computationally expensive, requiring topological changes in the mesh to account for
1.2. NUMERICAL MODELING OF QUASIBRITTLE FAILURE 7
the displacement discontinuities due to crack propagation, and the predicted crack paths are mesh
dependent, because fracture is considered to occur only along the interface zero-thickness cohesive
elements.
The cohesive crack models described above constitute numerical tools for modeling fracture
within the discrete approach, since they allow discontinuities in the displacement field (i.e. strong
discontinuities). Other models of softening damage and fracture of quasibrittle materials belonging
to the same approach, are the discrete models in the form of either lattice or particle models, which
consist of a collection of interconnected rigid bodies organized into a net structure. The first lattice
model was introduced by Hrennikoff [75] for the solution of elasticity problems, and in the 1970s
and 1980s lattice models were employed for the study of the fracture process of disordered media in
theoretical physics [61, 73, 123, 125, 139]. Further developments of the lattice discrete model for 2D
and 3D simulations of quasibrittle fracture were made by [28,47,48,80,81,127,141].
The discrete particle model was originated from the so-called distinct element model [45, 46],
which was developed for capturing the behavior of particulate materials such as cohesionless soils and
rock masses. The model was further extended to concrete materials [141], in which the computational
domain consisted of a set of rigid polyhedral elements. In further extensions of this approach, the
connections between elements were implemented by nonlinear springs of zero length [36,37,85]. An
important feature of discrete element models is the introduction of a characteristic length scale (e.g.,
particle size or lattice size), which prevents fracture from localizing into infinitesimally small regions
and thus ensures a finite amount of dissipated energy in the fracture process.
Recently, significant advances have been made towards realistic discrete models for various con-
crete materials. The lattice discrete particle model (LDPM) of Cusatis and coworkers [47, 48, 51],
which considered the random size distribution of concrete heterogeneities, has been experimentally
validated under many loading scenarios, and was shown to give very realistic predictions of cracking
damage and fracture. Discrete element models can also account for the randomness of mesoscale con-
stitutive properties, and they have been successfully applied in stochastic simulations of quasibrittle
fracture [59,69].
In general, discrete models have two major advantages: 1) they are based on a discontinuous
formulation, which avoids the singularity-related issues of continuum-based numerical simulation
methods, such as the finite element method (FEM), and 2) they allow for explicit representation
of the microstructure of the material [114]. Owing to the second advantage, lattice or particle
models are dominant in the simulation of the mesoscale or microscale of materials, where the effects
of material heterogeneity should be considered. A major disadvantage of discrete element models
1.2. NUMERICAL MODELING OF QUASIBRITTLE FAILURE 8
though is their high computational cost, which hampers their application for the simulation of large-
scale structures. For this reason, recent research efforts have focused on the development of robust
multiscale computational frameworks, which aim to bridge the discrete with the continuum approach
to fracture [120].
1.2.2 Crack Band Model
The crack band model [17] is a generalization of the discrete cohesive zone concept to a continuum
formulation by smearing the inelastic strains due to microcracking over the width of the FPZ. The
average inelastic strain ′′yy inside the FPZ is defined from the relationship w(x) =
∫ h
0
′′yy(x, y)dy,
where w(x) is the cohesive crack opening and h the crack band width (Fig. 1.3). Assuming that
the inelastic damage strain is uniformly distributed across the crack band width, its pointwise value
′′yy(x, y) coincides with the average value and thus we have ′′yy = ′′yy = w/h. This assumption also
implies that the fracture energy under mode I loading (opening mode) is given by:
GF = h
∫ f
0
σyyd′′yy (1.4)
where f = wf/h is the ultimate fracturing strain after which the stress component normal to the
localization plane decreases to zero. In contrast to the cohesive crack model, the crack band model
allows for description of the constitutive behavior of the FPZ by a triaxial nonlinear constitutive
law, which can be written in general form as:
σij = Fij (kl) (1.5)
where σij , kl are the stress and strain tensors, i, j, k, l are Cartesian tensorial subscripts and Fij is a
function of a history functional. Such constitutive law takes into account all other stress components
that can exist with or without the presence of σyy, which allows for a much richer description of
mechanical behavior. For instance, the shear components can lead to change of the direction of
propagation of the crack band, whereas the normal compressive stress components parallel to the
crack plane can precipitate failure and are able of causing splitting fracture, even in the absence of
the σyy stress component.
In finite element analysis, the crack band is modeled as a stack of finite elements undergoing
damage. The constitutive behavior of the elements inside the localization band can exhibit gradual
softening, but it has also been checked [11] that the model is applicable for brittle failure as well,
1.2. NUMERICAL MODELING OF QUASIBRITTLE FAILURE 9
since the stress intensity factor KI calculated for the case of elements exhibiting a sudden post-peak
stress drop was very close to the LEFM prediction (within 1% error at most). On the other hand,
the cohesive zone model is represented as an interface element of zero thickness, and, in contrast to
the crack band model, requires topological changes in the mesh in order to capture the displacement
discontinuity arising from crack initiation at the element boundaries. It can be considered as an
approximation of the crack band model neglecting the effect of the finite FPZ width.
Use of the crack band model in finite element analysis, which is more convenient than use of
the cohesive crack model, requires the user to specify the constitutive behavior of an element of size
equal to the crack band width. Apart from being a numerical parameter, the size of the elements
forming the localization band, (crack band width), is physically related to the width h of the FPZ,
which is a material property. The crack band width can be measured experimentally as the minimum
possible spacing of parallel cracks in specimens that have been designed so as to fail with distributed
cracking [27]. However, if the user selects elements of size he = h∗ larger than the material property
of h (for instance for reasons of computational efficiency), a scaling of the post-localization part
of the input constitutive law according to the chosen element size h∗ is required, so as to always
dissipate the same amount of energy during failure of the elements inside the developed band of
localized strains, independently of their size. This mesh-dependent adjustment of the post-peak
part of the constitutive law is essential in order to ensure the mesh objectivity of finite element
simulations of strain-softening continua for the case of failure with strain localization.
The crack band model cannot address the issue of mesh-induced directional bias and thus its
efficiency in mitigating mesh dependence is maximized when the mesh lines are aligned with the
direction of crack band growth (Fig. 1.5a)). In case when the mesh lines are inclined with respect to
the direction of crack propagation, the crack band width h can be determined by applying correction
factors to the element dimensions (Fig. 1.5b)). In general, the method of estimate of the effective
band width in finite element simulations using the crack band model can significantly affect its
performance in mitigating mesh dependence [79], since the scaling of the constitutive law according
to this model is based on the width of the FPZ, which may not always coincide with the finite
element width (Fig. 1.5b)).
Finally, it should be noted that the crack band model regularizes the global load-displacement
behavior of the structure. However, the numerically resolved width of the FPZ decreases and tends
to zero upon mesh refinement [25]. For more complex constitutive laws, it becomes increasingly
hard to identify which parameters govern the global behavior of the structure and should thus be
adjusted so as to mitigate mesh dependence by regularizing the energy dissipation in the softening
1.2. NUMERICAL MODELING OF QUASIBRITTLE FAILURE 10
5
a) b)
a) Crack band parallel to the element sides in a square mesh
b) Crack band along the element diagonals in a square mesh
Figure 1.5: a) Crack band parallel to the element sides in a square mesh, and b) Crack band along
the element diagonals in a square mesh
region. Other limitations of the crack band model include its rather simplified assumption of the
failure mechanism of the structure. It assumes that damage localization occurs as soon as damage
initiates. Meanwhile, as discussed earlier, the crack band model is also inefficient in mitigating mesh
dependence of stochastic finite element simulations, since it does not take into the randomness in
the location of crack band formation inside an element, as originating from the underlying random
field of material strength. The present dissertation particularly focuses on improving the energy
regularization methods within the crack band concept so as to overcome these last two limitations.
1.2.3 Nonlocal Damage Models
The development of cavities in the microscopic, mesoscopic and macroscopic processes of fracture
in materials together with the resulting deterioration in their macroscopic properties is called dam-
age [83, 94, 108]. Continuum damage mechanics aims at the analysis of damage development in
mesoscopic and macroscopic fracture processes in the framework of continuum mechanics. More
specifically, the mechanical effects of progressive microcracking, void nucleation and growth are rep-
resented by a set of state variables which act on the elastic and/or plastic behavior of the material
at the macroscopic level [118].
Description of the softening behavior of quasibrittle materials through continuum damage me-
chanics models can lead to localization instability and the spurious mesh sensitivity of finite element
analysis. One of the simplest remedies to alleviate mesh dependence of the global response is the
crack band model, described in Section 1.2.2. However, because continuum damage models are for-
mulated in terms of a stress-total strain relation, rather than an additive decomposition of the strain
to its elastic and inelastic part, it is not easy to scale only the postlocalization part of the strain, as
1.2. NUMERICAL MODELING OF QUASIBRITTLE FAILURE 11
dictated by the crack band model, while keeping the unloading part unaffected [25]. Contrariwise,
smeared crack models [17] decompose the strain into elastic and inelastic parts and hence make
application of the crack band model more convenient.
Another, more rigorous method of mitigating mesh dependence, are nonlocal formulations, which
are able to fully regularize the problem in the sense that the resulting strain field is highly concen-
trated in narrow zones but remains continuous [78]. Not only are nonlocal models able to suppress
mesh sensitivity of numerical simulations, but they also prevent strain localization into arbitrarily
small regions, which is why they are often characterized as true localization limiters, differentiating
them from partial localization limiters such as the crack band model and rate-dependent models.
The general approach in nonlocal models is to replace a certain variable at a given point by its
nonlocal counterpart, which is obtained by weighted averaging in the spatial neighborhood of the
point under consideration. In the first formulations of nonlocal models [24], the stress tensor at a
point x is calculated from the nonlocal strain tensor at the same point ¯(x), which is defined by the
following relation:
¯(x) =
∫
V0
α(x,x′)(x′)dV (x′) (1.6)
where α(x,x′) is an empirical bell-shaped weighting function (Fig. 1.6a)), and V0 is the volume of
the neighborhood whose size is controlled by a characteristic length of the material.
Since one of the goals of the development of nonlocal models has been preventing strain localiza-
tion, which most often occurs after inelastic processes initiate, nonlocality may be introduced only
in the variables that control the inelastic behavior, e.g. damage variable, and in this case nonlocal
models can reduce to standard local models at low strain values. This prevents some problematic
issues related to the original formulations, that applied the nonlocal averaging (Eq. (1.6)) to the
elastic part of the strain as well, leading to dispersive wave propagation even in homogeneous elastic
continua [15].
The efficiency of nonlocal models in preventing strain localization into regions of diminishing
size can be explained as follows. If strains tend to localize in a very narrow region, the nonlocal
strain, by its definition (Eq. (1.6)), will become high not only in the material point with high local
strain, but also in its immediate neighborhood. The damage variable, which depends on the nonlocal
strain, will grow in the neighborhood of the material point as well. Considering the simple case of a
uniform stress field (e.g. static equilibrium of a bar under uniaxial tension with zero body forces),
damage growth in the close neighborhood will lead to increase of the local strains so as to keep the
stress constant, therefore localization of strain and damage in an infinitesimal region is prevented.
The region of localized strains has a finite size, which is dictated by the length scale involved in the
1.2. NUMERICAL MODELING OF QUASIBRITTLE FAILURE 12
6
a) Weighting function
b) Boundaries
Rescaling Dirac Delta at 
the boundary
Dirac Delta at the 
point of interest
a) b)
Figure 1.6: a) Nonlocal weighting function, and b) treatment of the weighting function near the
structural boundaries
nonlocal weighting function.
Numerical implementation of nonlocal models has to deal with difficulties arising near the
structure boundaries, where the averaging neighborhood tends to protrude outside the boundaries
(Fig. 1.6b), upper). Attempts to resolve this issue include rescaling of the weighting function within
the body [25] or placing a Dirac delta function either along the structural boundary, or at the center
point of the nonlocal integral [38] (Fig. 1.6b), lower). However, these approaches were not based on
physical arguments. The problem with the boundary can be eliminated by the nonlocal boundary
layer model [16]. In this model, a boundary layer of specified thickness related to the material char-
acteristic length is excluded from the body in which the nonlocal averaging is carried out. Other
approaches are stress-based nonlocal damage models [67], which consider a non-fixed size of the
averaging region, and have been extended to rate-dependent behavior as well [116]. The efficiency of
nonlocal models in preventing strain localization can be undesirable in dynamic loadings, where it
was shown that use of nonlocal models can lead to an overspread of the damaged region and inability
of description of localized failure [55]. As a remedy to this problem, models that transition from
nonlocal to local with increasing strain rate can be considered, as suggested in [55] and successfully
developed in [116].
1.2. NUMERICAL MODELING OF QUASIBRITTLE FAILURE 13
1.2.4 Phase-Field Models
Phase-field models are numerical tools for modeling fracture within the framework of continuum
mechanics, and they are closely related to the variational approach to brittle fracture [63], which was
introduced in order to overcome limitations of the theory of fracture mechanics. Fracture mechanics,
both the linear elastic [70] as well as the nonlinear cohesive zone models [6], is not a self-contained
theory in the sense that additional criteria are required to predict phenomena such as initiation of
crack propagation, to estimate the crack path and to determine the onset of branching instabilities
in dynamics [30]. Phase-field models, in contrast, being solely based on energy minimization, are
able to automatically determine crack nucleation, growth and coalescence, without requiring an
assumption of predefined cracks and without suffering from the pathological mesh dependence of
other continuum approaches to fracture.
As an introduction to the phase-field models, we consider an elastic solid Ω ⊂ Rndim with a set
of sharp cracks Γ ⊂ Rndim−1, where ndim = 1, 2, 3 is the number of spatial dimensions (Fig. 1.7a)).
The external boundary of the solid ∂Ω ⊂ Rndim−1, with outward unit normal vector n, consists of
two disjoint parts ∂Ωu and ∂Ωt, i.e., ∂Ωu ∩ ∂Ωt = ∅ and ∂Ωu ∪ ∂Ωt = ∂Ω, subjected to prescribed
displacements u∗(x) for x ⊂ ∂Ωu and tractions t∗(x) for x ⊂ ∂Ωt, where x are the spatial coordi-
nates of the material particles of the solid. Kinematics of motion of the solid are described by the
displacement field u(x) and the infinitesimal strain field (x) := ∇su(x) for the symmetric gradient
operator ∇s(·) with respect to x. The potential energy functional of the external loads acting on
the solid is given by:
P(u) =
∫
Ω
b∗ · u dV +
∫
∂Ωt
t∗ · u dA (1.7)
where b∗ are the prescribed volumetric body forces.
In phase-field models, the sharp topology of cracks is regularized by a diffuse damage band of
finite thickness, characterized by a scalar phase-field parameter φ, where φ = 1 represents the cracked
material and φ = 0 the intact one [39, 40], as shown in Fig. 1.7b). Considering quasi-static brittle
fracture of isotropic elastic solids and the regularized crack topology, the surface energy Ψc(Γ) is
approximated by:
Ψc(Γ) =
∫
Γ
Gc dA ≈
∫
Γ
Gcγ (φ;∇φ) dV, γ (φ;∇φ) = 1
2
[
1
l0
φ2 + l0 (∇φ · ∇φ)
]
(1.8)
where Gc [J/m2] is the critical energy release rate of the material, γ (φ;∇φ) is the crack surface
density function [107], and l0 is a length scale that determines the width of the diffuse crack. The
1.2. NUMERICAL MODELING OF QUASIBRITTLE FAILURE 14
7
a) b)
A solid body Omega with the crack set Gamma:
a) Sharp cracks
b) Approximated diffuse crack bands
Figure 1.7: A solid body Ω with the crack set Γ: a) sharp cracks, and b) approximated diffuse crack
bands
length scale l0, though initially treated as a purely numerical parameter, was later considered to be a
material property that can be measured experimentally [110]. When l0 → 0, the solution converges
to that of the original problem with the sharp crack topology according to the Γ convergence theorem
[42]. It has also been shown that the results of the phase-field model model can be made independent
of t e length scale l0 [136] provided there is sufficient resolution of the crack phase-field, which was
not the case in earlier formulations of the model. Because of the phase-field regularization of the
sharp crack Γ, the stored strain energy of the elastic solid Ψs is a function of both the displacement
field u as well as the crack phase field φ and is expressed as follows:
Ψs (u, φ) =
∫
Ω
g(φ)ψ0 ((u)) dV (1.9)
where g(φ) is an energetic degradation function, which leads to a decrease of the initial elastic
energy density ψ0() with increasing values of the crack phase-field parameter φ. The total energy
functional of the solid is thus given by:
E(u, φ) =
∫
Ω
g(φ)ψ0 ((u)) dV +
∫
Ω
Gc
1
2
[
1
l0
φ2 + l0 (∇φ · ∇φ)
]
dV − P(u) (1.10)
The displacement field u and the crack phase-field φ at every instant are found by minimizing the
total potential energy of the cracking solid (Eq. (1.10)) subject to appropriate boundary conditions.
Similar to Griffith’s energetic approach to fracture, crack propagation results from the competition
1.2. NUMERICAL MODELING OF QUASIBRITTLE FAILURE 15
between the elastic strain energy stored in the bulk material and the surface energy. On the one
hand, deformations of the elastic solid due to loading increase the elastic energy. Any excess of the
stored elastic energy beyond the energy required for fracture, which is the surface energy in the case
of Griffith’s perfectly brittle fracture considered here, will lead to crack propagation by a continuous
reduction of the strain energy and a corresponding increase of the surface energy (all other dissipative
mechanisms ignored) through increase of the crack phase-field φ towards unity. Thus, phase-field
models can be regarded as solving fracture mechanics problems using partial differential equations
(PDEs): basically one solves two PDEs; one for the vectorial displacement field u and another one
for the scalar phase-field φ, through alternating minimization of each field in Eq. (1.10).
Phase field models have also been developed for ductile [3,56] and cohesive fracture [41,104,111,
131,132]. In the case of ductile fracture, the stored energy functional is decomposed into elastic and
plastic parts. The mechanisms of plastic deformation and crack growth may be considered either
as occurring independently [3,56], or the interaction between them can modeled by considering the
plastic energy term to affect the evolution of the crack phase-field [106]. A review of various phase-
field models for ductile fracture can be found in [2]. For the case of cohesive fracture, modifications
of the surface energy term are introduced so as to approximate the gradual release of fracture
energy with increasing crack opening displacement, which is the essential feature of this type of
failure. Extension of phase-field models to dynamics can be straightforward if the kinetic and
fracture energy are assumed to be independent of the crack phase-field. Successful attempts of
predicting crack branching [33] and providing quantitative estimates of the crack speed [34] have
been made. However, accurate prediction of the paths of ultra-high-speed cracks requires relaxing
the aforementioned assumptions and considering velocity or rate-dependent fracture energy [1] and
phase-field degraded kinetic energy [44], which is a more challenging task.
Advantages of phase-field models for fracture include prediction of phenomena such as crack initi-
ation, growth and coalescence purely from energy minimization and without any further (sometimes
ad hoc) criteria required, ability to deal with the merging and branching of multiple cracks with no
additional heuristics, as well as computational implementation that is straightforward in any num-
ber of dimensions. However, phase-field models are rather computationally expensive, since accurate
resolution of the gradient term in Eq. (1.8) requires sufficiently fine mesh inside the damaged zone.
Besides, more efforts are still required in order to make phase-field models capable of providing
accurate predictions of crack tip locations. Finally, lack of representation of the displacement dis-
continuity caused by fracture is a drawback of phase-field models and of all continuum approaches
to fracture. Nevertheless, the high performance of phase-field models in predicting complex crack
1.3. PROBABILISTIC MODELING OF STRENGTH DISTRIBUTION OF QUASIBRITTLE
STRUCTURES 16
patterns, their connection to both fracture mechanics and continuum damage mechanics as well as
their rigorous variational structure render them worthy of further research efforts, as also concluded
in a review of phase-field models by [137].
1.3 Probabilistic Modeling of Strength Distribution of Qua-
sibrittle Structures
The foregoing discussion focuses on deterministic analysis of quasibrittle fracture. Due to material
heterogeneity, the response of quasibrittle structures is inevitably stochastic. Up to present, there
have been developed three analytically tractable probabilistic models of structural failure: a) the
weakest-link model, b) the fiber-bundle model, and c) the fishnet model (Fig. 1.8). The weakest-
link and fiber-bundle models represent failure with damage localized into a single crack and failure
with distributed cracking, respectively, whereas the recently developed fishnet model is transitional
between the weakest-link and fiber-bundle models.
8
...
a)
b)
c)
a) Chain model
b) Fiber-bundle model , c) Fishnet model
Figure 1.8: Probabilistic models of structural failure: a) chain model, b) fiber-bundle model, and c)
fishnet model
1.3. PROBABILISTIC MODELING OF STRENGTH DISTRIBUTION OF QUASIBRITTLE
STRUCTURES 17
1.3.1 Chain Model
The chain model or series coupling model (Fig. 1.8a)) statistically represents the structure as a chain
of representative volume elements (RVEs), and is based on the assumption that structural failure will
occur as soon as one RVE fails. It should be noted that, for strength statistics, the RVE is defined
as the smallest material element whose failure causes the entire structure to fail, which is different
from the definition of the RVE in the homogenization theory of continuum mechanics. Each RVE
in the weakest-link model represents one of the potential fracture process zones (FPZs) inside which
damage can localize and lead to global failure. In Weibull’s theory of perfectly brittle structures [134],
the size of the FPZ is negligible compared to the structural dimensions, therefore the number of
RVEs in the weakest-link model can be considered as infinite. This hypothesis of brittleness is
not obeyed in quasibrittle failure, in which the size of the fracture process zone (FPZ) becomes
nonnegligible with respect to the system size and thus the number of RVEs is finite, as was considered
in the finite weakest-link model [18, 90]. The assumption of the failure mechanism employed in the
aforementioned weakest-link representations of the structure is applicable to structures of positive
geometry, for which the stress intensity factor at constant load increases with the crack length, and
thus load control conditions would lead to unstable crack propagation.
Assuming that the RVE strengths are statistically independent random variables with cumulative
distribution function (cdf) P1, the probability distribution of the structural strength Pf can be
written as follows:
Pf (σN ) = 1− [1− P1(σN )]n (1.11)
where n is the number of RVEs in the structure, and σN is a maximum load parameter with
the dimension of stress. Based on Eq. (1.11), the survival probability of the structure is the joint
probability of survival of all its RVEs. Recent studies on the statistics of static strength of quasibrittle
structures [14, 15, 18, 90] have shown that the RVE can be statistically modeled by a hierarchical
model, which consists of a combination of bundles and chains (Fig. 1.9a)). Each RVE represents
a potential FPZ, comprising a number of subscale cracks, as shown in Fig. 1.9b). The subscale
cracks which are aligned with the macrocrack direction can be considered to act approximately in
parallel coupling, i.e. they are constrained to open approximately equally and simultaneously while
the forces applied on them are summed. Transversely to the FPZ, the adjacent subscale cracks
act approximately in series coupling, i.e. they are subjected to approximately the same force and
their transverse deformations, or openings, are summed [23]. The above considerations justify the
hierarchical model of the RVE. Based on the hierarchical model, it was shown that the probability
1.3. PROBABILISTIC MODELING OF STRENGTH DISTRIBUTION OF QUASIBRITTLE
STRUCTURES 18
distribution of the static strength of the RVE can be approximated by a Gaussian distribution onto
which a Weibull tail (or, equivalently, a power-law tail) is grafted at a fairly low probability.
9
a) b)
c)
a) Hierarchical model of the RVE
b) Idealization of the force transmission mechanisms inside the 
FPZ
c) Stress-strain diagrams of the coupled elements
Figure 1.9: a) Hierarchical model of the RVE, b) idealization of the force transmission mechanisms
inside the FPZ, and c) stress-strain diagrams of the coupled elements
The weakest-link model has the following mathematical properties which have been proven ana-
lytically [15,18,90]: If the cdfs of strength of all the elements (links) have a power-law tail of exponent
m, then the cdf of strength of the entire chain has also a power-law tail of the same exponent m; and
each ten-fold increase of the number of elements in the chain extends the reach of its power-law tail
in terms of failure probability by an order of magnitude, thus for sufficiently large number of links,
the cdf of structural strength approaches the Weibull distribution. The latter property reflects the
fact that quasibrittle structures, whose strength statistics have been realistically described by the
finite weakest-link model, become perfectly brittle at the large size limit.
1.3.2 Fiber-Bundle Model
In contrast to the weakest-link model, which represents the damage localization mechanism, the fiber-
bundle model represents structural failure with distributed cracking. The structure is statistically
modeled by a set of elements, also called fibers, which are loaded in parallel, i.e. all fibers are
subjected to the same deformation whereas the load carried by the bundle is the sum of the loads of
the individual fibers (Fig. 1.8b)). Once a fiber fails, its load is redistributed among other surviving
1.3. PROBABILISTIC MODELING OF STRENGTH DISTRIBUTION OF QUASIBRITTLE
STRUCTURES 19
fibers, this way preventing global failure. The maximum load of the bundle is reached when a certain
fraction of the fibers breaks. Various phenomenological rules have been proposed in the literature
[53, 89, 117] in order to describe the stress redistribution after the failure of a fiber. According to
the equal or global load-sharing rule, the load of the failed element is equally carried by all the
remaining unbroken fibers, whereas in the local load-sharing rule, only the surviving fibers in the
vicinity of the failed element carry the redistributed load. Another approach is to derive the load-
sharing rule by assuming a constitutive behavior (e.g. brittle, plastic, softening) of the elements in
the bundle [15,90], as shown in Fig. 1.9c).
Assuming that the strengths of the individual fibers are statistically independent random vari-
ables with probability density functions fi(x), where i denotes the ith fiber, the cdf of bundle
strength, which is defined as the maximum of the average bundle stress, is given by the following
expression
Gn(S) = n!
∫
Ωn(S)
n∏
i=1
fi(σi)dσ1dσ2 . . . dσn (1.12)
where n is the number of fibers in the bundle model, and Ωn(S) is the admissible region of stresses
in all the fibers.
The fiber-bundle model has the following mathematical properties, which hold regardless of the
constitutive behavior of its elements [126]: If the cdf of strength of each fiber has a power-law tail
of exponent mi, then the cdf of bundle strength will also have a power-law tail, with an exponent
that is equal to the sum of the exponents of the power-law tails of the cdfs of fiber strengths, i.e.∑n
i=1mi. The extent of the power-law tail of the cdf of bundle strength in terms of failure probability
decreases with the number n of elements coupled in parallel, as described by the following relations:
Ptn ∼ (Pt1/n)n − (Pt1/3n)n for brittle bundles and Ptn ∼ (Pt1/n)n for plastic bundles [26], where
Pt1 is the failure probability at the point where the power-law tail of the strength cdf of a single fiber
terminates. Another important asymptotic property of the fiber-bundle model is the convergence of
the cdf of strength of sufficiently large bundles to the Gaussian distribution. This has been shown
in [53] for the case of brittle bundles, whereas for plastic bundles it is a consequence of the central
limit theorem, since the strength of a plastic bundle is the sum of strengths of all the fibers. The
latter property leads to the conclusion that the cdf of strength of ductile structures should follow
the Gaussian distribution, since in ductile failure the structural strength is a weighted sum of the
strengths of all the material elements along the failure surface.
1.3. PROBABILISTIC MODELING OF STRENGTH DISTRIBUTION OF QUASIBRITTLE
STRUCTURES 20
1.3.3 Fishnet Model
Both the hierarchical model of the RVE strength used in the finite weakest-link model as well as the
fiber-bundle model are based on simplified assumptions of the load transmission and redistribution
mechanisms inside the FPZ, which are necessary for the mathematical formulation of the probability
distribution of structural strength. In the recently developed fishnet probabilistic model [99,101,102],
which aimed to describe the strength distribution of nacreous materials, but was also shown able
of capturing the lower tail of failure probability of quasibrittle materials such as concrete or fiber
composites, the force transmission mechanism of the structure is modeled as a fishnet with diagonally
pulled links (Fig. 1.8c)). The links were initially considered as perfectly brittle, though in subsequent
formulations the fishnet model was generalized to links exhibiting a gradual softening postpeak
behavior (Fig. 1.9c)).
In contrast to the weakest-link model, which assumes that global failure will occur as soon as
a single link, the weakest one, fails, the fishnet model considers that several links can fail prior
to reaching the maximum load. Thus, the survival probability of the structure 1 − Pf (σ) at a
prescribed load σ can be expressed as the probability of the union of disjoint events, each one
representing survival under different number of failed links, as follows:
1− Pf (σ) = Ps0(σ) + Ps1(σ) + Ps2(σ) + . . . (1.13)
where Psk(σ) represents the survival probability of the entire fishnet under load σ for the case
when there have been exactly k failed links. For a two-dimensional fishnet consisting of r rows
and q columns of links (n = r × q), whose strengths are independent and identically distributed
(i.i.d.) random variables with cdf P1(σ) following the Gauss-Weibull grafted distribution, the first
term in Eq. (1.13) is the same as the finite weakest-link model and equals [1− P1(σ)]n, whereas
the remaining series terms differentiate the fishnet model from the weakest-link model. The terms
Psk(σ) with k > 0 can be calculated as the joint probability of failure of k number of links at
load σ, and survival of the n − k links under the corresponding level of applied load, which takes
into account the stress redistribution caused by the damages. For instance, for k = 1, Ps1 has the
following expression:
Ps1 ' nP1(σ) [1− P1(σ)]n−ν−1 [1− P1(ησ)]ν (1.14)
where it is assumed that ν number of links in the vicinity of the failed link will carry the redistributed
stress ησ, whereas the remaining n− ν − 1 will survive under the prescribed load σ. The analytical
1.3. PROBABILISTIC MODELING OF STRENGTH DISTRIBUTION OF QUASIBRITTLE
STRUCTURES 21
expressions of the terms corresponding to k ≥ 2 are rather complicated, but for brittle fishnets
and for not unusually large coefficient of variation of the link strengths, these higher terms are not
required since the peak load is reached when only a few links have failed. However, in the case of
links exhibiting a general softening behavior, the number of failed elements prior to the peak load
increases with decreasing brittleness of the links and thus truncating the series of Eq. (1.13) up to
the first few terms because of the difficulty of calculating the stress redistribution history for the
higher-order terms, cannot provide a reliable approximation of the survival probability of the entire
fishnet.
For this reason, a realistic and more computationally efficient approximation of the failure prob-
ability of fishnets with softening links was introduced in [100], by using the order statistics. More
specifically, it was considered that the maximum load of the structure is determined by the kth
weakest link, where k is a discrete random variable characterizing the extent of damage prior to
structural failure. The failure probability of the structure is then expressed as:
Pf (σ) = P(σmax ≤ σ) =
n∑
k=0
P(nc = k)Wk(γkσ) (1.15)
where nc is the number of failed links just before the peak load is reached,Wk(x) is the cdf of strength
of the kth smallest minimum (kth order statistics), and γk is the average stress concentration factor,
which depends on k. It has also been shown in [102] that the series expansion of Eq. (1.13) and the
approach using order statistics, Eq. (1.15), are equivalent at the lower tail of the structural failure
probability.
The fishnet shape, i.e., the width-to-length aspect ratio, where the length is measured in the
direction of the applied tensile loading, is particularly important for the probability distribution
of its strength. Increasing the aspect ratio from 0 to ∞, the fishnet gradually transitions from
the weakest-link model to the fiber-bundle model, which represent its limiting cases. However, in
contrast to the finite weakest-link model, in the fishnet model there is no fixed-size RVE, and the
predicted size effect is similar, but different from the Type I quasibrittle size effect described by the
finite weakest-link model. The main drawback of the fishnet model is that it requires the calculation
of stress redistributions after breakages of the links, a calculation which becomes intractable as the
number of failed links prior to the peak load increases.
1.4. RESEARCH OUTLINE 22
1.4 Research Outline
Despite the development of several successful numerical techniques to mitigate the spurious mesh
sensitivity of FE simulations of softening continua, known as the localization limiters, the mesh
dependence issue continues to attract significant research interest, as an essential part of the overar-
ching goal to continue improving the computational efficiency and reliability of numerical models of
failure and aid the design of safer structures. The highly sophisticated localization limiters, such as
nonlocal integral or gradient-type models, have been mostly used for academic applications, whereas
commercial codes use predominantly the simplest localization limiter, that is the crack band model,
owing to its convenient implementation and low computational cost.
An important limitation of the widely-used crack band model, as already discussed, is that its
ability to recover mesh-objective structural response is limited to cases of failure with severe strain
localization. However, failure of quasibrittle structures can manifest a complex damage process,
which may include a transition between localized and diffused damage patterns. Moreover, conditions
that may suppress strain localization (e.g. high strain-rate loading) render the crack band model
and generally localization limiters less appropriate. It is thus clear that a general approach for
addressing the issue of mesh dependence should be tied to the actual spatial and temporal evolution
of damage pattern in a structure, which is one of the motivations of the current study.
In Chapter 2, we focus on FE simulations of softening continua in the deterministic setting, and
propose an energy regularization method that adjusts the constitutive law according to the current
level of damage localization, as described by a numerical parameter referred to as the localization
parameter. The method is cast into an isotropic damage model, and is further extended to rate-
dependent behavior. The model is shown to be efficient in mitigating mesh dependence of the global
structural response, by its application to FE simulations of both static and dynamic failures of
ceramic specimens. Even though the proposed energy regularization scheme is applied to an isotropic
continuum damage mechanics model suitable for tension-governed failure, extensions of the model
to account for anisotropy, mixed-mode failure, and more complicated constitutive behaviors can be
pursued in the future.
The present research is further extended in the stochastic setting, since it is known that the
properties of quasibrittle materials could exhibit a significant degree of spatial variability. Effective
discretization of the random fields of constitutive properties of a given material while preserving
their original probability distribution and correlation structure is an issue which has been extensively
studied in the literature [65,71,129]. In contrast, in this study, we shift our attention to a mechanistic-
1.4. RESEARCH OUTLINE 23
based representation of the randomness of finite element properties, which takes into account both
the inherent variability as well as the failure pattern of the considered element. Since focus is
placed on the use of mesh sizes larger than the autocorrelation lengths of the underlying random
fields, the apparent fracture properties of different finite elements can be considered as statistically
independent random variables, and can thus be fully defined by their input probability distributions.
This approach constitutes a more efficient method, since it allows use of large element sizes, and it
also offers the benefit of retaining the same spatial discretization for both the mechanical problem
as well as the modeling of randomness of fracture properties, in contrast to stochastic FE methods.
However, it is important to take into account the size effects on the sampling probability distributions
of fracture properties that can arise due to the failure mechanism of an element, otherwise localization
of strains into the weakest elements will lead to the output probabilistic structural response be
dependent on the spatial discretization, as explained in Section 1.1.
Different probabilistic models of failure have been developed for the various failure mechanisms
of quasibrittle structures (Section 1.3). The models capture the size effects on the strength statis-
tics. Even though the described failure models consider statistically independent strengths of their
elements, they are able to partially account for the presence of correlations in the system through
their assumptions or approximations of the load redistribution mechanism after local failures. How-
ever, given that the structural failure mechanism cannot generally be known a priori, the numerical
parameters describing the evolving damage pattern introduced in our energy regularization scheme,
are also used to determine the sampling distribution functions of the constitutive properties. In this
way, a general probabilistic model of apparent fracture properties is developed in Chapter 4, which
can transition from a weakest-link model representation of structural failure in cases of damage lo-
calization, to a fiber-bundle model in cases of distributed damage. Insight for the formulation of this
model is gained by stochastic discrete element simulations of the size effects on the statistics of struc-
tural strength of ceramic specimens experiencing different failure mechanisms, triggered by loading
with a wide range of applied tensile strain rates (Chapter 3). The probabilistic model is applied to
continuum stochastic FE simulations of quasibrittle structures experiencing different failure mech-
anisms. The simulation results show that the model can effectively mitigate mesh dependence of
the predicted first and second statistical moments of structural strength. The probabilistic model is
intimately tied with the proposed energy regularization scheme and is not based on any assumptions
for the structural failure mechanism. Therefore, the present research contributes to the development
of a general and unified framework for addressing the mesh dependence of continuum FE simulations
of quasibrittle failure, in both deterministic and stochastic settings.
Chapter 2
Energy Regularization in
Deterministic Numerical Simulations
of Strain-Softening Continua
As discussed in Chapter 1, the strain-softening behavior of quasibrittle materials can lead to localiza-
tion instability and the spurious mesh sensitivity of finite element simulations. Energy regularization
methods in the concept of the crack band model [17] are based on the assumption of a fully local-
ized damage pattern, which occurs as soon as the peak load of the structure is reached. However,
depending on the loading configuration and structural geometry, quasibrittle structures can exhibit
a complex damage process, which may involve both localized and diffused damage patterns. In this
chapter, a generalized energy regularization method, that considers the spatial and temporal evolu-
tion of the internal damage variable, is presented. The method introduces a localization parameter,
which describes the local spatial distribution of damage at every time step, and can capture the
transition from diffused to localized damage during the failure process. It is the localization param-
eter that governs the energy regularization of the constitutive model, obviating the need for any a
priori knowledge of the failure mechanism of the structure. The method is cast into an isotropic
damage model, and is further extended to rate-dependent behavior. The efficiency of the model in
mitigating mesh dependence is shown by its application in finite element simulations of both static
and dynamic tension-governed failure of quasibrittle specimens. The simulation results demonstrate
the essential role of mechanism-based energy regularization of constitutive relation in continuum
24
2.1. MODEL FORMULATION 25
finite element modeling of quasibrittle fracture.
2.1 Model Formulation
The proposed energy regularization scheme is incorporated into an isotropic continuum damage
mechanics model [83,94,119], which is represented in general form by
σ = f(ω)C :  (2.1)
where C denotes the elastic stiffness tensor, ω is the damage variable ranging from 0 (virgin state)
to 1 (fully damaged state),  is the infinitesimal strain tensor, and f(ω) is a damage function, which
describes stiffness degradation of the material. The function f(ω) has a maximum value of one
at ω = 0, and a minimum value of zero at ω = 1. Meanwhile, f(ω) must also be monotonically
non-increasing, i.e. df/dω ≤ 0. The damage parameter represents a state variable, which can be
related to the strain tensor. The core of continuum damage models largely lies in the formulation
of this relationship. In this study, the damage evolution equation is formulated by linking damage
mechanics with fracture mechanics through an energy regularization scheme.
Considering a finite element with linear shape functions and a given Gauss point inside its
domain, let ~ng denote a unit vector perpendicular to the crack plane. The crack plane is assumed
to be perpendicular to the direction of the maximum principal tensile strain at the onset of damage,
and is considered to have a fixed orientation throughout the damage process - fixed crack model -
even though the principal strain directions can change during loading. The dimension of the element
projected in the direction of ~ng is denoted as h1, and h2 denotes the corresponding dimension in
the direction orthogonal to h1. It is noted that, inside a given element, each Gauss point may
have a different ~ng, and therefore the corresponding length h1, which physically represents the
crack band width, may also differ among Gauss points of the same element. Various numerical
estimates of the crack band width have been proposed in the literature, and the performance of the
crack band model [17] in addressing the mesh dependence of finite element simulations when using
different estimates of the localization bandwidth can vary significantly [79]. The modified projection
method [79], which considers a single value of h1 for all Gauss points of an element based on the
principal values of its average strain tensor, could also have been used in the present model. Besides,
the estimate of the crack band width h1 should also take into account the order of approximation
of the displacement field, since for higher order shape functions, damage could localize into a single
2.1. MODEL FORMULATION 26
Gauss point [79], but in this study we limit our focus to linear elements.
In order to derive the damage evolution law by employing fracture mechanics, it is not unduly
restrictive to consider that ωh2 represents the crack length a, which is similar to a common defi-
nition of the damage variable in continuum damage mechanics [83]. The energy required (per unit
thickness) for advancement of a crack by distance δa is equal to Gfh2δω, where Gf represents the
mode I fracture energy, since this study focuses on tensile damage. Crack advancement is made
possible through the release of energy from the entire element during an increment of damage δω,
which can be written as h1h2δU , where U represents an energy potential governing tensile failure.
Considering quasi-static crack propagation, and fracture as the only energy dissipation mechanism,
the energy balance equation in the limit of δω → 0 can be written as follows:
∂U
∂ω
+
Gf
h1
= 0 (2.2)
Evidently, Eq. (2.2) is based on the assumption that the damage of the element is manifested by a
single crack (Fig. 2.1a)). This corresponds to the damage localization mechanism. Eq. (2.2) indicates
that the constitutive relation varies with the mesh size h1.
On the other hand, if the element experiences diffused damage (Fig. 2.1b)), there would be h1/h0
number of cracks propagating together, where h0 = width of the fracture process zone (FPZ). Note
that h0 is a material length scale, which, under quasi-static loading, is typically on the order of two
to three times the size of material heterogeneity. In this case, the energy balance equation becomes
∂U
∂ω
+
Gf
h0
= 0 (2.3)
In contrast to Eq. (2.2), Eq. (2.3) suggests a mesh-independent constitutive behavior. To further
explain this difference, we emphasize that the constitutive model of a single finite element describes
the behavior of a material element of finite size. In the case of localized damage, only the material
inside the damage zone experiences cracking, which is manifested by strain softening, whereas the
material outside the damage zone undergoes unloading. Since the size of the damage zone is a
material constant, a larger material element implies a larger unloading zone. Consequently, the
overall mechanical behavior of the material element is dependent on the element size. In the case of
diffused damage, the entire material element experiences the same softening damage, and therefore
the mechanical behavior of the element is independent of its size.
Based on the foregoing discussion, it is evident that Eqs. (2.2) and (2.3) correspond to two distinct
damage mechanisms. In some cases, such as a pre-notched specimen, Eq. (2.2) can directly be applied
2.1. MODEL FORMULATION 27
a) b)
Figure 2.1: Two distinct crack patterns in a material element represented by a finite element: a)
localized cracking, and b) diffused cracking
since damage localization occurs right away. However, in other cases, the damage process involves
the transition from diffused to localized damage, and the onset of this transition is not known
a priori. Therefore, a general energy regularization method is required, which can automatically
adjust the constitutive model as the damage pattern evolves. To formulate such a method, the
key step is to determine the localization level from the evolving damage pattern. In recent studies,
an estimate for the onset of localization was calculated based on the spatial fluctuations of the
damage variable [31, 92]. However, in these studies, the formulation is either in a one-dimensional
setting [31] or does not consider the direction of damage propagation [92]. Here we develop an
improved formulation based on the concept proposed in [92].
Considering a finite element labeled by i, let ~vij denote the vector that connects the centroid
of element i with the centroids of its neighboring elements j (Fig. 2.2). For any given Gauss
point of element i, we determine the connection vectors that maximize and minimize the product
(~ng · ~vij)/|~vij |, and denote them by ~vm1 and ~vm2, respectively. It is clear that these two vectors
represent the directions of the two adjacent elements that most closely align with the direction of
the maximum principal strain of the Gauss point. We then define a localization parameter δω for
2.1. MODEL FORMULATION 28
Fig. 2
Figure 2.2: Schematic of finite element mesh for defining the localization parameter
the Gauss point as
δω =

〈
1− ω¯i1 + ω¯i2
2ω¯i
〉2
when ω¯i > 0
0 when ω¯i = 0
(2.4)
where 〈x〉 = max(x, 0), ω¯i = average damage value of element i, and ω¯i1, ω¯i2 = average damage
values of the two adjacent elements that are along the directions of ~vm1 and ~vm2. Eq. (2.4) indicates
that the localization parameter is governed by the relative difference between the damage of element
i and the average damage of its adjacent elements along the directions of ~vm1 and ~vm2. It is evident
from Eq. (2.4) that δω has a minimum value of zero, which corresponds to the case of a uniform
distribution of damage or the case where the damage level of element i is lower than the average
damage level of its adjacent elements. δω has a maximum value of one, which corresponds to the
case of a fully localized damage pattern.
By introducing the localization parameter, we now propose a general energy regularization equa-
2.1. MODEL FORMULATION 29
tion for each Gauss point:
∂U
∂ω
+
Gf
h0
g(δω) = 0 (2.5)
where we assume function g(δω) to have the following form:
g(δω) = 1 +
(
h0
h1
− 1
)
δκω (2.6)
where κ = constant. It is clear that Eq. (2.5) converges to Eq. (2.2) and Eq. (2.3) for the case of
localized damage (δω = 1) and diffused damage (δω = 0) respectively, regardless of the selected value
for κ. Exponent κ governs the transition between these two fundamental failure mechanisms as a
function of δω.
In this study, we choose a form of the free energy function suitable for tension dominated failure.
A simple choice is to consider the concept of equivalent strain, which has been widely used in
tension-governed continuum damage models [92,105]:
U =
1
2
f(ω)E¯2 (2.7)
where ¯ =
√√√√ 3∑
k=1
〈k〉2 (2.8)
E = Young’s modulus of the material, ¯ = equivalent strain, and k = the principal values of the
strain tensor. Note that the damage function in Eq. (2.7) is a natural consequence of the present
isotropic damage model (Eq. (2.1)). Substitution of Eq. (2.7) into Eq. (2.5) yields an evolution
equation of the damage function that is compatible with fracture mechanics:
df(ω)
dω
= − 2Gf
h0E¯2
g(δω)H(¯− 0)H( ˙¯) (0 ≤ ω < 1) (2.9)
where H(x) = Heaviside function, 0 = threshold strain, and ˙¯ = d¯/dt. Eq. (2.9) gives the energy-
based evolution law of the damage function, in which the energy regularization is governed by the
localization parameter. Note that the two Heaviside functions in Eq. (2.9) indicate that 1) damage
would not grow if ¯ ≤ 0, and 2) damage would not grow if the material experiences tensile unloading.
Note that, based on Eq. (2.1), function f(ω) = 0 when ω = 1, and therefore we have df(ω)/dω = 0
when ω = 1.
The original crack band model [17] applicable for localized damage was anchored by the preser-
vation of fracture energy (i.e. the product of the fracture energy density times the element size is
2.1. MODEL FORMULATION 30
set equal to the material fracture energy) [8,15,19]. In some cases, though, it can be hard to assess
which constitutive parameters should be adjusted based on the mesh size so as to recover mesh
objective simulation results. In the present model, there is no need for such an insight, since the
energy regularization is carried out by adjusting the point-wise tangential stiffness tensor through
Eq. (2.9).
The constitutive model is completed by combining Eq. (2.9) with an explicit form of the damage
function f(ω), which in this study is assumed to be (e.g. [97, 135])
f(ω) =
1− ω
1− aω (2.10)
where a = model parameter that is directly related to the material tensile strength ft, not to be
confused with the crack length a for which the same notation was used. To reveal this relation-
ship, we consider the case of monotonic uniaxial tension with applied strain . By differentiating
Eq. (2.10) with respect to the damage variable and combining it with Eq. (2.9), we obtain the
following relationship between the damage variable and the strain for  ≥ 0:
ω =
1
a
− 
a
√
h0E(1− a)
2Gfg(δω)
(2.11)
The condition of threshold strain (ω = 0 at  = +0 ) yields:
a = 1− 2Gfg(δω)
h0E20
(2.12)
Threshold strain 0 can be determined by investigating the tangential stiffness at  = +0 :
dσ
d
∣∣∣∣
+0
=
(
1 +
df(ω)
d
∣∣∣∣
+0
0
)
E = E/a (2.13)
Here we consider that the stress-strain response exhibits a softening response right after +0 , which
implies a < 0. In this case, 0 must correspond to the material tensile strength ft, i.e. 0 = ft/E,
and Eq. (2.12) becomes:
a = 1− 2lchg(δω)
h0
(2.14)
where lch = EGf/f2t = Irwin characteristic length. Eq. (2.14) indicates that the damage function
(Eq. (2.10)) is directly related to the localization parameter. For the case of localized damage, the
damage function explicitly depends on the mesh size, h1, whereas for the case of diffused damage it
2.2. EXTENSION TO RATE-DEPENDENT BEHAVIOR 31
does not. Fig. 2.3 shows the uniaxial tensile stress-strain curves for material elements of different
sizes for the case of localized damage. As seen, in this case the present model recovers the original
crack band model. Note that, due to the limit a < 0, Eq. (2.14) suggests a maximum element
size that can be used in the FE simulation: hm = 2lch. If the element size is larger than hm, the
post-peak tangential stiffness would become positive indicating a snap-back behavior.

[M
P
a]
Fig. 3Figure 2.3: Uniaxial tensile stress-strain curves for different mesh sizes for the case of fully localized
damage
2.2 Extension to Rate-Dependent Behavior
We now extend the model to dynamic tensile failure of quasibrittle materials. One way to capture the
rate effect on the constitutive behavior is through a kinetics model of damage growth [86]. Within
the present modeling framework, it is essential to incorporate the energy regularization scheme into
the formulation of the damage kinetics model, which can be applied to any arbitrary loading path.
The simplest way is to derive an equation for the damage growth rate directly from Eq. (2.9). To
this end, we introduce the rate dependence of the fracture energy and crack band width.
It is well documented that many quasibrittle materials exhibit a rate enhancement of fracture
2.2. EXTENSION TO RATE-DEPENDENT BEHAVIOR 32
energy [32,121,128]. In this study, we consider the following rate dependence on the fracture energy
Gdf (ω˙) = Gf exp [a1(ω˙/ω˙0)] (2.15)
where a1 = model constant, ω˙ = dω/dt = damage growth rate, and ω˙0 = reference damage
growth rate, which is defined as ω˙0 = CR/lch, where CR = the Rayleigh wave speed = (0.862 +
1.14ν)
√
E/2ρ(1 + ν)3 (ν = Poisson ratio, and ρ = material density).
The consideration of the rate effect on the crack band width is motivated by recent discrete
element simulations of dynamic fracture of ceramic specimens. The simulations showed that the
spacing of distributed cracks, which is a measure of the crack band width, decreases with an in-
creasing strain rate [91]. Therefore, we consider the following rate dependence on the crack band
width:
hd0(ω˙) = h0 exp [−a2(ω˙/ω˙0)] (2.16)
where a2 = model constant. By considering ω˙0 as the maximum damage growth rate, the minimum
crack band width predicted by Eq. (2.16) equals h0 exp(−a2). Under quasi-static loading, the crack
band width is about three times the maximum size of material heterogeneities [19]. It is reasonable
to consider that, under high loading rate, the minimum crack band width is about the maximum size
of material heterogeneities. Therefore, we have a2 ≈ 1.1. It should be mentioned that the present
model does not consider comminution and fragmentation of material heterogeneities. Nevertheless,
the effect of material comminution and fragmentation can be taken into account by introducing an
additional energy dissipation term to the constitutive model [9, 10].
The rate-dependent fracture energy and crack band width are incorporated into the evolution
equation of the damage function (Eq. (2.9)) as well as the function g (Eq. (2.6)). As a result, we
obtain a rate-dependent evolution equation of the damage function:
df(ω)
dω
= −2Gf exp [(a1 + a2)(ω˙/ω˙0)]
h0E¯2
gd(δω, ω˙) for ¯ ≥ 0 and ˙¯ > 0 (2.17)
where gd(δω, ω˙) = 1+
[
hd0(ω˙)/h1 − 1
]
δκω. However, the model parameter a is considered to be affected
only by the rate dependence of the localization parameter, i.e. a(δω, ω˙) = 1− 2lchgd(δω, ω˙)/h0. By
combining Eq. (2.17) with the assumed damage function f(ω) (Eq. (2.10)), we have
1
(1− aω)2 = exp [(a1 + a2)(ω˙/ω˙0)]
(0
¯
)2
for ¯ ≥ 0 and ˙¯ > 0 (2.18)
When using an explicit numerical scheme with a small time step, we may consider that the
2.2. EXTENSION TO RATE-DEPENDENT BEHAVIOR 33
localization parameter is evaluated based on the damage state calculated in the previous step, and
therefore can be kept constant during the current time step. Based on Eq. (2.18), the damage growth
rate for the current time increment can then be calculated as
ω˙ =
ω˙0
a1 + a2
ln
{
1
[1− a(δ0ω, ω˙0)ω]2
¯2
20
}
H( ˙¯)H(¯− 0) (0 ≤ ω < 1) (2.19)
where δ0ω, ω˙0 = localization parameter and damage rate calculated in the previous time step, re-
spectively. Note that Eq. (2.19) indicates a zero damage growth rate at ¯ = +0 , which ensures the
continuity of the damage growth function. Eq. (2.19) is further supplemented by the condition that
ω˙ = 0 when ω = 1. The combination of Eq. (2.19) and Eq. (2.10) completes the rate dependent
constitutive model.

[M
P
a]
Fig. 4
Figure 2.4: Rate-dependent uniaxial tensile stress-strain behavior of a material element for the case
of fully diffused damage
To demonstrate the constitutive behavior described by the present model, we consider a set of
model parameters for dense alumina [140]: E = 370GPa, ν = 0.24, ρ = 3900 kg/m3, ft = 200MPa,
Gf = 50 J/m
2, a1 = 2.5 and a2 = 1.1. These parameters are also used in the subsequent simulations
of quasi-static and dynamic failure. Fig. 2.4 presents the uniaxial tensile stress-strain curves of a
material element with diffused damage (i.e. gd(δω, ω˙) = 1) under different uniform applied strain
rates ˙. It is seen that, for all strain rates, the stress-strain curve deviates from the initial linear elastic
2.2. EXTENSION TO RATE-DEPENDENT BEHAVIOR 34
Figure 2.5: Rate effect on material properties: a) apparent tensile strength, b) fracture energy, and
c) Irwin’s characteristic length
response at the threshold strain 0. At increasing strain rates, the stress-strain response exhibits a
more pronounced nonlinearity before the peak strength is reached. Meanwhile, the present model
predicts a rate enhancement of material tensile strength, as shown in Fig. 2.5a). Note that, for the
case of constant strain-rate loading, the damage growth rate varies with time, and therefore the
model does not directly yield the dependence of the fracture energy on the strain rate. Nevertheless,
we may approximate this strain rate effect by using the average damage growth rate ω˙a = ˙/(u−0),
where u = ultimate strain at which the stress drops to zero. Fig. 2.5b) plots this approximated
strain rate effect on the fracture energy, which as expected exhibits a strong rate enhancement. The
rate effects on tensile strength and fracture energy directly provide the rate dependence of Irwin’s
characteristic length, as shown in Fig. 2.5c). It is seen that Irwin’s characteristic length, which is
proportional to the FPZ length, decreases with an increasing strain rate. Such a trend has been
experimentally observed in quasibrittle materials [13,19].
2.3. SIMULATIONS OF QUASI-STATIC FAILURE 35
2.3 Simulations of Quasi-Static Failure
Implemented in an open-source FEM software (OOFEM) [115], the present model is first applied to
simulate quasibrittle fracture under quasi-static loading. We use dense alumina as a model system.
The simulations consider both notched and unnotched beams subjected to three-point bending under
displacement controlled loading (Fig. 2.6). The unnotched beam has a size of 7.5 mm (length)×
2.5 mm (depth)×0.0375 mm (width). The notched beam has a size of 15 mm (length)× 5 mm
(depth)×0.0375 mm (width), and has a notch of length equal to half of the beam depth placed at
the mid-span. The quasi-static crack band width is chosen as h0 = 37.5 µm, which is about three
times the maximum grain size [91].
Fig. 6
P,
P,
Figure 2.6: Notched and unnotched three-point bend beams considered in the simulation (all dimen-
sions are in mm)
In the simulation, the beam is discretized with a uniform mesh of linear isoparametric brick
elements, each of which contains eight Gauss integration points. The dimensions of the finite elements
in the reference case are hx = h0, hy = 2h0, and hz = h0, where subscripts x, y, z denote the
directions of the beam span, depth and width, respectively. In the present simulation cases, the
mesh size along the direction of the maximum principal strain is approximately equal to the quasi-
static crack band width h0. To investigate the mesh sensitivity, we consider other two mesh sizes:
hx = 2h0 and hx = 4h0 while hy and hz are kept the same as those in the reference case. It is
noted that, for the notched beam, the notch width is set to be equal to the current dimension hx of
the finite element. Therefore, the notch width varies with the element size. Nevertheless, the ratio
between the notch depth and the notch width is large (on the order of 10) for all three mesh sizes,
so it is expected that the effect of the notch width on the simulation results is insignificant.
2.3. SIMULATIONS OF QUASI-STATIC FAILURE 36
For the purpose of comparison, we use three different approaches of energy regularization for
each simulation case:
I) Direct use of Eq. (2.2), which assumes damage localization in all cases.
II) Direct use of Eq. (2.3), which assumes diffused damage in all cases.
III) The present model, which handles diffused damage, localized damage, as well as the transition
between them (Eq. (2.5)).
It is evident that approach I resembles the original crack band model [17]. Approach II implies a
mesh independent constitutive model. Hereafter, we refer to approaches I and II as the localized
damage model and the diffused damage model, respectively. Approach III (i.e. the present model)
requires an additional parameter κ for determining the transition between different damage patterns.
As will be discussed later, the κ value will be determined based on the simulation of the unnotched
beam.
For the notched beam, damage localization occurs right away at the notch tip due to the
stress concentration at the pre-existing notch. This is evidenced by the simulated damage pat-
tern (Fig. 2.7). Thus, it is expected that the localized damage model would effectively mitigate
the spurious mesh dependence. Meanwhile, damage localization implies that δω = 1 in the present
model, so regardless of the value of κ, the present model yields the same result as the localized
damage model. Therefore, we do not seek to determine the optimum value of κ from the simulation
of the notched beam.
Fig. 7
Figure 2.7: Evolution of the damage pattern of the notched beam
2.3. SIMULATIONS OF QUASI-STATIC FAILURE 37
Fig. 2.8a) presents the simulated load-displacement curves of the notched beam using the afore-
mentioned three approaches. It is found that the localized damage model and the present model
essentially give the same results, which are independent of the mesh size. However, the diffused
damage model yields a strong mesh sensitivity. When a large mesh size is chosen, the simulation
predicts a significant increase in both the peak load and the total area under the load-displacement
curve. This can be attributed to the fact that, for localized damage, the fixed stress-strain model
implies a higher fracture energy for larger finite elements. Fig. 2.8b) shows the simulated profile of
normal stress along the ligament using the diffused damage model. It is seen that, when using a
larger mesh size, the FPZ becomes longer. This explains why the simulated post-peak regime for
larger mesh sizes exhibits a more ductile behavior.
Fig. 8
0.0 2.5 5.0 7.5
displacement [µm]
0
1
2
lo
ad
[N
]
localized damage model
hx = 4h0
hx = 2h0
hx = h0
0.0 2.5 5.0 7.5
displacement [µm]
diÆused damage model
0.0 2.5 5.0 7.5
displacement [µm]
present model
°400 °200 0 200
stress [MPa]
3
4
5
ve
rt
ic
al
co
or
d
in
at
e
y
[m
m
]
notch tip
upper surface
hx = 4h0
hx = 2h0
hx = h0
a)
b)
Figure 2.8: Simulation results of the notched beam: a) load-displacement curves, and b) profile of
normal stress along the ligament for the diffused damage model
The failure process in the unnotched beam is more complicated. Present simulations show that
all three approaches give a qualitatively similar damage evolution, which is shown in Fig. 2.9. It
is clear that the beam first experiences a diffused damage pattern and gradually transitions to
2.3. SIMULATIONS OF QUASI-STATIC FAILURE 38
localized damage featured by macrocrack propagation. Therefore, it is expected that neither the
localized damage model nor the diffused damage model will be able to effectively mitigate the issue
of mesh dependence. This is confirmed by the simulated load-deflection curves of the unnotched
beam as shown in Fig. 2.10. It is seen that, when using the localized damage model, the simulation
predicts a lower peak load as the mesh size increases. On the other hand, when using the diffused
damage model, the simulation shows an opposite trend, where a larger mesh size leads to a slightly
higher peak load. The mesh dependence by using the localized damage model is considerably more
pronounced than that by using the diffused damage model. This indicates that the energy dissipation
of distributed cracking contributes significantly to the peak load capacity. In the localized damage
model, a larger mesh size underestimates the energy dissipation when a diffused damage pattern
is present, and therefore it leads to a lower peak load capacity. The fact that the diffused damage
model does not completely suppress the mesh sensitivity indicates that damage localization has
occurred at the peak load. Similar to the foregoing analysis of the notched beam, a larger mesh
would overestimate the energy dissipation for the case of localized damage and yield a higher peak
load capacity.
Fig. 9
Figure 2.9: Snapshots of the damage pattern of the unnotched beam during the loading process
The opposite trends of mesh dependence of the localized damage model and the diffused damage
model indicate that the present model should be able to mitigate the mesh dependence if an appro-
priate value of the parameter κ is chosen. It is found that κ = 1 leads to the best performance of the
present model (Fig. 2.10), where the simulation result is almost insensitive to the mesh size. The
essence of the proposed model is that, as the loading proceeds, the energy regularization gradually
switches from Eq. (2.3) to Eq. (2.2) as damage localization emerges from distributed damage. The
2.3. SIMULATIONS OF QUASI-STATIC FAILURE 39
0 2 4 6
displacement [µm]
0
2
4
6
lo
ad
[N
]
localized damage model
hx = 4h0
hx = 2h0
hx = h0
0 2 4 6
displacement [µm]
diÆused damage model
0 2 4 6
displacement [µm]
present model
Fig. 10
Figure 2.10: Simulated load-displacement curves of the unnotched beam
same κ value will be used in the subsequent analysis of dynamic fracture as well.
It should be mentioned that the simulations show a lack of convergent solution for part of the
post-peak region (colored in grey in Fig. 2.10), which indicates a snap-back behavior. The post-peak
region with the snap-back behavior includes the kinetic energy due to unstable crack propagation,
and therefore cannot be directly used to investigate the effect of different energy regularization
schemes on the calculated energy dissipation due to damage and fracture. Nevertheless, it is clear
from Fig. 2.9 that the snap-back behavior occurs only during the propagation of a single macrocrack,
and the energy dissipation of the specimen is primarily due to the energy expended for the macrocrack
propagation. Therefore, it is expected that, if we are able to trace the snap-back behavior using a
different loading scheme (e.g. crack-mouth opening displacement controlled test), the present model
would be able to correctly capture the overall energy dissipation during this loading stage, since the
localization parameter would approach the value of 1, corresponding to a localized damage pattern.
The foregoing analysis of quasi-static fracture of beams indicates that energy regularization
plays an important role in mitigating mesh dependence in FE simulations. Meanwhile, energy
regularization is intimately tied with the underlying damage mechanism. Neither the crack band
model nor the mesh-independent local constitutive model can provide a general approach for energy
regularization due to their underlying assumptions of the damage mechanism. Figs. 2.8a) and 2.10
show some promising features of the present model for overcoming the limitations of these previous
models.
2.4. SIMULATIONS OF THE DYNAMIC SPALLING TEST 40
2.4 Simulations of the Dynamic Spalling Test
In this section, we apply the model to simulations of dynamic tensile failure of quasibrittle struc-
tures. We consider the modified split Hopkinson bar test, also known as the dynamic spalling test
(Fig. 2.11a)). In this test, a compressive stress pulse is applied to one end of the bar (the left end
was considered here), whereas the other end is free. As soon as the compressive stress wave reaches
the free end, it will reflect as a tensile wave and travel back into the bar, causing tensile damage
(Fig. 2.11b)). Depending on the amplitude and duration of the compressive stress pulse, the spalling
test can produce high strain rates, e.g. up to the order of 102/s for concrete, or even up to the order
of 104/s for higher strength materials, such as ceramics [62]. One important measurement of the
spalling test is the time history of the velocity of the free end. It has been suggested that the differ-
ence between the peak velocity and the rebound velocity, also referred to as the pullback velocity,
can be used to estimate the intrinsic dynamic tensile strength of the material [60, 76,98,128].
a)
incident wave
reflected wave
free end
b)
c) d)
hx
Fig. 11Figure 2.11: a) Schematic of the modified split Hopkinson bar test, b) stress wave traveling in thebar, c) stress pulse used in the present simulations, and d) discretization of the bar
In this study, we simulate the spalling test on a dense alumina bar of square cross-section. The
bar has a size of 12 mm (length)× 0.225 mm (depth) × 0.225 mm (thickness). The constitutive
properties are chosen to be the same as those used to produce Fig. 2.4 except that here we set
2.4. SIMULATIONS OF THE DYNAMIC SPALLING TEST 41
Poisson’s ratio to zero in order to preserve one-dimensional wave propagation. A sinusoidal pressure
pulse of half period tp = 1 µs is applied to the left end of the bar as an initial condition (Fig. 2.11c)).
Different stress amplitudes σa = 0.5, 0.75 and 1.25 GPa are considered in the simulations.
The bar is discretized with one row of linear isoparametric brick elements (Fig. 2.11d)). To in-
vestigate mesh dependence, three different mesh sizes along the longitudinal direction are considered
in the simulations, i.e. hx = h0, 2h0, and 4h0. Nonlinear dynamic analysis is performed by using
an explicit time integration scheme. A convergence study is first conducted in order to determine
the required time step. In order to achieve a convergent result for all mesh sizes, the time step is
found to be 5 × 10−4 µs, which is about 12% of the minimum time that the sound wave needs to
propagate through the smallest mesh with size hx = h0. For each loading case, we consider the
same three energy regularization methods as in the foregoing simulations of quasi-static fracture
(i.e. approaches I-III). As mentioned, for the present model, we set κ value equal to 1 based on the
simulation results of the unnotched beam.
Fig. 12
1.0 1.5 2.0 2.5
time [µs]
0
10
20
fr
ee
en
d
ve
lo
ci
ty
[m
/s
] present
model
æa = 0.5GPa
0
10
20
fr
ee
en
d
ve
lo
ci
ty
[m
/s
] diÆused
damage
model
æa = 0.5GPa
0
10
20
fr
ee
en
d
ve
lo
ci
ty
[m
/s
] localized
damage
model
æa = 0.5GPa
1.0 1.5 2.0 2.5
time [µs]
0
10
20
30
40 æa = 0.75GPa
0
10
20
30
40 æa = 0.75GPa
0
10
20
30
40 æa = 0.75GPa
1.0 1.5 2.0 2.5
time [µs]
0
20
40
60
æa = 1.25GPa
hx = 4h0
hx = 2h0
hx = h0
0
20
40
60
æa = 1.25GPa
0
20
40
60
æa = 1.25GPaa)
b)
c)
Figure 2.12: Time histories of free-end velocity predicted by the localized damage model (first row),
the diffused damage model (middle row), and the present model (bottom row)
2.4. SIMULATIONS OF THE DYNAMIC SPALLING TEST 42
Fig. 2.12 presents the time histories of the free-end velocity simulated by the three approaches
together with the elastic solution represented by the dashed line. It is seen that both the diffused
damage model and the present model yield essentially the same mesh-insensitive results, whereas
the simulation by the localized damage model shows a strong mesh dependence. Clearly all three
methods give the same results up to the point of damage initiation (i.e. point where the result
deviates from the elastic solution). Note that the peak free-end velocity corresponds to the instant
when the peak of the incident stress wave has reached the free end of the bar. After tensile damage
occurs, the diffused damage model and the present model predict that the free-end velocity decreases
continuously over the time domain considered in the simulation, due to continuous energy dissipation
in the damaged zone. By contrast, the localized damage model predicts a local minimum followed by
a second ascending branch of the free-end velocity diagram for the large mesh size, which indicates
the arrival of a second compressive wave at the right end. This compressive wave arises from the
reflection of the previously formed tensile wave at the spalling plane in the completely damaged
elements. On the other hand, for smaller mesh sizes the velocity history does not show any local
minimum, since the values of the velocity continuously decrease due to energy dissipation at the
spall plane. This result implies that, for the localized damage model, elements of larger size are
more brittle, and therefore are more prone to damage as compared to small elements. This behavior
is closely tied to its energy regularization. The central idea of the localized damage model is that,
regardless of the mesh size hx, the energy expended (per thickness) for complete damage of an
element is equal to Gfhy. The aforementioned mesh-dependent results indicate that the localized
damage model underestimates the energy dissipation for damaging a single element.
To further analyze the energy dissipation, we plot in Fig. 2.13 the time history of the total energy
dissipation due to damage: i.e. Wd(t) = A0
∫
L
[W(x, t) − 12σij(x, t)ij(x, t)]dx, where A0 = cross-
sectional area of the bar, and W(x, t) = ∫ t
0
σij(x, t
′)˙ij(x, t′)dt′. Similarly to the aforementioned
results of the time-history of the free-end velocity, the diffused damage model and the present model
predict the same history of energy dissipation, which is mesh insensitive. However, the localized
damage model yields a strongly mesh-dependent result. For the three stress pulses considered here,
as the mesh size increases, the energy dissipation grows at a faster rate but the ultimate amount of
dissipated energy becomes smaller. This result further supports the foregoing conclusion that the
energy dissipation of each finite element is underestimated by the localized damage model.
The underestimation of the dissipated energy by the localized damage model indicates that the
failure pattern of the bar is not featured by damage localization. In fact, the agreement between
the results of the present model and the diffused damage model implies that the bar must rather
2.4. SIMULATIONS OF THE DYNAMIC SPALLING TEST 43
Fig. 13
1.0 1.5 2.0 2.5
time [µs]
0.0
0.2
0.4
model
present
æa = 0.5GPa
0.0
0.2
0.4
d
is
si
p
at
ed
en
er
gy
[£
10
°4
J]
model
damage
diÆused
æa = 0.5GPa
0.0
0.2
0.4
model
damage
localized
æa = 0.5GPa
1.0 1.5 2.0 2.5
time [µs]
0.0
0.5
1.0
æa = 0.75GPa
0.0
0.5
1.0
æa = 0.75GPa
0.0
0.5
1.0
æa = 0.75GPa
1.0 1.5 2.0 2.5
time [µs]
0
1
2
æa = 1.25GPa
hx = 4h0
hx = 2h0
hx = h0
0
1
2
æa = 1.25GPa
0
1
2
æa = 1.25GPaa)
b)
c)
Figure 2.13: Time histories of energy dissipation simulated by the localized damage model (first
row), the diffused damage model (middle row), and the present model (bottom row)
experience a diffused damage pattern. Figs. 2.14 and 2.15 show snapshots of the damage pattern at
different times, as captured from the simulations using the localized damage model and the present
model, respectively. It is seen that both models predict a large damage zone, the size of which is
governed by the length of the reflected tensile wave. Therefore, for each stress pulse, the size of the
damage zone is approximately constant regardless of the mesh size. Such damage pattern indicates
that each finite element should experience diffused cracking, and therefore the energy dissipation
due to damage is proportional to the mesh size hx. However, the localized damage model enforces
constant energy dissipation independent of the mesh size, which leads to underestimation of the
dissipated energy as the mesh size increases. Moreover, one can note that the localized damage
model predicts, at the same instant of time, more severe damage for elements of larger size, which
indicates that the crack growth rate is itself mesh dependent, with damage propagating faster at
larger elements (Fig. 2.14). This is again, characteristic of the more brittle behavior of larger
elements imposed by the crack band model. By contrast, the present model ensures constant crack
propagation velocity independent of the mesh size, as is verified by the overlapping damage profiles
of all mesh sizes at all instants of time (Fig. 2.15). The present model is able to correctly capture
2.5. CONCLUSIONS 44
Fig. 14
1.0 1.5 2.0 2.5
time [µs]
0
20
v
[m
/s
]
A B
æa = 0.5GPa
1.0 1.5 2.0 2.5
time [µs]
0
25 A
B
æa = 0.75GPa
1.0 1.5 2.0 2.5
time [µs]
0
50
A B
æa = 1.25GPa
8 10 12
x [mm]
0.0
0.5
1.0
d
am
ag
e
!
A
B
hx = 4h0
hx = 2h0
hx = h0
8 10 12
x [mm]
0.0
0.5
1.0
A
B
8 10 12
x [mm]
0.0
0.5
1.0
A
B
Figure 2.14: Snapshots of the damage pattern of the bar at different times simulated by the localized
damage model
the energy dissipation of the element based on the information of the evolving damage state, and
therefore predicts a mesh-insensitive damage pattern and overall energy dissipation. It should be
mentioned that the present model may exhibit mesh sensitivity if the mesh size is so large that there
is lack of sufficient resolution of the damage profile.
2.5 Conclusions
A general energy regularization method is developed for finite element (FE) simulations of quasib-
rittle fracture for both localized and diffused damage mechanisms [68]. The model is cast in the
framework of continuum damage mechanics, with an explicit relation to fracture mechanics. By
employing a localization parameter, the model adjusts the energy regularization of the material con-
stitutive relation based on the prevalent local damage mechanism. The model is further extended
to rate-dependent constitutive behavior, in which the formulation of the damage growth rate is
seamlessly combined with the energy regularization scheme. The model can easily be implemented
in finite element codes, and is suitable for large-scale FE simulations.
The present model is applied to analyze both static and dynamic fracture of quasibrittle struc-
tures. It is shown that the model recovers the crack band model for the case where the structure
experiences a purely localized damage pattern. For structures which exhibit a transition from dif-
fused damage to localized damage, the present model outperforms both the localized damage model
and the diffused damage model in terms of suppressing the spurious mesh sensitivity. In the dynamic
2.5. CONCLUSIONS 45
Fig. 15
1.0 1.5 2.0 2.5
time [µs]
0
20
v
[m
/s
]
A B
æa = 0.5GPa
1.0 1.5 2.0 2.5
time [µs]
0
25 A B
æa = 0.75GPa
1.0 1.5 2.0 2.5
time [µs]
0
50 A
B
æa = 1.25GPa
8 9 10 11 12
x [mm]
0.0
0.5
1.0
d
am
ag
e
!
A
B
hx = 4h0
hx = 2h0
hx = h0
8 9 10 11 12
x [mm]
0.0
0.5
1.0
A
B
8 9 10 11 12
x [mm]
0.0
0.5
1.0
A
B
Figure 2.15: Snapshots of the damage pattern of the bar at different times simulated by the present
model
spalling test, the specimen exhibits a diffused damage pattern, and the present model is shown to
yield a mesh-independent simulation result. By contrast, the localized damage model grossly un-
derestimates the energy dissipation for large mesh size. These different numerical examples clearly
demonstrate the advantage of the present mechanism-based energy regularization method for FE
analysis of quasibrittle structures, which may exhibit complex damage patterns under different load-
ing scenarios.
Chapter 3
Rate and Size Effects on Strength
Distribution of Quasibrittle
Structures
Chapter 2 presented an energy regularization method that effectively mitigates mesh dependence
in continuum finite element (FE) simulations of quasibrittle fracture under general loading scenar-
ios. This energy regularization scheme was applied solely in deterministic analysis of quasibrittle
structures. As explained in Section 1.1, stochastic FE simulations of quasibrittle fracture can lead
to mesh dependent output probabilistic structural response. One potential means to mitigate the
mesh dependence issue in the stochastic setting is to consider that the probability distributions of
constitutive properties of finite elements depend on the mesh size. The underlying idea is to treat
each finite element as a material element of finite size. This implies that the sampling distribution
functions of constitutive parameters should depend on both the inherent randomness of material
properties as well as the failure mechanism of the corresponding element. As shown in Chapter 2,
quasibrittle structures can exhibit a complex failure pattern, including both localized and diffused
damage states, which evolve during loading. Therefore, it is important to understand the size effects
on quasibrittle structures under different loading configurations. For this reason, this chapter stud-
ies the combined rate and size effects on the strength statistics of quasibrittle structures through
stochastic discrete element model simulations of tensile failure of ceramic specimens. An analytical
probabilistic model is developed, that can capture the rate and size effects on the first and second
46
3.1. STOCHASTIC DISCRETE MODEL SIMULATIONS OF DYNAMIC FAILURE 47
statistical moments of structural strength, as provided by the stochastic discrete element computa-
tions. The proposed model has important implications for mitigating mesh dependence of stochastic
FE simulations of strain-softening continua.
3.1 Stochastic Discrete Model Simulations of Dynamic Failure
3.1.1 Model Description
In order to explore the underlying physical mechanisms leading to the size and rate dependent
mechanical and probabilistic structural response, we perform a series of stochastic discrete element
model simulations of dynamic fracture of ceramic specimens. In the computational model, the
specimen is represented by a set of interconnected discrete rigid bodies. The geometry of the model
is random to avoid any bias in the crack propagation direction. A set of nuclei is first randomly placed
in the domain, where the mutual distance of two adjacent nuclei is controlled to be approximately
equal to the grain size lmin of the material. Once these nuclei are created, Voronoi tessellation is used
to discretize the domain (Fig. 3.1a)). The location of the ith nuclei is denoted by xi. Each nucleus
has three translational and three rotational degrees of freedom. Due to the imposition of the length
scale lmin, each Voronoi body can be considered to represent a material grain. The contact surface
(henceforth referred to as facet) between the bodies represents the grain boundary (Fig. 3.1b)). The
overall nonlinear macroscopic behavior of the material is determined by the mesoscale constitutive
behavior of the facets.
In this study, the formulation of the constitutive behavior of the facet follows the lattice discrete
particle model (LDPM), which was originally developed for concrete materials [4,47–51]. Since this
study focuses primarily on macroscopic tensile failure, we adopt a simplified version of the LDPM,
which is briefly summarized here. The detailed description of the model can be found in [58].
The relative motion of two adjacent Voronoi bodies directly results in displacement jumps on the
facet. For two adjacent bodies i and j, the displacement jump vector ∆ij on their common facet
can be calculated as [58]
∆ij = Aj
uj
θj
−Ai
ui
θi
 (3.1)
where uk,θk (k = i, j) = displacement and rotation vectors of the nucleus of body k and Ak(k =
3.1. STOCHASTIC DISCRETE MODEL SIMULATIONS OF DYNAMIC FAILURE 48
a)
b)
Voronoi body
lij
n
m
l
Figure 3.1: Representation of the discrete element model: a) domain discretization, and b) Voronoi
body and facet
i, j) = transformation matrices, which can be written as
Ak =

1 0 0 0 xc3 − xk3 xk2 − xc2
0 1 0 xk3 − xc3 0 xk1 − xc1
0 0 1 xc2 − xk2 xk1 − xc1 0
 (3.2)
where xc1, xc2, xc3 = coordinates of the facet centroid, and xk1 , xk2 , xk3 = coordinates of the nucleus of
body k.
The components of the strain vector can be calculated as en = n ·∆ij/lij , el = l ·∆ij/lij , em =
m ·∆ij/lij , where ∆ij = displacement jump vector, lij = distance between the nuclei of rigid bodies
i, j, and n, l,m = unit vectors in the normal and two orthogonal tangential directions of the facet,
respectively. The corresponding traction vector transmitted over the facet can then be related to
3.1. STOCHASTIC DISCRETE MODEL SIMULATIONS OF DYNAMIC FAILURE 49
the strain vector by using damage mechanics, i.e.
tn
tl
tm
 = E0(1− ω)

en
αel
αem
 (3.3)
where E0 = elastic stiffness, ω = damage parameter, α represents the ratio of tangential to nor-
mal stiffness, and tn, tl, tm = components of the traction vector in the directions of n, l and m,
respectively. The damage parameter is expressed in terms of the equivalent stress and strain, i.e.
ω = 1− seq
E0eeq
(3.4)
seq = feq exp
(
K
feq
〈
χ− feq
E0
〉)
; eeq =
√
e2n + α(e
2
m + e
2
l ) (3.5)
For tension and shear dominated loading, the equivalent strength feq can be defined as [49,58]
feq = ft
 (η2 + µ2) sinψ −
√
(η2 − µ2)2 sin2 ψ + 4αη2 cos2 ψ
2µ2 sin2 ψ − 2α cos2 ψ
 (3.6)
where ft = tensile strength, fs = shear strength, ψ = tan−1
[
en/
√
α(e2m + e
2
l )
]
, η = fs/ft, and µ =
constant. The history variable χ is expressed as
χ =
eeqψ/ψ0 + emax (1− ψ/ψ0) ψ0 ≤ ψ < 0emax ψ ≥ 0 (3.7)
where emax =
√
max (e2n) + αmax (e
2
m + e
2
l ), in which the maxima is calculated for the entire loading
history, and ψ0 = − arctan (ζ
√
α). For the present numerical simulations, the facets mainly expe-
rience tensile-shear damage. Therefore, it is not necessary to cover the full scenario of compressive
damage.
Finally, parameter K represents the initial slope of the nonlinear branch of the constitutive
law [49], which is given by
K = −Kt
[
1−
(
ψ − pi/2
ψ0 − pi/2
)nt]
; nt =
ln [Kt/(Kt −Ks)]
ln (1− 2ψ0/pi) (3.8)
where Kt,Ks are the values of K for pure tensile and shear loading, respectively. Kt and Ks are
3.1. STOCHASTIC DISCRETE MODEL SIMULATIONS OF DYNAMIC FAILURE 50
expressed as
Kt =
2E0
lt/lij − 1 Ks =
2αE0
ls/lij − 1 (3.9)
where lt = E0Gt/f2t , ls = E0Gs/f2s , and Gt, Gs = mode I and II fracture energies. It should be
pointed out that Eq. (3.9) involves two material length scales lt and ls so that the energy expended
for the fracture of a unit area of facet is independent of the distance between the nuclei. Once the
traction vectors on the facets are determined, the forces and moments of the nuclei can be calculated
based on the principle of virtual work [58].
It is worthwhile to mention that the present constitutive model is itself rate-independent, and
therefore the rate effect on macroscopic behavior of the specimen is largely due to inertia and its
influence on the interaction of facet failures. Such a modeling approach has been adopted in several
previous discrete simulations of dynamic fracture of ceramic materials [54, 87, 88, 138], which were
able to capture the rate-dependent failure behavior reasonably well. Meanwhile, it has also been
shown that the rate dependence of the constitutive behavior plays a secondary role in the overall
failure behavior under high strain-rate loading [112,113], which is of interest to the present study.
The aforementioned constitutive model contains the following parameters: 1) elastic modulus
E0 in the normal direction, 2) tangential to normal stiffness ratio α, 3) mesoscopic tensile and
shear strengths ft, fs, 4) mesoscopic mode I and II fracture energies Gt, Gs, and 5) constants µ
and ζ. It is well known that the properties of quasibrittle materials usually exhibit a considerable
degree of spatial variability. For studying the failure behavior, the most relevant parameters are the
strengths and fracture energies [69]. In this study, we characterize the spatial random distribution
of mesoscopic strengths and fracture energies by assigning a single random field h(x):
ft(x) = f¯th(x); fs(x) = f¯sh(x) (3.10)
Gt(x) = G¯t [h(x)]
2
; Gs(x) = G¯s [h(x)]
2 (3.11)
where the overbar denotes the mean of the corresponding random variable. Eqs. (3.10) and (3.11)
imply that, on each facet, the material length scales lt and ls are deterministic constants. The
underlying probability distribution function governing the random field h(x) follows the Gauss-
Weibull grafted distribution (Eqs. (3.16a) and (3.16b)) with a mean value equal to one. The spatial
autocorrelation of the field h(x) is characterized by a Gaussian function, i.e.
R(δx) = exp[−(δx/la)2] (3.12)
3.1. STOCHASTIC DISCRETE MODEL SIMULATIONS OF DYNAMIC FAILURE 51
where δx = |x − x′| and la = autocorrelation length. The generation of the random field h(x)
consists of three steps [59]: 1) a standard Gaussian random field ĥ(x) is generated on a regular
square grid with a spacing of la/4, 2) the value of ĥ(x) on each facet is determined by using the
optimal linear estimation method [95], and 3) the standard Gaussian field generated on the facets
is converted to the Gauss-Weibull random field h(x) by using the isoprobabilistic transformation.
Fig. 3.2 presents a typical realization of the random field h(x).
Figure 3.2: A typical realization of the random field h(x)
The aforementioned computational model is numerically implemented with the implicit Newmark
method [7, 109]. The equation of motion is solved incrementally in an implicit scheme, which is
unconditionally stable. The mass matrix is constructed by using the full inertia tensor of the
Voronoi bodies. The moment of inertia of these polyhedral bodies are calculated by dividing them
into tetrahedra [130]. Though a stable numerical solution scheme is chosen, a small time step is
needed to capture the dynamic failure of the facet. For each simulation case, several trial simulations
are performed to determine a desirable time step, which yields a consistent result.
3.1. STOCHASTIC DISCRETE MODEL SIMULATIONS OF DYNAMIC FAILURE 52
𝑛
a) b)
Figure 3.3: A rectangular specimen loaded by a prescribed strain rate: a) schematic of the specimen,
and b) weakest link model representation
3.1.2 Simulations of Dynamic Tensile Failure of Aluminum Nitride Spec-
imens
The aforementioned stochastic discrete element computational model is applied to simulate the
dynamic tensile failure of aluminum nitride (AlN) specimens. In the simulations, we consider square
specimens of different in-plane sizes D = L = 50 µm, 100 µm, 200 µm, 400 µm, 800 µm (Fig. 3.3a)),
whereas the out-of-plane thickness is set constant b = 10 µm. The average in-plane grain size lmin
of AlN material is taken as 6 µm [54]. Each specimen is subjected to a constant strain-rate loading,
which is applied by imposing a non-uniform velocity field described as vx1 = x1˙ (x1 denotes the
horizontal position of each nucleus). This velocity gradient is applied as an initial condition (t = 0),
whereas at subsequent times (t > 0) only the velocities of the sections at x1 = 0 and x1 = L are
controlled (equal to 0 and ˙L respectively). In order to investigate the rate-dependent behavior, we
consider a wide range of applied strain rates, i.e. ˙ = 1, 1000, 2500, 5000, 7500, 104, 3× 104, 5×
104, 105, 2× 105/s.
We choose the following mesoscale material parameters for AlN: E0 = 530GPa, α = 0.17,
f¯t = 150MPa, f¯s = 3f¯t, G¯t = 2 Jm−2, G¯s = 16G¯t, µ = 0.2, and ζ = 0.95. The elastic parameters
(E0 and α) are determined to match the macroscopic elastic response of AlN reported in [54]. The
inelastic mesoscale material parameters are chosen so that the model predicts the quasi-static tensile
strength of a 50 µm square AlN specimen being about 120 MPa, which is similar to the published
results [54]. It is admitted that mixed-mode loading scenarios will be needed to better calibrate
3.1. STOCHASTIC DISCRETE MODEL SIMULATIONS OF DYNAMIC FAILURE 53
the coupling between tensile and shear damage. However, such a detailed calibration procedure
is not necessary for the present purpose. For the random field h(x), we consider that it has an
autocorrelation length la of 24µm, and the underlying Gauss-Weibull distribution function has a
mean value of one, a coefficient of variation of 20%, a Weibull modulus of 30 and a grafting probability
of 5×10−4.
For each specimen size and strain rate, about 70 realizations of different random fields h(x)
and mesostructures are used to determine the mean and variance of the peak load capacity of the
specimen. As mentioned earlier, the peak load capacity of the specimen is expressed in terms of the
nominal strength, which in this case is equal to the maximum value of the nominal stress σa along
the x1−direction (Fig. 3.3a)). For the present discrete system, we first calculate the fabric stress
tensor of a single Voronoi body k as in [5, 77,122]:
σ¯k =
1
Vk
nk∑
p=1
tp ⊗ (xcp − xk)Ap (3.13)
where Vk = volume of the kth body, nk = number of facets of the body, tp = traction vector on the
pth facet expressed in the global coordinate system, xcp = position vector of the centroid of the pth
facet, xk = position vector of the nucleus of the body, and Ap = surface area of the pth facet. The
nominal stress σa can then be calculated as a volume average of the fabric stress components in the
x1−direction, i.e.:
σa =
1
V
N∑
k=1
Vkσ¯
k
11 (3.14)
where V = volume of the specimen, and N = number of Voronoi bodies in the specimen.
Based on the simulations, we obtain a set of random responses of the relationship between the
nominal stress σa and the average strain a, where a = u/D and u = displacement measured at the
free edge of the specimen. Fig. 3.4 presents the simulated average σa–a curves for all specimen sizes
and strain rates. It can be seen that, for a given specimen size, the average nominal strength σ¯N ,
which is the maximum value of the average nominal stress, increases with the applied strain rate.
This rate enhancement is primarily due to the inertia effect. Meanwhile, we also observe that, as
the strain rate increases, the post-peak softening behavior becomes less pronounced, which indicates
that the specimen exhibits a more ductile behavior. This rate dependence can directly be observed
from the simulated damage pattern at the peak load. As a demonstration, Fig. 3.5 shows the damage
patterns of specimens of D = 400 µm at the peak load for different applied strain rates. It is seen
that, for a given specimen size, the damage is fairly localized at low strain rates while the specimen
3.1. STOCHASTIC DISCRETE MODEL SIMULATIONS OF DYNAMIC FAILURE 54
Figure 3.4: Simulated average nominal stress-strain curves for different specimen sizes and strain
rates
3.1. STOCHASTIC DISCRETE MODEL SIMULATIONS OF DYNAMIC FAILURE 55
exhibits a diffused damage pattern as the strain rate increases.
Figure 3.5: Simulated damage patterns of specimen of D = 400µm at the peak load
It should be pointed out that the aforementioned brittle-to-ductile transition is based on the
current simulation setup with the prescribed initial velocity field. Recent computational studies of
ceramics and concrete specimens with the same setup showed similar rate-dependent failure patterns
[54,64]. Several experiments on the rate-dependent behavior of concrete also showed that the tensile
stress-strain curve exhibits a more gentle post-peak softening behavior at high strain rates [72,133].
However, there is experimental evidence showing that the compressive stress-strain response of
brittle composites and ceramics could become more brittle at higher strain rates [35, 84, 96]. The
setups of these experiments are different from the present simulations. We may expect that the
rate dependence of failure behavior of quasibrittle materials could depend strongly on the specimen
geometry and test configuration. This is consistent with the well accepted fact that, for quasibrittle
structures, brittleness is influenced not only by the material properties, but also by the structure
geometry and loading configuration. [8, 15,19].
3.2. RATE-DEPENDENT FINITE WEAKEST LINK MODEL 56
We can also compare the simulated post-peak responses of specimens of different sizes for a given
strain rate. It is seen that, at low strain rates, the slope of the post-peak softening curve becomes
much steeper as the specimen size increases. This implies that the specimen experiences a more
brittle failure as the specimen size increases. Such a size dependent failure behavior at low strain
rates has been well documented for many quasibrittle materials, such as concrete, rock, composites,
ceramics, etc. [8, 19]. As the strain rate increases, the size dependence of the post-peak behavior
starts to diminish. This is consistent with the present simulation results, which show that, at high
strain rates, all the specimens considered in this study exhibit a diffused cracking pattern.
3.2 Rate-Dependent Finite Weakest Link Model
3.2.1 Model Formulation
The simulation results presented in Section 3.1.2 clearly indicate the rate dependence of the failure
behavior of quasibrittle structures. Accurate probabilistic modeling of failure requires the assump-
tions of the corresponding mathematical models to be consistent with the actual failure mechanism
of a structure, therefore the rate effect, through its influence on the failure mechanism, can have
profound implications in the development of such models. In this section, we develop an analytical
probabilistic model that can describe the combined rate and size effects on the strength distribution
of quasibrittle structures. The statistics of static strength of quasibrittle structures of positive ge-
ometry have been described by the finite weakest link model [18, 90], which considers the structure
to be statistically represented by a chain of representative volume elements (RVEs), and assumes
that structural failure will occur as soon as a RVE fails. This assumption is realistic for quasibrittle
structures of positive geometry, since these structures fail right at the initiation of a macroscopic
crack. At the same time, the finite weakest link model takes into account the influence of the finite
size fracture process zone (FPZ) on the failure behavior of quasibrittle structures, and is thus more
general than Weibull’s statistical theory [134], which applies to brittle failure only.
Since in this study we focus on quasibrittle structures, we essentially build upon the finite weakest
link model, which will provide the limiting strength distribution for the case of very low strain rates
(quasi-static), and extend this model so as to account for the change in the failure mechanism with an
increasing strain rate. Considering a set of geometrically similar rectangular quasibrittle specimens,
which are loaded under displacement control with a prescribed strain rate (Fig. 3.3a)), we define
the nominal tensile strength, σN , as the maximum average stress that the specimen can sustain.
In the weakest link representation of the structure (Fig. 3.3b)), we consider that the RVE size is
3.2. RATE-DEPENDENT FINITE WEAKEST LINK MODEL 57
taken to be larger than the autocorrelation length of the underlying random strength field, thus
the RVE strengths can be modeled as statistically independent random variables with cumulative
distribution function (cdf) P1. In this case, the structural failure probability Pf (under quasi-static
loading conditions) can be expressed as follows:
Pf (σN ) = 1− [1− P1(σN )]n (3.15)
where n is the number of RVEs in the structure, which, for the case of two-dimensional scaling, is
equal to n = A/l20s, A = DL = area of the specimen, and l0s = RVE size under quasi-static loading.
a)
b)
Figure 3.6: Modeling of strength distribution of a RVE: a) hierarchical model of static RVE, and b)
fiber-bundle model of RVE under high strain-rate loading
In recent studies on statistics of static strength [14, 15, 18, 90], it was shown that the RVE can
be statistically modeled by a hierarchical model, which consists of a combination of bundles and
chains (Fig. 3.6a)). This hierarchical model represents damage localization and load redistribution
mechanisms at different material scales. The basic element in the model represents a nanoscale
structure, either an atomic lattice or a disordered system of nano-particles. The failure statistics
of the nanoscale structure was derived based on the transition rate theory and atomistic fracture
mechanics [14,15,90]. Based on the hierarchical model, it was shown that the probability distribution
3.2. RATE-DEPENDENT FINITE WEAKEST LINK MODEL 58
of the static strength of the RVE can be approximated by a Gaussian distribution onto which a
Weibull tail (or, equivalently, a power-law tail) is grafted at a fairly low probability. Mathematically,
the cdf of the RVE strength can be written as
P1(σN ) =

1− exp (−〈σN/s0〉m0) ≈ 〈σN/s0〉m0 (σN ≤ σ0gr) (3.16a)
P 0gr +
rf0
δG0
√
2pi
∫ σN
σ0gr
e−(x−µG0)
2/2δ2G0dx (σN > σ0gr) (3.16b)
where m0 and s0 are the Weibull modulus and the scale parameter of the Weibull tail, respectively,
and µG0 and δG0 are the mean and standard deviation of the Gaussian core if considered extended
to −∞; rf0 is a scaling parameter required to normalize the grafted cdf such that P1(∞) = 1,
σ0gr = grafting stress, and P 0gr = (σ0gr/s0)m = grafting probability. Note that here the subscript and
superscript “0” of all the parameters indicate the quasi-static case. The continuity of the probability
density function (pdf) at the grafting point requires that [dP1/dσN ]σ0+gr = [dP1/dσN ]σ0−gr . Therefore,
there are in total four independent parameters to define the grafted distribution of the RVE strength.
To incorporate the strain rate effect into the probabilistic model of structural failure, we start
by modifying the weakest link model (Eq. (3.15)) in order to account for the rate dependence of
the RVE size, denoted as l0 (˙). Several recent studies [54,64] have investigated the dynamic tensile
fracture behavior of ceramic and concrete materials by considering specimens similar to that depicted
in Fig. 3.3a). These simulation results consistently showed that, with an increasing strain rate, the
damage pattern becomes more diffusive, in agreement with the results of the present stochastic
discrete element model computations described in Section 3.1.2. Therefore, we expect that the RVE
size l0 (˙) will be an increasing function of the strain rate. Besides, the strain rate can also affect the
statistical parameters of the distribution function P1 = P1(σN , ˙). Based on the foregoing discussion,
we can express the structural failure probability Pf (σN , ˙) by a rate-dependent weakest-link model:
Pf (σN , ˙) = 1− [1− P1(σN , ˙)] [1− P1(σN , ˙)]〈n(˙)−1〉 (3.17)
where P1(σN , ˙) = strength distribution of either one RVE or the entire specimen, whichever has
a smaller size, n(˙) = A/[l0(˙)]2, 〈x〉 = max(x, 0) = Macaulay bracket, and l0(˙) = RVE size
under dynamic loading conditions. Since the RVE size increases with the rate, for a given specimen
size, an increase in the applied strain rate could possibly make the RVE larger than the specimen
size, i.e. n(˙) < 1. In this case, Eq. (3.17) reduces to Pf (σN , ˙) = P1(σN , ˙). In other words,
P1(σN , ˙) becomes the strength distribution of the entire specimen. This implies that the weakest-
link representation of structural failure vanishes since there is no damage localization in the specimen.
3.2. RATE-DEPENDENT FINITE WEAKEST LINK MODEL 59
In order to determine the rate dependence of the distribution function P1(σN , ˙), we propose a
fiber-bundle model for the RVE under dynamic loading, which consists of dnbe fibers, where dnbe =
the least integer that is greater than or equal to nb, and nb represents the equivalent number of
static RVEs in the dynamic RVE or in the entire structure, whichever is smaller:
nb(˙) = min[n(˙), 1]
[
l0(˙)
l0s
]2
(3.18)
where, it should be noted that, in the present context, we refer to an RVE of size equal to l0s as
the static RVE. The fiber-bundle statistical model represents a structure as a series of elements,
called fibers, which are loaded in parallel (Fig. 3.6b)). In contrast to the weakest link model, the
fiber-bundle model represents structural failure under distributed cracking, since, as soon as a fiber
fails, its load is redistributed among other surviving fibers, preventing global failure. Thus, it is
considered to be able to realistically represent the damage redistribution mechanism of the RVE
under high strain rates (or of the entire specimen if n(˙) < 1), given the diffused damage pattern
observed (Section 3.1.2).
Each fiber in the bundle model of the dynamic RVE represents a static RVE. For quasibrittle
materials, the mechanical behavior of the static RVE will exhibit some degree of strain-softening.
On the other hand, we also note that the static RVE should behave in a quasi-plastic manner since it
contains only a few material heterogeneities [18]. This quasi-plastic failure behavior is also manifested
by the dominant Gaussian distribution of RVE strength as indicated by Eq. (3.16b) [15, 18, 90].
Therefore, the stress-strain curve of the static RVE must have a gentle softening behavior even
under quasi-static loading. As the applied strain rate increases, the static RVE would experience
more diffusive micro-cracking, which alleviates the degree of strain softening [54, 64, 72, 133]. For
this reason, we consider that the elements in the fiber-bundle model exhibit a plastic stress-strain
behavior. Besides, the individual fiber strengths can be treated as statistically independent random
variables, since the autocorrelation length of the random strength field is considered to be smaller
than the size of the static RVE. The plastic bundle model of the RVE under dynamic loading
indicates that the bundle strength is a weighted sum of the strengths of the individual fibers, which
leads to the following relationship:
σbN =
1
nb
dnbe−1∑
i=1
σi + (nb + 1− dnbe)σdnbe
 (3.19)
where σi (i = 1, ..., dnbe) are the random strengths of the individual fibers. Since nb may not be
3.2. RATE-DEPENDENT FINITE WEAKEST LINK MODEL 60
an integer, in writing Eq. (3.19) we have considered that the bundle contains dnbe fibers, where
dnbe − 1 fibers have the same cross-sectional area A0, and the remaining fiber has a cross-sectional
area of (nb + 1 − dnbe)A0. The strength distribution of each fiber representing a static RVE is
predominantly Gaussian with a very short power-law tail (Eqs. (3.16a) and (3.16b)). Previous
studies have shown that, for such a strength distribution of fibers, the strength distribution of the
bundle will also contain a power-law tail, and meanwhile the core of the strength distribution of
the bundle is Gaussian [18]. Therefore, it is clear that the previously proposed grafted distribution
function (Eqs. (3.16a) and (3.16b)) provides the correct functional form for the strength distribution
function P1(σN , ˙), though the values of the statistical parameters must now depend on the applied
strain rate.
As mentioned earlier, the grafted distribution function of the RVE strength can be determined
by any four statistical parameters in Eqs. (3.16a) and (3.16b). In this study, we incorporate the
strain rate dependence into the strength distribution function P1(σN , ˙) through the following four
parameters:
• The Weibull modulus m measures the exponent of the power-law tail of the bundle. By using
the series expansion method, it was shown that the power-law exponent of the bundle is equal
to the sum of the power-law exponents of the tail distributions of the individual fibers [18]. In
fact, this property also holds for brittle and softening bundles, where the randomness of the
stress-strain curve of the fibers is described by an affine transformation [15,90,126]. Meanwhile,
the strain rate may also affect the Weibull modulus of the static RVE itself, which is described
by an empirical function f(˙) to be determined later. Following these considerations, we can
express the Weibull modulus of the distribution function P1(σN , ˙) as
m(˙) = m0f(˙)dnb(˙)e (3.20)
• The grafting probability Pgr determines the extent of the power-law tail. Previous studies have
shown that, as the number of fibers increases, the power-law tail of the strength distribution
of the bundle shortens at the rate of (Pgr,f/dnbe)dnbe (Pgr,f is the grafting probability of the
strength distribution of an individual fiber) [18]. By further considering the potential influence
of the strain rate on the grafting probability of the static RVE, the rate-dependent grafting
probability is written as
Pgr(˙) ≈
[
P 0grg(˙)/dnb(˙)e
]dnb(˙)e (3.21)
3.2. RATE-DEPENDENT FINITE WEAKEST LINK MODEL 61
where the term P 0grg(˙) describes the rate-dependent grafting probability of the static RVE.
Based on Eq. (3.21), it is clear that the Weibull tail of the RVE strength is shortened drastically
as nb increases. This means that, when the strain rate is sufficiently high, Eq. (3.20) becomes
practically unimportant.
• Parameter µG is approximately equal to the mean strength of the plastic bundle. The power-
law tail of the strength distribution of each fiber is too short to influence the mean behavior.
Based on the statistics of the weighted sum of independent Gaussian variables (Eq. (3.19)),
it is clear that µG is equal to the mean strength of each fiber, which leads to the following
expression
µG(˙) = µG0p(˙) (3.22)
where p(˙) is the rate enhancement function of the mean strength of the static RVE.
• Parameter δG measures the standard deviation of the bundle strength since the power-law tail
of the distribution function has a minimal influence on the second moment of statistics. The
statistics of the weighted sum of independent Gaussian variables gives
δG(˙) =
δG0q(˙)
nb
[
dnbe − 1 + (nb + 1− dnbe)2
]1/2
(3.23)
where function q(˙) describes the influence of the strain rate on the standard deviation of the
strength of the static RVE.
By substituting Eqs. (3.20)-(3.23) into the grafted distribution function (Eqs. (3.16a) and (3.16b)),
we obtain the rate-dependent strength distribution function P1(σN , ˙). Based on the rate-dependent
weakest-link model (Eq. (3.17)), we can then calculate the mean dynamic strength of the specimen
as:
σ¯N (˙) =
∫ 1
0
σNdPf (σN , ˙) =
∫ ∞
0
[1− P1(σN , ˙)] [1− P1(σN , ˙)]〈n(˙)−1〉 dσN (3.24)
Meanwhile, the model also yields the rate-dependent second moment of strength statistics. For
instance, the standard deviation of the dynamic strength can be calculated as:
δσN (˙) =
√∫ 1
0
σ2NdPf (σN , ˙)− [σ¯N (˙)]2 (3.25)
By considering specimens of different sizes, we can determine the effect of the specimen size on the
statistics of dynamic strength. Together with the aforementioned rate dependence, Eqs. (3.24) and
3.2. RATE-DEPENDENT FINITE WEAKEST LINK MODEL 62
(3.25) describe the combined size and rate effects on the first and second moments of the statistics
of nominal tensile strength of the specimen.
3.2.2 Numerical Verification
We now use the rate-dependent finite weakest-link model (Eqs. (3.24) and (3.25)) to compare with
the rate and size effects on the mean and standard deviation of structural strength as simulated by
the stochastic discrete element model. Fig. 3.7 presents the size effect curves of the mean nominal
strength at different applied strain rates. The most salient feature is that the mean size effect
diminishes as the strain rate increases. At high strain rates, the entire specimen is occupied by
diffused cracking as the peak load is reached (Fig. 3.5). In such a case, there does not exist any
characteristic length scale governing the nominal strength. Since the failure is ductile (or quasi-
plastic), the size effect on the nominal strength must be absent [8,15]. As the strain rate decreases,
the damage pattern transitions from diffused cracking to a localized cracking zone. The size of the
localized damage zone represents a characteristic length scale, which leads to a non-power law form
of the mean size effect curve. As seen from the predictions of the rate-dependent finite weakest link
model in Fig. 3.7, the influence of the strain rate on the mean size effect behavior can correctly be
captured.
Based on this fitting, we determine the rate effect on the RVE size and on the statistical pa-
rameters of the strength distribution of the static RVE. From Fig. 3.8 it is clear that the RVE size
l0 increases with the strain rate. More specifically, the RVE size is almost constant for relatively
low strain rates, and increases mildly over a range of intermediate strain rates 3× 104 − 5× 104/s.
When the strain rate exceeds 5 × 104/s, the RVE size increases significantly with the strain rate.
In fact, for the high strain-rate regime studied here, the analysis indicates that all the specimens
are smaller than the dynamic RVE. Therefore, the actual RVE size is not known. At high strain
rates, the RVE can become a purely mathematical concept if the specimen does not exhibit damage
localization. In this case, we may set l0(˙) to be a very large number. However, introducing the
concept of the RVE for the present theoretical framework is essential because it captures the rate de-
pendence of the failure statistics, which transitions from a weakest-link model signifying the damage
localization mechanism to a fiber-bundle model representing the more diffused damage mechanism
with increasing strain rate.
The concept of the RVE is central to the present weakest-link model. In a recent study [93], it
has been shown that the RVE size can be best determined from the optimum fitting of the mean size
curve of nominal structural strength, and a relationship between the RVE size and basic material
3.2. RATE-DEPENDENT FINITE WEAKEST LINK MODEL 63
Simulations
Model
✏˙ = 1/s
✏˙ = 2500/s ✏˙ = 5000/s✏˙ = 1000/s
✏˙ = 7500/s ✏˙ = 104/s ✏˙ = 3⇥ 104/s
✏˙ = 5⇥ 104/s ✏˙ = 2⇥ 105/s✏˙ = 105/s
Figure 3.7: Simulated mean size effect curves of nominal strength at different strain rates and the
optimum fits by the rate-dependent finite weakest-link model
length scales (Irwin’s characteristic length and crack band width), was established for the case of
quasi-static loading. Based on the results of this recent study, we expect that the RVE size will
increase with the average grain size. Meanwhile, the RVE size will also increase with an increase
in the mesoscale fracture energies. This implies that the RVE size is dependent on the loading
configuration. For example, if the applied loading induces more shear-dominant fractures along the
grain boundaries, based on the material properties used in the present simulation the specimen will
exhibit a more ductile failure behavior and therefore the RVE size will increase.
Fig. 3.9 presents the influence of the applied strain rate on the strength statistics of a RVE with
size equal to the corresponding static value l0s (herein referred as the static RVE). We observe that
Weibull modulus increases with the applied strain rate (i.e. f(˙) > 1 in Eq. (3.20)) and the grafting
3.2. RATE-DEPENDENT FINITE WEAKEST LINK MODEL 64
RVE size is 
greater than the 
specimen size
103 104 3⇥ 104 5⇥ 104
Figure 3.8: Effect of the strain rate on the RVE size
probability decreases with the strain rate (i.e. g(˙) < 1 in Eq. (3.21)). It is noted that, when the
strain rate is higher than 7500/s, the power-law tail of the strength distribution of the static RVE
is too short to be determined for the range of specimen sizes considered here. The observed rate
dependence of the strength statistics of the static RVE can be attributed to the fact that, as the strain
rate increases, the static RVE experiences more dense micro-cracking. In the context of the present
stochastic discrete element model, microcracking is represented by the failure of grain boundaries
(facets). The increase in the Weibull modulus can be qualitatively explained by the fishnet statistics
model (Section 1.3.3), which was recently developed to analyze a fishnet structure for predicting the
probability distribution of static strength of nacreous imbricated lamellar materials. Though the
specimens dealt in the present study are very different from the fishnet structure, due to the similar
failure mechanism which is described by diffuse microcracking, the same mathematical framework
can be applied to both cases.
Following the fishnet model [99], the strength distribution of a static RVE can be calculated as
Fs(σN ) = 1− [Ps,0(σN ) + Ps,1(σN ) + · · ·+ Ps,N (σN )] (3.26)
3.2. RATE-DEPENDENT FINITE WEAKEST LINK MODEL 65
Figure 3.9: Rate dependence of the statistical parameters of P1(σN )
where Ps,r = survival probability when r number of facets have experienced damage at stress level
σN applied to the static RVE, and N = total number of facets. For obtaining an approximation
to Eq. (3.26), we just need to retain several dominant terms (r = 0, ..., k), which contribute to the
overall failure probability. The aforementioned rate effect on microcracking density indicates that
the number of damaged facets k increases with the applied strain rate. To calculate the individual
probability component Ps,r of Eq. (3.26), one would need to determine the load redistribution
among the facets as damage proceeds. This requires performing the nonlinear stochastic discrete
element simulations, which is computationally intensive. Nevertheless, the recent results of the
fishnet analysis indicate that the exponent of the power-law tail of the strength distribution of the
system increases in proportion to the dominant number of damaged facets (i.e. k in Eq. (3.26)) [99].
This explains the observed in the current study increase in the Weibull modulus with an increasing
strain rate.
3.2. RATE-DEPENDENT FINITE WEAKEST LINK MODEL 66
The decrease in the grafting probability indicates that the strength distribution of the static
RVE approaches the Gaussian distribution. This can be physically related to the increasing number
of damaged facets, which implies that the failure process becomes more ductile. Since more facets
contribute to the failure of the static RVE, we expect that the standard deviation would decrease
considerably. However, it is also noted that the random mesostructure causes a considerable vari-
ability of the stress field, and this variability is further enhanced by the inertia effect. Therefore, the
stress field becomes more random at high strain rates, which leads to a larger standard deviation of
the overall nominal strength. This effect counteracts the aforementioned decrease in the standard
deviation due to the increasing number of damaged facets. As shown in Fig. 3.9, the calibration
results show that the overall standard deviation of the strength of the static RVE decreases mildly
with an increasing strain rate.
The calibrated rate-dependent weakest-link model is now used to predict the size effect curves of
the standard deviation for different applied strain rates. Fig. 3.10 shows that the model prediction
agrees well with the simulation results. It is seen that, in contrast to the diminishing size effect on
the mean nominal strength at high strain rates, the size effect on the standard deviation becomes
more pronounced as the strain rate increases. At low strain rates, the standard deviation exhibits
an intricate size effect ending with a power-law asymptote at the large-size limit. This theoretical
power-law asymptote is attributed to the fact that, according to the weakest-link model, at low strain
rates the strength distribution of large-size specimens would approach the Weibull distribution, i.e.
Pf (σN ) = 1− exp[−(D/l0)2(σN/s0)m]. The corresponding standard deviation of σN can be written
as
δσN = s0
(
l0
D
)2/m√
Γ
(
1 +
2
m
)
− Γ2
(
1 +
1
m
)
(3.27)
where Γ(x) = Eulerian gamma function. Eq. (3.27) indicates that, at low strain rates, the asymptotic
size effect on the standard deviation is very weak since m is usually a large number.
For the high strain rates considered in this study (˙ > 105/s), the fitting of the mean size effect
curve indicates that the specimen size is smaller than the RVE size. According to Eq. (3.18), we
have nb = (D/l0s)2. Therefore, the statistics of the nominal strength of the specimen is modeled by
a plastic bundle, and the scaling of the standard deviation is given by Eq. (3.23). It can be easily
shown that, as nb increases (i.e. nb > 4), Eq. (3.23) is approximately equivalent to:
δG(˙) =
δG0q(˙)√
nb
= δG0q(˙)l0sD
−1 (3.28)
We note that the size effect expressed by Eq. (3.28) is much stronger than that predicted by the
3.2. RATE-DEPENDENT FINITE WEAKEST LINK MODEL 67
Simulations
Model
N / D1 N / D1
✏˙ = 2500/s ✏˙ = 5000/s
✏˙ = 7500/s ✏˙ = 104/s
✏˙ = 3⇥ 104/s
✏˙ = 1/s
✏˙ = 1000/s
✏˙ = 5⇥ 104/s ✏˙ = 2⇥ 105/s✏˙ = 105/s
Figure 3.10: Simulated size effects on the standard deviation of nominal strength at different strain
rates and its comparison with the rate-dependent finite weakest-link model
Weibull distribution. From Fig. 3.10, we see that Eq. (3.28) agrees well with the simulated size effect
on the standard deviation at high strain rates.
The foregoing analysis elucidates the size and rate dependence of the functional form of the
strength distribution. At low strain rates, the strength distribution of the specimen transitions from
a Gaussian distribution modified by a far-left Weibull tail to a purely Weibull distribution as the
specimen size increases. At increasing strain rates, for the specimens considered in this study, the
strength distribution becomes predominantly Gaussian with a size independent mean value which
increases with the rate, but with a standard deviation that decreases strongly with the specimen
size. This transitional behavior is well captured by the proposed rate-dependent weakest-link model.
3.3. CONCLUSIONS 68
3.3 Conclusions
This study shows that the strain rate has a profound influence on the scaling of the nominal tensile
strength of quasibrittle structures. This rate dependence arises from the influence of the strain rate
on the failure behavior of the structure. The discrete element model simulations show that the
structure exhibits a more diffused cracking pattern as the strain rate increases, which explains the
observed diminishing mean size effect at high strain rates. Meanwhile, the rate-dependent failure
behavior also leads to different scaling behavior of the standard deviation of the nominal strength
at different strain rates. It is shown that, for a given strain rate, the standard deviation decreases
with an increasing specimen size. At high strain rates, this decreasing trend becomes stronger and
it exhibits a power-law scaling behavior. The result also indicates, that, for a given specimen size,
an increasing strain rate leads to a reduction in the standard deviation of nominal strength.
It is shown that the simulated rate and size effects on the mean and standard deviation of
the nominal tensile strength can be well captured by a rate-dependent finite weakest-link model.
The model employs a rate-dependent length scale, which physically represents the transition from
damage localization to diffused damage with an increasing strain rate. The model directly predicts
the combined rate and size effects on the probability distribution of nominal strength. At low strain
rates, the strength distribution varies from a predominantly Gaussian distribution to a Weibull
distribution as the specimen size increases, whereas at high strain rates, the strength distribution
is Gaussian for a significantly large range of specimen sizes. The developed rate-dependent finite
weakest-link model provides an analytical tool to model the rate and size effects on the probability
distribution of nominal tensile strength of quasibrittle structures, which could be implemented into
stochastic FE simulations of quasibrittle fracture under various strain rates. More specifically,
considering that each finite element represents a material element of finite size, the present model
indicates that the input probability distribution of tensile strength of finite elements should depend
on the failure mechanism, which was shown to be influenced by both the element size and applied
strain rate. The above is an essential consideration for mitigating mesh dependence of the output
probabilistic response in continuum finite element simulations of quasibrittle structures with random
material properties.
Chapter 4
Stochastic Numerical Simulations of
Strain-Softening Continua
Chapter 3 investigated the size effects on the statistics of structural strength of quasibrittle struc-
tures through stochastic discrete element model computations. The simulated set of quasibrittle
specimens of different sizes exhibited different failure mechanisms ranging from fairly localized to
purely diffused, which were triggered by application of various uniform tensile strain rates, from
quasi-static up to the high strain-rate regime. It is clear that the observed failure mechanism gov-
erns the mathematical formulation of probabilistic models of structural failure. For instance, the
well-established weakest-link and fiber-bundle models statistically represent the failure mechanisms
of damage localization and distributed damage respectively. The developed rate-dependent finite
weakest-link model transitions from a weakest-link to a fiber-bundle model representation of struc-
tural failure as the damage pattern changes from localized damage to distributed cracking with an
increasing strain rate. Such transition was achieved by introducing rate dependence of the statistical
parameters of the strength distribution. The model was shown to accurately capture the simulated
size and rate-dependent mean and standard deviation of structural strength.
The aforementioned study has important implications for mitigating the mesh sensitivity of the
output probabilistic structural response in stochastic FE simulations of strain-softening continua,
which is the topic of the present chapter. The stochasticity of peak load capacity arises from the
spatial variability of element tensile strength and fracture energy. Focus is placed on the case when
the finite element mesh size is larger than the crack band width of the material, as well as the
autocorrelation length of the underlying random field of constitutive properties. Considering that
69
4.1. MODEL FORMULATION 70
each finite element represents a material element of finite size, the sampling probability distribu-
tions of the input tensile strength and fracture energy are related to the failure mechanism of the
corresponding element. This is achieved by introducing a dependence of the distribution functions
of apparent fracture properties on proposed numerical parameters, some of which were introduced
in Chapter 2, that are able to capture the actual evolution of damage pattern including the stages
of damage initiation, localization and propagation. Therefore, the present formulation of the prob-
ability distributions of fracture properties and the mechanism-based energy regularization scheme
of Chapter 2 are interconnected and are always consistent with the actual failure mechanism of the
structure. Besides, the explicit dependence of the present probabilistic model on the damage mech-
anism makes it applicable to general loading scenarios, including rate effects, and thus resulting in
a more general approach than the rate-dependent weakest link model of Chapter 3, which was also
published in [91]. The model is applied to simulate the stochastic failure behavior of quasibrittle
structures of different geometries, subjected to different loading configurations. The numerical ex-
amples are selected such that various failure mechanisms are triggered and therefore the generality
of the method can be tested. It is shown that the proposed model is efficient in mitigating mesh
dependence of the output strength statistics in all considered cases, thus the present study provides
a general framework for addressing the issue of mesh dependence in stochastic FE simulations of
quasibrittle failure.
4.1 Model Formulation
4.1.1 Probabilistic Modeling of Element Fracture Properties
The constitutive law to which our probabilistic model is applied was formulated in Chapter 2 and
is briefly described here. We consider an isotropic continuum damage mechanics model, written in
general form as
σ = f(ω)C :  (4.1)
where C is the elastic stiffness tensor, ω is a scalar damage variable ranging from 0 to 1,  is the
infinitesimal strain tensor, and f(ω) is a function which describes stiffness degradation of the bulk
material due to damage. The constitutive law is defined by the following parameters: E = Young’s
modulus, ν = Poisson’s ratio, ft0 = tensile strength, Gf0 = mode I fracture energy, and h0 =
crack band width, from which we model only the tensile strength and fracture energy as statistically
independent random variables, since these are governing the failure behavior. The subscript “0” is
4.1. MODEL FORMULATION 71
used to denote the true material properties, which are defined for an element of size equal to the crack
band width h0. The current analysis focuses on the case of finite elements of size larger than the
crack band width, and aims to formulate the probabilistic models of apparent tensile strength ft and
fracture energy Gf of such elements. It should be noted that the proposed sampling distributions
of apparent tensile strength and fracture energy of a given finite element need to be related to both
the inherent randomness of the corresponding material properties, as captured by the probability
distributions of ft0 and Gf0, as well as the failure mechanism experienced by the element, which will
be quantified from proposed numerical parameters resembling the localization parameter introduced
in Chapter 2.
Following the notation of Chapter 2, let us consider a finite element of in-plane dimensions
h1 and h2, where h1 denotes the element dimension perpendicular to the direction of incipient
crack propagation, and h2 the dimension orthogonal to h1 (Fig. 4.1). According the fictitious crack
model [74], cracks are considered to initiate and propagate inside an element as soon as its tensile
strength is reached. Denoting as nf the number of potential damage zone locations from which
fracture can initiate inside the element, we can expect that, due to the spatial variability of material
strength, the location of damage zone formation inside an element is intrinsically random, and in
this case nf = h1/hd, where hd is the size of the damage zone, to be determined from the failure
pattern of the element. However, in the presence of significant stress concentration, for instance due
to a propagating crack from the neighboring elements into the element of interest, there is only one
possible damage zone location from which fracture can initiate in the considered element, dictated
by mechanics rather than material randomness.
1
A finite element with a damage zone consisting 
of nb number of crack bands
ℎ𝑑
Figure 4.1: A finite element with a damage zone consisting of nb number of crack bands
Based on the foregoing discussion, we consider a weakest-link statistical representation of element
4.1. MODEL FORMULATION 72
failure, which reflects the randomness in the location of damage zone formation inside the element.
The probability distribution Pft of the tensile strength of a finite element can thus be expressed as:
Pft (σ) = 1−
[
1− P bft (σ)
]nf (4.2)
where nf = h1/hd and P bft is the probability distribution of strength of each damage zone, whose
size can vary between h0 and h1 depending of the failure mechanism of the element. For instance,
hd = h0 for localized damage through the formation of a single crack band, whereas hd = h1 for
purely diffusive damage inside the element. Based on Eq. (4.2), the entire element will fail under
controlled loads as soon as a single damage zone is formed.
The damage zone, by definition, experiences diffusive microcracking, therefore we statistically
represent it by a fiber-bundle model, which is suitable to describe this failure mechanism. In the
fiber-bundle model, the structure is modeled as a set of elements, called fibers, which are coupled in
parallel. Once a fiber fails, its load is redistributed among other surviving fibers, this way preventing
global failure. In this study, we consider that each fiber in the bundle model of the damage zone
represents a potential fracture process zone (FPZ) of width h0, and denote as nb = hd/h0 the number
of fibers in the bundle (Fig. 4.1). The crack band width is usually a few times the size of material
heterogeneities, and thus the reference element of size h1 = h0 is expected to exhibit a gradual
softening stress-strain behavior. To simplify the mathematical formulation of the distribution of
element strength, we assume that each fiber in the bundle model of the damage zone exhibits an
elastic-perfectly plastic stress-strain behavior, for which case the bundle strength is a weighted
sum of the strengths of the individual fibers. It has been shown [18] that the FPZ, which in our
case is represented by the elements of the bundle model, has a strength distribution that can be
approximated by a Gauss-Weibull grafted cumulative distribution function (cdf)
P1(σ) =

1− exp (−〈σ/s0〉m0) ≈ 〈σ/s0〉m0 (σ ≤ σgr0) (4.3a)
Pgr0 +
rf0
δG0
√
2pi
∫ σ
σgr0
e−(x−µG0)
2/2δ2G0dx (σ > σgr0) (4.3b)
where m0 and s0 are the Weibull modulus and scale parameter of the Weibull tail, respectively, µG0
and δG0 are the mean and standard deviation of the Gaussian core if considered extended to −∞; rf0
is a scaling parameter required to normalize the grafted cdf such that P1(∞) = 1, σgr0 = grafting
stress, and Pgr0 = (σgr0/s0)m0 = grafting probability. The continuity of the probability density
function (pdf) at the grafting point requires that [dP1/dσ]σ+gr0 = [dP1/dσ]σ−gr0 . Therefore, there are
in total four independent parameters that define the grafted distribution of the fiber strength, which
4.1. MODEL FORMULATION 73
in this study are taken to be f t0 = µG0, ω0 = δG0/µG0,m0, and Pgr0.
Because the crack band width h0 can be considered to be larger than the autocorrelation length
of the underlying random field of material strength, the fiber strengths can be treated as statistically
independent random variables. The plastic fiber-bundle model of the damage zone implies that the
cdf of bundle strength also follows a Gauss-Weibull grafted cdf, P bft = P1
(
f t, ω,m, Pgr
)
, which
can be defined by the four statistical parameters f t, ω,m, and Pgr representing the mean and
coefficient of variation of the Gaussian core of bundle strength, and the Weibull modulus and grafting
probability of its power-law tail respectively. According to the statistics of the sum of independent
Gaussian variables and the mathematical properties of the plastic fiber-bundle model, the statistical
parameters of the Gauss-Weibull cdf of bundle strength are given from the following equations
relating them to the corresponding parameters of the individual fibers [91,126]
f t = f t0 (4.4)
ω =
ω0
nb
[
dnbe − 1 + (nb + 1− dnbe)2
]1/2
(4.5)
m = m0dnbe (4.6)
Pgr = (Pgr0/dnbe)dnbe (4.7)
where it was considered that the bundle contains dnbe − 1 number of fibers with area A0 and a
single fiber with area (nb + 1− dnbe)A0, in order to ensure that the number of fibers is always an
integer, even if the quantity nb, which will be determined from numerical parameters quantifying
the damage pattern of an element (Eqs. (4.22), (4.23)), is non-integer.
Even though the weakest-link representation of element strength is suitable for failures featured
by damage localization, the proposed model of Eq. (4.2) is able to represent different failure mech-
anisms, because the quantities nf and nb governing the scaling of the strength distribution are not
predetermined constants. Instead, they are calculated during the simulation and can evolve accord-
ing to the actual failure mechanism of the element. At this point, it can be illustrative to provide
the limiting values of nf and nb during the three stages that are generally present in any quasibrittle
failure process, which are: 1) damage initiation, characterized by a diffused damage pattern through
the formation of multiple microcracks, 2) damage localization, during which the previously formed
microcracks coalesce into a single macrocrack, and 3) damage propagation, characterized by propa-
gation of the developed macrocrack leading to ultimate failure. Based on the definitions of nf and
nb, for a given element we should thus have: 1) nf = 1 and nb = h1/h0 during damage initiation,
4.1. MODEL FORMULATION 74
since the zone of diffusive cracking (damage zone) occupies the entire element, 2) nf = h1/h0 and
nb = 1 during damage localization, since a single crack has been formed in the element out of the
h1/h0 possible cracks with random locations of formation, 3) nf = 1 and nb = 1 during propagation
of a band of localized strains from the neighboring elements into the considered element, since stress
concentration due to damage localization in the neighbors determines the location of crack band
formation in the considered element.
As will be shown later, the mathematical expressions based on which nf and nb will be calculated
satisfy these limits, and upon their substitution to Eq. (4.2), it can be seen that, in cases of diffusive
damage, the entire element is statistically represented by the fiber-bundle model, in the case of
damage localization, the proposed model transitions to the weakest-link representation of element
failure, and in the case of damage propagation, the distribution of element strength coincides with the
Gauss-Weibull strength cdf of the reference element of size equal to the crack band width (Eq. (4.3)).
For all cases lying between purely diffused and purely localized damage, the proposed model is
transitional among the fiber-bundle, the weakest-link, and the constant distribution models. It is
also important to note that the dependence of the probability distributions of constitutive properties
on the damage pattern, as introduced in the present model, requires an assumption for the cdf of
tensile strength ft at the initial state, if undamaged. Since the randomness of material properties
can precipitate damage localization, we assume a weakest-link representation of structural failure
prior to damage initiation, whereas after the onset of damage the proposed model evolves during
the simulation according to the actual damage pattern. Eqs. (4.2), (4.3) and (4.4)-(4.7) provide
the formulation of the probability distribution of tensile strength of a given finite element, which is
governed by its failure mechanism.
The probabilistic model of apparent fracture energy Gf , which is the second constitutive param-
eter modeled as random variable in the present study, will be derived from energetic considerations
during propagation of nb number of crack bands of width h0 inside the damage zone formed in the
corresponding element. The energy required (per unit thickness) for simultaneous advancement of
nb cracks by the same distance δa under quasi-static conditions and assuming fracture as the only
dissipative process, is equal to nbGfδa = nbGfh2δω, where Gf represents the apparent fracture
energy of the element under mode I conditions and δa = h2δω, based on our interpretation of the
damage variable in Chapter 2. Crack growth is made possible through the release of strain energy
from the entire element during the increment of the crack length(s) δa, which can be written as
h1h2δU , where U represents a free energy function that governs tensile failure. The energy balance
4.1. MODEL FORMULATION 75
equation in the limit of δa→ 0, and thus δω → 0, yields
∂U
∂ω
+
Gf
h1
nb = 0 (4.8)
The fracture energy is defined as the energy required to create a unit of new crack surface area.
For this general case considered, with 1 ≤ nb ≤ h1/h0 number of crack bands that could be developed
inside the element, the total energy dissipated in failure is
Ed =
nb∑
i=1
AsiGfi (4.9)
where Gfi and Asi denote the fracture energy and surface area of crack i, respectively, and it is
considered that the fracture energies of all crack bands are statistically independent and identically
distributed random variables. In case when nb is not an integer, we consider that inside the element
there have been formed dnbe − 1 number of cracks with surface area As0 and a single crack with
area (nb + 1− dnbe)As0. The apparent fracture energy of the element can be obtained by dividing
Eq. (4.9) with the total surface area of all developed cracks, thus
Gf =
∑dnbe−1
i=1 As0Gfi + (nb + 1− dnbe)As0Gf0
nbAs0
(4.10)
Eq. (4.10) represents a weighted sum of statistically independent random variables, thus, following
the Central Limit Theorem (CLT), the distribution of apparent fracture energy Gf will approach
the Gaussian distribution when sufficiently large number of cracks have been developed, thus when
failure is characterized by a purely diffusive damage pattern.
In the present study, we consider that the fracture energies of elements of size equal to the crack
band width h0, which represent the actual material properties, are Gaussian random variables with
the same mean Gf0 and coefficient of variation ω′0, i.e. Gfi ∼ N
(
Gf0, ω
′
0
)
. Therefore, the apparent
fracture energy of the element will also follow the Gaussian distribution even if nb is not large enough
for the CLT to hold, with the following parameters:
Gf ∼ N
(
Gf0,
ω′0
nb
[
dnbe − 1 + (nb + 1− dnbe)2
]1/2)
(4.11)
Based on Eq. (4.11), in cases when failure of an element is featured by damage localization due
to a single crack that has either been formed inside the element or is about to propagate from its
neighbors, nb = 1, and so the distribution of the apparent fracture energy of the element coincides
4.1. MODEL FORMULATION 76
with that of the fracture energy of a single crack band, Gf ∼ N
(
Gf0, ω
′
0
)
. In cases of multiple
crack formation, the mechanics of energy dissipation lead to a size dependence of the coefficient of
variation of apparent fracture energy which decreases with the element size, and a size independent
mean.
Having completed the formulation of our probabilistic model of apparent fracture properties,
we can now combine it with our mechanism-based energy regularization scheme (Chapter 2). We
thus focus on Eq. (4.8), which is a general energy balance equation applicable to different failure
mechanisms depending on the value of nb. Eq. (4.8) is equivalent to Eq. (2.5) which is also reproduced
below, if the energy regularization function g is related to quantity nb as in Eq. (4.13).
∂U
∂ω
+
Gf
h0
g = 0 (4.12)
g =
(
h0
h1
)
nb (4.13)
Similarly to Chapter 2, function g transitions from 1 in cases of a purely diffused damage pattern
to h0/h1 in the case of damage localization. It has also been shown (Chapter 2 and [68]), that
satisfaction of the general energy balance equation (Eqs. (4.8) or (4.12)) by the damage mechanics
constitutive law is essential in order to preserve correct energy dissipation during failure regardless
of the spatial discretization, and thus ensure the mesh-objectivity of finite element simulations of
strain-softening continua in the deterministic setting. In the case of failure with damage localization,
i.e. nb = 1 and Eq. (4.8) implies a size-dependent input constitutive behavior, which is essential
for mitigating the otherwise mesh-dependent simulated structural response [17, 68]. Additionally,
the size dependence of the input sampling distributions of ft and Gf arising due to the failure
mechanism of a given element and described by the current probabilistic model (Eqs. (4.2), (4.4)-
(4.7) and (4.11)), will be shown to be effective in mitigating the mesh dependence of the output first
and second statistical moments of structural strength in the present case of stochastic analysis. From
Eqs. (4.12) and (4.13), it can be seen that the energy regularization scheme and the probabilistic
models of fracture properties are interconnected, since the term nb is present in both the general
energy balance equation as well as the distribution functions of ft and Gf .
Finally, we reproduce below the expressions for the damage function f(ω) and the damage
evolution law ω(¯) as taken from Chapter 2, which are required along with Eq. (4.1) to fully define
our constitutive model:
f(ω) =
1− ω
1− aω (4.14)
4.1. MODEL FORMULATION 77
ω =
1
a
− ¯
a
√
h0E(1− a)
2Gfg
(4.15)
a = 1− 2EGfg
h0f2t
(4.16)
where a is a parameter which depends on the material properties as well as the energy regularization
function g, not to be confused with the crack length a for which the same notation was used, and ¯ =
Mazar’s equivalent strain (Eq. 2.8). Damage is allowed to initiate only after the equivalent strain
¯ exceeds the threshold value 0 = ft/E related to the element tensile strength, and the damage
variable is irreversible, i.e. ω˙ ≥ 0. Eqs. (4.1), (4.14), (4.15) as well as Eqs. (4.2) and (4.11) complete
the description of mechanical and probabilistic response of a material element of finite size, which
are intimately tied to its failure mechanism.
4.1.2 Description of the Local Damage Pattern
The probabilistic and constitutive models described in Section 4.1.1 require the calculation of quan-
tities nf and nb, representing the number of potential damage zones inside an element and the
number of crack bands inside the damage zone respectively. It is clear that these quantities should
be related to the damage pattern of the corresponding finite elements. In this section, we introduce
two numerical parameters which capture the evolving damage pattern around each element, and
continue to describe the details of their implementation within the finite element procedure. Quan-
tities nf and nb can then be determined by the introduced numerical parameters based on proposed
expressions.
The first numerical parameter is called the localization parameter δω1 and has the same physical
meaning as parameter δω introduced in the energy regualization scheme of Chapter 2. δω1 aims
to describe the level of damage localization inside a given finite element. The value of δω1 ranges
from 0, representing the case when the considered element is either under a diffused damage state or
when damage localization cannot occur, to δ1m, representing a fully localized damage pattern. The
second parameter, δω2, is called the propagation parameter. Its minimum value is 0, corresponding
to the case when it is unlikely that a crack from the neighboring finite elements will propagate into
the element of interest, and its maximum value is δ2m, corresponding to the case when damage
propagation will most certainly occur.
In the finite element implementation of our constitutive model, we follow the procedure outlined
below to determine parameters δω1 and δω2, which is suitable for the case when linear shape functions
are used for the interpolation of the displacement field. For a given finite element i, we initially
4.1. MODEL FORMULATION 78
compute its average strain tensor from the corresponding strain values at all the element Gauss
points. We denote as ~ne the unit vector with the direction of maximum principal average strain
at the onset of damage. It is considered that cracks will form in the direction perpendicular to
~ne, which is kept constant after damage initiation, as assumed in the fixed crack model. Length
h1 can now be calculated by projecting the element onto the direction of ~ne and has therefore the
same value for all Gauss points of a given element, while length h2 will correspond to the element
dimension orthogonal to h1. At every Gauss point of element i, we construct the unit vectors ~v
g
ij
connecting the current Gauss point with the centroid of neighboring element j (Fig. 4.2), and among
all possible connection vectors, we identify those that are most closely aligned with the vector ~ne,
and denote them as ~vgm1 and ~v
g
m2. The localization parameter of every Gauss point δω1, can then
be calculated as follows
δω1 =
〈
1− ωˆi1 + ωˆi2
2ωˆi
〉s
(4.17)
where 〈x〉 = max(x, 0), ωˆi = characteristic damage level of element i, and ωˆi1, ωˆi2 = characteristic
damage levels of the two adjacent elements that are along the directions of ~vgm1 and ~v
g
m2. From the
definition of δω1 in Eq. (4.17), it is observed that δω1 is equal to 0 in the case of uniform distribution
of the damage ωˆ along the direction perpendicular to the crack path, as well as in the case when the
characteristic damage level of element i is lower than the average of the characteristic damage levels
of its adjacent elements most closely aligned with ~ne, since localization of damage in element i is
unlikely to occur under the conditions of the second case. δω1 attains its maximum value δ1m when
a fully localized damage pattern is observed, i.e. only the element of interest i experiences damage,
whereas its neighboring elements along the directions of ~vgm1 and ~v
g
m2 are undamaged.
The second parameter introduced in the present model is the propagation parameter δω2, which
is calculated for every Gauss point of the considered element i as follows:
δω2 = maxj
[
δω1
p
j
(
1− ‖~ne · ~vgij‖
)q] (4.18)
where δω1j is the average value of the localization parameter at a neighboring element j, and the
exponents p and q are set equal to 1 and 2 respectively. δω2 is calculated as the maximum of the
inside the brackets quantity in Eq. (4.18) among all neighboring elements j. Based on Eq. (4.18),
the propagation parameter is taken equal to the average localization parameter of the neighboring
element that is most closely aligned with the direction of crack propagation inside element i, since it
is reasonable to consider that, after damage localizes in this specific neighbor, the developed crack
will propagate into the element of interest. Here one could note that, in contrast to δω1, which
4.1. MODEL FORMULATION 79
aims to describe the level of damage localization inside a given element and is affected by both the
damage state of this element as well as that of its neighbors in the direction perpendicular to the
crack path, parameter δω2 does not take into account the damage pattern of the element of interest,
but only that of its neighbors along the crack direction.
2
Schematic of finite element mesh for calculation of 
the localization and propagation parameters
𝑛𝑒𝑖
𝑗
Ԧ𝜐𝑚1
𝑔Ԧ𝜐𝑚2
𝑔
Ԧ𝜐𝑖𝑗
𝑔
Figure 4.2: Schematic of finite element mesh for calculation of the localization and propagation
parameters
The characteristic damage levels ωˆk employed for the calculation of the localization and propa-
gation parameters, where k denotes the index of the corresponding element, are defined as follows:
ωˆk = αwω¯k + (1− αw)
(
ω¯k−1 + ω¯k+1
2
)
(4.19)
where ω¯k denotes the average damage level of element k, as calculated by using the damage values at
all Gauss points of the element, ω¯k−1, ω¯k+1 are the average damage levels of its neighboring elements
along the direction perpendicular to the direction of crack propagation inside element k, and αw is
a selected weighting coefficient in the interval (0, 1]. In the case of αw < 1 considered in the present
study, the characteristic damage level of a given element k is a nonlocal quantity, since it additionally
contains information of the average damage levels of elements k − 1 and k + 1. Introduction of a
nonlocal damage quantity for the determination of parameters δω1 and δω2 is essential in the case
of random properties considered here, since, due to the spatial variability of fracture properties, a
diffusive damage pattern of the structure may not be characterized by a homogeneous distribution
of damage, even for loadings that create a uniform strain field. Thus, calculation of the localization
parameter δω1 from the nonlocal damage levels of the corresponding elements ensures that it will
not reach its maximum value indicating full localization, when there exist weak spatial variations of
4.1. MODEL FORMULATION 80
the damage variable inside a predominantly diffused damage pattern. The maximum value of the
localization parameter can now be calculated from Eqs. (4.17) and (4.19) applied to a fully localized
damage pattern, i.e. only the element of interest i is damaged, whereas its neighbors 1 and 2 are
undamaged. This yields:
δ1m =
〈
3αw − 1
2αw
〉s
(4.20)
Based on Eq. (4.18), we can also calculate the maximum value of the propagation parameter δω2
δ2m = δ
p
1m if p > 0 (4.21)
where p = 1 in the present study, so the maximum values of the localization and propagation
parameters coincide. We should also note that in the initial, undamaged state, δ0ω1 = δ1m and
δ0ω2 = 0. As will be shown (Eqs. (4.22)-(4.24)), these initial conditions ensure that the tensile
strength ft is sampled from a weakest-link model of h1/h0 number of links prior to damage initiation,
consistent with our previously mentioned assumption in Section 4.1.1.
Having proposed expressions for the localization and propagation parameters δω1 and δω2 in
terms of the characteristic damage levels, we can quantitatively capture the failure mechanism
of any element, and use this information inside the energy regularization scheme as well as for
the determination of quantities nb and nf . We introduce the following equation for the energy
regularization function g, which scales the post-peak part of the input constitutive law according to
the failure mechanism of the element, so as to ensure correct amount of energy dissipation:
g =
h0
h1
+
(
1− h0
h1
){
1−
[ 〈δω1 − δt1〉
δ1m − δt1
]κ}{
1−
[ 〈δω2 − δt2〉
δ2m − δt2
]λ}
(4.22)
where δt1, δt2, κ, λ are parameters to be determined from numerical simulations. For δω2 = 0,
δt1 = 0 and αw = 1, Eq. (4.22) is equivalent to g = 1 +
(
h0
h1
− 1
)
δκω1, a functional form which,
when combined with a similar definition of the localization parameter δω1 (Eq. (4.17)), has been
shown (Chapter 2) to be efficient in mitigating mesh dependence of deterministic finite element
simulations of structures made of strain-softening materials under general failure mechanisms [68].
Solving Eq. (4.13) with respect to nb, we get:
nb =
h1
h0
g (4.23)
4.2. NUMERICAL EXAMPLES 81
and below, we introduce the relation for the number of potential damage zones nf inside an element
nf =
h1
nbh0
+
(
1− h1
nbh0
)[ 〈δω2 − δt2〉
δ2m − δt2
]λ
(4.24)
Eqs. (4.22)-(4.24) recover the limiting behaviors of function g and of quantities nb and nf ,
including damage initiation (δω1 = 0, δω2 = 0), localization (δω1 = δ1m, δω2 = 0), and propagation
(δω2 = δ2m). Section 4.1 has provided a complete description of the probabilistic and constitutive
model, which are interconnected through the energy regularization scheme. The remaining sections
of the current chapter present the application of the model in stochastic FE simulations of quasibrittle
structures experiencing different failure mechanisms, and evaluate its efficiency in mitigating mesh
dependence of the output strength statistics.
4.2 Numerical Examples
4.2.1 Description of Analysis
The performance of the proposed model in mitigating mesh dependence of the output first and
second statistical moments of structural strength is evaluated by its application in stochastic FE
simulations of quasibrittle specimens of different geometries under different loading configurations.
The simulations are performed in the open-source FEM software OOFEM [115], where the present
constitutive and probabilistic models were implemented. A two-dimensional state of plane stress
is assumed. The specimens are subjected to the three-point and four-point bending tests, under
displacement control conditions. For the three-point bending test, two different structure geometries
are considered; an unnotched beam of in-plane dimensions 7.5 mm (length) × 2.5 mm (depth)
(Fig. 4.3), and a notched beam of the same dimensions, with relative notch depth 0.5 and notch
width w = h0/2, where h0 = crack band width (Fig. 4.6). Only the unnotched beam is subjected to
the four-point bending test (Fig. 4.8). The beams are made of a dense alumina ceramic, with elastic
properties treated as deterministic quantities with values: E = 370GPa and ν = 0.24. The crack
band width is set constant and equal to h0 = 37.5 µm, which is a few times the average grain size
of the material.
The tensile strength ft and fracture energyGf of different finite elements in the discretization vary
spatially and are modeled as statistically independent random variables. Their values are calculated
by inverting the input sampling distributions at probability levels which are randomly generated at
the beginning of every simulation, and are the same among all Gauss points of a given element. The
4.2. NUMERICAL EXAMPLES 82
probability levels remain constant during each simulation, but the cdfs of fracture properties evolve
with the damage pattern of an element, following the current probabilistic model. For comparison
purposes, we consider another two probabilistic models of the fracture properties: the weakest-link
model and the assumption of constant input distribution functions. More specifically, the three
different probabilistic models employed in the simulations are:
I) The proposed model, where the tensile strength and fracture energy are sampled from the
distributions of Eqs. (4.2) and (4.11) respectively, which are related to the actual failure mech-
anism of the structure through quantities nf and nb, and can lead to a transition from a
fiber-bundle model statistical representation of failure to a weakest-link model as damage lo-
calization emerges from diffusive damage.
II) The weakest-link model, where the distribution of tensile strength takes the form Pft =
1 − (1− P1)h1/h0 with P1 sampled from Eq. (4.3), and the fracture energy follows a Gaussian
distribution with size independent statistical parameters corresponding to the properties of an
element of size equal to the crack band width, Gf ∼ N
(
Gf0, ω
′
0
)
. The weakest link model is
based on the assumption of structural failure with damage localization.
III) The assumption of constant distributions, where the material properties of tensile strength and
fracture energy are sampled from a Gauss-Weibull grafted cdf, Pft = P1
(
f t0,m0, ω0, Pgr0
)
,
and a purely Gaussian distribution Gf ∼ N
(
Gf0, ω
′
0
)
, respectively, with statistical parameters
which do not depend on the element size.
The input statistical parameters of the Gauss-Weibull cdf of material tensile strength are: f t0 =
200MPa, ω0 = 0.15, m0 = 30, Pgr0 = 5 × 10−3, and of the Gaussian distribution of fracture en-
ergy are Gf0 = 50 J/m2, ω′0 = 0.15. These values characterize the inherent randomness of material
properties of an element of size equal to the crack band width. It should be noted that the en-
ergy regularization scheme described in Section 4.1.2, Eq. (4.22), is applied to all aforementioned
probabilistic models, so that the corresponding deterministic problems are regularized and thus any
observed mesh dependence of the output statistics of structural strength when using approaches
I)-III) can be mainly attributed to incorrect scaling of the corresponding probability distributions
of constitutive properties.
All beams are discretized by a uniform mesh consisting of quadrilateral bilinear elements, using
full Gauss integration. To investigate the mesh dependence, three different discretizations of the
beams are considered: the reference case, with dimensions of finite elements hx = h0 and hy = h0,
4.2. NUMERICAL EXAMPLES 83
where subscripts x, y, denote the directions of the beam span and depth respectively, and two other
coarser discretizations, with hx = 2h0 and hx = 4h0 whereas hy is kept constant. The mesh lines
are aligned with the direction of overall crack band growth to avoid mesh-induced directional bias,
so the element size along the direction of the maximum principal strain h1 is approximately equal
to hx.
For each loading case and structure geometry, we calculated the mean and standard deviation of
the peak load of the structure, as captured from the output load-midpoint displacement diagram.
A sufficient number of random simulations was performed, ranging between 300-500, such that
convergence of both statistics of peak load is achieved. In all three different probabilistic models
used, the parameters of the energy regularization function δt1, δt2, κ, λ (Eq. (4.22)) were calibrated
such that the mesh dependence of the output mean peak load is alleviated, so the models can then
be compared in terms of their ability to mitigate mesh dependence of the second statistical moment
of structural strength. The remaining parameters of the model required in Eqs. (4.17)-(4.21) are
constant and set equal to: s = 0.5, p = 1, q = 2, aw = 0.5. Parameters δt1, δt2, κ, and λ, once
calibrated for a given material and a selected probabilistic model (I-III), retain the same values
among the different structure geometries and loading configurations.
4.2.2 Unnotched beams under three-point bending
In this section we investigate the performance of all three different cases of sampling of the ran-
dom variables of element tensile strength ft and fracture energy Gf (approaches I-III), in terms of
their ability to mitigate mesh dependence of the statistics of peak load of three-point bend beams
(Fig. 4.3). Fig. 4.4 presents the output mean and standard deviation of peak load, as calculated
from 500 random simulations, for all probabilistic models and element sizes hx = h0, 2h0, 4h0. The
parameters δt1, δt2, κ, and λ have been calibrated separately for each probabilistic model such that
the mean size effect is mitigated, as shown in Fig. 4.4 (left). The optimum values of these parameters
were found to be: δt1 = 0, δt2 = 0, κ = 1, λ = 1 for the proposed model and the case of constant
distributions, whereas δt1 = 0.7, δt2 = 0.7, κ = 1, λ = 1 for the weakest-link model. Given that all
models perform excellent in terms of the mean behavior, we can focus on their ability to mitigate the
size dependence of the second statistical moment. In Fig. 4.4 (right), we observe that the optimum
result is achieved by the proposed model, which gives a maximum error for the standard deviation
of the coarser mesh compared to the reference mesh of approximately 11%, followed by the case of
constant distributions with maximum error around 16%. The mesh dependence issue of the output
standard deviation when using the weakest-link model is exacerbated.
4.2. NUMERICAL EXAMPLES 84
3
Unnotched beam under three-point bending considered in 
the simulation (all dimensions are in mm)
Figure 4.3: Geometry of unnotched beams under three-point bending considered in the simulation
(all dimensions are in mm)
TPB unnotched (Gf=50J/m2)
1
Mean peak load of the unnotched beam under three-point 
bending with probabilistic models I)-III) (Gf0=50J/m2)
Standard deviation of peak load of the unnotched beam under three-
point bending with probabilistic models I)-III) (Gf0=50J/m2)
Figure 4.4: Mean (left) and standard deviation (right) of peak load of the unnotched beam under
three-point bending with probabilistic models I)-III) (Gf0 = 50J/m2)
The resulting size dependence of the second statistical moment of structural strength can be ex-
plained by the failure mechanism of this unnotched beam under three-point bending. Even though
the randomness of material properties can lead to different number of damaged elements as well as
different locations of crack initiation between random simulations, when the peak load is reached,
the most prevalent failure mechanism of this beam is failure where multiple cracking contributes
significantly to the peak load capacity. This means that the maximum load is governed by the con-
tributions of all softening elements inside the damage zone, and is not solely dictated by the strength
of the weakest element. Such averaging effect can lead to a decrease in the standard deviation of
structural strength with increasing number of elements inside the damage zone. Considering the
overall damage zone size being approximately constant, since it is mainly governed by the struc-
4.2. NUMERICAL EXAMPLES 85
ture geometry and loading configuration, a coarser discretization implies fewer elements inside the
damage zone, and thus less averaging effect. The above considerations explain why the present
simulation results show an increasing trend of the standard deviation of peak load with the mesh
size for all three probabilistic models.
The present model has the ability to “feel” the averaging effect because its sampling distributions
of fracture properties depend on the localization parameter δω1, which quantifies the actual damage
pattern around each element. In cases of distributed cracking, δω1 << δ1m so the model will
switch from the initial assumption of weakest-link failure to the fiber-bundle model, if there is not
significant propagation of localized damage, i.e. δω2 << δ2m. According to the fiber-bundle model,
the input CoV of tensile strength and fracture energy is scaled such that larger elements have less
variable properties (Eqs. (4.5)-(4.6),(4.11)), and this scaling is shown to significantly mitigate the
mesh dependence of the output standard deviation of structural strength, so the increasing trend
of standard deviation with the mesh size is insignificant when the proposed model is used. In the
case of constant distributions, though, the size dependence of the standard deviation becomes more
pronounced, since this approach does not introduce any dependence of the statistics of fracture
properties on the element size.
The worst result is obtained when the weakest-link model is used, which leads to a drastic increase
of the standard deviation with the mesh size. The weakest-link assumption is inconsistent with the
actual failure mechanism of the structure, which justifies the poor performance of this model. More
specifically, the weakest-link model assigns lower mean strength to larger elements, in contrast to
the fiber-bundle model as well as the case of constant distributions, where the input mean properties
are size-independent. An initial set of trial simulations with the weakest-link model and parameters
of the energy regularization function g set equal to δt1 = 0, δt2 = 0, κ = 1, λ = 1, which were
the optimum values of the parameters calibrated for approaches I) and III), resulted in a quite
pronounced reduction of the mean peak load of the structure with the mesh size. To improve this
result, thresholds δt1 and δt2 were both increased to the value of 0.7 < δ1m = δ2m = 0.707, which
means that the energy regularization function g was set equal to 1 for a large part of the damage
process, thus no scaling of the post-peak part of the softening law was applied. In the presence of
some level of damage localization, absence of energy regularization implied by using g = 1, leads
to larger elements having higher fracture energy so the mean peak load of these elements increased
and counteracted the reduction of mean peak load with the mesh size imposed by the weakest-link
sampling, resulting in a mesh-independent mean behavior. However, the output standard deviation
exhibits a strong increasing trend with the mesh size, even though the weakest-link model introduces
4.2. NUMERICAL EXAMPLES 86
a scaling of the input CoV such that larger elements have less variable tensile strength, unless the
large-size limit is reached, for which the CoV depends only on Weibull’s modulus m0. The reason for
this significant increase in the standard deviation with the mesh size probably lies in the fact that
use of g = 1 over a large range of the damage process to improve the mean size effect, may be an
assumption not reflecting the actual failure mechanism of the beam, which is best quantified when
using the calculated values of the localization and propagation parameters δω1, δω2 inside the energy
regularization scheme, that is when δt1 << δ1m and δt2 << δ2m. We could also note that, given that
the present model, which is transitional among the weakest-link, the fiber-bundle, and the constant
distributions models, has the best performance followed by the case of constant distributions, the
most prevalent failure mechanism of the considered beam must lie between distributed cracking
which, at some locations, localizes and propagates into neighboring elements.
In order to further demonstrate the ability of the proposed model to switch between the different
limiting cases of failure with distributed cracking, weakest-link type of failure, and failure due to
propagation of localized damage, we now consider a beam of exactly the same geometry (Fig. 4.3),
but with reduced input mean material fracture energy from Gf0 = 50 J/m2 to 15 J/m2. Since in the
present study we aim to cover many different failure mechanisms, we now consider a nominally less
ductile beam to present a case for which the weakest-link assumption becomes more appropriate.
Due to the different material, the relevant parameters δt1, δt2, κ and λ need to be recalibrated for
each probabilistic model. Their optimum values in this case are: δt1 = 0.6, δt2 = 0.6, κ = 1, λ = 1
for the proposed model, δt1 = 0, δt2 = 0, κ = 1, λ = 1 for the case of constant distributions, and
δt1 = 0.7, δt2 = 0.7, κ = 1, λ = 1 for the weakest-link model.
Fig. 4.5 presents the mesh size dependence of the output mean and standard deviation of peak
load for all three probabilistic failure models, where convergence was achieved after 300 random
simulations. Even though the weakest-link model leads to a non-negligible error for the mean
peak load of around 6%, its overall performance in mitigating mesh dependence of both mean
and standard deviation of structural strength is significantly improved compared to its application
to the previously considered more ductile material. We can expect that further reduction of the
fracture energy or increase of the beam size such that failure becomes perfectly brittle will lead to
the weakest-link model having the best performance, and our model converging to the weakest-link
model. In this study, we considered the first alternative to change the failure type, since increase
of the beam size would increase the computational cost. However, the input fracture energy we
could use was limited by the snap-back condition, h1 < 2EGf/f2t , which is required by our energy
regularization scheme whenever a fully localized damage pattern occurs. Due to the randomness,
4.2. NUMERICAL EXAMPLES 87
the actual sampled element fracture energies can be less than the mean value, which is the initial
input in the model, so there is no a priori way to determine the minimum value of input mean
fracture energy to avoid snap-back. For the selected value of Gf0 = 15 J/m2, which led to very few
random simulations exhibiting snap-back that were discarded, it seems that failure of the considered
beam is still not perfectly brittle, though the contribution of distributed cracking to the peak load
is significantly reduced. The latter can also be verified by the fact that the mesh dependence of
the output standard deviation when using the constant distributions is much more pronounced,
indicating less averaging effect. Lastly, we can observe that the present model has now shifted such
that is closer to the weakest-link model than to the case of constant distributions, and has similar
performance to the weakest-link model in terms of the standard deviation, but significantly improved
performance in terms of the mean size effect.
TPB unnotched (Gf=15J/m2)
Mean peak load of the unnotched beam under three-point 
bending with probabilistic models I)-III) (Gf0=15J/m2)
Standard deviation of peak load of the unnotched beam under 
three-point bending with probabilistic models I)-III) (Gf0=15J/m2) 2
Figure 4.5: Mean (left) and standard deviation (right) of peak load of the unnotched beam under
three-point bending with probabilistic models I)-III) (Gf0 = 15J/m2)
4.2.3 Notched beam under three-point bending
In the preceding analysis of unnotched beams under three-point bending, we focused on two inter-
mediate cases of failure lying between the limits of purely diffused and purely localized damage, for
which the performance of the proposed model in mitigating mesh dependence of the output strength
statistics has been satisfactory. Most importantly, it was shown that the model is able to switch
from the fiber-bundle to the weakest-link model as the failure process transitions from predomi-
nantly diffused to more localized, which is a consequence of the imposed dependence of the sampling
distribution functions of fracture properties on the numerical parameters δω1 and δω2 capturing the
4.2. NUMERICAL EXAMPLES 88
actual damage pattern of an element.
In the current and following sections, we further evaluate the consistency of probabilistic and
mechanical behavior in the present model, by verifying it applies the correct scaling of the proba-
bility distributions as well as of the post-peak slope of the constitutive law in the limiting cases of
propagation of localized damage and of a purely diffused damage pattern. We consider a notched
beam of the same size as the unnotched beam of the previous section (7.5 mm × 2.5 mm), having
relative notch depth 0.5 (Fig. 4.6). All elastic properties as well as input statistical parameters of
fracture properties retain the values provided in Section 4.2.1. To investigate the mesh dependence
of the output statistics of peak load of the beam, we again consider the same set of element sizes
hx = h0, 2h0, 4h0 whereas hy = h0 = constant, and a regular mesh consisting of rectangular elements
aligned with the direction of overall crack growth. It should be noted that, in order to ensure a large
notch width over notch depth ratio, there is a single element just above the notch, of trapezoidal
shape, with lower side representing the notch width of dimension w = h0/2, and upper side of di-
mension equal to the current dimension hx of the subsequent rectangular finite elements. Use of a
trapezoidal element with bases w and hx at the notch tip allows us not only to specify a sufficiently
small and constant notch width w such that the resulting stress concentration is significant and the
energy regularization scheme is efficient, but also to always have a single column of elements above
the notch, regardless of the selected value of hx.
4
Notched beam under three-point bending considered in the 
simulation (all dimensions are in mm)
1
.2
5
Figure 4.6: Notched beam under three-point bending considered in the simulation (all dimensions
are in mm)
Stress concentration induced by the presence of the notch dictates the location of damage initi-
ation to be confined inside the single element above the notch tip, regardless of the randomness of
material properties. Consequently, the mean peak load is a deterministic quantity whereas material
randomness affects only the variability of structural strength, as is known in this case of Type II
4.2. NUMERICAL EXAMPLES 89
size effect [15, 59, 103]. Strain localization is already pronounced when the peak load is reached, so
the energy regularization scheme will switch to the crack band model [17], and will be able to fully
recover mesh objectivity of the mean peak load, as shown in the results of the proposed model in
Fig. 4.7 (left), obtained with the previously calibrated optimum set of parameters for this material
(δt1 = 0, δt2 = 0, κ = 1, λ = 1). Following strain localization at the element above the notch,
damage propagation into a single band of elements will occur, and this will lead to the propagation
parameter δω2 quickly reach its maximum value, so the present model will transition to the case of
constant distributions. This can be verified by Fig. 4.7, where it is clear that the output statistics
obtained by application of the present model and of the assumption of constant distributions almost
coincide. Additionally, these two approaches are the most efficient in mitigating mesh dependence
of the output statistics of peak load. Given that, at the peak load, there is always a single column
of damaged elements just above the notch, consisting of approximately the same number of failed
elements since the dimension hy is kept constant between the different discretizations, the output
variability of peak load should directly reflect the input randomness, and should thus be independent
of the element size if the fracture properties are sampled from size-independent distributions.
Based on the above discussion, we expect that application of the weakest-link model in the
considered case of the notched beam will have a poor performance in providing size-independent
statistics of peak load capacity, since this model always introduces a dependence of the statisti-
cal parameters on the element size. More specifically, the weakest-link model assigns lower mean
strength and coefficient of variation to larger elements, apart from the perfectly brittle case, where
the imposed size-effect is still present in the mean but not in the coefficient of variation. Given that
the simulated notched beam has uncracked ligament of depth comparable to Irwin’s characteristic
length, as calculated using the mean values of the relevant material properties, its failure is not
perfectly brittle, so the peak load is controlled by both the tensile strength and fracture energy, and
thus its statistics will depend on the element size hx as a consequence of the size dependence of the
sampled from the weakest-link model tensile strengths. As shown in Fig. 4.7, both the mean as well
as the standard deviation of peak load eventually decrease for the largest element size used. The
mesh dependence of the second statistical moment is more pronounced than the corresponding of the
mean peak load, since the parameters of the energy regularization scheme, which remain the same as
in the unnotched beam of the same material (δt1 = 0.7, δt2 = 0.7, κ = 1, λ = 1), were calibrated such
that the mean size effect is mitigated. However, the induced delay in the transition of the energy
regularization scheme to the crack band model in order to increase the mean peak load of larger
elements and thus alleviate the mean size effect, is totally inconsistent with the failure mechanism
4.2. NUMERICAL EXAMPLES 90
of the notched beam, which requires the scaling of the post-peak branch of the constitutive law such
that the fracture energy is preserved to be applied as soon as damage initiates. Only towards large
beam sizes, we expect to obtain mesh-objective strength statistics even when using the weakest-link
model of tensile strength, since in that case the peak load of the structure is solely governed by the
fracture energy, which will be preserved by the energy regularization scheme and will be sampled
from the reference size-independent Gaussian distribution (Gf ∼ N
(
Gf0, ω
′
0
)
), due to the purely
localized damage pattern (nb = 1 when g = h0/h1) of a notched beam.
TPB notched (Gf=50J/m2)
Mean peak load of the notched beam under three-point 
bending with probabilistic models I)-III) 3
Standard deviation of peak load of the notched beam under three-point 
bending with probabilistic models I)-III) 
Figure 4.7: Mean (left) and standard deviation (right) of peak load of the notched beam under
three-point bending with probabilistic models I)-III)
4.2.4 Unnotched beam under four-point bending
The last limiting case of failure for which we aim to evaluate the performance of the current prob-
abilistic model is purely diffusive damage. We thus consider the four-point bending test, which can
induce more diffusive cracking due to the development in the specimen of a region of constant bend-
ing moment between the concentrated loads. We simulate an unnotched beam of inner span 7.5 mm
(same as the total length of the three-point bent beams), overall span of 11.25 mm, and depth 2.5
mm. A regular mesh of quadrilateral elements aligned with the expected direction of crack growth is
used, similarly to the previous simulations. To investigate mesh dependence, we consider the three
different element sizes hx = h0, 2h0, 4h0 whereas hy = h0 = constant. The input elastic and fracture
properties of the beam remain the same and have been provided in Section 4.2.1. The experimental
set-up is mimicked in the finite element model by applying increasing levels of vertical displacement
at the midpoint of the top part of the beam, and imposing the sum of the corresponding displace-
4.2. NUMERICAL EXAMPLES 91
ments of the end points of the inner span to equal the applied midpoint displacement (Fig. 4.8). To
avoid damage initiation at the loaded regions, the entire upper layer of the beam as well as its cross
sections above the supports have been discretized with elements of extremely high tensile strength
such that they always remain elastic. Such elastic layers have been used in all simulations presented
in the current study.
11
Unnotched beam under four-point bending considered in 
the simulation (all dimensions are in mm)
Figure 4.8: Unnotched beam under four-point bending considered in the simulation (all dimensions
are in mm)
The probabilistic structural response of the beam is captured by calculating the mean and stan-
dard deviation of peak load, where 400 random simulations were required to achieve convergence of
the second statistical moment. The sampling of the fracture properties was performed by approaches
I)-III), with optimum values of parameters as provided in the results of the three-point bending test
for the same material, i.e. δt1 = 0, δt2 = 0, κ = 1, λ = 1 for the proposed model and the case of
constant distributions, and δt1 = 0.7, δt2 = 0.7, κ = 1, λ = 1 for the weakest-link model.
Fig. 4.9 provides the simulated first and second statistical moments of peak load of the beam
with respect to the element size hx, for all three probabilistic failure models. As in all previous
tests, we observe that all models are able to effectively mitigate the mean size effect, since their
parameters have been calibrated such that the mean behavior is optimized. Regarding the mesh
dependence of the second statistical moment, it can be seen that the present model is able to provide
an almost mesh-insensitive result. Due to the purely diffused damage pattern of the four-point bend
beam, the simplification of a plastic fiber-bundle model of the damage zone introduced in the present
4.3. CONCLUSIONS 92
model is indeed realistic, which explains why it exhibits the highest performance for this test. The
assumption of constant distributions is also acceptable, since the error in the standard deviation
is approximately 7%. On account of the larger size of the damage zone in the four-point bending
test, the averaging effect is significantly more pronounced than in the three-point bending, and is
able to decrease the standard deviation for coarser discretizations as well, to the extent that even
the simplest case of size-independent distributions can lead to output strength statistics that are
only mildly dependent on the element size. In fact, due to the significant damage zone size, the
performance of all three probabilistic models in the four-point bending test is improved compared to
the three-point bending test of the unnotched beam of the same material. However, the weakest-link
model still leads to an unacceptably high error of the standard deviation, which is shown in Fig. 4.9
to strongly increase with the element size.
FPB unnotched (Gf=50J/m2)
4
Mean peak load of the unnotched beam under four-point 
bending with probabilistic models I)-III) 
Standard deviation of peak load of the unnotched beam under 
four-point bending with probabilistic models I)-III) 
Figure 4.9: Mean (left) and standard deviation (right) of peak load of the unnotched beam under
four-point bending with probabilistic models I)-III)
4.3 Conclusions
The present study introduces a general probabilistic model of apparent fracture properties, which is
informed by both the inherent material randomness as well as the failure mechanism of a structure.
The mathematical formulation of the model (Section 4.1.1) ensures it can recover the well-established
weakest-link and fiber-bundle models in the limiting cases of failure with damage localization and
failure with distributed cracking respectively. Only the spatial variability of tensile strength and
fracture energy is considered, as parameters on which the failure behavior of quasibrittle structures
strongly depends. The formulated probability distributions of element strength and fracture en-
4.3. CONCLUSIONS 93
ergy are used as the sampling distribution functions of constitutive parameters in stochastic FE
simulations of quasibrittle specimens of different geometries, under different loading configurations.
The model does not require any a priori knowledge of the structural failure mechanism, since the
distribution functions are made dependent on numerical parameters that can describe the evolving
damage pattern in the structure, and lead to a transition of the probabilistic model towards the
failure model most suitable with the actual damage mechanism.
The considered numerical examples of notched and unnotched beams under the three-point bend-
ing test, as well as of an unnotched beam subjected to four-point bending, are selected such that a
wide range of structural failure mechanisms is triggered and therefore the generality of the model
can be tested. It is shown that in all aforementioned tests resulting in damage patterns transition-
ing from purely localized to purely diffused through intermediate complex states, the model is able
to switch to the failure model that is consistent with the damage mechanism experienced by the
structure. An important feature of the probabilistic model is that is coupled with our recently devel-
oped mechanism-based energy regularization scheme [68], which was shown to be an efficient way of
mitigating mesh dependence of the mechanical structural response in deterministic FE simulations
of strain-softening continua, experiencing different damage mechanisms. In the present context of
stochastic FE simulations, the size effects on the statistics of fracture properties that can arise dur-
ing the evolving between loading steps, damage pattern of an element, are explicitly modeled by
accordingly adjusting the input sampling distributions, and such adjustment plays an important role
in mitigating mesh dependence of the output probabilistic structural response. Indeed, comparison
of the present probabilistic model with other, predetermined distribution functions of constitutive
parameters, shows that its performance in mitigating mesh dependence of the output statistics of
structural strength is significantly higher than the performance of the prescribed probability dis-
tributions, which can be efficient only if the failure mechanism based on which they were derived
happens to be the actual failure mechanism of the structure. Other benefits of the present model
include its computational efficiency which makes it an attractive option for stochastic FE simulation
of large-scale structures, where most probably the selected element sizes will be larger than the au-
tocorrelation length of material properties, as well as simplicity, since the same spatial discretization
used for the mechanical problem is also used for modeling of the randomness of fracture properties.
Chapter 5
Conclusions
This research introduces a general framework to address the issue of spurious mesh sensitivity of finite
element (FE) simulations of strain-softening continua, for both deterministic and stochastic analyses.
In the deterministic setting, mesh sensitivity is successfully mitigated by a mechanism-based energy
regularization of the constitutive law, which is governed by the actual damage mechanism in a
structure and ensures consistent energy dissipation in failure regardless of the mesh discretization. In
the stochastic setting, the structural failure mechanism is additionally considered in the formulation
of the sampling probability distributions of apparent fracture properties. The randomness of element
tensile strength and fracture energy is described by a probabilistic model, that can be informed by
the inherent variability of material properties, but more importantly, is adjusted based on the actual
failure pattern of an element and can transition between the weakest-link and fiber-bundle models
for the cases of purely localized and purely diffused damage respectively. It should be emphasized
that the behavior of each finite element is considered to represent the behavior of a material element
of finite size, so the input constitutive and probabilistic models for the finite elements represent
both material as well as structural behavior of the corresponding material elements. The latter is
manifested through the failure mechanism, which is quantified by numerical parameters describing
the nonlocal damage pattern around each element. These parameters are used in both the energy
regularization scheme as well as the probability distributions of fracture properties. The introduction
of these descriptors allows us to capture the evolving during the simulation damage pattern, which
obviates the need for any a priori knowledge of the structural failure mechanism. Therefore, the
proposed method can be applied to quasibrittle structures under general loading scenarios.
The following conclusions can be drawn based on this study:
94
95
1. Energy regularization of the constitutive law constitutes an important tool for mitigating the
spurious mesh sensitivity of FE simulations of strain-softening continua, owing to its generality
since it is connected with the energetic approach to fracture. The developed mechanism-
based energy regularization scheme overcomes the limitations of the most prevalent energy
regularization method, i.e. the crack band model. The present model is able to mitigate mesh
dependence under various failure mechanisms, and recover the crack band model in cases
of a purely localized damage pattern. Additionally, the crack band model is tied with the
assumption that the fracture energy is the most important parameter governing the failure
behavior and which should thus be preserved. For complicated constitutive laws, selection of
the relevant parameters to preserve can be challenging. Such decision is not required though
with the present model, since it adjusts the point-wise tangential stiffness tensor based on the
current state of damage pattern.
2. Under dynamic loading conditions, energy regularization can be made possible by controlling
the damage growth rate based on the actual failure mechanism. More specifically, if the
applied strain rate and structure geometry are such that a purely diffused damage pattern is
observed, then material point behavior and structural behavior coincide so the damage should
propagate at the same rate regardless of the element size. However, in cases of a localized
damage pattern (e.g. notched specimen under high strain-rate loading), the resulting more
brittle failure behavior of larger elements implies that damage should propagate faster inside
them, so the damage growth rate should depend on the element size.
3. Stochastic discrete element model simulations of unnotched quasibrittle specimens subjected
to a wide range of tensile strain rates show the transition of the damage pattern from fairly
localized to purely diffused with the rate. This result has important implications for the scaling
of the distribution of structural strength. It is shown that, for low to intermediate strain rates,
the mean and standard deviation of specimen strength decrease with an increasing specimen
size. At high strain rates, the size effect on the mean diminishes whereas the decreasing trend of
standard deviation becomes stronger and it exhibits a power-law scaling behavior, as would be
predicted by a fiber-bundle model statistical representation of failure. The results also indicate
that, for a given specimen size, an increasing strain rate leads to a reduction in the standard
deviation of nominal strength due to the more diffusive cracking pattern. The simulated rate
and size effects on the strength statistics have been captured by a rate-dependent weakest-
link model, which transitions from a weakest-link to a fiber-bundle model representation of
96
structural failure as high strain rates lead to a more diffusive damage pattern. This part of the
research has provided insight into the formulation of a general probabilistic model of apparent
fracture properties, that is governed by the actual failure mechanism of a structure.
4. The present research proposes a new method to mitigate the mesh dependence of probabilistic
structural response in stochastic FE simulations of softening continua. Instead of making an a
priori assumption for the structural failure mechanism and sampling the input tensile strength
from a prescribed probability distribution, a common practice in the literature, it is suggested
to adjust the sampling distributions of constitutive properties according to the evolving damage
pattern in the structure. The proposed model is formulated in a general form so that it
can recover the well-established probabilistic models suitable for some limiting failure cases,
such as the weakest-link and fiber-bundle models. The model also considers the effect of the
failure mechanism on the distribution of fracture energy, and is coupled with the mechanism-
based energy regularization scheme, introduced in the deterministic analysis. Stochastic FE
simulations of quasibrittle structures of different geometries subjected to different loading
configurations, show that the present model can efficiently mitigate the mesh dependence of
the output first and second statistical moments of structural strength.
Bibliography
[1] V. Agrawal and K. Dayal. Dependence of equilibrium Griffith surface energy on crack speed
in phase-field models for fracture coupled to elastodynamics. Int. J. Fract., 207(2):243–249,
2017.
[2] R. Alessi, M. Ambati, T. Gerasimov, S. Vidoli, and L. De Lorenzis. Comparison of Phase-
Field Models of Fracture Coupled with Plasticity. In Advances in Computational Plasticity,
volume 46 of Computational Methods in Applied Sciences, pages 1–21. Springer International
Publishing, Cham, 2017.
[3] M. Ambati, T. Gerasimov, and L. De Lorenzis. Phase-field modeling of ductile fracture.
Comput. Mech., 55(5):1017–1040, 2015.
[4] S. Esna Ashari, G. Buscarnera, and G. Cusatis. A lattice discrete particle model for pressure-
dependent inelasticity in granular rocks. Int. J. Rock Mech. Min. Sci., 91:49–58, 2017.
[5] K. Bagi. Microstructural Stress Tensor of Granular Assemblies With Volume Forces. ASME
J. Appl. Mech., 66(4):934–936, 1999.
[6] G. I. Barenblatt. The Mathematical Theory of Equilibrium Cracks in Brittle Fracture. volume 7
of Adv. Appl. Mech., pages 55 – 129. Elsevier, 1962.
[7] K.-J. Bathe. Finite Element Procedures. Prentice Hall, Englewood Cliffs, N.J., 1996.
[8] Z. P. Bažant. Scaling of Structural Strength. Elsevier, London, 2005.
[9] Z. P. Bažant and F. C. Caner. Comminution of solids caused by kinetic energy of high shear
strain, with implications for impact, shock and shale fracturing. Proc. Natl. Acad. Sci. U.S.A.,
110(48):19291–19294, 2013.
[10] Z. P. Bažant and F. C. Caner. Impact comminution of solids due to local kinetic energy of
high shear strain rate: I. Continuum theory and turbulence analogy. J. Mech. Phys. Solids,
64:223–235, 2014.
[11] Z. P. Bažant and L. Cedolin. Blunt crack band propagation in finite element analysis. J. Eng.
Mech. Div. ASCE, 105(2):297–315, 1979.
[12] Z. P. Bažant and L. Cedolin. Stability of Structures: Elastic, Inelastic, Fracture and Damage
Theories. Oxford University Press, New York, 1991.
[13] Z. P. Bažant and R. Gettu. Rate Effects and Load Relaxation in Static Fracture of Concrete.
ACI Mater. J., 89(5):456–468, 1992.
97
BIBLIOGRAPHY 98
[14] Z. P. Bažant and J.-L. Le. Nano-mechanics based modeling of lifetime distribution of quasib-
rittle structures. Eng. Fail. Anal., 16:2521–2529, 2009.
[15] Z. P. Bažant and J.-L. Le. Probabilistic Mechanics of Quasibrittle Structures: Strength, Life-
time, and Size Effect. Cambridge University Press, Cambridge, U.K., 2017.
[16] Z. P. Bažant, J.-L. Le, and C. G. Hoover. Nonlocal boundary layer model: overcoming bound-
ary condition problems in strength statistics and fracture analysis of quasibrittle materials. In
The 7th International Conference on Fracture Mechanics of Concrete and Concrete Structures
(FraMCoS-7), pages 135–143, Jeju, Korea, 2010.
[17] Z. P. Bažant and B. H. Oh. Crack band theory for fracture of concrete. Mater. Struct.,
16:155–177, 1983.
[18] Z. P. Bažant and S. D. Pang. Activation energy based extreme value statistics and size effect
in brittle and quasibrittle fracture. J. Mech. Phys. Solids, 55:91–131, 2007.
[19] Z. P. Bažant and J. Planas. Fracture and Size Effect in Concrete and Other Quasibrittle
Materials. CRC Press, London, 1998.
[20] Z. P. Bažant. Softening Instability: Part I–Localization Into a Planar Band. In J. Appl. Mech.,
volume 55, pages 517–522, New York, NY, 1988. ASME.
[21] Z. P. Bažant. Softening Instability: Part II–Localization Into Ellipsoidal Regions. In J. Appl.
Mech., volume 55, pages 523–529, New York, NY, 1988. ASME.
[22] Z. P. Bažant. Stable States and Paths of Structures with Plasticity or Damage. J. Eng. Mech.,
114(12):2013–2034, 1988.
[23] Z. P. Bažant. Design of quasibrittle materials and structures to optimize strength and scaling
at probability tail: an apercu. Proc. R. Soc. A, 475(2224):20180617, 2019.
[24] Z. P. Bažant, T. B. Belytschko, and T.-P. Chang. Continuum Theory for Strain-Softening.
ASCE J. Eng. Mech., 110(12):1666–1692, 1984.
[25] Z. P. Bažant and M. Jirásek. Nonlocal Integral Formulations of Plasticity and Damage: Survey
of Progress. ASCE J. Eng. Mech., 128(11):1119–1149, 2002.
[26] Z. P. Bažant and S. D. Pang. Mechanics-based statistics of failure risk of quasibrittle structures
and size effect on safety factors. Proc. Natl. Acad. Sci. U.S.A., 103(25):9434–9439, 2006.
[27] Z. P. Bažant and G. Pijaudier-Cabot. Measurement of Characteristic Length of Nonlocal
Continuum. ASCE J. Eng. Mech., 115(4):755–767, 1989.
[28] Z. P. Bažant, M. R. Tabbara, M. T. Kazemi, and G. Pijaudier-Cabot. Random Particle Model
for Fracture of Aggregate or Fiber Composites. ASCE J. Eng. Mech., 116(8):1686–1705, 1990.
[29] Z. P. Bažant and Q. Yu. Size-Effect Testing of Cohesive Fracture Parameters and Nonunique-
ness of Work-of-Fracture Method. ASCE J. Eng. Mech., 137(8):580–588, 2011.
[30] T. Belytschko, H. Chen, J. Xu, and G. Zi. Dynamic crack propagation based on loss of
hyperbolicity and a new discontinuous enrichment. Int. J. Numer. Methods Eng., 58(12):1873–
1905, 2003.
[31] E. Berthier, V. Démery, and L. Ponson. Damage spreading in quasi-brittle disordered solids:
I. Localization and failure. J. Mech. Phys. Solids, 102:101–124, 2017.
BIBLIOGRAPHY 99
[32] H. S. Bhat, A. J. Rosakis, and C. G. Sammis. A micromechanics based constitutive model for
brittle failure at high strain rates. ASME J. Appl. Mech., 79:031016, 2012.
[33] J. Bleyer and J.-F. Molinari. Microbranching instability in phase-field modelling of dynamic
brittle fracture. Appl. Phys. Lett., 110(15):151903, 2017.
[34] J. Bleyer, C. Roux-Langlois, and J.-F. Molinari. Dynamic crack propagation with a variational
phase-field model: limiting speed, crack branching and velocity-toughening mechanisms. Int.
J. Fract., 204(1):79–100, 2017.
[35] W. R. Blumenthal. High Strain Rate Compression Testing of Ceramics and Ceramic Com-
posites. In J. J. Swab, editor, Advances in Ceramic Armor Ceramic Engineering and Science
Proceedings, volume 26, pages 89–96, 2005.
[36] J. E. Bolander, G. S. Hong, and K. Yoshitake. Structural Concrete Analysis Using Rigid-
Body-Spring Networks. J. Comput. Aided Civil Infrastruct. Eng., 15(2):120–133, 2000.
[37] J. E. Bolander and S. Saito. Fracture analyses using spring networks with random geometry.
Eng. Fract. Mech., 61(5):569–591, 1998.
[38] G. Borino, B. Failla, and F. Parrinello. A symmetric nonlocal damage theory. Int. J. Solids
Struct., 40(13):3621 – 3645, 2003.
[39] B. Bourdin, G. A. Francfort, and J-J. Marigo. Numerical experiments in revisited brittle
fracture. J. Mech. Phys. Solids, 48(4):797–826, 2000.
[40] B. Bourdin, G. A. Francfort, and J-J. Marigo. The Variational Approach to Fracture. Springer
Dordrecht, Netherlands, 2008.
[41] B. Bourdin, G. A. Francfort, and J.-J. Marigo. The Variational Approach to Fracture. J.
Elasticity, 91(1-3):5–148, 2008.
[42] A. Braides. Approximation of Free-Discontinuity Problems. Lecture Notes in Mathematics.
Springer-Verlag, Berlin and New York, 1998.
[43] N. Chandra, H. Li, C. Shet, and H. Ghonem. Some issues in the application of cohesive zone
models for metal–ceramic interfaces. Int. J. Solids Struct., 39(10):2827 – 2855, 2002.
[44] C.-H. Chen, E. Bouchbinder, and A. Karma. Instability in dynamic fracture and the failure
of the classical theory of cracks. Nat. Phys., 13(12):1186–1190, 2017.
[45] P. A. Cundall. A Computer Model for Simulating Progressive Large Scale Movements in Blocky
Rock Systems. In Proc. Int. Symp. Rock Fracture, volume 1, pages II–8, Nancy, France, 1971.
[46] P. A Cundall and O. D. L. Strack. A discrete numerical model for granular assemblies. Géotech-
nique, 29(1):47–65, 1979.
[47] G. Cusatis, Z. P. Bažant, and L. Cedolin. Confinement-Shear Lattice Model for Concrete
Damage in Tension and Compression: I. Theory. ASCE J. Eng. Mech., 129(12):1439–1448,
2003.
[48] G. Cusatis, Z. P. Bažant, and L. Cedolin. Confinement-Shear Lattice Model for Concrete
Damage in Tension and Compression: II. Computation and Validation. ASCE J. Eng. Mech.,
129(12):1449–1458, 2003.
BIBLIOGRAPHY 100
[49] G. Cusatis and L. Cedolin. Two-scale study of concrete fracturing behavior. Eng. Fract. Mech.,
74(1–2):3–17, 2007. Fracture of Concrete Materials and Structures.
[50] G. Cusatis, A. Mencarelli, D. Pelessone, and J. Baylot. Lattice Discrete Particle Model
(LDPM) for failure behavior of concrete. II: Calibration and validation. Cem. Concr. Compos.,
33(9):891–905, 2011.
[51] G. Cusatis, D. Pelessone, and A. Mencarelli. Lattice Discrete Particle Model (LDPM) for
failure behavior of concrete. I: Theory. Cem. Concr. Compos., 22:881–890, 2011.
[52] G. Cusatis and E. A. Schauffert. Cohesive crack analysis of size effect. Eng. Fract. Mech.,
76(14):2163 – 2173, 2009.
[53] H. E. Daniels. The statistical theory of the strength of bundles of threads. I. Proc. R. Soc.
Lond. A, 183(995):405–435, 1945.
[54] N. P. Daphalapurkar, K.T. Ramesh, L. Graham-Brady, and J.-F. Molinari. Predicting vari-
ability in the dynamic failure strength of brittle materials considering pre-existing flaws. J.
Mech. Phys. Solids, 59:297–319, 2011.
[55] R. Desmorat, M. Chambart, F. Gatuingt, and D. Guilbaud. Delay-active damage versus non-
local enhancement for anisotropic damage dynamics computations with alternated loading.
Eng. Fract. Mech., 77(12):2294 – 2315, 2010.
[56] F. P. Duda, A. Ciarbonetti, P. J. Sánchez, and A. E. Huespe. A phase-field/gradient damage
model for brittle fracture in elastic–plastic solids. Int. J. Plast., 65:269–296, 2015.
[57] D. S. Dugdale. Yielding of steel sheets containing slits. J. Mech. Phys. Solids, 8(2):100 – 104,
1960.
[58] J. Eliáš. Adaptive technique for discrete models of fracture. Int. J. Solids Struct., 100-101:376–
387, 2016.
[59] J. Eliáš, M. Vořechovský, J. Skoček, and Z. P. Bažant. Stochastic discrete meso-scale simula-
tions of concrete fracture: Comparison to experimental data. Eng. Fract. Mech., 135(1):1–16,
2015.
[60] B. Erzar and P. Forquin. An Experimental Method to Determine the Tensile Strength of
Concrete at High Rates of Strain. Exp. Mech., 50(7):941–955, 2010.
[61] S. Feng and P. N. Sen. Percolation on elastic networks: new exponent and threshold. Phys.
Rev. Lett., 52(3):216–219, 1984.
[62] P. Forquin and B. Lukić. On the Processing of Spalling Experiments. Part I: Identification of
the Dynamic Tensile Strength of Concrete. J. Dyn. Behav. Mater., 4:34–55, 2018.
[63] G. A. Francfort and J.-J. Marigo. Revisiting brittle fracture as an energy minimization prob-
lem. J. Mech. Phys. Solids, 46(8):1319–1342, 1998.
[64] F. Gatuingt, L. Snozzi, and J. F. Molinari. Numerical determination of the tensile response
and the dissipated fracture energy of concrete: role of the mesostructure and influence of the
loading rate. Int. J. Numer. Anal. Methods Geomech., 37:3112–3130, 2013.
[65] R. G. Ghanem and P. D. Spanos. Stochastic Finite Elements : A Spectral Approach. Springer-
Verlag, New York, 1991.
BIBLIOGRAPHY 101
[66] A. Giacomini. Size Effects on Quasi-Static Growth of Cracks. SIAM J. Math. Anal.,
36(6):1887–1928, 2005.
[67] C. Giry, F. Dufour, and J. Mazars. Stress-based nonlocal damage model. Int. J. Solids Struct.,
48(25-26):3431–3443, 2011.
[68] A. Gorgogianni, J. Eliáš, and J.-L. Le. Mechanism-Based Energy Regularization in Computa-
tional Modeling of Quasibrittle Fracture. ASME J. Appl. Mech., 87(9), 2020.
[69] P. Grassl and Z. P. Bažant. Random Lattice-Particle Simulation of Statistical Size Effect in
Quasi-Brittle Structures Failing at Crack Initiation. ASCE J. Eng. Mech., 135(2):85–92, 2009.
[70] A. A. Griffith and G. I. Taylor. VI. The Phenomena of Rupture and Flow in Solids. Trans.
R. Soc. Lond. A, 221(582-593):163–198, 1921.
[71] M. Grigoriu. Simulation of Stationary Non-Gaussian Translation Processes. J. Eng. Mech.,
124(2):121–126, 1998.
[72] D.L. Grote, S.W. Park, and M. Zhou. Dynamic behavior of concrete at high strain rates and
pressures: I. experimental characterization. Int. J. Impact Eng., 25:869–886, 2001.
[73] H. J. Herrmann, A. Hansen, and S. Roux. Fracture of disordered, elastic lattices in two
dimensions. Phys. Rev. B, 39(1):637–648, 1989.
[74] A. Hillerborg, M. Modéer, and P.-E. Petersson. Analysis of crack formation and crack growth
in concrete by means of fracture mechanics and finite elements. Cem. Concr. Res., 6(6):773 –
781, 1976.
[75] A. Hrennikoff. Solution of problems of elasticity by the framework method. J. Appl. Mech.,
8(4):169–175, 1941.
[76] Y. K. Hwang, J. Bolander, and Y. M. Lim. Simulation of concrete tensile failure under high
loading rates using three-dimensional irregular lattice models. Mech. Mater., 101:136–146,
2016.
[77] K. Iwashita and M. Oda. Mechanics of Granular Materials: An Introduction. CRC press,
Boca Raton, FL, 1999.
[78] M. Jirásek. Nonlocal Theories in Continuum Mechanics. Acta Polytech., 44(5-6):16–34, 2004.
[79] M. Jirásek and M. Bauer. Numerical aspects of the crack band approach. Comp. Struct.,
110-111:60–78, 2012.
[80] M. Jirásek and Z. P. Bažant. Macroscopic fracture characteristics of random particle systems.
Int. J. Fract., 69(3):201–228, 1994.
[81] M. Jirásek and Z. P. Bažant. Particle Model for Quasibrittle Fracture and Application to Sea
Ice. ASCE J. Eng. Mech., 121(9):1016–1025, 1995.
[82] M. Jirásek. Mathematical analysis of strain localization. REGC, 11(7-8):977–991, 2007.
[83] L. Kachanov. Introduction to Continuum Damage Mechanics. Springer, Netherlands, 1986.
[84] A. Kara, A. Tasdemirci, and M. Guden. Modeling quasi-static and high strain rate deformation
and failure behavior of a (±45) symmetric E-glass/polyester composite under compressive
loading. Mater. Des., 49:566–574, 2013.
BIBLIOGRAPHY 102
[85] T. Kawai. New discrete models and their application to seismic response analysis of struc-
tures. Nucl. Eng. Des., 48(1):207–229, 1978. Special Issue Structural Mechanics in Reactor
Technology - Smirt-4.
[86] O. Keita, C. Dascalu, and B. François. A two-scale model for dynamic damage evolution. J.
Mech. Phys. Solids, 64:170–183, 2014.
[87] J. Kimberley, K. T. Ramesh, and N. P. Daphalapurkar. A scaling law for the dynamic strength
of brittle solids. Acta Mater., 61:3509–3521, 2013.
[88] R. H. Kraft, J. F. Molinari, K. T. Ramesh, and D. H. Warner. Computational micromechanics
of dynamic compressive loading of a brittle polycrystalline material using a distribution of
grain boundary properties. J. Mech. Phys. Solids, 56:2618–2641, 2008.
[89] F. Kun, S. Zapperi, and H.J. Herrmann. Damage in fiber bundle models. Eur. Phys. J. B,
17(2):269–279, 2000.
[90] J.-L. Le, Z. P. Bažant, and M. Z. Bazant. Unified nano-mechanics based probabilistic theory
of quasibrittle and brittle structures: I. strength, crack growth, lifetime and scaling. J. Mech.
Phys. Solids, 59:1291–1321, 2011.
[91] J.-L. Le, J. Eliáš, A. Gorgogianni, J. Vievering, and J. Kveton. Rate-Dependent Scaling of
Dynamic Tensile Strength of Quasibrittle Structures. ASME J. Appl. Mech., 85(2):021003,
2018.
[92] J.-L. Le and J. Eliáš. A probabilistic crack band model for quasibrittle fracture. ASME J.
Appl. Mech., 83(5):051005, 2016.
[93] J.-L. Le, Z. Xu, and J. Eliáš. Internal Length Scale of Weakest-Link Statistical Model for
Quasi-Brittle Fracture. ASCE J. Eng. Mech., 144(4), 2017.
[94] M. Lemaitre and J. L. Chaboche. Aspect Phenomenologique de la Rupture par Endommage-
ment. J. Mec. Appl., 2:317–365, 1978.
[95] C.-C. Li and A. Der Kiureghian. Optimal Discretization of Random Fields. ASCE J. Eng.
Mech., 119(6):1136–1154, 1993.
[96] Z. Li and J. Lambros. Determination of the dynamic response of brittle composites by the use
of the split Hopkinson pressure bar. Compos. Sci. Technol., 59:1097–1107, 1999.
[97] E. Lorentz and V. Godard. Gradient damage models: towards full-scale computations. Com-
put. Methods Appl. Mech. Eng., 200:1927–1944, 2011.
[98] Y. B. Lu and Q. M. Li. About the dynamic uniaxial tensile strength of concrete-like materials.
Int. J. Impact Eng., 38(4):171–180, 2011.
[99] W. Luo and Z. P. Bažant. Fishnet statistics for probabilistic strength and scaling of nacreous
imbricated lamellar materials. J. Mech. Phys. Solids, 109:264–287, 2017.
[100] W. Luo and Z. P. Bažant. Fishnet model with order statistics for tail probability of failure
of nacreous biomimetic materials with softening interlaminar links. J. Mech. Phys. Solids,
121:281–295, 2018.
[101] W. Luo and Z. P. Bažant. Fishnet Statistical Size Effect on Strength of Materials With
Nacreous Microstructure. ASME J. Appl. Mech., 86(8), 2019.
BIBLIOGRAPHY 103
[102] W. Luo and Z. P. Bažant. General Fishnet Statistics of Strength: Nacreous, Biomimetic,
Concrete, Octet-Truss, and Other Architected or Quasibrittle Materials. ASME J. Appl.
Mech., 87(3):1–14, 2019.
[103] W. Luo, J.-L. Le, M. Rasoolinejad, and Z. P. Bažant. Coefficient of Variation of Shear Strength
of RC Beams and Size Effect. J. Eng. Mech., 147(2), 2021.
[104] S. May, J. Vignollet, and R. de Borst. A numerical assessment of phase-field models for brittle
and cohesive fracture: Γ-Convergence and stress oscillations. Eur. J. Mech. A Solids, 52:72–84,
2015.
[105] J. Mazars. Application de la mécanique de l’endommagement au comportement non linéaire
et à la rupture du béton de structure. PhD thesis, University of Paris VI, Paris, France, 1984.
[106] C. Miehe, M. Hofacker, L.-M. Schänzel, and F. Aldakheel. Phase field modeling of fracture
in multi-physics problems. Part II. Coupled brittle-to-ductile failure criteria and crack prop-
agation in thermo-elastic–plastic solids. Comput. Methods Appl. Mech. Eng., 294:486–522,
2015.
[107] C. Miehe, F. Welschinger, and M. Hofacker. Thermodynamically consistent phase-field models
of fracture: Variational principles and multi-field FE implementations. Int. J. Numer. Methods
Eng., 83(10):1273–1311, 2010.
[108] S. Murakami. Notion of Continuum Damage Mechanics and its Application to Anisotropic
Creep Damage Theory. J. Eng. Mater. Technol., 105(2):99–105, 04 1983.
[109] N. M. Newmark. A Method of Computation for Structural Dynamics. University of Illinois,
Urbana, 1959.
[110] T. T. Nguyen, J. Yvonnet, M. Bornert, C. Chateau, K. Sab, R. Romani, and R. Le Roy.
On the choice of parameters in the phase field method for simulating crack initiation with
experimental validation. Int. J. Fract., 197(2):213–226, 2016.
[111] T.T. Nguyen, J. Yvonnet, Q.-Z Zhu, M. Bornert, and C. Chateau. A phase-field method for
computational modeling of interfacial damage interacting with crack propagation in realistic
microstructures obtained by microtomography. Comput. Methods Appl. Mech. Eng., 312:567–
595, 2016.
[112] J. Ožbolt, J. Bošnjak, and E. Sola. Dynamic fracture of concrete compact tension specimen:
Experimental and numerical study. Int. J. Solids Struct., 50(25-26):4270–4278, 2013.
[113] J. Ožbolt, A. Sharma, and H.-W. Reinhardt. Dynamic fracture of concrete – compact tension
specimen. Int. J. Solids Struct., 48(10):1534–1543, 2011.
[114] Z. Pan, R. Ma, D. Wang, and A. Chen. A review of lattice type model in fracture mechanics:
theory, applications, and perspectives. Eng. Fract. Mech., 190:382–409, 2018.
[115] B. Patzák. OOFEM - An Object-Oriented Simulation Tool for Advanced Modeling of Materials
and Structures. Acta Polytech., 52(6):59–66, 2012.
[116] L. F. Pereira, J. Weerheijm, and L. J. Sluys. A new effective rate dependent damage model
for dynamic tensile failure of concrete. Eng. Fract. Mech., 176:281–299, 2017.
[117] S. L. Phoenix. Stochastic strength and fatigue of fiber bundles. Int. J. Fract., 14(3):327–344,
1978.
BIBLIOGRAPHY 104
[118] G. Pijaudier-Cabot and J. Ludovic. Continuum damage modelling and some computational
issues. Revue Française de Génie Civil, 6(6):991 – 1017, 2002.
[119] Y. N. Rabotnov. Creep rupture. In M. Hetényi and W. G. Vincenti, editors, Applied Mechan-
ics, pages 342–349, Berlin, Heidelberg, 1968. International Union of Theoretical and Applied
Mechanics, Springer.
[120] R. Rezakhani and G. Cusatis. Asymptotic expansion homogenization of discrete fine-scale
models with rotational degrees of freedom for the simulation of quasi-brittle materials. J.
Mech. Phys. Solids, 88:320–345, 2016.
[121] A. J. Rosakis. Explosion at the Parthenon: Can We Pick Up the Pieces? Technical Report
GALCIT SM Report 99-3, California Institute of Technology, 1999.
[122] L. Rothenburg and R. J. Bathurst. Analytical study of induced anisotropy in idealized granular
materials. Géotechnique, 39(4):601–614, 1989.
[123] S. Roux and E. Guyon. Mechanical percolation : a small beam lattice study. J. Physique
Lett., 46(21):999–1004, 1985.
[124] J. W. Rudnicki and J. R. Rice. Conditions for the localization of deformation in pressure-
sensitive dilatant materials. J. Mech. Phys. Solids, 23:371–394, 1975.
[125] M. Sahimi and J. D. Goddard. Elastic percolation models for cohesive mechanical failure in
heterogeneous systems. Phys. Rev. B, 33(11):7848–7851, 1986.
[126] M. Salviato and Z. P. Bažant. The asymptotic stochastic strength of bundles of elements
exhibiting general stress-strain laws. Probabilistic Eng. Mech., 36:1–7, 2014.
[127] E. Schlangen and J.G.M van Mier. Experimental and numerical analysis of micromechanisms
of fracture of cement-based composites. Cem. Concr. Compos., 14(2):105–118, 1992.
[128] H. Schuler, C. Mayrhofer, and K. Thoma. Spall experiments for the measurement of the
tensile strength and fracture energy of concrete at high strain rates. Int. J. Impact Eng.,
32:1635–1650, 2006.
[129] M. Shinozuka and G. Deodatis. Simulation of Multi-Dimensional Gaussian Stochastic Fields
by Spectral Representation. Appl. Mech. Rev., 49(1):29–53, 1996.
[130] F. Tonon. Explicit Exact Formulas for the 3-D Tetrahedron Inertia Tensor in Terms of its
Vertex Coordinates. J. Math. Stat., 1(1):8–11, 2005.
[131] C. V. Verhoosel and R. de Borst. A phase-field model for cohesive fracture. Int. J. Numer.
Methods Eng., 96(1):43–62, 2013.
[132] J. Vignollet, S. May, R. de Borst, and C. V. Verhoosel. Phase-field models for brittle and
cohesive fracture. Mecc., 49(11):2587–2601, 2014.
[133] J. Weerheijm, I. Vegt, and K. van Breugel. Research developments and experimental data
on dynamic concrete behaviour. In Conference on Advances in Construction Materials, pages
765–773, Stuttgart, Germany, 2007.
[134] W. Weibull. A Statistical Distribution Function of Wide Applicability. ASME J. Appl. Mech.,
153(18):293–297, 1951.
BIBLIOGRAPHY 105
[135] J.-Y. Wu. A unified phase-field theory for the mechanics of damage and quasi-brittle failure.
J. Mech. Phys. Solids, 103:72–99, 2017.
[136] J.-Y. Wu and V. P. Nguyen. A length scale insensitive phase-field damage model for brittle
fracture. J. Mech. Phys. Solids, 119:20–42, 2018.
[137] J. Y. Wu, V. P. Nguyen, C. T. Nguyen, D. Sutula, S. Sinaie, and S. Bordas. Phase-field
modeling of fracture. Adv. Appl. Mech., 2019.
[138] F. Zhou and J. F. Molinari. Stochastic fracture of ceramics under dynamic tensile loading.
Int. J. Solids Struct., 41:6573–6596, 2004.
[139] J. M. Ziman and D. Thouless. Models of disorder: The theoretical physics of homogeneously
disordered systems. Phys. Today, 33(6):57–58, 1980.
[140] J. L. Zinszner, B. Erzar, P. Forquin, and E. Buzaud. Dynamic fragmentation of an alumina
ceramic subjected to shockless spalling : An experimental and numerical study. J. Mech.
Phys. Solids, 85:112–127, 2015.
[141] A. Zubelewicz and Z. P. Bažant. Interface Element Modeling of Fracture in Aggregate Com-
posites. ASCE J. Eng. Mech., 113(11):1619–1630, 1987.