This page was generated from unit-6.3-plasticity/plasticity.ipynb.
6.3 Elasto-Plasticity¶
This tutorial considers Hencky-type [1] plasticity with isotropic hardening in 2d (plane strain). It showcases the classes IntegrationRuleSpace, NewtonCF and MinimizationCF.
We begin with a classical formulation with an explicit yield surface that is probably more familiar to engineers. In a second stage, we consider a (constrained) minimization formulation which leads to a more streamlined implementation. Both formulations correspond to the isotropic hardening model discussed in [2].
The essential model ingredients are given in terms of the dissipation potential \(\Phi\) and the stored energy density contribution \(\Psi\), see eg, [3]. For convenience, we first introduce the strain energy density \(\Psi^\text{e}\)
\begin{align*} \Psi^\text{e}(\epsilon) = \frac{1}{2}\| \epsilon\|^2_{M} = \frac{E}{2(1+\nu)(1-2\nu)}\left((1-2\nu)\epsilon:\epsilon + \nu \,\text{tr}(\epsilon)^2\right) = \frac{1}{2} \epsilon : \mathbf{C} : \epsilon. \end{align*}
Moreover, we consider the full strain \(\epsilon\) to be the sum of the elastic \(\epsilon^e\) and the plastic strain contribution \(p\). Conversely, the elastic strain is the given as \begin{equation} \epsilon^\text{e} = \epsilon - p ,\end{equation} which is the strain quantity that will be actually employed as argument to \(\Psi^\text{e}\).
The isotropic hardening variable will be denoted by \(\alpha\) throughout this tutorial. The hardening potential density is simply chosen as \(\Psi^\text{h} = \frac{1}{2}\alpha^2\) such that the total energy density reads \begin{equation} \Psi(\epsilon, p, \alpha) = \Psi^\text{e}(\epsilon-p) + \Psi^\text{h}(\alpha) = (\epsilon - p) : \mathbf{C} : (\epsilon - p) + \frac{1}{2} \alpha^2. \end{equation}
The yield function \(f\) that represents the elastic limit is given as
\begin{equation} f(\sigma, \beta) = \left|\mathrm{dev} \sigma\right| - \sqrt{2/3}\sigma_Y(1 - H \beta), \end{equation}
where we note that the hardening parameter \(H\) enters in \(f\) rather than in \(\Psi\) for consistency with [2]. Deviating from [2], we scale the yield stress \(\sigma_Y\) by \(\sqrt{2/3}\) to obtain the von-Mises yield surface radius. We remark that by virtue of \(\mathrm{dev}\sigma\) appearing in \(f\), the plastic strain (rate) \(p\) (\(\dot p\)) will be deviatoric.
The boundary value problem solved in this tutorial is in plain strains. However, this does not mean that the plastic strains are plain and neither are the elastic ones. Only their sum has vanishing out-of-plane components. Therefore, for simplicity, the material model is implemented in 3d whereas displacements are in 2d.
Formulation with explicit yield surface¶
Based on the principle of maximum dissipation we can abstractly formulate the dissipation potential as
:nbsphinx-math:`begin{align*} Phi(dot{p}, dot{alpha}) =
- sup_{sigma} sup_{beta} inf_{Lambda geq 0}
: sigma : dot{p} + beta dot{alpha} - Lambda f(sigma, beta)
end{align*}`
The stationary conditions of the above functional read as
- :nbsphinx-math:`begin{align*}
&dot{p} = Lambda frac{partial left|mathrm{dev} sigmaright|}{partial sigma} \ &dot{alpha} = Lambda,sqrt{2/3}sigma_YH \ &{Lambda geq 0 : land : (left|mathrm{dev} sigmaright| - sqrt{2/3}sigma_Y(1 - Hbeta)) leq 0 : land : Lambda(left|mathrm{dev} sigmaright| - sqrt{2/3}sigma_Y(1 - Hbeta)) = 0}
end{align*}`
which are regarded as evolution equations for \(p\) and \(\alpha\). The actual values of \(\sigma\) and \(\beta\) are obtained through the constitutive relations provided by the framework of generalized standard materials, we have
\begin{align*} &\frac{\partial\Psi((\epsilon - p), \alpha)}{\partial p} + \frac{\partial\Phi}{\partial \dot{p}} = 0 \quad \Rightarrow \quad \sigma = \mathbf{C} : (\epsilon - p) \\ &\frac{\partial\Psi((\epsilon - p), \alpha)}{\partial \alpha} + \frac{\partial\Phi}{\partial \dot{\alpha}} = 0 \quad \Rightarrow \quad \beta = - \alpha . \end{align*}
The time-discrete form is simply obtained by replacing the time rates with the approximations \(\dot p \approx \frac{p_{k+1} - p_k}{\Delta t}\) and \(\dot \alpha \approx \frac{\alpha_{k+1} - \alpha_k}{\Delta t}\), whereby the value of \(\Delta t\) has no effect for the present rate-independent model.
Constrained minimization formulation¶
where \(\Psi\) denotes the augmented stored energy density and \(\|\bullet \|_{\varepsilon}\) indicates a perturbed norm in order to avoid divisions by zero in the evaluation of the time-discrete evolution equation derived from an incremental variational principle.
Notes in the spatial discretization¶
The state \(u\) is spatially discretized by continuous Lagrange elements, the internal states \(p\) and \(\alpha\) as well as the multiplier \(\Lambda\) reside at quadrature points, for which NGSolve has the notion of an integration rule space. Thus, the problem above can be decomposed into a “global” problem for the displacements and local problems for the internal variables. The latter can be solved individually for each quadrature point. However, there is a coupling between
the local problems and displacements, which has to be accounted for by what is called “algorithmically consistent linearization”. The latter is automatically obtained by using an appropriate BilinearForm object representing the boundary value problem and the internal evolution problems as demonstrated for both implementations.
References 1. W. Han and B.D. Reddy: Mathematical Theory and Numerical Analysis, Springer 1999 2. C. Carstensen, Domain Decomposition for a Non-smooth Convex Minimization Problem and its Application in Plasticity, 1997 3. K. Hackl, F.D. Fischer, On the relation between the principle of maximum dissipation and inelastic evolution given by dissipation potentials, PRSA, 2008
Implementation¶
General setup¶
[1]:
import numpy as np
from ngsolve import *
from ngsolve.comp import IntegrationRuleSpace
from ngsolve.fem import MinimizationCF, NewtonCF
from netgen.geom2d import CSG2d, Circle, Rectangle
from ngsolve.webgui import Draw
SetNumThreads(4)
The problem domain is the classical plate with hole. Thanks to symmetry of the problem to be solved, we actually use only one quarter of the domain together with additional symmetry boundary conditions. Parameters for the geometry and the material employed are taken from [4].
[2]:
# polynomial order for geometry and displacements
order = 3
geo = CSG2d()
circle = Circle(center=(100,100), radius=10.0, bc="curve").Maxh(1)
rect = Rectangle(pmin=(0,100), pmax=(100,200), bottom="bottom", left="left", top="top", right="right")
geo.Add(rect-circle)
mesh = Mesh(geo.GenerateMesh(maxh=5))
mesh.Curve(order)
Draw(mesh)
# points of interest for post-processing
node_A = mesh(100,200)
node_B = mesh(0,200)
Strain energy density and generic helper functions
Model parameters
[3]:
# Shear modulus
mu = 80193.8
# Bulk modulus
kappa = 164206
# Young's modulus
Emod = Parameter((9 * kappa * mu) / (3 * kappa + mu))
# Poisson's ratio
nu = Parameter((3 * kappa - 2 * mu)/(2 * (3 * kappa + mu)))
print("E =", Emod, "\nnu =", nu)
# Hardening parameter
H = Parameter(10)
# Yield stress
sigma_Y = Parameter(450)
# Perturbation parameter for tensor norms: this value is decisive for the accuracy of the result!
# for higher values of H is must be considerably smaller
norm_pert = Parameter(1e-16)
E = ParameterCF, val = 206900
nu = ParameterCF, val = 0.29
[4]:
def strain(u):
_strain = Sym(Grad(u))
return CoefficientFunction((_strain[0,0], _strain[0,1], 0,
_strain[1,0], _strain[1,1], 0,
0 , 0, 0), dims=(3, 3))
def elastic_strain_energy(eps, Emod, nu):
return 1/2 * Emod / ((1 + nu) * (1 - 2 * nu)) * ((1 - 2 * nu) * InnerProduct(eps, eps) + nu * Trace(eps)**2)
def tensor_norm(a, pert=0):
return sqrt(InnerProduct(a, a) + pert)
def dev(a):
return a - Trace(a) * Id(3) / 3
Common FE spaces¶
We employ a \(H^1\) space for the displacements and a “integration rule space” (IR space) as basis space for all internal variables.
[5]:
fes_u = VectorH1(mesh, order=order, dirichletx="right", dirichlety="bottom")
fes_ir = IntegrationRuleSpace(mesh, order=order-1)
Scalar internal states can be defined based on fes_ir as needed.
The “measure” corresponding to fes_ir:
[6]:
# extract integration rules from IntegrationRuleSpace
irs_dx = dx(intrules=fes_ir.GetIntegrationRules())
Load definitions¶
[7]:
force = 450
# For loadsteps
loadfactor = Parameter(1)
Some data structures for postprocessing¶
[8]:
# for postprocessing
drawfes = MatrixValued(L2(mesh, order=order-1), dim=3, symmetric=True, deviatoric=False)
pd, qd = drawfes.TnT()
pdraw = GridFunction(drawfes)
ad = BilinearForm(InnerProduct(pd, qd)*irs_dx, symmetric=True).Assemble()
invad = ad.mat.Inverse()
drawafes = L2(mesh, order=order-1)
pda, qda = drawafes.TnT()
adraw = GridFunction(drawafes)
ada = BilinearForm(InnerProduct(pda, qda)*irs_dx, symmetric=True).Assemble()
invada = ada.mat.Inverse()
Implementation with explicit yield surface¶
In this approach we employ a generic symmetric-matrix-valued space for \(p\).
[9]:
fes_p = MatrixValued(fes_ir, symmetric=True, deviatoric=False, dim=3)
fes_int = fes_p * fes_ir * fes_ir # p x alpha x Lambda
# We need a coupled space as well.
fes = fes_u * fes_int
u, p, alpha, Lambda = fes.TrialFunction()
# GridFunction for full solution
gfsol = GridFunction(fes)
gfu = gfsol.components[0]
# GridFunction for internal states only
gfint = GridFunction(fes_int)
gfp, gfalpha, gfLambda = gfint.components
# For history states
fes_hist = fes_p * fes_ir
gfhist = GridFunction(fes_hist)
gfp_k, gfalpha_k = gfhist.components
# "Trial" states for plastic evolution (these are not trial functions in FE parlance)
gftrial = GridFunction(fes_hist)
gfsigma_trial, gfbeta_trial = gftrial.components
# Time increment
Delta_t = Parameter(1.0)
# Not really needed for a rate-independent model. Value could be used for scaling the equations though.
# Short-hand for storing history variables. Note that the history value of Lambda is not relevant and thus omitted.
def store_internal():
gfp_k.Interpolate(gfp)
gfalpha_k.Interpolate(gfalpha)
Next, we setup functions for the energy density \(\Psi\) and the dissipation potential \(\Phi\)
[10]:
# Energy density
def Psi(strain, p, alpha):
strain_energy = elastic_strain_energy(strain-p, Emod, nu)
hardening = 1/2 * alpha**2
return strain_energy + hardening
# The yield condition: result > 0 ==> plastic evolution
def yield_condition(sigma, beta, pert=0):
return tensor_norm(dev(sigma), pert=pert) - sqrt(2/3) * sigma_Y * (1 - H*beta)
# The objective function defining the dissipation potential
def Phi_obj(sigma, beta, p_dot, alpha_dot, Lambda):
return InnerProduct(sigma, p_dot) + InnerProduct(beta, alpha_dot) - \
Lambda * yield_condition(sigma, beta, pert=norm_pert)
# The evolution equations
def evolution_eqs(strain, p, alpha, Lambda, evol_tol):
p_dot = (p - gfp_k) / Delta_t
alpha_dot = (alpha - gfalpha_k) / Delta_t
sigma = -Psi(strain, p, alpha).Diff(p)
beta = -Psi(strain, p, alpha).Diff(alpha)
_Phi_obj = Phi_obj(sigma, beta, p_dot, alpha_dot, Lambda)
dsigma = _Phi_obj.Diff(sigma)
# only active when there is hardening
dbeta = IfPos(H - 1e-16, _Phi_obj.Diff(beta), alpha_dot)
dLambda = _Phi_obj.Diff(Lambda)
# Only solve the evolution equations if the elastic trial state is beyond the yield surface.
# This essentially leads to the fulfillment of Lambda >= 0
return IfPos(yield_condition(gfsigma_trial, gfbeta_trial) - evol_tol,
CoefficientFunction((dsigma, dbeta, dLambda)),
CoefficientFunction((p_dot, alpha_dot, Lambda - gfLambda)))
Setup the handling of the evolution of internal variables. The evolution equations are a system of nonlinear equations per quadrature point. NewtonCF solves such equations with given starting points, tolerance and maximum interations. On failure, it returns NaN.
[11]:
# Parameters for evolution
tol = 1e-10
maxiter = 20
evol_tol = Parameter(tol)
eqs = evolution_eqs(strain(gfu), p, alpha, Lambda, evol_tol).Compile(realcompile=False)
evolution = NewtonCF(eqs, gfsol.components[1:], tol=tol, maxiter=maxiter)
Warning: differentiationg by a variable not marked as Variable,
might be optimized out. Call MakeVariable for differentiation CF
Warning: differentiationg by a variable not marked as Variable,
might be optimized out. Call MakeVariable for differentiation CF
Warning: differentiationg by a variable not marked as Variable,
might be optimized out. Call MakeVariable for differentiation CF
Warning: differentiationg by a variable not marked as Variable,
might be optimized out. Call MakeVariable for differentiation CF
Warning: differentiationg by a variable not marked as Variable,
might be optimized out. Call MakeVariable for differentiation CF
The “trial state” that indicates whether there will be a plastic evolution or is computed via an exact (interpolating) projection.
[12]:
# We employ a BilinearForm to evaluate the weak form of the projection problem.
# Thus, we will only call L_trial.Apply
L_trial = BilinearForm(fes_hist, nonassemble=True)
p_trial, alpha_trial = fes_hist.TrialFunction()
p_test, alpha_test = fes_hist.TestFunction()
Psi_trial = Psi(strain(gfu), p_trial, alpha_trial)
L_trial += InnerProduct(-Psi_trial.Diff(p_trial), p_test) * irs_dx
L_trial += InnerProduct(-Psi_trial.Diff(alpha_trial), alpha_test) * irs_dx
# The mass matrix of the projection.
M_trial = BilinearForm(fes_hist, symmetric=True, diagonal=True)
for _trial, _test in zip(*fes_hist.TnT()):
M_trial += InnerProduct(_trial, _test).Compile() * irs_dx
M_trial.Assemble()
M_trial_inv = M_trial.mat.Inverse()
# RHS vector
vec_trial = gftrial.vec.CreateVector()
def compute_trial_state():
L_trial.Apply(gfhist.vec, vec_trial)
gftrial.vec.data = M_trial_inv * vec_trial
Warning: differentiationg by a variable not marked as Variable,
might be optimized out. Call MakeVariable for differentiation CF
Warning: differentiationg by a variable not marked as Variable,
might be optimized out. Call MakeVariable for differentiation CF
The actual solution of the equations happens during an again exact (interpolating) projection that sets gfint.
[13]:
# A Linearform for which the evluation triggers the solution of the evolution equations
testfunc = CoefficientFunction(tuple(fes_int.TestFunction()))
L_evol = LinearForm(fes_int)
L_evol += InnerProduct(CacheCF(evolution), testfunc) * irs_dx
# Mass-Matrix for the internal variable space as generic interpolation tool
M_evol = BilinearForm(fes_int, symmetric=True, diagonal=True)
for _trial, _test in zip(*fes_int.TnT()):
M_evol += InnerProduct(_trial, _test).Compile() * irs_dx
M_evol.Assemble()
M_evol_inv = M_evol.mat.Inverse()
def evolve_internal():
compute_trial_state()
L_evol.Assemble()
gfint.vec.data = M_evol_inv * L_evol.vec
for i, comp in enumerate(gfint.components):
gfsol.components[i+1].vec.data = gfint.components[i].vec
What remains is to setup the BilinearForm that describes the global boundary value problem.
[14]:
testfunc = fes.TestFunction()
u_test = testfunc[0]
test_int = CoefficientFunction(tuple(testfunc[1:]))
# Create a strain variable for partial differentiation
strain_var = strain(u)
internal_virtual_work_density = InnerProduct(Psi(strain_var, p, alpha).Diff(strain_var), strain(u_test))
# This "constraint" ensures the consistent linearization
evolution_constraint = InnerProduct(evolution_eqs(strain_var, p, alpha, Lambda, evol_tol), test_int)
loading_potential = -u[1] * force * loadfactor * ds("top")
realcompile = False
# Possible optimization: use static condensation for p
a = BilinearForm(fes, symmetric=False)
a += internal_virtual_work_density.Compile(realcompile=realcompile) * irs_dx
a += evolution_constraint.Compile(realcompile=realcompile) * irs_dx
a += Variation(loading_potential).Compile(realcompile=realcompile)
Warning: differentiationg by a variable not marked as Variable,
might be optimized out. Call MakeVariable for differentiation CF
Warning: differentiationg by a variable not marked as Variable,
might be optimized out. Call MakeVariable for differentiation CF
Warning: differentiationg by a variable not marked as Variable,
might be optimized out. Call MakeVariable for differentiation CF
Warning: differentiationg by a variable not marked as Variable,
might be optimized out. Call MakeVariable for differentiation CF
Warning: differentiationg by a variable not marked as Variable,
might be optimized out. Call MakeVariable for differentiation CF
Warning: differentiationg by a variable not marked as Variable,
might be optimized out. Call MakeVariable for differentiation CF
The nonlinear PDE is solved with a Newton-Raphson scheme. Thereby, in each Newton iteration we solve a linearized PDE for the displacements and a nonlinear evolution problem at each quadrature point via ``NewtonCF``
In order to run the problem, we need to define vectors that store respective solution data, define loadsteps and setup a loop that iterates through the latter.
[15]:
# A set of functions/vectors for keeping data
w = gfsol.vec.CreateVector()
rhs = gfsol.vec.CreateVector()
# Postprocess internal states
fd = LinearForm(InnerProduct(gfp, qd)*irs_dx)
fda = LinearForm(InnerProduct(gfalpha_k, qda)*irs_dx)
# Load steps (chosen based on experience)
loadsteps = [0.1,0.3,0.5,0.7,0.8,0.9,0.95,1.0]
# Set solution to zero initially
gfsol.vec[:] = 0
gfint.vec[:] = 0
gfhist.vec[:] = 0
gftrial.vec[:] = 0
# Iterate through load steps
for ls in loadsteps:
loadfactor.Set(ls)
# Update old solution at time t = t_k
store_internal()
with TaskManager():
for i in range(20):
a.AssembleLinearization(gfsol.vec)
a.Apply(gfsol.vec, rhs)
# Static condensation could be employed to speedup solving
w.data = a.mat.Inverse(freedofs=fes.FreeDofs(False), inverse="umfpack") * rhs
gfsol.vec.data -= w
# Solve the evolution equastion (quadrature-point-wise)
evolve_internal()
if np.isnan(np.max(np.abs(gfp.vec.FV().NumPy()))):
raise Exception("Evolution solver failed")
# NOTE: this problem is not a minimization problem. Therefore, we cannot just sum up (W * R).
err = np.max(np.abs(w.FV().NumPy() * rhs.FV().NumPy()))
print("step ", i, "err = ", err)
# Check convergence
if err < 1e-6: break
print("force = ", ls * force, ", uy_A =", gfu(node_A)[1], ", ux_B =", gfu(node_B)[0],\
", int u2 =", Integrate(gfu[1] * ds("top"),mesh))
called base class apply, type = N5ngfem29SymMatrixDifferentialOperatorE
called base class apply, type = N5ngfem29SymMatrixDifferentialOperatorE
called base class apply, type = N5ngfem29SymMatrixDifferentialOperatorE
called base class apply trans, type = called base class apply trans, type = N5ngfem29SymMatrixDifferentialOperatorEN5ngfem29SymMatrixDifferentialOperatorE
called base class apply trans, type = N5ngfem29SymMatrixDifferentialOperatorE
step 0 err = 4.712438806739467
step 1 err = 5.3724340411825175e-28
force = 45.0 , uy_A = 0.020951376702133328 , ux_B = 0.007675790037017989 , int u2 = 2.0403444126180386
step 0 err = 18.84975522695788
step 1 err = 3.982601620357466e-27
force = 135.0 , uy_A = 0.0628541301064 , ux_B = 0.02302737011105396 , int u2 = 6.121033237854117
step 0 err = 18.849755226957644
step 1 err = 6.935654304789211e-07
force = 225.0 , uy_A = 0.10475742852374946 , ux_B = 0.03837870133701924 , int u2 = 10.201748410165306
step 0 err = 18.85048574839874
step 1 err = 2.2717339787904215e-06
step 2 err = 3.797428655664862e-17
force = 315.0 , uy_A = 0.14666331464887977 , ux_B = 0.053728844278525686 , int u2 = 14.282591195628013
step 0 err = 4.712957660941089
step 1 err = 3.4582901589029536e-07
force = 360.0 , uy_A = 0.1676176492296069 , ux_B = 0.06140328097636566 , int u2 = 16.323081028534332
step 0 err = 4.713548847956908
step 1 err = 5.563671983451368e-06
step 2 err = 1.7258537474016867e-18
force = 405.0 , uy_A = 0.1885782876717374 , ux_B = 0.06907508555342708 , int u2 = 18.363848034284345
step 0 err = 1.1788967604087703
step 1 err = 9.357895610238343e-07
force = 427.5 , uy_A = 0.19906271893531674 , ux_B = 0.07290934469916595 , int u2 = 19.384411357322264
step 0 err = 1.1793159351211673
step 1 err = 1.1171362918238072e-06
step 2 err = 1.7928029436945363e-19
force = 450.0 , uy_A = 0.20955157350721115 , ux_B = 0.07674190172012159 , int u2 = 20.40517812520739
[16]:
fd.Assemble()
pdraw.vec.data = invad * fd.vec
Draw(Norm(pdraw), mesh, "p")
Draw(gfu,mesh, "u")
if H.Get() > 1e-16:
fda.Assemble()
adraw.vec.data = invada * fda.vec
Draw(Norm(adraw), mesh, "alpha")
[17]:
# confirm that the trace vanishes
p1 = np.array(pdraw(mesh(90, 100, 0))).reshape((3,3))
np.einsum('ij,ij', p1, np.eye(3))
[17]:
-1.6940658945086007e-21
[ ]:
Implementation based on constrained minimization¶
We again start with the (coupled) function spaces. This time, we won’t need a TrialFunction “variable” for \(\alpha\) but only a GridFunction to store its value. Also, there is no need for \(\Lambda\) in this formulation.
Moreover, we use a deviatoric space for \(p\).
[18]:
fes_p = MatrixValued(fes_ir, symmetric=True, deviatoric=True, dim=3)
fes = fes_u * fes_p
u, p = fes.TrialFunction()
# GridFunction for solution
gfsol = GridFunction(fes)
gfu, gfp = gfsol.components
# Save previous solution
gfsol_k = GridFunction(fes)
gfu_k, gfp_k = gfsol_k.components
# alpha_0 = 0
alpha_k = GridFunction(fes_ir)
The energy desnity \(\Psi\) and the dissipation potential \(\Phi\)
[19]:
# The full energy density (incl. hardening)
def Psi(p, strain, p_k, alpha_k):
strain_energy = elastic_strain_energy(strain-p, Emod, nu)
alpha = alpha_k + IfPos(H - 1e-16, sqrt(2/3) * sigma_Y * H * tensor_norm(p - p_k, pert=norm_pert), 0)
hardening = 1/2 * alpha**2
return strain_energy + hardening
# Dissipation potential: For rate-independent plasticity (homogeneous of degree 1 in the rate of p)
# this is the dissipation! Also, because of rate independence, the parameter Delta_t is not needed for
# the evolution.
def Phi(p, p_k):
return sqrt(2/3) * sigma_Y * tensor_norm(p - p_k, pert=norm_pert)
Finally, we have all ingredients for the variational formulation:
Find for given \((u^k,p^k)\) a solution \((u^{k+1},p^{k+1})\) such that \begin{align*} \Psi(\epsilon(u^{k+1}), p^{k+1},p^{k}) + \Phi(p^{k+1},p^k)-\int f^{k+1}\cdot u^{k+1}\,dx \rightarrow \min! \end{align*}
[20]:
realcompile = False
# Possible optimization: use static condensation for p
a = BilinearForm(fes, symmetric=True)
a += Variation( Psi(p, strain(u), gfp_k, alpha_k) * irs_dx ).Compile(realcompile=realcompile)
a += Variation( Phi(p, gfp_k) * irs_dx ).Compile(realcompile=realcompile)
a += Variation( -u[1] * force * loadfactor * ds("top") ).Compile(realcompile=realcompile)
Since Phi(p, gfp_k) is included above, we automatically obtain an “algorithmically consistent linearization”, correspondig to the local evolution
[21]:
evolution_objective = (Psi(p, strain(gfu), gfp_k, alpha_k) + Phi(p,gfp_k)).Compile(realcompile=False)
which will be used to solve for p in a nested way, ie. for each “guess” for u^{k+1}.
In order to run the problem, we need to define vectors that store respective solution data, define loadsteps and setup a loop that iterates through the latter.
The nonlinear PDE is solved with a Newton-Raphson scheme. Thereby, in each Newton iteration we solve a linearized PDE for the displacements and a nonlinear evolution problem at each quadrature point via ``MinimizationCF``
[22]:
# a set of vectors for keeping data
vec_k = gfsol.vec.CreateVector()
w = gfsol.vec.CreateVector()
rhs = gfsol.vec.CreateVector()
# postprocess internal states
fd = LinearForm(InnerProduct(gfp, qd)*irs_dx)
fda = LinearForm(InnerProduct(gfalpha_k, qda)*irs_dx)
# load steps (chosen based on experience)
loadsteps = [0.1,0.3,0.5,0.7,0.8,0.9,0.95,1.0]
# set solution to zero initially
gfsol.vec[:] = 0
alpha_k.vec[:] = 0
# evolution params (default values of MinimizationCF at the time of writing)
tol = 1e-6
maxiter = 20
# iterate through load steps
for ls in loadsteps:
loadfactor.Set(ls)
# update old solution at time t = t_k
gfsol_k.vec.data = gfsol.vec
with TaskManager():
for i in range(20):
energy_k = a.Energy(gfsol.vec)
#print ("energy(t_k) = ", energy_k)
a.AssembleLinearization(gfsol.vec)
a.Apply(gfsol.vec, rhs)
# static condensation could be employed to speedup solving
w.data = a.mat.Inverse(freedofs=fes.FreeDofs(False), inverse="sparsecholesky") * rhs
# linesearch ( with damping)
vec_k.data = gfsol.vec
scale = 1
while scale > 1e-7:
gfsol.vec.data -= scale * w
energy1 = a.Energy(gfsol.vec)
# Here we solve the evolution equation as each quadrature point through evaluation of
# the MinimizationCF. This function takes the "objective", the initial guess (most
# likely the previous solution), tolerances and the maximum number of iterations.
gfp.Interpolate(MinimizationCF(evolution_objective, gfp, tol=tol, maxiter=maxiter))
energy2 = a.Energy(gfsol.vec)
if energy2 < energy_k + 1e-12:
break
scale *= 0.5
gfsol.vec.data = vec_k
err = sqrt(InnerProduct(w, rhs))
print("step ", i, "err = ", err)
# check convergence
if err < 1e-6: break
alpha_k.Interpolate(alpha_k + sqrt(2/3) * sigma_Y * H * sqrt(InnerProduct(gfp - gfp_k, gfp - gfp_k)))
print("force = ", ls * force, ", uy_A =", gfu(node_A)[1], ", ux_B =", gfu(node_B)[0],\
", int u2 =", Integrate(gfu[1] * ds("top"),mesh))
step 0 err = 9.582056006020684
step 1 err = 1.5781049432801721e-07
force = 45.0 , uy_A = 0.020951445177746744 , ux_B = 0.00767584635396324 , int u2 = 2.040351079705119
step 0 err = 19.164150980208845
step 1 err = 4.699416302057107e-06
step 2 err = 5.408170861741613e-14
force = 135.0 , uy_A = 0.06285441203217194 , ux_B = 0.02302760034959514 , int u2 = 6.121060641138886
step 0 err = 19.164275660610738
step 1 err = 0.0022147679157972677
step 2 err = 5.398996436151884e-09
force = 225.0 , uy_A = 0.1047580956190586 , ux_B = 0.03837924164836408 , int u2 = 10.201813157477611
step 0 err = 19.16446250714923
step 1 err = 0.00845588235670084
step 2 err = 4.333438684226842e-08
force = 315.0 , uy_A = 0.1466646044235403 , ux_B = 0.0537298767581776 , int u2 = 14.282716139638739
step 0 err = 9.583035961812161
step 1 err = 0.006162720625781453
step 2 err = 1.8553093960574856e-08
force = 360.0 , uy_A = 0.1676197537801811 , ux_B = 0.06140492759953414 , int u2 = 16.323283726940634
step 0 err = 9.58329884526734
step 1 err = 0.010615530771988626
step 2 err = 4.158501933676034e-08
force = 405.0 , uy_A = 0.1885812908915135 , ux_B = 0.06907760492794891 , int u2 = 18.364146502989453
step 0 err = 4.793815916797371
step 1 err = 0.007957838920054968
step 2 err = 2.5204653441212257e-08
force = 427.5 , uy_A = 0.1990667443701601 , ux_B = 0.07291289861381442 , int u2 = 19.384822032081203
step 0 err = 4.794109779013835
step 1 err = 0.010140602182098007
step 2 err = 3.952430637757814e-08
force = 450.0 , uy_A = 0.20955661620043797 , ux_B = 0.07674679669147616 , int u2 = 20.40571690418738
[23]:
fd.Assemble()
pdraw.vec.data = invad * fd.vec
Draw(Norm(pdraw), mesh, "p")
Draw(gfu,mesh, "u")
if H.Get() > 1e-16:
fda.Assemble()
adraw.vec.data = invada * fda.vec
Draw(Norm(adraw), mesh, "alpha")
[24]:
# confirm that the trace vanishes
p1 = np.array(pdraw(mesh(90, 100, 0))).reshape((3,3))
np.einsum('ij,ij', p1, np.eye(3))
[24]:
-6.776263578034403e-21
[ ]: