Unsteady Free Convection Boundary Layer Flows of a Bingham Fluid in Cylindrical Porous Cavities

We consider two unsteady free convection flows of a Bingham fluid when it saturates a porous medium contained within a vertical circular cylinder. The cylinder is initially at a uniform temperature, and such flows are then induced by suddenly applying either a new constant temperature or a nonzero heat flux to the exterior surface. As time progresses, heat conducts inwards and this may or may not overcome the yield threshold for flow. For the constant temperature case, flow begins immediately should the parameter, Rb, which is a nondimensional yield parameter, be sufficiently large. The ultimate fate, though, is full immobility as the cylinder eventually tends towards a new constant temperature. For the constant heat flux case, the fluid remains immobile but will begin to flow eventually should Rb be sufficiently large. The two cases have different critical values for Rb.


Introduction
The topic of the flow of a Bingham fluid when it is not saturating a porous matrix is a wellestablished field of study, and the literature is now quite mature. Efforts continue to develop different numerical methods which will allow an accurate computation of such flows, the chief difficulty being the presence of a yield surface which needs to be computed as part of the overall problem. In some simplified cases, the analysis can proceed analytically (such as for plane-Poiseuille and Hagen-Poiseuille flows), but in other one dimensional problems the analysis has to be completed by using a simple Newton-Raphson iteration equation on a transcendental equation which then allows the positions of the yield surfaces to be found. Examples of such works include those by Yang and Yeh (1965) and Bayazitoglu et al. (2007) who studied steady free convection in a vertical channel which is heated from the side. Steady convection will only ensue once the Rayleigh number is sufficiently large that the yield stress is overcome by the buoyancy force. For this free convection problem, there are two moving zero-shear regions (i.e. ones which act like solids) and three fluid regions undergoing shear. By contrast, Patel and Ingham (1994) considered mixed convection with the combination of buoyancy and a driving pressure gradient, while Barletta and Magyari (2008) studied a free convection version of vertical Couette flow. In both cases the number, the size and the locations of the plug-flow regions are dependent on further nondimensional parameters. One paper whose topic is closer to that of the present is that of Kleppe and Marner (1972) who considered a sudden change in the temperature of one the sidewalls and who then determined the evolution with time of the resulting velocity profile.
When a Bingham fluid saturates a porous medium there are some changes to the modelling of the equations of motion which are caused by the presence of the porous matrix. Just as the Navier-Stokes equations, which apply for a clear Newtonian fluid, are replaced by Darcy's law, which applies in a porous medium, so it is that the usual yield-stress form of the Navier-Stokes equations, which gives zero shear when stresses are less than the yield stress, needs to be replaced by a Darcy-Bingham law. The first papers which presented such a law were those by Pascal (1979Pascal ( , 1981. This law appears later, see Eq. (1), but it provides a piecewise-linear dependence of the flux velocity on the applied pressure gradient; thus, there is a yield pressure gradient. When the pressure gradient (or, equivalently, body forces) is less than a threshold value then there is no flow. Further realism is obtained by considering the porous medium as an assembly of identical channels or pores, which gives the well-known Buckingham-Reiner law (1921). In this case the initial rise in the flow is quadratic immediately post-threshold, as opposed to linear in Pascal's model. We also note that more sedate transitions to flow were found by Nash and Rees (2017) who considered distributions of channels/pores. In the present paper, we will adopt Pascal's piecewise-linear form.
In this short paper, we will consider the flow of a Bingham fluid in a vertical circular cylinder which is filled with a porous medium. The flow is induced by means of applying a sudden heating to the outer surface, and therefore the paper may be categorised in the same way as the works of Kleppe and Marner (1972) and of Bassom (2015, 2016). The present paper extends the analyses of Bassom (2015, 2016) into a cylindrical domain, and therefore the work presented has been made as concise as possible. Therefore, the heating considered here takes two forms: (1) a new temperature and (2) a nonzero heat flux. Using Pascal's piecewise-linear Darcy-Bingham law, the analysis proceeds analytically and again, a Newton-Raphson scheme has to be used in order to locate both the yield surfaces and, given that the flow domain is finite horizontally, the corresponding change to the initial hydrostatic pressure gradient.

