EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
Vol. 16, No. 2, 2023, 983-996
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

Mathematical Model for the effect of buoyancy forces
on the stability of a fluid flow

Maysoon M. Aziz

Department of Mathematics, College of Computer Sciences and Mathematics,
University of Mosul, Mosul, Iraq

Abstract. Stability analysis of heat transfer (by Conduction, Convection, and Radiation) has
been found for a model of flow between two horizontal plates, one of them is thermally insulated.
The stability measure through the neutral curve for this model shows that an increase of the
buoyancy represented by Gershoff number (Gr) leads to an increase in stability region, especially
for the considerable value of Reynolds number R. The results indicated that the effects of buoyancy
forces have significant contribution to the field profiles

2020 Mathematics Subject Classifications: 76E06, 80A20, 35B35, 76M10

Key Words and Phrases: Stability, Heat Transfer, Buoyancy, Gershoff Number

1. Introduction

The field of mathematical physics known as stability theory allows for the derivation
of equation-solving statements from the principle of critical values that distinguishes the
many regimes of flow as well as the types of fluid motion that occur in each of these
regimes. Chaotic behavior and heat transfer performance are represented very accurately
[4–6].

Any system whatsoever can be disturbed, and the question is (does that disturbance
die down or grows up with time-lapse) [7, 8].

The presence of buoyancy forces convection along a flat plate causes a coupling between
the momentum and energy equations and, at the same time, causes the boundary layer to
become non-similar [2, 3, 9, 10, 13].

Mori, Y. [11] has discussed the effects of buoyancy forces in a forced laminar convection
flow over a horizontal plate, it is shown that the solution may be expanded into power
series in Gr / Re whose first terms express solution for the purely forced convection flow.

A bu-Mulaweh H.I. and others [1] have discussed the measurements of velocity and
temperature distributions and reported for buoyancy – opposing, laminar, mixed convec-
tion flow over a vertical back ward – facing step.

DOI: https://doi.org/10.29020/nybg.ejpam.v16i2.4751

Email addresses: aziz maysoon@uomosul.edu.iq (M. M. Aziz)

https://www.ejpam.com 983 © 2023 EJPAM All rights reserved.



M. M. Aziz / Eur. J. Pure Appl. Math, 16 (2) (2023), 983-996 984

Maulucci S.V. [12] has discussed the characteristics of spatial type waves in the bound-
ary layer of incompressible fluid on a flat plate at high Reynolds number R.

The stability of a model of heat transfer by conduction and convection through a
porous medium has been investigated by Aziz M. al [13]. The resulting analysis shows
that the imposed disturbances are dying down, and the system under consideration is
always stable.

The present paper is to investigate the effect of buoyancy forces one of stability of fluid
flow between two plates on them is thermally insulated.

2. Formulation of the model and solution

Considering a fluid moving laminarly between two parallel horizontal wall separations
of 2h distance, one of them is thermally insulated.

Choosing the coordinate framework. The x-axis is parallel to the channel and along the
direction of the flow, and the y-axis is taken as the vertical coordinate measured position
upward, wilts the z-axis is the direction mutually orthogonal to the other two axes.

The basic differential equation governing such a model can be summarized as follows:

2.1. Stat Equation

ρ = ρ0[1− β1(T − T1)]. (1)

Where β1 is the coefficient of thermal expansion, T1 is the characteristic temperature,
and ρ is the density.

2.2. Continuity Equation

A small size component by taking and calculating the force balance condition from
which the continuity equation can be deduced in the general case in which the flow is
compressive and unsteady flow in three -dimensions:

∂ρ

∂t
+

∂(ρu)

∂x
+

∂(ρv)

∂y
+

∂(ρw)

∂z
= 0. (2)

Where is the density and u, v, w, are the velocity components in the x, y, and z direc-
tions, respectively, which are functions of the coordinates x, y, z and time t.

2.3. Motion Equation

From the second law of motion, the motion Equation

ρ(v⃗ grad v⃗ + g⃗) = µ∇v⃗– grad p. (3)

Where µ viscosity, p pressure and g gravitational acceleration.
From the law of conservation of energy, the energy equation can be formulated as

follow:



M. M. Aziz / Eur. J. Pure Appl. Math, 16 (2) (2023), 983-996 985

