Hi guys,
I wrote a code to solve a time dependent thermoelasticity problem. Unfortunately, I can't just copy the stress formula from the static problem to the dynamic problem because of arity mismatches. I just copied my complete code and on the bottom you'll find the error message I get. It would be really really nice, if someone could help me!!!!!
from dolfin import *
Form compiler options
parameters["form_compiler"]["cpp_optimize"] = True
params={"optimize":True, "eliminate_zeros":True,
"precompute_basis_const":True, "precompute_ip_const":True}
! /usr/bin/env python
-- coding: utf-8 --
Time definitions
t = 0.0 # initial time
T = 10.0 # final time
num_steps = 50 # number of time steps
dt = T / num_steps # time step size
Define geometry and mesh
p = Point(0.0, 0.0, 0.0)
q = Point(2.0, 1.0, 1.0)
mesh = BoxMesh(p, q, 12, 8, 8)
Define boundaries
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left, right, bottom, top = 1, 2, 3, 4
CompiledSubDomain("near(x[0], side) && on_boundary", side = p[0]).mark(boundaries, left)
CompiledSubDomain("near(x[0], side) && on_boundary", side = q[0]).mark(boundaries, right)
CompiledSubDomain("near(x[1], side) && on_boundary", side = p[1]).mark(boundaries, bottom)
CompiledSubDomain("near(x[1], side) && on_boundary", side = q[1]).mark(boundaries, top)
Surface integral element
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
Define function space
P = VectorElement("Lagrange", mesh.ufl_cell(), 1, 3)
Q = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, P*Q)
Define trial and test functions
(u, theta) = TrialFunctions(V)
(delta_u, delta_theta) = TestFunctions(V)
Collapse mixed function space
V_u = V.sub(0).collapse()
V_theta = V.sub(1).collapse()
Interpolate displacement, velocity and temperature
u_n = Function(V_u)
u_n.interpolate(Expression(("0.0", "0.2x[0]x[1]", "0.0"), degree = 2))
v_n = Function(V_u)
v_n.interpolate(Constant((0.0, 0.0, 0.0)))
theta_n = Function(V_theta)
theta_n.interpolate(Constant((0.0)))
Volume force/ heat source and prescribed tractions/ prescribed heat fluxes
b = Constant((0.0, 0.0, 0.0))
t_p = Constant((0.0, 0.0, 0.0))
r = Constant(0.0)
q_p = Constant(0.0)
Dirichlet boundary conditions
bcs = [DirichletBC(V.sub(0), Constant((0.0, 0.0, 0.0)), boundaries, left),
DirichletBC(V.sub(0), Constant((0.0, 0.0, 0.0)), boundaries, right),
DirichletBC(V.sub(1), Constant(0.0), boundaries, left)]
Material parameters
E = 215e9
nu = 0.3
mu = E/(2.0(1.0 + nu))
lmbda = Enu/((1.0 + nu)(1.0 - 2.0nu))
kappa = 50.0
alpha = 11.1e-6
temp_0 = 273.0
c_E = 461.0
roh = 7850.0
def beta(E, alpha):
return Ealphatr(Identity(3))*Identity(3)
Stress tensor (linear isotropic elasticity)
def sigma(u, theta):
return lmbdatr(sym(grad(u)))Identity(3) + 2.0musym(grad(u)) - alpha(3.0lmbda + 2.0mu)Identity(3)*theta
Weak form a==l
a = temp_0dtinner(beta(E, alpha), sym(grad(u)))delta_thetadx + dtc_Ethetadelta_thetadx + dtdtkappainner(grad(theta), grad(delta_theta))dx + dtdtinner(sym(grad(delta_u)), sigma(u, theta))dx + rohinner(delta_u, u)dx
l = dtdtinner(r, delta_theta)dx + dtdtinner(q_p, delta_theta)ds(top) + dtdtinner(b, delta_u)dx + dtdtinner(t_p, delta_u)ds(top) + rohinner(delta_u, u_n)dx + dtrohinner(delta_u, v_n)dx + dtc_Einner(theta_n, delta_theta)dx + temp_0dtinner(beta(E, alpha), sym(grad(u_n)))delta_theta*dx
Create displacement and temperature file
u_n.rename("u","displacement")
file_u = File("thermoelasticitydynamic_displacement.pvd", "compressed")
file_u << (u_n, t)
file_theta = File("thermoelasticitydynamic_temperature.pvd", "compressed")
theta_n.rename("theta","temperature")
file_theta << (theta_n, t)
Project stress field and create stress file
def dev(s):
return s-tr(s)Identity(3)/3.0
def von_mises(s):
return sqrt(3.0/2.0inner(dev(s), dev(s)))
Z = FunctionSpace(mesh, "Lagrange", 1)
stress = project(von_mises(sigma(u, theta)), Z)
stress.rename("stress","vonMises")
file_stress = File("thermoelasticitydynamic_stress.pvd", "compressed")
file_stress << (stress, t)
Time integration and function solving
x = Function(V)
for n in range(num_steps):
t += dt
solve(a == l, x, bcs)
(u, theta) = x.split(deepcopy=True)
v_n.assign(FunctionAXPY([(1/dt, u), (-1/dt, u_n)]))
u_n.assign(u)
file_u << (u_n, t)
theta_n.assign(theta)
file_theta << (theta_n, t)
stress = project(von_mises(sigma(u, theta)), Z)
file_stress << (stress, t)
ERROR MESSAGE
File "thermoelasticitydynamic.py", line 113, in
stress = project(von_mises(sigma(u, theta)), Z)
raise ArityMismatch("Applying nonlinear operator {0} to expression depending on form argument {1}.".format(o._ufl_class_.name, t))
ufl.algorithms.check_arities.ArityMismatch: Applying nonlinear operator Sqrt to expression depending on form argument v_1.
Aborted
