In the notebook [6.3](../tour6_energy/6_3_localized_energy_dissipation.ipynb#top) we have analyzed the correspondence between the damage function inducing softening and the propagation of discontinuity macroscopically observable as debonding in a pullout test. We mentioned the correspondence between the locally and globally energy evaluated energy dissipation indicating the energy dissipated at a material point can be used to estimate the total amount of dissipated energy by relating it to a localized crack of a known known surface area. In case of an interface area between two material components, like CFRP sheet and concrete, the area is well defined.
Let us now apply the same concept to a crack propagation through a concrete specimen. Even though additional aspects, like two-dimensional stress and strains states, need to be considered, the general concept of energy dissipation introduced previously remains the similar.
In the present notebook we apply the damage model inducing material softening to a two-dimensional case of stable crack propagation. In particular, we are going to use the three-point bending test as a means of rendering an elementary crack-propagation scenario, namely a straight crack propagating through a cross section. After discussing the test setup, we provide the test results in terms of load-deflection curve. Then, we will construct a numerical model that can reproduce such a behavior. Similarly to the previous lecture, we will consider the link between fracture and damage concepts that are both included in any finite element computation involving strain-softening material behavior
An isolated tensile crack propagation can be initiated using a notched specimen. The most common configurations used to study the cracking behavior for a tensile crack denoted as mode I are the wedge splitting test and notched, three-point bending test. Both these tests aim to characterize the material behavior in terms of the softening law describing the relation between the tensile stress transmitted across the localization zone and the corresponding crack opening.
Due to its simplicity, three-point-bending test on a notched concrete has become a standard (RILEM) to determine the fracture energy $G_\mathrm{f}$ characterizing the cracking behavior of concrete. The test induces a single crack in the notched section propagating straight upwards from the notch in a stable manner. The energy is dissipated in a local region in the crack vicinity so that it can be directly ascribed to the area of emerging crack surface.
The numerical simulation of this model can be readily performed using the material model with the damage function derived in previous notebooks. An example of the geometry and boundary conditions of the three-point bending test is provided by the [Petersson (1982)](https://portal.research.lu.se/portal/files/4705811/1785208.pdf) test series using the following setup.
The numerical model simulating this test is assuming a plain stress described by the stress tensor with the enumeration of spatial dimensions using the index variables $a, b$ = [1,2]
$$
\sigma_{ab}
=
\left[
\begin{array}{cc}
\sigma_{xx} & \sigma_{xy} \\
\sigma_{yx} & \sigma_{yy}
\end{array}
\right]
$$
and $\sigma_{zz} = \sigma_{xz} = \sigma_{yz} = 0$. The finite element discretization in this model applies the symmetry condition at the middle section of the beam. Upon loading loading, the damage will localize
The exponential damage function with the two parameters $\kappa_0$ defining the onset of damage and $\kappa_\mathrm{f}$ controlling the slope of the exponential softening branch at the onset of damage is defined as follows.
To achieve a fast convergence of the nonlinear simulation, the tangential algorithmic stiffness will be activated by setting the parameter `D_alg = 1`. Note that `D_alg = 0` corresponds to the secant stiffness which requires large number of iterations but can cope with large stress redistribution within the simulated domain.
The important feature of the provided damage model is the possibility to evaluate the fracture energy, i.e. the amount of energy to achieve a zero stress state. This value can be accessed as a `G_f` property of the material model.
The scalar damage material model can combine different measures of equivalent strain with a different type of damage function.
The task of the equivalent strain norm is to transform the two-dimensional strain tensor
$$
\varepsilon_{ab}
=
\left[
\begin{array}{cc}
\varepsilon_{xx} & \varepsilon_{xy} \\
\varepsilon_{yx} & \varepsilon_{yy}
\end{array}
\right]
$$
to a scalar value $\kappa$. The equivalent strain norms already mentioned in the notebook [5.1](../tour5_damage_bond/5_1_Introspect_Damage_Evolution_Damage_initiation.ipynb#strain_norm) play the role of an elastic limit. Their visual inspection
can be performed using the model component of the material model `MATSScalarDamage`.
Further strain norms can be chosen from the options
- Rankine strain norm
- Masar's strain norm
- Energy norm
These norms have to be combined with an elastic threshold to distinguish the elastic and inelastic domains.
Inspect the visual representation of the equivalent strain measure by changing the selector in the material model app. (It is necessary to select the respective tree node to render the visualization of the currently selected option.)
where $\varepsilon_I$ denotes the principal strains. The second version of the expression uses the index notation with a summation over an index that repeats in a product of two terms (see [Einstein summation rule](https://en.wikipedia.org/wiki/Einstein_notation)). The symbol $\delta_{II}$ represents an identity mapping called Kronecker delta, with $\delta_{IJ} = 1$ if $I=J$ and $\delta_{IJ} = 0$ otherwise.
The rate of convergence strongly depends on the quality of prediction of the stress state for the next step. The parameter `D_alg` can be used ton control the calculation by choosing between the secant and algorithmic, tangential stiffness. Derivation of the material stiffness might appear tedious at the first sight but in fact, it only requires the evaluation of the chain derivatives involved in the constitutive law.
Let us exemplify the derivation of the algorithmic material stiffness on the example of the isotropic damage model. In contrast to the version presented in notebook [5.1]() for the bond-slip law, we are applying the damage model to a 2D continuum described by the strain and stress tensor tensors $\boldsymbol{\varepsilon} = \varepsilon_{ab}$
and $\boldsymbol{\sigma} = \sigma_{ab}$, where letter indexes $a,b,c,d = [1,2]$ enumerate the orthogonal directions of the cartesian coordinate system. For example, $\varepsilon_{11} = \varepsilon_{xx}$. Then the stress-strain relation prescribed by the damage model is given as
where $\kappa(\varepsilon_{cd})$ represents the equivalent strain described above, and $\omega(\kappa)$ represents the particular damage function. The algorithmic stiffness is derived as the derivative of the above stress tensor with respect to a strain tensor. This means that every stress component must be differentiated with respect to every strain component. The result of this expression is a rank 4 tensor $D_{abcd}$ describing the instantaneous stiffness of a material point undergoing damage. The resulting expression delivers the material stiffness in the following form
The parameter `D_alg` is denoted here as $\theta_{\mathrm{alg,stiff}}$ and is provided as an input parameter in the material model. By setting this parameter to zero, the second term in the above expression is canceled. Then, only the secant stiffness is available. The importance of the algorithmic stiffness can be studied in the simulation of the bending test. The default value $\theta_\mathrm{alg,stiff} = 1$ is activating the algorithmic stiffness. By setting $\theta_\mathrm{alg,stiff} = 0$ one can see the strong reduction of the performance.
However, in some situation, for example when simulating simultaneous crack propagation in reinforced specimens the tangential stiffness might induce divergence. The reason for such behavior is the fact that the bond stress can rapidly change the sign. In these situations, advanced control strategies are required.
If a standard finite element discretization is combined with a softening material model,
1. it will always induce localization into a band of finite elements at the propagating stress concentration.
2. the area below the stress strain curve of a material model represents the energy dissipated by a unit volume of a material.
As a consequence, the total amount of the dissipated energy can be calculated as a product of the specific dissipation energy $G_\mathrm{f}$ $\mathrm{[J/mm^3]}$ and the volume of elements constituting the crack band $V_\mathrm{diss}$ $\mathrm{[mm^2]}$. In general,
$$
G_\mathrm{total} = G_\mathrm{f} V_\mathrm{diss}
$$
In a notched bending test, the crack band is running from the notch straight through the cross section. Therefore, it is possible to directly calculate the dissipative volume $V_\mathrm{diss}$ as
The total dissipated energy is available in the `hist` model component as a last element of the recorded array, so that we can pick it using the index `-1`
By changing the element size within the crack band we change the volume of material which dissipates energy. As observed, a larger element size leads to a larger amount of energy dissipation which makes the structural response more ductile. Why does it happen? The displacement approximated within a finite element is continuous function. Using low-order elements, with, e.g. bilinear shape functions, means that the strain within the element is of even lower order, i.e. constant or linear. Therefore, strain cannot localize within a smaller size than the size of the element through which the crack propagates. As a consequence, upon softening, an element exhibits a nearly uniform increasing strain which results in energy dissipation over its whole volume. The obtained results are therefore biased by the mesh size which documents an ill-posed model behavior.
There are several techniques how to avoid this mesh dependency. The principle idea of these so called regularization methods is to ensure a prescribed amount of energy dissipation upon localization. The most simple technique is sketched and demonstrated in the following parametric study.
The finite element codes solve this problem by adjusting the slope of the softening branch to keep the energy dissipation within a crack band invariant for a changed size of an element. This is the most common technique used in commercial finite element codes.
The above examples explain that the correspondence between the local value of $G_\mathrm{f}$ and the total energy $G_\mathrm{total}$ is established by realizing that energy dissipation is linked with the crack propagation. Indeed, if $G_\mathrm{f}$ is related to a unit area of the crack surface, the total energy dissipation in the whole test with a localized crack of the area $V_\mathrm{diss}$ is
In the above example, the crack was represented by a band of elements with the width specified by the $L_\mathrm{cb}$ model parameter. Since the crack developed over the whole notched cross section, we can express the dissipative volume as
To ensure that the amount of dissipated energy for a changed $V_\mathrm{diss}$ due to different mesh size does not change, we need to correspondingly adjust effective value of dissipated energy in depending on $L_\mathrm{cb}$. In other words, we need to introduce
This scaling must be achieved by adjusting the material parameters of the model depending on the size of the discretization element $L_\mathrm{cb}$. In present model we use the exponential damage function with the two parameters $\kappa_0$ and $\kappa_\mathrm{f}$ controlling the onset of damage and the slope of the descending branch, respectively. The above specified scaling of energy dissipation can be achieved by scaling the parameter
To verify, that this scaling correspondingly affects the effective amount of energy dissipation by evaluating the value of $G_{\mathrm{f}}^{L_\mathrm{cb}}$ for the values of $L_\mathrm{cb} = [1,2,3,4]$. These values should lead to the corresponding fractions [$1, \frac{1}{2}, \frac{1}{3}, \frac{1}{4}$] of $G_\mathrm{f}$. Assuming the reference value $\kappa_f$ corresponding to the length $L_\mathrm{cb} = 1$ we obtain
Even though the difference in the predicted peak of the load deflection curve shown left is still apparent, the amount of dissipated energy up to the control displacement of $w=2$ mm does not change with the element size any more. Let us remark that the applied scaling of $\kappa_\mathrm{f}$ was only approximate. A more sophisticated scaling and a smaller size of load increment would be required to achieve even a better match of the load deflection curves.
The described kind of scaling can be found in the most common finite-element codes. In the present example, the position and profile of the crack was known a-priori. Finite element codes implementing material models with softening behavior must implement a general procedure for the described crack band regularization. This regularization is performed in each material point by including the size of the containing element and the orientation of principle strains to reflect the width of the softening zone correctly.
The described procedure can be regarded as a _quick fix_. It makes the standard finite-element method applicable to simulation cracking for problems with straight crack evolution, i.e. tensile cracks. It is, however, insufficient for crack propagation with non-negligible shear stress within the localization zone. Indeed, most of the commercially available software based on continuous finite-element approximation of the displacement field cannot correctly predict the shear crack propagation. Pragmatic solutions, including an additional parameter like a shear retention factor, are provided in some software packages (ATENA).
More advanced methods address the problem of mesh independence using
- Nonlocal averaging of state variables over a defined size of a process zone
- Gradient methods of damage and plasticity including the gradient of state variables into the constitutive law
They are mentioned here for completeness and as an optional further reading.
**Note the dimensionality of $G_\mathrm{f}$:** To be precise, let us note that in contrast to the example shown in notebook [6.3](tour6_energy/6_3_localized_energy_dissipation.ipynb) for the CFRP pullout simulation, the damage function was applied to the bond-slip law, describing the behavior of a unit area of the interface. On the other hand, in the present notebook, we applied the damage function in a two-dimensional constitutive law.
