I am interested in solving the steady heat equation with a non-uniform source term and a convective Robin boundary condition as follows:

k \Delta T + q'''=0 \text{ on } \Omega = [0, 1] \times [0, 1]

-\dfrac{\partial T}{\partial n} = h(T-T_0) \text{ on } \partial \Omega

I am using a method similar to the one documented in the pure Neumann boundary condition example to come up with the missing constant. Because the coefficient k is a material property in this model, I want to make it dependent on T and therefore, the equation becomes non-linear. I have included minimum working example code below. The problem I am having is that when the iterative method converges, there are non-physical temperatures, namely T<0. I donât know if this is the method or my implementation in FEniCS and any suggestions, critiques, or advice would be welcome and appreciated!

```
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from dolfin import *
mesh = UnitSquareMesh(32, 32)
class Source(UserExpression):
def eval(self, values, x):
if abs(x[0]) <= 0.05 and abs(x[1]) <= 0.05:
values[0] = 1.0
else:
values[0] = 0.0
# Thermal conductivity
class thermal_conductivity(UserExpression):
def __init__(self, temp, **kwargs):
super().__init__(kwargs)
self.temp = temp
temp_array = np.linspace(0, 1000, 100)
cond_array = np.linspace(10, 30, len(temp_array))
self.tc = interp1d(temp_array, cond_array, fill_value="extrapolate")
def eval(self, value, x):
value[0] = self.tc(self.temp(x))
def value_shape(self):
return ()
# Select finite element space
CG1_elem = FiniteElement("CG", mesh.ufl_cell(), 1)
R_elem = FiniteElement("R", mesh.ufl_cell(), 0)
element = MixedElement([CG1_elem, R_elem])
V = FunctionSpace(mesh, element)
# Define Test and Trial functions
(u, c) = TrialFunctions(V)
(v, d) = TestFunctions(V)
# Define initial guess at trial Function
u0 = interpolate(Constant(20.0), V.sub(0).collapse())
c0 = interpolate(Constant(0.0), V.sub(1).collapse())
w_k = Function(V)
assign(w_k, [u0, c0])
(u_k, c_k) = split(w_k)
kT = thermal_conductivity(u_k, degree=1)
F = kT*(inner(nabla_grad(u), nabla_grad(v)) + c*v + u*d)*dx
f = Source(degree=1)
F -= f*v*dx
htc=10.0
u_air = Constant(20.0)
g = htc*(u - u_air)
F -= g*v*ds
a = lhs(F)
L = rhs(F)
# Setup Picard iteration
w = Function(V)
eps = 1.0
tol = 1.0e-8
iterations = 0
maxiter = 25
while eps > tol and iterations < maxiter:
iterations += 1
# Setup iterative solver
solve(a == L, w)
u, c = w.split()
e = u - u_k
eps = assemble(e**2*dx)
print(f'iter={iterations}: norm={eps}')
kT.temp = u
w_k.assign(w) # update for next iteration
ax = plt.subplot(1,1,1)
c = plot(u)
cax = plt.colorbar(c, ax=ax)
plt.savefig(f'T_h{int(htc)}.png', dpi=300)'''
```