Simulation does not converge when material properties inversely proportional to temperature

Hello all.

I tried setting the temperature-dependent material properties for HeatTransferProblem using the polynomial expressions from this paper.

Unfortunately, the simulation model instantly diverges with these expressions. I suppose the issue appears, since the polynomials include terms ∝T-2.

MWE:

import festim as F
import numpy as np

def thermal_cond_function(T):
    return 1 + T**2 + 1 / T

model = F.Simulation()

model.mesh = F.MeshFromVertices(vertices=np.linspace(0, 1, num=1000))

model.settings = F.Settings(
    absolute_tolerance=1e-10, relative_tolerance=1e-10, transient=False
)

mat = F.Material(id=1, D_0=1, E_D=0.2, thermal_cond=thermal_cond_function)
model.materials = mat

model.T = F.HeatTransferProblem(transient=False)

model.initialise()
model.run()

If the last term in thermal_cond_function is removed, the simulation runs without problems.

@VVKulagin the problem is that because of that division, if T=0 then the thermal conductivity is nan. It can be verified by adding model.log_level = 20.

One way to fix this is adding an epsilon to the division so that it’s never zero.

import festim as F
import numpy as np

import fenics as f


def thermal_cond_function(T):
    return 1 + T**2 + 1 / (T + f.DOLFIN_EPS)


model = F.Simulation()

model.mesh = F.MeshFromVertices(vertices=np.linspace(0, 1, num=1000))

model.settings = F.Settings(
    absolute_tolerance=1e-10, relative_tolerance=1e-10, transient=False
)

mat = F.Material(id=1, D_0=1, E_D=0.2, thermal_cond=thermal_cond_function)
model.materials = mat

model.T = F.HeatTransferProblem(transient=False, initial_condition=1)


model.boundary_conditions = [F.DirichletBC(surfaces=[1, 2], value=1 + F.x, field="T")]
model.initialise()
model.run()

Running the following code produces:

import matplotlib.pyplot as plt

f.plot(model.T.T)
plt.show()

image

Alternatively you can also do:

def thermal_cond_function(T):
    return 1 + T**2 + 1 / (T + 1e-12)
1 Like