Compute the Nusselt number (Boussinesq approximation)

General
1

Hello! I am simulating natural convection in a square cavity (two cases: bottom side heated, top side cooled or left side heated, right side cooled). I’ve used the non-dimensional Boussinesq approximation of the Navier Stokes eq-ns.

The questions is: how to calculate the Nusselt number for my case? Below is the variational formulation from my code (anything else is identical to the code in the Navier-Stokes tutorial of FeniCsX).

def epsilon(u):
    """Strain-rate tensor"""
    return sym(nabla_grad(u))

def sigma(u, p):
    """Stress tensor"""
    vis_c = 1.0 / np.sqrt(Gr)
    return 2 * vis_c * epsilon(u) - p * Identity(len(u))

# Crank-Nicolson average
U = 0.5 * (u_n + u)

# Step 1: Tentative velocity
F1 = dot((u - u_n) / k, v) * dx
F1 += dot(dot(u_n, nabla_grad(u_n)), v) * dx
F1 += inner(sigma(U, p_n), epsilon(v)) * dx
F1 += (+T_n) * dot(g_vec, v) * dx 

a1 = form(lhs(F1))
L1 = form(rhs(F1))

# Step 2: Pressure correction
a2 = form(dot(nabla_grad(p), nabla_grad(q)) * dx)
L2 = form(dot(nabla_grad(p_n), nabla_grad(q)) * dx - (1.0/k) * div(u_) * q * dx)

# Step 3: Velocity correction
a3 = form(dot(u, v) * dx)
L3 = form(dot(u_, v) * dx - k * dot(nabla_grad(p_ - p_n), v) * dx)

# Step 4: Temperature equation
Gr_half = np.sqrt(Gr)
T_mid = 0.5 * (T_n + T)
F4 = Gr_half * dot((T - T_n) / k, s) * dx
F4 += Gr_half * dot(dot(u_, grad(T_mid)), s) * dx
F4 += (1.0 / Pr) * dot(grad(T_mid), grad(s)) * dx

if orientation == 1:
    F4 -= dot(q_flux, s) * ds(1)      # left heat flux
    F4 += dot(Bi * T, s) * ds(2)      # right convective cooling
else:
    F4 -= dot(q_flux, s) * ds(3)      # bottom heat flux
    F4 += dot(Bi * T, s) * ds(4)      # top convective cooling

a4 = form(lhs(F4))
L4 = form(rhs(F4))