The Effect of Internal and External Heating on the Free Convective Flow of a Bingham Fluid in a Vertical Porous Channel

: We study the steady free convective ﬂow of a Bingham ﬂuid in a porous channel where heat is supplied by both differential heating of the sidewalls and by means of a uniform internal heat generation. The detailed temperature proﬁle is governing by an external and an internal Darcy-Rayleigh number. The presence of the Bingham ﬂuid is characterised by means of a body force threshold as given by the Rees-Bingham number. The resulting ﬂow ﬁeld may then exhibit between two and four yield surfaces depending on the balance of magnitudes of the three nondimensional parameters. Some indication is given of how the locations of the yield surfaces evolve with the relative strength of the Darcy-Rayleigh numbers and the Rees-Bingham number. Finally, parameter space is delimited into those regions within which the different types of ﬂow and stagnation patterns arise.


Introduction
This short paper considers the steady flow which is induced when a vertical porous channel is heated both externally and internally, and when that porous medium is saturated by a Bingham fluid. These fluids form, perhaps, the simplest model of a yield stress fluid where the fluid shears when the applied stresses are greater than a threshold value, the yield stress, but acts like a solid when the applied stresses are too small. Within the context of a porous medium the natural yield threshold is expressed in terms of a pressure gradient and, in the context of convective flows, this includes the buoyancy force. Pascal's piecewise-linear law (Pascal [1,2]) is itself the simplest relationship, and in this model there is no flow (as opposed to zero shear) when the threshold has not been met.
One of the earliest models for the flow of a yield stress fluid in a porous medium is by Gheorghitza [3] who, like Pascal, cites a number of authors reporting the presence of a yield threshold in filtration flows although Gheorghitza does not name the fluid as a Bingham fluid. Gheorghitza names this phenomenon as an initial gradient. Wu and Pruess [4] adds the practical observation that heavy oils in reservoirs and water within clay soils also exhibit a threshold gradient. Bingham fluids also arise elsewhere and these include drilling mud, ceramic pastes, yoghurt, mayonnaise, sewage sludges and magma.
Other models for the flow of a Bingham fluid in a porous medium also exist. One of the most frequently used is the Buckingham-Reiner model (Buckingham [5] and Reiner [6]) which, strictly speaking, corresponds to the Hagen-Poiseuille flow of a Bingham fluid, but it may be applied to a porous medium by assuming a unidirectional flow within a medium consisting of identical pores. More complicated scenarios may be envisaged, and it is worth mentioning that Nash and Rees [7] performed an analytical study of the effect of having distributions of pore diameters. In general this was found to 'soften' the initial stages of flow post-threshold, and each distribution of pores has its own analogue of a Buckingham-Reiner law. The work of Malvault et al. [8] relaxes the assumptions of [7] that the pores have a uniform cross-section.
The present paper considers what might be regarded as a very straightforward steady free convection problem in a vertical porous channel. The channel is heated via a uniform heat generation mechanism, but the walls are also held at different temperatures from one another. The resulting temperature distribution is simple to determine, and should the porous medium be saturated by a Newtonian fluid, most of the resulting analysis centres around the analytical determination of a suitable reference temperature. When the porous medium is saturated by a Bingham fluid then the locations of yield surfaces need to be found, and allowance has to be made for when a yield surface attaches to a bounding surface. In more detail, we find that four different flow regimes arise where the first corresponds to complete stagnation, the second to having one stagnant region, and two more to having two stagnant regions. Much depends on the competing effects of the internal and external heating mechanisms and how these interact with the yield threshold. This work is part of a study of different aspects of porous channel and boundary layer flows involving Bingham fluids; see Rees and Bassom [9][10][11]. These earlier works consider the unsteady evolution of yield surfaces in such flows.

