### Plasticity python code

299
views
2
4 months ago by
Hi,

I am solving a plasticity problem in Fenics python code, which I write the return mapping algorithm and use conditional for to determine yield or not, I post my code below. But it seems not converge (close but not exact)to the cpp version of Fenics solid mechanics code, any one interested in plasticity could found out the potential problems that we could discuss? I use the same mesh file Fenics Solid Mechanics use for beam3d.
#Filename: plasticity.py
from dolfin import *
import numpy as np
parameters["form_compiler"]["cpp_optimize"]= True
parameters["form_compiler"]["optimize"]= True
parameters["form_compiler"]["representation"]= "uflacs"
mesh=Mesh("beam10000.xml.gz")
#Define the scalar space
DGCoeff=FunctionSpace(mesh,"DG",1)
#Define the displacement space
Space=VectorFunctionSpace(mesh, "Lagrange", 2)
#Define the stress space
Stress=TensorFunctionSpace(mesh,"DG",1)

#define space and function------------------------------

#current displacement
u=Function(Space, name='disp')
#previous displacement
u0=Function(Space)
#define trial and test
u_trial=TrialFunction(Space)
u_test=TestFunction(Space)
# local yield stress, scalar
Sigma_Y0=Function(DGCoeff, name='yield')
# local yield stress
Sigma_Y=Function(DGCoeff)
# initilize
Sigma_Y0.vector()[:]=200
Sigma_Y.vector()[:]=200

Sigma0=Function(Stress)

#----------------------------------------------
# define boundary-------------------------------
def DirichletXY(x, on_boundary):
return near(x[0],0.) and near(x[2],0.5)
def DirichletZ(x, on_boundary):
return (near(x[0],0.) and near(x[2],0.5)) or (near(x[0], 5.0) and near(x[2], 0.5))

bc0=DirichletBC(Space.sub(0), Constant(0.), DirichletXY, method="pointwise")
bc1=DirichletBC(Space.sub(1), Constant(0.), DirichletXY, method="pointwise")
bc2=DirichletBC(Space.sub(2), Constant(0.), DirichletZ,  method="pointwise")
bc=[bc0,bc1,bc2]
#------------------------------------------------
#define material property---------------------
nu=0.3
E=Constant(20000.0)
#tangent E
E_t=Constant(20000.*0.3)
H=E_t/(1.-E_t/E)
mu=E/(2.*(1.+nu))
lamda=E*nu/(1.+nu)/(1.-2*nu)
#-------------------------------------------------
#external force
f = Constant([0., 0., 0.])

#total deps=deps_e+deps_p

dsigma=2*mu*deps+lamda*tr(deps)*Identity(3)

Sigma_tr=Sigma0+dsigma

J=derivative(F, u, u_trial)

#trial deviatoric
Shear_tr=Sigma_tr-1./3*tr(Sigma_tr)*Identity(3)
#trial equivalence stress
Sigma_eq_tr=sqrt(1.5*inner(Shear_tr, Shear_tr))
#yield function
F_tr=Sigma_eq_tr-Sigma_Y0
# when F_tr<0 return 0
deps_eq=conditional(F_tr>=0,F_tr/(H+3*mu), 0.)
Sigma=Sigma_tr-3*mu*deps_eq*Shear_tr/Sigma_eq_tr
#deviatoric stress
Shear=Sigma-1./3*tr(Sigma)*Identity(3)
#equivalent stress
Sigma_eq=sqrt(1.5*inner(Shear, Shear))

JJ=derivative(FF, u, u_trial)

F1=XDMFFile('u.xdmf')
F1.parameters['rewrite_function_mesh'] = False
F1.parameters["flush_output"] = True

T = 1.
nsteps = 20
dt = T/nsteps
t = 0.
F1.write(u, 0.)
F1.write(Sigma_Y0, 0.)
#------------------------------------
for n in range(nsteps):
t+=dt
f.assign(Constant([0., 0., -40.0*sin(t/2.*pi)]))
if mpi_comm_world().rank==0:
print("Time "+str(t))
print (-40*sin(t/2.*pi))
print("Solving the trial step")
if n==0:
A, rhs = assemble_system(J, -F, bc)
solve(A, u.vector(), rhs, "mumps")

solve(FF==0, u, bc, J=JJ, solver_parameters={"newton_solver":
{"linear_solver":"mumps",
"absolute_tolerance":1e-7,
"relative_tolerance":1e-7}})

Sigma_Y0.vector()[:]=project(Sigma_Y, DGCoeff, solver_type='gmres').vector()

Sigma0.vector()[:]=project(Sigma, Stress, solver_type="gmres").vector()

u0.vector()[:]=u.vector()
F1.write(u, t)
F1.write(Sigma_Y0, t)

​
Community: FEniCS Project
1
I find your code unreadable. Strip off unnecessary lines, add comment to specify what is the purpose of each line... format the code according to PEP8 specification. Thanks.
written 4 months ago by Michal Habera
Sure, I will modify the comment to make it more readable
written 4 months ago by Yuxiang Lin
Hi, I think I have comment the necessary part of the code, any advise would be appreciated.
written 4 months ago by Yuxiang Lin

0
4 months ago by
I am also interested in solve plasticity in Python.  Fortunately, SolidMechancis now has a python tutorial, I have not compiled it.

A quick compare,   DirichletXY could be different, you use near().
It might be worth of plotting boundary.

element type is also different.
element_t = VectorElement("Quadrature", mesh.ufl_cell(), degree=3, dim=36, quad_scheme=scheme)
Vt = FunctionSpace(mesh, element_t)​

https://bitbucket.org/fenics-apps/fenics-solid-mechanics/src/9309635e547ccaefd35850ce1b5b533c3ca05264/python/demo/beam3d/demo-beam3d.py?at=master&fileviewer=file-view-default
Yeah, near basically the same for setting the centerline of the end the dirichlet bc to zero displacement.

Quadrature element, in my understanding is to speed up the computing, avoid too many unnecessary stored on element. Using Lagrange element is basically not wrong, the convergence rate may be slow. Another reason I want python version is that the cpp version is a bit complicated to couple to other problems, in my project, it is coupled to a shell problem in pure python code...