Development of the ENDO_HETEROGENE law
=======================================

In *Code_Aster*, we made the choice to represent the formation and propagation of cracks by an integrated isotropic damage model. The damage-elasticity coupling law being given classically by:


.. math::
 :label: eq-3

    \ underline {\ sigma} = (1-d)\ underline {\ underline {C}}\ underline {\ epsilon}

With :math:`\stackrel{ˉ}{\sigma }` the stress tensor, :math:`\underline{\varepsilon }` the deformation tensor, :math:`\underline{\underline{C}}` the isotropic elasticity tensor, and :math:`d` the damage variable.

We place ourselves in the framework mentioned in section :ref:`2 <RefNumPara__31404055>`. The presence of cracks in the structure is modelled by lines of broken elements (i.e. for which the damage d is such as :math:`d=1`). For healthy elements, :math:`d=0`. The rupture of the elements can be caused by two phenomena that are distinct here; namely the initiation of a new crack from the element in question or the propagation through the element in question of a crack that has begun in another element. These two types of breakage are described by two different types of threshold functions in this approach.

Moreover, for laws of behavior corresponding to fragile materials, the introduction of a heterogeneity of the material projected onto the mesh has the effect of introducing a dependence on the size of the zone (even with the use of classical non-local models [:ref:`bib1 <bib1>`]). In addition, a local damage law cannot describe fragile weakest link behavior.

In order to remedy this problem, it is proposed here to:


1. modify local law to represent a weakest member-type fragile fracture in a heterogeneous material


2. introduce a regularization operator based on the stress gradient.

It is shown that a non-local stress gradient description is, on the one hand, very close to non-local deformation gradient models and on the other hand allows us to correctly calculate the stress intensity factor, which is the key parameter for describing the propagation of macro-cracks. A regularized constraint is introduced such as


.. math::
 :label: eq-3

    \ underline {\ stackrel {} {\ sigma}}} - {l} _ {c} ^ {2}\ Delta\ underline {\ stackrel {} {\ sigma}} =\ underline {\ sigma}} =\ underline {\ sigma}}

With :math:`\underline{\stackrel{ˉ}{\sigma }}` the regularized stress tensor and :math:`{l}_{c}` the characteristic length. We will consider the so-called natural conditions at the Neumann limits, namely


.. math::
 :label: eq-3

    (\nabla\ underline {\ stackrel {} {\ sigma}})\ mathrm {.} n=0

with :math:`n` the outgoing normal at the considered edge. We consider a perfectly fragile law of behavior of the type:


.. math::
 :label: eq-3

    D=h\ left (\ langle\ stackrel {tante} {{\ sigma} _ {1}} - {\ sigma} _ {y} (\ mathit {el})\ rangle\ right),\ dot {d}\ ge 0

Where :math:`{\sigma }_{y}` is a stress associated with the elastic limit.

The two thresholds will therefore be written according to a regularized constraint over a characteristic length :math:`{l}_{c}` which will be large compared to the size of the elements. This section will focus on justifying the choice of the two thresholds and the way in which they coordinate.




Priming cracks
---------------------

As said above, the initiation of new macrocracks is represented by a Weibull model. In order to describe the cracking step in a relevant manner, a dispersion on the parameters of the material (Weibull module, scale parameter) specific to each finite element is introduced. By first drawing lots on a uniform distribution for each finite element in order to derive a value of probability of failure for each element :math:`P(\mathrm{el})\in [:ref:`0;1 <0;1>`] `, we can describe the stress for which a crack will start in a finite element by the following equation:


.. math::
 :label: eq-3

    {\ sigma} _ {a} (\ mathrm {el}) =\ frac {{\ sigma} _ {0}} {{({\ lambda} _ {0} {Z}} _ {\ mathrm {el}}) =\ frac {{el}}) =\ frac {{\ sigma} _ {\ lambda} _ {0} {Z} _ {\ Z} _ {\ mathrm {el}} _ {\ mathrm {el}})} {\ mathrm {el}})} ^ {1/m}}\ mathrm {.} {\ left [-\ mathrm {ln} (1-P (\ mathrm {el}))\ right]} ^ {\ frac {1} {m}}

With :math:`{Z}_{\mathrm{el}}` the size (the volume or the area depending on the dimension) of the finite element, :math:`m` the Weibull modulus and :math:`\frac{{\sigma }_{0}^{m}}{{\lambda }_{0}}` the scale parameter.

