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

From Medusa: Coordinate Free Mehless Method implementation
Jump to: navigation, search
Line 14: Line 14:
 
<syntaxhighlight lang="c++" line>
 
<syntaxhighlight lang="c++" line>
 
//transport equations
 
//transport equations
        #pragma omp parallel for private(i) schedule(static)
+
  //transport equations
        for (i=0;i<interior.size();++i) {
+
#pragma omp parallel for private(i) schedule(static)
            int c=interior[i];
+
for (i=0;i<interior.size();++i) {
            //Navier-Stokes
+
    int c=interior[i];
            v2[c] = v1[c] + O.dt * ( - 1/O.rho    * op.grad(P1,c)
+
    //Navier-Stokes
                                    + O.mu/O.rho * op.lap(v1, c)
+
    v2[c] = v1[c] + O.dt * ( - 1/O.rho    * op.grad(P1,c)
                                    -              op.grad(v1,c)*v1[c]
+
                            + O.mu/O.rho * op.lap(v1, c)
                                    + O.g*(1 - O.beta*(T1[c] - O.T_ref))
+
                            -              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)
+
    //heat transport
                    -O.dt*O.rho*O.c_p * v1[c].transpose()*op.grad(T1,c) ;
+
    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)
+
//heat Neumann condition
        for (i=0;i<top.size();++i) {
+
#pragma omp parallel for private(i) schedule(static)
            int c = top[i];
+
for (i=0;i<top.size();++i) {
            T2[c] = op.neumann(T2, c, vec_t{0, 1}, 0.0);
+
    int c = top[i];
        }
+
    T2[c] = op.neumann(T2, c, vec_t{0, 1}, 0.0);
 +
}
  
        #pragma omp parallel for private(i) schedule(static)
+
#pragma omp parallel for private(i) schedule(static)
        for (i=0;i<bottom.size();++i) {
+
for (i=0;i<bottom.size();++i) {
            int c = bottom[i];
+
    int c = bottom[i];
            T2[c] = op.neumann(T2, c, vec_t{0, -1}, 0.0);
+
    T2[c] = op.neumann(T2, c, vec_t{0, -1}, 0.0);
        }
+
}
  
        //Mass continuity
+
//Mass continuity
        #pragma omp parallel for private(c) schedule(static)
+
#pragma omp parallel for private(c) schedule(static)
        for (i=0; i<domain.size(); ++i) {
+
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);
 
  
 +
    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);
 
</syntaxhighlight>
 
</syntaxhighlight>

Revision as of 14:21, 22 October 2017

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.

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

 1 //transport equations
 2   //transport equations
 3 #pragma omp parallel for private(i) schedule(static)
 4 for (i=0;i<interior.size();++i) {
 5     int c=interior[i];
 6     //Navier-Stokes
 7     v2[c] = v1[c] + O.dt * ( - 1/O.rho    * op.grad(P1,c)
 8                              + O.mu/O.rho * op.lap(v1, c)
 9                              -              op.grad(v1,c)*v1[c]
10                              + O.g*(1 - O.beta*(T1[c] - O.T_ref))
11     );
12     //heat transport
13     T2[c] = T1[c] + O.lam / O.rho / O.c_p * O.dt * op.lap(T1, c)
14             -O.dt*O.rho*O.c_p * v1[c].transpose()*op.grad(T1,c) ;
15 }
16 //heat Neumann condition
17 #pragma omp parallel for private(i) schedule(static)
18 for (i=0;i<top.size();++i) {
19     int c = top[i];
20     T2[c] = op.neumann(T2, c, vec_t{0, 1}, 0.0);
21 }
22 
23 #pragma omp parallel for private(i) schedule(static)
24 for (i=0;i<bottom.size();++i) {
25     int c = bottom[i];
26     T2[c] = op.neumann(T2, c, vec_t{0, -1}, 0.0);
27 }
28 
29 //Mass continuity
30 #pragma omp parallel for private(c) schedule(static)
31 for (i=0; i<domain.size(); ++i) {
32 
33     P2[i] = P1[i] - O.dl * O.dt * O.rho * op.div(v2, i) +
34             O.dl2 * O.dl * O.dt * O.dt * op.lap(P1, i);
35 }
36 //time step
37 v1.swap(v2);
38 P1.swap(P2);
39 T1.swap(T2);