Boundary-domain integral method and homotopy analysis method for systems of nonlinear boundary value problems in environmental engineering

Abstract In this paper we present the development of numerical methods for two problems in chemical engineering: adsorption of carbon dioxide into phenyl glycidyl ether and diffusion and reaction processes in a porous catalyst. Both problems are modelled by systems of two coupled nonlinear ordinary differential equations. Dirichlet and Neumann type boundary conditions are considered. We develop the standard homotopy analysis method and the boundary-domain element method to solve for the unknown steady-state concentrations of reactants. The validity of the proposed methods is checked and demonstrated by numerical examples. Convergence rates are compared and the advantages and drawbacks of the proposed methods are discussed and compared to the Mathematica NDSOLVE solver. We show, that the proposed homotopy analysis method features exponential convergence rate and is highly accurate and efficient. As low as 6 terms are needed to reach error norms of 10−5. Furthermore, the boundary-domain integral approach is also capable of solving the same problem very efficiently, using a very coarse computational mesh (21 nodes).


Introduction
Numerical experiments and simulations (AbdAlmjeed, 2017;Aswhad, 2011) are frequently used to model natural and engineering phenomena in the fields of physics, chemistry, biology and environmental sciences.The phenomena are usually modelled by a system of equations with appropriate boundary and initial conditions.One encounters ordinary or partial differential equations, which may be linear or nonlinear.Unfortunately, most of these equations do not have an analytical solution or closed form, which could be represented by wellknown elementary functions.Therefore, there is always a demand to develop reliable and efficient mathematical methods to obtain either a numerical or an approximate solution of these problems.
Carbon dioxide is one of the most important gases in nature (Subramaniam, Krishnaperumal, & Lakshmanan, 2012) as it is used in photosynthesis.Furthermore, several engineering systems also make use of carbon dioxide, such as for extinguishing fires, for powering of pneumatic systems, for the production of carbonated soft drinks and to remove caffeine from coffee, precipitation processes (Choe, Oh, Kim, & Park, 2010;Park, Park, Kim, Park, & Oh, 2004;Subramaniam et al., 2012;Zhao, Han, Jakobsson, Louhi-Kultanen, & Alopaeus, 2016).Ever since excess carbon dioxide has been found to cause global warming, its chemical fixation is an important research area in the field of environmental engineering (Chatterjee, Rayalu, Kolev, & Krupadam, 2016).
We model the adsorption of carbon dioxide (CO 2 ) into phenyl glycidyl ether (PGE) by relating the concentrations of both species using a system of two coupled nonlinear ordinary differential equations.This process has been studied by Park et al. (2004) and Choe, Oh, et al. (2010).The approximate analytical expression for the steady-state concentrations of carbon dioxide and phenyl glycidyl ether have been presented by Subramaniam et al. (2012) using the Adomian's decomposition method (ADM).The drawback of this approach is the need to calculate the Adomian polynomials, which is computationally demanding.The methods proposed in this paper, do not have such a drawback.
The Lane-Emden equation is used to model many phenomena in astrophysics and mathematical physics.The description of diffusion and reaction processes in a porous catalyst, which we consider in this paper, is one of them.Other examples include: the thermal behaviour of a spherical cloud of gas, the theory of stellar structure and the theory of thermionic currents.The systems of Lane-Emden equations arise when modelling many other physical phenomena as well (Khojasteh Salkuyeh & Tavakoli, 2016;Sun, Liu, & Keith, 2004;Wazwaz, 2001;Wazwaz & Rach, 2011).Wazwaz et al. ( Wazwaz, Rach, & Duan, 2013;Wazwaz, Rach, & Duan, 2014) used the Adomian decomposition method for solving the Volterra integral form of the Lane-Emden equations with initial values and boundary conditions.
The main advantage of the standard HAM is that in most cases its solution series converges very rapidly.Typically, only few iterations lead to a very accurate solution.HAM is a universal method and can be used to solve different types of nonlinear equations.Due to homotopy, HAM gives great freedom in the selection of initial approximations and auxiliary linear operators.The original HAM has been upgraded by introducing an auxiliary parameter h, which enables control of the convergence of the series.With this modification, the HAM is more efficient than the ADM (Adomian, 1994).
The domain-decomposition based boundarydomain integral method (BDIM) has been successfully used to solve various fluid flow and heat transfer problems (Ravnik, Skerget, & Zuni c, 2008;Ravnik, Skerget, & Zunic, 2009).It is based on the boundary element method (Haghighat & Binesh, 2014).In this work, we rework the method for the solution of two chemical engineering problems as an alternative to the homotopy analysis method.The proposed method is second order accurate and can, due to the inclusion of function derivatives into the system of equations, be a good alternative for problems, which feature sharp function profiles and large gradients.
In this work, HAM and BDIM are developed for the solution of two chemical engineering problems.Both methods are derived and adapted to work for the two chosen chemical engineering problems.The advantage of the proposed HAM method in this paper is that it features exponential convergence rate and it is not hindered by the nonlinearity of the problem, which is the case with the Adomian decomposition method, which was proposed by other authors.In the HAM method proposed in this paper, there is no need for calculation of the Adomian polynomials.Additionally, we present the boundary-domain integral approach and show that it is also capable of solving the same problem very efficiently, using a very coarse computational meshes.
The paper is organized as follows.The governing equations describing the chemical processes are presented in section 2. Boundary element based solution is shown in section 3. Homotopy analysis method is developed in section 4. Simulations and discussion is given in section 5, which is followed by conclusions in section 6.

