Flow of Granular Particles on an Inclined plane in LIGGGHTS EN 649 - - PowerPoint PPT Presentation

DEM simulation of Flow of Granular Particles

• n an Inclined plane

in LIGGGHTS

EN 649

Mohit Prateek Roll No. 09D02017

Results and Discussions Post processing Simulation Introduction to geometry and pre simulation values Introduction to simulation software

• LAMMPS is a classical molecular dynamics code,

and an acronym for Large-scale Atomic/Molecular Massively Parallel Simulator.

• LAMMPS has potentials for soft materials

(biomolecules, polymers) and solid-state materials (metals, semiconductors) and coarse- grained or mesoscopic systems. It can be used to model atoms or, more generically, as a parallel particle simulator at the atomic, meso, or continuum scale.

• LAMMPS also offers a "GRANULAR" package for

DEM simulations.

LAMMPS

• LIGGGHTS stands for LAMMPS Improved for General Granular

and Granular Heat Transfer Simulations.

• As this name implies, it is based on the Open Source MD code

LAMMPS.

• LIGGGHTS now brings these DEM features to a new level. The

following features have been implemented on top of the LAMMPS "GRANULAR" features:

 A re-write of the contact formulations, including the

possibility to define macroscopic particle cohesion

 Import and handling of triangular meshes from CAD  A moving mesh feature  Improved particle insertion  A model for heat generation and conduction between

particles in contact

LIGGGHTS

Geometrical Representation

*.STL File; Viewed by Paraview

• Region of Simulation: 10 cm x 4 cm x 4 cm
• SI units
• Inclined Plane:
• Base : 5 cm
• Angle: 13.92 degrees
• Gravity: 9.81 m/s/s

Pre-Simulation Values

• Insertion of 5000 atoms
• Diameter: 0.5 mm
• Volume Fraction: 0.7
• Density: 2500
• Initial velocity: 0,0,0
• Coefficient of Restitution: 0.9
• Young's Modulus: 5E6
• Poisson's ratio: 0.45

Granular Particle Properties:

• ITEM: ATOMS id type type x y z vx vy vz fx fy fz tqx tqy tqz omegax omegay omegaz radius
• First Timestep
• 254 1 1 -0.0455105 -0.0194732 0.0172402 0 0 -0.271634 0 0 -1.60516e-06 0 0 0 0 0 0 0.00025
• 346 1 1 -0.0449629 -0.0188075 0.016992 0 0 -0.280454 0 0 -1.60516e-06 0 0 0 0 0 0 0.00025
• … … …. … … … … …. … … … … …. … … … …. … … … … …. … … … … …. … … … 0.00025
• …. … … … … …. … … … … …. … … … … …. … … … … …. … … … … …. … … … … …. …
• Next Timestep
• 2141 1 1 -0.0435509 -0.018768 0.01721 0 0 -0.272721 0 0 -1.60516e-06 0 0 0 0 0 0 0.00025
• 4889 1 1 -0.0431077 -0.0196025 0.0169182 0 0 -0.283023 0 0 -1.60516e-06 0 0 0 0 0 0 0.00025
• … … …. … … … … …. … … … … …. … … … …. … … … … …. … … … … …. … … … 0.00025
• …. … … … … …. … … … … …. … … … … …. … … … … …. … … … … …. … … … … …. …

Simulation

*.DUMP File; Viewed by VMD

Simulation Video

*.DUMP File; Viewed by VMD

Post Processing

SCILAB

LEVEL 2 LEVEL 1

FUNCTIONS MAKE3D.SCI REMOVEX.SCI SORTBYCOLUMN.SCI

MAKE3D.SCI

function [M2]=make3d(M, division_size) [rows, column] = size(M); ds = division_size; M2 = M(1:ds,:,:); for i = 2:(rows/ds) M2(:,:,i) = M(((i*ds)-(ds- 1)):(i*ds),:,:); end return endfunction function [Mx2]=removex(Mx) len = length(Mx(:,4)); for i = 10:10:len if Mx(i,4) > 0.05 then Mx = Mx(i:len,:); i = i - 10; else break; end len = length(Mx(:,4)); end Mx2 = Mx; return endfunction

REMOVEX.SCI

Function Description

SORT_COLUMN_ROWWISE2D.SCI

function [A, k]=sort_column_rowwise2d(a, column_number) cs = column_number; [B,k]=gsort(a(:,cs),'g'); [r,c] = size(a); A = rand(r,c); for i = 1:r A(i,:) = a (k(i),:); end return endfunction function mohitplot() a=gca(); a.font_size=2; poly1= a.children.children(1); poly1.thickness = 3; a.title.font_size = 5; a.x_label.font_size = 3.5; a.y_label.font_size = 3.5; xgrid endfunction

