I am trying to solve the 1D polymer curing problem as stated in the COMSOL blog:
https://www.comsol.com/blogs/modeling-the-thermal-curing-process/
It is a coupled heat equation with a cure equation to find two scalar fields: one is the temperature and the other the degree of cure alpha
By modifying the Fenics poisson PDE example:
http://home.simula.no/~hpl/homepage/fenics-tutorial/release-1.0/webm/timedep.html
To do this I'm basically defining two scalar fields, alpha_1 and alpha (the degree of cure for two different time steps, and computing its derivative as a finite difference, which is then updated to the source term f. However when I evaluate this I get a long expression error.
Does it have to do with the degree value in the Expression?
The code is shown below:
"""
Heat transfer problem based on:
FEniCS tutorial demo program: Diffusion equation with Dirichlet
conditions and a solution that will be exact at all nodes.
This is for 1D polymer thermal curing, using same setup as shown below for COMSOL
https://www.comsol.com/blogs/modeling-the-thermal-curing-process/
"""
from __future__ import print_function
from dolfin import *
import numpy as np
#Constants (all in SI units)
rho = 1200 #density
C_p = 1000 #specific heat
k = 0.2 #thermal conductivity
H_r = 5.0e5 #total heat of reaction
A = 2.0e4 #frequency factor
E_a = 5.0e4 #activation energy
R = 8.314 #universal gas constant
n = 1.4 #order of reaction
l = 5.0e-3 #length of 1D interval
numElems = 20 #number of finite elements
Q_heat = 10.0e3 #10 kW/m^2 for one m^2 which is 10 kW of heating power
T_ambient = 20.0
# Get mesh and define function spaces
mesh = IntervalMesh(numElems, 0.0,l)
V = FunctionSpace(mesh, 'Lagrange', 2)
Q = FunctionSpace(mesh, 'Lagrange', 1)
class HeatFlux(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (x[0] > (l-DOLFIN_EPS))
#Initialize sub-domain instances
right = HeatFlux()
# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
right.mark(boundaries, 1)
ds = Measure("ds")[boundaries]
#Insulation boundary given already and not needed
plot(mesh, title="1D mesh")
interactive()
# Initial condition
u_1 = project(T_ambient, V)
g = Expression('Q_heat', degree=1, Q_heat=Q_heat)
alpha_1 = project(Expression('0.0', degree=1), V)
alpha = alpha_1
dt = 1.0 # time step
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = (alpha - alpha_1)/dt
a = rho*C_p*u*v*dx + dt*k*inner(nabla_grad(u), nabla_grad(v))*dx
L = (u_1 + dt*f)*v*dx + g*v*ds(1)
A = assemble(a) # assemble only once, before the time stepping
b = None # necessary for memory saving assemeble call
uFile = File("tempRise/temperatureDist.pvd")
# Compute solution
u = Function(V) # the unknown at a new time level
T = 8.0*60.0 # total simulation time: 8 minutes
t = dt
while t <= T:
print('time =', t)
b = assemble(L, tensor=b)
#bc.apply(A, b)
solve(A, u.vector(), b)
t += dt
u_1.assign(u)
alpha_1.assign(alpha)
alpha = project(Expression('-rho*H_r*A*exp(-E_a/R/u)*(1-alpha)^n', degree=2, rho=rho,H_r=H_r,A=A,E_a=E_a,R=R,n=n, u=u,alpha=alpha), V)
f.assign((alpha - alpha_1)/dt)
uFile << u
File("finalTemperatureDist.pvd") << u
plot(u, title="Temperature")
interactive()