The use of a Weibull model to describe the initiation of new cracks has a definite advantage in order to obtain results independent of the mesh. Let's consider two different partitions of domain :math:`{Z}_{\mathrm{do}}` each made up of sub-domains :math:`{Z}_{\mathrm{di}}^{1}` and :math:`{Z}_{\mathrm{di}}^{2}`. We can write:


.. math::
 :label: eq-3

    \ sum _ {i=1} ^ {\ mathrm {n 1}} {1}} {Z}} _ {\ mathrm {di}} ^ {1} = {Z} _ {\ mathrm {do}}} =\ sum _ {i=1}} {\ mathrm {1}}} {Z} _ {\ mathrm {di}}} =\ sum _ {i=1}} =\ sum _ {i=1}} ^ {1} ^ {2} ^ {2}} ^ {2}

where :math:`\mathit{n1}` and :math:`\mathit{n2}` mean subdomain numbers in both cases. The probability of initiating the first crack in a subdomain for a :math:`\sigma` request is written:


.. math::
 :label: eq-3

    {P} _ {\ mathrm {asdi}} =1-\ mathrm {exp}}\ left [- {Z} _ {\ mathrm {sd}} {\ lambda} _ {0} {{0}} {(\ frac {\ sigma}} {\ sigma} _ {0}})} ^ {m}\ right]

where :math:`{Z}_{\mathrm{sd}}` is the size (volume) of a subdomain. The probability of initiating the first crack on all sub-domains :math:`{P}_{\mathrm{aesd1}}` for a :math:`\sigma` stress can be deduced from this:


.. math::
 :label: eq-3

    {P} _ {\ mathrm {aesd1}} =1-\ prod _ {i=1} ^ {{n} _ {\ mathrm {sd}}}\ left [1- {P} _ {\ mathrm {asdi1}}} =1-\ mathrm {asdi1}}}\ right] =1-\ mathrm {exp}\ left [-\ sum _ {i=1} _ {\ mathrm {asdi1}}} {\ right] =1-\ mathrm {asdi1}}}\ right] =1-\ mathrm {exp}\ left [-\ sum _ {i=1} _ {1}} {Z} _ {\ mathrm {di}}} ^ {di}}} ^ {1} {1} {\ lambda} _ {0}})} ^ {m}\ right]} ^ {m}\ right] =1-\ mathrm {exp}\ right] =1-\ mathrm {exp}\ left [-\ sum _ {i=1}} ^ {{n} _ {2}} {2}\ right] =1-\ mathrm {exp}\ right] =1-\ mathrm {exp}\ left [-\ sum _ {i=1} ^ {{n} _ {2}} {2}}\ right] =1-\ mathrm {exp}\ right] =1-\ mathrm {exp}\\ mathrm {di}} ^ {2} {\ lambda} {\ lambda} _ {0} _ {\ sigma} {{\ sigma} _ {0}})} ^ {m}\ right] = {P}\ right] = {P} _ {\ mathrm {aesd2}}

Because:


.. math::
 :label: eq-3

    \ sum _ {i=1} ^ {{n} _ {1}} {1}} {1}} {Z}} _ {\ mathrm {di}} ^ {i=1} ^ {{n} _ {2}} {2}} {Z}}} {Z}} {2} = {Z} _ {\ mathrm {do}}}

The equation allows us to find a probability of the first crack starting on a domain of size :math:`{Z}_{\mathrm{do}}` independent of the division into sub-domains.

Finally, we can describe the law of evolution of damage d linked to the initiation of a new crack in an element:


.. math::
 :label: eq-3

    D=h (\ langle\ stackrel {} {} {{\ sigma} _ {1}} - {\ sigma} _ {a} (\ mathrm {el})\ rangle),\ dot {d}\ ge 0

where :math:`H` is the Heaviside step function and :math:`\stackrel{ˉ}{{\sigma }_{1}}` is the maximum principal regularized stress.




Fissure propagation
------------------------

In this part, we will demonstrate that the proposed non-local approach allows us to correctly calculate the stress intensity factor for cracks modelled as a series of broken elements. In fact, it is possible to link a criterion in terms of stress intensity factor, resulting from the mechanics of rupture to the regularized stress. It can be considered that the asymptotic solution of the stress field at the crack point is a good approximation of the fields that can be found at the tip of a propagating crack. The stress field given by the Westergaard solution can be written in a crack point coordinate system as shown in the Figure as a composition of elastic fields (for a flat problem):


