International Journal of Analysis and Applications Volume 19, Number 4 (2021), 604-618 URL: https://doi.org/10.28924/2291-8639 DOI: 10.28924/2291-8639-19-2021-604 A POSSIBLE GENERALIZED MODEL OF THE CHLORINE CONCENTRATION DECAY IN PIPES: EXACT SOLUTION YUSSRI M. MAHROUS∗ Department of Studies and Basic Sciences, Faculty of Community, University of Tabuk, Saudi Arabia ∗Corresponding author: y.mahrous@ut.edu.sa; mahrous1972@yahoo.com Abstract. This paper discusses a possible generalization of the transport model describing the chlorine concentration decay in pipes. The proposed generalized model is governed by a second-order fractional partial differential equation. The exact solution of the generalized model is obtained via the Laplace transform method and the method of residues. The exact solution reduces to the corresponding published one as the fractional order α tends to one. Analytical expression for the dimensionless cup-mixing average concentration is deduced. Influences of various parameters on the behavior of the dimensionless cup-mixing average concentration are discussed. It is shown that the physical interpretation of the dimensionless cup-mixing average concentration in view of the fractional calculus is completely different than its interpretation in the classical calculus. 1. Introduction Chlorine is the most commonly employed disinfectant in most countries and minimum levels of chlorine must be maintained to ensure the disinfection capacity of distributed water [1,2]. Studying the chlorine decay reflects its importance in engineering and industrial sciences [3]. In this paper, we propose a generalized Received April 18th, 2021; accepted May 25th, 2021; published June 24th, 2021. 2010 Mathematics Subject Classification. 35R11. Key words and phrases. fractional partial differential equation; Bessel function; boundary value problem; exact solution; Laplace transform. ©2021 Authors retain the copyrights of their papers, and all open access articles are distributed under the terms of the Creative Commons Attribution License. 604 https://doi.org/10.28924/2291-8639 https://doi.org/10.28924/2291-8639-19-2021-604 Int. J. Anal. Appl. 19 (4) (2021) 605 model of the chlorine transport in pipes. The standard model was formulated by Biswas et. al [4]. The proposed model, in dimensionless form, is governed by the fractional partial differential equation (FPDE): (1) ∂αu ∂xα = A0 r ∂ ∂r ( 1 r ∂u ∂r ) −A1u, α ∈ (0, 1], under the boundary conditions (BCs): u(0,r) = 1, 0 ≤ r ≤ 1,(2) ∂ ∂r u(x, 0) = 0, 0 ≤ x ≤ 1,(3) ∂ ∂r u(x, 1) + A2u(x, 1) = 0, 0 ≤ x ≤ 1,(4) where u(x,r) is the chlorine concentration, α is the order of the fractional derivative in Caputo sense. The dimensionless parameters A0, A1 and A2 are related to the chlorine decay. The parameter A0 stands for the radial diffusion. It depends on the pipe length, the effective diffusivity of chlorine, and the flow rate throughout the system. In addition, A1 depends on the reactivity of chlorine with species such as viable cells or chemical compounds in the bulk liquid phase and on the residence time in the system. The parameter A2 reflects the wall consumption and depends on the wall consumption rate V ∗ d , the pipe radius r ∗ 0 and the effective diffusivity of chlorine D, where A2 = V ∗ d r ∗ 0/D. The objective of this paper is to apply the Laplace transform (LT) method to obtain the exact solution of the system (1)-(4). The LT is a well-known method for solving ordinary differential equations (ODEs) and fractional differential equations (PDEs). Such LT method has been successfully applied on several models such as diffusions [5], heat transfer of nanofluids suspended with carbon-nanotubes [6], singular boundary value problems (SBVPs) related to fluid flow of carbon-nanotubes [7,8], and the MHD Marangoni convection over a flat plate [9]. Furthermore, the LT was successfully applied to solve the Ambartsumian delay equation [10]. Moreover, one can find in Refs. [11-24] other interesting applications of the LT. Int. J. Anal. Appl. 19 (4) (2021) 606 2. The LT method The application of LT, L(·), on both sides of Eq. (1) gives (5) L ( ∂αu ∂xα ) = L ( A0 r ∂ ∂r ( 1 r ∂u(x,r) ∂r )) −L (A1u(x,r)) , or (6) sαU(s,r) −sα−1u(0,r) = A0 r d dr ( 1 r dU(s,r) dr ) −A1U(s,r), From the BC (2) and Eq. (6), we have (7) d2U(s,r) dr2 + 1 r dU(s,r) dr − ( sα + A1 A0 ) U(s,r) = − sα−1 A0 . The solution of Eq. (8) is (8) U(s,r) = ρ1J0 ( i √ sα + A1 A0 r ) + ρ2Y0 ( i √ sα + A1 A0 r ) + sα−1 sα + A1 , where i = √ −1. Besides, J0 (·) and Y0(·) are Bessel functions and ρ1 and ρ2 are unknown constants. Physically, u(x,r) is bounded at r = 0, hence ρ2 must be vanishes and therefore, (9) U(s,r) = ρ1J0 ( i √ sα + A1 A0 r ) + sα−1 sα + A1 , and (10) dU(s,r) dr = −ρ1i √ sα + A1 A0 J1 ( i √ sα + A1 A0 r ) , where J′0(λr) = −λJ1(λr). The BC (4) gives (11) dU(s, 1) dr + A2U(s, 1) = 0. From Eqs. (9), (10) and (11), we obtain ρ1 as (12) ρ1 = − A2s α−1 (sα + A1) [ A2J0 ( i √ sα+A1 A0 ) − i √ sα+A1 A0 J1 ( i √ sα+A1 A0 )]. Int. J. Anal. Appl. 19 (4) (2021) 607 Substituting (12) into (9), yields (13) U(s,r) = − A2s α−1J0 ( i √ sα+A1 A0 r ) (sα + A1) [ A2J0 ( i √ sα+A1 A0 ) − i √ sα+A1 A0 J1 ( i √ sα−1+A1 A0 )] + sα−1 sα + A1 , which can be written as (14) U(s,r) = −A2H(s,r) + sα−1 sα + A1 , where (15) H(s,r) = sα−1J0 ( i √ sα+A1 A0 r ) (sα + A1) [ A2J0 ( i √ sα+A1 A0 ) − i √ sα+A1 A0 J1 ( i √ sα+A1 A0 )]. Applying the inverse LT on Eq. (14) leads to (16) u(x,r) = −A2h(x,r) + Eα (−A1xα) , where h(x,r) is the inverse LT of H(s,r) and (17) h(x,r) = L−1   sα−1J0 ( i √ sα+A1 A0 r ) (sα + A1) [ A2J0 ( i √ sα+A1 A0 ) − i √ sα+A1 A0 J1 ( i √ sα+A1 A0 )]   . 3. The exact solution In order to find the exact solution, it should be first evaluate the inverse LT of (17). We observe from Eq. (15) or Eq. (17) that the denominator has simple poles at sα = −A1 and i √ sα+A1 A0 = λ1,λ2, . . .λn, . . . Hence, we have simple poles at s = (−A1)1/α and s = (−A1 −A0λ2n)1/α, n = 1, 2, 3, . . . , where λn are the roots of (18) A2J0 (λn) −λnJ1 (λn) = 0. Int. J. Anal. Appl. 19 (4) (2021) 608 So, h(x,r) can be evaluated by applying theorem 1, in appendix A, by calculating the residues of esxH(s,r) at s = (−A1)1/α and s = (−A1 −A0λ2n)1/α, and then by taking their sum. At s = (−A1)1/α, we have (Res esxH)s =(−A1)1/α = lim s→(−A1)1/α ( s− (−A1)1/α ) esxsα−1J0 ( i √ sα+A1 A0 r ) (sα + A1) [ A2J0 ( i √ sα+A1 A0 ) − i √ sα+A1 A0 J1 ( i √ sα+A1 A0 )], = e(−A1) 1/αx lim s→(−A1)1/α sα−1J0 ( i √ sα+A1 A0 r ) A2J0 ( i √ sα+A1 A0 ) − i √ sα+A1 A0 J1 ( i √ sα+A1 A0 ) × lim s→(−A1)1/α s− (−A1)1/α sα + A1 , = e(−A1) 1/αx ( (−A1)(α−1)/αJ0 (0) A2J0 (0) − 0 ) . lim s→(−A1)1/α 1 αsα−1 , = e(−A1) 1/αx αA2 , where J0 (0) = 1.(19) At s = (−A1 −A0λ2n)1/α, we have (Res esxH)s=(−A1−A0λ2n)1/α = lim s→(−A1−A0λ2n)1/α ( s− (−A1 −A0λ2n)1/α ) esxsα−1J0 ( i √ sα+A1 A0 r ) (sα + A1) [ A2J0 ( i √ sα+A1 A0 ) − i √ sα+A1 A0 J1 ( i √ sα+A1 A0 )], or (Res esxH)s=(−A1−A0λ2n)1/α = lim s→(−A1−A0λ2n)1/α s− (−A1 −A0λ2n)1/α A2J0 ( i √ sα+A1 A0 ) − i √ sα+A1 A0 J1 ( i √ sα+A1 A0 ) × lim s→(−A1−A0λ2n)1/α esxsα−1J0 ( i √ sα+A1 A0 r ) sα + A1 , which can be written as (20) (Res esxH)s=(−A1−A0λ2n)1/α = e(−A1−A0λ 2 n) 1/αx ( −A1 −A0λ2n )(α−1)/α J0 (−λnr) −A0λ2n . lim s→(−A1−A0λ2n)1/α Q(s,r). Int. J. Anal. Appl. 19 (4) (2021) 609 Using the L’Hospital’s rule, we have lim s→(−A1−A0λ2n)1/α Q(s,r) = lims→(−A1−A0λ2n)1/α ( s− (−A1 −A0λ2n)1/α ) lims→(−A1−A0λ2n)1/α [ A2J0 ( i √ sα+A1 A0 ) − i √ sα+A1 A0 J1 ( i √ sα+A1 A0 )] = 0 0 , = lims→(−A1−A0λ2n)1/α d ds ( s + A1 + A0λ 2 n ) lims→(−A1−A0λ2n)1/α d ds [ A2J0 ( i √ sα+A1 A0 ) − i √ sα+A1 A0 J1 ( i √ sα+A1 A0 )], = 1 σ ,(21) where σ is defined by (22) σ = lim s→(−A1−A0λ2n)1/α d ds [ A2J0 ( i √ sα + A1 A0 ) − i √ sα + A1 A0 J1 ( i √ sα + A1 A0 )] . Assume that y = i √ sα+A1 A0 , then σ = lim s→(−A1−A0λ2n)1/α d ds [A2J0(y) −yJ1(y)] = lim s→(−A1−A0λ2n)1/α [ −A2J1(y) dy ds − y 2 (J0(y) −J2(y)) dy ds −J1(y) dy ds ] , = lim s→(−A1−A0λ2n)1/α dy ds [−A2J1(y) −yJ0(y)] ,(23) where the properties of Bessel functions are used, see appendix B. The magnitude dy ds is (24) dy ds = − αsα−1 2A0y . We also note at s = (−A1 −A0λ2n)1/α that (25) y = −λn, dy ds = α(−A1 −A0λ2n)(α−1)/α 2A0λn . Inserting Eqs. (25) into Eq. (23), noting that the functions J0 and J2 are even and J1 is odd, we find (26) σ = α(−A1 −A0λ2n)(α−1)/α 2A0λn (A2J1 (λn) + λnJ0 (λn)) . From (26) and (21), it then follows (27) lim s→(−A1−A0λ2n)1/α Q(s,r) = 2A0λn α(−A1 −A0λ2n)(α−1)/α (A2J1 (λn) + λnJ0 (λn)) . Int. J. Anal. Appl. 19 (4) (2021) 610 Substituting (27) into (20) and simplifying, gives (28) (Res esxH)s=(−A1−A0λ2n)1/α = − 2 α ∞∑ n=1 e(−A1−A0λ 2 n) 1/αxJ0 (λnr) λn [A2J1 (λn) + λnJ0 (λn)] . Hence, h(x,r) in Eq. (17) is given by h(x,r) = (Res esxH)s=(−A1)1/α + (Res e sxH)s=(−A1−A0λ2n)1/α , = e(−A1) 1/αx αA2 − 2 α ∞∑ n=1 e(−A1−A0λ 2 n) 1/αxJ0 (λnr) λn [A2J1 (λn) + λnJ0 (λn)] .(29) Inserting (29) into (14), and after simplifying, we obtain the solution u(x,r) as (30) u(x,r) = − 1 α e(−A1) 1/αx + 2 α ∞∑ n=1 A2e (−A1−A0λ2n) 1/αxJ0 (λnr) λn [A2J1 (λn) + λnJ0 (λn)] + Eα (−A1xα) , Eq. (18) implies (31) A2 = λnJ1 (λn) J0 (λn) . Substituting (31) into (30), yields (32) u(x,r) = − 1 α e(−A1) 1/αx + 2 α ∞∑ n=1 λnJ1 (λn) J0 (λnr) e (−A1−A0λ2n) 1/αx (A22 + λ 2 n) J 2 0 (λn) + Eα (−A1xα) , As α → 1, Eq. (32) reduces to (33) u(x,r) = −e−A1x + 2 ∞∑ n=1 λnJ1 (λn) J0 (λnr) e −(A1+A0λ2n)x (A22 + λ 2 n) J 2 0 (λn) + E1 (−A1x) , which can be simplified to (34) u(x,r) = 2 ∞∑ n=1 λnJ1 (λn) J0 (λnr) e −(A1+A0λ2n)x (A22 + λ 2 n) J 2 0 (λn) , where E1 (−A1x) = e−A1x. The solution (34) is identical to the same result obtained by Biswas et al. [4] for the chlorine decay model with classical partial derivative with respect to x. Int. J. Anal. Appl. 19 (4) (2021) 611 4. Results and Discussion According to Biswas et al. [4], the dimensionless cup-mixing average concentration is defined by (35) uav = 2 ∫ 1 0 u(x,r) rdr. Substituting (32) into (35), yields (36) uav = ( Eα (−A1xα) − 1 α e(−A1) 1/αx )∫ 1 0 2rdr + 4 ∞∑ n=1 λnJ1 (λn) e (−A1−A0λ2n) 1/αx (A22 + λ 2 n) J 2 0 (λn) ∫ 1 0 rJ0 (λnr) dr, or (37) uav = Eα (−A1xα) − 1 α e(−A1) 1/αx + 4 ∞∑ n=1 J21 (λn) (A22 + λ 2 n) J 2 0 (λn) e(−A1−A0λ 2 n) 1/α x. Implementing the relation (31) we have (38) uav = Eα (−A1xα) − 1 α e(−A1) 1/αx + 4 ∞∑ n=1 A22 λ2n (A 2 2 + λ 2 n) e(−A1−A0λ 2 n) 1/α x. If the pipe walls act as a perfect sink, i.e., V ∗0 →∞ or A2 →∞, then the cup-mixing average concentration is obtained from Eq. (38) by the limit: (39) uav = Eα (−A1xα) − 1 α e(−A1) 1/αx + 4 lim A2→∞ ( ∞∑ n=1 A22 λ2n (A 2 2 + λ 2 n) e(−A1−A0λ 2 n) 1/α x ) , which gives (40) uav = Eα (−A1xα) − 1 α e(−A1) 1/αx + ∞∑ n=1 4 λ2n e(−A1−A0λ 2 n) 1/α x, where λn’s are the roots of J0 (λn) = 0. Moreover, if V ∗ 0 → 0 or A2 → 0 (i.e., the pipe walls are inert and no chlorine consumption takes place at the walls), then u(x,r) in Eq. (19) reduces to (41) u(x,r) = Eα (−A1xα) , and accordingly, (42) uav = 2 ∫ 1 0 Eα (−A1xα) rdr = Eα (−A1xα) . Int. J. Anal. Appl. 19 (4) (2021) 612 Following Biswas et al. [4], we consider the first three terms of the series (38), hence, three roots λ1, λ2, and λ3 of Eq. (18) are to be used. Table 1 presents the three roots λ1, λ2 and λ3 of Eq. (18) at different values of A2 in the range 0.01 ≤ A2 < 1. The roots are calculated using the command “FindRoot” in MATHEMATICA. In Tables 2 and 3, the values of λ1, λ2 and λ3 are listed for selected values of A2 in the range 1 ≤ A2 < 10 and in the range 10 ≤ A2 < 1000, respectively. Table 1. The first three roots 1 , 2 , and 3 of Eq. (27) at different values of A2 in the range .101.0 2  A A2 1 2 3 0.01 0.141245 3.83431 7.01701 0.1 0.441682 3.85771 7.02983 0.2 0.616975 3.88351 7.04403 0.5 0.940771 3.95937 7.08638 Table 2. The first three roots 1 , 2 , and 3 of Eq. (27) at different values of A2 in the range .101 2  A A2 1 2 3 1 1.25578 4.07948 7.1558 2 1.59945 4.29096 7.28839 5 1.98981 4.71314 7.61771 Table 3. The first three roots 1 , 2 , and 3 of Eq. (27) at different values of A2 in the range .100010 2  A A2 1 2 3 10 2.1795 5.03321 7.95688 50 2.35724 5.4112 8.48399 100 2.3809 5.46521 8.56783 The curves of the cup-mixing average concentration uav are depicted in Figs. 1-4 versus A1, at the outlet x = 1 of a pipe, for several values of A0 and A2 when α = 1/3. Figure 1 indicates that the uav is a decreasing Int. J. Anal. Appl. 19 (4) (2021) 613 function in the parameter A1 in the absence of A2 (i.e., A2 = 0). However, the behavior of uav is different in the case A2 6= 0 where uav decreases in two subdomians of A1 and increases in a certain domain. This last conclusion can be also confirmed and seen in Figs. 2-4 for the curves of uav when A2 has a specified nonezero value, i.e., A2 doesn’t vanish. The influence of the fractional order α on the cup-mixing average concentration uav is displayed in Fig. 5. It can be seen from this figure that uav is a decreasing function in the full domain of the parameter A1 when α = 1 (classical derivative) of while uav is of different behavior when α = {1/3, 1/5, 1/7} (fractional derivative). The discussion above may give some lights about the modeling of chlorine decay in view of the fractional calculus. A2=0.0 A2=0.1 A2=0.3 A2=0.5 2 4 6 8 10 A1 -2.0 -1.5 -1.0 -0.5 0.5 1.0 uav Figure 1. The cup-mixing average concentration uav versus A1 at different values of A2 when α = 1/3 and A0= 1.4. Int. J. Anal. Appl. 19 (4) (2021) 614 A2=1 A2=3 A2=5 A2=7 2 4 6 8 10 A1 -1.0 -0.5 uav Figure 2. The cup-mixing average concentration uav versus A1 at different values of A2 when α = 1/3 and A0 = 1.4 × 10−3. A2=10 A2=20 A2=30 A2=40 2 4 6 8 10 A1 -1.5 -1.0 -0.5 uav Figure 3. The cup-mixing average concentration uav versus A1 at different values of A2 when α = 1/3 and A0 = 1.4 × 10−2. Int. J. Anal. Appl. 19 (4) (2021) 615 A2=100 A2=200 A2=300 A2=400 2 4 6 8 10 A1 -2.0 -1.5 -1.0 -0.5 uav Figure 4. The cup-mixing average concentration uav versus A1 at different values of A2 (higher values) when α = 1/3 and A0 = 1.4 × 10−2. α=1 α=1/3 α=1/5 α=1/7 2 4 6 8 10 A1 -6 -5 -4 -3 -2 -1 uav Figure 5. Influence of the fractional order α on the cup-mixing average concentration uav when A0 = 1.4 and A2 = 0.5. Int. J. Anal. Appl. 19 (4) (2021) 616 5. Conclusion A possible generalization of the transport model describing the chlorine concentration decay in pipes was analyzed. The exact solution of the generalized model was obtained using the LT and the method of residues. The obtained exact solutions reduced to the corresponding published solutions as the fractional order α tends to one. Analytical expression for the dimensionless cup-mixing average concentration was deduced. The effects of the impeded parameters on the dimensionless cup-mixing average concentration were discussed and analyzed. The results showed that the behavior of the dimensionless cup-mixing average concentration in view of the fractional calculus is completely different than its behavior using the classical calculus. Appendices: A. Residues method A basic theorem for obtaining the inverse LT using the method of residues is given below. Theorem 1: The inverse LT of a function H(s,r) using the method of residues is given by h(x,r) = Sum of residues of esxH(s,r) at all poles of H(s,r), see Ref. [25] for details. B. Properties of Bessel functions The Bessel functions J0(y), J1(y), and J2(y) are defined by J0(y) = ∞∑ k=0 (−1)k (k!) 2 (y 2 )2k ,(B.1) J1(y) = ∞∑ k=0 (−1)k k!(k + 1)! (y 2 )2k+1 ,(B.2) J2(y) = ∞∑ k=0 (−1)k k!(k + 2)! (y 2 )2k+2 ,(B.3) Int. J. Anal. Appl. 19 (4) (2021) 617 and satisfy the properties: d dy (J0(λy)) = −λJ1 (λy)(B.4) d dy (J1(λy)) = λ 2 (J0(λy) −J2(λy)) ,(B.5) yJ2(y) + yJ0(y) = 2J1(y).(B.6) Availability of data and materials: Not applicable. Funding: NA. Conflicts of Interest: The author(s) declare that there are no conflicts of interest regarding the publication of this paper. References [1] R.M. Clark, E.J. Read, J.C. Hoff, Analysis of Inactivation of Giardia Lamblia by Chlorine, J. Environ. Eng. 115 (1989), 80–90. [2] M.W. LeChevallier, C.D. Cawthon, R.G. Lee, Inactivation of biofilm bacteria, Appl. Environ. Microbiol. 54 (1988), 2492–2499. [3] B.F. Arnold, J.M. Colford, Treating Water With Chlorine at Point-of-Use to Improve Water Quality and Reduce Child Diarrhea in Developing Countries: A Systematic Review and Meta-Analysis, Amer. J. Trop. Med. Hyg. 76 (2007), 354–364. [4] P. Biswas, C. Lu, R.M. Clark, A model for chlorine concentration decay in pipes, Water Res. 27 (1993), 1715–1724. [5] J. Jakubowski, M. Wisniewolski, On matching diffusions, Laplace transforms and partial differential equations, Stoch. Proc. Appl. 125 (2015), 3663-3690. [6] A. Ebaid, M. Al Sharif, Application of Laplace transform for the exact effect of a magnetic field on heat transfer of carbon-nanotubes suspended nanofluids, Z. Naturforsch., A, 70 (2015), 471-475. [7] A. Ebaid, A.M. Wazwaz, E. Alali, B. Masaedeh, Hypergeometric Series Solution to a Class of Second-Order Boundary Value Problems via Laplace Transform with Applications to Nanouids, Commun. Theor. Phys. 67 (2017), 231. [8] A. Ebaid, E. Alali, H. Saleh, The exact solution of a class of boundary value problems with polynomial coefficients and its applications on nanofluids, J. Assoc. Arab Univ. Basi Appl. Sci. 24 (2017), 156-159. [9] S.M. Khaled, The exact effects of radiation and joule heating on magnetohydrodynamic Marangoni convection over a flat surface, Therm. Sci. 22 (2018), 63-72. [10] H.O. Bakodah, A. Ebaid, Exact solution of Ambartsumian delay differential equation and comparison with Daftardar-Gejji and Jafari approximate method, Mathematics, 6 (2018), 331. [11] S. Handibag, B.D. Karande, Laplace Substitution Method for Solving Partial Differential Equations Involving Mixed Partial Derivatives, Int. J. Comput. Eng. Res. 2 (2012), 1049-1052. Int. J. Anal. Appl. 19 (4) (2021) 618 [12] N. Dogan, Solution of the system of ordinary differential equations by combined Laplace transform-Adomian decomposition method, Math. Comput. Appl. 17 (2012), 203-211. [13] P. Rai, Application of Laplace Transforms to Solve ODE using MATLAB, J. Inform. Math. Sci. 7 (2015), 93-97. [14] S.S. Handibag, B.D. Karande, Laplace Substitution Method for nth Order Linear and Non-Linear PDE’s Involving Mixed Partial Derivatives, Int. Res. J. Eng. Technol. 2 (2015), 378-388. [15] A.A. Alshikh, M.M.A. Mahgob, A Comparative Study Between Laplace Transform and Two New Integrals “ELzaki” Transform and “Aboodh” Transform, Pure Appl. Math. J. 5 (2016), 145-150. [16] A. Atangana, B.S.T. Alkaltani, A novel double integral transform and its applications, J. Nonlinear Sci. Appl. 9 (2016), 424-434. [17] X. Lianga, F. Gao, Y.-N. Gao, X.-J. Yang, Applications of a novel integral transform to partial differential equations, J. Nonlinear Sci. Appl. 10 (2017), 528-534. [18] P.V. Pavani, U.L. Priya, B.A. Reddy, Solving Differential Equations by using Laplace Transforms, Int. J. Res. Anal. Rev. 5 (2018), 1796-1799. [19] B.M. Faraj, F.W. Ahmed, On the MATLAB technique by using Laplace transform for solving second order ODE with initial conditions exactly, Matrix Sci. Math. 3 (2019), 8-10. [20] A. Mousa, T.M. Elzaki, Solution of Volterra Integro-Differential Equations by Triple Laplace Transform, Irish Interdiscip. J. Sci. Res. 3 (2019), 67-72. [21] R.R. Dhunde, G.L. Waghmare, Double Laplace iterative method for solving nonlinear partial differential equations, New Trends Math. Sci. 7 (2019), 138-149. [22] D. Ziane, M.H. Cherif, C. Cattani, K. Belghaba, Yang-Laplace Decomposition Method for Nonlinear System of Local Fractional Partial Differential Equations, Appl. Math. Nonlinear Sci. 4 (2019), 489-502. [23] S. Mastoi, W.A.M. Othman, N. Kumaresan, Randomly generated grids and Laplace Transform for partial differential equations, Int. J. Disaster Recovery Bus. Contin. 11 (2020), 1694-1702. [24] H. Zhang, M. Nadeem, A. Rauf, Z. Guo Hui, A novel approach for the analytical solution of nonlinear time-fractional differential equations, Int. J. Numer. Meth. Heat Fluid Flow, 31 (2021) 1069–1084. [25] M.R. Spiegel, Laplace transforms, McGraw-Hill. Inc., New York, 1965. 1. Introduction 2. The LT method 3. The exact solution 4. Results and Discussion 5. Conclusion Appendices A. Residues method B. Properties of Bessel functions References