Vertical Sediment Concentration Distribution in High-Concentrated Flows:An Analytical Solution Using Homotopy Analysis Method

Manotosh Kumbhakar,Jitraj Saha,Koeli Ghoshal,Jitendra Kumar,and Vijay P.Singh

1Department of Mathematics,Indian Institute of Technology Kharagpur,Kharagpur 721302,India

2Department of Mathematics,National Institute of Technology,Tiruchirappalli,Tamil Nadu 620015,India

3Department of Biological and Agricultural Engineering&Zachry Department of Civil Engineering,Texas A&M University,College Station,TX 77843,USA

AbstractTransport of suspended sediment in open channel flow has an enormous impact on real life situations,viz.control and management of reservoir sedimentation,geomorphic evolution such as dunes,rivers,and coastlines etc.Transport entails advection and diffusion.Turbulent diffusion is governed by the concept of Fick’s law,which is based on the molecular diffusion theory,and the equation that represents the distribution of sediment concentration is the advection-diffusion equation.The study uses the existing governing equation which considers different phases for solid and fluid,and then couples the two phases.To deal with high-concentrated flow,sediment and turbulent diffusion coefficients are taken to be different from each other.The effect of hindered settling on sediment particles is incorporated in the governing equation,which makes the equation highly non-linear.This study derives an explicit closed-form analytical solution to the generalized one-dimensional diffusion equation representing the vertical sediment concentration distribution with an arbitrary turbulent diffusion coefficient pro file.The solution is obtained by Homotopy Analysis Method,which does not rely on the small parameters present in the equation.Finally,the solution is validated by comparing it with the implicit solution and the numerical solution.A relevant set of laboratory data is selected to check the applicability of the model,and a close agreement shows the potential of the model in the context of application to high-concentrated sediment-laden open channel flow.

Key words:open channel flow,diffusion equation,hunt equation,homotopy analysis method

1 Introduction

Transport of sediment in fluvial processes has been a long standing topic of research.[1−11]It plays a signi ficant role in the control and management of reservoir sedimentation,geomorphic evolution such as dunes,river,and coastlines,etc.that have immense importance in our day to day life.Better understanding of sediment transport mechanism in natural watercourses needs knowledge of open channel flow.Sediment transport in open channel flow occurs mainly in two different modes,namely,bedload and suspended load.[12]The part of the load that travels close to the bed is bed-load,while the other part travels in suspension as suspended load.In the study of suspended load,the spatial distribution of sediment concentration is of primary interest as it helps determine the discharge of suspended particles in watercourses,bed-load transport rate,and so on.The present study focuses on the vertical distribution of suspended sediment concentration in open channel turbulent flow.

Thegoverning equation forsedimentconcentration distribution,termed as advection-diffusion equation(ADE),has the origin in Fick’s law of diffusion.Rouse[7]was the pioneer in this field to provide an analytical expression for vertical concentration pro file by balancing the upward diffusion and the downward settling of sediment,and the solution provided by him is widely known as the Rouse equation.However,it has a major drawback that it cannot predict the sediment concentration values well near the channel bed or near the water surface,which might be due to the inclusion of Prandtl-von K´arm´an velocity pro file that does not model velocities well in sediment-laden flows.[13]Following the work of Rouse,[7]different modifications have been done by researchers to overcome those drawbacks by incorporating several turbulent mechanisms.[6,9,11,14−15]On the other hand,Hunt[3]analyzed the dynamics of particle suspension by treating the fluid phase and the solid phase separately,and then coupled the phases by decomposing the vertical velocity of the particle into the vertical flow velocity and the particle settling velocity in a quiescent fluid.The Hunt equation has more practical importance,as unlike the Rouse equation,it can be applied to a region where the volume occupied by the particles is appreciable.For simplicity,Hunt[3]assumed the turbulent diffusion coefficient to be the same as sediment diffusion coefficient,which is a poor approximation valid only in the case of small particles as mentioned in his work.Later on,different experimental investigations revealed the non-equality of these coefficients for different flow conditions,[16−18]and different models on these coefficients are reported in literature.[9,19]