Governing Equations
The present paper is concerned with a vertical cylindrical configuration and the subsequent evolution of the velocity and temperature fields once its bounding surface has its temperature raised from the reference temperature. The configuration we use is sketched in Fig. 1. In Case 1 the initial temperature is T 0 and at the time t = 0 the temperature of the outer boundary is raised suddenly to T 1 , while, by way of contrast, Case 2 considers a sudden change in the boundary heat flux.
For the isothermal and unidirectional flow of a Bingham fluid saturating a homogeneous and isotropic porous medium, we may follow Pascal (1981) and write the equation, where G is used to denote the threshold pressure gradient above which the fluid is able to flow. Other terms in Eq. (1) are the common ones used for flows in porous media and they are given in the Nomenclature. In the present paper, we shall considering a free convective flow, and therefore buoyancy forces may also to be included in such an equation and we obtain,  (2) We have taken w to be the vertical velocity here and have also assumed that the Boussinesq approximation applies. The reference temperature is T 0 .
The cases which we consider are regarded as being of infinite length in the vertical direction, although practically a long finite cylinder or annulus is highly likely to be essentially equivalent. This means that the thermal boundary layer which arises initially and its subsequent evolution will be unidirectional with the temperature and the vertical velocity being functions solely of the radius and time. In such cases, the equation of continuity is satisfied.
The governing equation for the heat transport is, where r is the radial coordinate. Here ς is the heat capacity ratio between the porous medium and the saturating fluid, and α is the thermal diffusivity of the porous medium. The governing equations, namely Eqs.
(2) and (3), may be nondimensionalised using the following scalings, where T 1 is the temperature of the hot surface and R is a chosen radius, and we obtain, and In the above equations, the Darcy-Rayleigh number is given by and the Rees-Bingham number is, This latter parameter is clearly a scaled version of the yield pressure gradient, G, and, given the presence of K and α, it could also be described as a porous convective Bingham number. The initial condition in nondimensional form is simply that θ = 0 at t = 0. The two cases which we consider now correspond to the following: Case 1: 0 ≤ r ≤ 1 θ = 1 at r = 1, Case 2: 0 ≤ r ≤ 1 ∂θ/∂r = −1 at r = 1.
Case 1 is considered in the following section, but when Case 2 is considered in Sect. 4, then the modified boundary condition at r = 1 requires minor changes in the scaling for T and for the definition of Ra; these will be introduced in the appropriate place in the text.

