tFirst commit - simple_DEM - a simple 2D Discrete Element Method code for educa… | |
git clone git://src.adamsgaard.dk/simple_DEM | |
Log | |
Files | |
Refs | |
LICENSE | |
--- | |
commit a17213846a2df2ba52051d8f9f9effcd8b58636d | |
Author: Anders Damsgaard Christensen <[email protected]> | |
Date: Sun, 14 Oct 2012 08:42:39 +0200 | |
First commit | |
Diffstat: | |
A Makefile | 19 +++++++++++++++++++ | |
A README.rst | 17 +++++++++++++++++ | |
A global_properties.h | 30 ++++++++++++++++++++++++++++++ | |
A grains.c | 119 +++++++++++++++++++++++++++++… | |
A header.h | 44 +++++++++++++++++++++++++++++… | |
A initialize.c | 28 ++++++++++++++++++++++++++++ | |
A main.c | 86 ++++++++++++++++++++++++++++++ | |
A output/blank.txt | 0 | |
A vtk_export.c | 134 +++++++++++++++++++++++++++++… | |
A walls.c | 109 +++++++++++++++++++++++++++++… | |
10 files changed, 586 insertions(+), 0 deletions(-) | |
--- | |
diff --git a/Makefile b/Makefile | |
t@@ -0,0 +1,19 @@ | |
+CC=g++ | |
+CFLAGS=-c -Wall -O2 | |
+LDFLAGS= | |
+#CFLAGS=-c -Wall -O2 -fopenmp -g -G | |
+#LDFLAGS=-fopenmp | |
+SOURCES=main.c | |
+OBJECTS=$(SOURCES:.c=.o) | |
+EXECUTABLE=simple_DEM | |
+ | |
+all: $(SOURCES) $(EXECUTABLE) | |
+ | |
+$(EXECUTABLE): $(OBJECTS) | |
+ $(CC) $(LDFLAGS) $(OBJECTS) -o $@ | |
+ | |
+.c.o: | |
+ $(CC) $(CFLAGS) $< -o $@ | |
+ | |
+clean: | |
+ rm -f $(EXECUTABLE) *.o output/* | |
diff --git a/README.rst b/README.rst | |
t@@ -0,0 +1,17 @@ | |
+simple_DEM | |
+========== | |
+ | |
+A basic CPU-based discrete element method algorithm, created for comparison ag… | |
+ | |
+Requirements | |
+------------ | |
+- GNU Make | |
+- GCC | |
+ | |
+Obtaining | |
+--------- | |
+`git clone https://github.com/anders-dc/simple_DEM.git` | |
+ | |
+Building | |
+-------- | |
+`make` | |
diff --git a/global_properties.h b/global_properties.h | |
t@@ -0,0 +1,30 @@ | |
+#ifndef GLOBAL_CONSTANTS_H_ | |
+#define GLOBAL_CONSTANTS_H_ | |
+ | |
+// Properties of sample | |
+const int ng = 5000; // Number of grains | |
+ | |
+// Size distribution | |
+const double rmin = 1e-3; // Min. radius | |
+const double rmax = 2e-3; // Max. radius | |
+ | |
+// Properties of grains | |
+const double kn = 1e5; // Normal stiffness | |
+const double nu = 30; // Normal damping | |
+const double rho = 1000; // Density of the grains | |
+const double mu = 0.5; // Sliding friction | |
+const double kt = kn; // Tangential stiffness | |
+ | |
+// Temporal variables | |
+const double dt = 1e-4; // Time step length | |
+const int maxStep = 3000; // Number of steps | |
+const int fileInterval = 20; // No. of steps between output | |
+ | |
+// Physical constants | |
+const double grav = 9.80; // Gravity | |
+ | |
+// Number of grains along the width in the initial configuration | |
+const int ngw = 100; | |
+ | |
+#endif | |
+ | |
diff --git a/grains.c b/grains.c | |
t@@ -0,0 +1,119 @@ | |
+void prediction(grain* g) | |
+{ | |
+ int i; | |
+ | |
+ #pragma omp parallel for shared(g) private (i) | |
+ for (i = 0; i < ng; i++) { | |
+ // Predict positions and velocities | |
+ g[i].x += dt * g[i].vx + 0.5 * dt * dt * g[i].ax; | |
+ g[i].y += dt * g[i].vy + 0.5 * dt * dt * g[i].ay; | |
+ g[i].vx += 0.5 * dt * g[i].ax; | |
+ g[i].vy += 0.5 * dt * g[i].ay; | |
+ | |
+ g[i].th += dt * g[i].vth + 0.5 * dt * dt * g[i].ath; | |
+ g[i].vth += 0.5 * dt * g[i].ath; | |
+ | |
+ // Zero forces | |
+ g[i].fx = 0.0; | |
+ g[i].fy = 0.0; | |
+ g[i].fth = 0.0; | |
+ g[i].p = 0.0; | |
+ } | |
+} | |
+ | |
+void interparticle_force(grain* g, int a, int b) | |
+{ | |
+ if (a > b) { // Use Newtons 3rd law to find both forces at once | |
+ | |
+ // Particle center coordinate component differences | |
+ double x_ab = g[a].x - g[b].x; | |
+ double y_ab = g[a].y - g[b].y; | |
+ | |
+ // Particle center distance | |
+ double dist = sqrt(x_ab*x_ab + y_ab*y_ab); | |
+ | |
+ // Size of overlap | |
+ double dn = dist - (g[a].R + g[b].R); | |
+ | |
+ if (dn < 0.0) { // Contact | |
+ double xn, yn, vn, fn; // Normal components | |
+ double xt, yt, vt, ft; // Tangential components | |
+ // Local axes | |
+ xn = x_ab / dist; | |
+ yn = y_ab / dist; | |
+ xt = -yn; | |
+ yt = xn; | |
+ | |
+ // Compute the velocity of the contact | |
+ double vx_ab = g[a].vx - g[b].vy; | |
+ double vy_ab = g[a].vy - g[b].vy; | |
+ vn = vx_ab*xn + vy_ab*yn; | |
+ vt = vx_ab*xt + vy_ab*yt - (g[a].R*g[a].vth + g[b].R*g[b].vth); | |
+ | |
+ // Compute force in local axes | |
+ fn = -kn * dn - nu * vn; | |
+ | |
+ // Rotation | |
+ if (fn < 0) | |
+ fn = 0.0; | |
+ ft = fabs(kt * vt); | |
+ if (ft > mu*fn) // Coefficient of friction | |
+ ft = mu*fn; | |
+ if (vt > 0) | |
+ ft = -ft; | |
+ | |
+ // Calculate sum of forces on a and b in global coordinates | |
+ g[a].fx += fn * xn; | |
+ g[a].fy += fn * yn; | |
+ g[a].fth += -ft*g[a].R; | |
+ g[a].p += fn; | |
+ g[b].fx -= fn * xn; | |
+ g[b].fy -= fn * yn; | |
+ g[b].p += fn; | |
+ g[b].fth += -ft*g[b].R; | |
+ | |
+ } | |
+ | |
+ } else { | |
+ return; | |
+ } | |
+} | |
+ | |
+void interact_grains(grain* g) | |
+{ | |
+ int a, b; | |
+ #pragma omp parallel for shared(g) private (a,b) | |
+ // Loop through particle a | |
+ for (a = 0; a < ng; a++) { | |
+ | |
+ // Loop through particle b | |
+ for (b = 0; b < ng; b++) { | |
+ interparticle_force(g, a, b); | |
+ } | |
+ | |
+ } | |
+ | |
+} | |
+ | |
+void update_acc(grain* g) | |
+{ | |
+ int i; | |
+ #pragma omp parallel for shared(g) private (i) | |
+ for (i = 0; i < ng; i++) { | |
+ g[i].ax = g[i].fx / g[i].m; | |
+ g[i].ay = g[i].fy / g[i].m - grav; | |
+ g[i].ath = g[i].fth / g[i].I; | |
+ } | |
+} | |
+ | |
+void correction(grain* g) | |
+{ | |
+ int i; | |
+ #pragma omp parallel for shared(g) private (i) | |
+ for (i = 0; i < ng; i++) { | |
+ g[i].vx += 0.5 * dt * g[i].ax; | |
+ g[i].vy += 0.5 * dt * g[i].ay; | |
+ g[i].vth += 0.5 * dt * g[i].ath; | |
+ } | |
+} | |
+ | |
diff --git a/header.h b/header.h | |
t@@ -0,0 +1,44 @@ | |
+#ifndef HEADER_H_ | |
+#define HEADER_H_ | |
+ | |
+//// Structure declarations | |
+struct grain | |
+{ | |
+ double m; // Mass | |
+ double R; // Radius | |
+ double I; // Inertia | |
+ double x, y, th; // Position | |
+ double vx, vy, vth; // Velocities | |
+ double ax, ay, ath; // Acceleration | |
+ double fx, fy, fth; // Sum of forces, decomposed | |
+ double p; // Pressure | |
+}; | |
+ | |
+ | |
+ | |
+//// Prototype functions | |
+ | |
+// initialize.cpp | |
+void triangular_grid(grain* g); | |
+ | |
+// grains.cpp | |
+void prediction(grain* g); | |
+void interparticle_force(grain* g, int a, int b); | |
+void interact_grains(grain* g); | |
+void update_acc(grain* g); | |
+void correction(grain* g); | |
+ | |
+// walls.cpp | |
+void compute_force_upper_wall(int i, grain* g, double wfy, double wupy); | |
+void compute_force_lower_wall(int i, grain* g, double wfy, double wloy); | |
+void compute_force_left_wall(int i, grain* g, double wfy, double wlex); | |
+void compute_force_right_wall(int i, grain* g, double wfy, double wrix); | |
+ | |
+ | |
+// vtk_export.cpp | |
+int vtk_export_grains(grain* g, int numfile); | |
+int vtk_export_walls(int numfile, double wlex, double wrix, double wloy, doubl… | |
+int vtk_export_forces(grain* g, int numfile); | |
+ | |
+ | |
+#endif | |
diff --git a/initialize.c b/initialize.c | |
t@@ -0,0 +1,28 @@ | |
+#include <stdlib.h> | |
+#include <math.h> | |
+#include "header.h" | |
+#include "global_properties.h" | |
+ | |
+void triangular_grid(grain* g) | |
+{ | |
+ // Initialize grain properties | |
+ for (int i = 0; i < ng; i++) { | |
+ g[i].R = (rand() / (double)RAND_MAX) * (rmax - rmin) + rmin; | |
+ g[i].m = M_PI * rho * g[i].R * g[i].R; | |
+ g[i].I = 0.5 * g[i].m * g[i].R * g[i].R; | |
+ g[i].p = 0.0; | |
+ } | |
+ | |
+ // Initialize grain positions in a triangular grid | |
+ for (int i = 0; i < ng; i++) { | |
+ int column = i%ngw; | |
+ int row = i/ngw; | |
+ | |
+ if (row%2 == 0) // Even row | |
+ g[i].x = rmax + 2*column*rmax; | |
+ else // Odd row | |
+ g[i].x = 2*rmax + 2*column*rmax; | |
+ g[i].y = rmax + 2*row*rmax; | |
+ | |
+ } | |
+} | |
diff --git a/main.c b/main.c | |
t@@ -0,0 +1,86 @@ | |
+#include <iostream> | |
+#include <cmath> | |
+ | |
+// Structure declarations and function prototypes | |
+#include "header.h" | |
+ | |
+// Global and constant simulation properties | |
+#include "global_properties.h" | |
+ | |
+// Functions for exporting data to VTK formats | |
+#include "vtk_export.c" | |
+ | |
+#include "initialize.c" | |
+#include "grains.c" | |
+#include "walls.c" | |
+ | |
+ | |
+ | |
+int main() | |
+{ | |
+ using std::cout; | |
+ using std::endl; | |
+ | |
+ cout << "\n## simple_DEM ##\n" | |
+ << "Particles: " << ng << endl | |
+ << "maxStep: " << maxStep << endl; | |
+ | |
+ | |
+ double time = 0.0; // Time at simulation start | |
+ | |
+ // Allocate memory | |
+ grain* g = new grain[ng]; // Grain structure | |
+ | |
+ | |
+ // Compute simulation domain dimensions | |
+ double wleft = 0.0; // Left wall | |
+ double wright = (ngw+1)*2*rmax; // Right wall | |
+ double wdown = 0.0; // Lower wall | |
+ double wup = (ng/ngw+1)*2*rmax; // Upper wall | |
+ | |
+ // Variables for pressures on walls | |
+ double wp_up, wp_down, wp_left, wp_right; | |
+ | |
+ // Initiailze grains inside the simulation domain | |
+ triangular_grid(g); | |
+ | |
+ | |
+ | |
+ // Main time loop | |
+ for (int step = 0; step < maxStep; step++) { | |
+ | |
+ time += dt; // Update current time | |
+ | |
+ // Predict new kinematics | |
+ prediction(g); | |
+ | |
+ // Calculate contact forces between grains | |
+ interact_grains(g); | |
+ | |
+ // Calculate contact forces between grains and walls | |
+ interact_walls(g, wleft, wright, wup, wdown, &wp_up, &wp_down, &wp_left, &… | |
+ | |
+ // Update acceleration from forces | |
+ update_acc(g); | |
+ | |
+ // Correct velocities | |
+ correction(g); | |
+ | |
+ // Write output files if the fileInterval is reached | |
+ if (step % fileInterval == 0) { | |
+ (void)vtk_export_grains(g, step); | |
+ (void)vtk_export_walls(step, wleft, wright, wdown, wup, wp_up, wp_down, … | |
+ (void)vtk_export_forces(g, step); | |
+ } | |
+ | |
+ } // End of main time loop | |
+ | |
+ | |
+ // Free dynamically allocated memory | |
+ delete[] g; | |
+ | |
+ | |
+ cout << "\nSimulation ended without errors.\n"; | |
+ return 0; // Terminate successfully | |
+} | |
+ | |
diff --git a/output/blank.txt b/output/blank.txt | |
diff --git a/vtk_export.c b/vtk_export.c | |
t@@ -0,0 +1,134 @@ | |
+#include <iostream> | |
+#include <cstdio> | |
+#include <fstream> | |
+#include "header.h" | |
+#include "global_properties.h" | |
+ | |
+int vtk_export_grains(grain* g, int numfile) | |
+{ | |
+ FILE* fout; | |
+ | |
+ char filename[25]; // File name | |
+ sprintf(filename, "output/grains%04d.vtk", numfile); | |
+ | |
+ if ((fout = fopen(filename, "wt")) == NULL) { | |
+ std::cout << "vtk_export error, cannot open " << filename << "!\n"; | |
+ return 1; | |
+ } | |
+ | |
+ // Write header | |
+ fprintf(fout, "# vtk DataFile Version 3.0\n"); | |
+ | |
+ // Write title | |
+ fprintf(fout, "'Time: %f s'\n", numfile*dt); | |
+ | |
+ // Write data type | |
+ fprintf(fout, "ASCII\nDATASET UNSTRUCTURED_GRID\n"); | |
+ | |
+ // Grain coordinates | |
+ fprintf(fout, "POINTS %d FLOAT\n", ng); | |
+ for (int i = 0; i < ng; i++) | |
+ fprintf(fout, "%f %f 0.0\n", g[i].x, g[i].y); | |
+ fprintf(fout, "POINT_DATA %d\n", ng); | |
+ | |
+ // Grain radii | |
+ fprintf(fout, "VECTORS Radius float\n"); | |
+ for (int i = 0; i < ng; i++) | |
+ fprintf(fout, "%f 0.0 0.0\n", g[i].R); | |
+ | |
+ // Grain velocities | |
+ fprintf(fout, "VECTORS Velocity float\n"); | |
+ for (int i = 0; i < ng; i++) | |
+ fprintf(fout, "%f %f 0.0\n", g[i].vx, g[i].vy); | |
+ | |
+ // Pressure | |
+ fprintf(fout, "SCALARS Pressure float 1\n"); | |
+ fprintf(fout, "LOOKUP_TABLE default\n"); | |
+ for (int i = 0; i < ng; i++) | |
+ fprintf(fout, "%e\n", g[i].p); | |
+ | |
+ // Angular velocity | |
+ fprintf(fout, "SCALARS Angvel float 1\n"); | |
+ fprintf(fout, "LOOKUP_TABLE default\n"); | |
+ for (int i = 0; i < ng; i++) | |
+ fprintf(fout, "%e\n", g[i].vth); | |
+ | |
+ fclose(fout); | |
+ return 0; | |
+} | |
+ | |
+// Save walls vtk file | |
+int vtk_export_walls(int numfile, double wlex, double wrix, double wloy, doubl… | |
+ double wp_up, double wp_down, double wp_left, double … | |
+{ | |
+ char fname[25]; // file name | |
+ sprintf(fname,"output/walls%04d.vtk",numfile); | |
+ | |
+ using std::endl; | |
+ | |
+ std::ofstream fow(fname, std::ios::out); | |
+ if (fow) | |
+ { | |
+ fow.precision(3); fow << std::scientific; | |
+ fow << "# vtk DataFile Version 3.0" << endl; | |
+ fow << "My walls" << endl; | |
+ fow << "ASCII" << endl; | |
+ fow << "DATASET POLYDATA" << endl; | |
+ fow << "POINTS 8 float" << endl; | |
+ // lower wall | |
+ fow << wlex << " " << wloy << " 0" << endl; | |
+ fow << wrix << " " << wloy << " 0" << endl; | |
+ // upper wall | |
+ fow << wlex << " " << wupy << " 0" << endl; | |
+ fow << wrix << " " << wupy << " 0" << endl; | |
+ // left wall | |
+ fow << wlex << " " << wloy << " 0" << endl; | |
+ fow << wlex << " " << wupy << " 0" << endl; | |
+ // right wall | |
+ fow << wrix << " " << wloy << " 0" << endl; | |
+ fow << wrix << " " << wupy << " 0" << endl; | |
+ fow << "LINES 4 12" << endl; | |
+ fow << "2 0 1" << endl; | |
+ fow << "2 2 3" << endl; | |
+ fow << "2 4 5" << endl; | |
+ fow << "2 6 7" << endl; | |
+ fow << "POINT_DATA 8" << endl; | |
+ fow << "SCALARS Rayon float" << endl; | |
+ fow << "LOOKUP_TABLE default" << endl; | |
+ fow << wp_up << endl; // No idea about the following sequence, … | |
+ fow << wp_up << endl; | |
+ fow << wp_left << endl; | |
+ fow << wp_left << endl; | |
+ fow << wp_right << endl; | |
+ fow << wp_right << endl; | |
+ fow << wp_down << endl; | |
+ fow << wp_down << endl; // ?? | |
+ } | |
+ return 0; | |
+} | |
+ | |
+int vtk_export_forces(grain* g, int numfile) | |
+{ | |
+ FILE* fout; | |
+ | |
+ char filename[25]; // File name | |
+ sprintf(filename, "output/forces%04d.vtk", numfile); | |
+ | |
+ if ((fout = fopen(filename, "wt")) == NULL) { | |
+ std::cout << "vtk_export error, cannot open " << filename << "!\n"; | |
+ return 1; | |
+ } | |
+ | |
+ // Write header | |
+ fprintf(fout, "# vtk DataFile Version 3.0\n"); | |
+ | |
+ // Write title | |
+ fprintf(fout, "'Time: %f s'\n", numfile*dt); | |
+ | |
+ // Write data type | |
+ fprintf(fout, "ASCII\nDATASET POLYDATA\n"); | |
+ | |
+ | |
+ fclose(fout); | |
+ return 0; | |
+} | |
diff --git a/walls.c b/walls.c | |
t@@ -0,0 +1,109 @@ | |
+#include <cmath> | |
+#include "header.h" | |
+#include "global_properties.h" | |
+ | |
+// Compute force between i and upper wall | |
+double compute_force_upper_wall(int i, grain* g, double wup) | |
+{ | |
+ double dn = wup - (g[i].y + g[i].R); // Overlap | |
+ | |
+ if(dn<0.0) { | |
+ double vn,fn; | |
+ // velocity (wall velocity = 0) | |
+ vn = g[i].vy; | |
+ // force | |
+ fn = -kn * dn - nu * vn; | |
+ // Update sum of forces on grains | |
+ g[i].fy -= fn; | |
+ // Add fn to pressure on grains i | |
+ g[i].p += fn; | |
+ // Update stress to the wall | |
+ return fn; | |
+ } | |
+ return 0.0; | |
+} | |
+ | |
+double compute_force_lower_wall(int i, grain* g, double wdown) | |
+{ | |
+ | |
+ double dn = g[i].y - g[i].R - wdown; // Overlap | |
+ | |
+ if(dn<0.0) { | |
+ double vn,fn; | |
+ // velocity (wall velocity = 0) | |
+ vn = g[i].vy; | |
+ // force | |
+ fn = - kn * dn - nu * vn; | |
+ // Update sum of forces on grains | |
+ g[i].fy += fn; | |
+ // Add fn to pressure on grains i | |
+ g[i].p += fn; | |
+ // Update stress to the wall | |
+ return fn; | |
+ } | |
+ return 0.0; | |
+} | |
+ | |
+double compute_force_left_wall(int i, grain* g, double wleft) | |
+{ | |
+ double dn = wleft + g[i].x - g[i].R; // Overlap | |
+ | |
+ if(dn<0.0) { | |
+ double vn,fn; | |
+ // velocity (wall velocity = 0) | |
+ vn = g[i].vx; | |
+ // force | |
+ fn = -kn * dn - nu * vn; | |
+ // Update sum of forces on grains | |
+ g[i].fx += fn; | |
+ // Add fn to pressure on grains i | |
+ g[i].p += fn; | |
+ // Update stress to the wall | |
+ return fn; | |
+ } | |
+ return 0.0; | |
+} | |
+ | |
+ | |
+double compute_force_right_wall(int i, grain* g, double wright) | |
+{ | |
+ double dn = wright - (g[i].x + g[i].R); // Overlap | |
+ | |
+ if(dn<0.0) { | |
+ double vn,fn; | |
+ // velocity (wall velocity = 0) | |
+ vn = -g[i].vx; | |
+ // force | |
+ fn = -kn * dn - nu * vn; | |
+ // Update sum of forces on grains | |
+ g[i].fx -= fn; | |
+ // Add fn to pressure on grains i | |
+ g[i].p += fn; | |
+ // Update stress to the wall | |
+ return fn; | |
+ } | |
+ return 0.0; | |
+} | |
+ | |
+void interact_walls(grain* g, double wleft, double wright, double wup, double … | |
+ double* wp_up, double* wp_down, double* wp_left, doubl… | |
+{ | |
+ double u = 0.0; | |
+ double d = 0.0; | |
+ double r = 0.0; | |
+ double l = 0.0; | |
+ | |
+ int i; | |
+ | |
+ #pragma omp parallel for shared(g, u, d, r, l) private (i) | |
+ for (i = 0; i < ng; i++) { | |
+ u += compute_force_upper_wall(i, g, wup); | |
+ d += compute_force_lower_wall(i, g, wdown); | |
+ r += compute_force_right_wall(i, g, wright); | |
+ l += compute_force_left_wall(i, g, wleft); | |
+ } | |
+ *wp_up = u; | |
+ *wp_down = d; | |
+ *wp_left = l; | |
+ *wp_right = r; | |
+} |