Next Article in Journal
On an Intuitionistic Fuzzy Form of the Goguen’s Implication
Previous Article in Journal
A Note on Some Definite Integrals of Arthur Erdélyi and George Watson
Previous Article in Special Issue
Calculation of Two Types of Quaternion Step Derivatives of Elementary Functions
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Mathematical Analysis of Oxygen Uptake Rate in Continuous Process under Caputo Derivative

by
Rubayyi T. Alqahtani
1,
Abdullahi Yusuf
2,3,* and
Ravi P. Agarwal
4
1
Department of Mathematics and Statistics, College of Science Al Imam Mohammad Ibn Saud Islamic University (IMSIU), Riyadh 11566, Saudi Arabia
2
Department of Computer Engineering, Biruni University, 34010 Istanbul, Turkey
3
Department of Mathematics, Federal University Dutse, Jigawa 7156, Nigeria
4
Department of Mathematics, Texas A M University-Kingsville, Kingsville, TX 78363, USA
*
Author to whom correspondence should be addressed.
Mathematics 2021, 9(6), 675; https://doi.org/10.3390/math9060675
Submission received: 25 February 2021 / Revised: 11 March 2021 / Accepted: 12 March 2021 / Published: 22 March 2021
(This article belongs to the Special Issue Nonlinear Equations: Theory, Methods, and Applications)

Abstract

:
In this paper, the wastewater treatment model is investigated by means of one of the most robust fractional derivatives, namely, the Caputo fractional derivative. The growth rate is assumed to obey the Contois model, which is often used to model the growth of biomass in wastewaters. The characteristics of the model under consideration are derived and evaluated, such as equilibrium, stability analysis, and steady-state solutions. Further, important characteristics of the fractional wastewater model allow us to understand the dynamics of the model in detail. To this end, we discuss several important analyses of the fractional variant of the model under consideration. To observe the efficiency of the non-local fractional differential operator of Caputo over its counter-classical version, we perform numerical simulations.

1. Introduction

Mathematical models as routes of interpreting real-life situations are included in the literature. These models display a vital update in the quantification and assessment of real-life problems and preventive measures [1,2,3]. Mathematical modeling has, in many ways, proven to be a very versatile and effective way of studying the dynamics of many situations. To build and test control measures that are persuasive, statistical analysis and numerical simulations may be used. It is also widely accepted that by using mathematical models, the prediction to control and solve real-life problems can be reached. A relevant argument here is that it is important to obtain acceptable criteria for the particular situation under consideration and to use these criteria to determine the effects of potential control measures. The dilemma, then, is how, from an optimal point of view, to implement such measures. Several noteworthy attempts have recently been made to establish results through mathematical modelling [4,5,6,7,8].
One of the key areas to see the application of mathematical models is the field of biochemical engineering, where oxygen is a key substrate in aerobic bioprocesses that is used for growth, maintenance, and other metabolic routes, including product synthesis [9]. Oxygen is continuously supplied via the gas phase due to its low solubility in broths, which are typically aqueous solutions, and thus, knowledge of the oxygen transfer rate is needed for the design and scale-up of the bioreactor, which can be addressed by a mathematical modelling approach. Both the transition of oxygen mass from the phase of liquid to gas and its microorganism intake are of decisive importance, as the transfer of oxygen is the control stage for microbial growth in many processes and may influence the evolution of bioprocesses. Early in the history of biochemical engineering, this fact was established and reviewed [10,11,12,13,14]. The oxygen uptake rate is one of the essential physiological features of culture growth, and has been used to optimize the fermentation method. The real oxygen uptake rate is usually estimated from the experimentally defined oxygen uptake rate [15,16].
We consider a wastewater treatment model with the Contois model for the growth rate. As the main feeds into a tank, the aerator is named [17,18,19], and the usage theory can be established. At this point, by biological oxidation of the substrate, one bacteria population degrades the pollutant. By absorbing oxygen, the reaction produces an aerobic environment. The mixture is transmitted to a settler tank in the second stage. Here, because of gravity, the solid pieces settle and collect at the bottom. Part of the sludge inside the aerator is recycled to promote oxidation. The topic of waste water treatment has also been studied from many points of view, such as complex observation [20,21] and control [22,23], although the anaerobic component is considered in the majority of situations [24,25,26].
However, this model has never been studied using the contemporary fractional operators with the Contois model. Additionally, it has been multifariously proved that natural phenomena can be more effectively enunciated by fractional models than by differential integer-order equations. The fractional calculus has taken on the value and popularity of modeling realistic cases, in particular those with memory effects, because of this advantage [27,28,29,30]. Fractional operators have been used in so many fields where [31,32,33,34,35] is used. We were inspired to study and evaluate a new fractional version of Caputo’s operator model given this value. This is the first time that the Caputo operator has been hired for the model being considered, to the best of our knowledge.

2. Formation of the Model

We present the model to be investigated in this research here. The model with dimensional form is given by:
V d S d t = F ( S 0 S ) V X μ ( S , X ) Y x / s , V d X d t = F β X + V X μ ( S , X ) K d V X .
The oxygen accumulation in the liquid phase can be expressed as
V d C d t = F β C + V K I a ( C 0 C ) V X μ ( S , X ) Y x / O 2 ,
where
-
μ ( S , X ) is the Contois model for the specific growth rate, which is expressed as μ ( S , X ) = μ m S K s X + S ;
-
The rate of oxygen transfer from the gas to the liquid phase is provided by
O T R = K I a ( C 0 C ) ;
-
Y x / O 2 is the overall yield of cells from oxygen; and
-
K I a is the oxygen mass transfer coefficient.
A detailed definition of the model can be found in Nomenclature.

The Non-Dimensional Model

