Parallel Plate Flow of a Third-Grade Fluid and a Newtonian Fluid With Variable Viscosity

2.1 Perturbation-Iteration Algorithm PIA(1,1) for Second-Order Coupled Differential Equations

The problems investigated are second-order nonlinear equations with coupling between the velocity and temperature. Hence, the theory is presented for such equations. The simplest iteration algorithm is the PIA(1,1) algorithm which is constructed by taking one correction term in the perturbation expansion and only the first-order derivate terms in the Taylor series expansion. For the general nonlinear second-order ordinary differential equations with two dependent variables

where both u and θ are functions of y, u=u(y), θ=θ(y), y being the coordinate vertical to the plates and ε is the perturbation parameter. If the equations do not possess a perturbation parameter, it can be inserted artificially. Broadly speaking, the perturbation parameter is usually inserted as a coefficient in front of the nonlinear terms. Taking just one correction term in the perturbation expansion

Substituting (3) and (4) into (1) and (2), and expanding them in a Taylor series with first-order derivatives leads to

where the total derivative is defined as

and evaluated at ε=0. Note that in this method, the dependent variables and their derivatives are considered to be independent of each other. Their contributions to the total derivative are calculated separately and added to each other as shown in (6). Using (5), one will have two coupled equations of (u_c)_n and (θ_c)_n to be solved in order. With two initial guesses, u₀ and θ₀, which satisfy the boundary conditions, (u_c)₀ and (θ_c)₀ can be calculated and later substituted into (3) and (4) to calculate u₁ and θ₁. This iteration procedure can be repeated until satisfactory results are obtained.

2.2 Perturbation-Iteration Algorithm PIA(1,2) for Second-Order Coupled Differential Equations

In this iteration algorithm, only one correction term in the perturbation expansion is taken and up to second-order derivatives in the Taylor series expansion are included, i.e. n=1, m=2. The algorithm is named as PIA(1,2). Again, the perturbation expansions with one correction terms

are substituted into (1) and (2) and expanded in a Taylor series up to second-order derivatives leading to

One will have two differential equations in terms of the unknowns (u_c)_n and (θ_c)_n again, which have to be solved and substituted into (7) and (8). The iterative procedure is set up in the same way as explained for the PIA(1,1).

3.1 Perturbation-Iteration Algorithm PIA(1,1)

If the reference viscosity is selected as the Newtonian viscosity, μ=1 for the constant viscosity case. Applying the constant viscosity, (10) and (11) are rewritten in the following form:

where ε is an artificially introduced small parameter. Substituting the perturbation expansions (3) and (4) into (15) and (16), then expanding in Taylor series with first-order derivatives (5) and rearranging yields

Starting with the initial functions that satisfy the boundary conditions

and using (17), (18), (3) and (4), one will obtain the first iteration solutions as follows:

Proceeding in the same way, the subsequent iterations are

Note that the solutions are of polynomial type for this specific algorithm and becomes lengthy as the number of iterations increase.

3.2 Perturbation-Iteration Algorithm PIA(1,2)

For this case, the problem will be solved using the perturbation-iteration algorithm PIA(1,2). Substituting the perturbation expansions (7) and (8) into (15) and (16), then expanding in Taylor series up to second-order derivatives (9), the iteration formulas for the PIA(1,2) reduces to

With the initial choices

the first and the second iteration solutions are

For the sake of brevity, the lengthy solution calculated for θ₂ using a symbolic manipulation program is not written here. The PIA(1,2) generated more complex solutions in functional forms as compared to the polynomial form solutions of the PIA(1,1).

3.3 Perturbation Solution

This particular problem was solved using regular perturbation methods with one correction term by Yürüsoy, Pakdemirli, and Yilbas [22]. The solutions they found under constant viscosity assumption are

It is well-known from the perturbation theory that the correction terms must be much smaller than the leading terms for the solutions to be valid. Hence, the validity criterion of these two solutions is

It should be noted that this solution is equivalent to the second iteration of PIA(1,1) when ε=1.

3.4 Comparisons With the Numerical Solution

