Difference between revisions of "De Vahl Davis natural convection test"

From Medusa: Coordinate Free Mehless Method implementation
Jump to: navigation, search
 
(88 intermediate revisions by 3 users not shown)
Line 1: Line 1:
−
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.
+
Click here to return back to [[Fluid Mechanics]]
 +
 
 +
{{Box-round|title= Related papers |
 +
[https://e6.ijs.si/ParallelAndDistributedSystems/publications/992507.pdf G. Kosec, B. Šarler; Solution of thermo-fluid problems by collocation with local pressure correction, International journal of numerical methods for heat & fluid flow, vol.18, 2008]
 +
 
 +
[https://e6.ijs.si/ParallelAndDistributedSystems/publications/1905659.pdf G. Kosec, M. Založnik, B. Šarler, H. Combeau; A meshless approach towards solution of macrosegregation phenomena, Computers, materials & continua : CMC, vol. 22, 2011 ]
 +
 
 +
[http://comms.ijs.si/~gkosec/data/papers/27339815.pdf G. Kosec, M. Depolli, A. Rashkovska, R. Trobec; Super linear speedup in a local parallel meshless solution of thermo-fluid problem, Computers & Structures, vol. 133, 2014]
 +
 
 +
[https://e6.ijs.si/ParallelAndDistributedSystems/publications/3218939.pdf G. Kosec, B. Šarler; Simulation of macrosegregation with mesosegregates in binary metallic casts by a meshless method, Engineering analysis with boundary elements]
 +
 
 +
[https://e6.ijs.si/ParallelAndDistributedSystems/publications/32388135.pdf M. Maksić, V. Djurica, A. Souvent, J. Slak, M. Depolli, G. Kosec; Cooling of overhead power lines due to the natural convection, International journal of electrical power & energy systems, 2019]
 +
}}
 +
 
 +
=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.  
 
[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.  
 
[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.  
 
[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.
 
[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.
  
−
Problem is schematically presented as
+
<figure id="fig:devahl_scheme">
 
[[File:image.png]].
 
[[File:image.png]].
 +
<caption>Scheme of the de Vahl Davis benchmark test </caption>
 +
</figure>
 +
 +
=Code=
 +
Full examples can be found under the examples in the [https://gitlab.com/e62Lab/medusa code repository].
 +
 +
==Explicit ACM method with CBS looks==
 +
<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>
 +
 +
== Explicit pressure correction ==
 +
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) {
 +
        // Explicit Navier-Stokes computed on whole domain, including boundaries
 +
        // without pressure
 +
        for (int c:all) {
 +
            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
 +
        for (int i:interior) rhs_pressure(c) = O.rho / O.dt * op.div(v_2, c);
 +
        for (int i: boundary) 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);
 +
        for ( int i = interior) v_2[c] -=  O.dt / O.rho * op.grad(P_c, c);
 +
        v_2[boundary] = 0; // force boundary conditions
 +
        //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>
 +
 +
== Implicit ==
 +
<syntaxhighlight lang="c++" line>
 +
for (int step = 0; step <= O.t_steps_i; ++step) {
 +
        time_1 = std::chrono::high_resolution_clock::now();
 +
        // NAVIER STOKES
 +
        M_velocity = mat_t(2 * N, 2 * N);
 +
        // system
 +
        M_velocity.reserve(Range<int>(2 * N, O.n));
 +
        for (int i : interior) {
 +
            op.valuevec(M_velocity, i, 1 / O.dt_i);
 +
            op.lapvec(M_velocity, i, -O.mu / O.rho);
 +
            op.gradvec(M_velocity, i, v_1[i]);
 +
        }
 +
        for (int i : boundary) op.valuevec(M_velocity, i, 1); //sets velocity to 0
 +
 +
        M_velocity.makeCompressed();
 +
        solver_v.compute(M_velocity);
 +
        // solution
 +
        Range <vec_t> rhs_vec(domain.size(), 0);
 +
        for (int i : interior) {
 +
            rhs_vec[i] = -1 / O.rho * op.grad(P, i) +
 +
                        v_1[i] / O.dt_i
 +
                        + O.g * (1 - O.beta * (T[i] - O.T_ref));
 +
        }
 +
        // for (int i:top) rhs_vec[i] = vec_t{0,1};
 +
 +
        v_2 = reshape<2>(solver_v.solveWithGuess(reshape(rhs_vec), reshape(v_1)));
 +
        // END OF NAVIER STOKES
 +
 +
        // 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
 +
        double alpha;
 +
        for (int i : interior) rhs_pressure(i) = O.rho / O.dt_i * op.div(v_2, i);
 +
        for (int i : boundary) rhs_pressure(i) = O.rho / O.dt * v_2[i].dot(domain.normal(i));
 +
 +
        VecXd solution = solver_p.solve(rhs_pressure);
 +
        alpha = solution[N];
 +
        VecXd P_c = solution.head(N);
 +
        // apply velocity correction
 +
        for (int i : interior) {
 +
            v_2[i] -= O.dl * O.dt_i / O.rho * op.grad(P_c, i);
 +
        }
 +
        P += O.dl * P_c;
 +
        // enforce velocity BC
 +
        // v_2[boundary] = 0;
 +
        // END OF PRESSURE CORRECTION
 +
 +
        // HEAT TRANSPORT
 +
        M_temperature = mat_t(N, N);
 +
        Range<int> per_row(N, O.n);
 +
        M_temperature.reserve(per_row);
 +
        // outer boundary dirichlet BC
 +
        for (int i : top) op.neumann_implicit(M_temperature, i, domain.normal(i), 1);
 +
        for (int i : bottom) op.neumann_implicit(M_temperature, i, domain.normal(i), 1);
 +
        for (int i : left) op.value(M_temperature, i, 1.0);
 +
        for (int i : right) op.value(M_temperature, i, 1.0);
 +
        // heat transport in air
 +
        for (int i : interior) {
 +
            op.value(M_temperature, i, 1.0);                          // time dependency
 +
            op.lap(M_temperature, i, -O.dt_i * O.lam / O.rho / O.c_p);  //laplace in interior
 +
            op.grad(M_temperature, i, O.dt * v_2[i]);
 +
        }
 +
        M_temperature.makeCompressed();
 +
        solver_T.compute(M_temperature);
 +
 +
        VectorXd rhs = VectorXd::Zero(N);
 +
        for (int i : interior) rhs(i) = T(i);
 +
        for (int i : top) rhs(i) = 0;
 +
        for (int i : bottom) rhs(i) = 0;
 +
        for (int i : left) rhs(i) = O.T_hot;
 +
        for (int i : right) rhs(i) = O.T_cold;
 +
        T = solver_T.solveWithGuess(rhs, T);
 +
        // END OF HEAT TRANSPORT
 +
        v_1.swap(v_2);
 +
}
 +
</syntaxhighlight>
 +
 +
=Results=
 +
== Comparison of MLSM solution with reference data  ==
 +
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>

Latest revision as of 19:54, 22 October 2022

Click here to return back to Fluid Mechanics

edit 

Related papers

G. Kosec, B. Šarler; Solution of thermo-fluid problems by collocation with local pressure correction, International journal of numerical methods for heat & fluid flow, vol.18, 2008

G. Kosec, M. Založnik, B. Šarler, H. Combeau; A meshless approach towards solution of macrosegregation phenomena, Computers, materials & continua : CMC, vol. 22, 2011

G. Kosec, M. Depolli, A. Rashkovska, R. Trobec; Super linear speedup in a local parallel meshless solution of thermo-fluid problem, Computers & Structures, vol. 133, 2014

G. Kosec, B. Šarler; Simulation of macrosegregation with mesosegregates in binary metallic casts by a meshless method, Engineering analysis with boundary elements

M. Maksić, V. Djurica, A. Souvent, J. Slak, M. Depolli, G. Kosec; Cooling of overhead power lines due to the natural convection, International journal of electrical power & energy systems, 2019


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.

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

Code

Full examples can be found under the examples in the code repository.

Explicit ACM method with CBS looks

    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);
    }

Explicit pressure correction

The solution of heat equation is the same as in above example 

    for (int step = 0; step <= O.t_steps; ++step) {
        // Explicit Navier-Stokes computed on whole domain, including boundaries
        // without pressure
        for (int c:all) {
            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
        for (int i:interior) rhs_pressure(c) = O.rho / O.dt * op.div(v_2, c);
        for (int i: boundary) 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);
        for ( int i = interior) v_2[c] -=  O.dt / O.rho * op.grad(P_c, c);
        v_2[boundary] = 0; // force boundary conditions
        //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);
}

Implicit

for (int step = 0; step <= O.t_steps_i; ++step) {
        time_1 = std::chrono::high_resolution_clock::now();
        // NAVIER STOKES
        M_velocity = mat_t(2 * N, 2 * N);
        // system
        M_velocity.reserve(Range<int>(2 * N, O.n));
        for (int i : interior) {
            op.valuevec(M_velocity, i, 1 / O.dt_i);
            op.lapvec(M_velocity, i, -O.mu / O.rho);
            op.gradvec(M_velocity, i, v_1[i]);
        }
        for (int i : boundary) op.valuevec(M_velocity, i, 1); //sets velocity to 0

        M_velocity.makeCompressed();
        solver_v.compute(M_velocity);
        // solution
        Range <vec_t> rhs_vec(domain.size(), 0);
        for (int i : interior) {
            rhs_vec[i] = -1 / O.rho * op.grad(P, i) +
                         v_1[i] / O.dt_i
                         + O.g * (1 - O.beta * (T[i] - O.T_ref));
        }
        // for (int i:top) rhs_vec[i] = vec_t{0,1};

        v_2 = reshape<2>(solver_v.solveWithGuess(reshape(rhs_vec), reshape(v_1)));
        // END OF NAVIER STOKES

        // 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
        double alpha;
        for (int i : interior) rhs_pressure(i) = O.rho / O.dt_i * op.div(v_2, i);
        for (int i : boundary) rhs_pressure(i) = O.rho / O.dt * v_2[i].dot(domain.normal(i));

        VecXd solution = solver_p.solve(rhs_pressure);
        alpha = solution[N];
        VecXd P_c = solution.head(N);
        // apply velocity correction
        for (int i : interior) {
            v_2[i] -= O.dl * O.dt_i / O.rho * op.grad(P_c, i);
        }
        P += O.dl * P_c;
        // enforce velocity BC
        // v_2[boundary] = 0;
        // END OF PRESSURE CORRECTION

        // HEAT TRANSPORT
        M_temperature = mat_t(N, N);
        Range<int> per_row(N, O.n);
        M_temperature.reserve(per_row);
        // outer boundary dirichlet BC
        for (int i : top) op.neumann_implicit(M_temperature, i, domain.normal(i), 1);
        for (int i : bottom) op.neumann_implicit(M_temperature, i, domain.normal(i), 1);
        for (int i : left) op.value(M_temperature, i, 1.0);
        for (int i : right) op.value(M_temperature, i, 1.0);
        // heat transport in air
        for (int i : interior) {
            op.value(M_temperature, i, 1.0);                          // time dependency
            op.lap(M_temperature, i, -O.dt_i * O.lam / O.rho / O.c_p);  //laplace in interior
            op.grad(M_temperature, i, O.dt * v_2[i]);
        }
        M_temperature.makeCompressed();
        solver_T.compute(M_temperature);

        VectorXd rhs = VectorXd::Zero(N);
        for (int i : interior) rhs(i) = T(i);
        for (int i : top) rhs(i) = 0;
        for (int i : bottom) rhs(i) = 0;
        for (int i : left) rhs(i) = O.T_hot;
        for (int i : right) rhs(i) = O.T_cold;
        T = solver_T.solveWithGuess(rhs, T);
        // END OF HEAT TRANSPORT
        v_1.swap(v_2);
}

Results

Comparison of MLSM solution with reference data

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


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.
  • DVD cont T Ra3.png
  • DVD cont T Ra4.png
  • DVD cont T Ra5.png
  • DVD cont T Ra6.png
  • DVD cont T Ra7.png
  • DVD cont T Ra8.png
  • DVD cont v Ra3.png
  • DVD cont v Ra4.png
  • DVD cont v Ra5.png
  • DVD cont v Ra6.png
  • DVD cont v Ra7.png
  • DVD cont v Ra8.png
  • DVD Nu convRa3.png
  • DVD Nu convRa4.png
  • DVD Nu convRa5.png
  • DVD Nu convRa6.png
  • DVD Nu convRa7.png
  • DVD Nu convRa8.png
  • DVD vMax convRa3.png
  • DVD vMax convRa4.png
  • DVD vMax convRa5.png
  • DVD vMax convRa6.png
  • DVD vMax convRa7.png
  • DVD vMax convRa8.png
Retrieved from "https://e6.ijs.si/medusa/wiki/index.php?title=De_Vahl_Davis_natural_convection_test&oldid=3457"