Pulsating flow in a channel filled with a porous medium under local thermal non-equilibrium condition: an exact solution

The present work investigates analytically the problem of forced convection heat transfer of a pulsating flow, in a channel filled with a porous medium under local thermal non-equilibrium condition. Internal heat generation is considered in the porous medium, and the channel walls are subjected to constant heat flux boundary condition. Exact solutions are obtained for velocity, Nusselt number and temperature distributions of the fluid and solid phases in the porous medium. The influence of pertinent parameters, including Biot number, Darcy number, fluid-to-solid effective thermal conductivity ratio and Prandtl number are discussed. The applied pressure gradient is considered in a sinusoidal waveform. The effect of dimensionless frequency and coefficient of the pressure amplitude on the system’s velocity and temperature fields are discussed. The general shape of the unsteady velocity for different times is found to be very similar to the steady data. Results show that the amplitudes of the unsteady temperatures for the fluid and solid phases decrease with the increase in Biot number or thermal conductivity ratio. For large Biot numbers, dimensionless temperatures of the solid and fluid phases are similar and are close to their steady counterparts. Results for the Nusselt number indicate that increasing Biot number or thermal conductivity ratio decreases the amplitude of Nusselt number. Increase in the internal heat generation in the solid phase does not have a significant influence on the ratio of amplitude-to-mean value of the Nusselt number, while internal heat generation in the fluid phase enhances this ratio.

Height of the channel (m) j Counter K Permeability of the porous medium ( m 2 ) k The ratio of fluid effective thermal conductivity to that of the solid k f Thermal conductivity of the fluid ( W m −1 K −1 ) k f, eff Effective thermal conductivity of the fluid, k f k s Thermal conductivity of the solid ( W m −1 K −1 ) Dimensionless steady bulk mean temperature of the fluid

