CHAPTER 3

MATHEMATICAL FORMULATION

3.1 Introduction

This chapter presents the complete mathematical formulation of the problem of

steady, laminar, two-dimensional boundary-layer flow and heat transfer of a second-grade

fluid over a permeable stretching sheet embedded in a porous medium. The formulation

proceeds in six stages. First, the physical configuration and the modelling assumptions are

described precisely. Second, the governing system of partial differential equations (PDEs) —

the continuity, momentum, and energy equations — is stated in dimensional form, together

with the constitutive model for the non-uniform internal heat source/sink. Third, the

dimensional boundary conditions are stated. Fourth, a similarity transformation is introduced

that reduces the three PDEs to a coupled system of two ordinary differential equations

(ODEs): a fourth-order nonlinear ODE for the dimensionless stream function f(η) and a

second-order linear ODE for the dimensionless temperature θ(η). Fifth, the dimensionless

parameters governing the problem are identified and their physical significance is discussed.

Sixth, the exact closed-form analytical solution for the velocity field is derived by direct

1
substitution, and the corrected cubic algebraic equation for the decay constant s is obtained

without approximation. The derivation explicitly identifies and resolves the algebraic

discrepancy present in the published literature for the case λ₁ = 0.5.

3.2 Physical Configuration and Modelling Assumptions

3.2.1 Geometry and Flow Configuration

Consider a flat, impermeable elastic sheet that occupies the plane y = 0 and stretches

along the positive x-direction from a fixed origin at x = 0 (the extrusion slot). The sheet

stretches with velocity u_w(x) = cx, where c > 0 is the constant stretching rate (units: s ⁻¹)

and x is the distance from the slot along the sheet. The half-space y > 0 is occupied by a

second-grade fluid that is in steady, laminar, two-dimensional, incompressible flow driven by

the stretching motion of the sheet. The fluid-saturated region y > 0 is modelled as a

homogeneous, isotropic porous medium of constant permeability k'. A uniform suction

velocity v₀ > 0 is applied at the sheet surface (y = 0), drawing fluid toward the sheet. The

wall is maintained at a constant temperature T_w, and the ambient fluid temperature at y →

∞ is T_∞, with T_w > T_∞ (heated wall). A non-uniform internal heat source or sink is

present within the fluid. The physical configuration is illustrated schematically in Figure 3.1.

The coordinate system is Cartesian: x is measured along the sheet in the direction of

stretching, and y is measured perpendicular to the sheet into the fluid domain. The velocity

2
components are u(x, y) in the x-direction and v(x, y) in the y-direction. The temperature field

is T(x, y). The problem is posed on the domain {(x, y) : x > 0, y > 0} with y = ∞ representing

the quiescent far-field region.

3.2.2 Modelling Assumptions

The following assumptions are adopted throughout this study:

(1) Steady state — all flow and thermal quantities are independent of time.

(2) Two-dimensional flow — all quantities are independent of the z-coordinate; the flow

is confined to the x-y plane.

(3) Incompressibility — the fluid density ρ is constant; the divergence of the velocity

field is zero.

(4) Laminar flow — no turbulence modelling is required; the Reynolds number based on

the stretching rate is assumed to be below the transition threshold.

(5) Second-grade constitutive model — the fluid obeys the Rivlin-Ericksen constitutive

equation with material constants μ ≥ 0 and α₁ ≥ 0 satisfying the Dunn-Fosdick

thermodynamic conditions.

3
 (6) Boundary-layer approximation — the streamwise length scale greatly exceeds the

boundary-layer thickness, so that ∂²u/∂x² ≪ ∂²u/∂y², ∂p/∂y ≈ 0, and v ≪ u within the

boundary layer.

(7) Darcy porous medium model — the porous medium resistance is represented by a

