Steady, nonlinear heat equation with pure neumann BCs

mathematics
1

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)'''

T_h10
T_h101920×1440 115 KB

2

I believe there are two mistakes in the formulation of the weak problem. First, there is a sign error in the boundary term. If we derive a weak problem via integration by parts, we get the following:

\int_\Omega -\nabla\cdot(k\nabla u)v\,d\mathbf{x} = \int_\Omega fv\,d\mathbf{x}
\Rightarrow \int_\Omega k\nabla u\cdot\nabla v\,d\mathbf{x} ~-~\int_{\partial\Omega}(k\nabla u\cdot\mathbf{n})v\,d\mathbf{s}~=~\int_\Omega f v\,d\mathbf{x}
\Rightarrow \int_\Omega k\nabla u\cdot\nabla v\,d\mathbf{x} ~+~\int_{\partial\Omega}h(u - u_\text{air})v\,d\mathbf{s}~=~\int_\Omega f v\,d\mathbf{x}\text{ ,}

under the assumption that the desired BC is

-k\nabla u\cdot\mathbf{n} = h(u - u_\text{air})\text{ ,}

which differs from the equation in the post by a factor of k on the left. (I have followed the notation/conventions used in the code.) Notice that the boundary term enters with a + sign, unlike the code, where it is subtracted.

Second, with a Robin BC like this, there is no need to introduce an auxiliary mechanism to pin down an arbitrary constant. Unlike the pure Neumann demo, if u satisfies the weak problem, u + C (for some nonzero constant C) does not satisfy it, because the boundary term would change.

Another minor point, comparing the code and equations in the post: If k varies in space, the weak form term \int_\Omega k\nabla u\cdot\nabla v\,d\mathbf{x} corresponds to the strong form term -\nabla\cdot(k\nabla u) (as in the derivation above), not -k\Delta u; the first one is what you want, based on deriving the PDE from a physical conservation law.

2 Likes
3

Ah–thanks for the help and pointing out the sign errors! And for the explanation of the Robin BC; I’ve been looking for a reference for that, but hadn’t found one–thanks! For the weak form–I think the code is still correct, even if the typeset formulation isn’t, no? The code below seems to work as I’d expect. Thanks again!

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.5 and abs(x[1]) <= 0.5:
            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)
V = FunctionSpace(mesh, CG1_elem)

# Define Test and Trial functions
u = TrialFunction(V)
v = TestFunction(V)

# Define initial guess at trial Function
u_k = interpolate(Constant(20.0), V)
kT = thermal_conductivity(u_k, degree=1)
F = kT*(inner(nabla_grad(u), nabla_grad(v)))*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
u = 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, u)
    e = u - u_k
    eps = assemble(e**2*dx)
    print(f'iter={iterations}: norm={eps}')

    kT.temp = u
    u_k.assign(u)   # update for next iteration

# Plot solution
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)
4

If you’re referring specifically to the issue of -\nabla\cdot(k \nabla u) vs. -k\Delta u issue, then the first one is what you typically want, so \int_\Omega k\nabla u\cdot\nabla v\,d\mathbf{x} (as in the code) is the “correct” weak form.

For the question of whether or not to include k in the statement of the Robin BC, that depends on what the physical interpretation of h is. (It might be good to do a sanity check of the units.)

1 Like
5

Excellent–thanks again for your help!

6

Hi , I try to run the (corrected) code from that post in Fenics because corresponds to one of my needs. And here is the error message I have:

ModuleNotFoundError: No module named ‘scipy’

Apparently, this line at the beginning of the code is not correct.
from scipy.interpolate import interp1d

How to do to get hte scipy module correctly ?
I am using Fenics through a server + Jupyter Notebooks.

Thanks

7

Scipy is a python package, and can be obtained with for instance.
pip3 install scipy on your server

1 Like