Transient Natural Convection in Partially Porous Vertical Annuli

For thermo-fluid analysis of problems involving the interaction between a layer assimilable to a porous medium and a free fluid which saturates it, two different approaches can be followed. A first one, called two domains approach, consists in solving two different equations for porous region and free fluid domain by defining properly coupling conditions for the interface. The formulation of these boundary conditions at the porous region- free fluid interface is not simple as often the different sub-regions are described by differential equations of different orders, e.g. Navier-Stokes and Darcy’s equation. Alternatively, a second approach more general can be introduced by defining a mathematical model based on a single equation suitable for application for both fluid and porous region. This model, called generalized porous medium model, is based on the concept of Representative Elementary Volume, whose size should be so that the values of the quantities of interest are independent of it; by averaging Navier-Stokes equations over this elementary unit using the volume averaging procedure introduced by Whitaker [30-31] it is possible to derive a generalized model that allows to define a single computational domain with different properties and therefore only matching conditions must be assumed at the interface, and there is no need for any empirical coefficient. Considering a uniform porous medium, under the hypothesis of Local Thermal Equilibrium, the non-dimensional form of the generalized model in cylindrical coordinates for the transient natural convection can be postulated as:

$\frac{1}{r} \frac{\partial\left(r \cdot u_{r}\right)}{\partial r}+\frac{\partial u_{z}}{\partial z}=0$(1)            

$\frac{1}{\varepsilon} \frac{\partial u_{r}}{\partial t}+\frac{1}{\varepsilon^{2}}\left(u_{r} \frac{\partial u_{r}}{\partial r}+u_{z} \frac{\partial u_{r}}{\partial z}\right)=-\frac{\partial p_{f}}{\partial r}+$$\sqrt{\frac{\operatorname{Pr}}{R a}} \frac{J}{\varepsilon}\left[\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial u_{r}}{\partial r}\right)-\frac{u_{r}}{r^{2}}+\frac{\partial^{2} u_{r}}{\partial z^{2}}\right]-\left(\frac{1}{D a} \sqrt{\frac{\operatorname{Pr}}{R a}}+\frac{F o}{D a}|\vec{u}|\right) u_{r}$(2)

$\frac{1}{\varepsilon} \frac{\partial u_{z}}{\partial t}+\frac{1}{\varepsilon^{2}}\left(u_{r} \frac{\partial u_{z}}{\partial r}+u_{z} \frac{\partial u_{z}}{\partial z}\right)=-\frac{\partial p_{f}}{\partial z}+$$\sqrt{\frac{\operatorname{Pr}}{R a}} \frac{J}{\varepsilon}\left[\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial u_{z}}{\partial r}\right)+\frac{\partial^{2} u_{z}}{\partial z^{2}}\right]-\left(\frac{1}{D a} \sqrt{\frac{\operatorname{Pr}}{R a}}+\frac{F o}{D a}\right) u_{z}+T$(3)

$\sigma \frac{\partial T}{\partial t}+u_{r} \frac{\partial T}{\partial r}+u_{z} \frac{\partial T}{\partial z}=\frac{\lambda}{\sqrt{\operatorname{Ra} \operatorname{Pr}}}\left(\frac{1}{r} \frac{\partial}{\partial r}\left(\frac{\partial(r \cdot T)}{\partial r}\right)+\frac{\partial^{2} T}{\partial y^{2}}\right)$(4) 

Here, The Boussinesq approximation is introduced to model buoyancy effects.

The scales and parameters used in the non-dimensional form are:

\(\begin{align}  & r=\frac{{{r}^{*}}}{L}{{;}^{{}}}{{^{{}}}^{{}}}T=\frac{{{T}^{*}}-T_{m}^{*}}{T_{h}^{*}-T_{c}^{*}}{{;}^{{}}}{{^{{}}}^{{}}}T_{m}^{*}=\frac{T_{h}^{*}+T_{c}^{*}}{2}{{;}^{{}}}{{^{{}}}^{{}}} \\ & Ra=\frac{g{{B}_{f}}\left( T_{h}^{*}-T_{c}^{*} \right){{L}^{3}}}{{{\nu }_{f}}{{\alpha }_{f}}}{{;}^{{}}}{{^{{}}}^{{}}}\Pr =\frac{{{\nu }_{f}}}{{{\alpha }_{f}}}; \\ & t={{t}^{*}}\frac{\sqrt{g\beta L\Delta T}}{L}{{;}^{{}}}{{^{{}}}^{{}}}p=\frac{{{p}^{*}}}{\rho g\beta L\Delta T}{{;}^{{}}}{{^{{}}}^{{}}}J=\frac{{{\mu }_{eff}}}{{{\mu }_{f}}}{{;}^{{}}}{{^{{}}}^{{}}} \\ & \lambda =\frac{{{\lambda }_{eff}}}{{{\lambda }_{f}}}{{;}^{{}}}{{^{{}}}^{{}}}Da=\frac{\kappa }{{{L}^{2}}};\sigma =\frac{\varepsilon {{\left( \rho {{c}_{p}} \right)}_{f}}+\left( 1-\varepsilon  \right){{\left( \rho {{c}_{p}} \right)}_{s}}}{{{\left( \rho {{c}_{p}} \right)}_{f}}}{{;}^{{}}}{{^{{}}}^{{}}} \\ & {{\nu }_{f}}=\frac{{{\mu }_{f}}}{{{\rho }_{f}}}{{;}^{{}}}{{^{{}}}^{{}}}{{\alpha }_{f}}=\frac{{{\lambda }_{f}}}{{{\left( \rho {{c}_{p}} \right)}_{f}}}{{;}^{{}}}{{^{{}}}^{{}}}{{u}_{i}}=\frac{u_{i}^{*}}{\sqrt{g\beta L\Delta T}} \\ \end{align}\) (5)

The Navier-Stokes equations can be recovered by just setting $\varepsilon \rightarrow 1$ and $D a \rightarrow \infty$ in the above generalized model equations, while for $\varepsilon \rightarrow 0$ and $D a \rightarrow 0$ Eq. (4) describes the pure heat conduction through a solid. As the generalized model allows to simplify the numerical solution of highly coupled problems it is particularly suitable for the analysis of real life interface problems.

In this section, the numerical results obtained for unsteady free convection in partially porous tall annuli with ARs of 4:1 and 8:1 are presented. In order to investigate the influence of porous layer position on thermo-fluid behavior of free convection in a partially porous vertical annulus, several analyses have been carried out by changing the position of the porous insert. In particular, two different configurations have been analyzed: a first one with the porous layer is on the left side of the computational domain (case A), and a second one obtained by exchanging the sub-domains (case B). The computational domain and the reference boundary conditions are represented in Figure 1. In particular, no slip conditions are assumed for all the walls of the cavity. The horizontal walls are assumed to be adiabatic, while Dirichlet temperature boundary conditions are imposed on vertical ones: in particular, the inner vertical boundary is kept at a dimensionless temperature Th=0.5 and the outer wall is at a fixed dimensionless temperature Tc=-0.5. Zero temperature and velocity fields are assumed as initial conditions. Moreover, the fluid is assumed to be incompressible and Newtonian, with a Prandtl number equal to 0.71, and the porous matrix is assumed to be isotropic and homogeneous, with a uniform porosity equal to 0.5.

figure_1.png

Figure 1. Computational domain with a detail of computational mesh and boundary conditions employed

Figure 1 also shows a detail of the computational mesh used for the 4:1 cavity: it is composed of 14, 478 elements and 7, 679 nodes, while the mesh selected for the taller annulus consists of 30, 054 elements and 15, 592 nodes. Both grids are refined near all the boundaries of the computational domain and have been chosen after a mesh sensitivity analysis. In order to analyze the transient behavior of the natural convection, five points have been identified to analyze time variation of temperature: these points are represented in Figure 1 and their coordinates are reported in Table 1.