Case 1: Constant Temperature
We consider the development of the temperature field within a vertical porous cylinder the initial temperature of which is θ = 0 and where the temperature of the boundary at r = 1 is raised suddenly to θ = 1. The evolution of the temperature field is governed solely by conduction and is unaffected by the induced flow. We may solve Eq. (6) using separation of variables; in this way it is a standard textbook exercise to show that the final solution is and where the λ n values are the positive zeros of the zeroth order Bessel function of the first kind, J 0 . Using standard results involving Bessel functions, this solution may be written in a simpler form: Figures 2 and 3 are concerned with the temperature field. Although such curves are well known and appear in many textbooks and lecture notes, they nevertheless set the context for the resulting flow and are included for that reason. Figure 2 shows the temperature profile at a selected set of times. Clearly a very distinct boundary layer forms at very early times, and, when t 1, it may be shown to be given by the equivalent Cartesian solution, , until curvature effects become significant. The temperature at the origin will also be important in what follows, and therefore Fig. 3 shows this as a function of time. The origin remains at the original temperature until t 0.03 after which the temperature rises towards 1. By the time t = 1 the temperature at the origin is θ = 0.9951. At this point in time, the temperature of the cylinder is almost uniform and therefore there are no buoyancy forces available to drive a flow. That this should be so is evident if the present model problem forms the central portion of a very tall closed cylindrical cavity where the mean flow up the layer be zero. If we consider Eq. (5) in its Newtonian form, i.e. with w = Ra θ − p z , then the initial state, θ = w = 0, means that p z = 0. However, once the cylinder has heated up completely to θ = 1, then the velocity field must again be zero, and therefore p z has now risen to the value given by Ra. This simply corresponds to an adjustment in the hydrostatic pressure gradient because the new uniform temperature is different from the reference value. Therefore, the variation of p z with time will also need to be computed as part of the solution procedure for w; this idea was also used in Bassom (2015, 2016).
The velocity field is given by solving Eq. (5). The chief difficulty is the determination of the location of the two yield surfaces between which the fluid is stationary. If the inner and outer yield surfaces are defined to be r = r 1 and r = r 2 , respectively, then Eq. (5) tells us that and at all times. The third condition which needs to be satisfied is that the overall vertical flow is zero. This means that must be satisfied. After some lengthy manipulations, this condition may be written in the form,  (12), (13) and (15) form three nonlinear equations for the three unknowns, r 1 , r 2 and p z as functions of t. As with our previous papers Bassom (2015, 2016), the solutions were obtained using a multi-dimensional Newton-Raphson scheme which utilises numerical differentiation to obtain the iteration matrix. Figure 4 shows some example velocity profiles for the case, Rb/Ra = 0.1. When t = 0.001 the temperature field, which is also shown in Fig. 4, has penetrated only a small distance from r = 1. The resulting velocity profile consists of a narrow region of upward flow near the other surface, another narrow region of stagnant medium, and a very wide region of low amplitude negative velocity which includes the origin. The closer one gets to t = 0, the weaker the downflow is and the narrower are the regions of upflow and stagnation, as shall be seen in Fig. 5.
Once time has reached t = 0.01 both the upflow and stagnant regions have widened substantially, and the magnitude of downflow in the centre of the cylinder has also increased. As t increases further the magnitude of the downflow and the width of the upflow region both reach a maximum and then decrease once more. The primary reason for this is the heating up of the centre of the cylinder which brings about a decrease in the available buoyancy force. Just before the fluid becomes completely stagnant, the upflow and downflow have both In Fig. 5, we give a summary of how the yield surfaces vary in time for different values of Rb/Ra. For each value of Rb/Ra two curves are presented and the fluid is stagnant in the region between those curves. Our previous comment about the narrowness of the stagnant region and its proximity to r = 1 at early times is clearly evident in that Figure. In addition, we also see that the width of the stagnant region increases as Rb/Ra increases; that this should be so is because of the reduced effective buoyancy when Rb is large. We also see how, in all cases, the stagnant region eventually expands to fill the cylinder, but stagnation happens earlier for larger values of Rb. Figure 6 shows what might be called the stagnation time as a function of Rb/Ra. This is computed easily by first substituting r 1 = 0, r 2 = 1 into Eqs. (12) and (13) to obtain, and by solving for t in terms of Rb/Ra using a single-unknown Newton-Raphson scheme. Figure 6, however, is created easily by finding Rb/Ra as a function t. Given the form of Eqs. (12) and (13), it is clear that it is impossible to have any flow at any time should Rb/Ra > 1 /2, and therefore flow cannot be initiated. On the other hand, the behaviour of the curve shown in Fig. 6 as Rb → 0 may be found using Eq. (11). If we assume that Rb/Ra is small and that t is large, then it is sufficient to retain one term in the summation in (11) For example, if we were to have a very weak yield threshold, one for which Rb/Ra = 10 −6 , then Eq. (17) gives a stagnation time of t = 2.351. This time increases to 4.739 when Rb/Ra = 10 −12 .
Finally, Fig. 7 shows that the hydrostatic pressure gradient always increases in time from a value which is identical to Rb (i.e. p z /Ra increases from Rb/Ra). It is only in the case of a Newtonian fluid that p z /Ra rises to a unit value.

