Natural convection in a square cavity: heat flux on one wall and convective cooling on the opposite, Boussinesq approximation

variational formulation
1

Hello! This is my variational formulation for the Boussinesq approximation (Navier-Stokes). Though I am getting somewhat a realistic simulation - I wanted to make sure whether I have correctly defined the boundary conditions in the weak form. Please I need your your help!

F4 -= dot(q_flux, s) * ds(hot_id)             # left or bottom heat flux
F4 += Bi * (T - T_ref) * s * ds(cool_id)      # right or top convective cooling

The full formulation here:

#=============================== VARIATIONAL FORMULATION ===============================#
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

F4 -= dot(q_flux, s) * ds(hot_id)             # left or bottom heat flux
F4 += Bi * (T - T_ref) * s * ds(cool_id)      # right or top convective cooling

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