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.

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