Summer Research Fellowship Programme of India's Science Academies

Granular flow through a vertical channel-steady state solution using kinetic model

Sai Sree Ramachandruni

Department of Chemical Engineering, National Institute of Technology, Warangal, Telangana, India, 506004

Prof K Kesava Rao

Professor, Department of Chemical Engineering, Indian Institute of Science, Bengaluru, Karnataka, 560012


Granular materials are the most processed materials on earth finding their applications in various industries such as pharmaceutical, agricultural and food processing, construction-based, and automotive industries. Hence, a better knowledge of granular flow will help in improving the efficiency of various industrial operations. Granular materials are essentially clusters of discrete macroscopic particles in contact. They exhibit the behaviour of both solids and fluids thereby increasing the complexity of the flow. They resist deformation on the application of shear stress like solids and flow like liquids on account of gravity. Unlike liquids, the shear stress in granular solids is independent of the velocity at low velocities, and results in a non-zero shear stress even when they are at rest. In this work, steady, fully developed flow through a vertical channel is examined. As the purely frictional model results in a physically unrealizable solution, the kinetic model or the frictional-kinetic models have to be used. Balance laws for mass, linear momentum, and energy were formulated for continuum models while ignoring the effect of the interstitial fluid and any localized torques. The intention was to predict the velocity profile for the kinetic model using the galerkin finite element method. However, due to a lack of time, only the solution of a system of two coupled nonlinear ordinary differential equations was examined using this method.

Keywords: kinetic model, fully developed flow, continuum model, balance laws


ODEOrdinary Differential Equation
DEMDiscrete Element Method
FEM Finite Element Method 



Granular materials are classified into three types based on the fluid present in the voids: Dry granular materials (fluid present in the voids is gas), saturated granular materials (fluid present in the voids is liquid) and partially saturated granular materials (some parts of the voids are filled by a gas and the rest by a liquid)[1]. Saturated and partially saturated granular materials are together called wet granular materials. Though granular materials exhibit some properties of fluids and solids, they differ from them in many aspects. Like hydraulic jump, there exists a phenomenon called granular jump. Shear stresses in granular materials arise due to friction which depends on thenormal force, unlike fluids where the stresses are due to viscous friction and depend on the shear rate. This implies that the stress remains unchanged with an increase in velocity. Also, dense granular solids dilate upon deformation. These considerations give rise to a model which is mathematically more complex than the Navier-Stokes’ model. This friction can be attributed to the solid-like behavior of the granular materials. Grains have no intrinsic thermal motion and hence the flow of grains cannot be approximated to that of the flow of dense gases [2]. They have a non-zero slip velocity and the mass flow rate from storage vessels has been found to be independent of the head, which makes them distinct from liquids.


The system under study consists of a channel having a rectangular cross section as shown in Fig.1. The channel is completely filled with coarse granular material and allowed to flow freely on account of gravity. The analysis is done in the region where the flow is fully developed, far from the channel ends. A steady flow may be obtained by continuously feeding material to top of the channel.

    Geometry of the channel


    The regimes of slow, dense flow and loose, rapid flow of grains are idealizations of grain interactions. In reality, the mode of contact in grains depends on elastic response of the material comprising the grains, the shape and surface roughness of the grains, the impact velocity, and various other factors. The total stress comprises of the contribution from kinetic and frictional stresses. The shear layer in the vicinity of the walls of the channel experiences maximum kinetic stress. The frictional stress becomes important in the central region of the flow. In the problem considered, the particles are coarse (of the order of 1 mm), and the drag arising due to the interstitial fluid is negligible compared to the other forces. But this cannot be the case with fine particles where the drag effects are significant. The interparticle collisions are assumed to be slightly inelastic. Taking this in to account, Lun et al. (1984) has developed a model by extending the kinetic theory of dense gases to granular flow. Granular flow can be modeled using discrete or continuum model. In discrete models, Newton’s laws were applied for each particle and solved numerically using Discrete Element Method (DEM). DEM becomes complex due to limited under-standing of contact forces between particles and between particles and wall. Though discrete models seem to be more realistic in case of granular materials, due to the difficulty in finding analytical or semi-analytical solutions, continuum models preferred over discrete models for large systems.


    Finite element method

    The solution for any second order ordinary differential equation (ODE) can be found using analytical or numerical methods. But, analytical methods are restricted to particular forms of ODEs only. Hence, we seek for an approximate solution with the numerical techniques which helps in finding the numerical values for unknowns in the domain. Finite element method, boundary element method and finite difference method are some of the widely used numerical techniques for solving ODEs. In this project, we have chosen finite element method to work with. In this method, the entire domain is divided into finite elements. An algebraic equation is written for each element as a linear combination of interpolation functions. The resulting equations from all the elements are then assembled after applying appropriate boundary conditions in such a way that the continuity is not broken at the nodes. For a domain divided into ‘N’ elements, the solution of the below form is assumed.

     uUNj=1Ncjϕj(x)+ϕo(x)  u\approx U_N\equiv\sum_{j=1}^{N}c_j\phi_j(x)+\phi_o(x) 

    Where ϕj(x)\phi_j(x) are the known interpolation function describing the jthj^{th} element in the domain and cjc_j are the coefficients to be found. ϕj(x)\phi_j(x) and ϕo(x)\phi_o(x) are selected in such a way that the given boundary conditions are satisfied by the N-parameter approximate solution UNU_{N}( Reddy, J.N. (1988) ​)

    The approximate solution U should satisfy the given ODE in weighted-integral sense


    where RR is called the residual and ww is called a weight function.

    Any approximate solution for an ODE does not satisfy the differential equation exactly and hence results in a residual RR. The value obtained on the right hand side of the differential equation on substituting the assumed solution, in place of zero is the residual RR. The residual of each element is multiplied with an appropriate weight function ‘w’ and integrated over the domain and equated to zero. This procedure is followed for every element of the domain which results in a set of algebraic equations. These equations are then then assembled by applying the continuity conditions at the boundary nodes of each element. This results in a matrix which is to be solved to find the coefficients cjc_j.



    1. Steady state flow

    2. 1-D flow

    3. Fully developed flow

    4. Drag due to interstitial fluid is negligible

    5. Channel thickness in Z direction>>W

    Governing equations

    The continuity and momentum equations were written for a 2-D flow and simplified based on the assumptions

    Continuity equation

    (ρvx)x+(ρvy)y=0{\displaystyle\frac{\partial(\rho v_x)}{\partial x}}+{\displaystyle\frac{\partial(\rho v_y)}{\partial y}}=0

    Momentum balance

    xx component:

    ρvxvxx+ρvyvxy+Σxxx+Σxyy=0\rho v_x{\displaystyle\frac{\partial v_x}{\partial x}}+\rho v_y{\displaystyle\frac{\partial v_x}{\partial y}}+{\displaystyle\frac{\partial\Sigma_{xx}}{\partial x}}+{\displaystyle\frac{\partial\Sigma_{xy}}{\partial y}}=0

    yy component:

    ρvxvyx+ρvyvyy+Σxxx+Σyyy+ρg=0\rho v_x{\displaystyle\frac{\partial v_y}{\partial x}}+\rho v_y{\displaystyle\frac{\partial v_y}{\partial y}}+{\displaystyle\frac{\partial\Sigma_{xx}}{\partial x}}+{\displaystyle\frac{\partial\Sigma_{yy}}{\partial y}}+\rho g=0

    Σxx\Sigma_{xx}, Σxy \Sigma_{xy} and Σyy\Sigma_{yy} are the components of the symmetric stress tensor(defined in the compressive sense), and g is the acceleration due to gravity. The equations for the kinetic stresses are taken from ​Lun et al. (1984)​ Considering only the kinetic contribution of the stress,​

    Σxx=Σxxk=ρpTh1(ν)ρpdpT(h2(ν)23h3(ν))(vxx+vyy)2ρpdpTh3(ν)vxx{\Sigma}_{xx}={\Sigma}^k_{xx}=\rho_pTh_1\left(\nu\right)-\rho_pd_p\sqrt{T}(h_2\left(\nu\right)-\frac{2}{3}h_3(\nu))\left(\dfrac{\partial{v_{x}}}{\partial{x}}+ \dfrac{\partial{v_{y}}}{\partial{y}}\right)-2\rho_pd_p\sqrt{T}h_3\left(\nu\right)\dfrac{\partial{v_{x}}}{\partial{x}}
    Σyy=Σyyk=ρpTh1(ν)ρpdpT(h2(ν)23h3(ν))(vxx+vyy)2ρpdpTh3(ν)vyy{\Sigma}_{yy}={\Sigma}^k_{yy}=\rho_pTh_1(\nu)-\rho_pd_p\sqrt{T}(h_2(\nu)-\frac{2}{3}h_3(\nu))\left(\dfrac{\partial{v_{x}}}{\partial{x}}+ \dfrac{\partial{v_{y}}}{\partial{y}}\right)-2\rho_pd_p\sqrt{T}h_3\left(\nu\right)\dfrac{\partial{v_{y}}}{\partial{y}}
    Σxy=Σxyk=ρpdpTh3(ν)(vyx+vxy){\Sigma}_{xy}={\Sigma}^k_{xy}=-\rho_pd_p\sqrt{T} h_3(\nu)\left(\dfrac{\partial{v_{y}}}{\partial{x}}+\dfrac{\partial{v_{x}}}{\partial{y}}\right)

    Here, dpd_p is the diameter of the particles, ρp\rho_p is the intrinsic density of particles and

    h1h3h_1-h_3 are functions of ν\nu. T is the grain temperature defined as​ 32T=12<v2>\frac{3}{2}T=\frac{1}{2}<{v'}^2>. Here 12<v2>\frac{1}{2}<{v'}^2> is the specific mean fluctuational kinetic energy of translation of a collection of particles.

    Energy balance

    The dependency of transport properties on the grain temperature results in the need for the momentum balance and the kinetic energy balance to be solved simultaneously

    ~q~σ~:~u~J~=0-\widetilde\nabla\cdot\widetilde q-\widetilde\sigma:\widetilde\nabla\widetilde u-\widetilde J=0
    dq~dx~+σxy~dvy~dx~+J~=0{\displaystyle\frac{d\widetilde q}{d\tilde x}}+\tilde{\sigma_{xy}}{\displaystyle\frac{d\tilde{v_y}}{d\tilde x}}+\widetilde J=0


    q~=ρpdpf4T~12dT~dx~ρpdpf4hT~32dνdx~\widetilde q=-\rho_pd_pf_4\widetilde T^\frac12{\displaystyle\frac{d\widetilde T}{d\widetilde x}}-\rho_pd_pf_{4h}\widetilde T^\frac32{\displaystyle\frac{d\nu}{d\widetilde x}}
    J~=ρpdpf5T~32\widetilde J={\displaystyle\frac{\rho_p}{d_p}}f_5\widetilde T^\frac32

    where J~\tilde J is the rate of dissipation of pseudo thermal energy per unit volume due to inelasticity of inter-particle collisions. h1h7h_1-h_7are defined as follows [5].​

    h1=ν(1+4νg(ν))h_1=\nu(1+4\nu g(\nu))
    h3=23π(5π96(1+85νg(ν))(1g(ν)+85ν)+85ν2g(ν))h_3=\frac{2}{3\sqrt{\pi}}(\frac{5\pi}{96}(1+\frac{8}{5}\nu g(\nu))(\frac{1}{g(\nu)}+\frac{8}{5}\nu)+\frac{8}{5}\nu^2g(\nu))
    h4=25π128(1+125ηνg(ν))(1g(ν)+125ην)+4πην2g(ν)h_4=\frac{25\sqrt{\pi}}{128}(1+\frac{12}{5}\eta\nu g(\nu))(\frac{1}{g(\nu)}+\frac{12}{5}\eta\nu)+ \frac{4}{\sqrt{\pi}}\eta\nu^2g(\nu)
    h7=πνg(ν)23νmaxh4(ν)h_7=\frac{\pi\nu g(\nu)}{2\sqrt{3}\nu_{max}h_4(\nu)}

    Fully developed flow equations

    Based on the assumptions made ( vy=vy(x),vx=0v_y=v_y(x), v_x=0and T=T(x), the above equation reduces to

    Σxx=ρpf1(ν)T~=σ^xx\Sigma_{xx} = \rho_pf_1(\nu)\tilde{T} = {\hat{\sigma}_{xx}}
    Σxy=ρpdpf2(ν)T~12dvy~dx~=σxy~\Sigma_{xy} = -\rho_pd_pf_2(\nu){\tilde{T}}^{\frac{1}{2}}\dfrac{d\tilde{v_y}}{d\tilde{x}} = \tilde{\sigma_{xy}}


    f1=ν(1+4νg(ν))f_1=\nu(1+4\nu g(\nu))
    f2=2+α3π(5π96η(2η)(1+85νηg(ν))(1g(ν)+85η(3η2)ν)+85ν2g(ν))f_2=\frac{2+\alpha'}{3\sqrt\pi}(\frac{5\pi}{96\eta(2-\eta)}(1+\frac85\nu\eta g(\nu))(\frac1{g(\nu)}+\frac85\eta(3\eta-2)\nu)+\frac85\nu^2g(\nu))
    f4=25π16η(4133η)(1+125ηνg(ν))(1g(ν)+125ην2(4η3)ν)+4πην2g(ν)f_4=\frac{25\sqrt\pi}{16\eta(41-33\eta)}(1+\frac{12}5\eta\nu g(\nu))(\frac1{g(\nu)}+\frac{12}5\eta\nu^2(4\eta-3)\nu)+\frac4{\sqrt\pi}\eta\nu^2g(\nu)
    f4h=25π16η(4133η)(1νg(ν)+12η5).12η5(2η1)(η1)(d(ν2g(ν))dν)f_{4h}=\frac{25\sqrt{\pi}}{16\eta(41-33\eta)}\left(\frac{1}{\nu g(\nu)}+\frac{12\eta}{5}\right).\frac{12\eta}{5}(2\eta-1)(\eta-1)\left(\frac{d(\nu^2g(\nu))}{d\nu}\right)
    f6=3πνg(ν)2νmaxf4(ν)f_6=\frac{\sqrt{3}\pi\nu g(\nu)}{2\nu_{max}f_4(\nu)}
    f7=πνg(ν)23νmaxf2(ν)f_7=\frac{\pi\nu g(\nu)}{2\sqrt3\nu_{max}f_2(\nu)}


    g(ν)=11(ννmax)13g(\nu)=\frac{1}{1-({\frac{\nu}{\nu_{max}})}^\frac{1}{3}} and η=12(1+ep)\eta=\frac{1}{2}(1+e_p)

    Substituting ep=1e_p=1, η=1\eta=1 except for the terms containing 1η1-\etaand taking the magnitude of (d(ν2g(ν))dν) \left(\frac{d(\nu^2g(\nu))}{d\nu}\right) to be zero since it is of negligible magnitude, f1f7f_1-f_7 becomes

    f1=ν(1+4νg(ν))f_1=\nu(1+4\nu g(\nu))
    f2=2+α3π(5π96(1+85νg(ν))(1g(ν)+85ν)+85ν2g(ν))f_2=\frac{2+\alpha'}{3\sqrt\pi}\left(\frac{5\pi}{96}\left(1+\frac85\nu g(\nu)\right)\left(\frac1{g(\nu)}+\frac85\nu\right)+\frac85\nu^2g(\nu)\right)
    f4=25π128(1+125νg(ν))(1g(ν)+125ν3)+4πν2g(ν)f_4=\frac{25\sqrt\pi}{128}\left(1+\frac{12}5\nu g(\nu)\right)\left(\frac1{g(\nu)}+\frac{12}5\nu^3\right)+\frac4{\sqrt\pi}\nu^2g(\nu)
    f4h=25π16η(4133η)(1νg(ν)+12η5).12η5(2η1)(η1)(d(ν2g(ν))dν)f_{4h}=\frac{25\sqrt{\pi}}{16\eta(41-33\eta)}\left(\frac{1}{\nu g(\nu)}+\frac{12\eta}{5}\right).\frac{12\eta}{5}(2\eta-1)(\eta-1)\left(\frac{d(\nu^2g(\nu))}{d\nu}\right)
    f6=3πνg(ν)2νmaxf4(ν)f_6=\frac{\sqrt{3}\pi\nu g(\nu)}{2\nu_{max}f_4(\nu)}
    f7=πνg(ν)23νmaxf2(ν)f_7=\frac{\pi\nu g(\nu)}{2\sqrt3\nu_{max}f_2(\nu)}
    f1=ν(1+4νg(ν))f_1=\nu(1+4\nu g(\nu))
    f2=2+α3π(5π96(1+85νg(ν))(1g(ν)+85ν)+85ν2g(ν))f_2=\frac{2+\alpha'}{3\sqrt\pi}\left(\frac{5\pi}{96}\left(1+\frac85\nu g(\nu)\right)\left(\frac1{g(\nu)}+\frac85\nu\right)+\frac85\nu^2g(\nu)\right)
    f4=25π128(1+125νg(ν))(1g(ν)+125ν3)+4πν2g(ν)f_4=\frac{25\sqrt\pi}{128}\left(1+\frac{12}5\nu g(\nu)\right)\left(\frac1{g(\nu)}+\frac{12}5\nu^3\right)+\frac4{\sqrt\pi}\nu^2g(\nu)
    f4h=25π16η(4133η)(1νg(ν)+12η5).12η5(2η1)(η1)(d(ν2g(ν))dν)f_{4h}=\frac{25\sqrt{\pi}}{16\eta(41-33\eta)}\left(\frac{1}{\nu g(\nu)}+\frac{12\eta}{5}\right).\frac{12\eta}{5}(2\eta-1)(\eta-1)\left(\frac{d(\nu^2g(\nu))}{d\nu}\right)
    f6=3πνg(ν)2νmaxf4(ν)f_6=\frac{\sqrt{3}\pi\nu g(\nu)}{2\nu_{max}f_4(\nu)}
    f7=πνg(ν)23νmaxf2(ν)f_7=\frac{\pi\nu g(\nu)}{2\sqrt3\nu_{max}f_2(\nu)}

    On substituting (20) and (21) in the momentum balance equations, we got

    dΣxxdx=dσ^xxdx~=0\frac{d\Sigma_{xx}}{dx}= \frac{d{\hat{\sigma}}_{xx}}{d\tilde{x}}=0


    (x,y)=1w~(x~,y~)(x,y)=\frac{1}{\tilde{w}}(\tilde{x},\tilde{y}), (vx,vy)=1uc(vx~,vy~)(v_x,v_y)=\frac{1}{u_c}(\tilde{v_x},-\tilde{v_y}), T=T~(dpw~)2uc~2T=\frac{\tilde{T}}{{(\frac{d_p}{\tilde{w}})^2}{\tilde{u_c}^2}},

    ​, ​​ σ=σ~(dpw~)2uc~2ρp\sigma=\frac{\tilde{\sigma}}{{(\frac{d_p}{\tilde{w}})^2}{\tilde{u_c}^2\rho_p}} and ​ q=q~(dpw~)4uc~3ρpq=\frac{\tilde{q}}{{(\frac{d_p}{\tilde{w}})^4}{\tilde{u_c}^3\rho_p}}

    For fully developed flow, the dimensionless equations are given by




    ddx(f2(ν)T~12dvydx)=νw2\frac d{dx}\left(f_2(\nu)\widetilde T^\frac12\frac{dv_y}{dx}\right)=-\nu w^2

    Energy balance


    On substituting for σxy~\tilde{\sigma_{xy}} and J~\tilde{J} we get,

    ddx(f4T12dTdx+f4hT32dνdx)w2[f5T32f2T12(duydx)2]=0{\displaystyle\frac d{dx}}\left(f_4T^\frac12{\displaystyle\frac{dT}{dx}}+f_{4h}T^\frac32{\displaystyle\frac{d\nu}{dx}}\right)-w^2\left[f_5T^\frac32-f_2T^\frac12\left({\displaystyle\frac{du_y}{dx}}\right)^2\right]=0


    Final equations

    ddx(μduydx)=νw2-\frac d{dx}\left(\mu\frac{du_y}{dx}\right)=\nu w^2
    ddx(KdTdx)+w2(f2T12(duydx)2f5T32)=0\frac d{dx}\left(K\frac{dT}{dx}\right)+w^2\left(f_2T^\frac12\left(\frac{du_y}{dx}\right)^2-f_5T^\frac32\right)=0
    ddx(f1T)=0\frac d{dx}\left(f_1T\right)=0

    Boundary conditions

    At x=0, duydx=0\frac{du_y}{dx}=0​; At x=0, dTdx=0\frac{dT}{dx}=0

    At x=1, duydx=K3uyww\frac{du_y}{dx}=K_3u_{yw}w; At x=1, dTdx=K4uywK5wT\frac{dT}{dx}=K_4{u_{yw}}K_5wT​​


    μ=f2T12\mu=f_2T^\frac{1}{2}, K=f4T12K=f_4T^\frac{1}{2}, K1=f2T12K_1=f_2T^\frac12, K2=f5=0K_2=f_5=0, K3=ϕf7K_3=-\phi'f_7​,

    K4=f63ϕK_4=\frac{f_6}{3}\phi' and K5=f62(1ew)2K_5=\frac{f_6}{2}(1-e_w)^2


    The following equations were taken from https://www.youtube.com/watch?v=rJhlZCo0Gi4

    d2udt2+(dvdt)2+uv=2+et(et+t2)\frac{d^2u}{dt^2}+\left(\frac{\displaystyle dv}{\displaystyle dt}\right)^2+uv=2+e^t(e^t+t^2)

    Boundary conditions

    u(0)=0, u(1)=1, v(0)=1, v(1)=e


      u Vs x
        v Vs x

          Code explanation

          Firstly the residual of each ODE was multiplied with the appropriate weight functions and their weak forms were formulated as shown below: 

          [w1dudt]0101dudtdw1dt+01(dvdt)2dt+01uvw1dt=012w1dt+01w1(e2t+ett2)\left[w_1\frac{du}{dt}\right]_0^1-\int_0^1\frac{du}{dt}\frac{dw_1}{dt}+\int_0^1\left(\frac{dv}{dt}\right)^2dt+\int_0^1uvw_1dt=\int_0^12w_1dt+\int_0^1w_1\left(e^2t+e^t\ast t^2\right)

          [L11(n+1)(n+1)L12(n+1)(n+1)L21(n+1)(n+1)L22(n+1)(n+1)][U(n+1)1V(n+1)1]=[R1(n+1)1R2(n+1)1]\begin{bmatrix} L11_{(n+1)*(n+1)} & L12_{(n+1)*(n+1)}\\ L21_{(n+1)*(n+1)} & L22_{(n+1)*(n+1)} \end{bmatrix} \quad \begin{bmatrix} U_{(n+1)*1} \\ V_{(n+1)*1} \end{bmatrix}=\begin{bmatrix} R1_{(n+1)*1}\\ R2_{(n+1)*1} \end{bmatrix} \quad

          Coupled ODEs have been solved using Galerkin finite element method. The domain was discretized in to 'n' equal elements. The solutions of u and v were assumed to be a linear combination of interpolation function of each element. The above matlab code was written to compute the values of u and v at every node. L is the global stiffness matrix, X is the unknown variable matrix and R is the right hand side matrix consisting of essential boundary conditions.  L11, L12, L21 and L22 are square matrices of order n+1 and together make the global stiffness matrix.  R1 and R2 are matrices of the order (n+1)*1 and contains the essential boundary conditions of equations (45) and (46) respectively. The natural boundary conditions have been employed in the stiffness matrix already. In the above code k11, k12, k21, k22 are the local stiffness matrices of L11, L12, L21, L22 respectively.  k11, k12,  k21, k22 are calculated for every element and slipped into the global matrix at appropriate positions. This is known as the assembly of the element equations. The stiffness matrix was further simplified by imposing the boundary conditions. Using left division command in matlab the inverse of the stiffness matrix L has been multiplied with the right hand side matrix to find the values of unknowns at each node. The final matrix will take the form LX=R{L}{X}={R}


          In this work, an initiative has been taken to solve coupled ODEs using Galerkin finite element method. The equations considered were only slightly non-linear which enabled us to solve them using our basic understanding on linearization. The same algorithm can be extended for solving the equations governing the granular flow. The above code uses linear interpolation function. Higher order interpolation function can be employed for better accuracy. Owing to the strong non linearity of the equations, it requires advanced linearization techniques in solving them.


          I would like to acknowledge all those people who have made this research work possible, an experience that I will cherish forever. 

          First of all, I am grateful to the Indian Academy of Sciences, Bangalore for giving me this opportunity to work at Indian Institute of Science, Bangalore, one of the top institutes for scientific research and higher education in India. 

          I owe my sincere thanks and the deepest sense of gratitude to my guide, Prof. K Kesava Rao, who has been a source of inspiration and support throughout the period of the project, for his guidance and constant encouragement. I have been extremely fortunate to have worked under his supervision.

          I would like to extend my heartiest thanks to Prof. K Kesava Rao’s research scholar Mr. Bhanjan Debnath and my co-intern Ms. Vaishnavee Manivannan, BITS Goa for their inspirational guidance, timely help, valuable suggestions and discussions throughout my project. 

          Further, I thank Prof. Satya Sai, Visiting faculty, Department of Chemical engineering, National Institute of Technology, Warangal for recommending me for this internship.

          Finally, I thank all the supporting staff of AuthorCafe for providing such an excellent space to write my report and collaborate with my guide. 


          1. function coupledFEM()
          2. n=10;
          3. nn=n+1;
          4. l=1;
          5. h=l/n;
          6. x=(0:h:l);
          7. AC=0.00005;
          8. U=zeros(nn,1);
          9. V=zeros(nn,1);
          10. U(1)=0;
          11. V(1)=1;
          12. U(nn)=1;
          13. V(nn)=exp(1);
          14. c=1.0;
          15. count=0;
          16. tic
          17. while(c>0)
          18. %disp('iteration %d starts\n',count)
          19. [U1,V1]=assembly(U, V, n, h);
          20. c=0.0;
          21. for i=1:nn
          22. if(abs(U(i)-U1(i))>AC)
          23. c=c+1;
          24. break;
          25. end
          26. if(abs(V(i)-V1(i))>AC)
          27. c=c+1;
          28. break;
          29. end
          30. end
          31. U=U1;
          32. V=V1;
          33. %disp('iteration %d ends\n',count)
          34. count=count+1
          35. end
          36. disp('Hence solution:');
          37. %output for primary and secondary variables
          38. diff=abs(U-(x.^2)');
          39. diff1=abs(V-exp(x)');
          40. fprintf('No.of elements=%d\n',n)
          41. disp(' x u(FEM) u(Exact) Error v(FEM) v(Exact) Error')
          42. disp([x',U,(x.^2)',diff,V,exp(x)', diff1])
          43. fprintf('No of iterations=%d\n',count)
          44. %plotting of primary variable
          45. figure(1)
          46. plot(x,U,'--rs','Linewidth',2)
          47. xlabel('x')
          48. ylabel('u(x)')
          49. title('Solution plot to given BVP')
          50. figure(2)
          51. plot(x,V,'--rs','Linewidth',2)
          52. xlabel('x')
          53. ylabel('v(x)')
          54. title('Solution plot to given BVP')
          55. toc
          56. end
          57. %derivation of element matrix and assembly
          58. function [U1,V1]=assembly(U, V, n, h)
          59. nn=n+1;
          60. L11=zeros(nn,nn);L12=zeros(nn,nn);R1=zeros(nn,1);
          61. L21=zeros(nn,nn);L22=zeros(nn,nn);R2=zeros(nn,1);
          62. lmm=[];
          63. for i=1:n
          64. lmm=[lmm;[i i+1]];
          65. end
          66. for i=1:n
          67. lm=lmm(i,:);
          68. %generation of element matrix(k11) and RHS matrix(f1)
          69. %integration using Gaussian Quadrature.
          70. [k111,k112,k113,k311,k312,k313,k314,f11,f12]=linearelement(h,i);
          71. k11=-k111
          72. k12=k313*V(lm(1))+k314*V(lm(2))+(k311*U(lm(1))+k312*U(lm(2)))
          73. k21=k112
          74. k22=-k111-k113
          75. f1=f11;
          76. f2=f12;
          77. %assembly according to connectivity matrix
          78. L11(lm,lm)=L11(lm,lm)+k11;
          79. L12(lm,lm)=L12(lm,lm)+k12;
          80. L21(lm,lm)=L21(lm,lm)+k21;
          81. L22(lm,lm)=L22(lm,lm)+k22;
          82. R1(lm)=R1(lm)+f1;
          83. R2(lm)=R2(lm)+f2;
          84. end
          85. %imposition of boundary conditions
          86. L11(1,:)=0.0; L12(1,:)=0.0;
          87. L11(nn,:)=0.0; L12(nn,:)=0.0;
          88. L21(1,:)=0.0; L22(1,:)=0.0;
          89. L21(nn,:)=0.0; L22(nn,:)=0.0;
          90. L11(1,1)=1; L11(nn,nn)=1;
          91. L22(1,1)=1; L22(nn,nn)=1;
          92. %L21(1,1)=1; L21(nn,nn)=1;
          93. %L12(1,1)=1; L12(nn,nn)=1;
          94. R1(1,1)=U(1); R1(nn,1)=U(nn);
          95. R2(1,1)=V(1); R2(nn,1)=V(nn);
          96. %solution of system of equations(F1)
          97. K=[L11,L12;L21,L22]
          98. R=[R1;R2]
          99. d=K\R
          100. U1=d(1:nn)
          101. V1=d(nn+1:2*nn)
          102. end
          103. function [k111,k112,k113,k311,k312,k313,k314,f11,f12]=linearelement(h,i)
          104. gL=[-sqrt(3/5);0;sqrt(3/5)];
          105. gW=[5/9,8/9,5/9];
          106. k111=zeros(2);k112=zeros(2);k113=zeros(2); k221=zeros(2);k222=zeros(2);
          107. k311=zeros(2);k312=zeros(2);k313=zeros(2); k314=zeros(2);
          108. f11=zeros(2,1);f12=zeros(2,1);
          109. for k=1:length(gW)
          110. s=gL(k);w=gW(k);
          111. n=[(1-s)/2,(1+s)/2];
          112. dns=[-1/2,1/2];
          113. coord=[(i-1)*h;i*h];
          114. x=n*coord(:,1);
          115. J=h/2;
          116. dx=(1/J)*dns;
          117. k111=k111+J*w*(dx)'*dx;
          118. k112=k112+J*w*n'*dx;
          119. k113=k113+J*w*(n)'*n;
          120. k221=k221+J*w*n'*dx*n(1);
          121. k222=k222+J*w*n'*dx*n(2);
          122. k311=k311+J*w*(dx)'*dx*n(1);
          123. k312=k312+J*w*n'*dx*n(1);
          124. k313=k313+J*w*(n)'*n*n(2);
          125. k314=k314+J*w*n'*dx*n(2);
          126. f11=f11+J*w*n'*(2+exp(2*x)+exp(x)*x^2);
          127. f12=f12+J*w*n'*(2*x);
          128. end
          129. end
          MATLAB code for a system of nonlinear coupled equations


          • Rao, K. K., Nott, P. R. (2008) An Introduction to Granular Flow, Cambridge U. Press.

          • Jaeger, H.M., Nagel, S. R., and Behringer, R. P. (1996). Granular solids, liquids, and gases. Rev. Mod. Phys. 68, 1259–1273.

          • Lun, C. K. K., Savage, S. B., Jeffrey, D. J. and Chepurniy, N. (1984) Kinetic theories for granular flow: inelastic particles in Couette flow and slightly inelastic particles in a general flow field. J. Fluid Mech. 140, 223-256.

          • Reddy, J. N. (1988) An Introduction to the Finite Element Method, Lecture Notes in Engineering, Finite Element Analysis for Engineering Design, 41-70

          • Mohan, L. S., Nott, P. R., and Rao, K. K. (1997). Fully developed flow of coarse granular materials through a vertical channel. Chem. Eng. Sci. 52, 913–933.

          Written, reviewed, revised, proofed and published with