Non-dimensionalization is very important in mathematical modeling. Here, we aim to non-dimensionalize (1) and (2) to study various mathematical properties of the model. Equations (1) and (2) can be written in non-dimensional form by the following equations using non-dimensional variable C * = C S 0 for oxygen concentration, S * = S S 0 for the concentration of substrate, X * = K s X S 0 as the microorganism concentration, and t * = μ m t as time
d S * d t = 1 τ ( 1 S * ) S * X * α * ( S * + X * ) , d X * d t = β X * τ * + S * X * ( S * + X * ) K d * X * , d C * d t = β C * τ * + K I a * ( C 0 * C * ) S * X * γ ( S * + X * ) .
In Equation (4), τ * is the main experimental control parameter. It is assumed for the investigation of the mathematical model that the substrate concentration entering into the bioreactor is positive S 0 > 0 and the rate of change of death is positive K d * > 0 . In addition, it is assumed that 0 < K d * < 1 and
A1:
All the parameters are nonnegative.
A2:
0 μ ( S , X ) 1 , μ ( 0 , X ) = 0 .
In Equation (4),
  • K d * = K d μ m represents the dimensionless death rate,
  • τ * = V μ m F represents the dimensionless residence time,
  • K I a * = K I a μ m represents the dimensionless oxygen mass transfer coefficient, and
  • C 0 * = C 0 S 0 , α * = K s Y x / s , and γ = K s Y x / O 2 are dimensionless parameters.
In the rest of the analysis of Model 4, we drop the star symbols from Equation (4).

3. Analysis of the Model in Classical Sense

In this section, our aim is to analyze the properties of the model in a classical sense.

3.1. Equilibria and Stability Analysis

We consider the number of equilibrium solutions of the model (4). It is obvious that the model has a branch of the washout given by
E 0 = ( S , X , C ) = 1 , 0 , C 0 K Ia τ K Ia τ + β .
One can obtain the model’s non-free steady-state (4) by setting the right side to zero. From the model (4), we have,
S = α τ ( 1 K d ) + α β , X = τ ( K d 1 ) + β α τ ( K d 1 ) α + β K d τ + β .
Substituting (6) into the third equation of (4), we obtain
C = C 0 K Ia τ K Ia τ + β τ ( K d 1 ) + β α γ τ ( K d 1 ) α + β K Ia τ + β .
Next, we focus on analyzing the stability of the equilibrium of (4). The stability of the two forms of the equilibrium, E 0 and E 1 , is studied.
Lemma 1.
The steady-state solution E 0 is locally asymptotically stable when τ < β 1 K d , and is unstable when τ > β 1 K d .
Proof. 
The Jacobian matrix of system (4) at E 0 is given by
J ( E 0 ) = τ 1 α 1 0 0 β τ + 1 K d 0 0 γ 1 [ β τ + K Ia ] .
The eigenvalues of (8) is expressed as
λ i = τ 1 τ ( 1 K d ) β τ K Ia τ + β τ
Note that λ 1 , 3 are negative. After some algebra, we notice λ 2 < 0 when τ < β 1 K d . □
Lemma 2.
The solution E 1 is locally asymptotically stable for a physically meaningful solution.
Proof. 
For the stability of the solution E 1 , it is possible to write the Jacobian of the system as
J ( E 1 ) = J 11 J 12 0 J 21 J 22 0 J 31 J 32 J 33 ,
where,
J 11 = 1 τ + X 2 α S + X 2 , J 12 = S 2 α S + X 2 , J 21 = X 2 S + X 2 , J 22 = X S S + X 2 J 31 = X 2 γ S + X 2 , J 32 = S 2 γ S + X 2 , J 33 = K Ia τ + β τ .
As a simple analysis, J ( E 1 ) characteristics are expressed as
λ 3 + D 1 λ 2 + D 2 λ + D 3 , where D 1 = J 33 + J 22 + J 11 , D 2 = J 12 + J 22 + J 33 J 11 + J 33 J 22 , D 3 = J 33 J 11 J 12 + J 22 .
We use the Rough–Hurtwiz criterion, that is, D 1 > 0 , D 3 > 0 , and D 1 D 2 D 3 > 0 , to illustrate the stability of the solution E 1 . We have
D 1 D 2 D 3 = D 11 J 11 2 + D 22 J 11 + B 33 > 0 . D 11 = J 12 + J 22 + J 33 , D 22 = J 22 J 12 + J 22 2 + 2 J 33 J 22 + J 33 2 , D 33 = J 33 + J 22 J 33 J 22 .
Thus, the condition of the Routh–Huewitz theorem for the stability of E 1 is satisfied. □

3.2. Numerical Simulations

The model’s steady-state solution is evaluated in this section with respect to parameters using the experimental values in Table 1. Figure 1 shows stable curves of the oxygen concentration, the substrate concentration, and the microorganism concentration as functions of residence time for different values of the parameter β . Note that the parameter β represents the configuration of the reactor. β = 1 corresponds to the continuous flow reactor, while β = 0 corresponds to the membrane reactor.
Figure 1a shows that the oxygen concentration decreases as a function of residence time to minimum value, called the minimum residence time. After this minimum value, the oxygen concentration increases as the value of the residence time increases. Note that the minimum residence time for β = 1 and β = 0.1 are 1.30025 and 0.201052, respectively. Thus, the value of parameter β must be taken lower to minimise the oxygen concentration.
In Figure 1b, the substrate concentration decreases as the value of the residence time increases for different values of β . In Figure 1c, the microorganism concentration increases to the maximum value, then decreases as the value of the residence time.
Figure 2 is an unfolding diagram for the minimum residence time in Figure 1a as a function of β and the influent concentration ( S 0 ). We note that the minimum residence time increases as the value of these parameters increases.

4. Analysis of the Model under Caputo Fractional Operator

Here, we aim to extend the model (4) by incorporating the fractional order in the sense of a Caputor fractional derivative. The model (4) in the sense of a Caputo fractional operator is given by
C D ψ S = 1 τ ψ ( 1 S ) S X α ψ ( S + X ) , C D ψ X = β ψ X τ ψ + S X ( S + X ) K d ψ X , C D ψ C = β ψ C τ ψ + K I a ψ ( C 0 ψ C ) S X γ ψ ( S + X ) .

4.1. Results for Existence of Solutions under Caputo