MOHITPLOT.SCI

Function Description

Read File using fscanfMat()

• Convert it into a hypermatrix
• Extract different values from the matrix like id, velocity, omega

Calculate Quantities required

• V = sqrt(vx.*vx + vy.*vy + vz.*vz);
• F_atom = sqrt(fx.*fx + fy.*fy + fz.*fz);
• T_atom = sqrt(tx.*tx + ty.*ty + tz.*tz);
• KE_atom = (1/2)*2500*(4/3*%pi*(0.00025^3))*(vx.*vx + vy.*vy + vz.*vz);
• RE_atom = (1/2)*(2/5)*2500*(4/3*%pi*(0.00025^3))*(0.00025^2)*(ox.*ox + oy.*oy + oz.*oz);
• KE_RE_atom = KE_atom + RE_atom

READING FILE

stacksize('max'); // To increase the limit in Scilab ! M_raw = fscanfMat('Default_edited.flow' ); M_raw = M_raw(:,1:18); // Removing radius column ! M = make3d(M_raw,5000); // There is a loss of data after if the division size is not a multiple

• f division size !

[row, column, rc] = size(M);

// Reading file and naming it ! id = M(:,1,:); x = M(:,4,:); y = M(:,5,:); z = M(:,6,:); vx = M(:,7,:); vy = M(:,8,:); vz = M(:,9,:); fx = M(:,10,:); fy = M(:,11,:); fz = M(:,12,:); tx = M(:,13,:); ty = M(:,14,:); tz = M(:,15,:);

• x = M(:,16,:);
• y = M(:,17,:);
• z = M(:,18,:);

// File reading done !

EXTRACTING VALUES

Code Description

FOR SINGLE PARTICLE

// Calculations start ! v = sqrt(vx.*vx + vy.*vy + vz.*vz); F_atom = sqrt(fx.*fx + fy.*fy + fz.*fz); T_atom = sqrt(tx.*tx + ty.*ty + tz.*tz); KE_atom = (1/2)*2500*(4/3*%pi*(0.00025^3))*(v x.*vx + vy.*vy + vz.*vz); RE_atom = (1/2)*(2/5)*2500*(4/3*%pi*(0.00025 ^3))*(0.00025^2)*(ox.*ox + oy.*oy +

• z.*oz);

KE_RE_atom = KE_atom + RE_atom;

for i =1:rc vtotal(i) = sum(v(:,:,i)); KE(i) = sum(KE_atom(:,:,i)); RE(i) = sum(RE_atom(:,:,i)); KE_RE(i) = sum(KE_RE_atom(:,:,i)); F(i) = sum(F_atom(:,:,i)); T(i) = sum(T_atom(:,:,i)); end

FOR ALL PARTICLES

Code Description

Graphs

SCILAB

LEVEL 4 LEVEL 3 LEVEL 2 LEVEL 1

GRAPHS WRT TIME

FOR ANGLE 13

FOR ANGLE 20 WRT AXIS FOR Z’-AXIS FOR ANGLE 13 FOR ANGLE 20 FOR X’-AXIS FOR ANGLE 13 FOR ANGLE 20

SAMPLE CODE FOR GENERATING A GRAPH

t = 1:1:rc; l = 1100; b = 750; scf(1); f=gcf(); // Create a figure f.figure_size= [l,b]; plot(t,vtotal); mohitplot(); xtitle("Variation of Velocity with Time", "Time (s)", "Velocity (m/s)");

Code Description

Variation of Velocity with time

Variation of Translational Kinetic Energy with time

Variation of Rotational Kinetic Energy with time

Variation of Total Kinetic Energy with time

Variation of Force with time

Variation of Energy with time

Rotate axis using the rotation matrix

• Sort it in increasing or decreasing order of X’ or Z’ axis
• Add these coordinates to the hypermatrix
• Group it into bins of 1000 and calculate a mean for each of them

Calculate Quantities required

