### Plasticity python code

299

views

2

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.

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"]["quadrature_degree"] = 4
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
deps=sym(grad(u-u0))
dsigma=2*mu*deps+lamda*tr(deps)*Identity(3)
Sigma_tr=Sigma0+dsigma
F=inner(Sigma_tr, grad(u_test))*dx-dot(f, u_test)*dx
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.)
# radial return mapping
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))
FF=inner(Sigma, grad(u_test))*dx-dot(f, u_test)*dx
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

### 2 Answers

0

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.

https://bitbucket.org/fenics-apps/fenics-solid-mechanics/src/9309635e547ccaefd35850ce1b5b533c3ca05264/python/demo/beam3d/demo-beam3d.py?at=master&fileviewer=file-view-default

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...

Any advise would be appreciated!

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...

Any advise would be appreciated!

written
4 months ago by
Yuxiang Lin

have you tried with quadrature degree = 2 or may be 3 ?

written
4 months ago by
Ovais

Haven't yet. You mean quadrature element for stress tensor and yield stress?

written
4 months ago by
Yuxiang Lin

0

Please login to add an answer/comment or follow this question.