In this section, for different values of the pressure gradient C, the non-Newtonian coefficient Λ and the Brinkman number Γ, solutions of the PIA(1,1), PIA(1,2), and the regular perturbation will be compared with each other as well as with the numerical solutions. The numerical solutions are calculated with Wolfram Mathematica’s built-in command, NDSolve.

Figure 1 shows that all the solutions match well enough, when small values of the parameters satisfying the validity criterion are chosen, as expected. Especially, the non-Newtonian behaviour is assumed to be small in the calculations of Figure 1.

Figure 1:

Comparison of the approximate velocity solutions with the numerical solution. ε=1, C=–2, Λ=0.025, Γ=5, Cr=0.2.

When the values of the parameters are increased further, the validity criterion does not hold anymore. The perturbation solution predicts backflow adjacent to the plates which are unphysical for a constant pressure drop between the plates (Fig. 2). The PIA(1,2) solution predicts the velocity profile better than the PIA(1,1) solution.

Figure 2:

Comparison of the approximate velocity solutions with the numerical solution. ε=1, C=–2, Λ=0.6, Γ=5, Cr=4.8.

When the physical parameters are considerably larger and far beyond the validity criterion (Cr=8), the perturbation solution predicts a backflow within all the domain which is not correct at all (Fig. 3). The PIA(1,1) solution does not predict the velocity profiles quantitatively well with an exaggerated maximum velocity. The PIA(1,2) velocity profiles are in acceptable agreement with the numerical ones.

Figure 3:

Comparison of the approximate velocity solutions with the numerical solution. ε=1, C=–2, Λ=1, Γ=5, Cr=8.

As far as the temperature profiles are considered, the matches are better among the solutions. For the criterion of 0.2, the temperatures are indistinguishable as shown in Figure 4. One reason might be that the first non-linear term involving only the Brinkman number is taken as an O(1) term in the energy equation, i.e. (11).

Figure 4:

Comparison of the approximate temperature solutions with the numerical solution. ε=1, C=–2, Λ=0.025, Γ=5, Cr=0.2.

In Figure 5, the perturbation and the PIA(1,1) solutions start diverging as the criterion gets larger. PIA(1,1) seems to perform slightly better than the perturbation solution. PIA(1,2), on the other hand, agrees well with the numerical solution. Nevertheless, all the solutions are within the acceptable percentage error.

Figure 5:

Comparison of the approximate temperature solutions with the numerical solution. ε=1, C=–2, Λ=0.6, Γ=5, Cr=4.8.

The last case for the temperature solutions are presented in Figure 6. The criterion is 8, a fairly large number compared to 1. The perturbation and the PIA(1,1) solutions diverge from the numerical solution while the PIA(1,2) solution still matches with the numerical solution reasonably. For this specific case, the perturbation solution performed slightly better than the PIA(1,1) solution.

Figure 6:

Comparison of the approximate temperature solutions with the numerical solution. ε=1, C=–2, Λ=1, Γ=5, Cr=8.

In addition to the comparisons of the velocity and temperature profiles, the comparisons of the dimensionless shear stress and the Nusselt number are also given since these dimensionless parameters are technologically important. The dimensionless equation defining the shear stress of a third-grade fluid is

where u is the dimensionless velocity and Λ is the dimensionless parameter related to the non-Newtonian behaviour as defined before. The dimensionless shear stress is related to the dimensional one (with over tilde) as follows

Using (37), Table 1 is computed for the comparison of the shear stresses calculated by the numerical and approximate solutions for different non-Newtonian parameter values.

An interesting feature in the numerical results is that the shear stresses remained almost unchanged by the increase in the non-Newtonian parameter. Increasing Λ decreases the slope of the velocity due to the non-Newtonian effects and their combination in the stress equation produces similar values. However, approximate solutions give diverging results near boundaries as Λ increases. This is due to the higher slopes of the approximate solutions than the numerical ones near the boundaries as it can also be verified visually from the graphs. The shear stress decreases as y is increased till the midline where the flow reaches its maximum speed and the lowest value of shear stress is observed there. When the three approximate solutions are compared, PIA(1,2) complies best with the numerical results.

For the heat transfer, the local Nusselt number is important

