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.

$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}$

$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)
  • 2
    $\begingroup$ Please provide the differential equations and what you have already done. $\endgroup$
    – Paddy
    Commented Jun 21, 2021 at 18:31

1 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.