Numerical Solutions of Fractional Differential Equations Arising in Engineering Sciences

: This paper deals with the numerical solutions of a class of fractional mathematical models arising in engineering sciences governed by time-fractional advection-diffusion-reaction (TF–ADR) equations, involving the Caputo derivative. In particular, we are interested in the models that link chemical and hydrodynamic processes. The aim of this paper is to propose a simple and robust implicit unconditionally stable ﬁnite difference method for solving the TF–ADR equations. The numerical results show that the proposed method is efﬁcient, reliable and easy to implement from a computational viewpoint and can be employed for engineering sciences problems. in the characterization and validation of solution for the packed bed columns model governed by a time-fractional advection-diffusion-reaction equation.


Introduction
In the last years, the fractional differential equations (FDEs) have attracted considerable interest by the numerous researchers. Although the FDEs have an intricate mathematical background, they are widely applied in various fields of applied sciences and engineering to describe nonlinear phenomena which are of great importance in scientific and technological fields. A wide and detailed collection of real world applications of the FDEs calculus in science and engineering is proposed in [1].
The great potential and applicability of the FDEs have motivated the development of new and efficient tools for the search of exact and numerical solutions of both stationary and dynamic fractional differential equations. For a better understanding of the historical aspects of fractional calculus, the analytical properties of the fractional differential equations, and the difficulties associated with them, the readers are referred to the books [2,3]. Because of the impossibility of exactly solving most of these problems, approximate solutions are of academical and practical importance. Therefore, new methods for obtaining approximate solutions of fractional nonlinear partial and ordinary differential equations have been developed. Recently, several methods have drawn special attention such as Adomian decomposition method [4][5][6], variational iteration method [7,8], homotopy analysis method [9][10][11][12], homotopy perturbation method [13][14][15][16], wavelet methods [17,18], finite difference methods [19,20], finite volume methods [21], finite element methods [22,23] and spectral methods [24]. Recently, a new approach has been developed [25][26][27][28][29] to find exact and numerical solutions of the fractional advection-diffusion-reaction equations involving Riemann-Liouville derivative by means the combination of the Lie symmetries theory and finite difference methods.
In this paper, we consider a class of mathematical models that, in recent years, has been a great deal of interest in fractional differential equations. These equations have important applications in various areas such as nonlinear hydrodynamics, population dynamic, biophysics, engineering, neurosciences, polymer physics, laser physics, plasma physics, surface physics, pattern formation, psychology, and marketing. In particular, we are interested in the models that link chemical and hydrodynamic processes. Combining physical and chemical processes into the same mathematical model, we obtain the time-fractional advection-diffusion-reaction (TF-ADR) equations used to describe a wide variety of nonlinear processes arising in many fields of the applied sciences. The governing partial differential equations are characterized by an advective transport term, a parabolic dissipative diffusive (dispersive) term and, possibly, an additional decay/reaction mechanism. The aim is to propose a simple and efficient implicit unconditionally stable finite difference method for solving a class of TF-ADR equations arising in engineering sciences. The proposed numerical method was applied for solving TF-ADR models defined on nonuniform grids and its consistency, unconditional stability and convergence properties were discussed [20]. In this paper, a mathematical model to predict the dynamics in a packed bed column is presented. A packed bed reactor is an assembly of usually uniformly sized particles, which are randomly arranged within a vessel or tube. The bulk fluid flows through the voids of the bed and the reactants are transported at first from the bulk of the fluid to the catalyst surface, then through catalyst pores, where the reactants adsorb on the surface of the pores and then undergo chemical transformation. When a fluid is flowing through a bed of inert particles, one observes the dispersion of the fluid, consequence of the combined effects of convection in the spaces between particles and molecular diffusion. Several works have been done on the description of the principles of solute transport in porous media in packed bed reactors and different models of various levels of complexity have been proposed in the specialized literature. An exact mathematical description of packed bed reactors is practically impossible because of the interactions of fluid mechanics, heat, and mass transport with chemical reactions. Then, determining the exact solution to characterize the flowing fluid through one of these structures becomes impossible. For this reason, numerical solutions of mathematical models have received a great attention over the years. The numerical solutions obtained reveal that the approach is easy to implement and accurate. The reported results show that the proposed numerical method is very simple and effective to perform for engineering problems.
The manuscript is organized as follows. In Section 2, the mathematical model is presented; in Section 3, the numerical scheme is discussed and its consistency, unconditionally stability, and convergence property are presented; in Section 4, an application to a model widely used in the applied sciences is given to illustrate the effectiveness of the proposed numerical method; and, finally, in Section 5, conclusions are given.

