Skip to content
Snippets Groups Projects
Commit a15f314a authored by Rostislav Chudoba's avatar Rostislav Chudoba
Browse files 

plasticity
parent a37bff6c
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
bmcs_course/2_2_PO_LF_LM_EL.ipynb
+ 3952
− 386
View file @ a15f314a
Source diff could not be displayed: it is too large. Options to address this: view the blob.
bmcs_course/4_2_BS_EP_SH_IK_A.ipynb 0 → 100644
+ 1207
− 0
View file @ a15f314a
Source diff could not be displayed: it is too large. Options to address this: view the blob.
bmcs_course/4_2_BS_elasto_plastic.ipynb deleted 100644 → 0
+ 0
− 694
View file @ a37bff6c
%% Cell type:markdown id: tags:
# 4.2 Bond behavior governed by plasticity
In the first part of this notebook, we rephrase the basic framework of
elasto-plastic models showing which conditions are used to find out how during yield
to describe the material behavior, once it crosses the elastic limit.
Perfect plasticity with constant level of yielding stress is considered first.
In the second part, hardening variable is included which allows for to expand
the elastic range
Define a bond-slip law governed yielding with hardening or softening.
%% Cell type:code id: tags:
``` python
%matplotlib inline
import matplotlib.pyplot as plt
import sympy as sp
sp.init_printing()
```
%% Cell type:markdown id: tags:
## The onset of inelasticity
%% Cell type:code id: tags:
``` python
s_el = sp.symbols('s_el')
E_b = sp.symbols('E_\mathrm{b}', positive=True)
tau_ = E_b * s_el
tau_
```
%% Output
$\displaystyle E_\mathrm{b} s_{el}$
E_\mathrm{b}⋅sₑₗ
%% Cell type:code id: tags:
``` python
s, s_pl = sp.symbols('s, s_pl')
s_el_ = s - s_pl
tau_.subs(s_el, s_el_)
```
%% Output
$\displaystyle E_\mathrm{b} \left(s - s_{pl}\right)$
E_\mathrm{b}⋅(s - sₚₗ)
%% Cell type:code id: tags:
``` python
tau = sp.symbols(r'\tau')
tau_bar = sp.symbols(r'\bar{\tau}')
f_tau = sp.sqrt( tau * tau ) - tau_bar
f_tau
```
%% Output
$\displaystyle - \bar{\tau} + \sqrt{\tau^{2}}$
_______
╱ 2
-\bar{\tau} + ╲╱ \tau
%% Cell type:code id: tags:
``` python
f_s = f_tau.subs(tau, tau_).subs(s_el, s_el_)
f_s
```
%% Output
$\displaystyle E_\mathrm{b} \sqrt{\left(s - s_{pl}\right)^{2}} - \bar{\tau}$
____________
╱ 2
E_\mathrm{b}⋅╲╱ (s - sₚₗ) - \bar{\tau}
%% Cell type:markdown id: tags:
The question is, how does the stress develop in the inelastic regime?
This means when $f = 0$?
%% Cell type:code id: tags:
``` python
s_pl_f0 = sp.solve(f_s, s_pl)
s_pl_f0[0]
```
%% Output
$\displaystyle s - \frac{\bar{\tau}}{E_\mathrm{b}}$
\bar{\tau}
s - ────────────
E_\mathrm{b}
%% Cell type:code id: tags:
``` python
tau_.subs(s_el, s_el_).subs(s_pl, s_pl_f0[1])
```
%% Output
$\displaystyle - \bar{\tau}$
-\bar{\tau}
%% Cell type:markdown id: tags:
So far, so good. We obtained a trivial result telling us that the stress cannot grow over the
introduced elasticity limit. But in this case we can just reproduce the constant bond-slip
law - now extended with a elastic stiffness. How to provide a general framework suitable
for an algorithmic treatment? How to distinguish between loading and unloading?
%% Cell type:markdown id: tags:
So how much plastic deformation was consumed? What was the direction of flow?
There might be several different processes going on the material structure?
Let us postulate, that they can be clustered by a positive flow variable $\lambda$. Then
we assume that the rate of change of all state variables representing yielding
can be related to the common positive parameter $\lambda$ and that the yielding direction
is normal to the yield surface $f$
%% Cell type:code id: tags:
``` python
lambda_ = sp.symbols(r'\lambda', nonnegative=True)
dot_s_pl_ = lambda_ * f_tau.diff(tau)
dot_s_pl_
```
%% Output
$\displaystyle \frac{\lambda \sqrt{\tau^{2}}}{\tau}$
_______
╱ 2
\lambda⋅╲╱ \tau
──────────────────
\tau
%% Cell type:markdown id: tags:
**Relate increment of yielding to the increment of primary kinematic state variables**
%% Cell type:markdown id: tags:
But what is the amount of yielding controlled now by $\lambda$? Solving for $f = 0$ is not
of much help any more. Can the elastic range change during the yielding process? An abstract
distinction between elastic and inelastic loading processes can be provided by Kuhn-Tucker conditions
\begin{align}
\lambda \dot{f} = 0, \; \lambda > 0,\; \dot{f} \le 0
\end{align}
%% Cell type:code id: tags:
``` python
dot_s, dot_s_pl = sp.symbols(r'\dot{s}, \dot{s}_\mathrm{pl}')
dot_tau_ = E_b * (dot_s - dot_s_pl)
dot_tau_
```
%% Output
$\displaystyle E_\mathrm{b} \left(\dot{s} - \dot{s}_\mathrm{pl}\right)$
E_\mathrm{b}⋅(\dot{s} - \dot{s}_\mathrm{pl})
%% Cell type:code id: tags:
``` python
dot_f = f_tau.diff(tau) * dot_tau_.subs(dot_s_pl, dot_s_pl_)
```
%% Cell type:code id: tags:
``` python
lambda_solved = sp.solve( dot_f, lambda_)[0]
lambda_solved
```
%% Output
$\displaystyle \frac{\dot{s} \tau}{\sqrt{\tau^{2}}}$
\dot{s}⋅\tau
────────────
_______
╱ 2
╲╱ \tau
%% Cell type:code id: tags:
``` python
dot_s_pl_.subs(lambda_, lambda_solved)
```
%% Output
$\displaystyle \dot{s}$
\dot{s}
%% Cell type:code id: tags:
``` python
dot_tau_.subs(dot_s_pl, dot_s_pl_).subs(lambda_, lambda_solved)
```
%% Output
$\displaystyle 0$
0
%% Cell type:markdown id: tags:
[Convex mathematical programming literature]
%% Cell type:markdown id: tags:
## Can the elastic range expand?
%% Cell type:code id: tags:
``` python
z = sp.symbols('z')
K = sp.symbols('K', positive=True )
Z = sp.symbols('Z')
f_tau = sp.sqrt(tau**2) - (tau_bar + Z)
f_tau
```
%% Output
$\displaystyle - Z - \bar{\tau} + \sqrt{\tau^{2}}$
_______
╱ 2
-Z - \bar{\tau} + ╲╱ \tau
%% Cell type:code id: tags:
``` python
dot_s_pl_ = lambda_ * f_tau.diff(tau)
dot_s_pl_
```
%% Output
$\displaystyle \frac{\lambda \sqrt{\tau^{2}}}{\tau}$
_______
╱ 2
\lambda⋅╲╱ \tau
──────────────────
\tau
%% Cell type:code id: tags:
``` python
dot_z_ = - lambda_ * f_tau.diff(Z)
dot_z_
```
%% Output
$\displaystyle \lambda$
\lambda
%% Cell type:markdown id: tags:
\begin{align}
\dot{f} =
\frac{\partial f}{\partial \tau} \dot{\tau}
+
\frac{\partial f}{\partial Z} \dot{Z}
\end{align}
%% Cell type:code id: tags:
``` python
dot_z = sp.symbols(r'\dot{z}')
dot_Z_ = K * dot_z
dot_Z_
```
%% Output
$\displaystyle K \dot{z}$
K⋅\dot{z}
%% Cell type:code id: tags:
``` python
dot_f = f_tau.diff(tau) * dot_tau_ + f_tau.diff(Z) * dot_Z_
dot_f
```
%% Output
$\displaystyle \frac{E_\mathrm{b} \left(\dot{s} - \dot{s}_\mathrm{pl}\right) \sqrt{\tau^{2}}}{\tau} - K \dot{z}$
_______
╱ 2
E_\mathrm{b}⋅(\dot{s} - \dot{s}_\mathrm{pl})⋅╲╱ \tau
─────────────────────────────────────────────────────── - K⋅\dot{z}
\tau
%% Cell type:code id: tags:
``` python
dot_f_lambda = dot_f.subs(dot_s_pl, dot_s_pl_).subs(dot_z, dot_z_)
dot_f_lambda
```
%% Output
$\displaystyle \frac{E_\mathrm{b} \left(\dot{s} - \frac{\lambda \sqrt{\tau^{2}}}{\tau}\right) \sqrt{\tau^{2}}}{\tau} - K \lambda$
⎛ _______⎞
⎜ ╱ 2 ⎟ _______
⎜ \lambda⋅╲╱ \tau ⎟ ╱ 2
E_\mathrm{b}⋅⎜\dot{s} - ──────────────────⎟⋅╲╱ \tau
⎝ \tau ⎠
────────────────────────────────────────────────────── - K⋅\lambda
\tau
%% Cell type:code id: tags:
``` python
lambda_solved = sp.solve(dot_f_lambda, lambda_)[0]
```
%% Cell type:code id: tags:
``` python
sp.simplify(dot_tau_.subs(dot_s_pl, dot_s_pl_).subs(lambda_, lambda_solved))
```
%% Output
$\displaystyle \frac{E_\mathrm{b} K \dot{s}}{E_\mathrm{b} + K}$
E_\mathrm{b}⋅K⋅\dot{s}
──────────────────────
E_\mathrm{b} + K
%% Cell type:markdown id: tags:
## Construct the bond slip model
%% Cell type:markdown id: tags:
Given an increment of slip, calculate the corresponding amount of stress
regarding the current state of the material
%% Cell type:markdown id: tags:
Let us now solve this problem numerically
\begin{align}
\end{align}
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
bmcs_course/4_3_BS_EP_SH_IK.ipynb deleted 100644 → 0
+ 0
− 1123
View file @ a37bff6c
%% Cell type:markdown id: tags:
# Bond-slip elasto-plstic model with isotropic and kinematic hardening
This notebook is a work in progress on an abstract and general implementation of time integration algorithm for general damage-plasticity modes. It serves for the development of a package that can be configured by specifying the ingredients of thermodynamically based model
- Vector of state variables $\boldsymbol{\mathcal{E}}$
- Vector of streses $\boldsymbol{\mathcal{S}}$
- Yield condition $f(\boldsymbol{\mathcal{S}},\boldsymbol{\mathcal{E}})$
- Flow potential $\varphi(\boldsymbol{\mathcal{S}},\boldsymbol{\mathcal{E}})$
as symbolic equations using the sympy package. The time-stepping algorithm gets generated automatically within the thermodynamically framework. The derived evolution equations and return-mapping to the yield surface is performed using Newton-Raphson scheme.
%% Cell type:code id: tags:
``` python
%matplotlib notebook
import sympy as sp
sp.init_printing()
import matplotlib.pyplot as plt
import numpy as np
```
%% Cell type:markdown id: tags:
## Material parameters
%% Cell type:code id: tags:
``` python
E_b = sp.Symbol('E_b', real=True, nonnegative=True)
gamma = sp.Symbol('gamma', real=True, nonnegative=True)
K = sp.Symbol('K', real=True)
tau_bar = sp.Symbol(r'\bar{\tau}', real=True, nonnegative=True)
```
%% Cell type:code id: tags:
``` python
py_vars = ('E_b', 'gamma', 'K', 'tau_bar')
map_py2sp = {py_var : globals()[py_var] for py_var in py_vars}
sp_vars = tuple(map_py2sp[py_var] for py_var in py_vars)
```
%% Cell type:markdown id: tags:
## State variables
%% Cell type:code id: tags:
``` python
s = sp.Symbol('s', real=True)
s_pi = sp.Symbol(r's_pi', real=True)
alpha = sp.Symbol('alpha', real=True)
z = sp.Symbol('z', real=True)
```
%% Cell type:code id: tags:
``` python
Eps = sp.Matrix([s_pi, z, alpha])
Eps.T
```
%% Cell type:markdown id: tags:
## Thermodynamic forces
%% Cell type:code id: tags:
``` python
tau = sp.Symbol('tau', real=True)
X = sp.Symbol('X', real=True)
Z = sp.Symbol('Z', real=True, nonnegative=True)
```
%% Cell type:code id: tags:
``` python
Sig = sp.Matrix([tau, Z, X])
Sig.T
```
%% Cell type:code id: tags:
``` python
Sig_Eps = sp.Matrix([[E_b * (s - s_pi), K * z, alpha * gamma]]).T
Sig_Eps
```
%% Cell type:markdown id: tags:
**Executable code for** $\boldsymbol{\mathcal{S}}(s,\boldsymbol{\mathcal{E}})$
%% Cell type:code id: tags:
``` python
get_Sig = sp.lambdify(
(s, Eps) + sp_vars, Sig_Eps.T, 'numpy'
)
```
%% Cell type:markdown id: tags:
To derive the time stepping procedure we will need also the matrix of derivatives of the generalized stresses $\boldsymbol{\mathcal{S}}$ with respect to the kinematic variables $\boldsymbol{\mathcal{E}}$
\begin{align}
\partial_\boldsymbol{\mathcal{E}} \boldsymbol{\mathcal{S}}
\end{align}
%% Cell type:code id: tags:
``` python
dSig_dEps = sp.Matrix([Sig_Eps.T.diff(eps) for eps in Eps]).T
dSig_dEps
```
%% Cell type:markdown id: tags:
**Executable Python code generation** $\partial_\boldsymbol{\mathcal{E}} \boldsymbol{\mathcal{S}}(s,\boldsymbol{\mathcal{E}})$
%% Cell type:code id: tags:
``` python
get_dSig_dEps = sp.lambdify(
(s, Eps) + sp_vars, dSig_dEps, 'numpy'
)
```
%% Cell type:markdown id: tags:
## Threshold function
%% Cell type:markdown id: tags:
To keep the framework general for different stress norms and hardening definitions let us first introduce a general function for effective stress. Note that the observable stress $\tau$ is identical with the plastic stress $\tau_\pi$ due to the performed sign switch in the definition of the thermodynamic forces.
%% Cell type:code id: tags:
``` python
f_Sig = sp.sqrt((tau - X)*(tau - X)) - (tau_bar + Z)
f_Sig
```
%% Cell type:markdown id: tags:
**Executable code generation** $f(\boldsymbol{\mathcal{E}}, \boldsymbol{\mathcal{S}})$
Note that this is a function of both the forces and kinematic state variables
%% Cell type:code id: tags:
``` python
get_f_Sig = sp.lambdify(
(Eps, Sig) + sp_vars, f_Sig, 'numpy'
)
```
%% Cell type:markdown id: tags:
The derivative of $f$ required for time-stepping $\partial_\boldsymbol{\mathcal{S}} f$ is obtained as
%% Cell type:code id: tags:
``` python
df_dSig = f_Sig.diff(Sig)
sp.simplify(df_dSig).T
```
%% Cell type:markdown id: tags:
**Executable code generation** $\partial_\boldsymbol{\mathcal{S}}f(\boldsymbol{\mathcal{E}}, \boldsymbol{\mathcal{S})}$
%% Cell type:code id: tags:
``` python
get_df_dSig = sp.lambdify(
(Eps, Sig) + sp_vars, df_dSig, 'numpy'
)
```
%% Cell type:markdown id: tags:
## Flow direction
Evolution equations have the form
\begin{align}
\dot{\boldsymbol{\mathcal{E}}} = \lambda \, \boldsymbol{\Phi}
\end{align}
with the vector $\Phi$ representing the flow direction within
the stress space. Assuming the normality condition, i.e. that
the flow direction coincides with the vector normal to the
yield condition we can write
\begin{align}
\boldsymbol{\Phi} = \boldsymbol{\Upsilon} \frac{\partial f}{\partial \boldsymbol{\mathcal{S}}}
\end{align}
The sign matrix $\boldsymbol{\Upsilon}$ is used to direct all the derivatives in the outward direction from the elastic range.
%% Cell type:code id: tags:
``` python
Sig_signs = sp.diag(-1,1,1)
Phi = -Sig_signs * df_dSig
Phi
```
%% Cell type:code id: tags:
``` python
get_Phi = sp.lambdify(
(Eps, Sig) + sp_vars, Phi, 'numpy'
)
```
%% Cell type:markdown id: tags:
# Time integration scheme
%% Cell type:markdown id: tags:
## Summary of the backward Euler scheme
The derived expressions can be now plugged-in into a generic return mapping algorithm that efficiently identifies a state that satisfies the discrete consistency condition. The general structure of an implicit integration scheme reads
\begin{align}
\boldsymbol{\mathcal{E}}_{n+1} &= \boldsymbol{\mathcal{E}}_{n} +
\lambda_\Delta \, \boldsymbol{\Phi}_{n+1} \\
f(\boldsymbol{\mathcal{E}}_{n+1}; \lambda_\Delta) &= 0
\end{align}
To reach an admissible state let us linearize the threshold function at an intermediate state $k$ as
\begin{align}
f(\boldsymbol{\mathcal{E}}^{(k)}; \lambda_\Delta^{(k)} )
& \approx
f^{(k)}
+
\left. \frac{\partial f}{\partial \lambda} \right|^{(k)}
\Delta \lambda
\end{align}
%% Cell type:markdown id: tags:
Thus, by rewriting the linearized equation as a recurrence formula, the iteration algorithm is obtained
\begin{align}
&\left. \frac{\partial f}{\partial \lambda} \right|^{(k)}
\Delta \lambda =
- f ^{(k)}
\\
&\lambda_{\Delta}^{(k+1)} = \lambda_{\Delta}^{(k)} + \Delta \lambda \\
& \boldsymbol{\mathcal{E}}^{(k+1)} =
\boldsymbol{\mathcal{E}}^{(k)} +
\lambda_\Delta \, \boldsymbol{\Phi}^{(k)}
\\
&k = k + 1
\end{align}
%% Cell type:markdown id: tags:
To define a generic return mapping we need to construct the derivatives of the flow rule $f$ with respect to $\lambda$. The dependency of $f$ on $\lambda$ is intermediated via the stresses $\boldsymbol{\mathcal{S}}$ and state variables $\boldsymbol{\mathcal{E}}$
\begin{align}
f(\boldsymbol{\mathcal{S}}(\boldsymbol{\mathcal{E}}(\lambda))).
\end{align}
To correctly resolve the dependencies in the derivative $\frac{\partial f}{\partial_\lambda}$, we need to apply rules for chaining of derivatives. Let us start with the derivative with respect to $\boldsymbol{\mathcal{E}}$ in the form
\begin{align}
\frac{\partial f(\boldsymbol{\mathcal{S}}(\boldsymbol{\mathcal{E}}))}{\partial \boldsymbol{\mathcal{E}}}
&=
\frac{\partial f(\boldsymbol{\mathcal{S}})}{\partial \boldsymbol{\mathcal{S}}} \,
\frac{\partial \boldsymbol{\mathcal{S}}(\boldsymbol{\mathcal{E}})}{\partial \boldsymbol{\mathcal{E}}}.
\end{align}
By expanding the derivatives of $\boldsymbol{\mathcal{E}}$ with respect to $\lambda_\Delta$ that will be abbreviate in index position as $\lambda$ for brevity we obtain
\begin{align}
\frac{\partial f(\boldsymbol{\mathcal{S}}(\boldsymbol{\mathcal{E}}(\lambda)))}{\partial \lambda}
&=
\frac{\partial f}{\partial \boldsymbol{\mathcal{S}}} \,
\frac{\partial \boldsymbol{\mathcal{S}}}{\partial \boldsymbol{\mathcal{E}}}
\frac{\partial \boldsymbol{\mathcal{E}}}{\partial \lambda}.
\end{align}
%% Cell type:markdown id: tags:
The last term $\frac{\partial \boldsymbol{\mathcal{E}} }{ \partial \lambda}$ can be obtained from the evolution equations
\begin{align}
\boldsymbol{\mathcal{E}} = \lambda \, \boldsymbol{\Phi} \; \implies
\frac{\partial \boldsymbol{\mathcal{E}} }{\partial \lambda} =
\boldsymbol{\Phi}
\end{align}
%% Cell type:markdown id: tags:
**Summarizing**: the algorithm can be written in a compact way as follows:
\begin{align}
&
\left(
\frac{\partial f}{\partial \boldsymbol{\mathcal{S}}}
^{(k)} \,
\frac{\partial \boldsymbol{\mathcal{S}}}{\partial \boldsymbol{\mathcal{E}}}
^{(k)} \,
\boldsymbol{\Phi}^{(k)}
\right)
\Delta \lambda = -
f^{(k)}\\
&\lambda_{\Delta}^{(k+1)} = \lambda_{\Delta}^{(k)} + \Delta \lambda \\
& \boldsymbol{\mathcal{E}}^{(k+1)} = \boldsymbol{\mathcal{E}}^{(k)} +
\lambda_\Delta \, \boldsymbol{\Phi}^{(k)}
\\
&k = k + 1
\end{align}
%% Cell type:markdown id: tags:
## Implementation concept
The gradient operators needed for the time-stepping scheme have been derived above and are now available for the implementation of the numerical algorithm both in `Python` and `C89` languages
<table style="width:50%">
<tr>
<th>Symbol</th>
<th>Python</th>
<th>C89</th>
</tr>
<tr>
<td>$\mathcal{S}(\boldsymbol{\mathcal{E}}) $
</td>
<td>get_Sig</td>
<td>get_Sig_C</td>
</tr>
<tr>
<td>$ f(\boldsymbol{\mathcal{S}}, \boldsymbol{\mathcal{E}})$</td>
<td>get_f</td>
<td>get_f_C</td>
</tr>
<tr>
<td>
$\displaystyle{\frac{\partial f}{\partial \boldsymbol{\mathcal{S}}}}(\boldsymbol{\mathcal{S}}, \boldsymbol{\mathcal{E}})$
</td>
<td>get_df_dSig</td>
<td>get_df_dSig_C</td>
</tr>
<tr>
<td>
$\displaystyle{\frac{\partial \boldsymbol{\mathcal{S}}}{\partial \boldsymbol{\mathcal{E}}}}(\boldsymbol{\mathcal{E}})$</td>
<td>get_dSig_dEps</td>
<td>get_dSig_dEps_C</td>
</tr>
<tr>
<td>$\boldsymbol{\Phi}(\boldsymbol{\mathcal{S}}, \boldsymbol{\mathcal{E}}) $</td>
<td>get_Phi</td>
<td>get_Phi_C</td>
</tr>
</table>
%% Cell type:markdown id: tags:
To avoid repeated calculation of the same expressions, let us put the evaluation of $f$ and $\frac{\partial f}{\partial \lambda}$ into a single procedure. Indeed, the iteration loop can be constructed in such a way that the predictor $\frac{\partial f}{\partial \lambda}$ for the next step is calculated along with the residuum $f$. In case that the residuum is below the required tolerance, the overhead for an extra calculated derivative is negligible or, with some care, it can be even reused in the next time step.
%% Cell type:code id: tags:
``` python
def get_f_df(s_n1, Eps_k, *margs):
Sig_k = get_Sig(s_n1, Eps_k, *margs)[0]
f_k = np.array([get_f_Sig(Eps_k, Sig_k, *margs)])
df_dSig_k = get_df_dSig(Eps_k, Sig_k, *margs)
Phi_k = get_Phi(Eps_k, Sig_k, *margs)
dSig_dEps_k = get_dSig_dEps(s_n1, Eps_k, *margs)
df_dSigEps_k = np.einsum(
'ik,ji->jk', df_dSig_k, dSig_dEps_k)
dEps_dlambda_k = Phi_k
df_dlambda = np.einsum(
'ki,kj->ij', df_dSigEps_k, dEps_dlambda_k)
df_k = df_dlambda
return f_k, df_k, Sig_k
```
%% Cell type:markdown id: tags:
The update of state variables $\boldsymbol{\mathcal{E}}^{(k+1)}$ for an newly obtained $\lambda_\Delta^{(k+1)}$ is performed using the evolution equations.
%% Cell type:code id: tags:
``` python
def get_Eps_k1(s_n1, Eps_n, lam_k, Eps_k, *margs):
Sig_k = get_Sig(s_n1, Eps_k, *margs)[0]
Phi_k = get_Phi(Eps_k, Sig_k, *margs)
Eps_k1 = Eps_n + lam_k * Phi_k[:,0]
return Eps_k1
```
%% Cell type:markdown id: tags:
The double loop over the time increments and over the return mapping iteration. The inner loop represents the material point level in a standard finite element calculation. The input is the maximum slip value, the number of time steps, the maximum number of iterations and a load function which can define cyclic loading as shown below. The procedure returns the record of $\boldsymbol{\mathcal{E}}(t)$ and $\boldsymbol{\mathcal{S}}(t)$
%% Cell type:code id: tags:
``` python
def get_response(margs, s_max=3, n_steps = 10, k_max=20, get_load_fn=lambda t: t):
Eps_n = np.zeros((len(Eps),), dtype=np.float_)
Eps_k = np.copy(Eps_n)
Sig_record, Eps_record, iter_record = [], [], []
t_arr = np.linspace(0,1,n_steps+1)
s_t = s_max * get_load_fn(t_arr) + 1e-9
for s_n1 in s_t:
lam_k = 0
f_k, df_k, Sig_k = get_f_df(s_n1, Eps_k, *margs)
f_k_norm = np.linalg.norm(f_k)
f_k_trial = f_k[-1]
k = 0
while k < k_max:
if f_k_trial < 0 or f_k_norm < 1e-8:
Eps_n[...] = Eps_k[...]
Sig_record.append(Sig_k)
Eps_record.append(np.copy(Eps_k))
iter_record.append(k+1)
break
dlam = np.linalg.solve(df_k, -f_k)
lam_k += dlam
Eps_k = get_Eps_k1(s_n1, Eps_n, lam_k, Eps_k, *margs)
f_k, df_k, Sig_k = get_f_df(s_n1, Eps_k, *margs)
f_k_norm = np.linalg.norm(f_k)
k += 1
else:
print('no convergence')
Sig_arr = np.array(Sig_record, dtype=np.float_)
Eps_arr = np.array(Eps_record, dtype=np.float_)
iter_arr = np.array(iter_record,dtype=np.int_)
return t_arr, s_t, Eps_arr, Sig_arr, iter_arr
```
%% Cell type:markdown id: tags:
# Support functions
To run some examples, let us define some infrastructure including a more complex loading history and postprocessing
%% Cell type:markdown id: tags:
## Loading history
This implementation uses the symbolic machinery which is not necessary a simpler data point based implementation with `numpy.interp1d` would be better ... later
%% Cell type:code id: tags:
``` python
t, theta = sp.symbols(r't, \theta')
n_cycles = 5
A = 2
ups = np.array([((theta-2*cycle)*A+(1-A), theta-2*cycle<=1)
for cycle in range(n_cycles)])
downs = np.array([((1-(theta-(2*cycle+1)))*A+(1-A),(theta-(2*cycle+1))<=1)
for cycle in range(n_cycles)])
ups[0,0] = theta
updowns = np.einsum('ijk->jik',np.array([ups, downs])).reshape(-1,2)
load_fn = sp.Piecewise(*updowns).subs(theta,t*n_cycles)
get_load_fn = sp.lambdify(t, load_fn,'numpy')
t_arr = np.linspace(0,1,600)
plt.plot(t_arr, get_load_fn(t_arr));
```
%% Cell type:markdown id: tags:
## Plotting functions
To simplify postprocessing examples, here are two aggregate plotting functions, one for the state and force variables, the other one for the evaluation of energies
%% Cell type:code id: tags:
``` python
def plot_Sig_Eps(t_arr, s_t, Sig_arr, Eps_arr, iter_arr, ax1, ax11, ax2, ax22, ax3, ax33):
colors = ['blue','red', 'green', 'black', 'magenta' ]
s_pi_, z_, alpha_ = Eps_arr.T
sig_pi_, Z_, X_ = Sig_arr.T
n_step = len(s_pi_)
ax1.plot(s_t, sig_pi_, color='black',
label='n_steps = %g' % n_step)
ax1.set_xlabel('$s$'); ax1.set_ylabel(r'$\tau$')
ax1.legend()
if ax11:
ax11.plot(s_t, iter_arr, '-.')
ax2.plot(t_arr, z_, color='green',
label='n_steps = %g' % n_step)
ax2.set_xlabel('$t$'); ax2.set_ylabel(r'$z$')
if ax22:
ax22.plot(t_arr, Z_, '-.', color='green')
ax22.set_ylabel(r'$Z$')
ax3.plot(t_arr, alpha_, color='blue',
label='n_steps = %g' % n_step)
ax3.set_xlabel('$t$'); ax3.set_ylabel(r'$\alpha$')
if ax33:
ax33.plot(t_arr, X_, '-.', color='blue')
ax33.set_ylabel(r'$X$')
```
%% Cell type:code id: tags:
``` python
from scipy.integrate import cumtrapz
def plot_work(ax, t_arr, s_t, Eps_arr, Sig_arr):
W_arr = cumtrapz(Sig_arr[:,0], s_t, initial=0)
U_arr = Sig_arr[:,0] * (s_t-Eps_arr[:,0]) / 2.0
G_arr = W_arr - U_arr
ax.plot(t_arr, W_arr, lw=2, color='black', label=r'$W$')
ax.plot(t_arr, G_arr, color='black', label=r'$G$')
ax.fill_between(t_arr, W_arr, G_arr, color='green', alpha=0.2)
ax.set_xlabel('$s$'); ax3.set_ylabel(r'$E$')
ax.legend()
```
%% Cell type:code id: tags:
``` python
def plot_dissipation(ax, t_arr, s_t, Eps_arr, Sig_arr):
colors = ['blue','red', 'green', 'black', 'magenta' ]
E_i = cumtrapz(Sig_arr, Eps_arr, initial=0, axis=0)
c = 'black'
ax.plot(t_arr, E_i[:,0], '-.', lw=1, color=c)
ax.fill_between(t_arr, E_i[:,0], 0, color=c, alpha=0.1)
c = 'black'
ax.plot(t_arr, E_i[:,0], color=c, lw=1)
ax.fill_between(t_arr, E_i[:,0], E_i[:,0],
color=c, alpha=0.2);
c = 'blue'
ax.plot(t_arr, E_i[:,1], '-.', lw=1, color='black')
ax.fill_between(t_arr, E_i[:,1], 0, color=c, alpha=0.1)
c = 'blue'
ax.plot(t_arr, E_i[:,1] + E_i[:,2], color='black', lw=1)
ax.fill_between(t_arr, E_i[:,1] + E_i[:,2], E_i[:,1],
color=c, alpha=0.3);
```
%% Cell type:markdown id: tags:
# Examples
%% Cell type:code id: tags:
``` python
material_params = {
E_b:1, gamma: 0.0, K:0.1, tau_bar:1,
}
margs = [material_params[map_py2sp[name]] for name in py_vars]
```
%% Cell type:markdown id: tags:
## Monotonic load
Let's first run the example with different size of the time step to see if there is any difference
%% Cell type:code id: tags:
``` python
fig, ((ax1,ax2,ax3)) = plt.subplots(1,3,figsize=(14,4), tight_layout=True)
ax11, ax22, ax33 = ax1.twinx(), ax2.twinx(), ax3.twinx()
for n_steps in [20, 40, 200, 2000]:
t_arr, s_t, Eps_arr, Sig_arr, iter_arr = get_response(
margs=margs, s_max=8, n_steps=n_steps, k_max=10
)
plot_Sig_Eps(t_arr, s_t, Sig_arr, Eps_arr, iter_arr, ax1, ax11, ax2, ax22, ax3, ax33)
```
%% Cell type:code id: tags:
``` python
fig, ax = plt.subplots(1,1,figsize=(9, 5))
plot_work(ax, t_arr, s_t, Eps_arr, Sig_arr)
plot_dissipation(ax, t_arr, s_t, Eps_arr, Sig_arr)
```
%% Cell type:markdown id: tags:
## Cyclic loading
%% Cell type:code id: tags:
``` python
fig, ((ax1,ax2,ax3)) = plt.subplots(1,3,figsize=(14,4), tight_layout=True)
ax11,ax22,ax33 = ax1.twinx(), ax2.twinx(), ax3.twinx()
t_arr, s_t, Eps_arr, Sig_arr, iter_arr = get_response(
margs, s_max=2, n_steps=20000, k_max=20, get_load_fn=get_load_fn
)
plot_Sig_Eps(t_arr, s_t, Sig_arr, Eps_arr, iter_arr, ax1, ax11, ax2, ax22, ax3, ax33);
```
%% Cell type:code id: tags:
``` python
fig, ax = plt.subplots(1,1,figsize=(9, 5))
plot_work(ax, t_arr, s_t, Eps_arr, Sig_arr)
plot_dissipation(ax, t_arr, s_t, Eps_arr, Sig_arr)
```
%% Cell type:markdown id: tags:
## Interactive application
%% Cell type:code id: tags:
``` python
def init():
global Eps_record, Sig_record, iter_record, t_arr, s_t, s0, t0, Eps_n
s0 = 0
t0 = 0
Sig_record = []
Eps_record = []
iter_record = []
t_arr = []
s_t = []
Eps_n = np.zeros((len(Eps),), dtype=np.float_)
def get_response_i(s1, margs, n_steps = 60, k_max=20):
global Eps_record, Sig_record, iter_record, t_arr, s_t, s0, t0, Eps_n
Eps_k = np.copy(Eps_n)
t1 = t0+n_steps+1
ti_arr = np.linspace(t0, t1, n_steps+1 )
si_t = np.linspace(s0,s1,n_steps+1)
for s_n1 in si_t:
lam_k = 0
f_k, df_k, Sig_k = get_f_df(s_n1, Eps_k, *margs)
f_k_norm = np.linalg.norm(f_k)
f_k_trial = f_k[-1]
k = 0
while k < k_max:
if f_k_trial < 0 or f_k_norm < 1e-6:
Eps_n[...] = Eps_k[...]
Sig_record.append(Sig_k)
Eps_record.append(np.copy(Eps_k))
iter_record.append(k+1)
break
dlam = np.linalg.solve(df_k, -f_k)
lam_k += dlam
Eps_k = get_Eps_k1(s_n1, Eps_n, lam_k, Eps_k, *margs)
f_k, df_k, Sig_k = get_f_df(s_n1, Eps_k, *margs)
f_k_norm = np.linalg.norm(f_k)
k += 1
else:
print('no convergence')
t_arr = np.hstack([t_arr, ti_arr])
s_t = np.hstack([s_t, si_t])
t0 = t1
s0 = s1
return
import ipywidgets as ipw
fig, ((ax1,ax2,ax3)) = plt.subplots(1,3,figsize=(14,4), tight_layout=True)
ax11 = ax1.twinx()
ax22 = ax2.twinx()
ax33 = ax3.twinx()
axes = ax1, ax11, ax2, ax22, ax3, ax33
axes = ax1, None, ax2, ax22, ax3, ax33
def update(s1):
global Eps_record, Sig_record, iter_record, t_arr, s_t, s0, t0, Eps_n, axes
global margs
get_response_i(s1, margs)
Sig_arr = np.array(Sig_record, dtype=np.float_)
Eps_arr = np.array(Eps_record, dtype=np.float_)
iter_arr = np.array(iter_record,dtype=np.int_)
for ax in axes:
if ax:
ax.clear()
plot_Sig_Eps(t_arr, s_t, Sig_arr, Eps_arr, iter_arr, *axes)
init()
s1_slider = ipw.FloatSlider(value=0,min=-4, max=+4, step=0.1,
continuous_update=False)
ipw.interact(update, s1 = s1_slider);
def reset(**material_params):
global margs
init()
s1_slider.value = 0
margs = [material_params[name] for name in py_vars]
n_steps = 20
margs_sliders = {
name : ipw.FloatSlider(description=name, value=val,
min=minval, max=maxval, step=(maxval-minval) / n_steps,
continuous_update=False)
for name, val, minval, maxval in [('E_b', 50, 0.5, 100),
('gamma', 1, -20, 20),
('K', 1, -20, 20),
('tau_bar', 1, 0.5, 20)]
}
ipw.interact(reset, **margs_sliders);
```
bmcs_course/4_3_BS_EP_SH_IK_N.ipynb 0 → 100644
+ 6015
− 0
View file @ a15f314a
Source diff could not be displayed: it is too large. Options to address this: view the blob.
index.ipynb
+ 6
− 5
View file @ a15f314a
%% Cell type:markdown id: tags:
# BMCS applications
Collection of notebooks accompanying the course on brittle-matrix composite structures.
@author: rosoba
## Rule of mixtures for elastic composites
- [J0101 - Elastic mixture rule](bmcs_course/1_1_elastic_stiffness_of_the_composite.ipynb)
# BOND
## Constant bond-slip law
- [J0201 - Pull-out of long fiber from rigid matrix](bmcs_course/2_1_PO_LF_LM_RG.ipynb)
- [J0202 - Pull-out of long fiber from long elastic matrix](bmcs_course/2_2_PO_LF_LM_EL.ipynb)
- [J0203 - Pull-out of short fiber from rigid matrix](bmcs_course/2_3_PO_SF_M_RG.ipynb)
- [J0204 - Comparison of several models](bmcs_course/2_4_PO_comparison.ipynb)
## Nonlinear bond-slip law
- [J0301 - Pull-out with softening and hardening](bmcs_course/3_1_PO_LF_LM_EL_FE_CB.ipynb)
- [J0302 - EXTRA - Newton iterative scheme](extras/newton_method.ipynb)
- [J0303 - EXTRA - Nonlinear finite-element solver for 1d pullout](extras/pullout1d.ipynb)
## Unloading, reloading and inelasticity
- [J0401 - Unloading with multi-linear bond-slip law](bmcs_course/4_1_PO_multilinear_unloading.ipynb) (PO_BS_ML)
- [J0402 - Plasticity framework](bmcs_course/4_2_BS_EP_SH_IK_analytical.ipynb) (BS_EP_SH_IK_A)
- [J0403 - Iterative numerical simulation](bmcs_course/4_3_BS_EP_SH_IK_numerical.ipynb) (BS_EP_SH_IK_N)
- [J0402 - Plasticity framework](bmcs_course/4_2_BS_EP_SH_IK_A.ipynb) (BS_EP_SH_IK_A)
- [J0403 - Iterative numerical simulation](bmcs_course/4_3_BS_EP_SH_IK_N.ipynb) (BS_EP_SH_IK_N)
- [J0404 - Damage framework](bmcs_course)
- [J0405 - Iterative numerical simulation](bmcs_course)
- [J0405 - Iterative numerical simulation](bmcs_course)/4_3_BS_EP_SH_IK_N
# CRACK
## Energy dissipation
# CRACK
%% Cell type:code id: tags:
``` python
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment