Granular flow through a vertical channelsteady state solution using kinetic model
Abstract
Granular materials are the most processed materials on earth finding their applications in various industries such as pharmaceutical, agricultural and food processing, constructionbased, 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 nonzero 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 frictionalkinetic 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.
Abbreviations
ODE  Ordinary Differential Equation 
DEM  Discrete Element Method 
FEM  Finite Element Method 
INTRODUCTION
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 NavierStokes’ model. This friction can be attributed to the solidlike 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 nonzero 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
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 understanding 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 semianalytical solutions, continuum models preferred over discrete models for large systems.
METHODOLOGY
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.
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 Nparameter approximate solution $U_{N}$( Reddy, J.N. (1988) )
The approximate solution U should satisfy the given ODE in weightedintegral sense
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$.
MODEL
Assumptions
1. Steady state flow
2. 1D 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 2D flow and simplified based on the assumptions
Continuity equation
Momentum balance
$x$ component:
$y$ component:
$\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,
Here, $d_p$ is the diameter of the particles, $\rho_p$ is the intrinsic density of particles and
$h_1h_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
Where
where $\tilde J$ is the rate of dissipation of pseudo thermal energy per unit volume due to inelasticity of interparticle collisions. $h_1h_7$are defined as follows ^{[5]}.
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
Where
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_1f_7$ becomes
On substituting (20) and (21) in the momentum balance equations, we got
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
$y$momentum
Energy balance
On substituting for $\tilde{\sigma_{xy}}$ and $\tilde{J}$ we get,
FINAL EQUATIONS AND BOUNDARY CONDITIONS
Final equations
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}(1e_w)^2$
SOLVING COUPLED EQUATIONS
The following equations were taken from https://www.youtube.com/watch?v=rJhlZCo0Gi4
Boundary conditions
u(0)=0, u(1)=1, v(0)=1, v(1)=e
RESULTS AND DISCUSSION
Code explanation
Firstly the residual of each ODE was multiplied with the appropriate weight functions and their weak forms were formulated as shown below:
$\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 nonlinear 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 cointern 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.
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;
tic
while(c>0)
%disp('iteration %d starts\n',count)
[U1,V1]=assembly(U, V, n, h);
c=0.0;
for i=1:nn
if(abs(U(i)U1(i))>AC)
c=c+1;
break;
end
if(abs(V(i)V1(i))>AC)
c=c+1;
break;
end
end
U=U1;
V=V1;
%disp('iteration %d ends\n',count)
count=count+1
end
disp('Hence solution:');
%output for primary and secondary variables
diff=abs(U(x.^2)');
diff1=abs(Vexp(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 variable
figure(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')
toc
end
%derivation of element matrix and assembly
function [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:n
lmm=[lmm;[i i+1]];
end
for i=1:n
lm=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=k111
k12=k313*V(lm(1))+k314*V(lm(2))+(k311*U(lm(1))+k312*U(lm(2)))
k21=k112
k22=k111k113
f1=f11;
f2=f12;
%assembly according to connectivity matrix
L11(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 conditions
L11(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\R
U1=d(1:nn)
V1=d(nn+1:2*nn)
end
function [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=[(1s)/2,(1+s)/2];
dns=[1/2,1/2];
coord=[(i1)*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);
end
end
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, 223256.

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

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.
Post your comments
Please try again.