2.4. Energy Equation

ρ

(
cv
∂T

∂t
+ v⃗ grad T

)
= k∇2T– div q⃗ + φ. (4)

Where cv specific heat at a constant temperature, T temperature, k coefficient of
thermal conductivity and φ is viscous dissipation, and q⃗ = i⃗ qx + j⃗ qy + k⃗qz is Radiative
flux vector and its Component are qx , qy , qz.

cpρ0
∂T

∂t
= k

∂2T

∂y2
+ 4ασ

(
TT 3

1 − 3T 4
1

)
+ µ

(
∂u

∂y

)2

. (5)

ρ0
∂u

∂t
= µ

∂2u

∂y2
+ ρg. (6)

Where ρ = ρ0[1−β1(T −T1)], cp specific heat at constant pressure, α is the obsorption
coefficient and σ Stefan constant.
Boundary condition of the model are:

u = 0 y = h
dT

dy
= 0

u = 0 y = −h T = T1

. (7)

To convert the governing partial differential equations into non-dimensional, we intro-
duce the following non-dimensional quantities:

T = θT1, u = u0∪, y = hη, t = t0t
∗. (8)

Pr =
µcp
k

(Prandtl number),

Ec =
u20
cpT1

(Eckert number),

Re =
hρ0u0
µ

(Reynolds number),

Bo =
ρ0u

3
0

T 4
1 σ

(Boltzmann number),

W = hα (Bouguer number),

Gr =
gρ20B1h

3

µ2
(Gershoff number),

Fr =
u20
gh

(Frood number).

and

A =
1

Pr. R
,C =

W . Ec

Bo
,D =

Ec

R
,E =

1

R
,F =

Gr

R2
, H =

1

Fr
+ F

Equations (5) and (6) become)

∂θ

∂t∗
= A

∂2θ

∂η2
–4Cθ +D(

∂u

∂η
)
2

–12C. (9)



M. M. Aziz / Eur. J. Pure Appl. Math, 16 (2) (2023), 983-996 986

∂u

∂t∗
= E

∂2u

∂η2
–Fθ +H. (10)

The boundary conditions become

∪ = 0 η = −1 θ = 1

∪ = 0 η = 1
dθ

dη
= 0

. (11)

3. Steady state solution

To find the steady state solution, we shall neglect the terms which depend on the time
and equations (9) and (10) become

A
∂2θ

∂η2
–4Cθ +D(

d∪
dη

)
2

–12C = 0. (12)

E
d2∪
dη2

–Fθ +H = 0. (13)

From equations (12) and (13) we get the following equation

d4∪
dη4

+
4C

A
(
d2∪
dη2

) + 4
CE

AE
(
H

F
− 3) +

DF

AE
(
dU

dη
)
2

= 0. (14)

and the boundary conditions become

∪ = 0 θ = 1 η = −1

,∪ = 0
dθ

dη
= 0 η = 1

,
d2∪
dη2

=
F −H

E
η = −1

,
d3u

dη3
= 0 η = 1.


(15)

Now we shall use perturbation method [8] and assume the solution of (14) to be as
follows.

∪(η) =
∞∑
n=0

∈n ∪n(η). (16)

Since D = EC
R and the values of EC from 0.1 to 0.01 and R is greater than 100, then

the maximum value of D is 0.001, thus we assume ∈= D we shall neglect terms which in
it Dn where n > 2 and thus the solution has the form

∪(η) = ∪0(η) +D ∪1 (η). (17)



M. M. Aziz / Eur. J. Pure Appl. Math, 16 (2) (2023), 983-996 987

Differentiate equation (17) many times and compensate it in equation (14) we get.

d4∪0

dη4
+

4C

A
(
d2∪0

dη2
) + 4

CF

AE
(
H

F
− 3) = 0. (18)

d4∪1

dη4
+

4C

A
(
d2∪1

dη2
) +

F

AE
(
d∪0

dη
)
2

= 0. (19)

and Boundary conditions

∪0 = 0 η = −1
,∪0 = 0 η = 1

,
d2∪0

dη2
=

F −H

E
η = −1

,
d3∪0

dη3
= 0 η = 1.


(20)

∪1 = 0 η = −1
,∪1 = 0 η = 1