Governing Equations
In the present paper we shall adopt Pascal's model (Pascal [1,2] to describe how a Bingham fluid flows within a homogeneous and isotropic porous medium. This is a piecewise linear relationship between the fluid velocity and the applied pressure gradient and it may be expressed as follows, where the quantity, G, denotes the threshold pressure gradient above which the fluid yields and flows but below which fluid is stagnant. Here K is the permeability, µ is the plastic viscosity, p is the pressure, z the vertical coordinate and w the vertical velocity. Here we assume that the plastic viscosity is a constant which corresponds to a piecewise-linear stress/strain relationship, and therefore Equation (1) models a Bingham fluid, rather than a Herschel-Bulkley fluid. More complicated expressions than Equation (1) exist, such as those derived in Nash and Rees [7] for different tube bundle distributions. The simplest expression is the Buckingham-Reiner model [5,6] which softens the discontinuous gradient in Pascal's law. Other expressions derived in [7] provide for slower approaches to the eventual linear relationship between the flow and the pressure gradient and depend on the probability distribution associated with the pore diameters. In the present paper we restrict ourselves solely to Pascal's model; the adoption of anything more detailed will only result in small quantitative differences in the solutions presented below. We are concerned with free convective flows for which buoyancy also acts as a body force. Therefore Equation (1) becomes, In this equation ρ f is the reference density of the fluid, g gravity, β the coefficient of cubical expansion and T the temperature. In the above we have assumed that the Boussinesq approximation applies when writing down the buoyancy term, and have taken T c to be the reference temperature which will also be the coldest temperature experienced by the medium. The full two dimensional heat transport equation is where C is the heat capacity, α = k pm /(ρC) pm is the thermal diffusivity and q is the uniform rate of heat generation. Both the temperature and the resulting velocity fields are independent of z and are functions solely of x, the horizontal coordinate. Thus the heat transport equation reduces to the ordinary differential equation, Finally, the boundary conditions for the temperature are that and where the channel itself lies in the range 0 ≤ x ≤ d.
The governing equations, namely Equations (2) and (4), and the boundary conditions, (5), may be nondimensionalised using the following scalings, where T h is the temperature of the hot surface. We obtain, and subject to, θ = 0 at x = 0, and θ = 1 at x = 1.
In the above the external and internal Darcy-Rayleigh numbers are given by respectively, and the Rees-Bingham number, which may be interpreted as being a convective porous Bingham number or as a nondimensional threshold gradient for the fluid. We note that it is possible to reduce the set of three nondimensional parameters to two and therefore it is also possible to give a comprehensive account of the present problem using any two independent ratios of the three. However, this degeneracy is removed once convection becomes nonlinear. The computations undertaken by Rees [12] consider convection in a rectangular porous cavity with the same external heating as here but without internal heating (Ra i = 0). The upper and lower surfaces of the cavity are insulated. The flow in [12] is two-dimensional as compared with the present one-dimensional flow. It is clear from the results of [12] that solution curves (such as the variation of the Nusselt number) do not map onto a single curve when plotted against Ra/Rb, which they will do when the channel is infinitely tall. Therefore it was decided to retain the present three parameters because of the context provided by [12].

Numerical Solutions
The solution for θ may be written down easily: and it is this which is used in Equation (7). We note that, when there is no internal heating (Ra i = 0), then the temperature profile is odd about x = 1 /2, and when there is no sidewall heating (Ra = 0) then the temperature profile is even about x = 1 /2. Later we shall see that the present system admits situations where either the whole channel is stagnant, or else that there exists either one or two stagnant regions bounded by flowing regions. An alternative viewpoint is that, when flow exists, then there may be two, three or four yield surfaces. One example of a situation in which there are four yield surfaces and two stagnant regions is shown in Figure 1 and this corresponds to the case, Ra = 3, Ra i = 50 and Rb = 1. The temperature profile has no symmetry and the locations of the four yield surfaces are denoted by the values, x 1 , x 2 , x 3 and x 4 . These values satisfy the following equations, Ra and are obtained by setting w = 0 into the first or third options in Equation (7). While the present flow is a free convection flow in an infinitely long channel, it also approximates very well the flow which will arise in a tall channel sufficiently far from the upper and lower boundaries, and therefore it is necessary to apply a zero upward flux condition. Thus we have, which is equivalent to, After integration this translates into, In the above we have used x 5 = 1 solely to show how the constant which arises when evaluating Equation (17) is related to the rest of the integrand.
Equations (13)-(16) and (19) form a set of five nonlinear equations for the four yield surfaces and the value of the vertical pressure gradient correction, p z . This correction needs to be found because T c was chosen to be the reference temperature and therefore the corresponding hydrostatic pressure gradient is incorrect if an overall zero mean flux is to be maintained. To illustrate this, if we choose to consider a Newtonian fluid (Rb = 0) in a channel with external heating only (Ra i = 0), then p z = 1 2 Ra corresponds to a zero mean flow, and w = Ra(x − 1 2 ). A similar analysis for solely internal heating yields p z = 1 12 Ra i . Both of these values may be found in the final line of Equation (19) when Rb set to be equal to zero, noting again that x 5 = 1. For general cases these five equations were solved using a standard Newton-Raphson method and therefore our solutions are essentially exact.
Showing the temperature and velocity profiles for the case Ra = 3, Ra i = 50 and Rb = 1. The locations of the yield surfaces are x 1 , x 2 , x 3 and x 4 . The two stagnant regions lie in the ranges, We note that when Ra is increased gradually for the parameter set shown in Figure 1, then the rightmost yield surface at x = x 4 moves towards the right hand boundary at x = 1 and eventually reaches it. In our code this was modelled by simply replacing Equation (16) by x 4 = 1. This corresponds to the final integral in Equation (18) being zero, and is modelled correctly in Equation (19) by having the last two bracketed terms cancelling one another. A further increase in Ra eventually leads to the x = x 3 yield surface reaching x = 1, which means that the flow now has only one stagnant region. In our code this was modelled by replacing Equation (15) by x 3 = 1. Therefore, for practical reasons, we always began a computation with four yield surfaces and followed their trajectory as one of the governing parameters was varied.

