View source for De Vahl Davis natural convection test

Click here to return back to [[Fluid Mechanics]]

=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 id="fig:devahl_scheme">
[[File:image.png]].
<caption>Scheme of the de Vahl Davis benchmark test </caption>
</figure>

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

<syntaxhighlight lang="c++" line>
    v2[boundary] = vec_t{0.0, 0.0};
    T2[left] = O.T_cold;
    T2[right] = O.T_hot;
    //Time stepping
    for (int step = 0; step <= O.t_steps; ++step) {
        for (int i_count = 1; i_count < _MAX_ITER_; ++i_count) {
            // Navier Stokes
            for (auto c : interior) {
                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)));
            }

            //Speed of sound
            Range<scal_t> norm = v2.map([](const vec_t& p) { return p.norm(); });
            scal_t C = O.dl * std::max(*std::max_element(norm.begin(), norm.end()), O.v_ref);
            // Mass continuity
            Range<scal_t> div_v;
            for (auto c:all) {
                div_v[c] = op.div(v2, c);
                P2[c] = P1[c] - C * C * O.dt * O.rho * div_v[c] +
                        O.dl2 * C * C * O.dt * O.dt * op.lap(P1, c);
            }
            P1.swap(P2);
        }

        //heat transport
        for (auto c : interior) {
            T2[c] = T1[c] + O.dt * O.lam / O.rho / O.c_p * op.lap(T1, c) -
                    O.dt * v1[c].transpose() * op.grad(T1, c);
        }
        for (auto c : top) T2[c] = op.neumann(T2, c, vec_t{0, -1}, 0.0);
        for (auto c : bottom) T2[c] = op.neumann(T2, c, vec_t{0, 1}, 0.0);
    }
</syntaxhighlight>

Snippet of code for explicit pressure correction scheme. Note that the solution of heat equation is the same as in above example

<syntaxhighlight lang="c++" line>

    for (int step = 0; step <= O.t_steps; ++step) {
        time_1 = std::chrono::high_resolution_clock::now();
        // Explicit Navier-Stokes computed on whole domain, including boundaries
        // without pressure -- Fraction step
        int i;
        #pragma omp parallel for private(i) schedule(static)
        for (i = 0; i < all.size(); ++i) {
            int c = i; //interior[i];
            v_2[c] = v_1[c] + O.dt (
                                   O.mu / O.rho * op.lap(v_1, c)
                                    - op.grad(v_1, c) * v_1[c]
                                    + O.g * (1 - O.beta * (T_1[c] - O.T_ref)));
        }
        // Pressure correction
        VecXd rhs_pressure(N + 1, 0); //Note N+1, +1 stands for regularization equation
        rhs_pressure(N) = 0; // = 0 part of the regularization equation
        #pragma omp parallel for private(i) schedule(static)
        for (i = 0; i < interior.size(); ++i) {
            int c = interior[i];
            rhs_pressure(c) = O.rho / O.dt * op.div(v_2, c);
        }
        #pragma omp parallel for private(i) schedule(static)
        for (i = 0; i < boundary.size(); ++i) {
            int c = boundary[i];
            rhs_pressure(c) = O.rho / O.dt * v_2[c].dot(domain.normal(c));
        }

        VecXd solution = solver_p.solve(rhs_pressure);
        alpha = solution[N];
        VecXd P_c = solution.head(N);

        // veclocity correction due to the pressure correction
        #pragma omp parallel for private(i) schedule(static)
        for ( int i = interior) v_2[c] -=  O.dt / O.rho * op.grad(P_c, c);
        
        v_2[boundary] = 0; // force boundary conditions

</syntaxhighlight>

=Results=

Following video shows evolution of temperature and velocity magnitude for the $Ra=10^8$ case. 

[[File:DeVahlDavisRa1e8.mp4|640px]]


In below galley you can find temperature contour plots, velocity magnitude contour plots, v_max and average hot side Nusselt number convergence behavior. The reference values are from:
*[a] 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. 
*[b]  Sadat H, Couturier S. Performance and accuracy of a meshless method for laminar natural convection. Numer Heat Transfer. 2000;B37:455-67. 
*[c] Wan DC, Patnaik BSV, Wei GW. A new benchmark quality solution for the buoyancy-driven cavity by discrete singular convolution. Numer Heat Transfer. 2001;B40:199-228. 
*[d] Šarler B. A radial basis function collocation approach in computational fluid dynamics. CMES-Comp Model Eng. 2005;7:185-93. 
*[e] Kosec G, Šarler B. Solution of thermo-fluid problems by collocation with local pressure correction. Int J Numer Method H. 2008;18:868-82. 

<gallery>
File:DVD_cont_T_Ra3.png
File:DVD_cont_T_Ra4.png
File:DVD_cont_T_Ra5.png
File:DVD_cont_T_Ra6.png
File:DVD_cont_T_Ra7.png
File:DVD_cont_T_Ra8.png
File:DVD_cont_v_Ra3.png
File:DVD_cont_v_Ra4.png
File:DVD_cont_v_Ra5.png
File:DVD_cont_v_Ra6.png
File:DVD_cont_v_Ra7.png
File:DVD_cont_v_Ra8.png
File:DVD_Nu_convRa3.png
File:DVD_Nu_convRa4.png
File:DVD_Nu_convRa5.png
File:DVD_Nu_convRa6.png
File:DVD_Nu_convRa7.png
File:DVD_Nu_convRa8.png
File:DVD_vMax_convRa3.png
File:DVD_vMax_convRa4.png
File:DVD_vMax_convRa5.png
File:DVD_vMax_convRa6.png
File:DVD_vMax_convRa7.png
File:DVD_vMax_convRa8.png
</gallery>