Surface-to-surface radiation in FENiCS

Hi everyone,

Recently I succesfully created a model for natural convection in a square cavity. I validated my code results with literature and it seems to work perfect. A schematic of this problem is shown in the following figure. The boundary conditions are also depicted. The left and right wall are set at a constant temperature while the top and bottom walls are considered as adiabatic walls. No-slip boundary condition is applied to all walls.
image

Now my final goal is to develop a convection-radiation model for air in this cavity. We will consider the working fluid as non-emitting and non-absorbing and thus it is radiatively non-participating. This in advance ensures that the governing equations of the pure convection problem are still valid. The only changes appear in the boundary conditions (surface-to-surface radiation). The walls are assumed to be diffuse-gray. This assumption is used often in literature too.

In equation form, the boundary conditions \frac{\partial\theta}{\partial Y}=0 will change to \frac{\partial\theta}{\partial Y}-N_rQ_r=0. So the sum of the convective and radiative heat fluxes on each of these two walls is zero. Here N_r is a dimensionless radiation-conduction number and Q_r is the net radiative heat flux along a the top or bottom surface, which is the sum of the net radiative heat flux of all elements on that surface. I want to use the following equations from a reference:

So it seems that I need to compute the radiative net heatflux for all elements which is dependent on the temperatures of all elements and transform this somehow in a Neumann boundary condition. Does anyone have experience with surface-to-surface radiation or have tips on how to implement this boundary condition. The full code which I want to transform to convection-radiation is shown at the bottom.

import fenics as fen
import matplotlib
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline

N = 40
mesh = fen.UnitSquareMesh(N,N)

P1 = fen.FiniteElement('P', mesh.ufl_cell(), 1) # Linear element
P2 = fen.VectorElement('P', mesh.ufl_cell(), 2) # Quadratic vector element

# Creating a mixed element in the order of pressure, velocity, temperature
element = fen.MixedElement([P1, P2, P1])

V = fen.FunctionSpace(mesh, element) # mixed finite element function space
q, v, s = fen.TestFunctions(V) # Test functions for the three components respectively


w = fen.Function(V) # system solution function
p, u, T = fen.split(w) # components of the solution: pressure, velocity vector and temperature

n = fen.FacetNormal(mesh)

# Define subsystems for boundary conditions 
V_p = V.sub(0)
V_u = V.sub(1)
V_T = V.sub(2)

# Define boundaries
left_wall = 'near(x[0],0)' 
right_wall = 'near(x[0],1)' 
adiabatic_walls = 'near(x[1],0) || near(x[1],1)' 
walls = left_wall + " | " + right_wall + " | " + adiabatic_walls

# Define hot and cold wall temperatures
Tc = fen.Constant(0)
Th = fen.Constant(1)

#Define boundary conditions
bcu = fen.DirichletBC(V_u, (0, 0), walls)
bcTh = fen.DirichletBC(V_T, Th, left_wall)
bcTc = fen.DirichletBC(V_T, Tc, right_wall)
bcs = [bcu, bcTh, bcTc]

w_n = fen.interpolate(
    fen.Expression(("0.", "0.", "0.", "T_h + x[0]*(T_c - T_h)"), 
                      T_h = Th, T_c = Tc,
                      element = element),
    V)

p_n, u_n, T_n = fen.split(w_n)


timestep_size = 0.001
Delta_t = fen.Constant(timestep_size)
u_t = (u - u_n)/Delta_t
T_t = (T - T_n)/Delta_t

Pr = fen.Constant(0.71)
Ra = fen.Constant(1000)

fBoussinesq = Ra*Pr*T*fen.Constant((0., -1.))

inner, dot, grad, div, sym, nabla_grad = \
    fen.inner, fen.dot, fen.grad, fen.div, fen.sym, fen.nabla_grad

mass = q*div(u)

momentum = dot(v, dot(grad(u), u)) - dot(grad(p), v) \
           + Pr*inner(nabla_grad(v), nabla_grad(u)) \
           + dot(v, u_t) + dot(v, fBoussinesq)

energy = s*T_t + dot(u, grad(T))*s + dot(grad(T), grad(s))

F = (mass + momentum + energy)*fen.dx

penalty_stabilization_parameter = 1.e-7

gamma = fen.Constant(penalty_stabilization_parameter)

F += -q*gamma*p*fen.dx

JF = fen.derivative(F, w, fen.TrialFunction(V))
problem = fen.NonlinearVariationalProblem(F, w, bcs, JF)
solver = fen.NonlinearVariationalSolver(problem)
w.vector()[:] = w_n.vector()
solver.solve()

for timestep in range(10):
    w_n.vector()[:] = w.vector()
    timestep_size *= 2.   
    Delta_t.assign(timestep_size)
    solver.solve()
    
p, u, T = w.split(deepcopy=True)