Issue with Thermal Source in Time Dependent Heat Equation

variational formulation
1

Hi everyone,

I am fairly new to Fenics and I am trying to model a 2D heat equation, accounting for radiation and conduction from a long cylindrical rod. The top face of the rod is illuminated by a Gaussian heat source of the form:
Q = \frac{2 A P_{max}}{\pi w^2} e^{\frac{-2(x-1.5)^2}{w^2}}
Where A is the absorptivity of the rod, P_{max} is the power of the Guassian heat source and w is the width of the Gaussian. The issue I am having is how to to implement this into a time dependent problem, as depending on the number of time steps I have, the final heat distribution also changes. For example, for t=900s with 150 time steps, I get this (in Kelvin):test2D

However, if I have 400 time steps. I get this:
test2D_2

I think this could just be a result of the coarseness of the timesteps, however the fact that I am getting negative Kelvin temperatures is definitely wrong, so there maybe an issue with my weak form differential or how my thermal source is implemented.

Any suggestions on how I could fix this? Thanks in advance

My code is shown below:

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

time = 900            # final time
num_steps = 400     # number of time steps
dt = time / num_steps # time step size

#Configure for dimensions of system
#build rectangular mesh to model 2D slice
nx = 15
ny = 50
x0 = 0
y0 = 0
x1 = 3  #mm
y1 = 10    #mm
pi = 3.141592653589793
T_am = Constant("0") #ambient vacuum temperature
T_a4 = T_am**4
epsilon = 0.05  # material emissivity
sigma = 5.67E-14 # W/(mm2.K4)
es = epsilon*sigma
Q = Expression('2*A*Pmax/(pi*w*w)*exp(-2*pow(x[0]-1.5, 2)/(w*w))',degree=2, A=0.4, Pmax=400,w=0.56) #power density
#----------THERMAL--PROPERTIES----------
kappa = 0.03        #conductivity [W/mm K] 
c = 0.132         #Specific Heat Capacity [J/gK]
rho = 0.0195            #Density [g/mm^3]
const = kappa /(c * rho)
#---------------------------------------
 
sw = Point (x0, y0)
ne = Point (x1,y1)

mesh = RectangleMesh(sw, ne, nx, ny, diagonal= "right")
V = FunctionSpace(mesh, 'P', 2)

# Define boundary condition
# build expression for power guassian (x0  the beam from the centre, A is the absorptivity of the target, Pmax is the maximum power of the beam, w is the width of the beam.)


#-----BOUNDARY-CONDITIONS--------------------------


#-----Define Markers for the different parts of the boundary-------
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
tol = 1E-14
class BoundaryX0(SubDomain):
    tol = 1E-14
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], x0, tol)
    
bx0 = BoundaryX0()
bx0.mark(boundary_markers, 0)

class BoundaryX1(SubDomain):
    tol = 1E-14
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], x1, tol)
    
bx1 = BoundaryX1()
bx1.mark(boundary_markers, 1)

class BoundaryY0(SubDomain):
    tol = 1E-14
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], y0, tol)
    
by0 = BoundaryY0()
by0.mark(boundary_markers, 2)

class BoundaryY1(SubDomain):
    tol = 1E-14
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], y1, tol)
    
by1 = BoundaryY1()
by1.mark(boundary_markers, 3)


# For the implimentation of these boundary conditions to be general, we can let the user specify what kind of boundary condition that applies to each of the four boundaries. We set up a Python dictionary for this purpose.


#ENCODE BOUNDARY CONDITIONS (MAY NEED ALTERING OR REMOVAL)
boundary_conditions = {0: {'Robin':     (es, T_a4)},
                       1: {'Robin':     (es, T_a4)},
                       2: {'Robin':     (es, T_a4)},
                       3: {'Robin':     (es, T_a4)}}


ds = Measure('ds', domain=mesh, subdomain_data = boundary_markers)
integrals_R = []
T = TrialFunction(V)
v = TestFunction(V)
T_ = Function(V)
T_n = interpolate(T_am, V)
#Account for Radiation
for i in boundary_conditions:
    if 'Robin' in boundary_conditions[i]:
        integrals_R.append(es*(T**4 - T_a4)*v*dt*ds(i)) #t-t_am will work for now to make problem solvable

#----------------------------------------
vtkfile = File('test2D/solution.pvd')
F = -T*v*dx - const*dt*dot(grad(T), grad(v))*dx + const*Q*v*dt*ds(3) - const*sum(integrals_R) -T_n*v*dx
T_= Function(V)   # the most recently computed solution
F = action(F,T_)
J = derivative(F, T_,T)
t = dt
bcs = []
while t <= time:
    print ('time =', t)
    problem = NonlinearVariationalProblem(F, T_, bcs, J)
    solver = NonlinearVariationalSolver(problem)
    solver.solve()
    prm = solver.parameters
    prm['newton_solver']['absolute_tolerance'] = 1E-9
    prm['newton_solver']['relative_tolerance'] = 1E-9
    prm['newton_solver']['maximum_iterations'] = 100
    prm['newton_solver']['relaxation_parameter'] = 1.0
    t += dt
    T_n.assign(T_)
    vtkfile << (T_, t)
    plot(T_)

p = plot(T_)
p.set_cmap("inferno")
#p.set_clim(0,100)
plt.colorbar(p)
plt.savefig("test2D_2.png")
plt.savefig("test2D_2.pdf")