,
d2∪1

dη2
= 0 η = −1

,
d3∪0

dη3
= 0 η = 1 .


(21)

Now we found the general solution of equation (18) by found complement solution at
homogenous equation of the following form.

∪op = A1 +A2 η +A3 cos(2S1η) +A4sin(2S1η). (22)

We found special solution of non-homogenous equation by method of variational pa-
rameters [7] which is ∪oc

∪oc(η) = Asη
2. (23)

Where A5 =
−F

2E
(
H

F
− 3), S =

C

A
,S1 = S

1
2

and using boundary conditions we found the constant A1, A2, A3, A4, where

A1 =
F

2E
(
H

F
− 3)− FA

2CE
[

cos2(2S1)

2cos2 (2S1)− 1
]

A2 =
FA

2CE
[

sin2 (2S1)

2cos2 (2S1)− 1
]

A3 =
FA

2CE
[

cos(2S1)

2cos2 (2S1)− 1
]

A4 =
FA

2CE
[

sin(2S1)

2cos2 (2S1)− 1
]

∪1p(η) = B1 +B2(η) +B3cos(2S1η) +B4sin(2S1η). (24)



M. M. Aziz / Eur. J. Pure Appl. Math, 16 (2) (2023), 983-996 988

and we found ∪1C(η) by method of variation which has the form

U1c(η) =
1

4
S−1 F

AE