.. math::
 :label: eq-3

    \ underline {\ underline {\ sigma}} (r,\ theta) = {K} _ {I}\ underline {\ underline {{f} _ {I}}} (r,\ theta) + {K} _ {\ mathit {II}}} (r,\ theta) = {K} _ {\ mathit {II}}}} (r,\ theta)

With :math:`{K}_{I},{K}_{\mathrm{II}}` the stress intensity factors in :math:`I` and :math:`\mathit{II}` mode, and :math:`\underline{\underline{{f}_{I},{f}_{\mathrm{II}}}}` the associated weight functions (equation and) in :math:`I` and :math:`\mathit{II}` mode.


.. math::
 :label: eq-3

    \ underline {\ underline {{f} _ {I}}}} (r,\ theta) =\ frac {1} {\ sqrt {2\ pi r}}\ left [\ begin {array} {ccc}\ mathrm {cos}}} (r,\ theta)}}} (r,\ theta) =\ frac {\ theta} {2}}\ left (1-\ mathrm {sin}\ frac {\ theta} {2}\ mathrm {sin}\ frac {3\ theta} {2}\ right) &\ mathrm {sin}\ frac {\ theta} {2}\ mathrm {cos}\ frac {\ theta} {2}\ frac {\ theta} {2}\ frac {\ theta} {2} & 0\\\ mathrm {sin}\ frac {sin} {2}\ frac {\ theta} {2} & 0\\\ mathrm {sin}\ frac {\ theta} {2} eta} {2}\ mathrm {cos}\ frac {\ theta} {2}\ mathrm {theta} {2}\ mathrm {cos}\ frac {\ theta} {2}\ frac {\ theta} {2}\ left (1+\ mathrm {sin}\ frac {\ theta} {2}\ mathrm {sin} {2}\ mathrm {sin}\ frac {sin} {2}\ mathrm {sin}\ frac {sin} {2}\ mathrm {sin}\ frac {sin} {2}\ mathrm {sin}\ frac {sin} {2}\ mathrm {sin}\ frac {sin} {3\ theta} {2}\ right) & 0\\ 0& 0& 0\ end {array}\ right]


.. math::
 :label: eq-3

    \ underline {\ underline {{f} _ {\ mathit {II}}}}}} (r,\ theta) =\ frac {1} {\ sqrt {2\ pi r}}\ left [\ begin {array} {{array} {ccc}} -\ mathrm {II}}}}} (r,\ theta) =\ frac {1} {\ sqrt {2\ pi r}}}\ left (2-\ mathrm {cos}}\ frac {cos} {ccc} -\ mathrm {ccc} -\ mathrm {cc}} -\ mathrm {\ sin}}}} (r,\ theta) =\ frac {1} {2}\ left (2-\ mathrm {cos}}\ frac {theta} {2}\ mathrm {cos}\ frac {3\ theta} {2}\ right) &\ mathrm {cos}\ frac {\ theta} {2}\ left (1-\ mathrm {sin}\ frac {sin}\ frac {\ theta} {2}\ right) & 0\ theta} {2}\ right) & 0\\\ mathrm {cos}\ frac {\ theta} {2} {2}\ left (1-\ mathrm {sin}\ frac {\ theta} {2}\ mathrm {sin}\ frac {3\ theta} {3\ theta} {2}\ right) &\ mathrm {\ theta} {2}\ mathrm {cos}\ frac {cos}\ frac {2}\ frac {\ theta} {2}\ mathrm {sin}\ frac {3\ theta} {2} & 0\\ 0& 0& 0\ end {array}\ right]


.. image:: images/10000200000002E7000002746F467E25E2D3604F.png
    :width: 2.7917in
    :height: 2.2307in


.. _RefImage_10000200000002E7000002746F467E25E2D3604F.png:

**Figure** 3.2-a **: Definition of a crack point location**

The regularized stress at the crack point can be estimated by


