Solving for two interconducting fluids in FEniCS

Question

I'm trying to model in FEniCS a steady-state situation in which a coolant fluid flows through a nuclear fluid. There is finite rate of conduction between the two fluids. The simulation should output the steady-state temperature profile for each fluid.

The issue I've run into is that I don't know how to get FEniCS to solve for u_coolant and u_fuel concurrently. I can solve for one after the other, but that isn't sufficient.

How do I solve for the two fluid temperature distributions concurrently in FEniCS? Please let me know if I can phrase this question better, or if there's additional information I should provide.

Equations:
$Q_{fuel-gen} = q''' V$
$Q_{fuel-conv} = h A(u_{fuel} - u_{coolant})$
$Q_{fuel-gen} = Q_{fuel-conv}$

$Q_{coolant-conv} = h A(u_{fuel} - u_{coolant})$
$Q_{coolant-adv} = C \dot{m} u_{coolant}$
$Q_{coolant-conv} = Q_{coolant-adv}$

Where:
$Q_{fuel-gen}$ is the rate at which heat is generated by the fuel ($W$).
$q'''$ is the volumetric heat generation rate ($W/m^3$).
$V$ is the volume of a differential shell ($m^3$).
$Q_{fuel-conv}$ is the rate at which heat is lost from the fuel via convection to the coolant ($W$).
$h$ is the convective heat transfer coefficient ($W/m^2)$.
$A$ is the area of coolant exposed to the fuel within a differential shell ($m^2$)
$Q_{coolant-conv}$ is the rate at which the coolant gains heat due to convection against the fuel ($W$).
$Q_{coolant-adv}$ is the rate at which the coolant within a differential shell loses heat due to advection ($W$).
$C$ is the specific heat of the coolant ($J/kg K$).
$\dot{m}$ is flow rate of coolant ($kg/s$).
$u_{fuel}$ is a function describing the temperature of the fuel as it varies with position $x$ ($K$).
$u_{coolant}$ is a function describing the temperature of the coolant as it varies with position $x$ ($K$).

q_genU= gen_rate*V_shell #heat deposition rate in the fuel
q_convU = h_conv*A_bubbles*(uU-uH)*vU #convection rate from fuel to coolant
LU = q_depoU*vU*dx + gU*vU*ds(1) + coolantBC_l
aU = (q_condU + q_convU)*dx + coolantBC_bl

q_convH = h_conv*A_bubbles*(uU-uH)
q_advH = C_H*cool_flow*uH
LH = gH*vH*ds
aH = (q_convH-q_advH)*vH*dx

uU = Function(V)
uH = Function(V)
solve(aU == LU, uU, bcs_U)
solve(aH == LH, uH, bcs_H)

Answer 1

You need to define your function space of your trial and test functions as the product space of the corresponding elements.

# CG element, degree 1
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)

# function space for all fields (here two)
W = FunctionSpace(mesh, P1*P1)

# trial and test functions
u_coolant, u_fuel  = TrialFunction(W)
p, q  = TestFunctions(W)

The declaration of W means $W=\{(u,v) \textrm{ such that } u\in \textrm{P1}, v \in \textrm{P1}\}$.

The return values of TrialFunction and TestFunction are split back into the two independent fields.

For a detailed example see here.