Rb
x (c) Ra = 0, Ra i = 120 When the heating is purely external then the system of equations for determining p z and the yield surfaces may be solved analytically. We find again that p z = 1 2 Ra, and that the fluid velocity is given by, It is clear from this expression that that there is flow only when 0 ≤ Rb ≤ 1 2 Ra and in this range of values of Rb the velocity profile is piecewise linear. The locations of the yield surfaces may be gleaned from Equation (20) and these linear functions of Rb are shown in Figure 2a. We have full stagnation when Rb > 1 2 Ra. When the heating is purely internal then it is not possible to find an analytical expression for p z or for the locations of the yield surfaces even though the symmetry of the system means that x 1 + x 4 = x 2 + x 3 = 1. Therefore we have resorted to purely numerical means to determine where the yield surfaces are as a function of Rb, and these are shown in Figure 2c. However, it is possible to find where the yield surfaces arise when Rb begins to increase from zero. Given that Ra = 0, Equations (13) and (14) with p z = 1 12 Ra i and Rb = 0 give both these values may be seen in Figure 2c. It is also possible to determine the value of Rb above which the channel becomes fully stagnant. This is achieved by setting x 1 = 0, x 2 = x 3 = 1 2 and x 4 = 1 into Equations (13)-(16), and then we find that stagnation corresponds to Rb > 1 16 Ra i . The intermediate range of cases is represented by the solutions shown in Figure 2b for which Ra = 10 and Ra i = 120. Once Rb has risen above zero, two narrow stagnant regions appear but these are not symmetrically placed about x = 1 2 . The right hand yield surface attaches onto the x = 1 boundary when Rb 3.831. Flow weakens as Rb increases further and full stagnation arises when Rb = 245 12 = 10.208333. Once more this value may be found analytically by first substituting x 1 = 0 into Equation (13), which yields p z = Rb, and then by noting that Equations (14) and (15) must represent a double zero since x 2 = x 3 . These two equations may be rearranged into the form, It is clear that there will be two different solutions for this equations (i.e., x 2 and x 3 ) when the right hand side is positive, but none when it is negative. Therefore incipient stagnation corresponds to when the solutions are equal, for which the right hand side must be zero. Hence, represents the critical value of Rb in general. For the example shown in Figure 2b we have Rb = 245 12 , as quoted above. Under these conditions x 1 is also zero and therefore full stagnation arises for larger values of Rb. The common values of x 2 and x 3 are now given by, and hence x 2 = x 3 = 7 12 for the case shown in Figure 2b. It is of interest to try to determine a general condition for stagnation to occur. The expression given in Equation (22) may be rearranged slightly to yield, However, this formula was derived by assuming that x 2 = x 3 , i.e., that the coaslescence of two yield surfaces takes place in the interior of the domain, and therefore it is essential to check if the corresponding value of x 2 = x 3 does indeed lie in the interior. Equation (21), with a zero right hand side, tells us that x 2 = x 3 = 1 2 + Ra Ra i and therefore the above analysis clearly applies only when Ra Ra i ≤ 1 2 . It is straightforward to check that this criterion is also the criterion that the maximum value of θ given in Equation (12) is an internal maximum. Therefore we have a simple delineation between two separate regimes, one with an interior maximum for θ and one where that maximum lies on the right hand boundary.
An illustration of the approach to stagnation as Rb increases for cases where Ra Ra i > 1 2 is shown in Figure 3 where Ra = 10 has been chosen and where the three separate values, 0, 10 and 20, have been taken for Ra i . In all cases stagnation occurs when Rb = 5, and more generally this will be when In Figure 3 we see that the curve for Ra i = 20, which is a transitional case because Ra Ra i = 1 2 , approaches x = 1 with a zero slope.

