0 Introduction

Flow and heat transport through porous media is widespread in natural system and engineering practice^[1-3]. Two types of models for the energy equation were commonly described in various heat transfer process, involving local thermal equilibrium (LTE) or non-equilibrium (LTNE) model^[4-6]. Using LTE and LTNE models respectively, Lee and Vafai^[7] obtained and compared the exact temperature solutions of thermal transfer in a porous media channel. The thermal bifurcation phenomenon in a forced convection porous medium channel by employing LTNE model was analyzed by Yang and Vafai^[8]. Considering the Darcy-Binkman-Forchheimer flow, Jiang and Ren^[9] used LTNE model to discuss forced convection process in a porous media channel. A metal foam microchannel heat sink was constructed by Li et al.^[10], and the heat transfer characteristics were numerically investigated using Brinkman extended Darcy and LTNE model. Based on fully developed flow and LTNE model, Xu^[11]theoretically investigated the flow inertia, velocity jumps, and thermal slips in a microchannel partially filled with micro-porous media. Using LTNE model, Shen et al.^[12]numerically studied the heat exchange characteristics of copper metal foam with different shapes in micro channels. Based on the Brinkman-Darcy momentum equation, considering viscous dissipation effect, Chen and Tso^[13] numerically simulated the forced convection process using LTNE model in a porous media channel, and gained the fluid-solid temperature distributions with invariant heat flux.

Most initial research is based on Newtonian fluid. In numerous practical engineering applications, many fluids belong to the category of non-Newtonian fluid, such as blood and cell fluid in human body, slurry, pigment, paint, and industrial grease in chemical production, milk, and beverage in food industry, showing rheological properties of non-Newtonian fluid^[14-15]. In the study of Thayalan and Huang^[16], the momentum boundary layer integral method was used to solve the velocity distribution and the temperature distribution using LTNE model. The effects of Da, n, Bi, and k on convection heat transfer were analyzed. However, the inertial effect and viscous dissipation effect of non-Darcy flow were not considered. Using Brinkman-Forchheimer extended Darcy flow model, Tian et al.^[17]compared and analyzed three different viscous dissipation terms on forced convection thermal characteristic of power-law fluid in a porous media channel, but did not take into consideration fluid-solid temperature difference.

To strengthen heat transfer of tube-fin heat exchangers^[18] or heat dissipation of microelectronic equipment, microchannel cooling technology using nanofluid is an effective solution. For example, for the convection heat transfer process of magnetic nanofluid^[19-20] in porous microchannels, the non-Newtonian characteristics of nanofluids^[21-22], LTNE effects, and other micro scale effects^[23] (such as viscous dissipation effect of nanofluids, nonlinear drag inertia effect, thermal dispersion effect, etc.) may need to be considered concurrently.

In summary, for the rheological characteristics of power-law fluid, there are few reports on the thermal transport in porous media including fluid-solid temperature difference, Darcy-Brinkman-Forchheimer flow, and viscous dissipation effect at the same time. In this paper, a convection model of power-law fluid in porous media channel is established. The effects of dimensionless parameters, such as Bi, k, Br, Da, F, and n on temperature distribution and heat transfer characteristics are analyzed numerically.

1 Physical Mathematical Model

In this study, as shown in Fig. 1, an infinite parallel channel of height 2H was investigated. Considering the fully developed steady state of two-dimensional laminar flow heat transfer, the z-direction flow heat transfer is not considered (the width of the model is much larger than the height, and the conditions in the width direction are uniform and symmetrical). Assume the fluid is a single-phase, incompressible power-law fluid, the porous skeletons are rigid and isotropic, and the fluid and skeletons are in local thermal non-equilibrium. Except for the fluid viscosity, the themo-physical properties of the fluid and solid phases are considered constant. The velocity gradient along flow direction (in x-direction)is infinitesimal, and the chemical reaction, radiation heat, and the axial conduction are ignored.

The Forchheimer-Brinkman extended Darcy's law by modifying the momentum equation for power-law fluid^[24] is expressed as

$ \frac{\mu}{\varepsilon^{n}} \frac{\mathrm{d}}{\mathrm{d} y}\left[\left|\frac{\mathrm{d} u}{\mathrm{~d} y}\right|^{n-1} \frac{\mathrm{d} u}{\mathrm{~d} y}\right]-\frac{\mu}{K^{*}} u^{n}-\rho b u^{2}+\frac{\mathrm{d} p}{\mathrm{~d} x}=0 $