where q˜ denotes the heat flux and the rest of the parameters are the same as defined in (12).

The Nusselt number and consequently the heat transfer rate decrease as we get closer to the upper plate (Table 2). PIA(1,2) solutions are closer to the numerical ones than the other solutions. As the non-Newtonian coefficient increases, the Nusselt numbers in the lower mid-part decrease, but a slight increase is observed for the upper mid-part.

In conclusion, for the velocity profiles, all solutions are in reasonable agreement within the range of validity of the perturbation solution. Beyond the range of validity, the PIA solutions agree better with the numerical ones and the PIA(1,2) performs better than the PIA(1,1). The perturbation solution, however, is not reliable as it gives quantitatively and qualitatively wrong behaviour for velocities if the criterion is not satisfied. For the temperature profiles, the qualitative behaviour is the same for all solutions although there are some small variations in the results quantitatively. For temperature profiles, the best agreement is observed for PIA(1,2). In conclusion, for non-Newtonian fluids, the perturbation iteration algorithms may be used to find approximate analytical solutions at least for some simple flow geometries.

4.1 Perturbation-Iteration Algorithm PIA(1,1)

Substituting the Reynolds model viscosity into the problem, (40) and (41) are rewritten in the following form:

One correction term perturbation expansions are substituted into (45) and (46) and then expanded in Taylor series with first-order derivatives. The following equations of iteration are obtained as a result

Simple initial functions that satisfy the boundary conditions can be selected as

The first iteration solutions are then

Finally, the second and third iteration solutions are calculated

Note again that the PIA(1,1) solutions are of polynomial type.

4.2 Perturbation-Iteration Algorithm PIA(1,2)

The iteration equations of PIA(1,2) are found by first substituting the perturbation expansions with one correction terms into (40) and (41) and then expanding in Taylor series up to second-order derivatives

Starting with the same initial guesses as in the previous section

the first iteration solutions are found as

From the definition of the logarithm, εm<1 must hold for the first iteration solutions to be valid. The next iteration is not calculated due to this limited validity of the solution and the complexities encountered in the equations.

4.3 Perturbation Solution

For the regular perturbation solution of the problem, expansions with two correction terms are taken into account

Reynolds model of viscosity is also expanded in a series

Substituting these expansions into (40) and (41) and separating at each order of approximation, one will have three sets of equations to be solved consecutively both for the velocity and temperature profiles. The equations at each order are solved and the solutions are substituted into the perturbation expansions. Thus, the two term perturbation solutions for the problem are

Imposing the condition that each correction term be much smaller than the previous term, the validity criteria are

The first criterion involves the viscosity decaying parameter and the second criterion involves a combination of the pressure drop as well as the Brinkman number.

4.4 Comparisons With the Numerical Solution

In this section, solutions of the PIA(1,1), PIA(1,2), and the regular perturbation will be compared with each other as well as with the numerical solutions for different values of the pressure drop, the viscosity parameter, and the Brinkman number. The numerical solutions are calculated with Wolfram Mathematica’s built-in command, NDSolve.

In Figure 7, the PIA(1,1) and the PIA(1,2) solutions are in good agreement with the numerical solution, while there is a small discrepancy for the perturbation solution especially at the fully developed flow region. Despite the fact that only the first criterion is satisfied, the perturbation solution does not introduce much error.

Figure 7:

Comparison of the approximate velocity solutions with the numerical solution. ε=1, m=0.2, C=–2, Γ=5, Cr₁=0.13, Cr₂=3.33.

In Figure 8, due to the increase in the viscosity parameter, the curves bent to the right and the maximum velocities occur closer to the upper plate where the temperature is greater compared to the lower one. The closest approximate solution to the numerical one is the PIA(1,1) solution. Since m is large and the PIA(1,2) solution consists of a logarithmic function for which εm<1 condition must hold, the solution over-predicts the location and magnitude of the maximum velocity.

Figure 8:

Comparison of the approximate velocity solutions with the numerical solution. ε=1, m=0.8, C=–2, Γ=5, Cr₁=0.53, Cr₂=3.33.

