Using a temperature dependent function in the heat equation

variational formulation
1

Hi everyone,

I am currently trying to model the cooling of a cylinder via radiation as a function of time using the heat equation, where the emissivity of the material epsilon_eff is dependent on the temperature. However, epsilon_eff seems not to be updating for each time step. The basic code is below:

from fenics import *
from dolfin import *

# Define the domain
r1 = 1.0   # inner radius
r2 = 2.0   # outer radius
L = 3.0    # length
N = 20     # number of mesh elements

# Create the mesh
mesh = CylinderMesh(Point(0, -L/2, 0), Point(0, L/2, 0), r1, r2, N)

# Define the function space
V = FunctionSpace(mesh, 'P', 1)

# emissivity function
emiss_in = [0.03,0.042,0.054,0.062,0.080,0.091,0.109,0.116,0.135,0.194,0.213,0.232,0.251]
temp_emiss = [300,400,500,600,700,800,900,1000,1100,1600,1800,2000,2200]
emiss_f = UnivariateSpline(temp_emiss, emiss_in, k = 2, ext = 3)

# Define the initial temperature
u0 = Expression('1000', degree=1)

# Define the boundary conditions
bc = DirichletBC(V, u0, 'on_boundary')

# Define the thermal conductivity and radiation coefficient
k = Constant(1.0)
sigma = Constant(5.67e-8)


# Define the variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0.0)

emiss_eff_array = emiss_f(u.vector().get_local())
emiss_space_eff = FunctionSpace(mesh, 'CG',1)
epsilon_eff = Function(emiss_space_eff)
epsilon_eff.vector()[:] = emiss_eff_array

a = k * inner(grad(u), grad(v)) * dx
L = f * v * dx - epsilon_eff * sigma * u**4 * v * ds

# Define the time-stepping method
T = 100.0    # total simulation time
dt = 0.1     # time step
t = 0.0      # initial time
u = Function(V)

# Create the file for storing the solution
file = File('cooling.pvd')

# Time-stepping loop
while t < T:
    print(t)
    t += dt
    solve(a == L, u, bc)
    print(max(u.vector()))
    print(max(epsilon_eff.vector()))
    file << (u, t)

I think that the issue maybe with the line, epsilon_eff.vector()[:] = emiss_eff_array but I am unsure about how this can be fixed.

Thanks in advance

2

You would have to call

emiss_eff_array = emiss_f(u.vector().get_local())
epsilon_eff.vector()[:] = emiss_eff_array

inside the loop for it to update.

DOLFIN cannot track changes to the vector (as in assigning data directly) and create a dependency tree. Thus you need to update this manually whenever u changes.

3

I see. When I copy said lines into the time loop, I get the following error:

RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 2 iterations (PETSc reason DIVERGED_INDEFINITE_PC, residual norm ||r|| = 2.786253e+00).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset:  ba376b6aebd7a9bc089be46b50bdd9f5c548fb91
*** -------------------------------------------------------------------------
4

Then you should inspect the values you set to epsilon_eff.vector() (for instance write epsilon_eff to file, and check if the values are sane.

5

The values of epsilon_eff are sane. I placed the lines after solve(a == L, u, bc) so I wonder if that affects it.

6

They should of course affect the solver. You wanted the epsilon_eff to update every time u changes, which is during solve. Thus after the solving the problem, you need to update epsilon_eff, such that this is reflected in the next solve.