(1)

Dimensionless parameters are defined as follows: n represents the power-law index, transverse coordinate Y =y/H; centerline velocity u_c=(K^*dp/(μdx))^1/n; characteristic velocity based on centerline velocity U=u/u_c; inertial parameter F= ρbK^*(u_c)^2－n/μ; Da number Da=K^*2/(n+1)/H²/ε^2n/(n+1). Then Eq. (1) with its boundary conditions is non-dimensionalized as

$ \frac{\mathrm{d}}{\mathrm{d} Y}\left[\left|\frac{\mathrm{d} U}{\mathrm{~d} Y}\right|^{n-1} \frac{\mathrm{d} U}{\mathrm{~d} Y}\right]=\frac{U^{n}+F U^{2}-1}{D a^{(1+n) / 2}} $

(2)

The corresponding energy equations are written as

$ k_{\mathrm{f}, \text { eff }} \frac{\partial^{2} T_{\mathrm{f}}}{\partial y^{2}}+h_{\mathrm{sf}} a_{\mathrm{sf}}\left(T_{\mathrm{s}}-T_{\mathrm{f}}\right)+\varOmega=\rho c_{P} u \frac{\partial T_{\mathrm{f}}}{\partial x} $

(3)

$ \varOmega=\frac{\mu u^{n+1}}{K^{*}}+\frac{\mu}{\varepsilon^{n}}\left|\frac{\mathrm{d} u}{\mathrm{~d} y}\right|^{n-1}\left(\frac{\mathrm{d} u}{\mathrm{~d} y}\right)^{2}+\rho b u^{3} $

(4)

$ k_{\mathrm{s}, \text { eff }} \frac{\partial^{2} T_{\mathrm{s}}}{\partial y^{2}}-h_{\mathrm{sf}} a_{\mathrm{sf}}\left(T_{\mathrm{s}}-T_{\mathrm{f}}\right)=0 $

(5)

where T_f, T_s, k_{s, eff}, k_{f, eff} are fluid-solid temperature and thermal conductivity, respectively; Ω is representative of different viscous dissipation, involving Darcy term (internal heating), Al-Hadhrami term (frictional heating), and Forchheimer term (nonlinear drag power); h_sf and a_sf are interface heat transfer coefficient and specific surface area, respectively.

According to the given conditions, the temperature gradient $\partial T/\partial x = \partial {T_{\rm{w}}}/\partial x = \partial {T_{\rm{m}}}/\partial x$. Define the wall temperature as T_w, the average temperature of fluid ${T_{\rm{m}}} = \int_0^H {uT{\rm{d}}y/H{u_{\rm{m}}}} $, and the average velocity${u_{\rm{m}}} = \int_0^H {u{\rm{d}}y/H} $. Then the following equation is obtained:

$ \frac{\mathrm{d} T_{\mathrm{m}}}{\mathrm{d} x}=\frac{q_{\mathrm{w}}+\int_{0}^{H} \varOmega}{\rho c_{p} H u_{\mathrm{m}}} $

(6)

$ q_{\mathrm{w}}=\left.k_{\mathrm{f}, \text { eff }} \frac{\partial T_{\mathrm{f}}}{\partial y}\right|_{y=0}+\left.k_{\mathrm{s}, \text { eff }} \frac{\partial T_{\mathrm{s}}}{\partial y}\right|_{y=0} $

(7)

Other dimensionless parameters are defined as follows: $\theta = {k_{{\rm{s}}, {\rm{ eff }}}}\left( {T - {T_{\rm{w}}}} \right)/\left( {{q_{\rm{w}}}H} \right)$; velocity based on the average velocity U'=u/u_m; modified inertial parameter ${F^\prime } = \frac{{\rho b{K^*}{{\left( {{u_{\rm{m}}}} \right)}^{2 - n}}}}{\mu } = F{\left( {\frac{{{u_{\rm{m}}}}}{{{u_{\rm{c}}}}}} \right)^{2 - n}}$; inter-phase heat transfer coefficient $Bi = \frac{{{h_{{\rm{sf}}}}{a_{{\rm{sf}}}}{H^2}}}{{{k_{{\rm{s}}, {\rm{ eff }}}}}}$; coefficient ratio $k = \frac{{{k_{{\rm{f}}, {\rm{ eff }}}}}}{{{k_{{\rm{s}}, {\rm{ eff }}}}}} = \frac{{\varepsilon {k_{\rm{f}}}}}{{(1 - \varepsilon ){k_{\rm{s}}}}}$; revised Brinkman number $B{r^\prime } = BrD{a^{(n + 1)/2}} = \mu {u_{\rm{m}}}^{n + 1}/\left( {{\varepsilon ^n}{q_{\rm{w}}}{H^n}} \right)$.

