# 7.2 Softening function and damage function


```python
%matplotlib widget
import sympy as sp
import numpy as np
import matplotlib.pylab as plt
sp.init_printing()
```

Let us consider a softening function of the form
\begin{align}
f(w) = c_1 exp( - c_2 w )
\end{align}
This is function should describe the decay of stress starting from the material tensile strength and continuously deminishing to zero.
The variable $w$ represents the crack opening and the parameters $c_1$ and $c_2$ are the material parameters.


```python
c_1, c_2 = sp.symbols('c_1, c_2', positive=True)
w = sp.symbols(r'w')
```


```python
f = c_1 * sp.exp(-c_2*w)
```

Let us plot the function to verify its shape for the material parameters set to the value 1


```python
f_w = f.subs({c_1:1, c_2:1})
f_w
```


```python
w_arr = np.linspace(0.01,10,100)
fig.canvas.header_visible=False
fig, ax = plt.subplots(1,1,figsize=(5,3), tight_layout=True)
ax.plot(w_arr, sp.lambdify(w, f_w)(w_arr));
```


    Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …


The function can be already used in this form. The question is however, how to set the material parameters $c_1$ and $c_2$. They can be directly associated to a particular type of material parameters - namely, to the tensile strength $f_\mathrm{t}$ and fracture energy $G_\mathrm{f}$ 


```python
f_t, G_f = sp.symbols('f_t, G_f', positive=True)
```

Softening starts at the level of the material strength so that we can set $f(w = 0) = f_\mathrm{t}$  to obtain the equation


```python
Eq1 = sp.Eq(f.subs({'w':0}), f_t)
Eq1
```

By solving for the $c_1$ we obtain the first substitution for our softening function.


```python
c_1_subs = sp.solve({Eq1}, c_1)
c_1_subs
```

Thus, if $w_0 = 0$, $c_1$ is equivalent to the tensile strength $f_\mathrm{t}$

The second possible mechanical interpretation is provided by the statement that softening directly represents the energy dissipation of a unit crack area. Thus, for large $w \rightarrow \infty$ it is equivalent to the energy producing a stress-free crack. This is the meaning of fracture energy.

We can thus obtain the fracture energy represented by the softening function by evaluating its integral in the range $w \in (0, \infty)$. 


```python
int_f_w = sp.integrate(f.subs(c_1_subs), w)
int_f_w
```

As $c_2 > 0$, only the second term matters.
The determinate integral
\begin{align}
\left[ - \frac{f_\mathrm{t}}{c_2} 
\exp(-c_2 w) \right]_0^{\infty}
\end{align}
is zero for $w = \infty$, so that the value in $w = 0$ delivers the result of the integral
\begin{align}
\frac{f_\mathrm{t}}{c_2} 
\end{align}

This integral is equal to the fracture energy $G_\mathrm{f}$.


```python
E_0 = sp.symbols('E_0', positive=True ) 
```


```python
Eq2 = sp.Eq(-int_f_w.subs({'w':0}), G_f)
Eq2
```

and the value of $c_2$ delivers the second substitution for the softening function


```python
c_2_subs = sp.solve({Eq2}, c_2)
c_2_subs
```

The softening function with strength and fracture energy as a parameter now obtains the form


```python
f_w = f.subs(c_1_subs).subs(c_2_subs)
sp.simplify(f_w)
```

Verify that the fracture energy is recovered at $w$ in infinity


```python
sp.integrate(f_w, (w,0,sp.oo))
```

## What if there is a finite elastic stiffness


```python
E_0, w_0 = sp.symbols(r'E_b, w_0', positive=True)
```


```python
f = c_1 * sp.exp(-c_2*w)
```


```python
f_w_0 = sp.Piecewise(
    (1, w < w_0),
    (f.subs(w, w-w_0), True)
)
f_w_0
```


```python
sig_w = E_0 * w * f_w_0
sig_w
```


```python
subs_c_1 = sp.solve( {sp.Eq( sig_w.subs(w, w_0), E_b * w_0 )}, {c_1} )
subs_c_1
```


```python
G_f_ = sp.simplify( sp.integrate(sig_w.subs(subs_c_1), (w,0,sp.oo)) ).factor()
G_f_
```


```python
_, subs_c_2 = sp.solve( sp.Eq( G_f, G_f_ ), {c_2} )
sp.simplify(subs_c_2)
```


```python
g1_w_ = sp.simplify( f_w_0.subs(c_2, subs_c_2).subs(c_1,1) )
g1_w_
```

\begin{align}
g(\tilde{s}) = \displaystyle{\frac{w_0}{\tilde{w}}}
\exp\left(\frac{2 E_\mathrm{b} w_0 ( \tilde{w} - w_0) }
{E_\mathrm{b} w^2_0 - 2 G_\mathrm{f}} 
 \right)
\end{align}


```python
g2_w_ = sp.Piecewise(
    (1, w < w_0),
    ( w_0 / w * sp.exp( 
    (2*E_b*w_0*(w-w_0))/(E_b*w_0**2-2*G_f) 
    ), True)
)
g2_w_
```


```python
get_g1_w_ = sp.lambdify( (w, w_0, E_b, G_f), g1_w_ )
get_g2_w_ = sp.lambdify( (w, w_0, E_b, G_f), g2_w_ )
```


