### P:lasticity using conditionals

436
views
2
8 months ago by
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.
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)
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

0
8 months ago by
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), pip3 -v install . --user
0
5 months ago by
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