Substituting Eqs. (4), (6), and (7) into Eq. (3), combining with Eq. (5), and integrating along the channel cross section, the following energy equation is obtained:

$ \begin{aligned} k &\frac{\partial^{2} \theta_{\mathrm{f}}}{\partial Y^{2}}+\frac{\partial^{2} \theta_{s}}{\partial Y^{2}}=U^{\prime}\left(1+B r \int_{0}^{1}\left(U^{\prime}\right)^{n+1} \mathrm{~d} Y+\right. \\ &\left.B r^{\prime} \int_{0}^{1}\left|\frac{\mathrm{d} U^{\prime}}{\mathrm{d} Y}\right|^{n-1}\left(\frac{\mathrm{d} U^{\prime}}{\mathrm{d} Y}\right) \mathrm{d} Y+B r F^{\prime} \int_{0}^{1}\left(U^{\prime}\right)^{3} \mathrm{~d} Y\right)- \\ &B r\left(U^{\prime}\right)^{n+1}-B r^{\prime}\left|\frac{\mathrm{d} U^{\prime}}{\mathrm{d} Y}\right|^{n-1}\left(\frac{\mathrm{d} U^{\prime}}{\mathrm{d} Y^{\prime}}\right)^{2}-B r F^{\prime}\left(U^{\prime}\right)^{3} \end{aligned} $

(8)

When the channel wall has finite thickness and the thermal conductivity is relatively high, under the fully developed steady state of two-dimensional flow heat transfer, θ_f=θ_s=θ_w, and considering axisymmetric boundary conditions, the corresponding dimensionless boundary conditions can be justified by

$ Y=0, U=U^{\prime}=0, \theta_{\mathrm{f}}=\theta_{\mathrm{s}}=0 $

(9a)

$ Y=1, \mathrm{~d} U / \mathrm{d} Y=\mathrm{d} U^{\prime} / \mathrm{d} Y=0, \partial \theta_{\mathrm{f}} / \mathrm{d} Y=\partial \theta_{\mathrm{s}} / \mathrm{d} Y=0 $

(9b)

It is noteworthy that Br characterizes the intensity of the viscous dissipation (Br=0, which indicates no viscous dissipation). It can be seen from the above expression that the viscous dissipation effect is closely related to Br, Da, F, and n. The Nusselt number^[25] can be evaluated as

$ N u=\frac{2 H q_{\mathrm{w}}}{k_{\mathrm{f}, \mathrm{eff}}\left(T_{\mathrm{w}}-T_{\mathrm{m}}\right)}=-\frac{2}{k \theta_{\mathrm{fb}}} $

(10)

where dimensionless average temperature of fluid phase $\theta_{\mathrm{fb}}=\frac{\int_{0}^{1} U \theta_{\mathrm{f}} \mathrm{d} Y}{\int_{0}^{1} U \mathrm{~d} Y}=\frac{\int_{0}^{1} U \theta_{\mathrm{f}} \mathrm{d} Y}{U_{\mathrm{m}}}$.

For LTE model, θ_f=θ_s. According to Eqs. (8) and (9), the dimensionless Nusselt number can be rewritten as follows:

$ N u_{\mathrm{LTE}}=\frac{2 H q_{\mathrm{w}}}{k_{\mathrm{f}, \mathrm{eff}}^{\prime}\left(T_{\mathrm{w}}-T_{\mathrm{m}}\right)}=-\frac{2}{k \theta_{\mathrm{b}}}=-\frac{2(k+1)}{k \theta_{\mathrm{fb}}} $

(11)

The difference between LTE and LTNE models can be evaluated as E=Nu_LTE-Nu

2 Numerical Methods and Solvers