.. math::
 :label: eq-3

    \ underline {\ underline {\ stackrel {249} {{\ sigma} {{\ sigma}} _ {p}}}} ({K} _ {\ mathit {II}}) =\ sum _ {i=1}} ^ {\ mathit {1}} ^ {1}} ^ {1} ^ {1} ^ {1} ^ {2} ^ {2}}\ left [\ frac {{K} _ {i}}} {\ pi {l}}}) =\ sum _ {i=1}}} {\ sum _ {i=1}} ^ 1} ^ {1}} ^ {2} ^ {2}} {\ int} _ {\ frac {-\ pi} {2}}} ^ {\ frac {\ pi} {2}} {\ int} _ {0} ^ {{\%} _ {c}}}\ underline {\ underline {\ pi}}}\ underline {\ pi}}}\ underline {\ underline {\ pi}}}\ underline {\ underline {{f} _ {f}}} (r,\ theta) r\ partial r\ partial\ theta\ right]

Which can finally be written in the following form:


.. math::
 :label: eq-3

    \ underline {\ underline {\ stackrel {ounty} {{\ sigma} {{\ sigma}} _ {p}}}} ({K} _ {\ mathit {II}}) =\ frac {{\ Gamma}}}) =\ frac {{\ Gamma}} ^ {2} ^ {2} (3/4)} {5\ pi\ sqrt {\ pi\ mathit {lc}}}) =\ frac {{\ Gamma}}}) =\ frac {{\ Gamma}} ^ {2} ^ {2} (3/4)} {5\ pi\ sqrt {\ pi\ mathit {lc}}}}\ left [\ begin {array} {ccc}}) =\ frac {{\ Gamma}}}) =\ frac {\} {\ mathrm {4K}} _ {I} & {\ mathrm {4K}}} _ {\ mathit {II}}} & 0\\ {\ mathrm {4K}} _ {\ mathit {II}}} & {\ mathit {II}}} & {\ mathrm {6K}}} _ {\ mathrm {6K}}} _ {\ mathrm {6K}}} _ {I}} & 0\\ 0&\ stackrel {249}} {{\ sigma}} _ {\ mathrm {6K}}} _ {\ mathrm {6K}}} _ {\ mathrm {6K}}} _ {\ mathrm {6K}}} _ {\ mathrm {6K}}} _ {m {33p}}}\ end {array}\ right]

where :math:`\stackrel{ˉ}{{\sigma }_{\mathrm{33p}}}=-\nu (\stackrel{ˉ}{{\sigma }_{\mathrm{11p}}}+\stackrel{ˉ}{{\sigma }_{\mathrm{22p}}})` in plane deformation and :math:`\stackrel{ˉ}{{\sigma }_{\mathrm{33p}}}=0` in plane stresses. Moreover, we recall that :math:`\Gamma` is the Gamma function, such as :math:`\Gamma (a)=\gamma (a,\infty )` and :math:`\gamma (a,b)={\int }_{0}^{b}{t}^{a-1}\mathrm{exp}(-t)\mathrm{dt}` incomplete gamma function. Concretely :math:`\Gamma (3/4)=\mathrm{1,225416}`.

Then, based on the expression of the regularized stress at the crack point, it is possible to propose a law of evolution associated with the propagation of a macro-crack relating to the maximum main stress of the type:


.. math::
 :label: eq-3

    D=h\ left (\ langle\ stackrel {tante} {\ sigma} {{\ sigma} _ {I}} - {\ sigma} _ {p} (\ mathit {el})\ rangle\ right),\ dot {d}\ ge 0

With:


.. math::
 :label: eq-3

    {\ sigma} _ {p} (\ mathit {el}) =\ frac {6 {\ Gamma} ^ {2} (3/4)} {5\ pi}\ frac {{K} _ {\ mathit {Ic}}} {\ mathit {Ic}}}} {\ sqrt {\ pi {l} _ {c}}}

:math:`{\sigma }_{p}(\mathrm{el})` is the breaking point of an element by the propagation of a crack at the macroscopic scale and :math:`{K}_{\mathrm{Ic}}` is the toughness of the material. For a load in pure I mode, the proposed threshold corresponds to a threshold of the stress intensity factor type resulting from fracture mechanics. The proposed threshold has in particular the advantage of being independent of the fact that one uses plane stresses or plane deformation for a 2D problem and it is similar to a criterion in terms of stress intensity factor in mode I for a 3D problem. It should be noted that in practice, the criterion presented is a mixed mode criterion, in fact, if we observe the threshold to break by propagation of a crack at the macroscopic scale under a load composed, for example, of modes I and II in plane stresses, we obtain at the crack point.