Numerous studies have been undertaken to investigate the vertical distribution of sediment concentration starting from the Hunt equation.But most of the researchers took the turbulent and sediment diffusion coefficients to be equal limiting the application of their studies.Very few researchers[6,20−21]considered the two coefficients to be different,and solved the equation numerically as the traditional analytical techniques did not succeed due to the non-linearity of the diffusion equation.To deal with the non-linear ODE-PDEs,classical perturbation[22]and asymptotic techniques[23−24]are widely applied to obtain analytic approximations of nonlinear problems in science and engineering.Unfortunately,perturbation and asymptotic techniques are too strongly dependent upon small/large physical parameters present in the governing equation,and thus are often valid only for weakly nonlinear problems.

The present study aims to derive an explicit analytical solution of the generalized diffusion equation with arbitrary turbulent diffusion coefficient,by means of a novel unified method,known as Homotopy Analysis Method(HAM).Liao[25]developed HAM to solve non-linear differential equations,which employs the concept of homotopy from topology to generate a convergent series solution for nonlinear systems.It has been proved[26]that the method is a unified one,which logically contains Lyapunov’s small artificial parameter method,[27]Adomian decomposition method,[28]Homotopy perturbation method(HPM),[29]the δ-expansion method,[30]and the Euler transform[31]as special cases.In summary,HAM distinguishes itself from various other analytical methods in three respects:(i)it does not directly rely on small physical parameters present in the governing equation;(ii)it assures the convergence of a non-linear differential equation;and(iii)it has flexibility regarding the choice of base functions and the auxiliary linear operator of the homotopy.Since its inception,HAM has successfully been applied to different fields of science and engineering.[32−44]The derived analytical solution obtained via HAM in the present study is validated with the implicit solution and with the 4th order Runge-Kutta(RK)method.The model is also compared with relevant experimental data on vertical sediment concentration of high-concentrated flow.

2 Mathematical Model

A sediment-laden flow of depth h in a wide open channel,in which the time-averaged suspended sediment concentration decreases monotonically from a maximum value Caat the reference level y=a to zero at the water surface,is presented in Fig.1.The formulation of the generalized diffusion equation representing sediment concentration and the methodology to obtain the analytical solution is described in what follows.

Fig.1 Schematic diagram of sediment concentration distribution along a vertical in open channel turbulent flow over a sediment bed;within the reference level y≤a,sediments move as bed load,and above that as suspended load.

2.1 Governing Equation

In the case of steady-uniform flow where sediment concentration varies only in the vertical direction,the generalized governing equation can be written as:[3]

where C is the volumetric suspended sediment concentration(dimensionless),y is the vertical axis,ωsis the settling velocity of sediment particles,and εsand εtare sediment and turbulent diffusion coefficients,respectively.The ratio of εsto εsis denoted by a proportionality factor β.In general,β is assumed to be a constant over the entire depth of flow.Hunt[3]assumed β ≈ 1(i.e.,εs≈ εs)for simplicity.However,experimental results revealed that β≈1 for finer sediment particles only,and β<1 for coarse sediment particles.[16−18]On the other hand,using the experimental results reported by Coleman,[45]van Rijn[9]suggested β to be greater than 1,which includes the in fluence of centrifugal acceleration enforcing the sediments to move outwards of an eddy,and results in an increase in the effective mixing length of eddies.Thus,the assumption β≈1 may not be physically realistic.

Experimental observations revealed that in the case of flow with suspended sediment particles,the flow around the neighboring settling particles exhibits a larger drag as compared to the clear water flow,known as hindered settling effect.This effect reduces the settling velocity in a sediment-laden flow from that in a clear fluid.The widely used expression for ωswas proposed by Richardson and Zaki,[46]which reads as:

where ω0is the settling velocity of sediment in a pure quiescent fluid,and nHis an exponent,which is termed as exponent of reduction of settling velocity.Therefore,to incorporate the real physics associated with flows carrying sediment,Eq.(2)is inserted in the governing equation(1).