[
h1η

2 + h2η
3 + h3η

4 + h4η cos (2S1η) + h5η sin (2S1η)+ h6 sin
2

(2S1η) + h7η
2 cos (2S1η) + h8η

2 sin (2S1η) + h9η sin
2 (2S1η) + h10 sin (2S1η) cos

(2S1η) + h11 sin (2S1η) cos
2 (2S1η) + h12 sin

2 (2S1η) cos (2S1η) + h13sin(2S1η)

cos2(2S1η) + h14ηsin
2 (2S1η) cos(2S1η).

(25)

Where

h1 =
1

2
A2

2 + 3SA2
4 − SA2

3 − S−1A2
5

h2 =
2

3
A2A5

h3 =
1

3
A2

5

h4 = 4A3A4 − 2S−1
1 A2A5 −A2A3 −

1

2
S−1
1 A4A5

h5 =
5

2
S
−1

1
A3A5 − 2S−1

1 A3A4 −A2A4

h6 =
2g

12
A

2

4
− 17

12
A

2

3
− S1A

2
4

h7 = −2A3A5

h8 = −A4A5

h9 = A2A5 + 2S1A3A4 − S−1
1 A2A5

h10 = −5

3
A3A4

h11 = −(S−1
1 A2A4 −

1

4
S−1
1 A4A5)

h12 = S−1
1 A2A4 −

1

2
S−1
1 A3A4 + S−1

1 A4A5

h13 = S−1
1 A3A5 − S−1

1 A4A5

h14 = −S−1
1 A4A5

We found ∪1(η) which has the form.
∪1(η) = ∪1p(η) + ∪1c(η)
By using boundary conditions, we found the constants B1, B2, B3, B4, where

B1 = −B1cos(2S1)−
1

8S
(
F

AE
)[ A2

2 + 6S A2
4–2S A2

3 − S A2
5 +

1

3
A2

5)− 4A3A5cos(2S1)

+ (S−1
1 A3A4 −A2A4)sin(2S1) + [

26

6
A2

4 −
17

6
A2

3 − 2S1A
2
4]sin

2(2S1)−
1

2S1
A4A5

sin(2S1)cos
2(2S1) + (2S−1

1 A2A4 − S−1
1 A3A4 + 2S−1

1 A4A5)cos(2S1)sin
2(2S1).

B2 = −B4 sin (2S1)−
1

8S

F

AE

[
−4

3
A2A5 −

(
8A3A4 − 4S−1

1 A2A4− 2A2A3 − S−1
1 A2A5 cos (2S1) +A4A5 sin (2S1)− (A2A5 + 4S1A3A4+ 5−1A2A5

)
sin2 (2S1) +

10

3
A3A4 sin (2S1) cos (2S1) +

(
2S−1

1 A3A4+ 4S−1A4A5 + 2S−1A3A5

)
sin (2S1) cos

2 (2S1)+ 2S−1
1 A4A5 sin

2 (2S1) cos (2S1)
]



M. M. Aziz / Eur. J. Pure Appl. Math, 16 (2) (2023), 983-996 989

B3 = B4 tan (2S1) +
1

16S2
.
F

AE
.

1

cos (2S1)

B4 =
1

32
S−5
1

1

1− 2 sin2 (2S1)

F

AE
(−32S1A3A4 − 32S1A4A5 − 12S1A2A5 − 12S1A2A3 + 56

S1A3A5 + 4S1A3A5 −8S3
1A2A4 − 8S3

1A4A5

)
+

(
2S1A

2
2 − 60S3

1A
2
4 +

56

3
S3
1A

2
3 + 16S2

1A
2
4

−4S−1
1 A2

5 + 8S3
1A2A5 + 16S2

1A3A4+ 8S1A
2
5 − 16S1A2A5

)
sin (2S1) (16S1A2A5+

136

3
S2
1A3A4 − 12A2A5 + 8A2

5

)
cos (2S1) + (64S1A3A4 − 48 S1A4A5 − 16S1A2A5 + 24S1

A2A3 + 72S1A3A5 + 16S1A3A5) sin
2 (2S1) + (76S1A4A5 − 126S1A3A4− 140S1A2A4 + 44

A3A5 − 44A4A5 + 6S1A3A5 + 64S3
1A3A4 − 16S3

1A2A3

)
sin (2S1) cos (2S1)−

(
232

3
S3
1A

2
4

− 136

3
S3
1A

2
3 − 96S2

1A
2
4 + 48S1A3A5 sin

2 (2S1) cos
2 (2S1) + (144S1A2A4 − 72S1A3A4+

144S1A4A5+ 36A4A5 − 36A3A5) sin
3 (2S1) cos (2S1)

]
.

Thus

∪ (η) = P1 + P2η + P3η
2 + P4η

3 + P5η
4 + P6 cos (2S1η) + P7 sin (2S1η) + P8η cos (2S1η)

+ P9η sin (2S1η) + P10 sin
2 (2S1η) + P11η

2 cos (2S1η) + P12η
2 sin (2S1η) + P13η sin

2 (2S1η)

+ P14 cos (2S1η) sin (2S1η) + P15 sin (2S1η) cos
2 (2S1η) + P16 sin

2 (2S1η) cos (2S1η)

+ P17η sin (2S1η) cos
2 (2S1η) + P18η sin

2 (2S1η) cos (2S1η)] .
(26)

Where

P1 = A1 +DB1, P2 = A2 +DB2, P3 = A5 +
1

8
S−1
1 kA2

2 +
3

4
kA2

4 − 1
4kA

2
3 −

1

4
S−2
1 A2

5

P4 =
1

6
S−1
1 kA2A5, P5 =

1

2
S−1
1 kA2A

2
5, P6 = A3 +DB3, P7 = A4 +DB4,

P8 = S−1
1 kA3A4 −

1

2
S−3
1 kA2A5 − 1

4S
−1
1 kA2A3 −

1

8
S−3
1 kA3A4,

P9 =
5

8
S−1
1 kA3A5 − 1

2S
−1
1 kA3A4 −

1

4
S−1
1 kA2A4,

P10 =
29

48
S−1
1 kA2

4 −
17

48
S−1
1 kA2

3 −
1

4
S−1
1 kA2

5, P11 =
1

2
S−1
1 kA3A5, P12 =

1

4
S−1
1 kA4A5,

P13 =
1

4
S−1
1 kA2A5 +

1

2
S−1
1 kA3A4 − 1

4
S−2
1 kA2A5, P14 = − 5

12S
−1
1 kA3A4,

P15 = −1

4
S−3
1 kA2A3−

1

16
S−2
1 kA4A5, P16 =

1

4
S−3
1 kA2A4−

1

8
S−2
1 kA3A4+

1

4
S−3
1 kA4A5,

P17 =
1
4S

−2
1 kA3A5 −

1

4
S−2
1 kA3A4, P18 = −1

4
S−1
1 kA4A5,K =

FD

AE



M. M. Aziz / Eur. J. Pure Appl. Math, 16 (2) (2023), 983-996 990

4. Stability analysis

In the standard stability, we assume that [10]

∪( η, t∗) = ∪s(η) + ∪d(η, t
∗). (27)

θ( η, t∗) = θs(η) + θd(η, t
∗). (28)

Where ∪d are small perturbations or departures from equilibrium, and Us , θs are the
velocity and temperature in steady state distribution. Substitution of the form (27), (28)
into (9), (10), we get the following equations.

∂θd
∂t∗

= A
∂2θd
∂η2

+ 4Cθd + 2D
∂∪
∂η

• d∪d

dη
. (29)

∂∪d

∂t∗
= E

∂2∪d

∂η2
− Fθd. (30)

d∪s

dη
which is itself

d∪
dη

which we get it from the derivative of equation (??)

The solutions of the system

θd (η, t∗) = C1e
αt∗+iβη. (31)

∪d = C2e
αt∗+iβη. (32)

Where C1, C2 are constant which represent the amplitudes, and β is a dimensionless
wave number, and ∝ is complex wave speed such that α =∝1 +i ∝2 is positive or negative
implies growth or decay of the disturbance.

By taking the derivative of relation (31) and (32) with respect to η, t∗ and then
substitute the result into equation (29) and (30), we get(

α+Aβ2 − 4C
)
c1 − 2iDβ dU

dη c2 = 0

, F c1 + (α+ Eβ)c2 = 0.

}
(33)

