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

Abstract

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

Abbreviations

Abbreviations
 ODE Ordinary Differential Equation DEM Discrete Element Method FEM Finite Element Method

Background

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.

PROBLEM DEFINITION

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

LITERATURE REVIEW

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.

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

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

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

$\int_{0}^{1}wRdx=0$

where $R$ is called the residual and $w$ is called a weight function.

Any approximate solution for an ODE does not satisfy the differential equation exactly and hence results in a residual $R$. The value obtained on the right hand side of the differential equation on substituting the assumed solution, in place of zero is the residual $R$. 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 $c_j$.

Assumptions

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

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

​ $x$ component:

$\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$

​ $y$ component:

$\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$

$\Sigma_{xx}$, $\Sigma_{xy}$and $\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,​

${\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}}$
${\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}}$
${\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, $d_p$ is the diameter of the particles, $\rho_p$ is the intrinsic density of particles and

$h_1-h_3$ are functions of $\nu$. T is the grain temperature defined as​ $\frac{3}{2}T=\frac{1}{2}<{v'}^2>$. Here $\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

$-\widetilde\nabla\cdot\widetilde q-\widetilde\sigma:\widetilde\nabla\widetilde u-\widetilde 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$

Where

$\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}}$
$\widetilde J={\displaystyle\frac{\rho_p}{d_p}}f_5\widetilde T^\frac32$

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

$h_1=\nu(1+4\nu g(\nu))$
$h_2=\frac{8}{3\sqrt{\pi}}\eta\nu^2g(\nu)$
$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))$
$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)$
$h_5=0$
$h_6=\frac{48}{\sqrt{\pi}}\eta(1-\eta)\nu^2g(\nu)$
$h_7=\frac{\pi\nu g(\nu)}{2\sqrt{3}\nu_{max}h_4(\nu)}$

Fully developed flow equations

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

$\Sigma_{xx} = \rho_pf_1(\nu)\tilde{T} = {\hat{\sigma}_{xx}}$
$\Sigma_{xy} = -\rho_pd_pf_2(\nu){\tilde{T}}^{\frac{1}{2}}\dfrac{d\tilde{v_y}}{d\tilde{x}} = \tilde{\sigma_{xy}}$

Where

$f_1=\nu(1+4\nu g(\nu))$
$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))$
$f_3=\frac8{3\sqrt\pi}\eta\nu^2g(\nu)$
$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)$
$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)$
$f_5=\frac{48}{\sqrt{\pi}}\eta(1-\eta)\nu^2g(\nu)$
$f_6=\frac{\sqrt{3}\pi\nu g(\nu)}{2\nu_{max}f_4(\nu)}$
$f_7=\frac{\pi\nu g(\nu)}{2\sqrt3\nu_{max}f_2(\nu)}$

where

$g(\nu)=\frac{1}{1-({\frac{\nu}{\nu_{max}})}^\frac{1}{3}}$ and $\eta=\frac{1}{2}(1+e_p)$

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

$f_1=\nu(1+4\nu g(\nu))$
$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)$
$f_3=\frac8{3\sqrt\pi}\nu^2g(\nu)$
$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)$
$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)$
$f_5=\frac{48}{\sqrt{\pi}}\eta(1-\eta)\nu^2g(\nu)$
$f_6=\frac{\sqrt{3}\pi\nu g(\nu)}{2\nu_{max}f_4(\nu)}$
$f_7=\frac{\pi\nu g(\nu)}{2\sqrt3\nu_{max}f_2(\nu)}$
$f_1=\nu(1+4\nu g(\nu))$
$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)$
$f_3=\frac8{3\sqrt\pi}\nu^2g(\nu)$
$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)$
$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)$
$f_5=\frac{48}{\sqrt{\pi}}\eta(1-\eta)\nu^2g(\nu)$
$f_6=\frac{\sqrt{3}\pi\nu g(\nu)}{2\nu_{max}f_4(\nu)}$
$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

$\frac{d\Sigma_{xx}}{dx}= \frac{d{\hat{\sigma}}_{xx}}{d\tilde{x}}=0$
$\frac{d\Sigma_{xy}}{dx}=-\rho{g}=\dfrac{d\tilde{\sigma_{xy}}}{d\tilde{x}}=\rho_p{\nu}g$

SCALING

$(x,y)=\frac{1}{\tilde{w}}(\tilde{x},\tilde{y})$, $(v_x,v_y)=\frac{1}{u_c}(\tilde{v_x},-\tilde{v_y})$, $T=\frac{\tilde{T}}{{(\frac{d_p}{\tilde{w}})^2}{\tilde{u_c}^2}}$,

​, ​​ $\sigma=\frac{\tilde{\sigma}}{{(\frac{d_p}{\tilde{w}})^2}{\tilde{u_c}^2\rho_p}}$ and ​ $q=\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

$x$-momentum

$\dfrac{d(f_1\tilde{T}\rho_p)}{d\tilde{x}}=0$
$\dfrac{d(f_1T)}{dx}=0$

$y$-momentum

$\frac d{dx}\left(f_2(\nu)\widetilde T^\frac12\frac{dv_y}{dx}\right)=-\nu w^2$

Energy balance

$\dfrac{d\tilde{q}}{d\tilde{x}}+\tilde{\sigma_{xy}}\dfrac{d\tilde{v_y}}{d\tilde{x}}+\tilde{J}=0$

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