The expression for εs(y)or εt(y)should be available to determine the pro file of vertical concentration distribution,as can be noticed from Eq.(1).To that end,two different expressions,(namely,linear and parabolic)available in the literature,are employed in this study.They can be presented as:

where κ is the von-K´arm´an constant,u∗is the shear velocity,and h is the maximum flow depth.The turbulent diffusion coefficient is mainly derived from the Boussenisq hypothesis based on the analogy of Newton’s law of viscosity for turbulent flow.[47]The parabolic pro file was obtained by incorporating the well-known logarithmic pro file for stream-wise velocity and a linear pro file for Reynolds shear stress.It was shown that under steady-uniform and unidirectional open-channel flows,[9]parabolic pro file estimates experimental data better than does the linear profile.It is noted that different pro files for turbulent diffusion coefficient are needed to characterize different types of flow,and thus both are considered in the present study.

Before proceeding further,we rearrange Eq.(1)in nondimensional form as follows:

in which

and A is the well-known Rouse number,i.e.,A= ω0/βκu∗.Equation(4)is a first-order highly non-linear ordinary differential equation,which represents vertical distribution of suspended sediment concentration.The solution methodology for solving Eq.(4)is described in what follows.

2.2 Implicit Solution Obtained for Eq.(4)without and with Hindered Settling Effect

Without the effect of hindered settling,i.e.,nH=0,Eq.(4)can be rewritten as:

Using the boundary conditionand integrating Eq.(5),we obtain

Further,with the effect of hindered settling through exponent nH,Eq.(4)can be rearranged as:

The volumetric sediment concentration C,which is the ratio of volume of sediment to volume of sediment- fluid mixture,takes on a maximum value 0.6 very near the channel bed.[48]Furthermore,the vertical concentration pro file C monotonically decreases from the channel bed to the water surface.The present study deals with the suspended sediment concentration distribution which starts from a certain height above the channel bed,i.e.,from the reference level.Thus,the reference concentration Câis much smaller than 0.6.Using this characteristic of C,the following approximation can be made:

Equation(8)is quite accurate due to the smallness of C in the region[a,h].Using the approximation Eq.(8)in Eq.(7),and then integrating using the boundary condition,one can obtain the implicit solution in the following form:

Equations(6)and(9)represent implicit solutions for the vertical sediment concentration distribution without and with the effect of hindered settling,respectively.The right hand side of Eqs.(6)and(9)can easily be integrated with the substitution of linear and parabolic pro files for K(z)mentioned previously.

However,Eqs.(6)and(9)are only implicit solutions.Our objective is to find an explicit analytical solution to Eq.(4).Explicit functions have many advantages over the implicit ones,for example,one can describe the character-istics of the dependent variable explicitly dependent on the independent variable,and an application to that would be to find out the suspended load transport rate,which is the integration of the product of sediment concentration and velocity over the flow depth.The present study,therefore,shows the potentiality of a novel unified mathematical method,called,Homotopy Analysis Method(HAM),for handling non-linear differential equations arising in various field of science and engineering.

2.3 Solution Obtained by HAM

We suggest this unified method,HAM,not only for solving the equation modeled in this study but for any strong non-linear problem arising in the field of hydraulics,especially sediment transport phenomenon;therefore,for the sake of completeness,here we brie fly present the methodology(as proposed by Liao[26])associated with it in a convenient way.For that purpose,let us consider the original non-linear Eq.(4)in the following form:

where N is a non-linear operator,C(z)is the unknown function,and z is the independent variable.One can construct the homotopy function as follows:

where q∈[0,1]is the embedding-parameter,is a nonzero auxiliary parameter,L is a linear operator,and H(z)is a non-zero auxiliary function.The underlying idea behind HAM is that a continuous mapping is described to relate the solution C(z)and the unknown function Φ(z;q),with the aid of the embedding-parameter q.[25]Thus,the initial approximation C0(z)of the solution C(z)is so chosen that Φ(z;q)varies from C0(z)to C(z),as q varies from 0 to 1.Representing this mathematically,one can write Φ(z;0)=C0(z)and Φ(z;1)=C(z).We now define:

where(∂mΦ(z;q)/∂qm)|q=0is called the m-th order deformation derivative.Using Taylor’s theorem,one can expand Φ(z;q)in a power series of the embedding parameter q as follows:

Assume that N,H(z),C0(z),andare so properly chosen that the series Eq.(13)converges at q=1.Then,at q=1,the series becomes

Equation(14)provides the relationship between the initial approximation C0(z)and the exact solution C(z).However,the higher order approximations Cm(z)for m≥1 are still unknown to us.Liao[25]showed that the higher order terms can be obtained by differentiating the zeroth order deformation Eq.(11)(with H(Φ(z;q);q,,H)=0)m times with respect to q,and setting q=0,and finally dividing by m!,as follows:

where

Liao[26]proposed some generalized rule of solution expressions,rule of coefficient ergodicity,and rule of solution existence in order to choose the initial approximation,linear operator,and the auxiliary function.The importance of the auxiliary parameterlies in the fact that it provides us a convenient way to control and adjust the convergence region of the obtained series solution.

The term(1−C)nH+1presents in the governing equation involves a non-integer exponent,which makes the equation highly non-linear and may take longer computational time when applying HAM.Thus,before proceeding further,we use some accurate approximation for that term similar to Eq.(8).Using the approximation,the governing Eq.(4)is rewritten as follows:

For Eq.(18),in accordance with Liao’s rule of solution expression,a single-term linear operator is chosen in order for simplicity of the subsequent derivations:

The nonlinear operator is selected as:

As pointed out by Vajravelu and Van Gorder,[49]without any loss of generality,the auxiliary function H(z)can be taken as 1.According to their analysis,the selection of H(z)is not necessary if one selects a proper or optimal linear operator L.Following them,we also have taken the function to be 1 to avoid any unnecessary difficulty in the computation.The key is to find out Rm’s from Eq.(17),the first few terms are calculated as follows(“′” represents first-order derivative with respect to the independent variable z):

Proceeding like above,one can find a closed-form for Rm(m≥1)corresponding to Eq.(18)as follows:

Now,applying the inverse of the linear operator Eq.(19),and using the boundary condition C(z=â)=Câ,one can obtain the m-th term as:

Finally,the M-th order HAM-based approximation of C(z)can be obtained explicitly as:

It can be observed from Eq.(26)that one can obtain the iteration as per wish once the initial approximation C0(z)is known.For subsequent derivations,the initial approximation C0(z)=Câis selected to achieve the HAM-based series solution.

3 Results and Discussion

In this section,we validate the HAM-based solution given by Eq.(26)with the implicit solutions(Eqs.(6)and(9))as well as with 4-th order Runge-Kutta(RK)method for two different pro files of K(z).Finally,the HAM-based solution is applied to a laboratory dataset of vertical concentration pro file for high-concentrated flow to check the applicability of the model.

3.1 Selected Expressions for the Parameters

It is to be noted that in order to apply the derived analytical model(Eq.(26)),the Rouse number A,proportionality parameter β,and the exponent of reduction of fall velocity nHneed to be computed.Many researches[46,50−51]have reported that nHdepends mainly on the particle Reynolds number.Out of several formulae available,we selected the widely used expression of Richardson and Zaki.[46]According to them,nHdepends on the particle Reynolds number Repby the following relation:

in which Rep=(ω0/νf)dp,dpbeing the particle diameter and νfbeing the kinematic viscosity of clear fluid.Here,the clear water settling velocity ω0is calculated from the well-known expression given by Cheng,[52]which reads as:

where d∗is the dimensionless particle diameter,defined as,,g being the gravitational acceleration.

On the other hand,investigations on the determination of expression for β have been limited since the start.[9,19]Recently,Pal and Ghoshal[19]studied the nature of β by analyzing a large number of data sets of different researchers,and established two relations of β for dilute and non-dilute flows.Their study claims that β values do not only depend upon normalized settling velocity,but also on reference level and reference concentration.They developed the following expressions:

In this study,we have used the expression given by Eq.(29).The Rouse number A:= ω0/βκu∗is computed using the expressions for ω0and β given by Eqs.(28)and(29),respectively.The parameter κ is a universal constant having the value 0.41.Calculation of the reference level z=â and the corresponding concentrationis difficult in sediment transport studies because the formulae available in the literature lead to an incorrect estimation of vertical concentration values.To get rid of it,a few researchers[6,20]computedâ as the lowest available height reported in the observed data set and the concentration thereat as.In this study,the reference level and reference concentration are selected in the same way.

3.2 Validation of the HAM-Based Solution

It can be observed that the HAM-based series solution given by Eq.(26)is basically a function of z and.The study of Liao[26]asserts that the auxiliary parameterplays the key role in controlling the convergence of the series solution Eq.(26).Thus,a suitable choice ofwill lead towards the convergence of the series Eq.(26)to the exact solution for all z in its domain of definition(for detailed proof,see Appendix A).As reported in Liao,[26]the region of flatness of-curves determines the values offor which the HAM-based series CSUMconverges to the exact solution C(z).Similar observation also follows for the derivative serieswhich converges to C′(z).The mathematical justification behind this observation is given by Abbasbandy et al.,[32]However,such-curves do not produce the best value offor which the series solution rapidly converges.To that end,the exact square residual error∆mat the m-th order approximation has been introduced(Liao[53])which is defined byFor a given order of approximation,the optimal value ofis given by the minimum of∆m.But the higher the order of approximation the greater is the computational time needed to obtain the exact residual error∆m.Therefore,one can use the average residual error Emto decrease the computational time,which is given as Em=In Figs.2(a)and 2(b),the average residual error Emis plotted for different orders of approximation,considering K(z)as linear and parabolic,respectively.It can be observed from the figures that for each order of approximation,there exists an optimalfor which the error Emis minimum.

Fig.2Average residual error Emversus the auxiliary parametertaking A=1,α=0.5,nH=4,â=0.08,and Câ=0.02.(a)K(z)as linear,(b)K(z)as parabolic.

Fig.3 Validation of the HAM-based solution(for=−1.585)with numerical solution(4th order RK-method)and implicit solution considering linear pro file for K(z).A=1,α=0.50,nH=4,a/h=0.08,Ca/h=0.02,K(z)=z.

Fig.4 Validation of the HAM-based solution(for=−3.445)with numerical solution(4th order RK-method)and implicit solution considering parabolic pro file for K(z). A=1,α=0.50,nH =4,a/h=0.08,Ca/h=0.02,K(z)=z(1−z).

Now,we validate the obtained HAM-based solution by comparing it with implicit solution given in Eqs.(6)and(9),and with 4th order Runge-Kutta method.For that purpose,a set of feasible values of the required parameters is chosen randomly to assess the solutions.Figsures 3 and 4 present three types of solutions:HAM-Based,implicit,and numerical solutions,for linear and parabolic K(z)pro files,respectively.The required parameter values are mentioned in the text boxes inside the figures.The auxiliary parameteris obtained by minimizing Em;the values obtained are−1.585 for the 10th order approximation with linear K(z)and−3.445 for the 9th order approximation with parabolic K(z),and HAM-based series solutions are close to both numerical and implicit solutions as seen from the figures.

3.3 Variation of Concentration Pro file with nH and A

The effect of hindered settling on vertical sediment concentration pro file is shown in Fig.5.It is observed from the above figure that the concentration values increase due to the effect of hindered settling.Due to the hindered settling mechanism,the downward settling of a particle is prevented by the particles present in the surrounding.

Fig.5 In fluence of hindered settling through nHon the vertical concentration pro file.A=1,α=0.25,a/h=0.05,Ca/h=0.10.

Fig.6 Vertical concentration pro files for different values of the Rouse number A.a/h=0.05,Ca/h=0.10,nH=4,α=0.25.