Variation with Ra i
The detailed Figure   When the external Darcy-Rayleigh number begins to rise from zero the pattern of flow and stagnation loses its symmetry. Flow begins at the right hand boundary at a larger value of Ra i /Rb than does flow at the left hand boundary and in the middle region both of which begin to flow simultaneously.
While Ra/Rb < 2, buoyancy forces remain too weak to cause flow when the heating is purely external, and therefore stagnation continues to be found at Ra i /Rb = 0. But when Ra/Rb passes through 2, the shapes of the yield surfaces change dramatically, and evolve from one continuous region to two.
When Ra/Rb = 2.1 very narrow regions of flow occur at the two boundaries when Ra i /Rb = 0 because buoyancy is only just in excess of what is required to overcome the yield threshold. But at larger values of Ra/Rb the two disconnected unyielded domains narrow, and, for Ra/Rb = 3 we see quite clearly the transition as Ra i /Rb increases from a flow pattern which is antisymmetric when Ra i /Rb = 0 to one which is symmetric when Ra i /Rb becomes large. Even when Ra/Rb = 50, the transitions shown for smaller values of Ra/Rb also occur but do so at much larger values of Ra i /Rb. Figure 5 summarises all of our discussion about when stagnation occurs, but we have reinterpreted the data in terms of the variation of Rb/Ra with Ra i /Ra. This figure delineates different regions of parameter space (ii) has a single stagnant region with flow either side, (iii) has two stagnant regions but the right hand one is attached to the right hand boundary, and (iv) has two stagnant and three flowing regions. These are indicated schematically on the figure itself for ease of interpretation.  Figure 5 It is possible to explain analytically many of the features exhibited in Figure 5 and for these purposes it is convenient to define the ratios

