tgrains.c - simple_DEM - a simple 2D Discrete Element Method code for educational purposes
(HTM) git clone git://src.adamsgaard.dk/simple_DEM
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tgrains.c (2840B)
---
1 #include <math.h>
2 #include "header.h"
3 #include "global_properties.h"
4
5
6 void prediction(grain* g)
7 {
8 int i;
9
10 #pragma omp parallel for shared(g) private (i)
11 for (i = 0; i < np; i++) {
12 /* Predict positions and velocities */
13 g[i].x += dt * g[i].vx + 0.5 * dt * dt * g[i].ax;
14 g[i].y += dt * g[i].vy + 0.5 * dt * dt * g[i].ay;
15 g[i].vx += 0.5 * dt * g[i].ax;
16 g[i].vy += 0.5 * dt * g[i].ay;
17
18 g[i].ang += dt * g[i].angv + 0.5 * dt * dt * g[i].anga;
19 g[i].angv += 0.5 * dt * g[i].anga;
20
21 /* Zero forces */
22 g[i].fx = 0.0;
23 g[i].fy = 0.0;
24 g[i].t = 0.0;
25 g[i].p = 0.0;
26 }
27 }
28
29 void interparticle_force(grain* g, int a, int b)
30 {
31 if (a > b) { /* Use Newtons 3rd law to find both forces at once */
32
33 /* Particle center coordinate component differences */
34 double x_ab = g[a].x - g[b].x;
35 double y_ab = g[a].y - g[b].y;
36
37 /* Particle center distance */
38 double dist = sqrt(x_ab*x_ab + y_ab*y_ab);
39
40 /* Size of overlap */
41 double dn = dist - (g[a].R + g[b].R);
42
43 if (dn < 0.0) { /* Contact */
44 double xn, yn, vn, fn; /* Normal components */
45 double xt, yt, vt, ft; /* Tangential components */
46 /* Local axes */
47 xn = x_ab / dist;
48 yn = y_ab / dist;
49 xt = -yn;
50 yt = xn;
51
52 /* Compute the velocity of the contact */
53 double vx_ab = g[a].vx - g[b].vx;
54 double vy_ab = g[a].vy - g[b].vy;
55 vn = vx_ab*xn + vy_ab*yn;
56 vt = vx_ab*xt + vy_ab*yt - (g[a].R*g[a].angv + g[b].R*g[b].angv);
57
58 /* Compute force in local axes */
59 fn = -kn * dn - nu * vn;
60
61 /* Rotation */
62 if (fn < 0)
63 fn = 0.0;
64 ft = fabs(kt * vt);
65 if (ft > mu*fn) /* Coefficient of friction */
66 ft = mu*fn;
67 if (vt > 0)
68 ft = -ft;
69
70 /* Calculate sum of forces on a and b in global coordinates */
71 g[a].fx += fn * xn;
72 g[a].fy += fn * yn;
73 g[a].t += -ft*g[a].R;
74 g[a].p += fn;
75 g[b].fx -= fn * xn;
76 g[b].fy -= fn * yn;
77 g[b].p += fn;
78 g[b].t += -ft*g[b].R;
79
80 }
81
82 } else {
83 return;
84 }
85 }
86
87 void interact_grains(grain* g)
88 {
89 int a, b;
90 #pragma omp parallel for shared(g) private (a,b)
91 /* Loop through particle a */
92 for (a = 0; a < np; a++) {
93
94 /* Loop through particle b */
95 for (b = 0; b < np; b++) {
96 interparticle_force(g, a, b);
97 }
98 }
99 }
100
101 void update_acc(grain* g)
102 {
103 int i;
104 #pragma omp parallel for shared(g) private (i)
105 for (i = 0; i < np; i++) {
106 g[i].ax = g[i].fx / g[i].m;
107 g[i].ay = g[i].fy / g[i].m - grav;
108 g[i].anga = g[i].t / g[i].I;
109 }
110 }
111
112 void correction(grain* g)
113 {
114 int i;
115 #pragma omp parallel for shared(g) private (i)
116 for (i = 0; i < np; i++) {
117 g[i].vx += 0.5 * dt * g[i].ax;
118 g[i].vy += 0.5 * dt * g[i].ay;
119 g[i].angv += 0.5 * dt * g[i].anga;
120 }
121 }
122