DEM simulation of Flow of Granular Particles on an Inclined plane in LIGGGHTS EN 649 Mohit Prateek Roll No. 09D02017
Introduction to simulation software Introduction to geometry and pre simulation values Simulation Post processing Results and Discussions
LAMMPS 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.
LIGGGHTS 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
Geometrical Representation *.STL File; Viewed by Paraview
Pre-Simulation Values 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
Granular Particle Properties: 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
Simulation *.DUMP File; Viewed by VMD • 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 Video *.DUMP File; Viewed by VMD
Post Processing SCILAB
LEVEL 1 LEVEL 2 MAKE3D.SCI FUNCTIONS REMOVEX.SCI SORTBYCOLUMN.SCI
Function Description MAKE3D.SCI REMOVEX.SCI function [ Mx2 ]=removex( Mx ) function [ M2 ]=make3d( M , len = length( Mx (:,4)); division_size ) [rows, column] = for i = 10:10:len size( M ); if Mx (i,4) > 0.05 then ds = division_size ; Mx = Mx (i:len,:); M2 = M (1:ds,:,:); i = i - 10; for i = 2:(rows/ds) else break; M2 (:,:,i) = M (((i*ds)-(ds- end 1)):(i*ds),:,:); len = length( Mx (:,4)); end end Mx2 = Mx ; return return endfunction endfunction
Function Description SORT_COLUMN_ROWWISE2D.SCI MOHITPLOT.SCI function [ A , function mohitplot() k ]=sort_column_rowwise2d( a , a=gca(); column_number ) a.font_size=2; cs = column_number ; poly1= a.children.children(1); [B, k ]=gsort( a (:,cs),'g'); poly1.thickness = 3; [r,c] = size( a ); a.title.font_size = 5; A = rand(r,c); a.x_label.font_size = 3.5; for i = 1:r a.y_label.font_size = 3.5; A (i,:) = a ( k (i),:); xgrid end endfunction return endfunction
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
Code Description READING FILE EXTRACTING VALUES stacksize('max'); // To increase the limit in Scilab ! // Reading file and naming it ! id = M(:,1,:); x = M(:,4,:); M_raw = fscanfMat('Default_edited.flow' y = M(:,5,:); ); z = M(:,6,:); M_raw = M_raw(:,1:18); // vx = M(:,7,:); Removing radius column ! vy = M(:,8,:); vz = M(:,9,:); fx = M(:,10,:); M = make3d(M_raw,5000); // There is a loss of data after if fy = M(:,11,:); the division size is not a multiple fz = M(:,12,:); of division size ! tx = M(:,13,:); [row, column, rc] = size(M); ty = M(:,14,:); tz = M(:,15,:); ox = M(:,16,:); oy = M(:,17,:); oz = M(:,18,:); // File reading done !
Code Description FOR SINGLE PARTICLE FOR ALL PARTICLES // Calculations start ! for i =1:rc v = sqrt(vx.*vx + vy.*vy + vz.*vz); vtotal(i) = sum(v(:,:,i)); F_atom = sqrt(fx.*fx + fy.*fy + fz.*fz); KE(i) = sum(KE_atom(:,:,i)); T_atom = sqrt(tx.*tx + ty.*ty + tz.*tz); KE_atom = RE(i) = sum(RE_atom(:,:,i)); (1/2)*2500*(4/3*%pi*(0.00025^3))*(v KE_RE(i) = x.*vx + vy.*vy + vz.*vz); sum(KE_RE_atom(:,:,i)); RE_atom = (1/2)*(2/5)*2500*(4/3*%pi*(0.00025 F(i) = sum(F_atom(:,:,i)); ^3))*(0.00025^2)*(ox.*ox + oy.*oy + oz.*oz); T(i) = sum(T_atom(:,:,i)); KE_RE_atom = KE_atom + RE_atom; end
Graphs SCILAB
LEVEL 1 LEVEL 2 LEVEL 3 LEVEL 4 FOR ANGLE 13 WRT TIME FOR ANGLE 20 GRAPHS FOR ANGLE 13 FOR Z’ -AXIS FOR ANGLE 20 WRT AXIS FOR ANGLE 13 FOR X’ -AXIS FOR ANGLE 20
Code Description 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)");
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 + oy_x.*oy_x + oz_x.*oz_x); • KE_RE_mean_x = KE_mean_x + RE_mean_x;
Code Description GROUPING AND CALCULATING ROTATING AND ADDING TO THE HYPERMATRIX [Mx, sort_index] = sort_column_rowwise2d(M_raw, 19); // Sorting in decreasing oreder of x' coordinate // Rotating the axis ! for i = 1:279 xmean(i) = mean(Mx(((i-1)*500+1):(i*500),19)); theta = 13.93*%pi/180; x1 = x*cos(theta) - z*sin(theta); vx_x(i) = mean(Mx(((i-1)*500+1):(i*500),7)); vy_x(i) = mean(Mx(((i-1)*500+1):(i*500),8)); y1 = y; vz_x(i) = mean(Mx(((i-1)*500+1):(i*500),9)); fx_x(i) = mean(Mx(((i-1)*500+1):(i*500),10)); z1 = x*sin(theta) + z*cos(theta); fy_x(i) = mean(Mx(((i-1)*500+1):(i*500),11)); fz_x(i) = mean(Mx(((i-1)*500+1):(i*500),12)); // Done rotating 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)); ox_x(i) = mean(Mx(((i-1)*500+1):(i*500),16)); oy_x(i) = mean(Mx(((i-1)*500+1):(i*500),17)); oz_x(i) = mean(Mx(((i-1)*500+1):(i*500),18)); End M(:,19,:) = x1; M(:,20,:) = y1; 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); M(:,21,:) = z1; T_mean_x = sqrt(tx_x.*tx_x + ty_x.*ty_x + tz_x.*tz_x); KE_mean_x = M_raw(:,19) = x1; (1/2)*2500*(4/3*%pi*(0.00025^3))*(vx_x.*vx_x + vy_x.*vy_x + vz_x.*vz_x); M_raw(:,20) = y1; RE_mean_x = (1/2)*(2/5)*2500*(4/3*%pi*(0.00025^3))*(0.00025^2)*(ox M_raw(:,21) = z1; _x.*ox_x + oy_x.*oy_x + oz_x.*oz_x); KE_RE_mean_x = KE_mean_x + RE_mean_x;
LEVEL 1 LEVEL 2 LEVEL 3 LEVEL 4 FOR ANGLE 13 WRT TIME FOR ANGLE 20 GRAPHS FOR ANGLE 13 FOR Z’ -AXIS FOR ANGLE 20 WRT AXIS FOR ANGLE 13 FOR X’ -AXIS 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
Recommend
More recommend