In this chapter, we address the existence and uniqueness of the suggested model’s solution via fixed-point theorems. Let’s reformulate the proposed (13) model in the form below:
C D 0 + ψ S = ζ 1 ( t , S , X , C ) , C D 0 + ψ X = ζ 2 ( t , S , X , C ) , C D 0 + ψ C = ζ 3 ( t , S , X , C ) ,
with
ζ 1 ( t , S , X , C ) = 1 τ ( 1 S ) S X α ( S + X ) , ζ 2 ( t , S , X , C ) = β X τ + S X ( S + X ) K d X , ζ 3 ( t , S , X , C ) = β C τ + K I a ( C 0 C ) S X γ ( S + X ) .
Now, (13) becomes
C D 0 ψ ϖ ( t ) = Λ ( t , ϖ ( t ) ) ; t J = [ 0 , b ] , 0 < α 1 , ϖ ( 0 ) = ϖ 0 0 ,
only if
ϖ ( t ) = ( S , X , C ) T , ϖ ( 0 ) = ( S 0 , X 0 , C 0 ) T , Λ ( t , ϖ ( t ) ) = ( ζ i ( t , S , X , C ) ) T , i = 1 , 2 , 3 ,
( · ) T implies the transpose operation. Now, (16) becomes
ϖ ( t ) = ϖ 0 + J 0 + υ Λ ( t , ϖ ( t ) ) = ϖ 0 + 1 Γ ( υ ) 0 t ( t τ ) υ 1 Λ ( τ , ϖ ( τ ) ) d τ .
Let E = C ( [ 0 , b ] ; R ) be a Banach space for all the functions which are continuous from R [ 0 , b ] with the norm ϖ = sup t J | ϖ ( t ) | .
Theorem 1.
Let the function Λ C ( [ J , R ] ) and maps of the bounded subset of J × R 5 be compact subsets of R . Further, there is a constant L Λ > 0 whereby
( A 1 ) | Λ ( t , ϖ 1 ( t ) ) Λ ( t , ϖ 2 ( t ) ) | L Λ | ϖ 1 ( t ) ϖ 2 ( t ) | ; t J and all ϖ 1 , ϖ 2 C ( [ J , R ] ) . Thus, (18), which is equivalent to (13), has a unique solution whenever Ω L Λ < 1 , and
Ω = b ψ Γ ( ψ + 1 ) .
Proof. 
Consider that P : E E is defined by
( P ϖ ) ( t ) = ϖ 0 + 1 Γ ( ψ ) 0 t ( t τ ) ψ 1 Λ ( τ , ϖ ( τ ) ) d τ .
Now, P is the unique solution of (13) and is well-defined and depicts the fixed-point of P. Consider sup t J Λ ( t , 0 ) = M 1 and κ ϖ 0 + Ω M 1 . Now, it suffices to justify P H κ H κ , and H κ = { ϖ E : ϖ κ } , as convex and closed. Thus, each ϖ H κ , we have
| ( P ϖ ) ( t ) | | ϖ 0 | + 1 Γ ( ψ ) 0 t ( t τ ) ψ 1 | Λ ( τ , ϖ ( τ ) ) | d τ ϖ 0 + 1 Γ ( ψ ) 0 t ( t τ ) ψ 1 | Λ ( τ , ϖ ( τ ) ) Λ ( τ , 0 ) | + | Λ ( τ , 0 ) | d τ ϖ 0 + ( L Λ κ + M 1 ) Γ ( ψ ) 0 t ( t τ ) ψ 1 d τ ϖ 0 + ( L Λ κ + M 1 ) Γ ( ψ + 1 ) b ψ ϖ 0 + Ω ( L Λ κ + M 1 ) κ .
We justify the results. Further, for ϖ 1 , ϖ 2 E , one attains
| ( P ϖ 1 ) ( t ) ( P ϖ 2 ) ( t ) | 1 Γ ( ψ ) 0 t ( t τ ) ψ 1 | Λ ( τ , ϖ 1 ( τ ) ) Λ ( τ , ϖ 2 ( τ ) ) | d τ L Λ Γ ( ψ ) 0 t ( t τ ) ψ 1 | ϖ 1 ( τ ) ϖ 2 ( τ ) | d τ Ω L Λ | ϖ 1 ( t ) ϖ 2 ( t ) | .
This justifies that ( P ϖ 1 ) ( P ϖ 2 ) Ω L Λ ϖ 1 ϖ 2 . Thus, due to Banach contraction, the unique solution for (13) is reached. □
Now, we go for the existence of (13) solutions by Krasnoselskii’s fixed-point justification.
Lemma 3.
Let M be a closed, bounded, and convex subset of a Banach Space E . Let two operators that respect the given relation be P 1 , P 2 .
  • P 1 ϖ 1 + P 2 ϖ 2 M , provided that ϖ 1 , ϖ 2 M ;
  • P 1 is compact and continuous;
  • P 2 is a contraction mapping.
Then, there is u M as far as u = P 1 u + P 2 u .
Theorem 2.
Surmising Λ : J × R 5 R is continuous and holds for the condition ( A 1 ) . Further, let
( A 2 ) | Λ ( t , ϖ ) | ψ ( t ) , f o r a l l ( t , ϖ ) J × R 5 a n d ψ C ( [ 0 , b ] , R + ) .
Thus, Equation (13) has at least one solution whenever
L K ϖ 1 ( t 0 ) ϖ 2 ( t 0 ) < 1 .
Proof. 
Setting sup t J | ψ ( t ) | = ψ and η ϖ 0 + Ω ψ , we consider B η = { ϖ E : ϖ η } . Assume P 1 , P 2 operators on B η are expressed as
( P 1 ϖ ) ( t ) = 1 Γ ( ψ ) 0 t ( t τ ) ψ 1 Λ ( τ , ϖ ( τ ) ) d τ t J ,
and
( P 2 ϖ ) ( t ) = ϖ ( t 0 ) , t J .
Now, each ϖ 1 , ϖ 2 B η , gives
( P 1 ϖ 1 ) ( t ) + ( P 2 ϖ 2 ) ( t ) ϖ 0 + 1 Γ ( ψ ) 0 t ( t τ ) ψ 1 Λ ( τ , ϖ 1 ( τ ) ) d τ ϖ 0 + Ω ψ η < .
Thus, P 1 ϖ 1 + P 2 ϖ 2 B η .
Further, the contraction of P 2 will be proved.
For t J and ϖ 1 , ϖ 2 B η , one reaches
( P 2 ϖ 1 ) ( t ) ( P 2 ϖ 2 ) ( t ) ϖ 1 ( t 0 ) ϖ 2 ( t 0 ) .
Having Λ as a continuous function, P 1 will also be continuous. Moreover, for any t J and ϖ 1 B η ,
P 1 ϖ Ω ψ < + ,
implies that P 1 is uniformly bounded. At last, one can justify that P 1 is compact. Defining sup ( t , ϖ ) J × B η | Λ ( t , ϖ ( t ) ) | = Λ * , yields
| ( P 1 ϖ ) ( t 2 ) ( P 1 ϖ ) ( t 2 ) | = 1 Γ ( ψ ) | 0 t 1 [ ( t 2 τ ) ψ 1 ( t 1 τ ) ψ 1 ] Λ ( τ , ϖ ( τ ) ) d τ + t 1 t 2 ( t 2 τ ) ψ 1 Λ ( τ , ϖ ( τ ) ) d τ | Λ * Γ ( ψ ) 2 ( t 2 t 1 ) ψ + ( t 2 ψ t 1 ψ ) 0 , a s t 2 t 1 .
This shows that by the Arzelá Ascoli principle, (13) has at least one solution. □

