© O. Hireche et al., published by EDP Sciences 2020

1 Introduction

In thermoacoustic devices, dynamic and thermal interactions between acoustic oscillations of a working gas and solid walls generate either heat pumping for refrigerating purposes or mechanical work for engine operating purposes [1]. These devices represent an ecological alternative for industrial applications but they suffer from poor performance. This is generally caused by the existence of several nonlinear phenomena interfering with the process of energy transfer and penalizing the global performance of the device.

The energy conversion takes place in the vicinity of solid-fluid interfaces and in order to increase the exchange surface, a porous medium is used for this conversion. The porous medium is a significant component of thermoacoustic devices and its characteristics greatly influence the device performance. The pores are usually of various shapes and sizes. In order to maintain the temperature distribution along the stack stable, thermal conductivity of the solid matrix should be lower than that of the working gas. Arnott et al. [2] developed a general thermoacoustic linear theory for the elements of thermoacoustic devices that have a microstructure and compared their performances. The porous media considered were solid stacks of plates or solids with many capillary tubes of circular, square or triangular geometry. They concluded that heat and work flows are approximately 10 percent greater for the parallel plate stack geometry than for other pore geometries. This model is still used in DeltaEC software [3], which is the reference code for computing a complete thermoacoustic system.

In thermoacoustic systems several nonlinear phenomena can take place, reducing their performance, such as: jet streaming, acoustic streaming, end effects and natural convection. These effects are still under investigation. Among them, acoustic streaming has already been the subject of many studies during the last decades [4–6]. Natural convection is a phenomenon that has mostly been neglected until now in thermoacoustic applications, although it does exist because of the presence along the resonator of heat exchangers at different temperatures: two heat exchangers are generally placed at each end of the porous medium, in order to maintain a temperature gradient along the porous medium in thermoacoustic engines, or to benefit from heat pumping effects. Therefore, if the device is horizontal, there will be horizontal temperature gradients along the porous medium saturated with gas, as well as along the two cavities filled with gas on both sides of the porous medium. Natural convection flows can then develop as in differentially heated enclosures.

Natural convection in fluids and porous cavities has been treated extensively in the literature (see for example [7,8]). A more recent study has been conducted by Weisman et al. [9] in a geometry of a thermoacoustic engine. The goal was to predict the temperature distribution in a 2D differentially heated horizontal cavity, filled with a stack of plates which constitutes an anisotropic porous medium. This study showed that a darcean porous model gives correct orders of velocity scales and allows for a qualitative description of the solutions. Some recent experimental studies of natural convection in thermoacoustic configurations show the influence of natural convection on Rayleigh streaming [10,11].

In the present work, natural convection flows are investigated in a simplified geometry of a standing-wave thermoacoustic engine, in order to assess the flow and temperature distributions and their order of magnitude. A 3D numerical study of natural convection flow is performed using a finite volume code for solving Navier-Stokes equations under Boussinesq approximation in the fluid domain. In the porous medium, the equations are modified using a Darcean term, according to Darcy-Brinkman formulation [7,8]. The Darcean term becomes dominant as permeability decreases. The influence of the porous medium’s physical characteristics (permeability, thermal conductivity, anisotropy) on the flow and temperature fields is investigated. The results are then analysed in order to explore thermo-physical characteristics adequate for thermoacoustic devices [12].

2 Model description

2.1 Geometry and boundary conditions

A simplifiedgeometry for a standing-wave thermoacoustic engine is introduced: a horizontal cylindrical tube of length L and inner radius R, closed at both ends, is filled with air. A cylindrical porous medium of length L_p and radius R is inserted in the tube. In a standing-wave thermoacoustic engine, the fundamental mode is of acoustic wavelength λ = 2L. The porous medium length L_p is usually small compared to the wavelength, such that L_p ≪ L. The computational domain is shown in Figure 1.

