Multiphase flow in porous channels and heat transfer have been studied more and more recently. This research area has large-scale potential in engineering and geophysical applications. Some major applications of such flows are in agricultural engineering for studying surface and underground water flows [10], in nuclear engineering for designing pebble bed reactors, in the petroleum industry for studying the flow of hydrocarbons in reservoir rocks, and in geotechnical engineering for underground waste disposal. Other applications include infiltration of water, sewage, porous bearings, solid matrix heat exchangers, and bioconvection in porous media. Many researchers [5, 8, 21] have considered blood flow as a two-phase flow.

Vajravelu et al. [32] analysed the hydromagnetic unsteady flow of two conducting immiscible fluids between two permeable beds with different permeabilities. They obtained expressions for interface velocity, velocity distribution, and mass flow rate. Subsequently, they analysed the pulsatile flow of a viscous fluid between two permeable beds using Darcy's law [33]. Analytical solutions of a flow in a porous medium inside permeable beds under an exponentially decaying pressure gradient were obtained by Prasad and Kumar [18]. Jogie and Bhatt [13] studied the laminar flow of two immiscible and incompressible fluids through permeable channels by using the B-J slip condition. Panda et al. [17] studied the three-layer fluid flow in a channel with a small obstruction on an impermeable bottom. Umavathi et al. [31] studied the unsteady flow in a porous medium sandwiched between viscous fluids. In the case of flow through porous mediums, Beavers and Joseph [3] demonstrated experimentally that the usual no-slip boundary condition is no longer valid for a fluid-porous boundary. They postulated the existence of slip at the fluid-porous boundary resulted in a condition called B-J slip. According to this condition, the Poiseuille velocity in a channel and Darcy's velocity in the porous medium can be coupled using the following equation:

Here, the clear fluid region occupies the region (y>0), u_f is the fluid velocity, and u_f and ∂uf∂y are evaluated at y=0+. The Darcy velocity u_m is evaluated at some small distance below y=0. The Beaver-Joseph constant α is dimensionless and depends on the structure of the porous medium and independent of the fluid viscosity and permeability.

Flow rheology is the main principle of lubrication theory. Two main theories explain the flow of liquids: classical continuum mechanics, which neglects the effects of fluid particle size, and microcontinuum mechanism, which accounts for the intrinsic motion of material particles (e.g. polymer molecules in a polymer suspension). Further, to study the discriminative behaviour of the fluids containing substructures such as polymer fluids, many theories of microcontinuum have been developed [1, 2, 28, 29]. Stokes [28] presented a simple generalisation of the classical theory of fluids, which allows for polar effects such as the presence of an antisymmetric stress tensor, couple stresses, and body couples. Couple stresses arise in fluids containing additives with large molecules. Major applications of couple stress fluid models include pumping fluids such as colloidal fluids, synthetic lubricants, liquid crystals, and biofluids (e.g. blood [26, 27]). Many researchers have applied the couple stress fluid model on lubrication problems such as squeeze film bearings [4, 20], thrust bearings [9, 14], and journal bearings [6, 16].

Soundalgekar and Aranake [23] studied the effects of couple stress on the magnetohydrodynamic Couette flow of a conducting fluid between two parallel plates. They concluded that current density is significantly affected by small values of couple stress parameter, and skin friction is affected by couple stress in the presence of a magnetic field. A brief discussion on the theory of microstructure fluids can be found in [29]. Szen and Rajagopal [30] analysed the flow of a non-Newtonian fluid between heated plates and obtained expressions by using constant and variable viscosity models. Recently, Siddiqui et.al. [22] studied Couette and Poiseuille flow of two non-Newtonian fluids, namely a fourth-grade fluid and Sisko fluid and concluded that solutions of the Newtonian and fourth-grade fluids can be recovered by substituting n=0 and 3 in Sisko model solutions. Islam and Zhou [11] derived exact solutions for the two-dimensional flow of couple stress fluids by using the inverse method. Farooq et al. [7] studied the non isothermal flow of a couple stress fluid with variable viscosity between two parallel plates. Iyengar and Bitla [12] investigated the pulsating flow of a couple stress fluid between permeable beds through constant injection and suction. Subsequently, Bitla and Iyengar [19] studied the oscillatory flow of a couple stress fluid inside permeable beds in the presence of a uniform magnetic field. Srinivas et al. [24] analysed entropy generation by the flow of two immiscible couple stress fluids between porous beds. The effects of Hall currents on the thermal instability of a compressible couple stress fluid in the presence of a horizontal magnetic field were investigated by Mehta et al. [15]. Srinivas and Murthy [25] studied the flow of two immiscible couple stress fluids between permeable beds and concluded that the presence of the couple stress reduces the flow velocity.