4.2. Iterative and Stability Analysis via Caputo Operator

Let ( B , | | . | | ) be a Banach space and Y * be a self-map of B . Further, consider y n + 1 = h ( Y * , y n ) and G ( Y * ) as a fixed-point set of non-empty Y * . It can be interpreted that y n has converged to q * G ( Y * ) . Defining | | z n + 1 * h ( Y * , z n * ) | | with { z n * B } . The iterative scheme, y n + 1 = h ( Y * , y n ) is Y * -stable if lim n c n = 0 , that is, lim n c n * = p * . To get convergence for z n , it has to possess an upper bound. Having satisfied the aforementioned conditions for y n + 1 = Y * with n representing Picard’s iteration as in [40], then the iteration is Y * -stable. Thus:
Theorem 3
([40]). Let ( B , | | . | | ) be a Banach space and Y * be a self-map on B . Then, for all x , y B , the inequality as comes next holds
Y x * Y y * K x Y x * + k x y ,
with K 0 , k [ 0 , 1 ) . If we assume that Y * is Picard Y * -stable, we can write
S n + 1 ( t ) = S n ( t ) + L 1 1 s α L 1 τ ψ ( 1 S n ) S n X n α ψ ( X n + S n ) , X n + 1 ( t ) = X n ( t ) + L 1 1 s α L β ψ X n τ ψ + S n X n ( X n + S n ) K d ψ X n , C n + 1 ( t ) = C n ( t ) + L 1 1 s α L β ψ C n τ ψ + K I d ψ ( C 0 C n ) S n X n γ ψ ( X n + S n ) .
Theorem 4.
Let F be a self-map. Then, it is defined as follows:
F [ S n ( t ) ] = S n + 1 ( t ) = S n ( t ) + L 1 1 s α L 1 τ ψ ( 1 S n ) S n X n α ψ ( X n + S n ) , F [ X n ( t ) ] = X n + 1 ( t ) = X n ( t ) + L 1 1 s α L β ψ X n τ ψ + S n X n ( X n + S n ) K d ψ X n , F [ C n ( t ) ] = C n + 1 ( t ) = C n ( t ) + L 1 1 s α L β ψ C n τ ψ + K I d ψ ( C 0 C n ) S n X n γ ψ ( X n + S n ) ,
that is F -stable in the space of L 1 ( a , b ) , only if
1 + 1 τ ψ ( 1 κ 1 ( ρ ) ) κ 1 ( ρ ) κ 2 ( ρ ) α ψ ( P + Q ) < 1 , 1 β ψ κ 2 ( ρ ) τ ψ + κ 1 ( ρ ) κ 2 ( ρ ) ( P + Q ) K d ψ κ 2 ( ρ ) < 1 , 1 β ψ κ 3 ( ρ ) τ ψ + K I d ψ ( C 0 κ 3 ( ρ ) ) κ 1 ( ρ ) κ 2 ( ρ ) γ ψ ( P + Q ) < 1 .
Proof. 
Since F is a fixed point, thus, for each ( m , n ) N × N , one attains
F [ S n ( t ) ] F [ S m ( t ) ] = S n + 1 ( t ) = S n ( t ) + L 1 1 s α L 1 τ ψ ( 1 S n ) S n X n α ψ ( X n + S n ) L 1 1 s α L 1 τ ψ ( 1 S m ) S m X m α ψ ( X m + S m ) , F [ X n ( t ) ] F [ X m ( t ) ] = X n + 1 ( t ) = X n ( t ) + L 1 1 s α L β ψ X n τ ψ + S n X n ( X n + S n ) K d ψ X n L 1 1 s α L β ψ X m τ ψ + S m X m ( X m + S m ) K d ψ X m , F [ C n ( t ) ] F [ C m ( t ) ] = C n + 1 ( t ) = C n ( t ) + L 1 1 s α L β ψ C n τ ψ + K I d ψ ( C 0 C n ) S n X n γ ψ ( X n + S n ) L 1 1 s α L β ψ C m τ ψ + K I d ψ ( C 0 C m ) S m X m γ ψ ( X n + S n ) .
Imposing norms on both sides of the first equation in (29), we attain
F ( S n ( t ) ) F ( S m ( t ) ) = S n ( t ) S m ( t ) + L 1 1 s α L 1 τ ψ ( 1 S n ) S n X n α ψ ( X n + S n ) L 1 1 s α L 1 τ ψ ( 1 S m ) S m X m α ψ ( X m + S m ) , ,
and imposing a triangular inequality gives
F ( S n ( t ) ) F ( S m ( t ) ) = S n ( t ) S m ( t ) + L 1 1 s α L 1 τ ψ ( 1 S n ) S n X n α ψ ( X n + S n ) L 1 1 s α L 1 τ ψ ( 1 S m ) S m X m α ψ ( X m + S m ) , .
Thus, (31) becomes
F ( S n ( t ) ) F ( S m ( t ) ) S n ( t ) S m ( t ) + L 1 1 s α L 1 τ ψ ( 1 S n ) S n X n α ψ ( X n + S n ) L 1 1 s α L 1 τ ψ ( 1 S m ) S m X m α ψ ( X m + S m ) .
To achieve our goal, we assume that
X n ( t ) X m ( t ) S n ( t ) S m ( t ) , C n ( t ) C m ( t ) S n ( t ) S m ( t ) .
Inserting (33) into the relation (32), one can have
F ( S n ( t ) ) F ( S m ( t ) ) S n ( t ) S m ( t ) + L 1 1 s α L 1 τ ψ ( 1 S n ) S n X n α ψ ( X n + S n ) L 1 1 s α L 1 τ ψ ( 1 S m ) S m X m α ψ ( X m + S m ) .
Since S m ( t ) , X m ( t ) , and C m ( t ) get convergent and bounded, there exists P > 0 , Q > 0 , and R > 0 for every t. Thereby, we get
S m ( t ) < P , X m ( t ) < Q , C m ( t ) < R , ( m , n ) N × N .
From (34) and (35), one reaches
F ( S n ( t ) ) F ( S m ( t ) ) 1 + 1 τ ψ ( 1 κ 1 ( ρ ) ) κ 1 ( ρ ) κ 2 ( ρ ) α ψ ( P + Q ) ( S n ( t ) S m ( t ) ) ,
with κ 1 , κ 2 , and κ 3 being acquired functions by the inverse Laplace transform in (34). In a similar way, we achieve
F ( Q n ( t ) ) F ( Q m ( t ) ) 1 β ψ κ 2 ( ρ ) τ ψ + κ 1 ( ρ ) κ 2 ( ρ ) ( P + Q ) K d ψ κ 2 ( ρ ) Q n ( t ) Q m ( t ) ,
F ( E n ( t ) ) F ( E m ( t ) ) 1 β ψ κ 3 ( ρ ) τ ψ + K I d ψ ( C 0 κ 3 ( ρ ) ) κ 1 ( ρ ) κ 2 ( ρ ) γ ψ ( P + Q ) E n ( t ) E m ( t ) ,
where the condition (28) is valid. Now, F has a fixed point. To prove that F holds for Theorem 3.3, let (37) and (38) hold, and the following relations are satisfied:
k = ( 0 , 0 , 0 ) , K = 1 + 1 τ ψ ( 1 κ 1 ( ρ ) ) κ 1 ( ρ ) κ 2 ( ρ ) α ψ ( P + Q ) < 1 , 1 β ψ κ 2 ( ρ ) τ ψ + κ 1 ( ρ ) κ 2 ( ρ ) ( P + Q ) K d ψ κ 2 ( ρ ) < 1 , 1 β ψ κ 3 ( ρ ) τ ψ + K I d ψ ( C 0 κ 3 ( ρ ) ) κ 1 ( ρ ) κ 2 ( ρ ) γ ψ ( P + Q ) < 1 .
Hence, we reach the intended result. □