and thus, we transform the differential equation to an algebra equation and take the
determent and we equal to zero∣∣∣∣ α+Aβ2 − 4C −2iDβ dU

dη

, F α+ β2

∣∣∣∣ = 0. (34)

which gives

α2 + Lα+M +Ni = 0. (35)

L = Aβ2 − 4C + Eβ2

M = 4Eβ2–4ECβ2



M. M. Aziz / Eur. J. Pure Appl. Math, 16 (2) (2023), 983-996 991

N = 2Dβ
d∪
dη

Through the solution of equation (35), we get the relation between R and β for different
values of Gr, 100, 300, 400 and for different locations of η, 1, 0.5, 0, -0.5, -1 as given in
table (1),(2),(3) it was found that the location has a small effect on the stability of the
system and then if we take one location such as η=0 the system is then become stable for
α1 < 0 and unstable for α1 > 0.

The stable and unstable are shown in Figures (1), (2), (3), (4) from Table (1), (2), and
(3).

Stability curves for different values at Bo = 1, P r = 0.7, Ec = 0.1, W = 0.15, F r =
1000, Gr = 100

Table 1: Values of Neutral Stability Curve at Bo =1, Pr=0.7, Ec=0.1, W=0.15, Fr =1000, Gr=100

R/η -1 -0.5 0 0.5 1

100 2.05 2.06 2.06 2.06 2.06

200 3.08 3.08 3.09 3.08 3.10

300 3.64 3.64 3.64 3.67 3.67

400 4.08 4.12 4.08 4.12 4.08

500 6.10 6.11 6.09 6.08 6.07

600 5.74 5.68 5.72 5.68 5.74

700 5.41 5.43 5.42 5.41 5.42

800 6.11 6.15 6.15 6.15 6.11

900 9.18 9.18 9.18 9.18 9.13

1000 8.94 8.94 9.01 8.94 8.94



M. M. Aziz / Eur. J. Pure Appl. Math, 16 (2) (2023), 983-996 992

Figure 1: Neutral Stability Curve for Bo = 1, Pr = 0.7, Ec = 0.1, W = 0.15, Fr = 1000, Gr = 100

Table 2: Values of Neutral Stability Curve at Bo =1, Pr=0.7, Ec=0.1, W=0.15, Fr =1000, Gr=300

R/η -1 -0.5 0 0.5 1

100 2.31 2.28 2.13 2.46 2.48

200 2.92 2.93 2.92 2.94 2.94

300 3.55 3.55 3.5 3.57 3.57

400 4.11 4.10 4.12 4.11 4.10

500 5.14 5.17 5.13 5.14 5.13

600 5.19 5.19 5.20 5.19 5.19

700 5.44 5.41 5.44 5.41 5.44

800 5.86 5.86 5.86 5.86 5.86

900 7.44 7.43 7.47 7.44 7.44

1000 7.48 7.48 7.45 7.48 7.48



