Possible scaling numerical precision errors in heat equation 1D problem

Question

I'm trying to solve a Poisson problem with the same physical constants as the problem solved in COMSOL described here:
https://www.comsol.com/blogs/modeling-the-thermal-curing-process/

I get unphysical answers as well as 0 temperature or NaN values, depending on the scaling factor 'sf' or the time step 'dt'.

The program doesn't seem to be able to enforce the initial conditions; for example if i set it to 20 degrees and 100 degrees on the right as a bc, I get the solution after a few seconds curving and reaching 100 at the right as expected, but the value on the left can be lower than 20 degrees.

As I understand dolfin an insulating bc on the right means a 0 Newman bc which is 'already' enforced. Is this the case?

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt

#Constants (all in SI units)
sf = 1.0e-3 #length scale factor (m to cm)

rho = 1200.0/sf/sf/sf    #density
C_p = 1000.0*sf*sf    #specific heat
k   = 0.2/sf     #thermal conductivity
H_r = 5.0e5*sf*sf   #total heat of reaction

k_A = 20.0e3   #frequency factor
E_a = 50.0e3*sf*sf   #activation energy
R   = 8.314*sf*sf   #universal gas constant
n   = 1.4     #order of reaction

l   = 5.0e-3*sf  #length of 1D interval

numElems = 100 #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 = Expression('20+273.15', degree=2)  #In K!


# Get mesh and define function spaces
mesh = IntervalMesh(numElems, 0.0,l)

V = FunctionSpace(mesh, 'Lagrange', 2)
Q = FunctionSpace(mesh, 'Lagrange', 2)


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 = interpolate(T_ambient, V)

g = Expression('Q_heat', degree=2, Q_heat=Constant(Q_heat))

dt = 1.0      # time step

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

f = project(Expression('0.0', degree=2), V)

a = rho*C_p*u*v*dx + dt*k*inner(nabla_grad(u), nabla_grad(v))*dx
L = (u_1 - dt*rho*H_r*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 = 60*8.0     # total simulation time: 8 minutes
t = dt

k2 =  E_a/R

#dof_coordinates = np.transpose(mesh.coordinates())

while t <= T:
    print 'time = %f' % t
    b = assemble(L, tensor=b)

    #bc.apply(A, b)
    solve(A, u.vector(), b)

    t += dt

    u_1.assign(u)

    f = interpolate(Expression('A*exp(-k2/u)*pow(1-alpha, n)',degree=2,A=Constant(k_A),k2=Constant(k2),n=Constant(n),u=u,alpha=f), V)
    L = (u_1 - dt*rho*H_r*f)*v*dx + g*v*ds(1)

    #uArray = u.vector().array()


#File("finalTemperatureDist.pvd") << u
#plt.plot(dof_coordinates, uArray)
#plt.show()

plot(u, title="Temperature")
interactive()

Answer 1

Why do you even interpolate f? You should be able to use the Expression directly in the form.

More importantly, it looks like you are having a non-linear problem, so you should use a non-linear solver.

Comments

No, you need to use UFL functions, so no interpolation or Expression.

---

I looked at the UFL manual and wrote the implicit scheme like this:

F = alpha - alpha_1 + dt*k_A*exp(-k2/u)*(1-alpha)**n
J = derivative(F, alpha)
solve(F == 0, alpha, J)

To use a Newton method. But when I use the derivative I get on the derivative line: AttributeError: 'Sum' object has no attribute 'arguments'

And without it I get on the solve line:
NotImplementedError: Wrong number or type of arguments for overloaded function 'la_solve'.

---

You need to multiply with a test function and dx. Furthermore, you should use a mixed formulate, so that you can solve the coupled equations.

---

I'm not sure what I should be multiplying by a test function and dx. The issue is that the cure equation is not a variational equation in space. Is that just to let DOLFIN know it is a variational form? I could add + 0.0vdx at the end? I don't think using a mixed space is needed. I saw that being used in one of the old Navier Stokes examples but not in the new one: https://fenicsproject.org/pub/tutorial/html/._ftut1009.html#ftut1:NS

---

I used a mixed formulation like in the advection-diffusion example, but now I'm getting:
TypeError: unsupported operand type(s) for +: 'Sum' and 'Form'
I posted the updated code and this error in another question