Diffusion and reaction processes in porous catalyst
Diffusion and reaction processes in porous catalyst have been a topic of research of several authors.Sun et al. (2004) consider a planar nonlinear model of diffusion and reaction in porous catalysts using a decomposition method.Similarly, a spherical pellet is considered by Shi-Bin, Yan-Ping, and Scott (2003).Rach, Duan, and Wazwaz (2014) proposed the Adomian decomposition method for the solution of the catalytic diffusion reactions.
Considering mass balance on a differential volume element in a spherical porous medium we may write the following transient diffusion-reaction equation (Bird, Stewart, & Lightfoot, 2007) where c is the reactant concentration, t is time, D the effective diffusivity of the reactant and R is the rate of reaction per unit volume.We assume that the process is steady state and that the diffusivity is constant (temperature variations are small), thus we may simplify Let r measure the distance from the centre of the catalyst pellet (r 0 ¼ 0) to the surface (r 0 ¼ r 0 ).The process is subjected to the following boundary conditions where c s is the reactant concentration at the surface of the porous medium.Now we consider a nth order irreversible reaction and model the reaction term using power-law kinetics as R ¼ k v c0 n , where the reaction constant is k v and it is in general temperature dependent.
Introducing nondimenzionaliton r ¼ r 0 =r 0 and c ¼ c 0 =c 0 and considering spherical symmetry we get where / is the Thiele modulus / 2 ¼ r 2 0 kc nÀ1 0 =D: The derivation of the above expression is based on the validity of Fick's law, i.e., concentration independent diffusivity.
In multicomponent mixture the transport of species must be described by Maxwell-Stefan equations.The application of these to mass transport in porous medium leads to the Dusty Fluid Model (Higler, Krishna, & Taylor, 2000).Applying these, we obtain the following catalyst diffusion model (Flockerzi & Sundmacher, 2011) where D is a constant matrix of diffusion coefficients, N is the stoichiometry matrix of the considered reaction network and R is a vector of reaction kinetics expressions.Since chemical reactions are mostly based on nono-molecular or bi-molecular events, the rate expressions are linear, bilinear or quadratic expressions.For this kind of reaction kinetics, we can identify a sub-vector u containing concentrations of species which appear in linear and bilinear forms and another vector v containing concentrations of species which also appear in quadratic terms.An example of such a case is acid-catalysed dimerization of C4-olefines (Talwalkar et al., 2007).This class of catalyst diffusion problems can be reformulated as a 2 D Lane-Emden boundary value problem (Flockerzi & Sundmacher, 2011): in the domain r 2 ½0, 1 having the following boundary conditions Here, b 1 and b 2 are positive constants and k is a nonzero matrix.

