### P:lasticity using conditionals

436

views

2

Dear Fenics users,

I am new to Fenics and would like some clarifications on several aspects. I am attaching a python script which is supposed to handle 1-D plasticity based on linear hardening. This is implemented using conditionals. It works for the loading history presented here (there are some lose ends to fix) but I am doubtful of several things. If you are trying to run it please wait for about 15 seconds for it to finish running.

1) Firstly, are conditionals the best way to actually handle plasticity? Do these conditionals apply to all the elements individually. In this 1-D example, the stress everywhere is the same, so it will not make a difference. Suppose certain points were yielding while others were not, would conditionals catch them?

2) Is there a way to assign stress to quadrature point of each element and then loop over element to check if yielded and perform the updates? Like the traditional way? Does this basically defeat the purpose of fenics?

3) Lastly, I tried to install the fenics solid mechanics package from https://bitbucket.org/fenics-apps/fenics-solid-mechanics/src/9309635e547ccaefd35850ce1b5b533c3ca05264/python/?at=master. when I do pip3 -v install --user, the error says, pip3 must have at least one argument to install. Can anyone please help me with the instructions. I use ubuntu and have python3-dolfin and the python2 versions installed installed.

Thanks a lot

I am new to Fenics and would like some clarifications on several aspects. I am attaching a python script which is supposed to handle 1-D plasticity based on linear hardening. This is implemented using conditionals. It works for the loading history presented here (there are some lose ends to fix) but I am doubtful of several things. If you are trying to run it please wait for about 15 seconds for it to finish running.

