Adding thermal strain component to elasticity problem

Question

I am trying to solve simple 2D thermoelasticity problem on the rectangular beam. The following code doest work:

mesh = RectangleMesh(Point(0,0),Point(Length,Width),40,20)
E0 = 10e5
nu = 0.28
alpha = 5e-6
Tref = 293
T = 1273
delta_T = T - Tref
E  = E0
mu = E/(2*(1+nu))
lamb = E*nu/((1+nu)*(1-2*nu))
thermal_coef= alpha*delta_T*(3+2*mu/lamb)
# define rectangular beam
# define function space
V = VectorFunctionSpace(mesh,'CG',1)
# Boundary Conditions
tol = 1e-14
def clamped_boundary(x,on_boundary):
    return on_boundary and x[0]<tol
bc = DirichletBC(V,Constant((0,0)),clamped_boundary)
# define strain-displacement 
def epsilon(u):
    return 0.5*(nabla_grad(u)+nabla_grad(u).T)
# define stifness tensor
def sigma(u):
    return lamb*((nabla_div(u)-alpha*delta_T*(3+2*mu/lamb))*Identity(d))+2*mu*epsilon(u)
# elasticity problem
u = TrialFunction(V)
d = u.geometric_dimension()
v = TestFunction(V)
f = Expression(( '0','-20'))
a = inner(sigma(u),epsilon(v))*dx
L = dot(f,v)*dx
F = assemble(L)
F = F.array()
u = Function(V)
solve(a == L,u,bc)
plot(u,mode='displacement')

With an error:

ArityMismatch: Adding expressions with non-matching form arguments ()
vs
(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange',
triangle, 1), dim=2), 24721), VectorElement(FiniteElement('Lagrange',
triangle, 1), dim=2)), 1, None),).

Comments

Hi Antonoguy,
Can you give me the reference where you found the equation for the thermoelastic problems?
Thanks in advance!

Answer 1

Hi, the identity term in sigma has no TrialFunction and thus a is not a left hand side of the equation (the term leads to contribution to right hand side). Consider instead

form = inner(sigma(u),epsilon(v))*dx - dot(f, v)*dx
a = lhs(form)
L = rhs(form)

Comments

Mirok, thanks!