Attached to the left end of the porous medium is a hot heat exchanger (positioned at a distance L_L from the left end of the tube), which can be a coiled resistance modeled here as solid rings of width e_h, maintained at given (hot) temperature T_h. Usually the “hot” part of the tube (left of coiled resistance) is insulated and therefore in the present model, all walls in this part are considered adiabatic. There is no cold heat exchanger in the present model, although there usually is one in experiments. Here the “cold” temperature is the ambient air temperature outside the cylinder, noted T_c. In the present model this cold temperature is imposed on all walls at the right of the porous medium. Finally, a linear temperature profile is imposed on the tube wall around the porous medium. Thermal boundary conditions are shown in Figure 1 (bottom). At all solid-fluid interfaces, we assume continuity of temperature and heat flux as well as no slip conditions. The porous-fluid interface is not treated as such since the equations are solved in both porous and fluid domains (see next section), a phase variable allowing to select the appropriate equation. Continuity of velocity, temperature, heat flux and stresses are implicitely ensured with the use of a staggered mesh.

The porous medium is usually positioned inside the tube at about 1/8 of the total length from the left end. However, in this natural convection study, with the imposed boundary conditions, we have shown that the gas in the right half of the tube is almost stagnant at cold temperature. Therefore it is not necessary to extend the computational domain too much. In the following simulations the porous medium will be positioned at about 1/4 of the total tube length.

The porous medium, consisting of its solid matrix and the fluid pores, is modeled as an equivalent homogeneous porous medium, with porosity ɛ, permeability tensor K, thermal conductivitytensor k_p. Also, it is common in thermoacoustics to use porous media composed of a solid (ceramic) matrix with horizontal cylindrical pores. In that case the porous medium is strongly anisotropic and this possibility will also be considered.

Fig. 1 Computational domain (top) and thermal boundary conditions (bottom).

2.2 Governing equations

In the fluid domain the Navier-Stokes equations are classically written under Boussinesq assumptions. The reference gas density ρ, thermal conductivity k, heat capacity c_p, dynamic viscosity μ and volume expansion coefficient at constant pressure β = 1∕T for an idealgas are constant and equal to their value for cold temperature T_c.

In the porous medium, the momentum equations are modified using a Darcy term. In the resulting Darcy-Brinkman equation, two simplifications are made related to the fact that in the application, porosity ɛ is rather large (around 0.8): the 1∕ɛ factor is omitted in the unsteady term (which does not affect the steady solution) and the effective viscosity is close to the fluid viscosity μ_eff = μ∕ɛ ≃ μ [7]. The advection and Forchheimer terms are also omitted in the formulation since the pore Reynolds number Re_p is assumed to be very small (Re_p ≪ 1). This assumption will be verified on results. The energy equation in the porous medium is deduced from the fluid and solid energy equations averaged over an elementary volume [7]. The value of ${(\rho c_p)}_p$ in the porous medium, given by ${(\rho c_p)}_p=\varepsilon (\rho c_p)+(1-\varepsilon ){(\rho c_p)}_s$, does not affect steady state. A phase variable X is introduced to identify the nature of the medium (X = 1 for gas, X = 0 for the porous medium). Therefore the equations can be written in a unified form as: $$\begin{array}{lll}\hfill & \hfill & \nabla .V=0\hfill \\ \hfill & \hfill & \frac{\partial (\rho V)}{\partial t}+X\nabla .(\rho V\otimes V)=-\nabla p+\mu \text{Δ}V\hfill \\ \hfill & \hfill & +\rho g\beta (T-T_c)e_x-(1-X)\mu K^{-1}V\hfill \\ \hfill & \hfill & \left((1-X){(\rho c_p)}_p+X\rho c_p\right)\frac{\partial T}{\partial t}+\rho c_{p}V.\nabla T\hfill \\ \hfill & \hfill & =\nabla .(XkI+(1-X)k_p)\nabla T.\hfill \end{array}$$

In the above equations V is the velocity vector (equal to the seepage velocity in the porous medium), p is the nonhydrostatic pressure, T is the temperature, t is time, e_x is the unit vector in the ascending vertical direction x, g is the acceleration of gravity and I the identity tensor.

2.3 Numerical solution

Computations are performed with the code SUNFLUIDH, which has been developped at LIMSI for simulations of unsteady incompressible or low Mach number flows.

The previous governing equations are solved using a finite volume approach with a 2nd order discretization in space and time. An implicit treatment of the viscous and thermal conductive terms ensures better numerical stability in regard to the time step which only depends on the CFL criterion. The equation set is solved by means of an ADI method [13]. An incremental projection method is used in order to ensure the divergence-free velocity field [14,15]. This implies to solve a Poisson’s equation at each time step which is carried out with a multigrid-SOR algorithm [16]. The solution corresponds to the time-variation of pressure from which the correction term is estimated in order to get the divergence-free velocity field.