Absorption of carbon dioxide into phenyl glycidyl ether
Due to global warming, chemical fixation of carbon dioxide has become an important research topic and adsorption is on the processes capable to address this issue.The problem of absorption of carbon dioxide into phenyl glycidyl ether has been investigated by several authors.Choe, Oh, et al. (2010) experimentally studied absorption of carbon dioxide with the aid of a catalyst.(Duan, Rach, & Wazwaz, 2015;Duan & Rach, 2011) proposed the use of Adomian decomposition method for the simulation of the absorption of carbon dioxide into phenyl glycidyl ether.Theoretical analysis of mass transfer was given by Subramaniam et al. (2012).
The reaction consists of two consecutive steps: a reversible reaction between PGE and a catalyst THA-CP-MS41 (QX) to form an intermediate complex C 1 and a reaction between the complex and carbon dioxide to form the catalyst and a five-membered cyclic carbonate C. The reactions are k 1 , k 2 and k 3 are reaction rates.Governing equations (Bird et al., 2007) for mass balances of CO 2 and PGE obtained using film theory are (Choe, Oh, et al., 2010): where S t is the area of the catalyst, c A and c B are concentrations of CO 2 and PGE, K 1 is the reaction equilibrium constant, D A and D B are diffusivity of CO 2 and PGE, r A, cons is the steady state CO 2 chemical reaction rate and z is the distance.The boundary conditions are: These equations may be normalized and rewritten in terms of the non-dimensional distance r ¼ z=L : with boundary conditions: The a 1 , a 2 , b 1 and b 2 are normalized parameters and u and v are dimensionless concentrations of CO 2 and PGE.The ratio between CO 2 concentrations at both sides of the domain is k ¼ c A0 =c AL :