Detailed Analysis of
The dotted line separates the regions where the temperature field has an internal maximum ( Ra i > 2) and where it has a maximum at x = 1 ( Ra i < 2). The boundary of the fully stagnant region has already been shown to be given by Equation (24) respectively. The other two lines in Figure 5 were obtained by suitably modified Newton-Raphson solvers and some of the numerical data corresponding to the yield surfaces just on the point of attaching to x = 1 in Figure 4. The middle region and the lower right hand region are distinguished by the location of the right hand stagnant region; in the former case it is attached to x = 1 but is not in the latter case.
We can derive expressions for the two boundaries in Figure 5 that meet at Ra i = 6. The left one of these corresponds to the case when the single stagnant region that is present when Ra i and Rb are both relatively small evolves into a pattern for which another stagnant zone begins to form at x = 1. The boundary is defined by x 3 = 1 and then combining Equations (14) and (15) gives that x 2 = 2/ Ra i . We can then determine x 1 in terms of x 2 and Rb by eliminating the pressure term between Equations (13) and (14). Finally, the flux condition Equation (19) then simplifies so that if then X satisfies the cubic This equation does not possess a simple analytic solution (although such an equation has solutions that can be written down albeit in very complicated form), but we can confirm some elementary results. First, when Ra i = 6 we have X = 2 and then Rb = 0 as expected from Figure 5. Moreover, if Ra i = 6 − δ, with 0 < δ 1, then X = 2 + δ + · · · and Equation (28) gives Rb = 1 8 δ + · · · . This suggests that the boundary has slope, − 1 8 , at Rb = 0. Equation (29) also enables us to deduce the local behaviour near the other end of the boundary at ( Ra i , Rb) = (2, 1 2 ) where it joins with the lines Equation (27). Near this point, if Ra i = 2 + δ then Equation (29) leads to X = 2 + 1 2 δ − 1 2 √ 3 δ 3/2 + · · · and then to Rb = 1 2 − 1 4 √ 3 δ 3/2 + · · · . The presence of fractional powers is required to resolve an apparent contradiction at O(δ 3 ) in the expansion of Equation (29), but this matches perfectly the numerical data of the appropriate curve in Figure 5.
We now turn to the last boundary on Figure 5 that separates those flows with two stagnant regions with the right hand one attached to x = 1 and the situation when there are two stagnant and three flowing regions. We can pursue an analysis that in many ways parallels the argument just above. The boundary of interest arises when x 4 = 1; then Equations (13) and (16) lead quickly to the conclusion that x 1 = 2/ Ra i . Equation (16) yields the pressure gradient p z = 1 + Rb while Equation (14) shows that x 2 + x 3 = 1 + (2/ Ra i ). If these relationships are substituted into the flux condition and if we define (cf. Equation (28)) then we find that Z 3 = 3 2 Ra i − 1. Hence is an analytical description of the boundary. It follows from this that, for small values of ( Ra i − 6), we have Rb ≈ 1 16 ( Ra i − 6) + · · · . We can also understand the behaviour of the system for small values of Rb and illustrated in Figure 6. If we consider a fixed small value of Rb, which is taken to be 0.01 or 0.05 in Figure 6, then for values Ra i much less than 6 the flow contains a single stagnant zone. As Ra i approaches 6 so a second stagnant region forms on x = 1 and later it detaches so that the two stagnant zones separate three flowing regions. We can capture all this behaviour analytically by seeking solutions when Ra i = 6 + c Rb for values of c = O(1). We suppose that the points x 1 -x 4 are located at If we substitute these expressions into Equations (13)-(16) and Equation (19) then at O( Rb) we obtain four linear equations for x 1x 4 whose solution gives that We remark that there may appear to be a something of a contradiction here because for sufficiently negative values of c these predictions for x 3 and x 4 do not lie in the region 0 ≤ x ≤ 1. It turns out that the results quoted for x 1 and x 2 hold to the accuracy quoted irrespective of whether x 3 and/or x 4 lie inside the channel or not. We can then interpret the results encompassed by Equation (33) in conjunction with Figures 5 and 6 in the following way. For sufficiently large negative values of c only x 1 and x 2 lie within the channel so that there is a single isolated stagnant region corresponding to the lower left-hand part of Figure 5. As we slowly increase c then the first restructure of the flow occurs when the predicted value of x 3 becomes 1; so that c = −8. (We point out that at this stage Rb = 1 8 (6 − Ra i ) confirming that the slope of the boundary here is − 1 8 .) For c slightly larger than this value, then the flow consists of one stagnant zone around x = 1 3 while a second stagnant region next to the wall at x = 1. The signal for this second region detaching from the x = 1 is that x 4 = 1 which clearly occurs when c = 16. Now Rb = 1 16 ( Ra i − 6) confirming that the slope of this bounding curve in Ra i / Rb parameter space is indeed 1 16 . For values of c greater than 16 the flow consists of two stagnant zones each away from either boundary.
This analysis is all consistent with the results shown in Figure 6 for two relatively small values of Rb. When Ra i 6 there is the single stagnant zone isolated from the walls of the channel. As Ra i grows so the centre of the stagnant region drifts towards x = 1/3 until at a value of Ra i slightly less than 6 there is evidence of the second stagnant zone forming on x = 1. As Ra i increases further, this second stagnant region soon detaches itself and we are left with three flowing regions separated by two stagnant zones (corresponding to region (iv) in Figure 5). We may use also Figure 5 in a qualitative manner by, for example, choosing a case for which Rb/Ra = 1 4 while Ra i /Ra = 0, i.e., pure external heating. This case lies in that part of the Figure for which there is a single stagnant region in the interior. Once more we see that as the strength of the internal heating is increased so a new stagnant region is induced at the right hand boundary, and then this ultimately detaches. We also see that it is impossible to have a single interior stagnant region when Ra i /Ra > 6 irrespective of the value of Rb.

Conclusions
In this short paper we have attempted to provide a comprehensive description of the properties of the flow in a vertical porous channel which is subjected to both internal heating (i.e., uniform heat generation) and external heating (i.e., a temperature difference across the layer) and where a Bingham fluid saturates the porous medium. We have adopted what is, perhaps, the simplest model for such flows namely the Pascal model. Some analytical results have been provided but the core results have been obtained using a multi-dimensional Newton-Raphson solver.
We have found the criterion for flow to arise and this is given either by Equations (24) or (25) depending on the value of Ra i /Ra, and it is also shown as the uppermost curve in Figure 5. The number and locations of the resulting stagnant regions depend quite strongly on the balance between the values of Rb/Ra and Ra i /Ra. We find that there is only one stagnant region when external heating dominates, but that a gradual transformation to a state where there are two stagnant regions arises as the strength of the internal heating increases.
The qualitative nature of our results will not change should a more accurate form of Pascal's law for the flow of a Bingham fluid in a porous medium be replaced by other models, such as the Buckingham-Reiner law which models porous media which are composed of identical unidirectional pores, or many of the models given in [7] which correspond to more general distributions of pores or channels.

Conflicts of Interest:
The authors declare no conflict of interest.