De Vahl Davis natural convection test

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

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

//transport equations
        #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);