.. math::
 :label: eq-3

    \ stackrel {279} {{\ sigma} _ {\ mathit {Ip}}}} =\ frac {6 {\ Gamma} ^ {2} (3/4)} {5\ pi}\ frac {1} {\ frac {1} {\ sqrt {1}} {\ sqrt {\ sqrt {\ pi}} {\ sqrt {\ pi}} {1} {\ sqrt {\ pi}} {1} {\ sqrt {\ pi}} {\ sqrt {\ pi}} {\ sqrt {\ pi}} {l} _ {I}}}\ left [{\ mathrm {5K}}}\ left [{\ mathrm {5K}}} _ {I} +\ sqrt {{K}} _ {I}} {I}} ^ {2} + {\ mathrm {16K}}} _ {\ mathit {II}} ^ {2}}\ right]

which also depends on the level of demand in :math:`\mathit{II}` mode.




Coordination of thresholds
-----------------------

The two-parameter Weibull model describes crack initiation with a very conservative character. In fact, using a Weibull model with two parameters is equivalent to assuming the existence of a non-zero probability of crack initiation as soon as the maximum principal stress becomes positive. On the other hand, the use of a crack propagation criterion as a factor of stress intensity on a crack at the microscopic scale is equivalent to assuming the existence of a non-zero initiation threshold. For example, consider the case of an area that has a dimension of the order of :math:`{l}_{c}` and is included in a larger area that is stressed uniformly. If we assume the existence of a Weibull :math:`{\sigma }_{a}` initiation threshold, we can associate a crack length with it using a fracture mechanics criterion:


.. math::
 :label: eq-3

    {\ sigma} _ {a} (\ mathit {el}) =\ frac {{el}) =\ frac {{\ sigma} _ {0} {Z} _ {\ mathit {el}})}}} ^ {1/m})}} ^ {1/m}}} {[-\ mathrm {ln}} {{\ lambda}} _ {0} {Z} _ {Z} _ {\ mathit {el}}})} {\ mathit {el}}))} {\ mathit {el}}))}} ^ {el}}))} ^ {el}}))} ^ {el}})} ^ {el}})} ^ {el}})} ^ {el}})} ^ {el}})} ^ {el}})} {{K} _ {\ mathit {Ic}}}} {\ sqrt {\ pi a (\ mathit {el})}}

With :math:`a(\mathrm{el})` the half-length of the largest crack in the area under consideration. For the sake of consistency with the propagation threshold :math:`{\sigma }_{p}`, care should be taken to ensure that the use of a Weibull model to describe crack initiation does not involve assuming the presence in a zone under consideration of a crack larger than the size of the zone under consideration. Moreover, it seems reasonable to assume that a given element will be more easily broken by the propagation of a macro-crack than by the initiation of a new crack. In this way, we will also ensure that the priming sites that are in the vicinity of the tip of a crack will be screened by it. This hypothesis amounts to posing the following inequality:


.. math::
 :label: eq-3

    {\ sigma} _ {a} (\ mathit {el})\ ge {\ sigma} _ {p} (\ mathit {el})

which is equivalent to:


.. math::
 :label: eq-3

    a (\ mathit {el})\ le {\ left (\ frac {5\ pi} {6\ Gamma (3/4)}\ right)} ^ {2} {l} {l} _ {c}

to define a critical value for the initial defect size. The preceding inequality can be interpreted as follows: it cannot be assumed that an initial crack of comparable size at the macroscopic scale is contained in an area of comparable size at the microscopic scale. In practice, if there are initial cracks with a half-length greater than approximately :math:`{l}_{c}`, they cannot be taken into account by the proposed priming model; on the other hand, they can be represented by a line of broken finite elements. The Weibull model will therefore be used to represent the initiation of cracks from initial defects smaller than :math:`2{l}_{c}`. The characteristic length :math:`{l}_{c}` appears as a transition scale between the microscopic scale of crack initiation and the macroscopic scale of crack propagation. The characteristic length can then be interpreted as being of the order of magnitude of the largest cracks initially present at the microscopic scale.

Moreover, we note that if each finite element does not contain any initial defects larger than itself and if the characteristic length :math:`{l}_{c}` is greater than the size of the finite elements, then the inequality is automatically verified.


.. _Toc268168108:


.. _Toc242524034:

To coordinate the thresholds, it is necessary to identify the elements. To do this, a criterion relating to the neighborhood of a finite element has been developed. In this context, the propagation criterion will only be verified on elements that are in the direct vicinity of a crack that has already begun. This model therefore requires the knowledge of the neighbors of an element. The details of this tracking are explained in the next chapter.