Solution by boundary-domain integral method
Due to the diffusive nature of the considered transport processes, both boundary value problems, defined by Equations ( 6)-( 8) and ( 12)-( 14) feature a second order derivative, which describes diffusion at constant diffusivity.Mixed Dirichlet/Neumann boundary conditions are known for both cases.For a general 3 D setting, all equations may be written as where u is the unknown function and f is a nonlinear forcing term.We consider Dirichlet and/or Neumann type boundary conditions  15) into its integral representation (Wrobel, 2002): where u ?ðr, nÞ ¼ 1=4pjrÀnj is the fundamental solution for the diffusion operator, The free coefficient cðnÞ is determined using solid surface angle at the source point position.The domain integral term in Equation ( 17) results from the source term on the right hand side of Equation ( 15).
In order to write a discrete version of Equation ( 17) we must interpolate the unknown function and its flux over boundary and domain elements.We employ 27 noded hexahedral domain elements, which enable quadratic interpolation of function in the domain using shape functions U i : Flux is interpolated on the boundary only.We use four noded discontinuous linear interpolation functions / i : Interpolation of function at boundary elements is performed using nine nodes and quadratic shape functions u i : Finally, the discrete version of Equation ( 17 where b lists boundary elements, e domain elements and i shape functions.The integrals of Equation ( 18) depend solely on the mesh geometry and source point location.They are calculated in advance and reused during the iterative solution process.The Gaussian quadrature algorithm is used to calculate the integrals.We use indirect calculation of cðnÞ using a known solution of the rigid body movement.The accuracy of integration is checked using nodes with known value of cðnÞ: A collocation scheme is used to write a system of linear equations for the unknown values of function and flux.The source point n is placed into boundary nodes and the following matrices of integrals are thus calculated: Let the square brackets denote matrices and curly brackets vectors of nodal values of functions.Using this notation the discrete version of Equation ( 18) is As the boundary-domain integral method requires discretization of the domain and since the matrix of domain integrals ½B is full, we avoid excessive storage and computational time usage by using a domain decomposition technique (Ram sak, Skerget, Hriber sek, & Zuni c, 2005) Domain decomposition leads to a sparse system of equations.In this work, we consider the subdomains to be domain mesh elements.Connection between subdomains is obtained by considering that the function and the flux must be continuous across subdomains boundaries.The described procedure leads to a sparse and over-determined system of linear equations.We use (Paige & Saunders, 1982) least-squares solver with diagonal preconditioning to find the solution.
As the problems considered in this paper are nonlinear (the forcing function f in Eq. ( 15) depends on u), we set up an iteration procedure during which, we estimate f using values of u in the previous iteration.Underrelaxation of 0.1 was used to achieve convergence.Furthermore, as the problems considered are 1 D and the BDIM method is written in 3 D, we used appropriate (zero flux) boundary conditions on the side walls.

Solution by homotopy analysis method
To verify the basic idea of HAM for a nonlinear system of ordinary differential equations (Bataineh, Noorani, & Hashim, 2009;Gabriel, 2016;Liao, 2004;Wazwaz, 2009), let us consider this problem: where N i are the nonlinear operators, x is the independent variable and u i ðxÞ are the unknown functions.To clarify this method simplify, we will ignore the conditions.Liao (2004) and Wazwaz (2009) configured the so-called zeroth-order deformation equation where q 2 ½0, 1 represents an embedding parameter, h i represent the nonzero auxiliary functions, L is the auxiliary linear operator, u i, 0 ðxÞ represent the initial guesses of the function u i ðxÞ and h i ðx; qÞ are unknown functions.The advantage of HAM is that there is great freedom in the selection of the auxiliary objects, by means of h i and L.
It is obvious that when putting the parameter q ¼ 0 and q ¼ 1, Equation ( 22) will become h i ðx; 0Þ ¼ u i, 0 ðxÞ and h i ðx; 1Þ ¼ u i ðxÞ, (23) respectively.Therefore, when q increases from 0 to 1, the solution h i ðx; qÞ will vary from the initial guess u i, 0 ðxÞ to the solution u i ðxÞ: The expansion of h i ðx; qÞ in Taylor series with respect to q will be: where The auxiliary parameters h i are the important to ensure the convergence of the series (Liao, 1992).So, if the auxiliary parameters h i , the auxiliary linear operator and the initial guesses u i, 0 ðxÞ, are appropriately chosen, then the series (Liao, 1992) will converge at q ¼ 1.Thus must be one of the original nonlinear equation's solutions, which is proved by Liao (2004) and Wazwaz (2009).If h i ¼ À1, then (Khojasteh Salkuyeh & Tavakoli, 2016) will become This zeroth-order deformation equation is used in the homotopy perturbation method (HPM) [Rach et al., 2014).In accordance with the relationship (Liao, 2003), the governing equation could be derived from the zeroth-order deformation equation ( 22).Let us define a vector ũi, n ðxÞ ¼ fu i, 0 ðxÞ, u i, 1 ðxÞ, u i, 2 ðxÞ, u i, n ðxÞg, i ¼ 1, 2, , n: By differentiating equation ( 27) with respect to the embedding parameter q for m times, then substituting q ¼ 0 and dividing them by m!, the mth-order deformation equation is derived as: where and Implementing L À1 (an integral operator) on Equation ( 29) in both sides, we have in iterative way, it is easy to get u i, m for m ! 1 and for mth-order the following will be obtained An accurate approximate function for the general differential equation ( 21) will be obtained when j !þ1, if Equation ( 21) possesses a unique solution, then the HAM will present this unique solution, but is Equation ( 21) hasn't a unique solution then this method will provide a solution among another active solutions.In the following, each problem will be reviewed, and then we will apply the HAM on each of them.

The HAM for solving the Lane-Emden boundary value problem
The boundary value problem for the coupled equations of Lane-Emden with the product nonlinearities (Flockerzi & Sundmacher, 2011;Rach et al., 2014) are: with the following boundary conditions The above boundary value problem represents a model of catalytic diffusion reactions (Flockerzi & Sundmacher, 2011).The parameters of the system: b 1 , b 2 , k 11 , k 12 , k 21 and k 22 can be selected for the actual chemical reactions.Flockerzi and Sundmacher (2011) assumed that in the qualitative analysis for the solutions.In the following, we produce an approximate solution for the above nonlinear problem using the HAM (Bataineh et al., 2009;Gabriel, 2016;Liao, 2004;Wazwaz, 2009).In order to implement this method to the current system, we will solve the system introduced by Rach et al. (2014) which is used in the modified ADM.So, the system above will be reduced to the following form of two-coupled nonlinear Fredholm-Volterra integral equations To implement the HAM to the nonlinear system above, the basic idea of HAM will be applied.So, the zeroth approximations are The linear operators for Equations ( 37) and ( 38), respectively will be Now, defining a system of nonlinear operators in the following forms: According to the basic idea of HAM, the mthorder deformation equations will be where, for m ! 1, we have We get the following , and so on.When h ¼ À 1, we obtain Also, we get and so on.Proceeding in this way, the components were also calculated but for brevity, they are not listed here.
The used approximate functions of solutions will be as in Rach et al. (2014), in the following formulas: The used convenient error remainder functions will be also as in Rach et al. (2014), in the forms with the maximal error remainder parameters forms

The HAM for solving the absorption of CO2 by the PGE boundary value problem
A model for steady state concentrations of carbon dioxide CO2 and PGE has been proposed as in Duan et al. (2015) and Subramaniam et al. (2012).The following system of two second order nonlinear ordinary differential equations representing the two chemicals must be solved: u 00 ðxÞ ¼ a 1 f ðuðxÞ, vðxÞÞ, (49) v 00 ðxÞ ¼ a 2 f ðuðxÞ, vðxÞÞ, where with the following boundary conditions uð0Þ ¼ 1, uð1Þ ¼ k, v 0 ð0Þ ¼ 0 and vð1Þ ¼ 1: Here u(x) stands for the CO2 concentration and v(x) is the PGE concentration.The normalized system parameters are a 1 , a 2 , b 1 and b 2 : The carbon dioxide concentration at catalyst surface is denoted by k Finally, the dimensionless distance from the centre is x (Duan et al., 2015).Adomian (1994), Duan et al. (2015), Biazar, Babolian, and Islam (2004), have used the Adomian decomposition method to develop a solution of this problem.They transform the governing equations into a system of coupled integral equations.Due to the nonlinear nature of the problem, Adomian polynomials have to be calculated, which requires large computational resources.
In the following, we present an approximate solution of this nonlinear problem using the HAM (Bataineh et al., 2009;Gabriel, 2016;Liao, 2004;Wazwaz, 2009).In order to implement this method to the current system, we will solve the following system given in Flockerzi and Sundmacher (2011).So, By implementing the HAM to the nonlinear system above, the selected zeroth approximations are Consecutively we obtain the following terms by HAM: ! and so on.When h ¼ À 1, we get and so on.Proceeding in this way, the components were also calculated but for brevity, they are not listed here.
The used approximate functions of solutions will be as in (Wazwaz, 2005).The used convenient error remainder functions will be also as in Duan et al. (2015), in the forms and the maximal error remainder parameters are calculated as in equation (Wazwaz & Rach, 2011).

Diffusion and reaction processes in porous catalyst
In order to validate the proposed HAM method and examine the HAM convergence properties, we solve the problem of diffusion and reaction processes in porous catalyst using the parameter values suggested by Wazwaz et al. (2014) , and k 22 ¼ 1: We will evaluate the approximate solutions, the error remainder and the maximal error remainder.Furthermore, MER 1, n and MER 2, n are calculated, so that the values of n are increasing from 1 to 10.In the left panel of Figure 1 we display the convergence curve for MER 1, n : The convergence of MER 2, n is similar.

Absorption of carbon dioxide into phenyl glycidyl ether
In order to validate the proposed HAM method, we solve the problem of absorption of carbon dioxide into phenyl glycidyl ether using the parameter values suggested by Duan et al. (2015): and k ¼ 0:5: We will evaluate the approximate solutions, the error remainder and the maximal error remainder.Furthermore, MER 1, n and MER 2, n are calculated, so that the values of n are increasing from 1 to 10.In the right panel of Figure 1 the convergence curve for MER 1, n : The convergence of MER 2, n is similar.

Solution of the catalytic diffusion problem
In Figure 2 we compare concentration profiles for HAM and BDIM.The parameters used were: b Looking at HAM results, we observe that solution obtained using n 3 are different from the accurate n ¼ 10 solution, while for n ! 4 the convergence profiles do not differ significantly.The BDIM results have been obtained using different computational grids.Grids having 21 to 201 equidistantly places nodes were used.Due to the simple nature of the problem (diffusion only) and due to the use of quadratic interpolation, we detect virtually no difference between results obtained on different meshes.

Solution of absorption of carbon dioxide into phenyl glycidyl ether
Solution profiles for the absorption of carbon dioxide into phenyl glycidyl ether are shown in Figure 3.The   following parameters were used: a 1 ¼ 1, a 2 ¼ 2, b 1 ¼ 1, b 2 ¼ 3, k ¼ 0:5: In this case, only the n ¼ 1 HAM solution show significant difference from the converged n ¼ 10 solution, while all others (n ! 2) are very close to each other.Similarly to the catalytic diffusion problem, the BDIM is able to capture the physics of the phenomena using a very coarse mesh (21 nodes).

Comparison to Mathematica NDSOLVE solver
We solved both problems using Mathematica's NDSOLVE functions.Please refer to the appendix for the complete code.In order to asses the difference in solutions, we calculated the RMS difference between solution profiles as In Figure 4 we presented the RMS differences versus location x, while in Figure 5 average values are presented.We observe that the RMS error decreases with increasing n for the HAM method and with increasing number of nodes for the BDIM method.Looking at the RMS vs. x profile we see that the accuracy is approximately constant for all values of x.Accuracy of solution of the absorption problem is  higher than the accuracy for the catalytic diffusion problem for both methods.

Conclusions
In this work, the standard Homotopy analysis method and the boundary domain integral method has been developed for the solution of two process problems in chemical engineering: the simulation adsorption of carbon dioxide in phenyl glycidyl ether and simulation of catalytic diffusion reactions in porous catalysts.
We have shown, that the proposed HAM method features exponential convergence rate and is highly accurate and efficient.As low as 6 terms are needed to reach error norms of 10 À5 : On the contrary, due to the nonlinearity of problem, the application of the Adomian decomposition method, which was proposed by other authors needs the calculation of the Adomian polynomials.
Additionally, we have shown that the boundarydomain integral approach is also capable of solving the same problem very efficiently, using a very coarse computational mesh (21 nodes).
where C D is the Dirichlet part of the boundary and C N is the Neumann part.At the Neumann part of the boundary the flux qðrÞ is known and set to a known value of qðrÞ: At the Dirichlet part of the boundary, the function uðrÞ is set to a known value of uðrÞ: Using a source point n 2 C at the boundary we can recast the Poisson type equation (

Figure 1 .
Figure 1.Norms MER 1, n shown for HAM for n ¼ 1, . . ., 10 for the solution of the diffusion and reaction processes in porous catalyst problem (left) and the adsorption problem (right).Norms MER 2, n exhibit similar characteristics.

Figure 2 .
Figure 2. Comparison of solutions of the catalytic diffusion problem.u and v concentration profiles are presented for the HAM (left) and BDIM (right).

Figure 3 .
Figure 3.Comparison of solutions of the absorption problem.u and v concentration profiles are presented for the HAM (left) and BDIM (right).

Figure 4 .
Figure 4. RMS difference between u and v concentration profiles and the Mathematica NDSLOVE solver.Results for the catalytic diffusion problem are shown in the top panels and results of absorption of carbon dioxide are shown in the bottom panels.

Figure 5 .
Figure 5. Average RMS difference between u and v concentration profiles and the Mathematica NDSLOVE solver.Results for HAM are shown in the left panel and results for BDIM are shown in the right panel.Lines marked with CD are the catalytic diffusion problem and AB stands for absorption of carbon dioxide problem.