${\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

$-\frac d{dx}\left(\mu\frac{du_y}{dx}\right)=\nu w^2$
$\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$
$\frac d{dx}\left(f_1T\right)=0$

Boundary conditions

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

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

Where

$\mu=f_2T^\frac{1}{2}$, $K=f_4T^\frac{1}{2}$, $K_1=f_2T^\frac12$, $K_2=f_5=0$, $K_3=-\phi'f_7$​,

$K_4=\frac{f_6}{3}\phi'$ and $K_5=\frac{f_6}{2}(1-e_w)^2$

SOLVING COUPLED EQUATIONS

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

$\frac{d^2u}{dt^2}+\left(\frac{\displaystyle dv}{\displaystyle dt}\right)^2+uv=2+e^t(e^t+t^2)$
$\frac{d^2u}{dt^2}+\frac{du}{dt}-v=2t$

Boundary conditions

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

RESULTS AND DISCUSSION

Solution
u Vs x
v Vs x
Results

Code explanation

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

$\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)$
$\left[w_2\frac{dv}{dt}\right]_0^1-\int_0^1\frac{dv}{dt}\frac{dw_2}{dt}+\int_0^1w_2\frac{du}{dt}dt-\int_0^1vw_2dt=\int_0^12tw_2dt$

$\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 ${L}{X}={R}$

CONCLUSION AND RECOMMENDATIONS

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.

Acknowledgement

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.

Appendices

function coupledFEM()n=10;nn=n+1;l=1;h=l/n;x=(0:h:l);AC=0.00005;U=zeros(nn,1);V=zeros(nn,1);U(1)=0;V(1)=1;U(nn)=1;V(nn)=exp(1);c=1.0;count=0;ticwhile(c>0)%disp('iteration %d starts\n',count)[U1,V1]=assembly(U, V, n, h);c=0.0;for i=1:nnif(abs(U(i)-U1(i))>AC)c=c+1;break;endif(abs(V(i)-V1(i))>AC)c=c+1;break;endendU=U1;V=V1;%disp('iteration %d ends\n',count)count=count+1enddisp('Hence solution:');%output for primary and secondary variablesdiff=abs(U-(x.^2)');diff1=abs(V-exp(x)');fprintf('No.of elements=%d\n',n)disp(' x u(FEM) u(Exact) Error v(FEM) v(Exact) Error')disp([x',U,(x.^2)',diff,V,exp(x)', diff1])fprintf('No of iterations=%d\n',count)%plotting of primary variablefigure(1)plot(x,U,'--rs','Linewidth',2)xlabel('x')ylabel('u(x)')title('Solution plot to given BVP')figure(2)plot(x,V,'--rs','Linewidth',2)xlabel('x')ylabel('v(x)')title('Solution plot to given BVP')tocend%derivation of element matrix and assemblyfunction [U1,V1]=assembly(U, V, n, h)nn=n+1;L11=zeros(nn,nn);L12=zeros(nn,nn);R1=zeros(nn,1);L21=zeros(nn,nn);L22=zeros(nn,nn);R2=zeros(nn,1);lmm=[];for i=1:nlmm=[lmm;[i i+1]];endfor i=1:nlm=lmm(i,:);%generation of element matrix(k11) and RHS matrix(f1)%integration using Gaussian Quadrature.[k111,k112,k113,k311,k312,k313,k314,f11,f12]=linearelement(h,i);k11=-k111k12=k313*V(lm(1))+k314*V(lm(2))+(k311*U(lm(1))+k312*U(lm(2)))k21=k112k22=-k111-k113f1=f11;f2=f12;%assembly according to connectivity matrixL11(lm,lm)=L11(lm,lm)+k11;L12(lm,lm)=L12(lm,lm)+k12;L21(lm,lm)=L21(lm,lm)+k21;L22(lm,lm)=L22(lm,lm)+k22;R1(lm)=R1(lm)+f1;R2(lm)=R2(lm)+f2;end%imposition of boundary conditionsL11(1,:)=0.0; L12(1,:)=0.0;L11(nn,:)=0.0; L12(nn,:)=0.0;L21(1,:)=0.0; L22(1,:)=0.0;L21(nn,:)=0.0; L22(nn,:)=0.0;L11(1,1)=1; L11(nn,nn)=1;L22(1,1)=1; L22(nn,nn)=1;%L21(1,1)=1; L21(nn,nn)=1;%L12(1,1)=1; L12(nn,nn)=1;R1(1,1)=U(1); R1(nn,1)=U(nn);R2(1,1)=V(1); R2(nn,1)=V(nn);%solution of system of equations(F1)K=[L11,L12;L21,L22]R=[R1;R2]d=K\RU1=d(1:nn)V1=d(nn+1:2*nn)endfunction [k111,k112,k113,k311,k312,k313,k314,f11,f12]=linearelement(h,i)gL=[-sqrt(3/5);0;sqrt(3/5)];gW=[5/9,8/9,5/9];k111=zeros(2);k112=zeros(2);k113=zeros(2); k221=zeros(2);k222=zeros(2);k311=zeros(2);k312=zeros(2);k313=zeros(2); k314=zeros(2);f11=zeros(2,1);f12=zeros(2,1);for k=1:length(gW)s=gL(k);w=gW(k);n=[(1-s)/2,(1+s)/2];dns=[-1/2,1/2];coord=[(i-1)*h;i*h];x=n*coord(:,1);J=h/2;dx=(1/J)*dns;k111=k111+J*w*(dx)'*dx;k112=k112+J*w*n'*dx;k113=k113+J*w*(n)'*n;k221=k221+J*w*n'*dx*n(1);k222=k222+J*w*n'*dx*n(2);k311=k311+J*w*(dx)'*dx*n(1);k312=k312+J*w*n'*dx*n(1);k313=k313+J*w*(n)'*n*n(2);k314=k314+J*w*n'*dx*n(2);f11=f11+J*w*n'*(2+exp(2*x)+exp(x)*x^2);f12=f12+J*w*n'*(2*x);endend
MATLAB code for a system of nonlinear coupled equations

References

• 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.

More
Written, reviewed, revised, proofed and published with