```python
get_g1_w_(0.1, 0.1, 10, 0.1)
```




    array(1.)




```python
fig, ax = plt.subplots(1,1)
ax2 = ax.twinx()
w_range = np.linspace(0.0001,1,100)
ax.plot(w_range, get_g1_w_(w_range, 0.1, 10, 0.1), color='red')
ax.plot(w_range, get_g2_w_(w_range, 0.1, 10, 0.1), color='blue')
ax2.plot()
```


```python
#f_t = E_b * s_0
#g30_s_ =  f_t / E_b / s_equiv * sp.exp( - f_t / G_f_0 * (s_equiv - s_0) )
#g3_s_ = sp.simplify(g30_s_.subs(G_f_0, G_f - E_b*s_0**2/2))
```

## How to apply the softening function to a zone of length $L_s$?

If we wish to embed the softening behavior into a finite element simulation we encounter the problem that the deformation is not described by the crack opening $w$. The finite element discretization assumes by definition a smooth stran field $verepilon$ and there is no notion of discontinuity. 

As a consequence, we need to account for a length of the softening zone  $L_s$ which is actually equivalent to the size of the finite element.


```python
L_s, epsilon = sp.symbols('L_s, varepsilon')
```

The crack opening $w$ is simply substituted by
\begin{align}
 w = \varepsilon L_\mathrm{s}
\end{align}
to obtain the softening function in terms of strains


```python
f_epsilon = f_w.subs({'w' : epsilon * L_s})
f_epsilon
```

Note the important fact that the integral of the softening function over the strains scales the dissipated energy by the term $1/L_\mathrm{s}$. 


```python
sp.integrate(f_epsilon.subs(L_s,1), (epsilon,0,sp.oo))
```

As a consequence, by changing the size of the softening zone we also change the total amount of the energy dissipated!!! This feature of the softening function will be exploited in the implementation of the finite elements.

## How to convert a softening function to damage function?

The shape of the softening function describes the process of deterioration - starting from an undamaged state and ending in a fully damage state. Let us establish an equivalence between softening and damage evolution by requiring that they describe the same kind of stress decay.

Considering the softening zone of the length $L_s$ let us describe the damage within this zone using the state variable $omega$ and the elastic modulus $E_c$


```python
omega, E_c = sp.symbols('omega, E_c')
```

Then, the constitute law to be used in a finite element of the zone is given as


```python
sigma = (1 - omega) * E_c * epsilon
```

The stress decay described by the damage law and by the softening law should be the same. Let them set equal and solve for the damage variable $omega$. Using the sympy solver we obtain the algebraic solution as


```python
omega_Gf = sp.solve(sigma - f_epsilon, omega)[0]
omega_Gf
```


```python
sigma_v = (1- omega_Gf) * E_c * epsilon
sigma_v
```

This new damage function is defined using material parameters with a clear mechanical interpretation. Such kind of material law is attractive because it makes it possible to design tests that focus on an isolated phenomenon, i.e. the determination of the material strength, E modulus or fracture energy.  

Let us now visually verify the shape of the damage function.
The set of parameters is assembled in a dictionary and then
they are all substituted into the damage function.


```python
data_f = dict(E_c = 28000, f_t = 3, L_s = 1, G_f = 0.01)
omega_Gf_epsilon = omega_Gf.subs(data_f)
omega_Gf_epsilon
```

The damage function is only valid in an inelastic regime. Therefore, we have to quantify the onset of inelasticity first as
\begin{align}
\varepsilon_0 = \frac{f_t}{E_c}
\end{align}


```python
epsilon_0 = (f_t / E_c).subs(data_f)
epsilon_0
```

Then the damage function can be plotted as


```python
splot(omega_Gf_epsilon, (epsilon,epsilon_0,epsilon_0*100))
```

The corresponding stress strain curve has then the form


```python
sigma = sp.Piecewise( (E_c * epsilon, epsilon < epsilon_0),
    (( 1 - omega_Gf ) * E_c * epsilon, epsilon >= epsilon_0))
sigma
```


```python
splot(sigma.subs(data_f), (epsilon,0,epsilon_0*100))
```

Note how the stress strain function scales with the change of the fracture energy and of the zone length.

Larger fracture energy makes the stress-strain response more ductile, while smaller makes it brittle.

On the other hand, larger size of the softening zone makes the softening behavior more brittle and smaller size of the zone makes it more ductile. Why?


```python
g_f = sp.integrate(sigma, epsilon)
```


```python
g_f
```


```python

```

ecc3458a4182bfc62167edfd770ef9709b8b1ee6
ipynb
Jupyter Notebook
tour6_energy/7_2_Softening_law_with_embedded_fracture_energy.ipynb
bmcs-group/bmcs_tutorial
4e008e72839fad8820a6b663a20d3f188610525d
MIT
tour6_energy/7_2_Softening_law_with_embedded_fracture_energy.ipynb
bmcs-group/bmcs_tutorial
4e008e72839fad8820a6b663a20d3f188610525d
MIT
tour6_energy/7_2_Softening_law_with_embedded_fracture_energy.ipynb
bmcs-group/bmcs_tutorial
4e008e72839fad8820a6b663a20d3f188610525d
MIT
Qwen/Qwen-72B
1. YES
2. YES