The Mathematical Model
We consider the class of TF-ADR equations of the type on the domain x ∈ (a, b) for t > 0, with the boundary conditions and the initial condition of the form In the model, c = c(x, t) represents the variable concentration field of one space variable x at time t, and which time-fractional derivative ∂ α c ∂t α (x, t) is given by means of the α order Caputo fractional derivative defined by where Γ(·) is the well-known Gamma function. The functions u(x, t) and D(x, t) > 0, ∀x, t, represent the velocity field and the diffusion or dispersion term, respectively. In this context, we assume them, in general, as assigned functions of the time t and space x. The function f (c) is used in order to represent source and sink terms. φ(x), ϕ(t), and ψ(t) are known smooth functions. We assume the problem (1) has a unique and sufficiently smooth solution under the above initial and boundary conditions, and that this solution is non-negative solution since negative solutions have no physical interpretation. The class of TF-ADR models (1) includes a wide variety of equations. The Equation (1) reduces to the time-fractional nonlinear reaction-diffusion equation when u(x, t) = 0 and reduces to the time-fractional nonlinear gas dynamics equation when D(x, t) = 0. Moreover, if f (c) = 0, it can be easily seen that Equation (1) is a special case of the well-known Fokker-Planck equation, where the diffusion and the drift coefficients are functions of the space and time. When α = 1, the model (1) leads to the classical advection-diffusion-reaction model used in order to describe several phenomena of relevant interest in many fields of applied sciences [30][31][32][33].

An Implicit Finite Difference Method
In this section, we present the unconditionally stable implicit finite difference method, proposed for solving the TF-ADR model (1).
For the derivation of the proposed method, first, we construct a computational uniform grid in the x and t directions: we define the spatial size of the mesh ∆x = x j − x j−1 , for 1 ≤ j ≤ J, and the temporal size of the mesh ∆t = t n − t n−1 , for 1 ≤ n ≤ N, with J and N positive integers. We define the mesh points (x j , t n ) with x j = a + j∆x, j = 0, · · · , J and t n = n∆t n , for n = 0, · · · , N.
We denote by C n j the numerical approximation provided by the difference method of the exact solution c(x j , t n ) at the mesh points (x j , t n ), for j = 0, · · · , J and n = 0, · · · , N.
According to (2), we discretize the time-fractional Caputo derivative of the function c(x, t), evaluated at the mesh point (x j , t n ), by means of the following classical approximate formula This is the well-known so-called L1 formula defined in [34] with R n = O(∆t 2−α ) local truncation error (see [35,36]). From the truncation error estimate of the L1 formula, it is important to note that the accuracy is dependent on the value of the fractional order α. This is justified by a weakly singular kernel (t − s) −α that is contained in the integral.
We assume, as usual, that the solution is sufficiently smooth and discretize its first ∂/∂x and second order ∂ 2 /∂x 2 spatial derivatives by the second-order three-point central finite difference formula so that By replacing c(x j , t n ) with its numerical approximation C n j and by neglecting the local truncation errors, the time-fractional advection-diffusion Equation (1) is discretized in the following way , and finally, The initial and boundary conditions can be rewritten as It holds that [37], for any temporal meshes on [0, T] and for any n = 1, 2, · · · , N, Taking into account that T n,n = ∆t −α , we can write then, we obtain the following implicit finite difference method where we set In general, by the integration of fractional differential equations by using implicit numerical methods, at each step, a nonlinear equation could be determined, therefore it is necessary to use an algorithm to solve nonlinear equations, such as the classic Newton method. Newton's method is an iterative method and is one of the most used approaches. At each iteration it requires the evaluation of the Jacobian matrix that can be evaluated analytically but, often, it is calculated numerically, if too expensive by the computational point of view.
The Equation (8) can be written in vectorial form as where K is a tridiagonal matrix and the difference operator L is defined as follows, Here, and in the following, we assume the convention that, if the lower bound is larger that the upper bound, the summation is equal to zero. The obtained method (9) is implicit, to compute the numerical solution C n a system with the tridiagonal coefficients matrix K has to be solved.
Note that the operator L is a kind operator with memory, due to the non-local character of the fractional derivative, this means that the effect on C at time t n , C n , depends on all the previous values, C 0 , C 1 , · · · , C n−1 , evaluated at all the previous time t 0 , t 1 , · · · , t n−1 . The main difference compared to the non-fractional case is that, to evaluate L, the numerical solutions for all the n previous time values t 0 , t 1 , · · · , t n−1 are required, whereas for non-fractional equations, only the solution to the previous value t n−1 is used. The computational cost to obtain the solution at the time t n from the solution at the time t n−1 increases as n, that is, increases as the number of terms in the summation that compares in the second term of the (10). This implies that the computational effort to go from t 0 to t n grows as n 2 .