```
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
mesh = UnitIntervalMesh(1)
V = FunctionSpace(mesh,'Lagrange',1)
disp = Function(V) # Actual displacement
ddisp = Function(V) # Incremental displacement
w = TestFunction(V)
du = TrialFunction(V)
A = 1
p= Point(1)
#=====================================================================================
class NeumannBoundary(SubDomain):
def inside(self,x,on_boundary):
return abs(x[0]) >= 1.0 - DOLFIN_EPS and on_boundary
class DirichletBoundary(SubDomain):
def inside(self,x,on_boundary):
return x[0] <= DOLFIN_EPS and on_boundary
def CT(status):
Y=Constant(10000.0)
H=Constant(1000.0)
K=Constant(Y*H/(Y+H))
return project(conditional(eq(status,1),Y,K),V)
def getbeta(sigmay,s,stress):
print('SIGMAY IN FUNCTION IS=',project(sigmay,V)(p))
print('ABS Stress is=',project(abs(stress),V)(p))
print('Abs sigma is=',project(abs(s),V)(p))
return project((sigmay-abs(s))/(abs(stress)-abs(s)),V)
def calculatesigmay(sigmay,ep):
H = Constant(1000.0)
return project(500 + H*ep,V)
def computestress(s,beta,dstress,dstrain):
return project(s +beta*dstress+CT(0)*(1-beta)*dstrain,V)
def computeplasticstrain(ep,dstrain,beta):
H=Constant(1000.0)
Y=Constant(10000.0)
return project(ep + (1-beta)/(1+H/Y)*abs(dstrain),V)
#=====================================================================================
Gamma_F = NeumannBoundary()
Gamma_u = DirichletBoundary()
bc = DirichletBC(V,Constant(0.0),Gamma_u)
#=====================================================================================
#Initial values
g = Constant(1000.0)
sigmay=Constant(500.0)
H = Constant(1000.0)
Y = Constant(10000.0)
K=Constant(Y*H/(Y+H))
beta=1
iselastic2=Constant(1)
stat=Constant(1)
dstrain = Constant(0.0)
sigma=Constant(0.0)
dstress=Constant(0.0)
ep=Constant(0.0)
stress=Constant(0.0)
iselastic=Constant(1)
beta=Constant(1)
prevep=Constant(0.0)
prevsigma=Constant(0.0)
strain=Constant(0.0)
change=Constant(0)
J=[0.1,0.2,0.3,0.4,0.5,0.51,0.52,0.53,0.54,0.55,0.56,0.57,0.58,0.59,0.6,0.5,0.4,0.3,0.2,0.1,-0.1,-0.2,-0.3,-0.4,-0.5,-0.6,-0.61,-0.62,-0.63,-0.64,-0.65,-0.66,-0.67,-0.68,-0.69,-0.7,-0.69,-0.68,-0.67,-0.68,-0.69,-0.7,-0.71,-0.72,-0.73,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.1]
#J=[0.1,0.2,0.3,0.4,0.52,0.53,0.6,0.7,0.8]
for factor in J:
stat=iselastic
prevsigma=project(sigma,V)
PS=project(prevsigma,V)
print (' ')
print (' ')
print ('#==OUTER FOR LOOP====')
print ('Factor is=',factor)
print ('Prev sigma is=',PS(p))
prevep=project(ep,V)
sigmay=project(calculatesigmay(sigmay,prevep),V)
PS=project(prevsigma,V)
print ('#=======================')
print ('SigmaY is=',sigmay(p))
print ('#=======================')
conv=2
iter = 0
# iselastic=Constant(1)
while conv > 1E-4:
iter = iter +1
if iter==1:
CTL=project(CT(1),V)
else:
CTL=project(CT(iselastic),V)
PS=project(prevsigma,V)
IE=project(iselastic,V)
CTE=project(CT(iselastic),V)
print ('#======INNER WHILE=================')
print ('Prev sigma B4 Solver is 1=',PS(p))
print ('ISELASTIC is=',IE(p))
print ('SIGMA B4 SOLVER is=',project(sigma,V)(p))
print ('MODULUS US=',CTE(p))
a=CTL*A*Dx(du,0)*Dx(w,0)*dx
L=-A*sigma*Dx(w,0)*dx+factor*g*w*ds
solve(a==L,ddisp,bc)
dstrain = Dx(ddisp,0)
dstress=project(Y*dstrain,V)
stress=project(prevsigma+dstress,V)
PS=project(prevsigma,V)
print ('Prev sigma After solver is=',PS(p))
DStr=project(dstrain,V)
print ('dstrain is=',DStr(p))
DS=project(stress,V)
print('Stress is=',DS(p))
Sigma=project(sigma,V)
print('Sigma after solver=',Sigma(p))
kt=project(prevsigma*dstress,V)
cstatus=conditional(lt(project(abs(stress),V),sigmay),1,0)
cstatus2=conditional(lt(kt,0),1,0)
iselastic2=conditional(eq(stat,1),cstatus,cstatus2)
#==============================================================================================
# Compare previous status with new status
change=conditional(And(eq(stat,1),eq(iselastic2,0)),1,0)
CH=project(change,V)
print('Change is=',CH(p),'Changed from elastic to plastic')
cbeta=conditional(lt(project(abs(stress),V),sigmay),1,getbeta(sigmay,prevsigma,stress))
cbeta2=conditional(lt(kt,0),1,0)
beta=conditional(eq(stat,1),cbeta,cbeta2)
BE=project(beta,V)
IEL=project(iselastic,V)
IEL2=project(iselastic2,V)
print('THE BETA IN LOOPS IS=',BE(p))
print('ISELASTIC IN LOOPS IS=',IEL(p))
print('ISELASTIC2 IN LOOPS IS=',IEL2(p))
cstress=conditional(lt(project(abs(stress),V),sigmay),stress,computestress(prevsigma,beta,dstress,dstrain))
cstress2=conditional(lt(kt,0),stress,computestress(prevsigma,beta,dstress,dstrain))
sigma=conditional(eq(stat,1),cstress,cstress2)
cplasstrain=conditional(lt(project(abs(stress),V),sigmay),prevep,computeplasticstrain(prevep,dstrain,beta))
cplasstrain2=conditional(lt(kt,0),prevep,computeplasticstrain(prevep,dstrain,beta))
ep=conditional(eq(stat,1),cplasstrain,cplasstrain2)
iselastic=iselastic2
# sigma=computestress(prevsigma,beta,dstress,dstrain)
# ep=computeplasticstrain(prevep,dstrain,beta)
#============================================================================================
RI=-A*sigma*Dx(w,0)*dx
RE=factor*g*w*ds
TR = RI + RE
RIA=assemble(RI)
REA=assemble(RE)
TRA=assemble(TR)
bc.apply(RIA)
bc.apply(REA)
bc.apply(TRA)
print ('Internal force=',RIA.array())
print ('External force=',REA.array())
conv=np.dot(TRA.array(),TRA.array())/(1+np.dot(REA.array(),REA.array()))
print ('Conv=',conv)
#================================================================================================
strain=project(strain+dstrain,V)
S1=project(sigma,V)
E1=project(strain,V)
EP=project(ep,V)
B=project(beta,V)
PS1=project(prevsigma,V)
STAT=project(stat,V)
print ('===CONVERGED VALUES')
print ('Value of Stress is=',S1(p))
print ('Value of Strain is=',E1(p))
print ('Plastic strain=',EP(p))
print ('STATUS=',STAT(p))
print ('Number of iterations used is=',iter)
print ('#=============DONE LOAD============')
plt.plot(E1(p),S1(p),'bo')
plt.pause(1)
plt.show()
```