4.3. The Solution’s Positiveness

Let us decide the invariant region and indicate that for all t 0 , all solutions of the model’s fractional system are positive. The main objective is to implement the convenience of the analyzed model solutions by observing whether they join the Υ invariant field. Employing the benefits of the fractional derivative of Caputo, we assume that
Υ = ( S , X , C ) R + 3 , R + 3 = { n R + 3 : n 0 }
is a model solution with positive initial conditions.
Further, n = ( S ( t ) , X ( t ) , C ( t ) ) T . We must also show that on each hyperplane that is bound by the non-negative hyperoctant, the vector field points to R + 3 . Thus, one can write
C D ψ S ( t ) = 1 τ ( 1 S ) 0 , C D ψ Q ( t ) = S X S + X 0 , C D ψ E ( t ) = K I a ( C 0 C ) 0 .
Thus, we give the convenient region as follows:
Υ = { ( S , X , C ) : S 0 , X 0 , C 0 , S + X + C 1 } .
Thus, in the Υ area, the fractional model is biologically acceptable and mathematically well-defined when t > 0 . In addition, this area is positively invariant, that is, the underlying system’s solutions are positive for all t.

5. Numerical Dynamics

This is where, with the aid of the best-suited parameters as obtained in the table, we present different numerical results for the proposed model by using a fractional differential operator known as the Caputo operator. To serve the purpose of model simulations, a predictor-corrector numerical technique for the fractional type of dynamical systems as implemented and analyzed in [41,42,43] has been used on an effective software called Matlab 2019a. In terms of the Caputo differential operator of order ψ , the following form of the Cauchy ordinary differential equation is taken into consideration:
0 C D t ψ φ t = ϖ t , φ t , φ p 0 = φ 0 p , 0 < ψ 1 , 0 < t τ ,
with p = 0 , 1 , , n 1 , n = ψ . The Volterra equation version of the Equation (43) becomes:
φ t = p = 0 n 1 φ 0 p t p p ! + 1 Γ ψ 0 t t s ψ 1 ϖ s , φ s d s .
Applying the described predictor-corrector related with the Adam-Bashforth-Moulton solver [42], we assume h = τ / N , t j = j h , and j = 0 , 1 , , N Z + , by taking φ j φ t j , with φ j standing for the approximate solution and φ ( t j ) for the true one. We can now write
S k + 1 = j = 0 k 1 S 0 j t k + 1 j j ! + h ψ Γ ψ + 2 j = 0 k q j , k + 1 1 τ ( 1 S j ) S j X j α ( S j + X j ) + h ψ Γ ψ + 2 j = 0 k q k + 1 , k + 1 1 τ ( 1 S k + 1 p ) S k + 1 p X k + 1 p α ( S k + 1 p + X k + 1 p ) , X k + 1 = j = 0 k 1 X 0 j t k + 1 j j ! + h ψ Γ ψ + 2 j = 0 k q j , k + 1 β X j τ + S j X j ( S j + X j ) K d X j + h ψ Γ ψ + 2 j = 0 k q k + 1 , k + 1 β X k + 1 p τ + S k + 1 p X k + 1 p ( S k + 1 p + X k + 1 p ) K d X k + 1 p , C k + 1 = j = 0 k 1 C 0 j t k + 1 j j ! + h ψ Γ ψ + 2 j = 0 k q j , k + 1 β C j τ + K I a ( C 0 C j ) S j X j γ ( S j + X j ) + h ψ Γ ψ + 2 j = 0 k q k + 1 , k + 1 β C k + 1 p τ + K I a ( C 0 C k + 1 p ) S k + 1 p X k + 1 p γ ( S k + 1 p + X k + 1 p ) ,
where
q j , k + 1 = k ψ + 1 k ψ k + 1 ψ , if j = 0 , k j + 2 ψ + 1 + k j ψ + 1 2 k j + 1 ψ + 1 , if 1 j k , 1 , if j = k + 1 .
Here, ( p ) implies the predictor portion and it can be determined as:
S k + 1 = j = 0 k 1 S 0 j t k + 1 j j ! + h ψ Γ ψ + 2 j = 0 k w j , k + 1 1 τ ( 1 S j ) S j X j α ( S j + X j ) X k + 1 = j = 0 k 1 X 0 j t k + 1 j j ! + h ψ Γ ψ + 2 j = 0 k w j , k + 1 β X j τ + S j X j ( S j + X j ) K d X j C k + 1 = j = 0 k 1 C 0 j t k + 1 j j ! + h ψ Γ ψ + 2 j = 0 k w j , k + 1 β C j τ + K I a ( C 0 C j ) S j X j γ ( S j + X j )
where
w j , k + 1 = k + 1 j ψ k j ψ .
The numerical simulations of the Caputo model (13) using the parameters mentioned in Table 2 are presented in this section. In order to solve the nonlinear model shown in Equation (13), we used the predictor-corrector Adams technique and obtained graphical results based on the given parameters values. The concentration of substrate S ( t ) , of oxygen C ( t ) , and the microorganism concentration X ( t ) are investigated using different values of the fractional order ψ in order to see the dynamical behavior of the Caputo model (13) for all three compartments. It can be observed in Figure 3 that the behaviour of the curves changes as the value of fractional order increases. Moreover, the values of one of the most important and effective parameters in the model, which is the main experimental control parameter τ , have been varied in Figure 4 to practically see its influence in the model. Comparing Figure 3 and Figure 4, the changes in the curves can be noticed as a result of the effect of the main experimental control, τ .

