strain energy linear elastic

Question

Hi, I am solving a linear elastic problem, I have to solve and get strain energy. I have calculated but i am not getting the correct solution. It will be great if anybody can help to figure it out the possible error.

Thanks in advance.

from dolfin import *
n = 50        
mesh  = UnitSquareMesh(n,n)
print "number of unknown",mesh.num_vertices()
elements = mesh.num_cells()
print elements
plot(mesh)
interactive()


W = VectorFunctionSpace(mesh,'CG',1)
u , v = TrialFunction(W), TestFunction(W)
E = 20.8e3
u = 0.3

lmbda, mu = Constant(E*nu/((1.0 + nu )*(1.0-2.0*nu))) , Constant(E/(2*(1+nu)))

class top(SubDomain):
       def inside(self,x,on_boundary):
             tol = 1e-10
             return abs(x[1]-1) < tol and on_boundary

class bottom(SubDomain):
          def inside(self,x,on_boundary):
             tol = 1e-10
             return abs(x[1]) < tol and on_boundary

Top = top()
Bottom = bottom()

fix_b_bottom_x =  DirichletBC(W.sub(0),Constant(0.0),Bottom)
fix_b_bottom_y =  DirichletBC(W.sub(1),Constant(0.0),Bottom)

 disp_topy = DirichletBC(W.sub(0),Constant(0.1),Top)


 bc_disp = [fix_b_bottom_x, fix_b_bottom_y, disp_topy ]

def epsilon(v):
       return 0.5*(grad(v)  + grad(v).T)

 def sigma(u):
       return 2.0*mu*epsilon(u) + lmbda*tr(epsilon(u))*Identity(2)
def en_dens(u):
       str_ele = 0.5*(grad(u) + grad(u).T)
       IC = tr(str_ele)
       ICC = tr(str_ele * str_ele)
       return (0.5*lmbda*IC**2) + mu*ICC


 E_du = inner(grad(v),sigma(u))*dx

 u = Function(W)
 problem_disp = LinearVariationalProblem(lhs(E_du),rhs(E_du),u,bc_disp)
 solver_disp = LinearVariationalSolver(problem_disp)
 solver_disp.solve()

 plot(u, interactive  = True)

 strr = en_dens(u)

 print strr

Answer 1

To obtain the energy density function you need to project en_dens(u) to a scalar function space, let say:

F = FunctionSpace(mesh, 'DG', 0)
W = project(en_dens(u), F)
plot(W, interactive=True)

Comments

Thanks for answer