1. Introduction
Applications of the magnetic field in MHD pump, heat exchangers, etc., has attracted many researchers’ attention [
1]. Hartmann (1998) formulated the electromagnetic viscous force [
2]. However, unsteady MHD nanofluid compressed flow between parallel plates, observed by Dogonchi and Ganji [
3,
4], studied a similar effect with solar radiation. Ref. [
5] recently analyzed unsteady magnetohydrodynamics (MHD) nanofluid flow using natural conviction with entropy generation. Additionally, entropy formation in a non-Newtonian microflow regulated by conjugate heat transport was presented in [
6].
Nayak and Redd et al. [
7] studied the viscous dissipation effect and Maxwell liquid effect on MHD flow using thermal radiation. Refs. [
8,
9,
10] depict some of the work on MHD flow through the variable surface topology.
The industrial applications of heat and mass transport using chemical reactions are substantial for many scientists. Route et al. [
11] observed the chemical reaction, mass, and heat sources on MHD flow in an asymmetrical channel. Refs. [
12,
13] studied the effect of chemical reaction on convection flow of viscous nanofluid using a magnetic field and over a permeable sheet to solve the governing equation. In [
14,
15,
16], fluid flows with effects of chemical reactions were studied using Jeffery liquid flow due to radially stretching sheets. Refs. [
17,
18,
19] are some examples in which the chemical reaction effect was coupled with nanofluid, and resulting consequences were intimated on the industrial level. In [
20], heat and mass transfer of upper-convected Maxwell nanofluid with the magnetic field’s effects was studied past an inclined stretching surface. A shooting method consisting of the fourth-order Runge–Kutta method (RK4) was utilized to find the solution to ordinary differential equations (ODEs) obtained from PDEs.
One of the applications of Newtonian fluids is their use in lubrication. In lubrication, the fluid does not change viscosity as a function of shear, but it depends on the temperature. Motor oil is one example of this process. The oil does not become too thin at higher engine operating temperatures. Moreover, its important characteristic is that the viscosity does not change if this oil is shared between engine elements [
21].
Various authors studied the flow in a duct for heat and mass transfer effects. The flow in a rectangular cross-section was given in [
22], and the flow in a circular cross-section was given in [
23,
24]. The problem of flow in a rectangular duct with a constant streamwise pressure gradient was given in [
25]. In [
26], the authors focused on the entropy generation for viscoelastic fluid flow through a parallel plate microchannel under pressure gradient, interfacial slip, and conjugate heat transfer. The effect of viscous dissipation on convective heat transfer and entropy formation for Poiseuille flow in an asymmetrically heated slit microchannel was studied in [
27].
An extensive literature study revealed that the amount of education based on blood flow in a porous saturated stenotic artery has dwindled. Dash et al. [
28] studied blood flow in the tube via a homogenous porous medium. The central theme behind this study was that the distribution of fatty cholesterol and artery-clogging blood clots in the lumen of a coronary artery is equivalent to a fictitious porous medium in some pathological situations. This idea was manipulated by [
29]. The blood was considered a Newtonian flood and found its flow through a porous medium under magnetic field influence.
The present attempt proposes a new class of explicit numerical schemes for the temporal discretization of parabolic and first-order hyperbolic differential equations. These schemes are constructed at different time levels, and construction was made through the Taylor series expansions. A Von Neumann stability criterion was used to find the stability condition of the second-order scheme for a parabolic linear equation with diffusion term only. However, from a computation point of view, the scheme takes a small temporal step size to produce an accurate solution. In order to reduce the computational time, a modification is suggested. This suggested modification was used when the Du Fort–Frankel scheme was constructed from the second-order unconditionally unstable scheme. The modification in diffusion terms made the Du Fort–Frankel method unconditionally stable. A modification that made the Du Fort–Frankel method unconditionally stable is the average of two neighboring time levels to approximate the time at level “n”. Therefore, here, the same modification was implemented in the proposed scheme on the second-order central difference formula for diffusion term. Additionally, using this modification, the scheme permits a little larger time step size. Despite this, the scheme loses some of its temporal order as a result of this adjustment, and the application of a first-order approximation causes this loss of order. For obtaining second-order accuracy, one can use a third-order approximation to approximate the time at level “n,” so the order of the modified scheme will not be affected if one can find the solution using the modified approach in an implicit manner. The work gave the third-order modification for the Du Fort–Frankel method [
30]. Since the first-order approximation in the scheme permitted a little larger value of the diffusion number for the present flow problem, only the second-order proposed and the first-order modified scheme was used to produce the final results for the displayed graphs.
Since Adams–Bashforth’s methods are the class of methods used to solve the ordinary differential equations, in [
30], a class of methods was given for temporal discretization of 1D hyperbolic partial differential equations. This class of methods [
30] can give second to fifth-order accuracy for hyperbolic partial differential equations. In [
31,
32], a Fully Explicit and Compact Numerical Scheme of the Fourth Order for the Heat Transfer in Boundary Layer Flow was proposed and discussed.
The second-order Adams–Bashforth method for the temporal discretization of PDEs given in [
28] can be employed for the following parabolic partial differential equation
By using the existing numerical scheme, the discretized equation is given by
where
.
This method (2) is second-order accurate in space and time, which finds a solution at “” time level by using the previous two time levels, “” and “”. This existing explicit method (2) has an advantage of accuracy over an existing Du Fort–Frankel method.
2. Proposed Numerical Scheme
This contribution is consisted on a similar kind of explicit class of schemes as given in Equation (2). For developing a scheme with a better stability region than existing scheme (2), the following modification is considered in the scheme given in Equation (2)
and so the proposed scheme for discretizing Equation (1) is given as
In the initial stages of the this contribution, an explicit second-order numerical scheme of temporal discretization is proposed for time-dependent partial differential equations. Standard difference formulae or any other spatial discretization schemes can be used to discretize space variable(s). Additionally, stability and consistency for scalar equation and stability for the considered system of equations for heat and mass transfer of reactive chemical fluid flow in a rectangular duct are given.
The general difference form of discretized version of Equation (1) is given by
The Taylor series expansions for
,
and
terms are given in Equation (3) as
Substituting (6)–(8) into Equation (5) gives
On equating the coefficients of
,
and
on both sides of Equation (9) gives
Solving Equation (10) gives
Therefore, for second-order in time, an explicit scheme is given by
Similarly, a scheme with temporal order of three can be constructed, which is given by
where
The mentioned scheme (13) consists of two free parameters, and schemes with orders four and five can be constructed similarly with three and four free parameters, respectively. Therefore, the class of Adams–Bashforth schemes for time-dependent PDEs is just a special case for the proposed schemes. Moreover, the second-order scheme with gives the same stability region as produced by the second-order Runge–Kutta method, which is more than the standard second-order Adams–Bashforth method in PDEs.
For Equation (1), the scheme given in Equation (12) with
can be expressed by
If the second-order central difference formula for spatial discretization in Equation (14) is employed, then the resulting scheme is given as
In the next section, the stability of the numerical scheme is given using Von Neumann stability criteria for the following equation:
3. Scalar Stability
Discretize Equation (16) by using the constructed scheme in the previous section with second-order central spatial discretization. The discretize equation is given by
In order to employ the Von Neumann stability criteria for Equation (17), consider the following transformations
where
. Substituting transformations (18) into Equation (17) gives
Dividing both sides of Equation (19) by
gives
Equation (20) can be expressed as
Let ] and ]
Now, Equation (21) can be expressed as
One more equation can be constructed as
System of Equations (22) and (23) can be expressed in the following single matrix–vector of the form
The stability conditions can be imposed on the eigenvalues of the matrix given in (24) as
where
The most restrictive condition can be made by checking the expression
to be positive or negative or zero. Here, two cases were considered for extreme values of
and
. For the first case, when
, then the expression is
given by
Since the coefficient of is positive, its graph, which is a parabola, opens upwards, and by using optimization, it can also be proved that the vertex of the parabola (26) lies above to -axis; due to these two mentioned facts, the expression (26) gives a positive value for every value of the diffusion number for the considered case. For the second case, when , the expression becomes zero, and the remaining two cases can be discussed similarly to the first case. Thus the eigenvalues are real, and the stability conditions can be found in inequalities and . It is to be noted that the main contribution to stability condition can be found from two inequalities, and and these inequalities can be implied from the inequalities and .
Now, by using the inequality
. The following inequality can be constructed:
, which gives
where
.
Similarly, other two cases, when
and
or
and
, can be discussed. The case when
and
is straightforward, and the inequality
becomes true for every value of the diffusion number
. The remaining inequality to be checked is
, and from this inequality, one of the conditions that can be imposed on stability is that
which holds because
for
. Moreover, the inequality
implies
which gives
or
for four different cases depending upon the extreme values of
and
, and these inequalities are true for any choice of diffusion number
. The main contributing cases for stability were discussed, and the stability condition for the considered PDE is
.
4. Consistency of the Scheme
The consistency of the scheme is checked by employing Taylor series expansions for the terms in the scheme. Consider the scheme (17) for discretizing Equation (16) and the Taylor series expansions for the terms in Equation (17) are given as
Similarly, Taylor series expansions for the term
and
can be constructed. Now, some of the terms in the scheme (17) can be expressed by
Adding Equations (31) and (32) gives
Now, the Taylor series expansion for the terms
and
are given by
Similarly, Taylor series expansion for the terms
and
can be constructed. By using (34)–(36), some terms in the scheme (17) for time level “
” can be simplified as
Similarly, using the Taylor series, the terms
and
can be expanded and simplified as
By adding Equations (37) and (38), the resulting equation is given as
By using Equations (17), (33) and (39), and Taylor series expansions for the terms
and
, the resulting equation is given by
Simplifying the terms in Equation (40) yields
By using the consistency criteria, and applying the limits the resulting equation becomes the original parabolic partial differential Equation (16) evaluated at the grid point and at time level “”. Thus the scheme is consistent. Now using Equation (16) in Equation (41), the leading error term in the scheme is given by
5. Problem Formulation in Fluid Flow
Consider a fully developed laminar, incompressible flow in a rectangular duct with a streamwise constant pressure gradient. The cross-sectional area of a duct is shown in
Figure 1. Let
be a velocity component along
-axis,
is the temperature of the fluid, and let
is the concentration poured into the fluid. Two boundaries of the duct are kept at a hot temperature
, and an ambient temperature is used for the left and right walls of the rectangular cross-sectional region. The flow is generated due to the constant pressure gradient exerted on the fluid. Dirichlet-type boundary conditions for velocity, temperature, and concentration are used for all four sides of a duct. The concentration is poured into the specified duct location from its bottom wall.
The governing equations of the considered flow problem using the effect of the chemical reaction can be expressed as
subject to the boundary conditions
where
is the density of the fluid,
is the kinematic viscosity,
is the thermal diffusivity,
is the Brownian diffusion coefficient, and
is the reaction rate parameter.
To make Equations (42)–(47) dimensionless, the following transformations are employed:
where
is the average velocity of the fluid. By incorporating transformations (48) for Equations (42)–(47), the following set of partial differential equations are obtained
Subject to the dimensionless boundary conditions
where
is the Reynolds number,
and
are the thermal and solutal Péclet numbers, respectively, and
is the reaction rate parameter, which is defined by
To solve Equations (49)–(54), a proposed numerical scheme for temporal discretization is employed with second-order central differencing formulas are adopted for spatial discretization.
Using a proposed second-order scheme for temporal discretization with standard central difference formulas for the spatial discretization to Equations (49)–(51) yields the following equation
where
and similarly, the remaining quantities in Equations (56) and (58) can be defined. The summary for finding solutions of Equations (49)–(54) using Matlab software is given in the form of Algorithm 1.
Algorithm 1 Summary of the Matlab code |
1. Initialization | Begin The Code With Input Data. |
2. Temporal For Loop | Begin The “For Loop” for Temporal Variable First and Give Initial Conditions. |
3. Spatial For Loops | Begin “For Loops” and for Spatial Variables and Give Boundary Conditions Using if-else Statements. |
4. Processing | Find the Solution for at the Grid Points and at Time Level |
5. End | End the Conditional Statements if-else and also end the Spatial and Temporal Loops. |
6. Plots | Draw the Desire Plots Using Matlab Command Window. |
6. Vector Stability
In order to check the stability conditions for Equations (49)–(51), a Von Neumann stability criteria is adopted again for this case. Therefore, for doing so, first, convert the Equations (49)–(51) into a single following matrix–vector equation
where
.
Employing the proposed numerical scheme for temporal discretization with the second-order central scheme for spatial discretization for Equation (59) gives:
where
.
By inserting the following transformations
into Equation (60), it gives the following equation after dividing both sides of the resulting equation by
Equation (61) can be expressed as
Let
where
and
can be expressed by
with
,
,
and
is an identity matrix.
By substituting
and
into Equation (53), the resulting equation can be expressed by
Consider now an equation of the form
System of Equations (63) and (64) can be expressed in the form of a matrix–vector equation given as
The stability conditions can be imposed on the eigenvalues of the matrix given in Equation (65) as
where
and
and
.
In the expressions of and , denotes the maximum eigenvalues of the matrix .
Theorem 1. A presently proposed scheme tends to converge for the dimensionless conservative form of considered PDEs if the error made by the scheme at the first time level is bounded and the following condition holds
where
,
and
Proof. Non-dimensional form of Equations (43)–(45) can be expressed in the following matrix–vector form
By discretizing Equation (67) using the proposed scheme with second-order central difference formulas, the following equation is obtained:
Let the exact discretized scheme for discretizing Equation (67) be given as
By subtracting Equation (69) from Equation (68), denoting
,
, the resulting Equation is given by
By applying norms on both sides of Equation (70), the following inequality is obtained:
By using inequality
at time level
and applying
, the resulting inequality can be expressed as
Now, if
, then the Inequality (72) can be simplified as
Let
, then Inequality (73) can be expressed by
For
in Inequality (74), the following inequality is obtained
Since the initial condition is exact, so
therefore, Inequality (75) becomes
By using the hypothesis, the following inequality can be constructed
For
in Inequality (74), it gives
For
in Inequality (74), the following inequality is obtained
For
in Inequality (74), it yields
Similarly, this process continues, and for even
,
When
is odd, an inequality becomes
Since the series is a geometric series that converges for and similarly if , then it can be proved that and therefore theorem is proved. □
7. Results and Discussion
The set of governing non-dimensional partial differential equations was solved by employing the proposed second-order in time with second-order and combined second and fourth-order spatial discretization. The proposed second-order scheme was proved to be conditionally stable for the case of scalar parabolic partial differential equations. In addition to this, a modification was suggested that can be coupled with the proposed scheme, which permits a little larger stability region. However, it has the disadvantage of reduction in the order of accuracy. Since the modification is first-order accurate, the order of the resultant scheme is reduced, and the proposed scheme (4) and Du Fort–Frankel scheme become first-order accurate. However, since Richardson’s scheme was unconditionally unstable, the modification makes it unconditionally stable, and here, the same modification was used to enhance the stability region. The modification overcomes the computational time issue for a second-order scheme. Because a scheme with a very small temporal step size may take more time to converge compared to a scheme that may take a short time to converge, the disadvantage is dropping a temporal order of a scheme. However, one can use a third-order approximation instead of Du Fort–Frankel’s first-order approximation for the diffusion term. The third-order approximation [
29], as mentioned in the introduction of this work, can also be considered, and so in this way, the order of the scheme is not reduced in some or most of the cases. Moreover, the Du Fort–Frankel method consists of the time derivative on
time level. For this case, the time derivatives were used for two consecutive time levels,
and
, but the modification is suggested for the
time level and with a spatial discretization of a diffusion term(s). The stability of the system of equations was also given.
In drawing
Figure 2,
Figure 3,
Figure 4,
Figure 5,
Figure 6,
Figure 7,
Figure 8,
Figure 9,
Figure 10,
Figure 11 and
Figure 12,
were used.
Figure 1 shows the geometry of the cross-section of a rectangular duct. Dirichlet-type boundary conditions were implemented on all four sides of the rectangular cross-sectional region of the duct. The
coordinate is the streamwise coordinate of the flow, and the flow was generated due to different constant pressure gradients exerted on the fluid from the inlet. A small part of the bottom boundary was used for concentration. Therefore, in non-dimensional form, the concentration has a maximum value at the particular part of the boundary, and at the remaining boundaries, it is zero.
Figure 2 shows the flow velocity at different pressure gradients. Since the fully developed flow was considered, the parabolic behavior of the velocity can be visualized in
Figure 2.
Figure 3 shows velocity contours at different times. Since no-slip boundary conditions were imposed at the walls, the same color representing these conditions is shown in
Figure 3.
Figure 3 shows that more layers of fluid are disturbing with time.
Figure 4 shows the flow temperature with the variation in thermal Péclet number.
By decaying thermal Péclet number, the temperature of the flow de-escalates. This results from lower thermal conductivity due to a decrease in thermal diffusivity. The decay in thermal diffusivity is the consequence of the direct relationship between thermal diffusivity and thermal Péclet number.
Figure 5 shows the contours of the temperature of the flow.
Since the temperature has its maximum values at the top and bottom edges of the rectangular region, the temperature of the fluid near these walls is maximum. Then, due to transferring heat into adjacent layers of the fluid, the temperature of those layers also increases. The increase in thermal Péclet number produces lower heat transfer, which can be seen in
Figure 5. The decrease in the heat transfer by increment in thermal Péclet number is the consequence of decreasing thermal conductivity discussed earlier in this contribution.
Figure 6 shows the contours of the temperature of the flow at different value times.
Figure 6 deliberates that more layers of fluid become hotter with time. The temperature has its minimum values at the left and right walls of the rectangular cross-sectional region.
Figure 7 deliberates the temperature of the flow at different times and different
coordinates.
Due to transferring of heat, the temperature approaches its maximum values at both locations of chosen
coordinate by time. Effect of thermal Péclet number on the temperature of the flow is shown in
Figure 8.
These two-dimensional graphs show that temperature is decreasing by growing values of thermal Péclet number. The effect of solutal Péclet number on concentration is shown in
Figure 9.
The mass transfer is better at smaller values of the solutal Péclet number. This is the consequence of decreasing mass diffusivity by increasing solutal Péclet number because both are inversely proportional to each other.
Figure 10 shows the contours of concentration with time. The concentration is maximum at the specified location at
-axis due to imposed boundary conditions. The transfer of mass from this location to the adjacent layers of fluid flow can be seen in
Figure 10.
The effect of the reaction rate parameter on concentration can be seen in
Figure 11. The dual behavior of concentration can be seen in
Figure 11.
For verification of the code for the proposed scheme for a time-dependent partial differential equation, the following example is considered:
subject to the boundary conditions
and using zero as the initial condition. The problem is solved using the proposed scheme in time and fourth-order spatial discretization except at those grid points where fourth-order central spatial discretization cannot be employed. The exact solution of the problem (83) and (84) is available. The comparison of the proposed scheme with the exact solution is given in
Figure 12.
Figure 12 shows the absolute error produced by the proposed scheme, which is obtained by finding the absolute value for the difference between numerical and exact solutions.
Figure 13 demonstrates the usefulness of the proposed scheme (4), which is the combination of the existing second-order scheme (2) and a modification given in the existing DuFort method (3). The proposed scheme produces a better stability region, which is evidence of the proposed strategy in
Figure 13. The existing second-order scheme does not converge on chosen temporal step size. However, the proposed scheme (4) converges on the same temporal step size, but the temporal order of the proposed scheme (4) is reduced to one due to first-order modification (3). In the caption of
Figure 13,
,
, and
denote the number of time levels, grid points, and final time, respectively.