6. Conclusions

In the current analysis, the wastewater model has been studied from a classical viewpoint and modelled by one of the most robust non-local fractional operators, called the Caputo operator. Firstly, some important properties of the model in a classical sense were reported. The dynamical features of the model were depicted through the steady-state region in its classical sense. Since it has been proved in the literature that the fractional operators are especially suitable for studying the dynamical behaviour of real-world problems, it is of vehement importance to analyze more critically the dynamics of the model under consideration with the robust fractional operator called Caputo, and also make some comparisons with the classical variant. The fractionalized order is ψ . As a consequence, many important characteristics of the proposed fractional version of the model have been reported, such as the formation of the model, the existence and uniqueness of the solution through the fixed-point theorem, iterative solutions, and the solutions’ positivity. It should be observed that the fractional version of the model under investigation more correctly understands the behavior of the model than the variant of the integer order. This has been done by numerous numerical simulations by using an efficient numerical scheme.

Author Contributions

Conceptualization, A.Y.; Data curation, R.T.A. and A.Y.; Formal analysis, R.T.A. and A.Y.; Funding acquisition, R.T.A.; Investigation, R.P.A.; Methodology, R.T.A. and A.Y.; Visualization, R.P.A.; Writing–original draft, R.P.A.; Writing–review–editing, R.P.A. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

All data used in this work can be found within the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

SymbolDescription
FThe flow rate through the reactor.
SThe substrate concentration within the reactor.
S * The dimensionless substrate concentration within the reactor [ S * = S / S 0 ].
CThe dissolved oxygen concentration.
C * -The dimensionless dissolved oxygen concentration C S 0 .
K I a Volumetric oxygen mass transfer coefficient
C 0 The dissolved oxygen solubility limit
C 0 * The dimensionless dissolved oxygen solubility limit
K d The death coefficient
K d * The dimensionless death coefficient
K s The saturation constant
S 0 The concentration of substrate flowing into the reactor.
VThe volume of the reactor.
Y x l s The yield factor for substrate
α * The dimensionless yield factor for substrate
Y x l o 2 The yield coefficient for oxygen
γ * The dimensionless yield coefficient for oxygen.
μ m The specific growth rate parameter
β The reactor model for the reactor.
XThe biomass concentration within the reactor
X * The dimensionless biomass concentration within the reactor
tThe time.
t * The dimensionless time.
μ ( S , X ) The specific growth rate model.
τ The residence time.
τ * The dimensionless residence time.