Consistency, Stability and Convergence
In this Section, we discuss the consistency, the stability and the convergence of the implicit finite difference method.
Consistency. According to Equations (3)-(5), the local truncation error R n j of the finite difference method (6) is The implicit finite difference scheme defined by Equations (6) or (8) is consistent with the problem (1) of order O(∆t 2−α + ∆x 2 ).
Stability. For the stability analysis of the implicit finite difference method, we rewrite Equation (8) as follows LetC n j be another approximate solution of the finite difference scheme (8), and let be the corresponding round-off error. We let ρ n = (ρ n 0 , ρ n 1 , · · · , ρ n J ) T , and we consider the infinity norm ||ρ n || ∞ = max 0≤j≤J |ρ n j | = |ρ n i | .
The round-off error satisfies the following round-off equations, To check whether the finite difference method is stable, we study how the size of the round-off error ρ n evolves in time. We define Equation (13) can be written as From (13) and taking into account that T n,k+1 − T n,k > 0, we have where we have defined However, as ∑ n−1 k=0 (T n,k+1 − T n,k ) = ∆t −α and recalling that T n,n = ∆t −α and T n,0 = 0, it has Thus, we can conclude that i.e., ||ρ n || ∞ ≤ ||ρ 0 || ∞ , for 1 ≤ n ≤ N. This means that the proposed method is unconditionally stable. Convergence. Let, c(x j , t n ) be the exact solution of Equation (1) at mesh point (x j , t n ) for j = 0, 1, · · · , J and n = 0, 1, · · · , N. Denoting n j = c(x j , t n ) − C n j , we get the error equations with 0 j = 0, for j = 0, 1, · · · , J and n J = 0, for n = 1, 2, · · · , N. R n j is the local truncation error (11) of the numerical scheme. We introduce the following norm where R max = max n,i |R n i |. For n = 1, we have with C R positive constant dependent on T, α, and the exact solution c(x, t), but independent of ∆t and ∆x. Then, we prove that the solution of the finite difference scheme (8), with initial and boundary conditions given by (7), is convergent. The proposed numerical method was validated in terms of consistency, accuracy, and convergence in [20], where it was implemented on a nonuniform grid and used for solving mathematical models with known exact solutions.