In view of the various applications of couple stress fluid flow in natural systems (e.g. in ground water flow, where multiple layers of fluids are present between permeable soil layers), human systems, and many engineering problems (e.g. in the petroleum industry to study the flow of many immiscible hydrocarbons between porous rocks), we analysed an MHD pulsating flow and heat transfer of two immiscible, incompressible, and conducting couple stress fluids between two permeable beds. The effects of different flow parameters on velocity and temperature profiles are displayed graphically and those on shear stress and heat transfer rates at permeable beds are presented numerically through tables.

Consider the flow of two electrically conducting immiscible couple stress fluids in a channel of height 2h bounded by two permeable beds. Permeable beds have infinite thicknesses with varying permeabilities. The permeabilities of lower and upper permeable beds are K₁ and K₂, respectively. The flow geometry is presented in Fig. 1. The origin is the centre of the channel, and X and Y are horizontal and vertical coordinates, respectively. The region I (−h≺_y≺_0) is occupied by an electrically conducting couple stress fluid of density ρ1, viscosity μ1, electric conductivity σ1, and thermal conductivity k₁. The region II (0≺_y≺_h) is filled with an electrically conducting couple stress fluid of density ρ2(≺ρ1), viscosity µ ₂, electric conductivity σ ₂, and thermal conductivity k₂. The lower and upper permeable beds are held at different constant temperatures Tw1 and Tw2, respectively, with Tw2≻Tw1. The field equations describing a couple stress fluid are similar to Navier-Stokes equations and are given by:

Figure 1. Schematic of flow

Conservation of mass

Conservation of momentum

Conservation of energy

Following Stokes, the constitutive equations of couple stress fluid are:

where material constants λ and µ are the viscosity coefficients and η and η and η' are the couple stress viscosity coefficients, which satisfy the constraints μ≽0,3λ+2μ≽0,η≽0, and ∣η'∣≺_η. ημ is the length parameter, which represents the characteristic measure of the polarity of the couple stress fluid and is identically zero for non polar fluids.

The scalar quantities P and ρ denote pressure at any point in the fluid and density of the fluid, respectively. The vectors q→, ω, f→, and l→ are the velocity, vorticity, body force per unit mass, and body couple per unit mass, respectively. The tensors τij and m_ij are force stress tensor and couple stress tensor, respectively. δ _ij is the Kronecker symbol, d_ij is the component of the rate of shear strain, ρck is body couple vector, w_i,j is the spin tensor, and ϵi,jk is the Levi-Civita symbol. Drs is the deformation tensor, which is equal to the symmetric part of the velocity gradient. E is the internal energy density per unit area, h_i is the influx of energy per unit area, q is the internal energy source density per unit area, and the comma denotes covariant differentiation.

Flow in both regions is assumed to be one-dimensional, laminar, and driven only by a pulsatile pressure gradient applied at the inlet of the channel. Assuming the permeable beds to be homogeneous and of infinite thickness so that Darcy's law can be applied with B-J slip condition at the fluid-porous interface, the velocity and temperature distribution in the channel are given by:

2.1. Velocity Distribution

The governing equation in Region I is

The governing equation in Region II is

To determine velocities u₁ and u₂, we adopt the following boundary and interface conditions:

• 1. At the lower fluid--porous boundary, couple stress disappears (i.e. no spin exists) and B-J slip condition is considered:

  (2.8)∂2u1∂y2=0    at y=-h

  (2.9)u1=u1'   aty=-h

  (2.10)∂u1∂y=αK1(u1'-Q1)at y=-h.

• 2. At the fluid-fluid interface, velocity, rotation, shear stress, and couple stress are continuous:

  (2.11)u1=u2, ∂u1∂y=∂u2%∂y, μ1∂u1∂y−η1%∂3u1∂y3=μ2∂u2∂y−η2∂3u2∂y3,η1∂2u1∂y2=η2∂2u2∂y2                     at y=0.

• 3. At the upper fluid--porous boundary, couple stress disappears (i.e. no spin exists) and B-J slip condition is considered:

  (2.12)∂2u2∂y2=0    at y=h

  (2.13)u2=u2'   at y=h

  (2.14)∂u2∂y=-αK2(u2'-Q2)    at y=h,n

  where u₁ and u₂ are velocities in Regions I and II, respectively, and B_o is the applied magnetic field normal to the flow direction. Q1=−K1μ1∂p∂x and Q2=−K2μ2∂p∂x are Darcy velocities in Regions I and II, respectively. ∂p∂x is the pulsatile pressure gradient given by:

  −∂p∂x=∂p∂xs+∂p∂xoeiωt,

  where ∂p∂xs and ∂p∂xo represent the steady and oscillatory parts of pulsatile pressure gradient, respectively, and ω is the frequency.

2.2. Temperature Distribution

The governing equation in Region I is

The governing equation in Region II is

To determine temperature profiles T₁ and T₂, we adopt the following boundary and interface conditions:

• 1. At the lower permeable bed, the temperature is constant and is equal to the temperature of the lower permeable bed, that is,

  (2.17)T1=Tw1   at y=-h.

• 2. At the fluid-fluid interface, the temperature and heat flux are continuous, that is,

  (2.18)T1=T2,  k1∂T1∂y= k2∂T2∂y    aty=0.

• 3. At the upper permeable bed, the temperature is constant and is equal to the temperature of the upper permeable bed, that is,

  (2.19)T2=Tw2   at y=h,

  where C_p is the specific heat at constant pressure.

Introducing the following non dimensional quantities

and dropping the asterisks, the velocity and temperature distributions in the non dimensional form are given as:

3.1. Velocity Distribution

Equations (2.6) and (2.7), respectively, become:

The boundary and interface conditions become:

where −∂p∂x=∂p∂xs+∂p∂xoeiωt is the non dimensional pressure gradient.

3.2. Temperature Distribution

Equations (2.15) and (2.16), respectively, become:

The boundary and interface conditions become:

Regarding the pulsating pressure gradient, let us assume that velocities and temperatures are in the form:

By using Equation (4.1) into the governing equations for velocity and temperature distributions and neglecting higher order terms, partial differential equations can be reduced into ordinary differential equations as follows:

4.1. Velocity Distribution

4.1.1. Steady Part

The governing equations of steady flow are:

The boundary conditions to be satisfied are:

4.1.2. Oscillatory Part

The governing equations of oscillatory flow are:

The boundary conditions to be satisfied are:

4.2. Temperature Distribution

4.2.1. Steady Part

The governing equations for steady flow are:

The boundary conditions to be satisfied are:

4.2.2. Oscillatory Part

The governing equations for steady flow are:

The boundary conditions to be satisfied are:

4.3. Velocity Profile

4.3.1. Steady Flow Solution

The solution for the steady flow described in Section 4.1.1 is given by:

where constants A_i, i=1,2,3,4 are provided in the appendix, whereas constants C_i, i=1,2,3,4,5,6,7,8 are not reported to achieve brevity.

4.3.2. Oscillatory Flow Solution

The solution for the oscillatory flow described in Section 4.1.2 is given by:

where constants A_i, i=5,6,7,8 are provided in the appendix, whereas constants C_i, i=9,10,11,12,13,14,15,16 are not reported to achieve brevity.

4.3.3. Pulsatile Flow Solution

The solution for the pulsatile velocity profile is given by:

where u11[y],u12[y],u21[y], and u₂₂[y] are known from the Equations (4.30), (4.31), (4.32), and (4.33), respectively.

4.4. Temperature Profile

4.4.1. Steady Flow Solution

The solution for the steady flow described in Section 4.2.1 is given by:

where constants A_i, i=5,6,7,8 are provided in the appendix, whereas constants C_i, i=9,10,11,12,13,14,15,16 are not reported to achieve brevity.

4.4.2. Oscillatory Flow Solution

The solution for the oscillatory flow described in Section 4.2.2 is given by:

where constants C_i, i=17 to 24 are not reported to achieve brevity, and all N_i's and M_i's are provided in the appendix.

4.4.3. Pulsatile Flow Solution

The solution for the pulsatile temperature profile is given by:

where θ11[y], θ12[y], θ21[y] and θ22[y] are known from Equations (4.35), (4.36), (4.37), and (4.38), respectively.

4.5. Mass Flux

The instantaneous mass fluxes are given by:

where Q₁ and Q₂ are mass fluxes, respectively, in Regions I and II.

4.6. Shear Stress

The shear stress in non dimensional forms at the permeable beds are given by:

4.7. Rate of Heat Transfer

The rates of heat transfer through the permeable beds to the fluid in the non dimensional form are given by:

In Section 4, we present analytical solutions for not only the pulsating velocity, temperature profile, and mass flux but also shear stress and heat flux at both permeable beds. The analytical expressions were evaluated numerically for different flow parameters values. The results are depicted graphically in Figs. 2-8 and numerically in tables. In the numerical evaluation, we take R2=ρ′μ′R1, S2=μ′η′S1, and ϵ1=ϵ2=ϵ.

Figure 2 shows the pulsating velocity (Figs. 2a and 2b) and temperature (Figs. 2c and 2d) profiles in both regions. The velocity corresponds to slip velocities at the interfaces of the lower (y=-1) and upper permeable beds (y=1). With one parameter kept fixed, Figs. 3a and 3b} depict the variation in flow velocity with respect to couple stress parameter S₁ and Hartmann number M. At M=0.4 and above, with an increase in S₁, flow velocity in both the regions decreased, whereas at M=0.3, flow velocity increased with an increase in S₁. Temperature profiles with respect to S₁ and M are presented in Figs. 3c and 3d. An increment in M or a decrement in S₁ resulted in an increment in temperature in both the regions. This was attributed to the fact that as S₁ decreases, couple stress velocity increases, which in turn enhances the internal heat. The temperature increases with an increase in M, and the increment in temperature occurs due to the Joule heating effect. Velocity and temperature profiles with respect to the slip parameter α and porosity parameter ϵ are presented in Fig. 4. The velocity in both the regions increased with an increase in α or a decrease in ϵ (Figs. 4a and 4b). The temperature in both the regions increased with an increase in α or a decrease in ϵ (Figs. 4c and 4d). At ϵ=1, the temperature profile was linear (Fig. 4d), whereas the variation in α appeared to exert no effect on the temperature. The effects of viscosity ratio coefficient µ' and couple stress viscosity ratio coefficient η' are depicted in Fig. 5. An increase in µ' or a decrease in η' resulted in a decrease in flow velocity because as µ' increases, fluids become thicker. A similar trend was noted for the temperature profile.

Figure 2. Velocity and temperature profiles with time t at ρ'=0.6,µ'=0.8,R₁=1,S₁=1,k'=1.2,M=1,σ'=1.2,Pr=1,Ec=1,η'=0.8,ϵ=0.5,α=0.6,P_s=1,P_o=1,ω=1.

Figure 3. Velocity and temperature profiles with time M and S₁ at ρ'=0.6,µ'=0.8,R₁=1,k'=1.2,σ'=1.2,Pr=1,Ec=1,η'=0.8,ϵ=0.5,α=0.6,P_s=1,P_o=1,ω
t = π 4 .

Figure 4. Velocity and temperature profiles with time α and ϵ at
ρ ′ = 0.6 , μ ′ = 0.8 , R 1 = 1 , S 1 = 1 , k ′ = 1.2 , M = 1 , σ ′ = 1.2 , P r = 1 , E c = 1 , η ′ = 0.8 , P s = 1 , P o = 1 , ω t = π 4 .

Figure 5. Velocity and temperature profiles with time µ' and η' at
ρ ′ = 0.6 , R 1 = 1 , S 1 = 1 , k ′ = 1.2 , M = 1 , σ ′ = 1.2 , P r = 1 , E c = 1 , ϵ = 0.5 , α = 0.6 , P s = 1 , P o = 1 , ω t = π 4 .

Figure 6 presents the variation in temperature with respect to Prandtl number Pr, Eckert number Ec, electric conductivity ratio σ', and thermal conductivity ratio k'. From Figs. 6a and 6b, it can be observed that as Pr increased, the temperature profile also increased. This increase is attributed to an increase in viscous diffusion in the presence of viscous dissipation, which enhances internal heat generation. As Ec increased, fluid frictional effects increased and hence the temperature profile increased. However, the temperature also increased with an increase in σ' (Fig. 6c), whereas it decreased with an increase in k' (Fig. 6d). From Figs. 7a and 7b, we observed that as the frequency parameter ω t increased, the flow velocity and temperature in both the regions decreased. Figure 8 presents the variation in velocity and temperature profiles when fluids in both the regions were considered identical. Figures 8a and 8b depict the velocity and temperature profiles, respectively, when M=0. The velocity profile resembles the Hagen--Poiseuille flow, and the temperature profile is almost linear. Further, the velocity at permeable beds becomes zero for ϵ1=ϵ2=ϵ (Fig. 8c).

Figure 6. Temperature profiles with Pr, Ec, σ' and k' at
ρ ′ = 0.6 , μ ′ = 0.8 , R 1 = 1 , S 1 = 1 , M = 1 , η ′ = 0.8 , ϵ = 0.5 , α = 0.6 , P s = 1 , P o = 1 , ω t = π 4 .

Figure 7. Velocity and temperature profiles with ω t at
ρ ′ = 0.6 , μ ′ = 0.8 , R 1 = 1 , S 1 = 1 , k ′ = 1.2 , M = 1 , σ ′ = 1.2 , P r = 1 , E c = 1 , η ′ = 0.8 , ϵ = 0.5 , α = 0.6 , P s = 1 , P o = 1.

Figure 8. Some particular cases.

Non dimensional shear stresses at lower and upper permeable beds were evaluated numerically, and the results are presented in Table 1. As ω t increased, shear stress at the lower permeable bed increased, whereas at the upper permeable bed, as ωt increased through values 0 to π4, shear stress decreased, and a further increase in ωt resulted in an increase in shear stress. From Table 1, it is clear that as ρ′,α,R1, and S₁ increased, shear stress at both the permeable beds also increased, whereas as µ', η', and ϵ increased, shear stress at both the permeable beds decreased. As M increased, shear stress at the lower permeable bed decreased, whereas shear stress at the upper permeable bed increased with an increase in M.

Numerical values of rates of heat transfer through both the permeable beds to the fluid were calculated for different values of governing flow parameters and are given in Table 2. As ωt,ϵ, and μ′ increased, the rate of heat transfer through the lower permeable bed decreased, whereas the rate of heat transfer through the upper permeable bed increased. As k′,η′,α,R1,Pr, and Ec increased, rates of heat transfer through the lower and upper permeable beds increased and decreased, respectively. As σ' increased, the rate of heat transfer through the lower permeable bed decreased, whereas as σ' increased through values 0.8 to 1.0, the rate of heat transfer through the upper permeable bed increased, and a further increase in σ' resulted in a decrease in the rate of heat transfer through the upper permeable bed. As S₁ increased, the rate of heat transfer through the lower permeable bed decreased, whereas as S₁ increased through values 0.3 to 0.4, the rate of heat transfer through the upper permeable bed increased, and a further increase in S₁ resulted in a decrease in the rate of heat transfer through the upper permeable bed. The rate of heat transfer through the lower permeable bed first increased and then decreased with an increase in M, whereas the rate of heat transfer through the upper permeable bed increased with an increase in M.

We analysed the MHD pulsating flow and heat transfer of two immiscible, incompressible, and conducting couple stress fluids between two permeable beds. The analytical solutions of the velocity profile, temperature profile, mass flux, shear stress, and rate of heat transfer were obtained. Analytical solutions were numerically evaluated for the different values of governing flow parameters. The effects of the flow parameters on flow velocity and temperature variation are depicted through graphs, whereas variations in shear stress and rate of heat transfer at the permeable beds are presented through tables. The following conclusions can be drawn from the entire analysis:

• • The velocity profile decreases as M increases. At M=0.4 or above, velocity is a decreasing function of couple stress parameter S₁, whereas at M=0.3 or below, velocity is an increasing function of S₁.

• • Both velocity and temperature increase as the slip parameter α increases, whereas for porosity parameter ϵ, both exhibit an inverse behaviour.

• • An increase in viscosity ratio µ' or a decrease in couple stress viscosity ratio results in a decrease in both the velocity and temperature profile.

• • The temperature profile is an increasing function of Prandtl number Pr, Eckert number Ec, and electric conductivity ratio σ', whereas increases in the thermal conductivity ratio k' have the tendency to cool down the thermal state.

• • When M=0, the temperature profile is nearly linear.

• • An increase in the density ratio ρ', slip parameter α, Reynolds number R₁, and couple stress parameter S₁, the shear stress at both the permeable beds increases. However shear stress shows a reverse trend for viscosity ratio µ', couple stress viscosity ratio η', and porosity parameter ϵ.

• • The rates of heat transfer at both the permeable beds are positive for the considered set of parameters, which indicates that heat is transferred to the permeable beds by the fluids. This occurs due to internal heat generation by viscous dissipation effects.