# Simulation does not converge when material properties inversely proportional to temperature

**URL:** <https://festim.discourse.group/t/simulation-does-not-converge-when-material-properties-inversely-proportional-to-temperature/36>\
**Category:** User Support\
**Created:** [April 18, 2024, 5:26pm UTC](https://festim.discourse.group/t/simulation-does-not-converge-when-material-properties-inversely-proportional-to-temperature/36 "2024-04-18T17:26:51Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![VVKulagin](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/vvkulagin/32/51_2.png) [@VVKulagin](https://festim.discourse.group/u/VVKulagin)\
**Post date:** [April 18, 2024, 5:26pm UTC](https://festim.discourse.group/t/simulation-does-not-converge-when-material-properties-inversely-proportional-to-temperature/36/1 "2024-04-18T17:26:51Z")

</div>

Hello all.

I tried setting the temperature-dependent material properties for `HeatTransferProblem` using the polynomial expressions from [this](https://www.sciencedirect.com/science/article/pii/S2352179117300388) paper.

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

MWE:

```auto
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.

---

<div class="post-metadata">

**Author:** ![remidm](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/remidm/32/4_2.png) [@remidm](https://festim.discourse.group/u/remidm)\
**Post date:** [April 18, 2024, 6:11pm UTC](https://festim.discourse.group/t/simulation-does-not-converge-when-material-properties-inversely-proportional-to-temperature/36/2 "2024-04-18T18:11:25Z")

</div>

@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.

```python
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:

```python
import matplotlib.pyplot as plt

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

```

![image](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/c1ffb9d2a2fabb62f49b9f7b6744c6ecc9ee83cb.png)

Alternatively you can also do:

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

```