I am completely new to FEniCS and I am trying to solve linear elasticity equation using DG method:
−div σ = b on a square plate of which one side is fixed and other is pulled with force. But I'm not getting the proper solution. Any help is appreciated. program is as:
from dolfin import *
from mshr import *
import numpy
import math
#====================== generate mesh and function space
mesh = UnitSquareMesh(10, 10)
V = VectorFunctionSpace(mesh ,"DG",2) #function space generation
#====================== define material property
E = 200e9 #young's modulas
nu = 0.3 #poisson ratio
lmbda = Enu /((1.0 + nu)(1.0-2*nu)) #plain strain constant
mu = E / (2 * (1 + nu)) #plain strain constant
#====================== defination of stress and strain
def epsilon(v):
D=nabla_grad(v)
return 0.5 * (D+D.T)
def sigma(v):
return 2.0 * mu * epsilon(v) + lmbda * tr(epsilon(v)) * Identity (2)
#===================== defining boundary condition
tol=1E-14
class left(SubDomain):
def inside(self,x,on_boundary):
return abs(x[0])<=tol and on_boundary
class right(SubDomain):
def inside(self,x,on_boundary):
return abs(x[0]-1)<=tol and on_boundary
left=left()
right=right()
boundaries=FacetFunction("size_t",mesh)
boundaries.set_all(0)
left.mark(boundaries,1)
right.mark(boundaries,2)
Define normal vector and mesh size
n = FacetNormal(mesh)
h = CellSize(mesh)
h_avg = (h('+') + h('-'))/2
u0 = Constant(("0.0", "0.0"))
bc1=DirichletBC(V,u0,boundaries,1)
ds=Measure("ds")[boundaries]
bc=[bc1]
#======================== define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f=Constant(("0.0","0.0")) # body force
g=Expression(("100e8","0.0")) # traction
a = inner(nabla_grad(v), sigma(u))dx \
- inner(avg(sigma(v)n), jump(u))dS \
- inner(jump(v), avg(sigma(u)n))dS \
- inner(avg(sigma(v)n), jump(u))ds(1) \
- inner(jump(v), avg(sigma(u)n))*ds(1)
L= inner(v,f)dx + inner(v,g)ds(2)
#======================== compute solutions
u = Function(V)
solve(a==L, u)
plot(u,mode="displacement",axes=True)
plot(mesh)
interactive()