The mesh is body-fitted irregular, with 64 points regularly spaced along the radial direction, 48 points regularly spaced along the azimutal direction and 256 points placed along the axis, with refined spacing near the heat exchanger and in the porous medium. The minimum and maximum mesh sizes in the axial direction are Δz_min = 3.9 × 10⁻⁴m and Δz_max = 4.2 × 10⁻³m. In the radial direction, Δr = 3.9 × 10⁻⁴ m. This spatial resolution is adequate for the given problem configurations.

The time step is adjusted throughout the simulation in order to satisfy a CFL<0.3 criterion.

3 Results

3.1 Reference configuration

The geometrical and thermal characteristics chosen for all cases are shown in Table 1. The thermophysical characteristics for air (fixed) and for the porous medium (for the reference case) are shown in Table 2. These parameters correspond (mostly) to some existing experiments in thermoacoustics [11]. The porous medium is considered to be isotropic and the permeability and thermal conductivity tensors are reduced to scalar values (K = KI and k_p = k_pI). The permeability value 10⁻⁶ m² constitutes a high limit to practical values found in the literature [12]. If the porous medium is made of horizontal capillaries of square cross-section, the permeability is K = a²ɛ∕8 corresponding to a hydraulic radius a equal to 3 mm. The equivalent thermal conductivity k_p is calculated from that of the solid matrix constituting the porous medium (k_s), the gas thermal conductivity and the porosity: k_p = ɛk + (1 − ɛ)k_s.

The nondimensional parameter controlling natural convection is the Rayleigh number, which is classically calculated as $Ra=\frac{\rho g\beta \text{Δ}T{H}^3}{\mu \alpha }$, with H a reference height, ΔT the imposed temperature difference and $\alpha =\frac{k}{\rho c_p}$ the fluid thermal diffusivity [7]. The associated vertical velocity scale is $U_{\text{ref}}=\frac{\alpha }{H}R{a}^{1/2}$. For the reference configuration, with H = R and ΔT = 207 K, we find that Ra = 2 × 10⁵ and U_ref = 0.48 m/s. The solution should thus be stationary. However as we will see in the results section, several convection flows can be distinguished within the domain, and this choice of Rayleigh number may not be representative of the overall flow. Note that for this value of ΔT the ratio ΔT∕(2T_ref) = 0.3 is rather large and falls outside of the validity range for the Boussinesq approximation (for air). Non Boussinesq simulations would require the use of a low Mach number model which is harder to implement. Since the flow is steady, expected differences should be small [17]. Therefore Boussinesq approximation is used throughout this paper.

In a porous cavity, the flow slows down and therefore a modified Rayleigh number is introduced, $R{a}_p=\frac{\rho g\beta \text{Δ}THK}{\mu {\alpha }_p}$, with ${\alpha }_p=\frac{k_p}{\rho c_p}$ the porous medium equivalent thermal diffusivity [7]. The associated vertical velocity scale is $U_{\text{ref}p}=\frac{{\alpha }_p}{H}R{a}_p$. For the reference configuration, using the same values for ΔT and H =R, this yields Ra_p = 6.71, U_{ref p} = 0.35 m/s and we expect a stationary solution of the natural convection problem. Again these reference values are just indicators since the problem here is a mixed fluid and porous domain with prescribed temperature differences unevenly distributed.

3.2 Numerical results for the reference configuration

The flow and temperature fields are first described in the reference configuration presented in the previous section. They are fully three-dimensional. Figure 2 shows on the vertical middle plane section the velocity vectors with contours of axial velocity (Fig. 2a) and contours of temperature (Fig. 2b). White dotted lines show the extent of the porous medium. The flow and temperature fields on this plane are characteristic of natural convection in a differentially heated fluid cavity on the left of the heat exchanger (with adiabatic conditions on all boundaries except for the heat exchanger section). On the right of the heat exchanger, the flow and temperature field resemble that of natural convection in a differentially heated porous-fluid cavity (with prescribed temperature on all boundaries). However the heat exchanger is not a wall and therefore there is fluid passing from one cavity to the other. Also the heat exchanger is not ideal and the condition T = T_h is not imposed homogeneously on the cross-section.