References

  1. Brauer, F.; Castillo-Chavez, C.; Feng, Z. Mathematical Models in Epidemiology; Springer: New York, NY, USA, 2019. [Google Scholar]
  2. Qureshi, S.; Yusuf, A. Modeling chickenpox disease with fractional derivatives: From caputo to atangana-baleanu. Chaos Solitons Fractals 2019, 122, 111–118. [Google Scholar] [CrossRef]
  3. Qureshi, S.; Yusuf, A. Mathematical modeling for the impacts of deforestation on wildlife species using Caputo differential operator. Chaos Solitons Fractals 2019, 126, 32–40. [Google Scholar] [CrossRef]
  4. Jajarmi, A.; Yusuf, A.; Baleanu, D.; Inc, M. A new fractional HRSV model and its optimal control: A non-singular operator approach. Phys. A 2020, 547, 123860. [Google Scholar] [CrossRef]
  5. Abdulhameed, M.; Muhammad, M.M.; Gital, A.Y.; Yakubu, D.G.; Khan, I. Effect of fractional derivatives on transient MHD flow and radiative heat transfer in a micro-parallel channel at high zeta potentials. Phys. A 2019, 519, 42–71. [Google Scholar] [CrossRef]
  6. VDubey, P.; Kumar, R.; Kumar, D. Analytical study of fractional Bratu-type equation arising in electro-spun organic nanofibers elaboration. Phys. A 2019, 521, 762–772. [Google Scholar]
  7. Chang, A.; Sun, H.; Zhang, Y.; Zheng, C.; Min, F. Spatial fractional Darcy’s law to quantify fluid flow in natural reservoirs. Phys. A 2019, 519, 119–126. [Google Scholar] [CrossRef]
  8. Goulart, A.G.; Lazo, M.J.; Suarez, J.M.S. A new parameterization for the concentration flux using the fractional calculus to model the dispersion of contaminants in the Planetary Boundary Layer. Phys. A 2019, 518, 38–49. [Google Scholar] [CrossRef] [Green Version]
  9. Ochoaa, F.G.; Gomez, E.; Santos, V.E.; Merchuk, J.C. Oxygen uptake rate in microbial processes: An overview. Biochem. Eng. J. 2010, 49, 289–307. [Google Scholar] [CrossRef]
  10. Garcia-Ochoa, F.; Gomez, E.; Santos, V. Oxygen transfer and uptake rates during xanthan gum production. Enzyme Microb. Technol. 2000, 27, 680–690. [Google Scholar] [CrossRef]
  11. Garcia-Ochoa, F.; Gomez, E. Mass transfer coefficient in stirrer tank reactors for xanthan solutions. Biochem. Eng. J. 1998, 1, 1–10. [Google Scholar] [CrossRef]
  12. Ho, C.; Oldshue, J.Y. (Eds.) Biotechnology Processes: Scale-Up and Mixing; AIChE: New York, NY, USA, 1987. [Google Scholar]
  13. Andrew, S.P.S. Gas–liquid mass transfer in microbiological reactors. Trans. IChemE 1982, 60, 3–13. [Google Scholar]
  14. Aiba, S.; Humphrey, A.E.; Mills, N.F. (Eds.) Biochemical Engineering, 2nd ed.; Academic Press: New York, NY, USA, 1973. [Google Scholar]
  15. Zou, X.; Hang, H.; Chu, J.; Zhuang, Y.; Zhang, S. Oxygen uptake rate optimization with nitrogen regulation for erythromycin production and scale-up from 50 L to 372 m3 scale. Bioresour. Technol. 2009, 100, 1406–1412. [Google Scholar] [CrossRef] [PubMed]
  16. Palomares, L.A.; Ramirez, O. The effect of dissolved oxygen tension and the utility of oxygen uptake rate in insect cell culture. Cytotechnology 1996, 22, 225–237. [Google Scholar] [CrossRef] [PubMed]
  17. Bhattacharya, S.C.; Khai, P.Q. Kinetics of Anaerobic Cowdung Digestion. Energy 1987, 12, 497–500. [Google Scholar] [CrossRef]
  18. Serhani, M.; Gouze, J.L.; Raissi, N. Dynamical study and robustness of a nonlinear wastewater treatment problem. J. Nonlinear Anal. RWA 2011, 12, 487–500. [Google Scholar] [CrossRef]
  19. Haandel, A.V.; Lubbe, J.V. Handbook Biological Waste Water Treatment; IWA Publishing: Leidschendam, The Netherlands, 2007. [Google Scholar]
  20. Alqahtani, T.R.; Nelson, I.M.; Worthy, L.A. Analysis of a chemostat model with variable yield coefficient: Contois kinetics. ANZIAM J. 2012, 53, C155–C171. [Google Scholar] [CrossRef] [Green Version]
  21. Alqahtani, T.R.; Nelson, I.M.; Worthy, L.A. A fundamental analysis of continuous flow bioreactor models with recycle around each reactor governed by Contois kinetics. III. Two and three reactor cascades. Chem. Eng. J. 2012, 183, 422–432. [Google Scholar] [CrossRef]
  22. Koumboulis, F.N.; Kouvakas, N.D.; King, R.E.; Stathaki, A. Two-stage robust control of substrate concentration for an activated sludge process. ISA Trans. 2008, 47, 267–278. [Google Scholar] [CrossRef]
  23. Jourani, A.; Serhani, M.; Boutoulout, A. Dynamic and controllability of a nonlinear wastewater treatment problem. J. Appl. Math. Inform. 2012, 30, 883–902. [Google Scholar]
  24. Khan, K.; Zarin, R.; Khan, A.; Yusuf, A.; Al-Shomrani, M.; Ullah, A. Stability analysis of five-grade Leishmania epidemic model with harmonic mean-type incidence rate. Adv. Differ. Equ. 2021, 86. [Google Scholar] [CrossRef]
  25. Serhani, M.; Cartigny, P.; Raissi, N. Robust feedback design of wastewater treatment problem. Math. Model. Nat. Phenom. 2009, 4, 1139–1143. [Google Scholar] [CrossRef]
  26. Haegeman, B.; Lobry, C.; Harmand, J. Modeling bacteria flocculation as density-dependent growth. AIChE J. 2007, 53, 535–539. [Google Scholar] [CrossRef]
  27. Singh, J.; Kumar, D.; Baleanu, D. On the analysis of fractional diabetes model with exponential law. Adv. Differ. Equ. 2018, 2018, 231. [Google Scholar] [CrossRef] [Green Version]
  28. Yusuf, A.; Acay, B.; Mustafa, U.T.; Inc, M.; Baleanu, D. Mathematical modeling of pine wilt disease with Caputo fractional operatorChaos. Solitons Fractals 2021, 143, 110569. [Google Scholar] [CrossRef]
  29. Musa, S.S.; Qureshi, S.; Zhao, S.; Yusuf, A.; Mustapha, U.T.; He, D. Mathematical modeling of COVID-19 epidemic with effect of awareness programs. Infect. Dis. Model. 2021, 6, 448–460. [Google Scholar]
  30. Baba, I.A.; Yusuf, A.; Al-Shomrani, M. A mathematical model for studying rape and its possible mode of control. Results Phys. 2021, 22, 103917. [Google Scholar] [CrossRef]
  31. Acay, B.; Inc, M. Fractional modeling of temperature dynamics of a building with singular kernels. Chaos Solitons Fractals 2021, 142, 110482. [Google Scholar] [CrossRef]
  32. Ahmed, I.; Goufo, E.F.D.; Yusuf, A.; Kumam, P.; Chaipanya, P.; Nonlaopon, K. An epidemic prediction from analysis of a combined HIV-COVID-19 co-infection model via ABC fractional Operator. Alex. Eng. J. 2021, 60, 2979–2995. [Google Scholar] [CrossRef]
  33. Acay, B.; Bas, E.; Abdeljawad, T. Non-local fractional calculus from different viewpoint generated by truncated M-derivative. J. Comput. Appl. Math. 2020, 366, 112410. [Google Scholar] [CrossRef]
  34. Kirtphaiboon, S.; Humphries, U.; Khan, A.; Yusuf, A. Model of rice blast disease under tropical climate conditions. Chaos Solitons Fractals 2021, 143, 110530. [Google Scholar] [CrossRef]
  35. Noeiaghdam, S.; Sidorov, D.N. Caputo-Fabrizio Fractional Derivative to Solve the Fractional Model of Energy Supply-Demand System. Math. Model. Eng. Probl. 2020, 7, 359–367. [Google Scholar] [CrossRef]
  36. Hu, W.C.; Thayanithy, K.; Forster, C.F. A kinetic study of the anaerobic digestion of ice-cream wastewater. Process. Biochem. 2002, 37, 965–971. [Google Scholar] [CrossRef]
  37. Fikar, M.; Chachuat, B.; Latifi, M.A. Optimal operation of alternating activated sludge processes. Control Eng. Pract. 2005, 13, 853–861. [Google Scholar] [CrossRef] [Green Version]
  38. Henze, M.; Grady, C.P.L., Jr.; Gujer, W.; Marais, G.V.R.; Matsuo, T. A general model for single-sludge wastewater treatment systems. Water Res. 1987, 21, 505–515. [Google Scholar] [CrossRef] [Green Version]
  39. Yoon, S.-H.; Lee, S. Critical operational parameters for zero sludge production in biological wastewater treatment processes combined with sludge disintegration. Water Res. 2005, 39, 3738–3754. [Google Scholar] [CrossRef]
  40. Wang, Z.; Yang, D.; Ma, T.; Sun, N. Stability analysis for nonlinear fractional-order systems based on comparison principle. Nonlinear Dyn. 2014, 75, 387–402. [Google Scholar] [CrossRef]
  41. Diethelma, K.; Ford, N.J.; Freed, A.D. A Predictor-Corrector Approach for the Numerical Solution of Fractional Differential Equations. Nonlinear Dyn. 2002, 29, 3–22. [Google Scholar] [CrossRef]
  42. Diethelm, K. An algorithm for the numerical solution of differential equations of fractional order. Electron. Trans. Numer. Anal. 1997, 5, 1–6. [Google Scholar]
  43. Diethelm, K.; Ford, N.J.; Freed, A.D. Detailed error analysis for a fractional Adams method. Numer. Algorithms 2004, 36, 31–52. [Google Scholar] [CrossRef] [Green Version]