Case 2: Constant Heat Flux
In this section we shall change the boundary condition on the outside of the cylinder from one with a constant temperature to one with a constant heat flux. Therefore, we consider the development of the temperature field where the initial temperature is θ = 0, where the boundary condition at r = 1 is changed suddenly to ∂θ/∂r = −1. In dimensional terms, this means that we shall apply and that the Darcy-Rayleigh number, previously defined in Eq. (7) for Case 1, is now defined as, Once more, we may solve Eq. (6) using separation of variables, and standard techniques may be used to show that the final solution is where the σ n values are the positive zeros of J 0 , or, equivalently, the positive zeros of J 1 . Figure 8 shows how the temperature field evolves with time. Once more a thin boundary layer develops at early times, and it may easily be shown that θ ∼ 2t 1/2 ierfc[(1 − r )/2 √ t]; see Rees and Bassom (2016). Therefore, we expect the outer surface temperature to rise initially as 2(t/π) 1/2 before curvature effects become significant. But without the knowledge of the solution given in Eq. (20), it is clear that when there is uniform heat flux into the cylinder, then the overall heat content must rise linearly with time. Like Fig. 2, this Figure shows the presence of a narrow thermal boundary layer at small times which thickens until the origin is reached, but unlike in Fig. 2, the maximum temperature continues to rise. The evolution with time of both the boundary and central temperatures is shown in Fig. 9, together with the large-t asymptotic states, namely The exact solution and the large-t asymptotic solutions cannot be resolved visually once t is as large as 0.4. At this time, the asymptotic solution for the mid-point temperature is in error by roughly 0.001, while that for the circumference is in error by roughly 0.0004. At t = 0.5, the respective errors are smaller by a factor of 5. Figure 10 shows a selection of velocity profiles for Rb/Ra = 0.05 and 0.2. As for the constant temperature case, the yield locations and the hydrostatic pressure adjustment were computed first by solving Eqs. (12)-(14) using Eq. (20) as the temperature profile. Equation (14), the zero mean velocity condition, is found to be equivalent to The value, Rb/Ra = 0.05, represents a weak yield threshold, and therefore the induced velocities are strong relative to when Rb/Ra is larger, and the region over which there is no flow is quite small. We also see that the flow strengthens as time progresses, although it eventually tends towards a steady state. When Rb/Ra = 0.2, which is just below the critical value above which flow is not possible, then the induced flow is weaker and the region of stagnation is much larger.
Given that the buoyancy force which is available to induce fluid motion is proportional to the temperature difference across the cylinder, it is clear that, for any chosen value of Rb, there must be an interval of time before that force overcomes the yield threshold and flow is induced. Therefore, there will an initiation time for flow. Figure 11 Fig. 12 where it may also be seen that the initiation time tends to infinity as Rb/Ra → 1 4 . That this is the maximum value for Rb/Ra may be shown easily by noting first that the maximum temperature difference within the cylinder varies from zero at t = 0 to 0.5 [see Eq. (22); Fig. 9]. Then, the equivalent of Eq. (16) is for convection to arise, and hence the largest value of Rb for which convection arises is Rb = Ra/4. Figure 13 shows the variation with time of the hydrostatic pressure gradient adjustment which is required in order to maintain a zero mean velocity, and which accounts for the change in the density of the fluid due to increasing mean temperature. Little detail may be seen when plotting the variation of p z with time due to the magnitude of the 2t term in the expression for θ in (20), and therefore the second frame in Fig. 13 shows how Ra −1 p z − 2t varies with time. The beginning point of each curve represents the initiation time for convection. In all cases, this modified hydrostatic pressure gradient adjustment tends to zero as time progresses, suggesting that p z ∼ 2Ra t for large times. Therefore, it is worth determining if  20) and (22) tend exponentially to zero as t increases. It is then quite a short analysis to show that precisely when exponentially terms are neglected. The given locations of the yield surfaces agree to more than eight significant figures with our numerical work when t = 1.

Conclusions
In this paper, we have considered where the outer impermeable surface of a porous circular cylinder is subject either to a sudden change in the surface temperature or in the applied heat flux. The porous medium within the cylinder is saturated with a Bingham fluid. In both cases a thin thermal boundary layer is initiated at the outer surface which conducts inwards. When a constant temperature is applied, flow is initiated immediately when Rb/Ra < 1 2 , but this flow eventually ceases as the fluid at the centre of the cylinder heats up, thereby providing a buoyancy force which is too small to overcome the yield threshold of the Bingham fluid. On the other hand, when the surface is subject to a constant heat flux, there is no flow at early times, again due to have too small a buoyancy force. However, flow eventually initiates if Rb/Ra < 1 4 , and the flowfield tends towards a steady state. In undertaking this analysis, we have assumed that the flow and temperature fields are independent of z, the vertical coordinate, and this is the situation which is likely to pertain when the cylinder is very tall compared with its radius; such a tall but finite height means that the overall flow must satisfy a zero vertical mean value.
The present analysis may now be extended to flows which arise outside a heated cylinder which is embedded within a porous medium, or to cases of vertical annuli where the porous medium is confined between the two surfaces, or even to cases where a volumetric heating is initiated suddenly. We believe that it is extremely unlikely that these flows will undergo a convective instability because the Newtonian counterpart has already been shown to be unconditionally stable in various separate analyses; see Gill (1969), Rees (1988, 2011), Straughan (1988, Lewis et al. (1995) and Scott and Straughan (2013). The effect of the threshold pressure gradient is to reduce the strength of the basic fluid flow, which generally means that flows are less susceptible to unstable disturbances. We note, however, that such channels will admit instability, but only when quite substantial alterations have been made from the configuration considered by Gill (1969). Thus Barletta (2015) found that instability can arise when the impermeable surfaces are replaced by constant pressure surfaces. Shankar and Shivakumara (2017) allowed the porous medium to be saturated by an Oldroyd-B fluid, thereby extending the work of Rees (1988), and found that this situation may also be subject to instabilities. It may be argued that these two respective changes to the Gill's configuration allow the fluid more freedom to be destabilised, whereas the present configuration represents less freedom.