Figure 3 shows on several vertical cross-sections the contours of vertical velocity (left) and of temperature (right). The corresponding cross sections are shown on the sketch at the bottom of Figure 3. The flow and temperature fields at z = z₁ (Fig. 3ab1), on the left of the heat exchanger, are characteristic of natural convection flow in a differentially heated cylindrical cavity, with adiabatic boundaries. On a section cutting through the heat exchanger, at z = z₂ (Fig. 3ab2) the flow is constrained between the solid rings that constitute the heat-exchanger and hot temperature T_h is maintained on the rings. On a section cutting through the porous medium, at z = z₃ (Fig. 3ab3), the flow and temperature field are essentially characteristic of natural convection in a differentially heated porous cylindrical cavity, with fixed temperature on the tube wall. Finally, on the right part of the porous medium, at z = z₄ (Fig. 3ab4), natural convection flow within a fluid cavity with cold temperature on the boundaries is recovered.The flow is slower and the temperature field becomes more homogeneously cold. There is a clear symmetry of flow and temperature with respect to the vertical axis in all figures. The maximum vertical velocity U_max = 0.115 m/s is read on the scale of Figure 3. The maximum axial velocity W_max = 0.17 m/s is read on the scale of Figure 2a. These are of the same order of magnitude as the reference velocity scales U_{ref p} and U_ref calculated in the previous section. In the porous medium the average velocity amplitude is 0.01 m/s. Using the hydraulic radius $r\simeq \sqrt{(8K/\varepsilon )}=3$ mm, this yields the pore Reynolds number Re_p ≃ 1.7. This value is actually border line for justifying the use of Darcy-Brinkman model without the introduction of the quadradic Forchheimer drag and the advective terms. The transition region between models corresponds to Re_p between 1 and 10 [7].

For this reference configuration, the boundary layer thickness can be estimated for the porous domain using ${\delta }_{\nu }\simeq {\delta }_{\kappa }\simeq HR{a}_K^{-1/2}=9.6\times {10}^{-3}$m [7]. Therefore since the mesh is refined near the porous medium’s left and right sections there should be about 20 mesh points within the boundary layer, which is plenty. In the fluid domain on the right of the porous medium, the temperature difference is much lower (it can be estimated to be about 25 K, based on numerical results), which gives an estimate of a fluid Rayleigh number Ra_H ≃ 2.4 × 10⁴ and the boundary layer thickness ${\delta }_{\nu }\simeq {\delta }_{\kappa }\simeq HR{a}_H^{-1/4}=2\times {10}^{-3}$m [7]. There should be about 5 points within the boundary layer, which is sufficient.

The reference configuration being most restrictive (for all other cases considered in the next section the boundary layer thickness is larger) numerical results can be considered as converged.

Fig. 2 Reference configuration, vertical middle plane section: (a) Velocity vectors with contours of axial velocity (m/s). (b) Contours of temperature (in K), white dotted lines show the extent of the porous medium.

Fig. 3 Reference configuration, contours of (a) vertical velocity (m/s) and (b) temperature (in K), for (1) z1 = 0.05 m, (2) z2 = 0.101 m, (3) z3 = 0.12 m and (4) z4 = 0.2 m. Locations of z1 to z4 are shown on the sketch at the bottom of the figure.

3.3 Parametric study

In this section, several parameters characterizing the porous medium are varied, keeping a fixed porosity: permeability, anisotropy, equivalent thermal conductivity.

In order to analyse the influence of each parameter, we define several flow and temperature output parameters: The maximum vertical velocity amplitude in the tube is noted U_max and ${(x_{\text{max}},y_{\text{max}},z_{\text{max}})}_U$ is the corresponding location. The maximum axial velocity amplitude in the tube is noted W_max and ${(x_{\text{max}},y_{\text{max}},z_{\text{max}})}_W$ is the corresponding location. The total kinetic energy within the tube is noted E_c. Also, the ratio of vertical versus horizontal temperature gradient in the core of the porous medium is noted $\gamma =\left|\frac{\partial T}{\partial x}\right|/\left|\frac{\partial T}{\partial z}\right|$. It is a measure of the strength of natural convection: if conduction is dominant γ → 0 and if convection is dominant γ →∞.