From the expressions of momentum equation and energy equation, it can be deduced that when the inertia flow and viscous dissipation effect are considered, there is a strong nonlinear effect of non-Darcy flow in porous media, which must be solved by numerical method. To do this, Eqs. (2) and (8) are first recombined and reduced to first order differential Eqs. (12). Obviously, this is a boundary value problem of ordinary differential equations. Then, assuming the initial value and a reasonable step size (0.01), the classical Runge Kutta method is used to carry out variable alternating numerical iteration from 0 to 1 until the iteration accuracy (10^-5) and the boundary condition Eqs. (13) are satisfied. The parameters defined are as follows:

$ y_{1}=U(Y), y_{2}=(\mathrm{d} U / \mathrm{d} Y)^{n} $

$ y_{3}=\int\limits_{0}^{Y} U \mathrm{d} Y, y_{4}=\int\limits_{0}^{Y} U^{n+1} \mathrm{d} Y $

$ y_{5}=\int\limits_{0}^{Y} U^{3} \mathrm{d} Y, y_{6}=\int\limits_{0}^{Y}(\mathrm{d} U / \mathrm{d} Y)^{n+1} \mathrm{d} Y $

$ y_{7}=U_{\mathrm{m}}, y_{8}=\left(U_{\mathrm{m}}\right)^{n+1} $

$ y_{9}=\left(U_{\mathrm{m}}\right)^{3}, y_{10}=\int\limits_{0}^{1} U^{n+1} \mathrm{d} Y $

$ y_{11}=\int\limits_{0}^{1} U^{3} \mathrm{d} Y, y_{12}=\int\limits_{0}^{1}(\mathrm{d} U / \mathrm{d} Y)^{n+1} \mathrm{d} Y $

$ y_{13(\mathrm{f})}=\theta_{(\mathrm{f})}(Y), y_{13(\mathrm{s})}=\theta_{(\mathrm{s})}(Y) $

$ y_{14(\mathrm{f})}=\mathrm{d} \theta_{(\mathrm{f})} / \mathrm{d} Y, y_{14(\mathrm{s})}=\mathrm{d} \theta_{(\mathrm{s})} / \mathrm{d} Y $

and

$ \begin{array}{l} \begin{array}{l} y_{1}^{\prime}=\left(y_{2}\right)^{\frac{1}{n}}, y_{2}^{\prime}=\left(\left(y_{1}\right)^{n}+F\left(y_{1}\right)^{2}-1\right) / D a^{\frac{n+1}{2}}\\ y_{3}^{\prime}=y_{1}, y_{4}^{\prime}=\left(y_{1}\right)^{n+1}, y_{5}^{\prime}=\left(y_{1}\right)^{3}, y_{6}^{\prime}=\left(y_{2}\right)^{\frac{n+1}{n}}\\ y_{7}^{\prime}=y_{8}^{\prime}=y_{9}^{\prime}=y_{10}^{\prime}=y_{11}^{\prime}=y_{12}^{\prime}=0, y_{13}^{\prime}=y_{14} \end{array}\\ \ \ \ \begin{aligned} k y&_{14(\mathrm{f})}^{\prime}+B i\left(y_{13(\mathrm{~s})}-y_{13(\mathrm{f})}\right)= \\ &\frac{y_{1}}{y_{7}}\left(1+B r \frac{y_{10}}{y_{8}}+B r^{\prime} \frac{y_{12}}{y_{8}}+F^{\prime} B r \frac{y_{11}}{y_{9}}\right)- \\ &B r \frac{\left(y_{1}\right)^{n+1}}{y_{8}}-B r^{\prime} \frac{\left(y_{2}\right)^{(n+1) / n}}{y_{8}}-F^{\prime} B r \frac{\left(y_{1}\right)^{3}}{y_{9}} \cdot \\ &y_{14(\mathrm{~s})}^{\prime}=B i\left(y_{13(\mathrm{~s})}-y_{13(\mathrm{f})}\right) \end{aligned} \end{array} $

(12)

$ \begin{array}{l} y_{1}(0)=y_{2}(1)=y_{3}(0)=y_{4}(0)=y_{5}(0)= \\ \ \ \ \ y_{6}(0)=y_{13}(0)=y_{14}(1)=0 \\ \ \ \ \ \ \ y_{7}(0)=y_{3}(1), y_{8}(0)=\left(y_{3}(1)\right)^{n+1} \\ \ \ \ \ \ \ \ \ y_{9}(0)=\left(y_{3}(1)\right)^{3}, y_{10}(0)=y_{4}(1) \\ \ \ \ \ \ \ \ \ \ \ y_{11}(0)=y_{5}(1), y_{12}(0)=y_{6}(1) \end{array} $