Table 1.   Non dimensional coordinates of probe points

Notes: see Figure 1.

The numerical investigations have been carried out assuming the value of Rayleigh number at which oscillations start both in the case of free fluid in a partially porous tall rectangular enclosure [23] (3.4 ∙ 10⁶ for the 4:1 cavity and 3.4 ∙ 10⁵ for the 8:1 cavity).

Temperature contours at real time t=500 are shown in Figures 2-5; in particular, Figure 2 and 3 refer to 4:1 annulus while the others refer to the taller cavity. A comparison between the temperature fields related to case (B) with that obtained for the case (A) shows that the position of the porous layer strongly affects the thermal behavior of the cavity: in fact, for both cases the configuration (B) leads to a chaotic evolution of temperature contours for high value of Da numbers whereas for very low values of Da, the heat conduction becomes dominant in both porous and free fluid region with positive effects on heat transfer.

figure_2.png

Figure 2. Isotherms in a 4:1 partially porous annulus (r_o/r_i=2) for different Da and Ra=3.4 ∙ 10⁶ (case A)

figure_3.png

Figure 3. Isotherms in a 4:1 partially porous annulus (r_o/r_i=2) for different Da and Ra=3.4 ∙ 10⁶ (case B)

figure_4.png

figure_4.png

Figure 4. Isotherms in a 8:1 partially porous annulus (r_o/r_i=2) for different Da and Ra=3.4 ∙ 10⁵ (case A)

figure_5.png

Figure 5. Isotherms in a 8:1 partially porous annulus (r_o/r_i=2) for different Da and Ra=3.4 ∙ 10⁵ (case B)

Transient oscillations for two cases are quantitatively represented in Figures 6-7; in particular, Figure 6 refers to transient thermal behavior for different values of Da numbers for the shorter annulus while Figure 7 represents the same results for the taller cavity. The analysis of the results highlights that for high values of Da the presence of the porous layer does not damp non-periodic oscillations in the configuration B, whereas transient oscillations tend to disappear approaching steady state for the case A.

figure_6.png

Figure 6. 4:1 annulus: Temperature time evolution at five probe points (see Table 1) for case A and case B for different Da and Ra=3.4 ∙ 10⁶ (r_o/r_i=2)

figure_7.png

Figure 7. 8:1 annulus: Temperature time evolution at five probe points (see Table 1) for case A and case B for different Da and Ra=3.4 ∙ 10⁵ (r_o/r_i=2)

In order to investigate the behavior of the flow regime with the annulus radius, several analyses have been carried out by changing the radius ratio r_o/r_i in the range between 3 and 1.1.

As transient oscillations affect natural convection for high value of Darcy number, the reference value of Da=1 has been assumed for all the simulations. Figures 8 and 9 show the results obtained both for the case A and B for the two ARs.

figure_8.png

Figure 8. 4:1 annulus: Temperature time evolution at five probe points (see Table 1) for case A and case B for different r_o/r_i and Ra=3.4 ∙ 10⁶ (Da=1)

figure_9.png

Figure 9. 8:1 annulus: Temperature time evolution at five probe points (see Table 1) for case A and case B for different r_o/r_i and Ra=3.4 ∙ 10⁵ (Da=1)

In the first case, for all values of r_o/r_i, oscillations dampen and the flow is not oscillatory in time on all the considered points. As the value of radii ratio decreases, it is possible to observe an increase of the time necessary to reach the steady state. A different trend is observed for the results related to case B as transient oscillations do not disappear for all the values of radii ratios examined here. Finally, it can be noted that in both cases, the amplitude of these oscillations increases with the decrease of the value of radii ratio. It is interesting to note that transient oscillations are more evident and chaotic for the taller cavity.