Figure 9 shows the comparison of the approximate solutions for m=1.1. Only the PIA(1,1) and the perturbation solutions are plotted here, as this is an undefined region for the PIA(1,2) solution. The perturbation solution under-predicts the location and magnitude of the maximum velocity. However, PIA(1,1) solution still remains close to the numerical one.

Figure 9:

Comparison of the approximate velocity solutions with the numerical solution. ε=1, m=1.1, C=–2, Γ=5, Cr₁=0.73, Cr₂=3.33.

Regarding the temperature profiles, the PIA(1,1) and the perturbation solution matches well with the numerical solution in the case where one of the criteria is not satisfied (Fig. 10). On the other hand, the first iteration solution of PIA(1,2) being θ₁=y is a rough estimation and differs from the others.

Figure 10:

Comparison of the approximate temperature solutions with the numerical solution. ε=1, m=0.2, C=–2, Γ=5, Cr₁=0.13, Cr₂=3.33.

When the viscosity parameter m is increased, it is observed from Figure 11 that the solutions begin diverging. Here, the perturbation solution performed relatively better than PIA(1,1) especially in the region near the upper plate.

Figure 11:

Comparison of the approximate temperature solutions with the numerical solution. ε=1, m=0.6, C=–2, Γ=5, Cr₁=0.4, Cr₂=3.33.

The last case where m=1.1 is presented in Figure 12. Both criterion are higher than should be for a regular solution, and hence, the perturbation solution is further away from the numerical solution when compared with PIA(1,1) which agrees quite well up to y=0.4 and gradually diverges as y increases.

Figure 12:

Comparison of the approximate temperature solutions with the numerical solution. ε=1, m=1.1, C=–2, Γ=5, Cr₁=0.73, Cr₂=3.33.

Similar to the previous problem, the shear stresses and the Nusselt numbers of the flow along the y-axis will be contrasted. The dimensionless shear stress equation of the steady-state parallel plate flow of a Newtonian fluid is

with the dimensionless parameters defined in (12) and (38). Substituting the Reynolds model of viscosity into (69), the shear stress equation takes the below form

In Table 3, the shear stress values slightly increase as m gets larger. When the approximate results are compared, with few exceptions, the shear stress values calculated by PIA(1,1) comply better to the numerically calculated ones than the other methods.

The Nusselt number for the steady-state parallel plate flow of a Newtonian fluid with temperature-dependent viscosity is defined as

where q˜ is the heat flux. Table 4 is given for the comparison of the Nusselt numbers calculated by the different methods.

In Table 4, it is shown that increasing m has increased the Nusselt number in the region near the lower plate while decreasing it near the upper plate. When the Nusselt number is tracked along the y-axis for a fixed m value, it seems to decrease as we get close to the upper plate. In the table, PIA(1,2) is not presented as its first iteration solution is a rough estimation and the second iteration cannot be calculated due to a singularity in the velocity profile. When the other two approximate solutions are compared, both produce reasonable approximations to the numerical data with few exceptions.

Consequently, for the velocity profiles, all solutions agree well with each other when the validity criteria are satisfied. Exceeding the limit of the criteria, solutions start to diverge with perturbation solution inheriting too much error. In addition to that PIA(1,2) solutions agree well for the small values of the viscosity parameter m; however, they diverge due to the criterion caused by the logarithmic function as m gets larger. In general, the PIA(1,1) gives the most reliable results for the velocity profiles for the case of Newtonian flow. Calculating PIA(1,1) up to third iteration may be the significant cause, while due to limitations, PIA(1,2) algorithm is terminated at the first iteration. For the temperature profiles, PIA(1,1) and the perturbation solutions give similar results. If the criteria are too large, then the perturbation solution ceases to be valid.

The solutions of the previous section and this section can only be contrasted for the limiting cases of Λ=0 when the non-Newtonian parameter vanishes. Both PIA(1,1) and PIA(1,2) solutions given in (25) and (33) reduce to

For (33), when Λ approaches zero, the singularity should be removed by the L’Hopital’s rule. The above solution can also be retrieved for the Newtonian fluid with Reynolds model parameter m=0. For m=0 in (55) and (66), (72) is obtained.