Non-linear transient heat equation

[Dave]

Dear all,

I am trying to solve the non-linear transient heat equation.

When I execute the below MWE, the resultant temperature field is discontinuous (see picture attached). However, when the thermal properties (rho, Cv, conductivity) are set equal to 1, the resultant temperature field is continuous (see picture attached).

Could anyone please offer any advice as to why the PDE becomes ill-conditioned with these parameters? and how I could fix it?

Many thanks

Dave

from fenics import *
import mshr

#Define unstructured example mesh
domain = mshr.Box(Point(0,0,0), Point(1,1,1))
mesh = mshr.generate_mesh(domain, 40)

#Boundary condition for heat flux
class S1(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0],0)

S1 = S1()
Surfaces = MeshFunction('size_t', mesh, 2)
Surfaces.set_all(0)
S1.mark(Surfaces, 1)
ds = Measure("ds")(subdomain_data=Surfaces)

#Define X-functions
A = FunctionSpace(mesh, "CG", 1)

Tn = interpolate(Constant(0), A)
Tn_plus_1 = Function(A)
vT = TestFunction(A)
duT = TrialFunctions(A)

#Thermal parameters
rho = 3000
Cv = 800
conductivity = 20

timeStep = 1
for i in range(2):
    
    #Weak form of heat equation
    thermalForm = rho*Cv*Tn_plus_1*vT*dx - rho*Cv*Tn*vT*dx + Constant(timeStep)*inner(conductivity*grad(Tn_plus_1),grad(vT))*dx - Constant(timeStep)*inner(Constant(1e6),vT)*ds(1)
    
    #Construct non-linear variational problem
    J = derivative(thermalForm, Tn_plus_1, duT)
    problem = NonlinearVariationalProblem(thermalForm, Tn_plus_1, J=J)
    solver = NonlinearVariationalSolver(problem)
    solver.solve()
    
    #Assign current solution
    Tn.assign(Tn_plus_1)
    
    #Save time step to file
    File('Visualisation/T' + str(i) + '.pvd') << Tn_plus_1