The value of the permeability K depends on the geometry of the porous medium. Usually in thermoacoustic configurations, the porosity is fairly large (between 0.75 and 0.85), and the permeability is kept small to maximize surface exchange (between 10⁻⁶ m² and 10⁻⁹ m²). In the following paragraph, the permeability is varied between these values, considering the porous medium to be homogeneous and isotropic with constant porosity ɛ = 0.84. The hydraulic radius a corresponding to each value of K can be estimated from the relationship K = a²ɛ∕8 (for capillaries of square cross-section). The hydraulic radius varies between ≃ 3 mm (K = 10⁻⁶ m²) to ≃ 0.1 mm (K = 10⁻⁹ m²).

Table 3 presents the results for 4 different values of the porous medium permeability, when the porous medium is isotropic. As K decreases, the total kinetic energy and the γ ratio also decrease, which is to be expected because the porous medium acts more and more like a solid barrier and natural convection flow looses strength. It is interesting to note that the maximum vertical velocity decreases and that its location shifts from just right of the heat exchanger (on the central horizontal plane, facing the middle of the outer heat exchanger fluid space) to just outside the porous medium as permeability decreases (slightly higher than the tube center). The maximum axial velocity also decreases with K and its location (in the right fluid cavity, on the central vertical plane, facing the middle of the outer heat exchanger fluid space) shifts to just outside the porous medium as permeability decreases. This is consistent with natural convection flow developing outside the porous medium when K is very small. The orders of magnitude of the maximum velocities vary with the value of K, but not proportionally as in a simpler configuration with natural convection through a porous medium (see the reference velocity scales U_ref and U_refp presented inthe previous section).

Figure 4 shows on the vertical middle plane section, for K = 10⁻⁷ m², the velocity vectors with contours of axial velocity (Fig. 4a) and contours of temperature (Fig. 4b). White dotted lines show the extent of the porous medium. Although velocities are much smaller than in the reference case, the effect of natural convection on the temperature distribution is clear and a much lesser fraction of the porous medium is submitted to a horizontal temperature gradient. Therefore the thermoacoustic energy conversion should be lowered significantly. In this configuration the hydraulic radius is on the order of 1 mm, the average velocity amplitude estimated in the porous medium is about 5.5 mm/s and the pore Reynolds number is ≃0.3, which validates the use of Darcy-Brinkman equation.

Figure 5 shows on the vertical middle plane section, for K = 10⁻⁸ m² the velocity vectors with contours of axial velocity (Fig. 5a) and contours of temperature (Fig. 5b). White dotted lines show the extent of the porous medium. It is clear on this figure that the porous medium behaves almost like a conducting solid. On the right end of the porous medium, the temperature is only slightly higher than T_c, so that convection is greatly diminished and the maximum velocities U_max and W_max are much smaller than in the reference case (on the order of 0.05 m/s, which corresponds to the fluid velocity scale U_ref for H = R and ΔT = 2 K). On the left of the heat exchanger, the temperature is almost homogeneous (equal to T_hot) and there is nearly no flow. Similar observationscan be made for K = 10⁻⁹ m². The pore Reynolds number are ≃10⁻² for K = 10⁻⁸ m² and ≃3 × 10⁻⁴ for K = 10⁻⁹ m², which validates the use of Darcy-Brinkman equation.

The equivalent porous medium thermal conductivity is varied next, for homogeneous isotropic porous media, keeping all other parameters corresponding to the reference case. Four values were tested, from k_p = 0.028 to k_p = 1.4 W.m⁻¹.K⁻¹, characteristic of porous media used in thermoacoustics [12]. It was shown that the flow and temperature field variations are negligible within this range of values.

Anisotropy is considered next: in the case of pores that are horizontal capillary tubes, the porous medium is almost impermeable in the transverse direction, so that K_xx ≃ 0. The permeability in the axial direction K_zz is varied between 10⁻⁶ m² and 10⁻⁹ m². Throughout this study the equivalent porous medium thermal conductivity is assumed to be homogeneous and isotropic. Figure 6 shows on the vertical middle plane section the velocity vectors with contours of axial velocity (Fig. 6a) and contours of temperature (Fig. 6b), for the anisotropic case, K_zz = 10⁻⁶ m² (all other parameters being maintained the same as in the reference configuration). The flow structure is strongly modified inside the porous medium, since it is constrained horizontally and flows through the fluid space between the heat exchanger rings (compare with Fig. 2a). The temperature field is also modified (compare with Fig. 2b). The solution inside the porous medium is analogous to natural convection flow in a shallow horizontal cavity, with on the top half hot fluid flowing to the right and on the bottom half cold fluid flowing to the left. Then on the right of the porous medium, a more standard natural convection fluid cavity flow is recovered. For this configuration, the maximum vertical velocity is 0.12 m/s measured at (x = R∕2, y = 0, z = 0.151 m) (just on the right of the porous medium) and the maximum axial velocity is 0.17 m/s measured at (x = 9R∕10, y = 0, z = 0.165 m) (also on the right of the porous medium). The total kinetic energy is E_c = 2.84 × 10⁻⁴ J (smaller than in the reference configuration) but the ratio γ is equal to γ = 4.56 (larger than in the reference configuration).

