Hello All,
I am trying to implement the Thermoelastodynamics based on FAR91 on the thermoelasticity_transient.py
(FAR91: Farhat, C., Park, K. C., & Dubois-Pelerin, Y. (1991). An unconditionally stable staggered algorithm for transient finite element analysis of coupled thermoelastic problems. Computer Methods in Applied Mechanics and Engineering, 85(3), 349-365. https://doi.org/10.1016/0045-7825(91)90102-C)
The following is my formulation:
Here we have a multiplication of two test functions v (testfunction for vectorspace) and vT (testfunction for scalarspace). As the thermal stress is a function of Temperature and has to be accounted in this equation, and I fail to find reference to how its done.
The following is my code:
from __future__ import print_function
from dolfin import *
from mshr import *
import numpy as np
L = 1.
R = 0.1
N = 50 # mesh density
domain = Rectangle(Point(0.,0.), Point(L, L)) - Circle(Point(0., 0.), R, 100)
mesh = generate_mesh(domain, N)
T0 = Constant(293.)
DThole = Constant(10.)
E = 70e3
nu = 0.3
lmbda = Constant(E*nu/((1+nu)*(1-2*nu)))
mu = Constant(E/2/(1+nu))
rho = 2700. # density
alpha = 2.31e-5 # thermal expansion coefficient
kappa = Constant(alpha*(2*mu + 3*lmbda))
cV = Constant(910e-6)*rho # specific heat per unit volume at constant strain
k = Constant(237e-6) # thermal conductivity
d = 1 # interpolation degree
Vue = VectorElement('CG', mesh.ufl_cell(), d) # displacement finite element
Vte = FiniteElement('CG', mesh.ufl_cell(), d) # temperature finite element
MixedSpace = FunctionSpace(mesh, MixedElement([Vue, Vte]))
VectorSpace = FunctionSpace(mesh,Vue)
ScalarSpace = FunctionSpace(mesh,Vte)
def inner_boundary(x, on_boundary):
return near(x[0]**2+x[1]**2, R**2, 1e-3) and on_boundary
def bottom(x, on_boundary):
return near(x[1], 0, DOLFIN_EPS) and on_boundary
def left(x, on_boundary):
return near(x[0], 0, DOLFIN_EPS ) and on_boundary
bc1 = DirichletBC(VectorSpace.sub(0), Constant(0.), left) # For u(x)
bc2 = DirichletBC(VectorSpace.sub(1), Constant(0.), bottom) # For u(y)
bcT = DirichletBC(ScalarSpace, DThole, inner_boundary) # For Temperature
bcTd = DirichletBC(ScalarSpace, Constant(0.), inner_boundary) # For Temperature Rate
vu = TestFunction(VectorSpace)
vT = TestFunction(ScalarSpace)
u = Function(VectorSpace) #Displacement Current Step
T = Function(ScalarSpace) #Temperature Current Step
un = Function(VectorSpace) #Displacement Previous Step
Tn = Function(ScalarSpace) #Temperature Previous Step
u_dot = Function(VectorSpace) #Velocity Current Step
u_dotn = Function(VectorSpace) #Velocity Previous Step
ud_dot = Function(VectorSpace) #Acceleration Current Step
ud_dotn = Function(VectorSpace) #Acceleration Previous Step
T_dot = Function(ScalarSpace) #Temperature Rate Current Step
T_dotn = Function(ScalarSpace) #Temperature Rate Previous Step
def eps(v):
return sym(grad(v))
def sigma(v, T):
return (lmbda*tr(eps(v)) - kappa*T)*Identity(2) + 2*mu*eps(v)
dt = Constant(0.)
Traction = Constant((0.0,0.0))
f = Constant((0.0,0.0))
Flux = Constant(0.0)
#Equation 2 - 22-23
u = TrialFunction(VectorSpace)
#mech_form = (inner(sigma(u, Tn), eps(vu)))*dx #Static Formulation
bcU = [bc1, bc2]
mech_form = rho*dot(u,vu)*dx - rho*dot(un,vu)*dx - (dt)*(rho*dot(u_dotn,vu)*dx) - (dt*dt/4)*(rho*dot(ud_dotn,vu)*dx) - (dt*dt/4)*dot(Traction,vu)*ds + (dt*dt/4)*inner(sigma(u,Tn),eps(vu))*dx
au, Lu = lhs(mech_form), rhs(mech_form)
u = Function(VectorSpace)
mechproblem = LinearVariationalProblem(au, Lu, u, bcU)
mechsolver = LinearVariationalSolver(mechproblem)
#Equation 1 - 20-21
T = TrialFunction(ScalarSpace)
#therm_form = (cV*(T-Tn)/dt*vT + kappa*T0*tr(eps(u-un))/dt*vT + dot(k*grad(T), grad(vT)))*dx #Backward Euler Formulation
therm_form = (cV)*((2/dt)*(T-Tn)-T_dotn)*vT*dx + dot(k*grad(T), grad(vT))*dx + kappa*T0*tr(eps(u_dotn))*vT*dx
aT, LT = lhs(therm_form), rhs(therm_form)
T = Function(ScalarSpace)
thermalproblem = LinearVariationalProblem(aT, LT, T, bcT)
thermsolver = LinearVariationalSolver(thermalproblem)
#Equation 3 - 24a
ud_dot = TrialFunction(VectorSpace)
v_a = TestFunction(VectorSpace)
bc_ud_dot = [bc1, bc2]
ud_dot_form = rho*dot(ud_dot,v_a)*dx - dot(Traction,v_a)*ds + inner(sigma(un,Tn),(eps(v_a)))*dx - dot(f,v_a)*dx
aa,La = lhs(ud_dot_form), rhs(ud_dot_form)
ud_dot = Function(VectorSpace)
ud_dot_problem = LinearVariationalProblem(aa, La, ud_dot, bc_ud_dot)
ud_dot_solver = LinearVariationalSolver(ud_dot_problem)
#Equation 4 - 24b
u_dot = TrialFunction(VectorSpace)
v_v = TestFunction(VectorSpace)
bc_u_dot = [bc1, bc2]
u_dot_form = dot(u_dot,v_v)*dx - dot(u_dotn,v_v)*dx - (dt/2)*dot((ud_dot + ud_dotn),v_v)*dx
av,Lv = lhs(u_dot_form), rhs(u_dot_form)
u_dot = Function(VectorSpace)
u_dot_problem = LinearVariationalProblem(av, Lv, u_dot, bc_u_dot)
u_dot_solver = LinearVariationalSolver(u_dot_problem)
#Equation 5 - 24c
T_dot = TrialFunction(ScalarSpace)
v_T = TestFunction(ScalarSpace)
bc_T_dot = [bcTd]
T_dot_form = cV*T_dot*v_T*dx + Flux*v_T*ds + dot(k*grad(T), grad(vT))*dx + kappa*T0*tr(eps(u_dotn))*vT*dx
aTd,LTd = lhs(T_dot_form), rhs(T_dot_form)
T_dot = Function(ScalarSpace)
T_dot_problem = LinearVariationalProblem(aTd, LTd, T_dot, bc_T_dot)
T_dot_solver = LinearVariationalSolver(T_dot_problem)
fileResults = XDMFFile('FENICS_FILE/Dynamic.xdmf')
fileResults.parameters["flush_output"] = True
fileResults.parameters["functions_share_mesh"] = True
tt = 0
Nincr = 10
t = np.logspace(1, 4, Nincr+1)
Nx = 100
u.rename("Displacement", "label")
T.rename("Temperature", "label")
fileResults.write(u,tt)
fileResults.write(T,tt)
for (i, dti) in enumerate(np.diff(t)):
dt.assign(dti)
tt+=float(dt)
print("Increment " + str(i+1),"Time Increment",float(dt), "Total Time", tt)
thermsolver.solve() #Solve Equation 20-21
Tn.assign(T)
mechsolver.solve() #Solve Equation 22-23
un.assign(u)
ud_dot_solver.solve() #Update acceleration 24a
u_dot_solver.solve() #Update velocity 24b
T_dot_solver.solve() #Update Heat Rate 24c
ud_dotn.assign(ud_dot)
u_dotn.assign(u_dot)
T_dotn.assign(T_dot)
fileResults.write(u,tt)
fileResults.write(T,tt)
If i use
therm_form = (cV)*((2/dt)*(T-Tn)-T_dotn)*vT*dx + dot(k*grad(T), grad(vT))*dx + kappa*T0*tr(eps(u_dotn))*vT*dx
instead of
therm_form = (cV)*((2/dt)*(T-Tn)-T_dotn)*vT*dx + dot(k*grad(T), grad(vT))*dx + kappa*T0*tr(eps(dot(u_dotn,vu)*dx + (dt/2)*dot((ud_dotn + (dot(Traction,vu)*ds - inner(sigma(u,T),(eps(vu)))*dx + dot(f,vu)*dx)),vu)*dx))*vT*dx , I can solve the problem but the dependence of Temperature within the u_dotn calculation is ignored.
How can I solve this problem?