(13)

Then the expression of the Nu number can be obtained by the Eqs. (10), (11), and (12):

$ N u=-\frac{2 y_{7}}{k \int_{0}^{1} y_{1} y_{13(\mathrm{f})} \mathrm{d} Y} $

(14)

$ N u_{\mathrm{LTE}}=-\frac{2(k+1) y_{7}}{k \int_{0}^{1} y_{1} y_{13(\mathrm{f})} \mathrm{d} Y} $

(15)

To prove the validation of the numerical program, the present solution with indirect decoupling analytical solution by Wang et al.^[1] is compared in Fig. 2, which illustrates dimensionless fluid-solid temperature profile with varying Biot number and k number along the cross section when Da=0.1, Br=0, and F=0. It can be seen that the present numerical predictions and the work by Wang et al.^[1] are in good agreement.

In order to further verify numerical solution, suppose Br=F=0, θ_(f)=θ_(s)(assuming LTE model). When n=1 (Newton fluid), Da≤0.001, the value of Nu is 5.92, which is in line with the conclusion that when Darcy is a small number (i.e., piston flow, Nu is close to 6); when Da≥1000, the simulation result, Nu=4.121, is also consistent with the conclusion that when Darcy is very large (i.e., parallel plate Poiseuille flow, Nu tends to be 70/17), which is consistent with the research of Nield et al.^[26] When n≠1 (non-Newtonian fluid), Da≥100, the simulation result of non-Newton fluid (n = 0.5, and n=1.5) are 4.380 and 4.010, respectively, which are in line with the conclusion of Chen and Tso^[27].

3 Results and Analysis

To characterize the effects of dimensionless parameters on convection heat transfer characteristics, such as Bi, k, Br, Da, F, n, and so on, there are many variables of velocity field and temperature field that can be solved. The velocity field is mainly affected by dimensionless Darcy number and Forchheimer number, and the dimensionless velocity profile of fluids for different power-law indexes with varying Da and F number had been characterized and analyzed in Ref. [17]. The relationship of temperature distribution and the Nusselt number of power-law fluid of shear thinning are focused on here.

3.1 Impact of Bi

Fig. 3 describes the temperature profile of fluid and solid with change of Biot number and power-law indexes along the cross section in porous media. As depicted in Fig. 3, with the increase of Bi, the temperature difference of fluid and solid decreases gradually. When Bi =100, the change of temperature difference is not obvious. It can be inferred that with a lager Bi (>1000), the temperature of the two phases tends to be the same. This shows that a larger Bi reduces fluid-solid thermal resistance and enhances convection heat transfer rate. Meanwhile, it is found that the direction of heat flux near wall is opposite, and the phenomenon of heat flux divergence^[28]occurs. With the decrease of Bi, the phenomenon of heat flux divergence is more obvious. It can also be seen that with the decrease of n, the fluid-solid temperature difference increases gradually, indicating the fluid-solid heat exchange is gradually weakened. It can be inferred that with the decrease of n, the rheological property of pseudoplasticity is enhanced, the ability of intermolecular entanglement is strong, and the mesh structure formed between the molecules has larger resistance to the flow, which affects the velocity distribution. Meanwhile, the skeleton shear makes the viscosity dissipation between molecules further enhanced, reducing local convection heat transfer rate.

3.2 Impact of k

The temperature profile with variable of k number and power-law indexes in porous media along the cross section is described in Fig. 4. With the increase of k, both fluid and solid temperatures increase, but fluid temperature rises faster than that of solid, which makes the fluid-solid temperature difference slightly reduce. As shown in Fig. 4, it is found that with the decrease of k (especially when k < 0.1), similar to Bi, bifurcation phenomenon appears at the temperature gradient near the wall, which indicates that the heat flow divergence will occur under certain parameter combination. It can also be seen that with the decrease of n, the fluid-solid temperature difference increases gradually, indicating that convection heat transfer rate decreases.

3.3 Impact of Br