Introduction
Convective heat transfer in porous media has been a subject of intense studies due to its wide range of application in the industries such as oil recovery, geothermal engineering, thermal insulation, carbon storage, heat transfer augmentation, solid matrix or micro-porous heat exchangers and porous radiant burners [1,2]. Studies regarding the thermal characteristics of non-pulsating flow in porous media are more concentrated on heat transfer enhancement in domains filled with porous materials subjected to a heat source at the wall boundaries or internal heat generation. There were abundance of experimental (e.g., [3][4][5][6][7]), numerical (e.g., [8][9][10][11][12]) and theoretical (e.g., [13][14][15][16][17][18][19][20][21][22]) studies, which demonstrated the use of porous material as a promising passive technique in enhancing forced convection heat transfer in different industrial applications in micro-and large scales. Investigations regarding heat transfer of pulsating flow have mostly been conducted in empty (non-porous) channels and pipes (e.g., [23][24][25]). Understanding the fluid flow and heat transfer of pulsating flow in porous media has a pivotal role to play in biological applications such as blood flow in vessels due to heart beating and also industrial applications such as mesh-type regenerators used in the stirling cycle devices [26][27][28] and electronic cooling by utilizing oscillating flow [29,30].
In theoretical modeling of heat transfer in porous media, two primary models are generally used. The local thermal equilibrium (LTE) and local thermal non-equilibrium (LTNE) models. The LTE model holds when the heat exchange between the solid and fluid phases is high, and the temperature difference between the two phases is negligibly small [11]. This model requires solving only one energy equation to predict temperature field within the porous medium, which simplifies the analysis of heat transfer in porous media (e.g., [10,[31][32][33][34][35]). LTNE model however requires solving separate energy equations for the two phases in the porous regions, which are coupled through an internal heat exchange term. Hence, the LTNE model leads to more accurate prediction of the temperature fields in porous media (e.g., [8,20,[36][37][38]). Guo et al. [39] studied numerically heat transfer of pulsating flow based on LTE model in a tube partially filled with a porous medium attached to the pipe walls. They presented and discussed the relationship between the effective thermal diffusivity and thickness of the porous layer. Byun et al. [40] conducted an analytical characterization of heat transfer of oscillatory velocity flow in a large slab of a porous medium with a hot and a cold side using the LTNE model. They defined a criterion for validity of LTE model and reduced the solutions into four regimes of asymptotic solutions. Kuznetsov and Nield [41] provided analytical expressions for pulsating flow and forced convection heat transfer, produced by an applied oscillating pressure in a channel/pipe utilizing LTE model. They [41] found that the fluctuating part of the Nusselt number increases to a maximum value and then decreases with the increase in frequency. This observation was not in agreement with other studies that investigated pulsating flow in an empty channel or tube [24,25,42]. Forooghi et al. [42] performed a numerical investigation using the LTNE model for both steady and pulsating flow in a parallel-plate channel partially filled with porous layers attached to the channel walls. They found that an increase in the thermal conductivity ratio between the two phases, or amplitude of the pressure gradient, results in an enhancement in the dimensionless average Nusselt number. Yang and Vafai [36] discussed that the LTE model is not suitable to use for transient heat transfer in porous media. They [36] discussed that the temperature difference between the two phases is relatively small when the process reaches its steady condition. However, during the transient process the temperature difference between the two phases is considerable. Their [36] results revealed that utilizing LTE model for time-dependent problems of heat transfer in porous media induces certain inaccuracies in predicting temperature field [36]. From the previous studies, it could be noted that the problem of pulsating flow in a channel/ pipe filled with a porous medium under local thermal nonequilibrium condition has not yet been fully understood. The current work presents an analytical solution to investigate the effects of pulsating flow on the velocity, temperature distributions and Nusselt number in a channel filled with a porous medium under LTNE condition and considering internal heat generation in the fluid and solid phases. There are several practical examples involving internal heat generation in porous media such as electronic cooling, agricultural product storage [13,43] and metabolic reactions in biological media [21]. The problem of pulsating flow in a medium leads to certain time-dependent governing equations which need to be addressed theoretically or numerically. As discussed above, some attempts have been made to cope with such a problem but in different geometries. Some works addressed unsteady governing equations for stretching permeable sheets using HAM (Homotopy Analysis Method) [44][45][46][47]. Some other interesting analytical solution for different applicable mathematical problems could be found in [48][49][50][51][52][53][54][55][56].
This problem is worth investigating analytically, since the unsteady pulsating fluid flow and the heat transfer between particles and fluid in an unsteady pulsating flow is complex and clearly is expensive to study experimentally. Furthermore, for realistic porous systems, pore-scale modeling of porous systems is computationally prohibitive and hence deploying the volume-averaged method is recommended. In the volume-averaged method, the local thermal non-equilibrium (LTNE) model is deployed to calculate temperatures of the fluid and solid phases in the porous media (as considered in this study) by solving different energy equations for each phase in the media [57]. However, for the volume-averaged solution based on the LTNE model, the internal heat transfer coefficient between the particles and fluid has to be known a priori. This coefficient is required for the term, which couples the two energy equations of the particle and fluid in the porous media. Therefore, the present work aims to shed some light on the problem. The subject is also of interest as a basic research of unsteady forced convection problem.
The Darcy-Brinkman flow model is used to represent the fluid flow in the porous medium [16] and the LTNE model is employed to find exact solutions for the temperature distributions and Nusselt number in the system. The effect of parameters such as thermal conductivity ratio, Biot number, Darcy number, dimensionless frequency, coefficient of pressure amplitude and Prandtl number on the flow and heat transfer characteristics are presented and discussed.

Problem description and governing equations
The schematic diagram of the problem is shown in Fig. 1. Fluid flows through a parallel-plate channel filled with a porous medium subjected to a constant heat flux boundary condition. Heat generation in the solid phase, S s and the fluid phase, S f are considered uniform but different [13]. The distance between the plates is 2H and the heat flux q w is applied to the channel walls. The incoming flow has pressure gradient, which oscillates in time about a non-zero mean value. Following assumptions are considered to simplify the governing equations [13,16] Based on these assumptions, the governing equations are represented as [13,16]: Momentum [57] Momentum equation is the sum of unsteady Darcy equation f u t = − p x − K u and Brinkman term eff 2 u y 2 , which is unsteady Brinkman-extended Darcy model.
The pressure gradient is considered to vary with time in a sinusoidal waveform about a constant steady value: where ( P x ) st is the steady component of pressure gradient [23], is coefficient of pressure amplitude and is oscillation frequency. This is a known form for the pressure gradient to represent pulsating flow, which was also used in previous works (e.g., [25,58]).
Using Eq. (2), the momentum Eq. (1) is rewritten as: Energy equation for the fluid phase is expressed as: (1) Energy equation for the solid phase is written as: where subscripts f and s represent the fluid and solid phases, respectively. Subscript st refers to steady flow. is the dynamic viscosity of the fluid and eff = ∕ [59] is the effective viscosity of the porous medium. K is the permeability of the porous medium, ρ is density and c p is the specific heat. k f,eff and k s,eff are the effective thermal conductivity of the fluid and solid phases, respectively. a sf is the interfacial area per unit volume of the porous medium and h sf is the fluid-to-solid heat transfer coefficient [16].

Boundary conditions
The following boundary conditions are employed to solve the systems of the governing Eqs. (1)-(5) [13,16]: No-slip condition at the channel wall: Symmetry boundary condition is applied at the channel centerline: When the channel wall has a high thermal conductivity with a finite thickness attached to a porous medium, the temperature of the solid and the fluid phases are almost equal at the wall [60,61]. Using this assumption at the channel wall, Model A boundary condition is adopted at the interface between the porous medium and the channel wall [13,16]. This model assumes that the prescribed heat flux at the wall is split between two phases relative to the physical values of their effective thermal conductivities and temperature gradients. This model further assumes that the steady component of the temperature of the fluid and solid phases at the wall are equal to the wall temperature [13,16]: u| y=±H = 0.

Normalization
The following dimensionless variables are introduced to normalize the governing equations and boundary conditions [16,25,62]: is a characteristic velocity. = f is the kinematic viscosity of the fluid and Re is Reynolds number. w * is Womersly number and is the porosity of the porous medium. Using the dimensionless variables, the dimensionless form of momentum Eq. (3) and the associated boundary conditions Eqs. (6) and (7) are written as: Energy Eqs. (4) and (5) and their associated boundary conditions Eqs. (8) and (10) are also rewritten as: where T * f and T * s are defined as: The dimensionless variables used in Eqs. (15) and (16) are as: where Pr is Prandtl number and Ps is a dimensionless variable similar to Prandtl number appeared in the normalization process of solid-phase energy equation, Eq. (5). Bi is Biot number. k f,eff = k f and k s,eff = (1 − )k s are the effective thermal conductivity of the fluid and solid phases, respectively [16,59] and k is the thermal conductivity ratio.

Analytical solution for the momentum equation
To solve the momentum Eq. (12), it is divided into steady and unsteady components [25]: where subscripts st and un refer to steady and unsteady terms, respectively. Using Eq. (21) and considering U st t * = 0, the governing Eq. (12) and the boundary conditions (13) and (14) are written as: and The initial condition for the velocity is considered zero for convenience. However, the correct initial value is obtained by applying the fully developed assumption. The effect of initial condition on solution is discussed later in the paper. Therefore, the initial condition for the momentum equation is written as: The steady velocity Eq. (22) has been solved in the previous studies (e.g., [63]) and here only the final solution is presented, which is as: where The unsteady momentum Eq. (24) is a non-homogeneous equation with homogeneous boundary conditions (Eq. 25). Hence, it is solved using the method of eigenfunction expansion [64]. Therefore, the unsteady component of the velocity is given by Eq. (29). The procedure of solving the unsteady velocity using this method is explained in "Appendix" (Sect. 5.1). where

Analytical solution for the energy equations
The analytical solution of the energy equations is explained in "Appendix". The distribution of the temperatures in this section is presented in a form of = T * − T * w , which is used to calculate the Nusselt number. See "Appendix" (Sect. 5.2) for more details.

Steady components of the energy equations
The procedure of solving the steady energy equations and finding key parameters are demonstrated in "Appendix" (Sect. 5.2.1). The steady equations of the problem have already been studied and discussed in the previous studies. Kun and Vafai [13] solved the equivalent steady problem using the Darcian flow model. The focus of the present study is on the unsteady solutions and to prevent repetitions, solutions presented by Ref. [13] are provided here.
where (31) f,st = f,st and s,st are the dimensionless steady temperatures of the fluid phase and solid phase, respectively, and T w,st is the steady component of the wall temperature.

Unsteady components of energy equations
The unsteady components of the temperature are normalized using the unsteady wall temperature (in the form of ) to be comparable with the steady components: and where f,un and s,un are the dimensionless unsteady temperatures of the fluid phase and solid phase, respectively, and T w,un is the unsteady component of the wall temperature. F m,n (t * ) and R m,n (t * ) are time-dependent coefficients (Eqs. 95, 97 in "Appendix") and r 1 and r 2 are roots of a characteristics equation presented by Eq. (89) in "Appendix". (34) and (35) is a constant parameter obtained when solving the steady energy equations (see Sect. 5.2.1 of "Appendix" for more details):

Calculation of Nusselt number (Nu)
The wall heat transfer coefficient is defined as [13]: with the Nusselt number obtained as [13]: where 4H is the hydraulic diameter of the channel and f,b is the dimensionless bulk mean temperature of the fluid. Con- f,st,b is obtained using the following relation [16].
Mahmoudi [16] obtained an analytical solution for the equivalent steady equations considering the velocity slip and temperature jump at the channel walls. The solution for f,st,b and consequently, Nu st , for the present work can be obtained from the analytical solutions presented in [16] by setting the velocity slip coefficient and the temperature jump coefficient to zero.
Using Eqs. (29), (30) and (34), f,un,b is obtained as: where and where a j (t * ) is obtained by replacing subscript j instead of n in Eq. (30).

Validation
In this section, we present the validation of the unsteady velocity field in comparison with the analytical solutions of Siegel and Perlmutter [23] presented for pulsating flow in a channel without porous medium. According to Eq. (12) when the Darcy number is high enough the resulted momentum equation will be similar to that in a channel without a porous medium. The applied pressure gradient for the results obtained by Siegel and Perlmutter [23] was in the Cosine waveform [23]. By substituting t * + 2 instead of t * in Eq. (30), the unsteady velocity at the channel centerline (Y = 0) for = 0.1 , = 1 and M = 1 obtained from the present solutions is compared against those presented in [23] and shown in Fig. 2. An excellent agreement is observed between the two solutions. To the best of our knowledge, there is no closely relevant work in the literature on unsteady temperature field in a pulsating flow to be considered for validation of the temperature solution presented in this work. Figure 3 shows the unsteady velocity profile for = 0.5 , = 20 , M = 1.1 at Y = 0 . The velocity profile for the fully developed solution, i.e., the exponential term exp − (2n − 1) 2 2 M∕4 + 1∕Da t * in Eq. (30) is not considered. It is seen that the initial transient part decays after a short period of time [23] (here after t * > 0.2). The initial condition was considered zero for convenience. Considering  a different initial condition leads to a different exponential term, which decays in a short time. Since only the fully developed oscillatory part of the solution is of interest, for the rest of the results presented here, the initial transient part is not considered [23]. The effect of Darcy number (Da) on the dimensionless unsteady part of the velocity is shown in Fig. 4

Unsteady velocity profile
It is seen that the velocity amplitude has a direct relationship with Da number. Increasing the Da number by a factor ten increases the amplitude of the velocity field by almost the same factor. Increasing Da number is equivalent to increasing the permeability of the porous medium for a fixed channel height. This enhancement results in decrease in the resistance against flow (the term − U un Da in Eq. 24), which leads to a higher velocity magnitude. It can be seen from Eq. (27) that the dimensionless steady component of the velocity has also a direct relationship with Da number. Figure 5 shows the unsteady velocity versus dimensionless time for three values of dimensionless frequency ( ), for two Darcy numbers of 10 −3 and 10 −1 , with flow conditions of = 0.5 , M = 1 at Y = 0 . It is seen in Fig. 5a that for high Da numbers, the velocity amplitude decreases by increasing the frequency. While for low Da values, the effect of frequency on the velocity amplitude is negligible. If ( Da → 0 ) in Eqs. (29) and (30) the dimensionless unsteady velocity for low Darcy numbers can be written as: From Eq. (45), it can be seen that for low Darcy numbers a change in the value of only changes the period of oscillation (as sin ( t * ) ) and does not change the amplitude of unsteady velocity. Figure 6 shows variation of unsteady velocity component with Y for different dimensionless times during a dimensionless period of oscillation ( * = 2 ) for = 0.5 , = 5 , M = 1.1 and Da = 10 −3 along with the corresponding steady velocity. As expected, the velocity at the wall ( Y = ±1 ) is zero and the symmetry condition at the channel centerline is satisfied for all the profiles. The shapes of the unsteady velocity profiles are very similar to that of the steady component.

Unsteady temperature
The results presented in this section are obtained using the mean properties of water as the fluid phase and steel as the solid phase in the porous medium. The thermal properties used are M = 1.1 , k = 0.3 , Pr = 7.7 , Ps = 2.5 and the porosity of the porous medium is = 0.9, for all the results except for the cases mentioned in the text. In addition, similar to the unsteady velocity, here only the fully developed oscillating part of the temperature solution is of interest. Thus, the initial transient part of the temperature profile will not enter the solution. To achieve this, the exponential terms (of the coefficients F m,n (t * ) and R m,n (t * ) ) in Eqs. (34) and (35) are neglected to obtain the fully developed solutions. As an  Figure 8 shows the dimensionless steady and unsteady temperature for the fluid and solid phases as a function of Y at different times without internal heat generations. For both phases, the steady parts are negative meaning that f and s are lower than the wall temperature. The steady components of f and s are equal to the wall temperature at the channel wall. Away from the wall, the difference between the fluid and solid temperature with the wall temperature increases and reaches to its maximum at the channel centerline. Figure 8 further shows that the unsteady components of f and s oscillate around zero near the channel wall, while the unsteady solid temperature s has a more uniform distribution.
For the unsteady components of the dimensionless temperature distributions, the internal heat generation appears in constant C in Eq. (36). For steady-state Darcian flow  (31) and (32) show, only the internal heat generation in the solid phase ( S * s ) influences the dimensionless temperature distributions s,st and f,st , and S * f has no effect on them. While, for unsteady temperature components as shown in Eq. (36), the sum of the uniform internal heat generation in the solid and fluid phases ( S * ) influences the variation of s,un and f,un . Figure 9 shows the variation of It is seen that the amplitude of the unsteady dimensionless temperatures for the fluid phase is relatively bigger than the solid phase. Additionally, the dimensionless temperatures of the two phases of all cases decrease with increasing Bi number. Furthermore, the graphs show that for large value of Bi , which translates to a strong internal heat transfer between the fluid and solid phases, the difference between s,un and f,un is relatively small, which is in agreement with the results already presented for the steady flow [13]. It seems that the LTE model is valid for large Biot numbers in unsteady flows. For large Biot numbers, the amplitudes of unsteady temperatures are very small (close to zero). Therefore, the total value of the dimensionless temperature for the two phases is close to the steady flow. In fact, for large Biot numbers, the effect of unsteadiness on heat transfer decays largely. The figures also illustrate that by increasing the thermal conductivity ratio k , the amplitude of s,un and f,un decreases and also become similar. Figure  and Ps = 2.5 at ( Y = 0 ). The general trend for the two phases is similar. It is seen that the amplitudes of s,un and f,un increase with the increase of Pr number. Prandtl number is defined as the ratio of the momentum diffusivity to thermal diffusivity of the fluid [65]. Hence, increase in the fluid Pr number enhances the influence of heat convection relative to heat conduction in the fluid flow. On the other hand, the main factor of the heat transfer in this flow is the convective heat transfer between the walls and the fluid flow. Thus, it is expected that an increase in Pr number increases the magnitude of the temperature distributions.

Nusselt number (Nu)
Nusselt number obtained using Eq. (39) for different conditions are presented in this section. Similar to the discussion presented for the temperature and velocity fields, since the initial transient part of the Nu number decays shortly, here we only present the results of the fully developed flow. Figure 12    Present work Steady flow [16] Nu Present work Steady flow [16] Nu Present work Steady flow [16] Nu Present work Steady flow [16] Nu Present work Steady flow [16] Nu Present work Steady flow [16]  the steady-state Darcian flow (i.e.,Da → 0 ) only S * s has effect on the Nusselt number and S * f does not affect it [16]. Furthermore, it can be seen that the amplitude of the Nu number increases with the increase of S * f , since the value of S * increases. Figure 13 shows the effect of Darcy number ( Da ) on Nu number for S * s = 0 and S * f = 0 . It is seen that increasing the Da number decreases the mean value of Nu number [16], while increases significantly the amplitude of the oscillation. For example, as Da number increases from 10 −5 to 10 −2 , the mean value of Nu number decreases by about 10%, while the amplitude of oscillation increases by more than 15 times. An increase in Da number increases the amplitude of the unsteady velocity according to Fig. 4 and the amplitude of the f,un according to Fig. 10. Hence, from Eq. (42) this will increase the amplitude of oscillation of f,un,b , which results in amplifying the amplitude of oscillation of Nu number according to Eq. (39). For the case of Da = 10 −5 , the amplitude of U un and the amplitude of f,un are so small that hence the oscillation of Nu number is very small.
The effects of thermal conductivity ratio k and Biot number, on the Nu number, are shown in Figs. 14 and 15 for Da = 10 −5 and Da = 10 −2 , respectively. Results for each case are compared with Nu numbers for the steady flow with the same conditions. Similar to Fig. 13, it is seen that for low Da number of 10 −5 , the ratio of amplitude-to-mean value of the Nusselt number for all cases in Fig. 14 is less than 1%, and the pulsations are almost negligible. This figure indicates that variations in the values of k or Bi do not affect significantly the pulsation of Nu number. From Fig. 15, it is seen that for a high Da number, the amplitude of Nu number reduces with the increase in k or Bi.  [16] Steady flow [16] 124  The effect of the coefficient of pressure amplitude ( in Eq. 2) is shown in Fig. 16. It is seen that an increase in , increases significantly the amplitude of Nu number, which was also presented in previous works containing pulsatile flow in an empty channel or pipe [24,25,62]. This demonstrates that waveform amplitude has a deterministic role in controlling the rate of convective heat transfer between the walls and the fluid flowing in the channel. Figure 17 shows the impact of the dimensionless frequency ( ) on Nu number for = 0.1 and for other conditions similar to Fig. 16. It is seen that as the value of increases, the amplitude of pulsation of Nu number decreases. This is in agreement with the findings of previous works (e.g., [24,25,58]). The main influence of is on the period of oscillation, as expected.
For the steady-state fully developed flow in a channel filled with a porous material, Nu number is independent of the Prandtl (Pr) number [13,16]. While for the pulsating flow, Pr number has significant influence on Nu number [25,41]. Figure 18 shows the effect of Pr number on the Nu number under LTNE condition for Da = 10 −5 (Fig. 18a) and Da = 10 −2 (Fig. 18b). It is seen that for small Da number, the amplitude of Nu number increases slightly with the increase in Pr number. For large Da numbers, as Da increases the amplitude of Nu oscillation increases. Further increase in Da number decreases the amplitude of Nu number. Additionally, it is seen that the effect of Pr number on the Nu number is profound for high Da numbers. Comparing Fig. 18a   Present work Steady flow [16] Present work Steady flow [16] Present work Steady flow [16] Present work Steady flow [16]  results of Kuznetsov and Nield [41] obtained incorporating the LTE model. Kuznetsov and Nield [41] investigated the effect of the Pr number in a channel filled with a saturated porous medium under LTE model. They found that for large Da numbers, increasing Pr number decreases the amplitude of Nu number. While for small Da numbers, the effect of Pr number on Nu number is found to be negligible.

Conclusions
This paper studied analytically the problem of forced convection heat transfer of pulsating flow due to oscillatory applied pressure gradient in a channel filled with a porous medium subjected to a constant wall heat flux. By considering internal heat generations in the solid and the fluid phases in the porous region, energy equations were solved using a local thermal non-equilibrium (LTNE) model. Considering specific conditions at the wall interface, the approach known  as Model A for the thermal boundary conditions was used to solve the governing equations. Exact solutions for the unsteady velocity ( U un ), temperature of the solid phase ( s,un ) and temperature of the fluid phase ( f,un ) and Nusselt number (Nu) were obtained. The effect of different parameters is analyzed. These parameters are Darcy number (Da), Prandtl number (Pr), Biot number (Bi), fluid-to-solid thermal conductivity ratio ( k ), heat generation in the solid phase ( S * s ) and fluid phase ( S * f ), dimensionless frequency ( ) and the coefficient of pressure amplitude ( ). Important results of this study are summarized as follows: • The amplitude of the unsteady velocity increases with the increase of or Da , while decreases with the increase in dimensionless frequency . • The amplitude of the unsteady dimensionless temperatures for the fluid phase is relatively higher than that of the solid phase. •  All that remains is to determine a n (t * ) in Eq. (50) which solves the non-homogeneous partial differential equations of momentum Eq. (24). To find a n (t * ), the term-by-term differentiations are calculated as: Equations (53) and (54) are now substituted into Eq. (24), and the resultant equation can be written as: Due to orthogonality of the eigenfunctions ( � n (Y) ) obtaining from the Sturm-Liouville eigenvalue problem of Eq. (34), the following equation for the time-dependent coefficient of the term cos (2n − 1) 2 Y in Eq. (55) is obtained [64]: Equation (56) is a first-order ODE problem with respect to t * which can easily be solved to find a n (t * ) . Using initial condition presented in Eq. (52), a n (t * ) is obtained (Eq. 30).

Energy equations
Similar to the procedure deployed to solve for the velocity, the solutions for the temperatures of the two phases are divided into steady and unsteady components [25] as: (52) a n (0) = 0.
Substituting Eqs. (57) and (58) into Eqs. (15) and (16), the resultant equations are as: Since T * f,st t * = 0 and according to Eq. (57) T * f,un is not a function of X [25], each of the energy Eqs. (59) and (60) can be divided into two separate steady and unsteady equations and solved separately as explained in the following.

Steady energy equations
Using Eqs. (59) and (60), a new set of steady equations are derived as: The flow is considered to be fully developed, implying that T * f,st X = C is constant [66]. Using a dimensionless variable st = k s,eff (T−Tw,st) The boundary conditions for the steady equations are also derived utilizing a same procedure as for momentum equation by separating steady and unsteady components from Eqs. (17) and (18). The heat flux, q w , is a constant value, and therefore, it is applied to boundary condition of the steady components [24]; so employing the dimensionless variable , boundary conditions can be rewritten as: Furthermore, one more boundary condition is needed to solve steady set of energy equations that is derived from Eq. (10).
Equation (66)  Finally, the unsteady temperatures of the two phases can be also normalized in the form ( un = T * un − T * w,un ) using the unsteady wall temperature (see Eqs. 34,35).