• v_mean_x = sqrt(vx_x.*vx_x + vy_x.*vy_x + vz_x.*vz_x);
• F_mean_x = sqrt(fx_x.*fx_x + fy_x.*fy_x + fz_x.*fz_x);
• T_mean_x = sqrt(tx_x.*tx_x + ty_x.*ty_x + tz_x.*tz_x);
• KE_mean_x = (1/2)*2500*(4/3*%pi*(0.00025^3))*(vx_x.*vx_x + vy_x.*vy_x + vz_x.*vz_x);
• RE_mean_x = (1/2)*(2/5)*2500*(4/3*%pi*(0.00025^3))*(0.00025^2)*(ox_x.*ox_x +
• y_x.*oy_x + oz_x.*oz_x);
• KE_RE_mean_x = KE_mean_x + RE_mean_x;

ROTATING AND ADDING TO THE HYPERMATRIX

// Rotating the axis ! theta = 13.93*%pi/180; x1 = x*cos(theta) - z*sin(theta); y1 = y; z1 = x*sin(theta) + z*cos(theta); // Done rotating M(:,19,:) = x1; M(:,20,:) = y1; M(:,21,:) = z1; M_raw(:,19) = x1; M_raw(:,20) = y1; M_raw(:,21) = z1;

[Mx, sort_index] = sort_column_rowwise2d(M_raw, 19); // Sorting in decreasing oreder of x' coordinate for i = 1:279 xmean(i) = mean(Mx(((i-1)*500+1):(i*500),19)); vx_x(i) = mean(Mx(((i-1)*500+1):(i*500),7)); vy_x(i) = mean(Mx(((i-1)*500+1):(i*500),8)); vz_x(i) = mean(Mx(((i-1)*500+1):(i*500),9)); fx_x(i) = mean(Mx(((i-1)*500+1):(i*500),10)); fy_x(i) = mean(Mx(((i-1)*500+1):(i*500),11)); fz_x(i) = mean(Mx(((i-1)*500+1):(i*500),12)); tx_x(i) = mean(Mx(((i-1)*500+1):(i*500),13)); ty_x(i) = mean(Mx(((i-1)*500+1):(i*500),14)); tz_x(i) = mean(Mx(((i-1)*500+1):(i*500),15));

• x_x(i) = mean(Mx(((i-1)*500+1):(i*500),16));
• y_x(i) = mean(Mx(((i-1)*500+1):(i*500),17));
• z_x(i) = mean(Mx(((i-1)*500+1):(i*500),18));

End v_mean_x = sqrt(vx_x.*vx_x + vy_x.*vy_x + vz_x.*vz_x); F_mean_x = sqrt(fx_x.*fx_x + fy_x.*fy_x + fz_x.*fz_x); T_mean_x = sqrt(tx_x.*tx_x + ty_x.*ty_x + tz_x.*tz_x); KE_mean_x = (1/2)*2500*(4/3*%pi*(0.00025^3))*(vx_x.*vx_x + vy_x.*vy_x + vz_x.*vz_x); RE_mean_x = (1/2)*(2/5)*2500*(4/3*%pi*(0.00025^3))*(0.00025^2)*(ox _x.*ox_x + oy_x.*oy_x + oz_x.*oz_x); KE_RE_mean_x = KE_mean_x + RE_mean_x;

GROUPING AND CALCULATING

Code Description

LEVEL 4 LEVEL 3 LEVEL 2 LEVEL 1

GRAPHS WRT TIME

FOR ANGLE 13

FOR ANGLE 20 WRT AXIS FOR Z’-AXIS FOR ANGLE 13 FOR ANGLE 20 FOR X’-AXIS FOR ANGLE 13 FOR ANGLE 20

Variation of Velocity along the incline

Variation of Translational Kinetic Energy along the incline

Variation of Rotational Kinetic Energy along the incline

Variation of Total Kinetic Energy along the incline

Variation of Force along the incline

Variation of Energy along the incline

Variation of Velocity normal to the incline

Variation of Translational Kinetic Energy normal to the incline

Variation of Rotational Kinetic Energy normal to the incline

Variation of Total Kinetic Energy normal to the incline

Variation of Force normal to the incline

Variation of Energy normal to the incline

• Rheology and Segregation of Granular Mixtures –

Anurag Tripathi, 2010

• Velocity correlations in dense gravity-driven granular

chute flow - Oleh Baran, Deniz Ertas, and Thomas C. Halsey, 2006

• Flow rule of dense granular flows down a rough incline
• Tamás Börzsönyi and Robert E. Ecke, 2007
• Transition in a dense granular flow - V. Kumaran and S.

Maheshwari

• Discrete Element Modelling Of Granular Snow Particles

using LIGGGHTS - Vinodh Vedachalam, 2011

References

Special Thanks To

Prateek Maheshwari Aditya Telang Chaitanya Wadi

Thank You ! :)

Mohit Prateek Roll No. 09D02017