De Vahl Davis natural convection test

Intro

The classical de Vahl Davis benchmark test is defined for the natural convection of the air ($\Pr =0.71$) in the square closed cavity (${{\text{A}}_{\text{R}}}=1$). The only physical free parameter of the test remains the thermal Rayleigh number. In the original paper [1] de Vahl Davis tested the problem up to the Rayleigh number ${{10}^{6}}$, however in the latter publications, the results of more intense simulations were presented with the Rayleigh number up to ${{10}^{8}}$. Lage and Bejan [2] showed that the laminar domain of the closed cavity natural convection problem is roughly below $\text{Gr1}{{\text{0}}^{9}}$. It was reported [3, 4] that the natural convection becomes unsteady for $\text{Ra}=2\cdot {{10}^5}$. Here we present a MLSM solution of the case.

\begin{equation} \text{Ra}\text{=}\,\frac{\left| \mathbf{g} \right|{{\beta }_{T}}\left( {{T}_{H}}-{{T}_{C}} \right){{\Omega }_{H}}^{3}{{\rho }^{2}}{{c}_{p}}}{\lambda \mu } \end{equation} \begin{equation} \text{Pr}=\frac{\mu {{c}_{p}}}{\lambda } \end{equation}

[1] de Vahl Davis G. Natural convection of air in a square cavity: a bench mark numerical solution. Int J Numer Meth Fl. 1983;3:249-64.

[2] Lage JL, Bejan A. The Ra-Pr domain of laminar natural convection in an enclosure heated from the side. Numer Heat Transfer. 1991;A19:21-41.

[3] Janssen RJA, Henkes RAWM. Accuracy of finite-volume disretizations for the bifurcating natural-convection flow in a square cavity. Numer Heat Transfer. 1993;B24:191-207.

[4] Nobile E. Simulation of time-dependent flow in cavities with the additive-correction multigrid method, part II: Apllications. Numer Heat Transfer. 1996;B30:341-50.

. Figure 1: Scheme of the de Vahl Davis benchmark test

Code

The snippet of the openMP parallel MLSM code for an explicit ACM method with CBS looks like: (full examples, including implicit versions, can be found under the examples in the code repository Main Page).

//Navier-Stokes
#pragma omp parallel for private(i) schedule(static)
for (i=0;i<interior.size();++i) {
    int c=interior[i];
    //Navier-Stokes
    v2[c] = v1[c] + O.dt * ( - 1/O.rho    * op.grad(P1,c)
                             + O.mu/O.rho * op.lap(v1, c)
                             -              op.grad(v1,c)*v1[c]
                             + O.g*(1 - O.beta*(T1[c] - O.T_ref))
    );
    //heat transport
    T2[c] = T1[c] + O.lam / O.rho / O.c_p * O.dt * op.lap(T1, c)
            -O.dt*O.rho*O.c_p * v1[c].transpose()*op.grad(T1,c) ;
}
//heat Neumann condition
#pragma omp parallel for private(i) schedule(static)
for (i=0;i<top.size();++i) {
    int c = top[i];
    T2[c] = op.neumann(T2, c, vec_t{0, 1}, 0.0);
}
#pragma omp parallel for private(i) schedule(static)
for (i=0;i<bottom.size();++i) {
    int c = bottom[i];
    T2[c] = op.neumann(T2, c, vec_t{0, -1}, 0.0);
}
//Mass continuity
#pragma omp parallel for private(c) schedule(static)
for (i=0; i<domain.size(); ++i) {
    P2[i] = P1[i] - O.dl * O.dt * O.rho * op.div(v2, i) +
            O.dl2 * O.dl * O.dt * O.dt * op.lap(P1, i);
}
//time step
v1.swap(v2);
P1.swap(P2);
T1.swap(T2);