Figure 5 presents two sediment concentration pro files with and without the effect of hindered settling.It is seen from the figure that the difference in concentration values is quite prominent.We observe that up to a very small distance from the channel bed and near the water surface,this difference is small.However,in the main flow region where suspension occurs,this difference is visible and decreases with the increasing vertical distance.The explanation behind these behaviors of concentration pro files is conveyed as follows.The hindered settling mechanism is not very active at the channel bed or near water surface,since the particles are not in suspension or are not mostly present there,respectively;however,in the main flow region,a large number of particles come into suspension mode due to excess bed shear stress,and the concentration difference increases thereat.On the other hand,the entrainment of sediment decreases with increasing vertical height which reduces the impact of nHand as a result,the difference of concentration values almost vanishes.

In Fig.6,the vertical normalized concentration proif les are plotted for different values of the Rouse number A(=0.25,0.5,1,and 1.5).A smaller value of A indicates finer sediment particles,and thus for a given vertical height the concentration value decreases with an increase in A as can be seen from the figure.For coarser sediment particles,the concentration value at the free surface is vanishingly small;however,a finite value at the free surface is meaningful for finer sediment particles as observed from the figures.Moreover,interestingly,near the channel bed(for z<0.2),the concentration curves are converging in nature for different Rouse numbers.

3.4 Further Modification to the Solution and Application to Laboratory Data

According to Liao’s rule of solution expression,one can choose the initial approximation,auxiliary linear operator,and auxiliary function with great freedom.In the previous section,we obtained the HAM-based series solution by choosing these functions.However,the solution obtained was seen to be close to the numerical solution for terms considering up to 9th order or so.Here,we further re fined our choices of these functions,in accordance with the rule of solution expression of HAM,to search for a reliable 3rd–4th order analytical solution to the problem.To that end,{((1−z)/z)n;n≥0}was chosen as the set of base functions to represent the solution C(z)for the problem.Accordingly,the auxiliary linear operator L and the initial approximation C0(z)were selected as follows:

The auxiliary function H(z)was chosen as 1 for simplicity,and the non-linear operator N was chosen the same as that of Eq.(20).Therefore,the final solution was obtained,considering parabolic pro file for K(z),as follows:

Fig.7 Validation of the HAM-Based Solution(Eq.(26)with Eq.(31)):(a)with numerical solution(4th order RK method)(A=2.212,α=−0.230,a/h=0.04,Ca/h=0.022,nH=2.624,h=−0.575)and(b)with Run S1 of Einstein and Chien[54]data.

The present study considered the most generalized governing equation by taking into account the physics of sediment-laden flow through the incorporation of hindered settling effect,arbitrary K(z)pro file,and the non-identity of the ratio of sediment diffusion coefficient to turbulent diffusion coefficient.Specifically,the equation modeled here is applicable for high-concentrated flows.To check the applicability of the model,a relevant set of laboratory data was chosen for verification.Experimental data of Einstein and Chien[54]was selected for verification as it is the most cited data for high-concentrated flows.They carried out their experiments in a recirculating flume of 12.2 m length and 30.7 cm width for three different particle sizes,namely,0.274 mm,940 mm,and 1.3 mm.The sediment concentration distribution near the channel bed is always a difficult task to measure,as it is under the in fluence of both sediment concentration and the channel bed,and hence,is of great importance.Einstein and Chien[54]measured concentration values very close to the channel bed,and reported 16 different measurement runs,named Run S1 to S16.We selected three test cases having three different particle diameters:S1(1.3 mm),S6(0.940 mm),and S12(0.274 mm)for validation.

Figures 7(a)and 7(b)show the validation of HAM-based solution(Eq.(26)together with Eq.(32))with the numerical solution(the 4th order RK method)and with the experimental data set S1 from Einstein and Chien,[54]respectively.The required parameter values were computed from the formulae mentioned in the previous section and are reported in Table 1.

Table 1 Summary of the selected experimental data from Einstein and Chien[17]for highconcentrated sediment-laden flow.

Fig.8 Validation of the HAM-Based Solution(Eq.(26)with Eq.(32)):(a)with numerical solution(the 4th order RK method)(A=0.202,α= −0.078,a/h= −0.043,Ca/h=0.0106,nH=2.562,h= −0.618)and(b)with Run S6 of Einstein and Chien[54]data.