Figure 1. A steady-state diagram of the oxygen concentration, the substrate concentration, and the microorganism concentration varying the residence time τ for different values of parameter β . The value of the parameters stated in Table 1. β = 0.1 (Black color), β = 0.5 (Blue color), β = 1 (Red color).
Figure 1. A steady-state diagram of the oxygen concentration, the substrate concentration, and the microorganism concentration varying the residence time τ for different values of parameter β . The value of the parameters stated in Table 1. β = 0.1 (Black color), β = 0.5 (Blue color), β = 1 (Red color).
Mathematics 09 00675 g001aMathematics 09 00675 g001b
Figure 2. Unfolding diagram for the minimum residence time in Figure 1a as a function of β and the influent concentration (S0). The value of the parameters stated in Table 1. (a) Minimum residence time versus parameter S0. (b) Minimum residence time versus parameter S0.
Figure 2. Unfolding diagram for the minimum residence time in Figure 1a as a function of β and the influent concentration (S0). The value of the parameters stated in Table 1. (a) Minimum residence time versus parameter S0. (b) Minimum residence time versus parameter S0.
Mathematics 09 00675 g002
Figure 3. Numerical simulation of the oxygen saturation concentration, the substrate concentration, and the microorganism concentration for different values of fractional order. The value of the parameters is stated in Table 1. (a) Dynamic of the concentration of substrate. (b) Dynamic of the substrate concentration. (c) Dynamic of the microorganism concentration.
Figure 3. Numerical simulation of the oxygen saturation concentration, the substrate concentration, and the microorganism concentration for different values of fractional order. The value of the parameters is stated in Table 1. (a) Dynamic of the concentration of substrate. (b) Dynamic of the substrate concentration. (c) Dynamic of the microorganism concentration.
Mathematics 09 00675 g003
Figure 4. Numerical simulation of the oxygen saturation concentration, the substrate concentration, and the microorganism concentration for different values of residence time. The value of the parameters is stated in Table 1. (a) Dynamic of the concentration of substrate with different values of τ. (b) Dynamic of the substrate concentration with different values of τ. (c) Dynamic of the microorganism concentration with different values of τ.
Figure 4. Numerical simulation of the oxygen saturation concentration, the substrate concentration, and the microorganism concentration for different values of residence time. The value of the parameters is stated in Table 1. (a) Dynamic of the concentration of substrate with different values of τ. (b) Dynamic of the substrate concentration with different values of τ. (c) Dynamic of the microorganism concentration with different values of τ.
Mathematics 09 00675 g004
Table 1. Parameters’ values.
Table 1. Parameters’ values.
ParametersValuesUnitReference
μ max 0.9297 day 1 [36]
K d 0.0131 day 1 [36]
Y x / s 0.2116 g g 1 [36]
K s 0.4818 g / L [36]
K L a 108 day 1 [37]
Y x / O 2 0.67 g g 1 [38]
S 0 200 m g / L [39]
C 0 10 m g / L [37]
β 0–1-Range value of the parameter
Table 2. Parameters’ values.
Table 2. Parameters’ values.
ParametersValues
K d 0.0140905
α 0.10194888
K I a 116.1665
γ 0.322806
C 0 1 20
β 0–1
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Alqahtani, R.T.; Yusuf, A.; Agarwal, R.P. Mathematical Analysis of Oxygen Uptake Rate in Continuous Process under Caputo Derivative. Mathematics 2021, 9, 675. https://doi.org/10.3390/math9060675

AMA Style

Alqahtani RT, Yusuf A, Agarwal RP. Mathematical Analysis of Oxygen Uptake Rate in Continuous Process under Caputo Derivative. Mathematics. 2021; 9(6):675. https://doi.org/10.3390/math9060675

Chicago/Turabian Style

Alqahtani, Rubayyi T., Abdullahi Yusuf, and Ravi P. Agarwal. 2021. "Mathematical Analysis of Oxygen Uptake Rate in Continuous Process under Caputo Derivative" Mathematics 9, no. 6: 675. https://doi.org/10.3390/math9060675

APA Style

Alqahtani, R. T., Yusuf, A., & Agarwal, R. P. (2021). Mathematical Analysis of Oxygen Uptake Rate in Continuous Process under Caputo Derivative. Mathematics, 9(6), 675. https://doi.org/10.3390/math9060675

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop