De Vahl Davis natural convection test
Click here to return back to Fluid Mechanics
Contents
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 testCode
Full examples can be found under the examples in the code repository Main Page).
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
}Full 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.
Quick comparison of explicit ACM, explicit pressure correction, and implicit methods
In below images you can see quick comparison of explicit stepping with ACM p-v coupling (EACM), explicit stepping with pressure correction p-v coupling (EPC), and fully implicit stepping with pressure correction p-v coupling (IPV). EACM is naturally the less computationally expensive, followed by EPC, and finally the IPC is much more costly due to the required sparse system solution at every step.