Fig. 4 Isotropic porous medium, K = 10−7 m2, vertical middle plane section: (a) Velocity vectors with contours of axial velocity (m/s) and (b) contours of temperature (in K), white dotted lines show the extent of the porous medium.

Fig. 5 Isotropic porous medium, K = 10−8 m2, vertical middle plane section: (a) Velocity vectors with contours of axial velocity (m/s) and (b) contours of temperature (in K), white dotted lines show the extent of the porous medium.

Fig. 6 Anisotropic porous medium, Kzz = 10−6, on the vertical middle plane section: (a) Velocity vectors with contours of axial velocity (m/s) and (b) contours of temperature(in K), white dotted lines show the extent of the porous medium.

3.4 Consequences for thermoacoustic applications

The previous simulations of natural convection flows were performed on configurations (geometry, materials, applied temperature differences) applicable to thermoacoustic engines, even though there was no acoustics.

Based on these simulations and observations, steady natural convection flows should not be neglected in the design of standing-wave thermoacoustic engines. The onset of an acoustic wave is based on the existence of a temperature gradient regularly distributed along the porous medium. Since natural convection flows induce very distorted temperature fields, they can delay or even prevent the onset of acoustic waves. Natural convection flows could interfere with acoustics and also with other nonlinear phenomena occuring in such thermoacoustic devices, such as streaming flow. It would be interesting to run some simulations of acoustics and natural convection together.

Finally such flows may have applications for other devices involving thermoacoustics, such as thermoacoustic heat pumps or traveling-wave thermoacoustic engines, as long as there is a horizontal temperature difference between the constituting elements.

4 Conclusion

A three-dimensional numerical study of natural convection flow was performed using a finite volume code for solving mixed Navier-Stokes and Darcy-Brinkman equations under Boussinesq approximation, in a configuration applicable to thermoacoustic devices. Theinfluence of the porous medium’s physical characteristics (permeability, thermal conductivity, anisotropy) on the flow and temperature fields was investigated. It was shown that such flows are fully three-dimensional and can modify significantly starting as well as steady operating conditions of thermoacoustic engines. Simulations of both acoustics and natural convection in thermoacoustic configurations are in progress, in order to analyse their interaction.

Nomenclature

Acknowledgements

This work was partly financed by ANR project TACOT (Projet-ANR-17-CE06-0007).

References

All Tables

All Figures

	Fig. 1 Computational domain (top) and thermal boundary conditions (bottom).
In the text

	Fig. 2 Reference configuration, vertical middle plane section: (a) Velocity vectors with contours of axial velocity (m/s). (b) Contours of temperature (in K), white dotted lines show the extent of the porous medium.
In the text

	Fig. 3 Reference configuration, contours of (a) vertical velocity (m/s) and (b) temperature (in K), for (1) z1 = 0.05 m, (2) z2 = 0.101 m, (3) z3 = 0.12 m and (4) z4 = 0.2 m. Locations of z1 to z4 are shown on the sketch at the bottom of the figure.
In the text

	Fig. 4 Isotropic porous medium, K = 10−7 m2, vertical middle plane section: (a) Velocity vectors with contours of axial velocity (m/s) and (b) contours of temperature (in K), white dotted lines show the extent of the porous medium.
In the text

	Fig. 5 Isotropic porous medium, K = 10−8 m2, vertical middle plane section: (a) Velocity vectors with contours of axial velocity (m/s) and (b) contours of temperature (in K), white dotted lines show the extent of the porous medium.
In the text

	Fig. 6 Anisotropic porous medium, Kzz = 10−6, on the vertical middle plane section: (a) Velocity vectors with contours of axial velocity (m/s) and (b) contours of temperature(in K), white dotted lines show the extent of the porous medium.
In the text