Darcy drag term −(ν/k')u, valid for low-Reynolds-number flow through the pores.

(8) Constant thermophysical properties — the kinematic viscosity ν, thermal diffusivity

α, and specific heat c_p are independent of temperature.

(9) Viscous dissipation included — the work done by viscous stresses is converted to

thermal energy and appears as a source term in the energy equation.

(10) Non-uniform heat source/sink — the internal heat generation rate follows the

model of Vajravelu and Hadjinicolaou (1993) as adopted by Abel et al. (2010).

(11) Negligible radiation and chemical reaction — thermal radiation and species

transport effects are not considered.

3.3 Governing Partial Differential Equations

4
 Under the boundary-layer approximations and the modelling assumptions of Section

3.2, the governing equations for the conservation of mass, momentum, and energy for the

second-grade fluid are as follows (Liu, 2005; Abel et al., 2010).

3.3.1 Continuity Equation

The conservation of mass for an incompressible fluid requires that the velocity field

be divergence-free:

∂u/∂x + ∂v/∂y = 0 (3.1)

Equation (3.1) guarantees that no fluid is created or destroyed within the domain. It is

automatically satisfied by introducing a stream function ψ(x, y) such that u = ∂ψ/∂y and v =

−∂ψ/∂x, as described in Section 3.5.

3.3.2 Momentum Equation

The x-momentum equation for a second-grade fluid within the boundary-layer

approximation, including the Darcy drag term for the porous medium, is (Liu, 2005):

u ∂u/∂x + v ∂u/∂y = ν ∂²u/∂y² + (α₁/ρ)[u ∂³u/∂x∂y² − (∂u/∂y)(∂²u/∂x∂y) (3.2)

+ v ∂³u/∂y³] − (ν/k')u

5
 The first term on the right-hand side, ν ∂²u/∂y², represents the standard Newtonian

viscous diffusion of momentum in the wall-normal direction. The bracketed group of terms,

multiplied by α₁/ρ, represents the elastic contribution of the second-grade fluid: the three

terms within the bracket arise from the projection of the upper-convected derivative of the

first Rivlin-Ericksen tensor onto the x-momentum balance, and they introduce the third-order

derivative ∂³u/∂y³, which raises the order of the ODE from three (Newtonian case) to four

(second-grade case). This higher order necessitates an additional boundary condition, which

in the present problem is supplied by the far-field condition f''(∞) → 0. The final term,

−(ν/k')u, is the Darcy drag force exerted by the porous medium on the fluid in the streamwise

direction; it opposes the flow and scales with the local streamwise velocity u.

3.3.3 Energy Equation

The energy equation for the second-grade fluid, incorporating convective heat

transport, thermal diffusion, viscous dissipation, and the non-uniform internal heat source, is

(Abel et al., 2010):

u ∂T/∂x + v ∂T/∂y = α ∂²T/∂y² + (ν/c_p)(∂u/∂y)² + (1/ρc_p)(∂u/∂y)[u (3.3)

∂u/∂x + v ∂u/∂y] + q'''/ρc_p

In Equation (3.3), α = κ/(ρc_p) is the thermal diffusivity of the fluid (where κ is the

thermal conductivity), and the four terms on the right-hand side represent: (i) thermal

diffusion in the wall-normal direction; (ii) viscous dissipation due to the primary shear ∂u/∂y;

6
(iii) the additional viscous dissipation term arising from the interaction of the elastic stress

and the convective acceleration (sometimes called the viscoelastic cross-term); and (iv) the

volumetric heat generation/absorption rate q''', to be specified in Section 3.3.4. The energy

equation (3.3) is coupled to the momentum equation (3.2) through the velocity field (u, v).

3.3.4 Non-Uniform Internal Heat Source/Sink Model

The non-uniform internal heat generation or absorption term q''' is modelled

following Vajravelu and Hadjinicolaou (1993) as adopted by Abel et al. (2010). The

volumetric heat source rate is decomposed into two physically distinct contributions:

q''' = (κu_w / νx) [ A*(T_w − T_∞) f'(η) + B*(T − T_∞) ] (3.4)

The first term in Equation (3.4), proportional to the dimensionless velocity gradient

f'(η) with coefficient A*, is called the space-dependent (or velocity-gradient-dependent) heat

source/sink. Since f'(η) is large near the sheet and decays exponentially to zero as η → ∞,

this term represents heat generation or absorption that is concentrated near the wall —

mimicking, for example, chemical reactions or metabolic processes that are activated by the

local shear rate. The second term, proportional to the local temperature excess (T − T_∞)

with coefficient B*, is called the temperature-dependent heat source/sink. Since it scales with

the local temperature, it is distributed throughout the boundary layer in proportion to the

thermal excess, analogous to a volumetric reaction whose rate depends on the local

temperature. Positive values of A* and B* correspond to heat generation (sources), which

7
elevate the fluid temperature; negative values correspond to heat absorption (sinks), which

depress it.

3.4 Dimensional Boundary Conditions

The boundary conditions for the velocity field are prescribed as follows. At the sheet

surface y = 0, the no-slip condition requires that the tangential velocity equal the sheet

velocity, and the wall-normal velocity equals the (inward) suction velocity:

u(x, 0) = u_w = cx, v(x, 0) = −v₀, T(x, 0) = T_w (3.5)

In Equation (3.5), the negative sign on −v₀ reflects the convention that v₀ > 0 denotes

suction (fluid drawn into the wall). In the far field, the velocity must decay to the quiescent

ambient value and the temperature must return to the ambient temperature:

u(x, y) → 0 as y → ∞, T(x, y) → T_∞ as y → ∞ (3.6)

The additional condition required by the fourth-order nature of the momentum

equation (due to the second-grade term) is that the shear-rate gradient must vanish in the far

field:

∂²u/∂y² → 0 as y → ∞ (3.7)

8
 Physically, Equation (3.7) states that the viscous stress gradient must vanish at the

edge of the boundary layer, ensuring that the second-grade elastic forces do not persist into

the quiescent far field. This condition is automatically satisfied by the exponential velocity

solution, as shown in Section 3.7.

3.5 Similarity Transformation

3.5.1 Motivation and Stream Function

The governing PDEs (3.1)–(3.3) constitute a system in three unknown fields: u(x, y),

v(x, y), and T(x, y). A direct numerical solution of this PDE system on the semi-infinite

domain {y > 0} would require discretisation in both x and y and the specification of far-field

boundary conditions at large but finite y. The existence of a similarity structure in the

problem — arising from the linear stretching velocity law u_w = cx — permits a dramatic

reduction: the three PDEs can be collapsed into a coupled system of two ODEs in a single

similarity variable η, which is far more tractable both analytically and numerically.

To satisfy the continuity equation (3.1) identically, a stream function ψ(x, y) is

introduced such that u = ∂ψ/∂y and v = −∂ψ/∂x. Dimensional analysis of the stretching-sheet

boundary-layer problem suggests that the stream function takes the self-similar form:

9
 ψ(x, y) = √(cν) · x · f(η), where η = y√(c/ν) (3.8)

Here η is the dimensionless wall-normal coordinate (the similarity variable), and f(η)

is the dimensionless stream function. The scaling √(c/ν) ensures that η is of order unity

within the boundary layer, whose thickness scales as √(ν/c). From Equation (3.8), the

velocity components are obtained by differentiation:

u = ∂ψ/∂y = cx f'(η), v = −∂ψ/∂x = −√(cν) f(η) (3.9)

where primes denote differentiation with respect to η. The boundary condition u(x, 0)

= cx is satisfied if f'(0) = 1, and the suction condition v(x, 0) = −v₀ is satisfied if −√(cν) f(0)

= −v₀, i.e., f(0) = v₀/√(cν) = R, where R is the dimensionless suction parameter defined in

Equation (3.12) below.

3.5.2 Dimensionless Temperature

The dimensionless temperature is defined as:

θ(η) = (T − T_∞) / (T_w − T_∞) (3.10)

so that θ(0) = 1 at the wall and θ(∞) = 0 in the far field. The definition (3.10) is

consistent with the assumption T_w > T_∞ (heated wall); for a cooled wall (T_w < T_∞), the

same definition applies with θ becoming negative.

10
3.6 Dimensionless Governing Parameters

Substituting the similarity variables (3.8)–(3.10) into the momentum and energy

equations yields five dimensionless parameters that govern the solution. Their definitions and

physical interpretations are as follows:

λ₁ = α₁c/(ρν), λ₂ = ν/(ck'), Pr = ν/α, Ec = c²x²/[c_p(T_w − T_∞)], (3.11)

R = v₀/√(cν)

Note that the Eckert number in (3.11) involves x², which would make it position-

dependent unless the wall-temperature difference is also position-dependent. In the similarity

solution framework, the temperature difference is taken to scale as c²x²/c_p, so that Ec is a

global constant characterising the overall importance of viscous dissipation relative to

convective heat capacity. The five parameters are described in Table 3.1 below.

Table 3.1 Dimensionless governing parameters: definitions, physical interpretation, and

representative ranges

Parameter Sym Definition Physical Meaning Typical

bol Range

Second-grade λ₁ α₁c/(ρν) Ratio of elastic to viscous 0.5 – 2.0

(viscoelastic) forces; λ₁ = 0 recovers
Newtonian flow

Porosity λ₂ ν/(ck') Ratio of Darcy resistance to 0.5 – 3.0

11
 Parameter Sym Definition Physical Meaning Typical
bol Range

stretching inertia; larger λ₂

means lower permeability

Prandtl number Pr ν/α Ratio of momentum to thermal 0.72 – 3.0

diffusivity; governs relative
thickness of velocity and
thermal BLs

Eckert number Ec u_w²/ Ratio of kinetic energy to 0.0 – 2.0

[c_p(T_w−T_∞) thermal energy; measures
] viscous dissipation strength

Suction R v₀/√(cν) Dimensionless wall suction 0.5 – 3.0

velocity; R > 0 thins both
boundary layers

Heat source A* Defined in (3.4) Controls wall-concentrated heat −0.5 – 1.0

(space) generation/absorption (positive
= source)

Heat source B* Defined in (3.4) Controls temperature- −0.5 – 1.0

(temp.) proportional distributed heat
generation/absorption

3.7 Reduced Ordinary Differential Equation System

Substituting the similarity variables (3.8)–(3.10) and the definitions (3.11) into the

governing PDEs (3.2)–(3.3) and performing the algebra of the boundary-layer transformation

yields the following coupled ODE system.

3.7.1 Fourth-Order Momentum ODE

12
 After substitution and cancellation of the common factor cx from all terms, the

momentum equation (3.2) reduces to the following fourth-order nonlinear ODE for f(η):

(f')² − ff'' + λ₂f' = f''' + λ₁(2f'f''' − (f'')² − ff'''') (3.12)

Equation (3.12) is a fourth-order nonlinear ODE in f(η). The left-hand side contains

the inertial terms (f')² − ff'' (Crane's Newtonian terms) plus the Darcy drag λ₂f'. The right-

hand side contains the Newtonian viscous term f''' augmented by the three second-grade

elastic terms: 2f'f''' (the interaction of local velocity and third-order shear), −(f'')² (the square

of the shear rate), and −ff'''' (the interaction of the stream function and the fourth derivative).

These three terms arise from the projection of the upper-convected derivative of A₁ onto the

boundary-layer momentum balance.

It is important to note that Equation (3.12) is of order four (due to the f'''' term),

whereas the corresponding Newtonian equation is of order three. This requires one additional

boundary condition beyond those in the Newtonian problem; the condition f''(∞) → 0 serves

this purpose, and it is automatically satisfied by the exact exponential solution derived in

Section 3.8.

3.7.2 Second-Order Energy ODE

13
 After substitution of the similarity variables into the energy equation (3.3) and the

heat source model (3.4), the energy equation reduces to the following second-order ODE for

θ(η):

θ'' + Pr·f·θ' + A*·θ = −[ B*·f' + Ec·Pr·{ (f'')² + λ₂(f')² } ] (3.13)

On the left-hand side of Equation (3.13): θ'' is the thermal diffusion term; Pr·f·θ' is the

convective transport of heat in the wall-normal direction by the flow f(η); and A*·θ is the

temperature-dependent heat source/sink. On the right-hand side: B*·f' is the space-dependent

(velocity-gradient-driven) heat source/sink; and the term Ec·Pr·{(f'')² + λ₂(f')²} represents the

total viscous dissipation — comprising (f'')² from the primary shear and λ₂(f')² from the

Darcy dissipation in the porous medium (following Al-Hadrami et al., 2002). The energy

equation (3.13) is coupled to the momentum equation (3.12) through the velocity field f(η)

and its derivatives.

3.7.3 Transformed Boundary Conditions

In terms of the similarity variables, the boundary conditions (3.5)–(3.7) become:

For the momentum equation:

f(0) = R, f'(0) = 1, f'(η) → 0 as η → ∞, f''(η) → 0 as η → ∞ (3.14)

14
For the energy equation:

θ(0) = 1, θ(η) → 0 as η → ∞ (3.15)

The boundary conditions (3.14) provide four conditions for the fourth-order ODE

(3.12): two at η = 0 (the wall value f(0) = R and the velocity slope f'(0) = 1) and two in the

far field (f'(∞) = 0 and f''(∞) = 0). The two conditions at η = 0 are used as initial values for

forward integration, while the two far-field conditions are enforced through a shooting

procedure, as described in Chapter 4. For the energy ODE (3.13), the wall condition θ(0) = 1

is known, and the far-field condition θ(∞) = 0 is enforced by determining the unknown initial

slope θ'(0) = β through the shooting algorithm.

3.8 Exact Analytical Solution for the Velocity Field

3.8.1 Proposed Ansatz

Following Liu (2005), we propose that the dimensionless velocity f'(η) takes the

exponential form:

f'(η) = e^{−sη} (3.16)

15
 where s > 0 is a real positive constant to be determined. The ansatz (3.16) satisfies the

far-field boundary condition f'(∞) = 0 for any s > 0, and satisfies f'(0) = 1 exactly. Integrating

(3.16) and applying the initial condition f(0) = R gives:

f(η) = R + (1 − e^{−sη}) / s (3.17)

The higher derivatives of f required for substitution into Equation (3.12) are:

f'(η) = e^{−sη}, f''(η) = −s e^{−sη}, f'''(η) = s² e^{−sη}, f''''(η) = (3.18)

−s³ e^{−sη}

Note that f''(∞) = lim_{η→∞}(−se^{−sη}) = 0 for any s > 0, so the additional far-

field condition f''(∞) = 0 is automatically satisfied by the ansatz.

3.8.2 Derivation of the Corrected Cubic Equation

We now substitute the expressions (3.17)–(3.18) directly into the momentum ODE

(3.12), computing each term explicitly. This direct substitution — without any intermediate

approximation — yields the corrected cubic equation for s.

Step 1: Evaluate each term on the left-hand side of Equation (3.12).

(f')² = e^{−2sη}

16
 ff'' = [R + (1−e^{−sη})/s](−se^{−sη}) = −sRe^{−sη} − e^{−sη} + e^{−2sη}

λ₂f' = λ₂e^{−sη}

Therefore the left-hand side of (3.12) is:

LHS = e^{−2sη} − (−sRe^{−sη} − e^{−sη} + e^{−2sη}) + λ₂e^{−sη}

= e^{−2sη} + sRe^{−sη} + e^{−sη} − e^{−2sη} + λ₂e^{−sη}
= (sR + 1 + λ₂) e^{−sη}

Step 2: Evaluate each term on the right-hand side of Equation (3.12).

f''' = s²e^{−sη}

2f'f''' = 2e^{−sη} · s²e^{−sη} = 2s²e^{−2sη}

(f'')² = (−se^{−sη})² = s²e^{−2sη}

ff'''' = [R + (1−e^{−sη})/s](−s³e^{−sη}) = −s³Re^{−sη} − s²e^{−sη} + s²e^{−2sη}

Therefore the second-grade bracket in the right-hand side of (3.12) is:

2f'f''' − (f'')² − ff'''' = 2s²e^{−2sη} − s²e^{−2sη} − (−s³Re^{−sη} − s²e^{−sη} +

s²e^{−2sη})
= 2s²e^{−2sη} − s²e^{−2sη} + s³Re^{−sη} + s²e^{−sη} − s²e^{−2sη}
= (2s² − s² − s²)e^{−2sη} + (s³R + s²)e^{−sη}

17
 = 0 · e^{−2sη} + s²(sR + 1)e^{−sη}
= s²(sR + 1)e^{−sη}

This is a crucial intermediate result: all e^{−2sη} terms cancel exactly in the second-

grade bracket, leaving only a term proportional to e^{−sη}. Therefore the right-hand side of

(3.12) is:

RHS = s²e^{−sη} + λ₁ · s²(sR + 1)e^{−sη}

= s²[1 + λ₁(sR + 1)]e^{−sη}
= s²(1 + λ₁ + λ₁sR)e^{−sη}

Step 3: Set LHS = RHS and cancel the common factor e^{−sη} > 0:

sR + 1 + λ₂ = s²(1 + λ₁ + λ₁sR)

Expanding the right-hand side and rearranging all terms to one side:

λ₁Rs³ + (1 + λ₁)s² − Rs − (1 + λ₂) = 0

This is the corrected cubic equation for the decay constant s. Dividing through by R

(which is valid since R > 0 for suction) gives the equivalent form:

λ₁s³ + [(1 + λ₁)/R] s² − s − (1 + λ₂)/R = 0 (3.19)

18
 Equation (3.19) is a cubic in s with real coefficients. For the physically relevant

parameter ranges λ₁ ≥ 0, λ₂ > 0, R > 0, it can be shown that there exists exactly one positive

real root s > 0, which is the physically meaningful solution ensuring f'(η) → 0 as η → ∞.

This positive root is computed numerically using Cardano's formula or Newton's method for

each combination of (λ₁, λ₂, R). The exact values of s for the parameter sets used in Chapter

5 are tabulated in Table 3.2.

3.8.3 Correction of the Literature Discrepancy

The cubic equation reported by Abel et al. (2010) and reproduced by Rashidi et al.

(2011) for the case λ₁ = 0.5, λ₂ = 2, R = 1 yields s ≈ 1.8026, whereas the corrected Equation

(3.19) yields s = √2 ≈ 1.41421. The source of this discrepancy is a subtle algebraic

inconsistency in the derivation adopted by those authors, which can be traced as follows.

In Abel et al. (2010) and Rashidi et al. (2011), the stream function is assumed to take

the form f(η) = A(1 − e^{−sη}), and the boundary condition f'(0) = Ase^0 = As = 1 is used to

set A = 1/s. This approach, however, yields f(0) = 0, which is inconsistent with the correct

boundary condition f(0) = R. The suction parameter R is effectively dropped from the stream

function, so the initial condition f(0) = R is not enforced. When this incomplete stream

function f(η) = (1 − e^{−sη})/s (with A = 1/s and f(0) = 0) is substituted into the ODE, the

f(0) term is absent from the expressions for ff'' and ff'''', producing a different set of

coefficients in the resulting algebraic equation. The cubic that emerges from this incomplete

19
substitution is not the correct characteristic equation for the ODE with f(0) = R, and its

positive root does not generally equal the physically correct value of s.

The present derivation avoids this inconsistency entirely by including the suction

parameter R in the stream function from the outset: f(η) = R + (1 − e^{−sη})/s, which

correctly satisfies both f(0) = R and f'(0) = 1. Direct substitution of this complete expression

into (3.12) — as carried out in Steps 1–3 of Section 3.8.2 — yields the corrected cubic (3.19)

without any approximation. The correctness of (3.19) can be verified by substituting f'(η) =

e^{−sη} back into Equation (3.12) after computing s from (3.19); the residual is zero to

machine precision for all parameter combinations tested.

Table 3.2 compares the corrected values of s obtained from Equation (3.19) with the

values reported by Abel et al. (2010) and Rashidi et al. (2011) for selected parameter

combinations. For the case λ₁ = 1.0, the two approaches agree (because both give the same

cubic when R = 1), but for λ₁ = 0.5 the discrepancy is clearly evident: the corrected value s =

√2 ≈ 1.41421 differs significantly from the reported value s ≈ 1.8026.

Table 3.2 Comparison of the decay constant s: corrected Equation (3.19) vs. Abel et al.

(2010) (λ₂ = 2, R = 1)

λ₁ λ₂ R s (Corrected, Eq. s (Abel et al., f''(0) Corrected f''(0) Abel et al.

3.19) 2010)

0.5 2.0 1.0 √2 ≈ 1.41421 ≈ 1.8026 −1.41421 −1.8026

1.0 2.0 1.0 1.14789904 1.14789904 −1.14789904 −1.14789904

1.5 2.0 1.0 1.00000000 1.00000000 −1.00000000 −1.00000000

20
 λ₁ λ₂ R s (Corrected, Eq. s (Abel et al., f''(0) Corrected f''(0) Abel et al.
3.19) 2010)

2.0 2.0 1.0 0.90129398 0.90129398 −0.90129398 −0.90129398

The fact that all rows with λ₁ ≥ 1.0 agree between the two approaches arises because,

for R = 1, the dominant λ₁ term correctly dominates the cubic and the suction-related

asymmetry is less influential. The discrepancy is most pronounced at λ₁ = 0.5 precisely

because the omission of the R-dependent terms in the Abel et al. derivation has a relatively

larger effect on the cubic when the λ₁ coefficient is small.

3.8.4 Verification by Direct Substitution

To verify Equation (3.19) independently, we substitute f'(η) = e^{−sη} — with s the

positive root of (3.19) — back into the momentum ODE (3.12) and compute the residual

R(η) = LHS − RHS at η = 0, 1, 2, ..., 8. For the base parameter set (λ₁ = 1, λ₂ = 2, R = 1, s =

1.14789904), the residual is:

R(η) = (sR + 1 + λ₂)e^{−sη} − s²(1 + λ₁ + λ₁sR)e^{−sη}

= [(sR + 1 + λ₂) − s²(1 + λ₁ + λ₁sR)] e^{−sη}
= [λ₁Rs³ + (1+λ₁)s² − Rs − (1+λ₂)] · (−1/R) · e^{−sη} = 0

since s satisfies Equation (3.19). The residual is identically zero for all η, confirming

that the exact solution (3.16)–(3.17) satisfies the momentum ODE (3.12) exactly — not

merely approximately — at every point in the domain. This property makes it an exact

21
benchmark in the strict mathematical sense: it is not a truncated series, a perturbation

approximation, or a numerical solution, but a genuine closed-form solution.

3.8.5 Skin-Friction Coefficient

The local skin-friction coefficient C_f is defined as:

C_f Re_x^{1/2} = 2(1 + 3λ₁) f''(0) (3.20)

From the exact solution (3.18), f''(0) = −s, so:

C_f Re_x^{1/2} = −2(1 + 3λ₁) s (3.21)

Since s > 0, the skin-friction coefficient is negative, consistent with the convention

that the wall exerts a retarding force on the fluid in the direction opposite to the stretching

motion. The magnitude |f''(0)| = s is an exact, parameter-dependent quantity that serves as the

primary benchmark for the velocity solution throughout Chapter 5. Table 3.3 below tabulates

the exact values of s and f''(0) for the parameter combinations investigated in the parametric

study.

Table 3.3 Exact values of the decay constant s and skin-friction parameter f''(0) = −s for

all parameter combinations (λ₂ = 2, R = 1 base)

22
 λ₁ λ₂ R s (positive root of Eq. f''(0) = −s (Exact)
3.19)

0.5 2.0 1.0 1.41421356 −1.41421356

1.0 2.0 1.0 1.14789904 −1.14789904

1.5 2.0 1.0 1.00000000 −1.00000000

2.0 2.0 1.0 0.90129398 −0.90129398

1.0 1.0 1.0 1.00000000 −1.00000000

1.0 3.0 1.0 1.26953100 −1.26953100

1.0 2.0 2.0 1.11208500 −1.11208500

1.0 2.0 3.0 1.09072200 −1.09072200

3.9 Temperature Equation: Absence of Closed-Form Solution

Unlike the momentum equation (3.12), the energy equation (3.13) does not admit a

known exact closed-form solution for general values of the governing parameters Pr, Ec, A*,

and B*. To understand why, rewrite (3.13) in the standard form of a second-order linear ODE

with variable coefficients:

θ'' + P(η)θ' + Q(η)θ = G(η) (3.22)

where the coefficient functions are:

P(η) = Pr·f(η) = Pr·[R + (1−e^{−sη})/s] (3.23)

23
 Q(η) = A* (3.24)

G(η) = −B*·e^{−sη} − Ec·Pr·[ s²e^{−2sη} + λ₂e^{−2sη} ] = (3.25)

−B*e^{−sη} − Ec·Pr·(s² + λ₂)e^{−2sη}

The non-constant coefficient P(η) = Pr·f(η), which involves the exponential function

e^{−sη} embedded in f(η), prevents the application of standard constant-coefficient solution

techniques. The equation is not of Euler-Cauchy type, Bessel type, or Legendre type, and it

does not reduce to any standard form in the classical catalogues of exactly solvable ODEs

(Polyanin and Zaitsev, 2003). Moreover, the non-homogeneous forcing term G(η) contains

two exponential components at different decay rates (e^{−sη} and e^{−2sη}), which further

complicates the search for a particular integral.

For the special case A* = 0, B* = 0, Ec = 0 (no heat source and no viscous

dissipation), the energy equation (3.22) reduces to θ'' + Pr·f·θ' = 0, which is related to the

confluent hypergeometric equation and can be expressed as an incomplete gamma function

integral (Grubka and Bobba, 1985). However, for general non-zero values of A*, B*, and Ec

— as considered in the present study — no such reduction is available, and the energy

equation must be solved by numerical or semi-analytical methods.

The approach adopted in the present work is twofold. For validation purposes, a high-

precision numerical solution is obtained using the fourth-order Runge-Kutta scheme with

adaptive step-size control (DOP853, rtol = 10⁻¹², atol = 10⁻¹⁴) combined with a Brent

bracketing shooting algorithm to enforce θ(η_max) = 0. This numerical solution serves as the

24
temperature benchmark θ_RK4(η). For the primary semi-analytical solution, the Multistep

Differential Transform Method (MsDTM) is applied to Equation (3.13), as described in detail

in Chapter 4. The MsDTM produces a piecewise polynomial approximation to θ(η) that

matches θ_RK4(η) to better than 10⁻¹³ throughout the domain [0, 8], as demonstrated in

Chapter 5.

3.10 Physical Interpretation of the Exact Velocity Solution

The exact velocity field f'(η) = e^{−sη} has a number of physically important

properties that are worth noting before proceeding to the solution methodology in Chapter 4.

First, the exponential decay of f'(η) with η confirms that the velocity boundary layer

has an exponential rather than algebraic thickness. The effective boundary-layer thickness

δ_v, defined as the value of η at which f'(η) = 0.01 (i.e., the velocity has decayed to 1% of the

wall value), is δ_v = −ln(0.01)/s ≈ 4.605/s. For the base parameters (λ₁ = 1, λ₂ = 2, R = 1, s =

1.14790), this gives δ_v ≈ 4.01, consistent with the computational domain η_max = 8 being

more than sufficient.

Second, since s depends on λ₁, λ₂, and R through the cubic (3.19), the velocity

profile is parametrically sensitive to all three parameters. Specifically: (a) increasing λ₁

(more elastic fluid) decreases s, producing a thicker velocity boundary layer with slower

decay — the elastic stresses resist velocity-gradient changes and retard momentum diffusion;

25
(b) increasing λ₂ (lower porous permeability) increases s, thinning the velocity boundary

layer — the Darcy drag dissipates momentum closer to the wall; (c) increasing R (stronger

suction) has a relatively modest effect on s for the parameter ranges of interest, as

demonstrated in Chapter 5.

Third, the exact solution shows that the full stream function is f(η) = R + (1 −

e^{−sη})/s, which approaches the constant f(∞) = R + 1/s as η → ∞. This finite limiting

value of the stream function is consistent with the far-field condition f'(∞) = 0 (since f'(η) →

0 as η → ∞) and confirms that the normal velocity component v = −√(cν)f(η) approaches the

constant −√(cν)(R + 1/s) in the far field, representing the entrainment of fluid into the

boundary layer.

3.11 Chapter Summary

This chapter has presented the complete mathematical formulation of the second-

grade fluid stretching-sheet problem. The governing PDEs — continuity (3.1), momentum

(3.2), and energy (3.3) — together with the non-uniform heat source model (3.4) and the

boundary conditions (3.5)–(3.7), were reduced via the similarity transformation (3.8)–(3.10)

to the coupled ODE system comprising the fourth-order nonlinear momentum ODE (3.12)

and the second-order energy ODE (3.13), subject to the boundary conditions (3.14)–(3.15).

Seven dimensionless governing parameters were identified and their physical meanings

26
explained: the second-grade parameter λ₁, the porosity parameter λ₂, the Prandtl number Pr,

the Eckert number Ec, the suction parameter R, and the heat source coefficients A* and B*.

The exact closed-form velocity solution f'(η) = e^{−sη} was derived by direct

substitution of the complete stream function f(η) = R + (1 − e^{−sη})/s into the momentum

ODE (3.12). The derivation yielded the corrected cubic equation (3.19) for the positive

constant s, which was shown to differ from the cubic reported by Abel et al. (2010) and

Rashidi et al. (2011) due to the omission of the suction parameter R from the stream function

in those works. The corrected cubic was verified by direct substitution — the residual of the

ODE is identically zero — and the discrepancy was quantified: for λ₁ = 0.5, λ₂ = 2, R = 1,

the corrected value is s = √2 ≈ 1.41421, compared with s ≈ 1.8026 in the prior literature. For

all other parameter combinations considered, the exact values of s and the skin-friction

parameter f''(0) = −s were tabulated.

The energy equation (3.13) was shown to have no known exact closed-form solution

for general parameter values. Its solution requires semi-analytical or numerical methods, and

the MsDTM approach adopted in this thesis is described fully in Chapter 4. The physical

properties of the exact velocity solution were discussed, revealing the parametric dependence

of the boundary-layer thickness on λ₁, λ₂, and R through the decay constant s.

27