Application to the Chemical Engineering
In this section, we show the applicability and reliability of the proposed numerical method by solving a TF-ADR model arising in chemical engineering. All the results are calculated by using the software MatLab.
Dispersion in porous media has been studied by a relevant number of researchers for its important role in describing of many phenomena, such as, for example, in contaminant transport in ground water flows, in miscible displacement of oil and gas and in reactant, product transport in packed bed reactors, and so on. The phenomenon of dispersion (transverse and longitudinal) in packed beds is an important phenomenon of dynamics in porous media.
A packed bed reactor is an assembly of usually uniformly sized particles, which are randomly arranged within a vessel or tube. The bulk fluid flows through the voids of the bed. The reactants are transported at first from the bulk of the fluid to the catalyst surface, then through catalyst pores, where the reactants adsorb on the surface of the pores and then undergo chemical transformation.
When a fluid is flowing through a bed of inert particles, one observes the dispersion of the fluid largely caused by the complex flow patterns in the reactor induced by the presence of the packing. Dispersion effects are the consequence of the combined effects of convection in the spaces between particles and molecular diffusion. Generally, the dispersion coefficient in radial direction is inferior to the dispersion coefficient in longitudinal direction. Several works have been done on the description of the principles of solute transport in porous media in packed bed reactors and different models of various levels of complexity have been proposed in the specialized literature. An exact mathematical description of packed bed reactors is practically impossible because of the interactions of fluid mechanics, heat and mass transport with chemical reactions. In fact, in the study of dispersion in packed beds there are several variables that must be considered, such as the structure of a porous medium, the length of the packed column, viscosity and density of the fluid, effect of fluid velocity and effect of temperature, ratio between particle diameter and column diameter, particle size distribution particle shape, and ratio of column length to particle diameter. Determining the exact solution to characterize the flowing fluid through one of these structures becomes then impossible. For this reason, numerical solutions of mathematical models have received a great attention over the years. In this paper, we propose an efficient unconditionally stable implicit finite difference method (8) for solving a dispersion model in a packed bed reactor.
We consider a packed bed, with uniform porosity , of length L and diameter d, filled with porous solid absorbent particles of uniform radius R. At time t = 0, a liquid, generally an aqueous solution, enters the packed bed, at an interstitial velocity u and initial concentration of solute c 0 , in which a tracer, with concentration of solute c and continuously injection, is dispersed in axial and radial directions. In this context, for a small ratio of column diameter to length (d/L) and large fluid velocity, the transverse dispersion is neglected in comparison with axial dispersion. Moreover, we assume fluid viscosity and density constant, isothermal operation, constant fluid flow rate, and mass transfer rate in the adsorption process described by linear driving force rate term. The species is transferred convectively from the bulk of the flowing solution to the external surface of the particles. Species then diffuses occupying the void fraction of the particles, where it is progressively adsorbed. In Figure 1, we schematically illustrate the behavior of a packed bed reactor [38]. To describe the temporal variation of order α of the concentration of the solute at any position in the bulk solution occupying the void space v of the adsorption packed bed, we propose the following mathematical model [39,40], Equation (17) describes the mass balance for a single component, as, for example, ethanol, glucose, glycerol, acetic acid, and so on, taking into account the axial dispersion, the bulk flow through the column, and convective mass transfer from the bulk to the surface of the particles. x represents the axial position, c(x, t) is the liquid-phase concentration of the species, and q(x, t) is the average adsorbed-phase concentration and ρ the particle density. The first term on the left hand side describes the fractional order time-change of the concentration of the species c within the interstitial space, the second term represents the advection of the species in the pores of spherical adsorbent particles, whereas the last term on the left-hand side includes an effective dispersion (not diffusion) coefficient D L . Finally, the term of the right-hand side incorporates the reactions occurring on the surface.
The main difference, between the mathematical model under study (17) and the classical formulation by the classical PDE, lies in using a time fractional partial differential operator of α order in substitution of the integer partial differential operator. By using the fractional derivative, we introduce a memory formalism in our model, a good instrument to describe nonlocal effects. As known, the fractional modeling allows an adequate representation of the dynamics of anomalous diffusion.
Since the rate of accumulation of solute in the solid surface is equal to the rate of transfer of solute across the liquid film [41], the mass balance for the single species is with c * the liquid-phase concentration at the particle surface. Now, model (17) becomes The reaction term is described by linear driving force rate term and describes a mass transfer rate in the adsorption process. Here, k c is the mass transfer coefficient of species and a v is the volumetric surface area of the packed bed. c * represents the hypothetical concentration of the solution that would be in equilibrium with the concentration of the solute inside the packed bed, i.e., solubility. In this context, we assume that the uptake rate of a species is proportional to the linear difference between the concentration of that species at the outer surface of the particle (equilibrium adsorption amount) and its average concentration within the particle. This means that species is adsorbed instantaneously when it reaches the surface of the particles in a way that it becomes immediately in equilibrium with the concentration of species in the solution at the surface of the particle. See the work in [42] for more details. The coupling equation between the liquid and solid concentration is the so-called equilibrium isotherm.
It is assumed that, at time t = 0, the concentration c(x, t) in the liquid occupying the void is equal to zero, ∀x > 0, and that, at the same time, the concentration is subjected to a step change from zero to c 0 at the entrance x = 0 of the packed bed The model under study is a boundary value problem and requires two boundary conditions both for the inlet and the outlet of the reactor. Then, at the entrance x = 0 of the packed bed, we assume that there is an assigned prefixed value of the species and, at the exit x = L, there is no diffusion, so that the two boundary conditions can be expressed as follows, according to [43] where c in is a specified quantity. Note that, for the example under study, being f (c) a linear function respect to the field variable c, by using the numerical method (8), we obtain a system of linear equations for which solution the Newton's iteration method does not need. Moreover, we set c 0 = c in = 1.7. The data used for numerical computation is v = 0.7, u = 1.3, D L = 1, k c = 1.5, a v = 1.2, and c * = 2.5. In this context, the coefficients that appear in Equation (18) are assumed constant. However, the numerical method impose no limitations in the choice of coefficients. In Figure 2, the simulated concentration profiles at different time steps are presented. The concentration slowly growths with increasing time and it approaches its equilibrium adsorption value. In Figure 3, we report the numerical solutions obtained with different parameter values α. Note that, for values of α approaching 1, the numerical solution tends to the equilibrium value, it reaches the exact equilibrium values only for α = 1. In Figure 4, the breakthrough curves obtained by plotting the outflow/inflow concentration ratio (C/C 0 ) for different values of fractional α at final time t N are reported. The solution behavior is in agreement with the results available in the specialized literature [38,40,41,44].  Example 2: In this second example, to validate the reliability of the proposed numerical method, we compute the numerical solution by setting the parameters values according to the test presented in [41]: c 0 = c in = c * = 0.2 (mM), v = 0.5196, D L = 3.8 × 10 −8 (m 2 /s), u = 1.061 × 10 −4 (m/s), k c = 10 −3 (s −1 ) and, finally, a v = 10 −2 (m). To simulate the breakthrough curves, it is required to calculate the column exit concentration of species as a function of time t on a computational domain [0, T] × [0, L], with T = 4000 (s) and L = 6.7 × 10 −2 (m) with J = N = 100 grid points in both directions. The predicted breakthrough curves are plotted in Figure 5, for different values of the fractional order α. Note the effect of the fractional order α on the behavior of the solution: the smaller the value of α, the longer is the time needed to reach the equilibrium state. Moreover, we observe that what is seen in the fractional model takes longer to be seen respect to the integer order model. The simulated concentration profiles predict accurately the breakthrough curves of solute obtained by Escudero et al. [41] (see Figure 1). It is important to note that our results are in agreement with the numerical and experimental ones obtained by Escudero et al. when α = 1. Example 3: Finally, a numerical test is considered to simulate the hydrogen penetration process in the flow bed reactor, filled with a Zr 2 Fe powder. According to the model studied by Jiangfeng et al. in [44], we compute the numerical solution by setting the parameters values as follows, c 0 = c in = c * = 1.1 (mM), v = 0.8, D L = 1.63 × 10 −4 (m 2 /s), u = 0.02 (m/s), k c = 10 −3 (s −1 ), and a v = 0.17 (m). We simulate the hydrogen concentration for different values of fractional order α = 0.5, 0.8, 1. In Figures 6-8    By comparison with the experimental and numerical results available in the specialized literature [38,40,41,44], we can conclude that the proposed computational method is an efficient and accurate tool in the characterization and validation of solution for the packed bed columns model governed by a time-fractional advection-diffusion-reaction equation.

Concluding Remarks
In this paper, we propose an unconditionally stable implicit finite difference method for solving a class of mathematical models governed by time-fractional partial differential equations. The main aim of this paper is to show how the proposed numerical method is an efficient and reliable tool for solving TF-ADR models characterized by an advective transport term, a parabolic diffusive (dispersive) term, and an additional reaction/decay mechanism. In particular, we present a numerical study of a TF-ADR model used to predict and simulate the dynamics of a single component in the packed bed column. The numerical results provide valuable and highly accurate approximations to the breakthrough curves for designing packed bed adsorption processes for field applications [38,40,41,44].
Funding: This research received no external funding.