Fig. 5 describes the temperature profile with the change of Br and different power-law indexes in porous media along the cross section. It is shown that without the viscous dissipation terms (Br = 0), the temperature distribution of the two phases corresponding to three different fluids is almost equal, and the fluid-solid temperature difference is also small. When considering viscous dissipation effect (Br≠0), the fluid-solid temperature difference changes obviously with the increase of Br, which can be explained that when numbers of Br, Da, and F are not zero, the viscous dissipation effect generated by Darcy term, Alhdharami term, and Forchheimer term will overlap with each other. With the increase of Br, the friction heat dissipation effect corresponding to Alhdharami term is more intense in the vicinity of the wall and decreases towards the channel center. The fluid temperature increases by heating effects^[29], and the temperature gradient bifurcation phenomenon also appears near the wall. On the contrary, the internal heat dissipation generated by Darcy term plays a major role near the channel center, which reduces fluid temperature and increases solid temperature. It should be noted that there is a critical value of Br, where the sign of the dimensionless temperature changes. It can also be seen that with the decrease of n, the rheological property of pseudoplasticity is enhanced, and the temperature difference increases gradually, which shows that the convection heat transfer between two phases decreases somewhat, and the effects of Br on the pseudoplastic fluid is more obvious.

3.4 Impact of Da

Fig. 6 plots the temperature profile with change of Darcy number and different power-law indexes along the cross section. It is observed that the fluid-solid temperature difference increases obviously with the increase of Da. It can be explained that the larger Da makes the flow deviate from Darcy flow, and the non-Darcy flow area increases gradually. The velocity profile tends to be a parabola shape, the contact surface area of heat transfer decreases for the increased porosity, leading to the decrease of convection heat transfer rate, and the fluid-solid temperature difference gradually increases. Meanwhile, it is presented that when Da is small, the temperature difference increases with decrease of n, indicating the convection heat transfer is weakened. On the contrary, with a large Da, a smaller fluid-solid temperature difference is obtained with the decrease of n, and the convection heat transfer is enhanced. The above results show that the range of influence of Da on pseudoplastic fluid is small, which is consistent with Ref. [17].

3.5 Impact of F

The temperature profile along the cross section with the change of F number and different power-law indexes is depicted in Fig. 7. It is observed that increasing F induces a slight increase of fluid and solid temperatures, and the phenomenon of heat flux divergence occurs near the wall. The flow speed increases with a larger F, and the flow resistance of solid matrix to the fluid increases. However, it can produce a larger velocity gradient near the wall, and the velocity field will be modified by Drag forces dissipation effect. The influence of the change of F on shear thinning fluid is also more obvious.

3.6 Impact of Parameters on Nu

The convection heat transfer characteristics with the change of related parameters (Bi, k, Br, Da, F, and n) is described in Fig. 8. The difference between LTE model and LTNE model is also evaluated. As described in Fig. 8(a), with the increase of Bi and n, Nu increases gradually, while the difference of E decreases gradually with the increase of Bi and the decrease of n. Contrary to Fig. 8(a), with the increase of k, Nu decreases quickly as depicted in Fig. 8(b), while the difference of E decreases with the increase of k and the decrease of n. Fig. 8(a) and Fig. 8(b) show that the local thermal equilibrium is approaching with a larger number of Bi and k.

As illustrated in Fig. 8(c), because of the reverse of fluid temperature, Nu has negative value, which shows that for a certain combination of other parameters, there is a maximum value of Nu when an appropriate Br number is selected. The general trend of the difference decreases with the increase of Br. Fig. 8(d) shows that with the increase of Da, the number of Nu and E decrease gradually, but the decrease is becoming slower for pseudoplastic fluid. As depicted in Fig. 8(e), with the increase of F and n, the number of Nu and E increases gradually. It can be explained that the larger velocity makes the fluid flow deviate from the Darcy region, and LTNE effect is more significant.

4 Conclusions

1) The general model of power-law fluid forced convection considering LTNE and viscous dissipation effects in a porous media channel was established. The calculation expressions of dimensionless fluid-solid temperature profiles and Nusselt number were derived with a clear relationship between physical parameters. The classical fourth-order Runge Kutta method was applied to numerically solve self-programmed solution under a constant heat flux boundary condition.

2) With a bigger Bi and k, the temperature profile of solid was close to that of fluid, which means a local thermal equilibrium was reached. While with the decrease of Bi, k, Br, Da, and the increase of F, the difference of temperature profile between fluid and solid increased gradually, intensifying the thermal non-equilibrium effect, and the fluid with different power-law indexes also had a non-negligible impact on the convection heat transfer characteristics.