1) Firstly, are conditionals the best way to actually handle plasticity? Do these conditionals apply to all the elements individually. In this 1-D example, the stress everywhere is the same, so it will not make a difference. Suppose certain points were yielding while others were not, would conditionals catch them?

2) Is there a way to assign stress to quadrature point of each element and then loop over element to check if yielded and perform the updates? Like the traditional way? Does this basically defeat the purpose of fenics?

3) Lastly, I tried to install the fenics solid mechanics package from https://bitbucket.org/fenics-apps/fenics-solid-mechanics/src/9309635e547ccaefd35850ce1b5b533c3ca05264/python/?at=master. when I do pip3 -v install --user, the error says, pip3 must have at least one argument to install. Can anyone please help me with the instructions. I use ubuntu and have python3-dolfin and the python2 versions installed installed.

Thanks a lot

Community: FEniCS Project

### 2 Answers

0

1.) Meaning of conditional depends on the way it is evaluated, i.e. depends on quadrature degree, quadrature scheme and possibly finite element space that is integrated.

2.) For this purpose FEniCS has finite element family called "Quadrature". Have a look at the fenics solid mechanics package to see how it is used for plasticity. To my best knowledge, it assembles your matrix without any integration - just quadrature points evaluated forms.

3.) I guess there should be a path (e.g. dot),

2.) For this purpose FEniCS has finite element family called "Quadrature". Have a look at the fenics solid mechanics package to see how it is used for plasticity. To my best knowledge, it assembles your matrix without any integration - just quadrature points evaluated forms.

3.) I guess there should be a path (e.g. dot),

`pip3 -v install . --user`

0

Hey:

I am doing plasticity problem also using python version, initially the way I use to determine yield or not is to write a form (u+abs(u))/2, which assumes it equals to zero when part of u is smaller than zero and u when part of u is positive. But it seems when project this form to functionspace, it will smooth out however, when part of u is positive and part of is not. I am looking at quadrature element now to see if it solves the problem.

Best,

Yuxiang

I am doing plasticity problem also using python version, initially the way I use to determine yield or not is to write a form (u+abs(u))/2, which assumes it equals to zero when part of u is smaller than zero and u when part of u is positive. But it seems when project this form to functionspace, it will smooth out however, when part of u is positive and part of is not. I am looking at quadrature element now to see if it solves the problem.

Best,

Yuxiang

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