discrepancy between numerical solution and analytical solution for linear elastic material

Is whereZero(x[0]-.5*lx) non-zero?

Hi Lutz,

whereZero(x[1])*whereZero(x[0]-.5*lx)*[1,1] and whereZero(x[1])*whereZero(x[0]-.5*lx)*[0,0] are used to fix the DOFs at the center of the bottom.

 Even if the conditions are set as
# Dirichlet BC positions, smooth at bottom and top
Dbc = whereZero(x[1])*[0,1]+whereZero(x[1]-ly)*[0,1]
# Dirichlet BC values
vel = -0.0001
Vbc = whereZero(x[1])*[0,0]+whereZero(x[1]-ly)*[0,vel]

the results are the as # Dirichlet BC positions, smooth at bottom and top, fixed at the center of bottom
Dbc = whereZero(x[1])*[0,1]+whereZero(x[1]-ly)*[0,1]+whereZero(x[1])*whereZero(x[0]-.5*lx)*[1,1]
# Dirichlet BC values
vel = -0.0001
Vbc = whereZero(x[1])*[0,0]+whereZero(x[1]-ly)*[0,vel]+whereZero(x[1])*whereZero(x[0]-.5*lx)*[0,0]

There still exists the discrepancy.

Best regards,
Jiadun

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

Is still a relevant ?

Hi Lutz,

Yes, it is still relevant.

Best regards,
Jiadun

the problem you are solving in 2D is a plain strain problem. For this case you get:

stressYY = E*vel*t/ly * 1/(1-nu**2)) (nu is the Poisson ration)

To get this answer you need to set the horizontal velocity at x[0]=0 to zero. the facxe x[0]=lx needs to remain free to get zero horizontal stress.

Note: you should set integrationOrder=-1.

Hi Lutz,

I still didn't get the answer.

I think the problem is that the strain in XX direction is not included in the analytical expression.

By the way, the expression stressYY = E*vel*t/ly is one from material mechanics rather than elastic one.

from esys.escript import *
from esys.escript.linearPDEs import LinearPDE
from esys.finley import Rectangle
from esys.weipa import saveVTK

lam = 3.462e9
mu = 3.462e9

# E to be used in analytical expression
E = mu*(3*lam+2*mu)/(lam+mu)

rho=1154.

lx = 0.1
ly = 1.0

dom = Rectangle(l0=lx,l1=ly,n0=1,n1=1,order=1,integrationOrder=-1)

C=Tensor4(0.,Function(dom))
for i in range(dom.getDim()):
 for j in range(dom.getDim()):
  C[i,i,j,j] += lam
  C[i,j,i,j] += mu
  C[i,j,j,i] += mu

pde=LinearPDE(dom,numEquations=dom.getDim(),numSolutions=dom.getDim())
pde.setValue(A=C)

stress=Tensor(0.,Function(dom))
stress[0,0]=-0.0e5
stress[1,1]=-0.0e5

x = dom.getX()
# Dirichlet BC positions, smooth at bottom and top, fixed at the center of bottom
Dbc = whereZero(x[0])*[1,0]+whereZero(x[1])*[0,1]+whereZero(x[1]-ly)*[0,1]+whereZero(x[1])*whereZero(x[0]-0.5*lx)*[1,1]
# Dirichlet BC values
vel = -0.0001
Vbc = whereZero(x[0])*[0,0]+whereZero(x[1])*[0,0]+whereZero(x[1]-ly)*[0,vel]+whereZero(x[1])*whereZero(x[0]-0.5*lx)*[0,0]

u = Vector(0.,Solution(dom))

xInitial = dom.getX()

t = 1

tEnd = 3

while t<tEnd:
 print "t= ",t

 pde.setValue(X=-1.0*stress,q=Dbc, r=Vbc)

 u = pde.getSolution()
 g=grad(u)
 sig=mu*(g+transpose(g))+lam*trace(g)*kronecker(dom)
 stress += sig

 dom.setX(xInitial+u)

 print "Numerical stressYY ", stress[1,1]

 #analytical expression
 stressYY = E*vel*t/(ly*(1-mu**2))
 print "Analytical stressYY " ,stressYY

 print "Ratio : ", stress[1,1]/stressYY

 saveVTK("deform_%d.vtu"%t, disp=u,stress=stress)
 t += 1

Best regards,
Jiadun

Thanks Lutz Gross, that solved my question.

Hi Lutz,

I'm just curious about how to define a plane stress problem in Escript.

Best regards,
Jiadun

The script you listed is for the plane stress case.

Hi Lutz Gross,

I checked that the stressYY = E*vel*t/(ly*(1-nu**2)) (This nu = lam/(2*(lam+mu))!!! not mu) you given is for plane strain.
and the tangent stiffness C given in my scripts is also for plane strain.

One more question is that I don't how to derive the expression

u[0]= -vel*t*x[0] * nu/(1-nu)/ly

Could you please explain it in more detail?

Thank you very much.

Best regards,
Jiadun

Thanks Lutz Gross, that solved my question.