Asymptotic approach to the generalized Brinkman’s equation with pressure dependent viscosity and drag coefficient

Journal of Applied Fluid Mechanics, Vol. 9, No. 6, pp. 3101-3107, 2016. Available online at www.jafmonline.net, ISSN 1735-3572, EISSN 1735-3645. Asymptotic Approach to the Generalized Brinkman’s Equation with Pressure-Dependent Viscosity and Drag Coefficient I. Pažanin1†, M. C. Pereira2 and F. J. Suárez-Grau3 1 Department of Mathematics, Faculty of Science, University of Zagreb, Bijenička 30, 10000 Zagreb, Croatia 2 Departamento de Matemática Aplicada, Instituto de Matemática e Estatística, Universidade de São Paulo, Rua do Matão, 1010, CEP 05508-090, São Paulo, Brazil 3 Departamento de Ecuaciones Diferenciales y Análisis Numérico, Facultad de Matemáticas, Universidad de Sevilla , Avenida Reina Mercedes S/N, 41012 Sevilla, Spain †Corresponding Author Email: pazani[email protected] (Received October 26, 2015; accepted May 4, 2016) ABSTRACT In this paper we investigate the fluid flow through a thin (or long) channel filled with a fluid saturated porous medium. We are motivated by some important applications of the porous medium flow in which the viscosity of fluids can change significantly with pressure. In view of that, we consider the generalized Brinkman’s equation which takes into account the exponential dependence of the viscosity and the drag coefficient on the pressure. We propose an approach using the concept of the transformed pressure combined with the asymptotic analysis with respect to the thickness of the channel. As a result, we derive the asymptotic solution in the explicit form and compare it with the solution of the standard Brinkman’s model with constant viscosity. To our knowledge, such analysis cannot be found in the existing literature and, thus, we believe that the provided result could improve the known engineering practice. Keywords: Brinkman’s equation; Pressure-dependent viscosity; Pressure-dependent drag coefficient; Transformed pressure; Asymptotic analysis. NOMENCLATURE F dimensionless flux H thickness of the channel l length of the channel ˆ k pressure-viscosity coefficient k dimensionless pressure-viscosity coefficient ε M Characteristic number ˆ P pressure P dimensionless pressure dP dimensionless pressure drop 0 P referent pressure q transformed pressure ˆ u velocity u dimensionless velocity 0 V referent velocity ˆ β drag coefficient β dimensionless drag coefficient ε small parameter ˆ  viscosity  dimensionless viscosity σ auxiliary parameter 1. INTRODUCTION The steady-state flow of an incompressible, viscous fluid through a porous media is given by the conservation of mass and conservation of linear momentum principles. While the conservation of mass is described by the standard continuity equation ˆ div 0,u (1) different laws have been used to describe the conservation of the linear momentum. Without doubt, Darcy law (Darcy 1856) is the most commonly used model stating that the filtration I. Pažanin et al. / JAFM, Vol. 9, No. 6, pp. 3101-3107, 2016. 3102 velocity is proportional to the applied pressure gradient. However, it is based on several (restrictive) assumptions and, thus, its range of applicability is limited. In particular, expressed as a first order PDE for the velocity, Darcy law cannot sustain the no-slip boundary condition for the fluid velocity imposed on an impermeable wall. That inspired H. Brinkmann (Brinkman 1947) to propose the correction of the Darcy law in order to be able to impose such condition on an obstacle submerged in porous medium. In the absence of the external force, Brinkman equation can be written as ˆ ˆˆ ˆ ˆ β0,p   uu (2) Note that the second-order Eq. (2) can handle the presence of a boundary on which (physically relevant) no-slip condition is imposed. Nevertheless, in many geophysical problems the variations of the viscosity with pressure cannot be ignored if the flow is subjected to very high pressure drops. Such situation naturally occurs in petroleum engineering, namely in problems such as enhanced oil recovery and CO2 sequestration. In view of that, the Brinkman equation needs to be generalized in order to be able to capture the effects of the pressure-dependent viscosity. The notion that the fluid viscosity can depend on the pressure goes back to the celebrated work by Stokes (1845). Since then, numerous experimental investigations (see e.g. Binding et al. (1998), Goubert et al. (2001), Del Gaudio and Behrens (2009)) confirmed that, as the pressure is increased by several orders of magnitude, the variations of the viscosity with pressure should be taken into account while the flow is still incompressible. The viscositypressure relation is most commonly described by the Barus law (Barus (1893)) stating that the viscosity increases exponentially with pressure: ˆˆ 00 ˆ ˆˆ ˆ ( ) , , const. 0. kp pe k    (3) In the context of Eq. (2), the pressure-dependent viscosity implies that the porous medium parameter ˆ β (drag coefficient) also varies with the pressure (see e.g. Srinivasan et al. (2013)). More precisely, the exponential dependence (3) leads to a drag coefficient of the form ˆˆ 00 ˆ ˆˆ ˆ ( ) , , const. 0. kp pe k    (4) As a consequence, we arrive at the following generalized version of the Brinkman’s equation:     ˆ ˆˆ ˆ ˆ ˆˆ div 2 ( ) β0,ppp      Du u (5) where is   1 ˆˆˆ 2 T     Du u u the symmetric part of the velocity gradient tensor. The aim of this paper is to study the flow governed by Eqs. (1) and (5) from the analytical point of view. Note that when ˆ0   and ˆ β is constant, Eq. (5) reduces to simple Darcy equation, while for ˆ β0 and ˆ  constant, we get the Stokes equation. Finally ˆ  and ˆ β are assumed to be constants (i.e. k = 0), we have the Brinkman's equation (2). Introducing the exponential dependence (3), (4) into the Brinkman’s equation makes the problem highly nonlinear and, thus, we cannot expect to derive the exact solution of the full system (1), (5) even in the case of the simple two-dimensional channel (i.e. fracture with plane-parallel walls). Therefore, we introduce the small parameter ε (denoting the ratio between thickness of the channel and its length) and propose the asymptotic approach as ε → 0. After rewriting the Eqs. (1), (5) in the non-dimensional form, we apply the concept of the transformed pressure originally proposed by Marušić-Paloka and Pažanin (2013). By doing that, we transform the momentum Eq. (5) into the equation with small nonlinear perturbation that we can control. In such transformed system, we compare the characteristic (non-dimensional) number ε M describing the flow with small parameter ε trying to identify the critical case in which all the effects we seek for are balanced at the main order. We compute the asymptotic solution of the transformed problem in the most interesting (critical) case and then recover the solution of the governing problem by applying the inverse transformation. As a result, we obtain the effective model described by the explicit expressions for the velocity and pressure distribution. It enables us to clearly observe the influence of the viscositypressure dependance and porous structure on the effective flow. In particular, we can easily compare our asymptotic solution with the solution of the standard Brinkman’s model with constant viscosity. To conclude the Introduction, we provide some bibliographic remarks on the subject. For the flow through porous media, as governed by the Brinkman equation and its various generalizations, analytical treatments can be found only in the constant viscosity case ( ˆ0k ). We refer the reader to Durlofsky and Brady (1987), Kuznetsov (1998), Malashetty et al . (2001), Merabet et al . (2008), Marušić-Paloka et al . (2012), Khan et al . (2014)). In the variable viscosity case ( ˆ0k ), the numerical approach has been developed and we refer the reader to Naskhatrala and Rajagopal (2009), Srinivasan et al . (2013). In those papers, the results based on numerical simulations have been provided clearly suggesting that the corresponding solutions exhibit significantly different characteristics than the solutions of the classical Darcy or Brinkman’s equation. We confirm those findings in the present paper using analytical (asymptotic) approach and without making any simplifying assumptions on the governing system of PDEs and its solution. If we assume that û= û (ˆx) i (in fact, û (ˆx) = const., due to (1)), then the Eq. (5) is reduced to a simple ODE for the pressure which can be easily solved by the separation of variables (see Subramanian and Rajagopal (2007)). As far as we know, this is the first attempt to carry out such an I. Pažanin et al. / JAFM, Vol. 9, No. 6, pp. 3101-3107, 2016. 3103 analysis and, thus, we believe that it could improve the known engineering practice. 1If we assume that ˆu = ˆu (xˆ)i (in fact, ˆu(xˆ) = const., due to (1)), then the Eq. (5) is reduced to a simple ODE for the pressure which can be easily solved by the separation of variables (see Subramanian and Rajagopal 2007). 2. FORMULATION OF THE PROBLEM We consider the flow in a simple two-dimensional domain    ˆˆˆ ˆ ˆ ,:0,0. x yxlyh      2 R (6) The ratio ε = h/l is assumed to be small meaning that the channel under consideration is either very thin or very long. Such situation appears naturally in the applications we aim to address. As explained in the Introduction, the channel is filled with a fluid saturated porous medium and the flow is modeled by the generalized Brinkman’s equation with pressure-dependent viscosity and drag coefficient. Thus, we have the following system satisfied by the unknown fluid velocity ˆ u and pressure ˆ p :     ˆ ˆˆ ˆ ˆ ˆˆ div 2 ( ) β0,ppp      Du u (7) ˆ div 0,u (8) where the dependence of the functions ˆ  and ˆ β upon the pressure is given by Eqs. (3) and (4) respectively. To complete the problem, suitable boundary conditions have to be added. We impose a standard no-slip boundary condition for the velocity on the channel walls and we assume that the flow is governed by the prescribed pressure drop 0 ˆˆ l dP p p . In view of that, the boundary conditions read: ˆˆ 0 for 0, , y hu (9) ˆˆ 0 for 0, , x l ui (10) ˆˆ ˆ for , 0, . i pp xii l (11) Our goal is to investigate the asymptotic behavior of the flow described by (7)–(11), as ε → 0. 3. ANALYSIS 3.1 Equations in Non-Dimensional Form The first step is to write the governing problem in the non-dimensional form representing the appropriate framework for our analysis. We do it in a standard way by introducing 00 ˆˆ ˆ ˆ , , , , x yp xy p llVP     u u (12) 00 0 ˆˆ ˆβ , , β=ˆ 1ˆβ k k P    (13) where V0, P0 denote referent values of the velocity and pressure. Setting 0 0 0 β P Vl  (14) we obtain ε2() , kp M epMk p          ε uuDu (15) ε div 0.u (16) Note that the above equations are posed in the εindependent domain Ω = (0,1)2. The nondimensional characteristic number Mε appearing in the Eq. (15) represents the ratio between the frictional effect due to drag and the frictional effect due to viscosity, namely 2 0 0 ˆ ˆ l M     (17) Finally, hereinafter we use the following notation for the partial differential operators: 12 2 2 ε222 11 div , v , vv v v xy x y          v  εεε φ1φ1 , ( ) 2 T v xy          ε ijDv 11 1 ε22 1 , 1 2 vv xy vv +v vv xy                vij 3.2 Transformed Pressure Following the idea from Marušić-Paloka and Pažanin (2013), we introduce a new function qε, called the transformed pressure, such that εε . kp epq  (18) From (18) we deduce   1, . kkp qe e k     R (19) For the sake of the further analysis, it is important to observe that parameter σ in (19) can be chosen in an arbitrary way. Since εε ε 1, kp k p eq q ekq      (20) from (15)-(16) we get the following system satisfied by the velocity u and transformed pressure q: εεεε ε 2() , k k M qM q ekq           uuDu (21) I. Pažanin et al. / JAFM, Vol. 9, No. 6, pp. 3101-3107, 2016. 3104 ε div 0.u (22) The transformed Eq. (21) is still nonlinear but we are able control the nonlinearity appearing on the right-hand side. Indeed, the liberty in choice of parameter σ (see (19)) enables us to choose it small enough so that 2 lim 0. k k ekq       (23) Such choice of the parameter σ will be justified in the sequel by the fact that the effective pressure does not depend on σ at all (see (35)). That means that, throughout the analysis, σ plays the role of just an auxiliary parameter, i.e. by choosing σ such that (23) holds, we do not impose any additional constraints in the process. 3.3 Asymptotic Solution Now we have to construct the asymptotic solution of the transformed system (21)-(22). Before proceeding, it is important to make the following observation: if we kept characteristic number ε M constant (i.e. independent of ε), then a simply calculation would yield the model with no contribution of the porous structure at the main order. Thus, we need to compare ε M with small parameter ε trying to identify the critical case in which those effects remain in the macroscopic model. It can be easily verified that the critical case takes place when 2 ε(ε)MO  . Indeed, if we take 2 ε(ε)MO  we would obtain a simple Darcy law not accounting the effects of the Brinkman (viscous) term. On the other hand, the assumption 2 ε(ε)MO  would yield the situation already described above. In view of that, we set ε2, (1) ε M MMO (24) and perform a careful analysis under this assumption. Note that the assumption (24) suggest that the frictional effects between the fluid and the pores in the solid are dominant over the frictional effects within the fluid due to viscosity. Expanding the solution as 01 ε...  uu u 01 ε...qq q  and taking into account (23), we first conclude that 00 ()qqx. Next, we obtain 0 200 10 22 2 2 0 1: 0 ε 1: 0 , ε 0 for 0,1. dq dq MMM dx dy y uin y y               uiju u (25) We deduce 1 00 11 (, ) , ()uxy q qxui , leading to 21 1 00 0 2= in . udq Mu M dx y      (26) We can treat (0,1)xas a parameter and solve the above equation as a linear second-order ODE with respect to y. We get 10 01 2 (, ) () () . My My dq uxy Cxe Cxe dx    (27) The functions 12 (), ()CxC xcan be computed by taking into account that 1 00u for 0,1.y We obtain 00 12 11 () , () . MM MM MM dq dq ee Cx C x dx dx ee ee      (28) It remains to determine the transformed pressure q0. From the divergence equation we deduce 12 01 1 : 0 in . uu xy     (29) Integrating the above equation with respect to y gives 11 11 001 00 0 const.udy udy D x        (30) Taking into account (27), we obtain 1 01 2 2 () . 22 MM MM Me e M qx D x D ee           (31) In view of (19), the boundary conditions for q0 are given by   0 1 () , 0,1, i kp k qi e e i k      (32) with 0 01 00 ˆˆ ,l p p pp p p  . Consequently, we arrive at     00 1 0 11 () . kp kp kp k qx e e x e e kk     (33) In order to reconstruct the effective pressure eff p related to the governing (dimensionless) problem2, we apply the inverse transformation (see (19)) to obtain 0 11 () ln . () eff k pxkekqx        (34) In view of (33) we deduce  00 1 11 () ln . eff kp kp kp pxkee ex         (35) The effective velocity is provided as I. Pažanin et al. / JAFM, Vol. 9, No. 6, pp. 3101-3107, 2016. 3105 1 0 () (), eff yuyui where 1 0 u is given by (27).Since 01 ,pp it follows 0 1 kp kp ee   implying that eff p is well-defined. Moreover, observe that the effective pressure does not depend on the parameter σ at all. This fact justifies the above transformation procedure, namely the choice of the parameter σ such that (23) holds. 4. D ISCUSSION Let us write our (dimensionless) asymptotic solution. In view of (27), (28), (33), and (35), it has the following form: 01 1 0 1 0 () (), () [(1 ) () +( 1) ( )], eff kp kp MMy MM MMyMM yuy ee uy e e ke e ee ee         ui (36) 00 1 11 ln , () eff kp kp kp pkeeex        (37) for ,(0,1)xy On the other hand, it can be easily verified that the solution corresponding to the standard Brinkman model (1)-(2) with constant viscosity reads: 10 () (), () [(1 ) () +( 1) ( )], BB MMy BMM MMyMM yuy pp uy e e ee ee ee         ui (38) 10 0 () ( ) . b px p px p  (39) The above velocity and pressure profiles have been plotted in Figs. 1-2 for different values of the pressure-viscosity coefficient k. According to the experimental results that one can find in the literature (see Srinivasan et al. (2013) for details), it is reasonable to choose values  0.01,0.02,0.035k in the case of the referent pressure 6 0 10 P  . In the following we also take M = 1 and put 0 300p and 1 1p ensuring the setting with high (dimensional) pressure drop. From Fig. 1 we conclude that the velocity profiles exhibit different characteristics in constant ( k = 0) and variable viscosity case ( k > 0). Indeed, though the symmetry (around y = 0.5) is preserved in both settings, the values of the velocity for k > 0 are more than 200 percent lower than those for k = 0. Furthermore, as the pressure-viscosity coefficient k increases, the velocity profile becomes more flattened indicating that the velocity, for higher values of k, does not vary too much with respect to y . Fig. 1. Plots of (0) B uk and 1 0 u for various values of k > 0. Fig. 2. Plots of (0) B pk and eff p for various values of k > 0. Fig. 3. Plots of F B (k = 0) and F for various values of k > 0. The difference between constant and variable viscosity case is even more apparent when we compare the pressures (see Fig. 2). While for k = 0 we have a simple linear function describing the pressure distribution, for k > 0 the outcome is completely different. We observe that the pressure exhibits a huge decrease within a short distance near the left end of the channel. After that it does not vary too much with respect to x. That suggests the appearance of the pressure boundary layer in the vicinity of x = 0 in case of high pressure drops between the channel’s ends. Note that the decrease near x = 0 becomes more significant as the pressureviscosity coefficient k increases. Last but not least, it is interesting to investigate how the corresponding flow rate (flux) changes with respect to the pressure drop 01 (1, 300)dP p p By a simple integration, for Brinkman’s model (k = 0), we obtain that the flux increase linearly with the I. Pažanin et al. / JAFM, Vol. 9, No. 6, pp. 3101-3107, 2016. 3106 change of the pressure difference across the channel, namely 111 01 0 () () ( 3 4). B dP FdP uydy e e ee      (40) On the other hand, in the variable viscosity case, from (36), we deduce 11 1 0 00 300 (300 ) 1 1 () . (), = ( 3 4). () eff kkdP FdP u dy u ydy ee ee ke e        i (41) As we can observe from Fig. 3., the values of the flux are significantly lower for k > 0 than those for k = 0 so we may conclude that the standard Brinkman’s model in some sense overpredicts the flux. Moreover, the flux in variable viscosity case increase very slowly, with the increase of dP. In fact, the increase in the flux is negligible up to some (critical) value of the pressure drop dP (being rather high even for moderate values of k). Such phenomena becomes more dominant, as the value of the pressure-viscosity coefficient k goes up. For instance, for k = 0.035, the critical value of the pressure drop is almost 250, i.e. the variations of the flux are negligible for pressure drops dP < 250. 5. CONCLUDING REMARKS In the present paper we address one important application of the porous medium flow: the fluid flow governed by the high pressure drop through a thin channel (i.e. thin fracture with plane-parallel walls) filled with fluid saturated porous medium. Such flows appear naturally in the industrial applications, in particular in petroleum engineering. The framework with high values of pressure induce the significant variations of the viscosity with respect to pressure making standard Brinkman’s model with constant viscosity inappropriate for describing such flows. Thus, we consider the generalized Brinkman’s model with pressuredependent viscosity and drag coefficient. The viscosity and drag-pressure relation is assumed to be exponential (Barus law) and the effective flow is found using asymptotic approach with respect to the thickness of the channel. The key idea is to introduce the notion of the transformed pressure enabling us to control the nonlinearity in the momentum equation. The obtained result clearly indicates that the asymptotic solution exhibits utterly different characteristics than the solution corresponding to the classical constant viscosity case. To our best knowledge, this is the first attempt to solve this problem using analytical approach and without making any simplifications on the starting problem. Finally, one of the benefits of the analysis presented here lies in the fact that it can be straightforwardly generalized in two directions. Instead of addressing simple fracture with planeparallel walls, we can consider more complex (and realistic) domain, namely the fracture with constrictions (see e.g. Gipouloux and Marušić- Paloka (2002)). Besides the effects of the porous structure and pressure-dependent viscosity, in such setting we are particularly interested to detect the effects of the shape functions describing the constrictions. Another possible generalization would be to consider general viscosity-pressure dependence ˆˆˆ () p p  satisfied by Barus law and other empiric laws (see e.g. Marušić -Paloka and Pažanin (2013)). Instead of (19), we would simply introduce the transformed pressure q as ξ () (ξ) pd Bp  and continue with the same procedure as above. ACKNOWLEDGMENTS The first author of this work has been supported by the Croatian Science Foundation (scientific project 3955: Mathematical modeling and numerical simulations of processes in thin or porous domains. The second author has been partially supported by FAPESP 2013/22275-1, CNPq 302960/2014-7 and 471210/2013-7, Brazil. The third author has been partially supported by the projects MTM2014- 53309-P of the Ministerio de Economı´a y Competitividad and FQM309 of the Junta de Andalucı´a. The authors would like to thank the referees for their helpful comments and suggestions that helped to improve the paper. REFERENCES Barus, C. (1893). Isothermals, isopiestics and isometrics relative to viscosity. AM. J. Sci. 45, 87–96. Binding, D. M., M. A. Chouch and K. Walters (1998). The pressure dependence of the shear and elongational properties of polymer melts. J. Non-Newtonian Fluid Mech. 79, 137-155. Brinkman, H. (1947). A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles. Appl. Sci. Res. A1, 27–34. Darcy, H. (1856). Les fontaines publiques de la ville de Dijon, Victor Darmon, Paris. Del Gaudio, P. and H. Behrens (2009). An experimental study on the pressure dependence of viscosity in silicate melts. J. Chem. Phys. 13. Durlofsky, L. and J. F. Brady (1987). Analysis of the Brinkman equation as a model for flow in porous media. Phys. Fluids 30, 3329-3341. Gipouloux, O. and E. Marušić-Paloka (2002). Asymptotic behavior of the incompressible Newtonian flow through thin constricted fracture. Multiscale Problems in Science and Technology, Springer 189–202. Goubert, A., Vermant, J., P. Moldenaers, A. Gotffert and B. Ernst (2001). Comparison of measurement techniques for evaluating the pressure dependence of viscosity. Appl. Rheol. 11, 26-37. I. Pažanin et al. / JAFM, Vol. 9, No. 6, pp. 3101-3107, 2016. 3107 Khan, W., F. Yousafzai, M. I. Chohan, A. Zeb, G. Zaman and I. H. Jung (2014). Exact solutions of Navier-Stokes equations in porous media. Int. J. Pure Appl. Math. 96, 235–247. Kuznetsov, A. S. (1998). Analytical study of fluid flow and heat transfer during forced convection in a composite channel partly filled with a Brinkman Forchheimer porous medium. Flow, Turbulence and Combustion 60, 173–192. Malashetty, M. S., J. C. Umawathi, J. P. Kumar (2001). Convective flow and heat transfer in an inlined composite porous medium. J. Porous Media 4, 8. Marušić -Paloka, E. and I. Pažanin (2013). A note on the pipe flow with a pressure-dependent viscosity. J. Non-Newtonian Fluid Mech. 197, 5–10. Marušić-Paloka, E., I. Pažanin and S. Marušić (2012). Comparison between Darcy and Brinkman laws in a fracture. Appl. Math. Comput. 218, 7538–7545. Merabet, N., H. Siyyam and M. H. Hamdan (2008). Analytical approach to the Darcy-Lapwood- Brinkman equation. Appl. Math. Comput. 196, 679–685. Naskhatrala, K. B. and K. R. Rajagopal (2009). A numerical study of fluids with pressure dependent viscosity flowing through a rigid porous media. Int. J. Num. Meth. Fluids 67, 342-368. Srinivasan, S., A. Bonito and K. R. Rajagopal (2013). Flow of a fluid through a porous solid due to high pressure gradients. J. Porous Media 16, 193-203. Stokes, G. G. (1845). On the theories of the internal friction of fluids in motion, and of the equillibrium and motion of elastic solids. Trans. Camb. Philos. Soc. 8, 287-305. Subramanian, S. C. and K. R. Rajagopal (2007). A note on the flow thorugh porous solids at high pressures. Comput. Math. Appl. 53, 260–275.