Hello @yschapira,
thank you for the suggestion. I will look into the package.
Here is a MWE with just the necessary functions from the tutorials. The geometry is the same from the vonMises_Plasticity tutorial.
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
parameters["form_compiler"]["representation"] = 'quadrature'
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)
E = Constant(70e3)
nu = Constant(0.3)
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2./(1+nu)
sig0 = Constant(250.) # yield strength
Et = E/100. # tangent modulus
H = E*Et/(E-Et) # hardening modulus
Re, Ri = 1.3, 1. # external/internal radius
mesh = Mesh("thick_cylinder.xml")
facets = MeshFunction("size_t", mesh, "thick_cylinder_facet_region.xml")
ds = ds(subdomain_data=facets)
deg_u = 2
deg_stress = 2
V = VectorFunctionSpace(mesh, "CG", deg_u)
We = VectorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, dim=4, quad_scheme='default')
W = FunctionSpace(mesh, We)
W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, quad_scheme='default')
W0 = FunctionSpace(mesh, W0e)
sig = Function(W)
sig_old = Function(W)
n_elas = Function(W)
beta = Function(W0)
p = Function(W0, name="Cumulative plastic strain")
u = Function(V, name="Total displacement")
du = Function(V, name="Iteration correction")
Du = Function(V, name="Current increment")
v = TrialFunction(V)
u_ = TestFunction(V)
u_old = Function(V)
v_old = Function(V)
a_old = Function(V)
bc = [DirichletBC(V.sub(1), 0, facets, 1), DirichletBC(V.sub(0), 0, facets, 3), DirichletBC(V, Constant((0., 0.)), facets, 2)]
rho = Constant(3000.e-6)
eta_m = Constant(0.)
eta_k = Constant(0.)
alpha_m = Constant(0.2)
alpha_f = Constant(0.4)
gamma = Constant(0.5+alpha_f-alpha_m)
beta_ga = Constant((gamma+0.5)**2/4.)
T = 0.01
Nsteps = 100
dt = Constant(T/Nsteps)
cutoff_Tc = 2e-3
p0 = 200.
n = FacetNormal(mesh)
q_lim = float(2/sqrt(3)*ln(Re/Ri)*sig0)*1.2
loading = Expression("-p*tanh(10*t/tc)", p=p0, tc=cutoff_Tc, t=0, degree=2)
def eps(v):
e = sym(grad(v))
return as_tensor([[e[0, 0], e[0, 1], 0],
[e[0, 1], e[1, 1], 0],
[0, 0, 0]])
def sigma(eps_el):
return lmbda*tr(eps_el)*Identity(3) + 2*mu*eps_el
def as_3D_tensor(X):
return as_tensor([[X[0], X[3], 0],
[X[3], X[1], 0],
[0, 0, X[2]]])
ppos = lambda x: (x+abs(x))/2.
def proj_sig(deps, old_sig, old_p):
sig_n = as_3D_tensor(old_sig)
sig_elas = sig_n + sigma(deps)
s = dev(sig_elas)
sig_eq = sqrt(3/2.*inner(s, s))
f_elas = sig_eq - sig0 - H*old_p
dp = ppos(f_elas)/(3*mu+H)
n_elas = s/sig_eq*ppos(f_elas)/f_elas
beta = 3*mu*dp/sig_eq
new_sig = sig_elas-beta*s
return as_vector([new_sig[0, 0], new_sig[1, 1], new_sig[2, 2], new_sig[0, 1]]), \
as_vector([n_elas[0, 0], n_elas[1, 1], n_elas[2, 2], n_elas[0, 1]]), \
beta, dp
def sigma_tang(e):
N_elas = as_3D_tensor(n_elas)
return sigma(e) - 3*mu*(3*mu/(3*mu+H)-beta)*inner(N_elas, e)*N_elas-2*mu*beta*dev(e)
metadata = {"quadrature_degree": deg_stress, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)
def m(u, u_):
return rho*inner(u, u_)*dxm
def k(u_):
return inner(eps(u_), as_3D_tensor(sig))*dxm
def c(u, u_):
return eta_m*m(u, u_) + eta_k*k(u_)
def Wext(u_):
return loading*dot(n, u_)*ds(4)
def update_a(u, u_old, v_old, a_old, ufl=True):
if ufl:
dt_ = dt
beta_ = beta_ga
else:
dt_ = float(dt)
beta_ = float(beta_ga)
return (u-u_old-dt_*v_old)/beta_/dt_**2 - (1-2*beta_)/2/beta_*a_old
def update_v(a, u_old, v_old, a_old, ufl=True):
if ufl:
dt_ = dt
gamma_ = gamma
else:
dt_ = float(dt)
gamma_ = float(gamma)
return v_old + dt_*((1-gamma_)*a_old + gamma_*a)
def update_fields(u, u_old, v_old, a_old):
u_vec, u0_vec = u.vector(), u_old.vector()
v0_vec, a0_vec = v_old.vector(), a_old.vector()
a_vec = update_a(u_vec, u0_vec, v0_vec, a0_vec, ufl=False)
v_vec = update_v(a_vec, u0_vec, v0_vec, a0_vec, ufl=False)
v_old.vector()[:], a_old.vector()[:] = v_vec, a_vec
u_old.vector()[:] = u.vector()
def avg(x_old, x_new, alpha):
return alpha*x_old + (1-alpha)*x_new
a_new = update_a(v, u_old, v_old, a_old, ufl=True)
v_new = update_v(a_new, u_old, v_old, a_old, ufl=True)
res = m(avg(a_old, a_new, alpha_m), u_) + c(avg(v_old, v_new, alpha_f), u_) \
+ k(u_) - Wext(u_) + inner(eps(v), sigma_tang(eps(u_)))*dxm
a_form = lhs(res)
L_form = rhs(res)
def local_project(v, V, u=None):
dv = TrialFunction(V)
v_ = TestFunction(V)
a_proj = inner(dv, v_)*dxm
b_proj = inner(v, v_)*dxm
solver = LocalSolver(a_proj, b_proj)
solver.factorize()
if u is None:
u = Function(V)
solver.solve_local_rhs(u)
return u
else:
solver.solve_local_rhs(u)
return
u_tip = np.zeros((Nsteps+1, 3))
file_results = XDMFFile("plasticity_results.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
P0 = FunctionSpace(mesh, "DG", 0)
p_avg = Function(P0, name="Plastic strain")
time = np.linspace(0, T, Nsteps+1)
Nitermax, tol = 200, 1e-7 # parameters of the Newton-Raphson procedure
results = np.zeros((Nsteps+1, 2))
for (i, t) in enumerate(np.diff(time)):
t = time[i+1]
print("Time: ", t)
loading.t = t-float(alpha_f*dt)
A, Res = assemble_system(a_form, L_form, bc)
nRes0 = Res.norm("l2")
nRes = nRes0
Du.interpolate(Constant((0, 0)))
print("Increment:", str(i+1))
niter = 0
while nRes/nRes0 > tol and niter < Nitermax and nRes > 1e-8:
solve(A, du.vector(), Res, "mumps")
update_fields(u, u_old, v_old, a_old)
Du.assign(Du+du)
deps = eps(Du)
sig_, n_elas_, beta_, dp_ = proj_sig(deps, sig_old, p)
local_project(sig_, W, sig)
local_project(n_elas_, W, n_elas)
local_project(beta_, W0, beta)
A, Res = assemble_system(a_form, L_form, bc)
nRes = Res.norm("l2")
print(" Residual:", nRes, "niter:", niter)
niter += 1
u.assign(u+Du)
p.assign(p+local_project(dp_, W0))
sig_old.assign(sig)
file_results.write(u, t)
p_avg.assign(project(p, P0))
file_results.write(p_avg, t)
p.t = t
if MPI.comm_world.size == 1:
u_tip[i+1, 0] = u(Ri, 0)[0]
if MPI.comm_world.size == 1:
plt.figure()
plt.plot(time, u_tip[:, 0])
plt.xlabel("Time")
plt.ylabel("Tip displacement_0")
plt.show()
The code would give some results but I am not sure if it’s correct. I am a bit skeptical about the residual as well. Hope that helps!