Fig.9 Validation of the HAM-Based Solution(Eq.(26)with Eq.(32)):(a)with numerical solution(the 4th order RK method)(A=1.444,α=0.978,a/h=0.03,Ca/h=0.077,nH=3.562,h=−0.662)and(b)with Run S12 of Einstein and Chien[54]data.

It is observed that only the 3rd order approximate solution obtained by HAM given in Eq.(26)was close to the numerical one.Moreover,the corresponding solution can predict data well near the channel bed for highconcentrated sediment-laden flows as can be seen from the figure.In a similar manner,other two sets of data,namely,S6 and S12 are validated in Figs.8 and 9.In all cases,good agreement shows potential of the present model in the context of application to particle- fluid mixture flows where the volume of particles is appreciable.

4 Conclusions

The following conclusions can be drawn from this study:

(i)The study derives an explicit analytical solution to the most generalized governing equation for suspended sediment transport in open channels using a powerful mathematical tool,known as Homotopy Analysis Method(HAM).The method derives the solution for an arbitrary turbulent diffusion coefficient pro file.Unlike the classical perturbation techniques,HAM does not rely on the presence of small/large parameters in the governing equation.Besides,the method is a generalized one,which contains several other analytical methods such as ADM,HPM,classical perturbation method etc.as special cases,and it paves the way for monitoring the convergence of the series solution through an auxiliary parame-ter.In all,the method can be applied to the field of hydraulics where non-linear differential equations often arise.

(ii)The governing equation for vertical sediment concentration distribution is formulated considering the different phases of fluid and solid,and then coupling both of the phases:a concept which was originally developed by Hunt,and the corresponding equation is widely known as the Hunt equation. To deal with high-concentration flow which is of great importance,sediment and turbulent diffusion coefficients are assumed to be different from each other.The effect of hindered settling on sediment particles is included in the model,and the concentration proif le is analyzed due to this effect,which explains the real physics associated with high-concentrated flows.

(iii)To assess the derived analytical solution,required parameters are computed from the recent studies of the corresponding parameters.An implicit solution to the governing equation is also obtained for validating the HAM-based solution.Furthermore,the implicit solution,HAM-based solution,and the numerical solution(the 4th order RK method)are compared for different test cases considering two different pro files for eddy viscosity,namely,linear and parabolic pro file.

(iv)At the end,the HAM-based solution is further achieved for different choices of the linear operator and the initial approximation to the solution.The derived solution is then seen to be close to the numerical one only for three terms of the series.Finally,to check the applicability of the derived model,a widely cited laboratory dataset on vertical concentration pro file in high-concentrated flow is selected.A close agreement is found between the observed and computed values,which shows the potential of the model in the context of application to the highconcentrated free surface sediment-laden flows.

Acknowledgements

The authors are thankful to the Science and Engineering Research Board(SERB),Department of Science and Technology(DST),Govt.of India for providing the financial support through the research project with No.EMR/2015/002434.

The authors are thankful to Prof.Subhasish Dey(Dept.of Civil Engineering,IIT Kharagpur),the advisor of this project,for his valuable advices.The first author is thankful to Mr.Suvra Kanti Chakraborty,Senior Research Scholar,Dept.of Mathematics,IIT Kharagpur for discussion regarding computation.

Appendix A Convergence Analysis

The detailed proof of the convergence of the HAM-based series solution for the generalized governing equation given by Eq.(18)is presented here. In this regard,we recall Theorems 2.1 and 2.2 given in the book of Liao,[55]and prove them according to our governing equation Eq.(18).

Theorem 1 If the homotopy seriesandconverge,thendefined by the relation Eq.(24)satis fies

Proof The auxiliary linear operator is defined as follows

According to Eq.(15),we obtain

Adding all the above terms,we get

Proof We first recall the following definition:

De finition 1 Cauchy product of two in finite series:Letandbe two in finite series of real/complex terms.Then the Cauchy product of the above two series is defined by the discrete convolution as follows:

Therefore,using the above rule in relation to Eq.(24),we get

which is basically the original governing equation Eq.(18).Furthermore,subject to the initial conditionand the conditions for the higher-order deformation equation=0,for m ≥ 1,we easily obtain

Hence,the convergence result follows.