M. M. Aziz / Eur. J. Pure Appl. Math, 16 (2) (2023), 983-996 993

Figure 2: Neutral Stability Curve at Bo =1, Pr=0.7, Ec=0.1, W=0.15, Fr =1000, Gr=300

Table 3: Values of Neutral Stability Curve at Bo =1, Pr=0.7, Ec=0.1, W=0.15, Fr =1000, Gr=400

R/η -1 -0.5 0 0.5 1

100 2.23 2.21 2.19 2.22 2.20

200 2.91 2.91 2.90 2.95 2.98

300 3.57 3.57 3.60 3.55 3.56

400 4.12 4.08 4.09 4.10 4.13

500 5.00 4.98 5.02 4.98 5.00

600 5.12 5.10 5.15 5.10 5.13

700 5.41 5.40 5.45 5.40 5.44

800 5.82 5.85 5.82 5.82 5.85

900 7.14 7.14 7.14 7.14 7.14

1000 7.23 7.17 7.22 7.23 7.20



M. M. Aziz / Eur. J. Pure Appl. Math, 16 (2) (2023), 983-996 994

Figure 3: Neutral Stability Curve at Bo =1, Pr=0.7, Ec=0.1, W=0.15, Fr =1000, Gr=400

Figure 4: Stability curves for different values of Gershoff number at B0=1, Pr=0.7, Ec=0.1, W=0.15,Fr=1000

5. Conclusion

The stability of the mathematical model of disturbed fluid flow between two parallel
plates has been investigated.

Analysis of the results showed the effect of buoyancy force on heat distribution and
velocity, and it turns out that when the Gershoff number is increased, the speed decreases.



REFERENCES 995

In the figures (1), (2), (3), and (4), the boundaries between stable and unstable states
were examined and found for a perturbed system and for different values of the Gershoff
numbers through the Neutral stability curve.

It has been shown through this, that the stability of the model increases with the
increase in the Gershoff numbers (Gr), where the gravitational forces resist the viscous
forces towards stability.

The results indicated that the effects of buoyancy forces have significant contribution
on the basic variable of the model.

References

[1] H I Abu-Mulaweh, B F Armaly, and T S Chen. Measurements in buoyancy-
opposing laminar flow over a vertical backward-facing step. Journal of Heat Transfer,
116(1):247–250, 1994.

[2] Maysoon M Aziz. Applications of mathematics: Dispersion of soluble matter in
solvent flowing through a tube under a steady pressure gradient. Periodicals of En-
gineering and Natural Sciences (PEN), 9(3):302–308, 2021.

[3] Maysoon M Aziz. Applied mathematical model for heat transfer. Periodicals of
Engineering and Natural Sciences (PEN), 9(2):940–945, 2021.

[4] MM Aziz. Stability analysis of mathematical model. International Journal of Science
and Research (IJSR), 7:147–148, 2018.

[5] MM Aziz and DM Merie. Stability and chaos with mathematical control of 4-d dy-
namical system. Indonesian Journal of Electrical Engineering and Computer Science,
20(3):1242–1251, 2020.

[6] D Joseff Daniel. Stability of fluid motions 1 berlin heidelberg, new york. 1979.

[7] L Elsgolts. Differential equations and Calculus of variations. 2016.

[8] H Grab Muller. Singular perturbation techniques applied to integro-differential equa-
tions. Great Britain at Biddles of Guild Ford, 1978.

[9] OM Jihad and MM Aziz. Stability & chaos tests of 2d discrete time dynamical system
with hidden attractors. Open Access Library Journal, 8(6):1–11, 2021.

[10] JD Logan. Applied mathematics. John Wiley & Sons, 2013.

[11] S.V. Maulucci. Spatial type instability wave in the boundary layer at reynolds num-
bers. Journal Fluid Mechanics Soviet Research, 11(6):25–32, 1989.

[12] Y. Mori. Buoyancy effects in forced laminar convection flow over a horizontal flat
plate. J. Heat Transfer, 83(1):479, 1961.



REFERENCES 996

[13] Konstantin G Shavarts and Abd alaziz. Effect of rotation on stability of a defec-
tive flow in horizontal liquid layer with free upper boundary. Journal of Physics:
Conference Series, 216(1):012002, 2010.