tgrain.c - granular - granular dynamics simulation
(HTM) git clone git://src.adamsgaard.dk/granular
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
tgrain.c (11985B)
---
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4 #include <err.h>
5 #include "granular.h"
6 #include "util.h"
7 #include "grain.h"
8 #include "arrays.h"
9 #include "contact.h"
10
11 #define FLOATPREC 17
12
13 struct grain
14 grain_new(void)
15 {
16 struct grain g;
17 grain_defaults(&g);
18 return g;
19 }
20
21 void
22 grain_defaults(struct grain *g)
23 {
24 int d;
25
26 g->diameter = 1.0;
27 for (d = 0; d < 3; d++) {
28 g->pos[d]
29 = g->vel[d]
30 = g->acc[d]
31 = g->force[d]
32 = g->angpos[d]
33 = g->angvel[d]
34 = g->angacc[d]
35 = g->torque[d]
36 = g->disp[d]
37 = g->forceext[d]
38 = g->contact_stress[d]
39 = 0.0;
40 g->gridpos[d]
41 = g->acc_lock[d] = 0;
42 }
43 g->density = 2.5e3;
44 g->fixed = 0;
45 g->rotating = 1;
46 g->enabled = 1;
47 g->youngs_modulus = 1e9;
48 g->poissons_ratio = 0.2;
49 g->friction_coeff = 0.4;
50 g->damping_n = 0.0;
51 g->damping_t = 0.0;
52 g->tensile_strength = 1e8;
53 g->shear_strength = 1e8;
54 g->fracture_toughness = 1e8;
55 g->ncontacts = 0;
56 for (d = 0; d < MAXCONTACTS; d++)
57 contact_defaults(&g->contacts[d]);
58 g->thermal_energy = 0.0;
59 g->color = 0;
60 }
61
62 static void
63 print_padded_nd_double(FILE *stream, const double *arr)
64 {
65 int d;
66
67 for (d = 0; d < 3; d++)
68 fprintf(stream, "%.*g\t", FLOATPREC, arr[d]);
69 }
70
71 static void
72 print_padded_nd_int(FILE *stream, const int *arr)
73 {
74 int d;
75
76 for (d = 0; d < 3; d++)
77 fprintf(stream, "%d\t", arr[d]);
78 }
79
80 static void
81 print_padded_nd_uint(FILE *stream, const size_t *arr)
82 {
83 int d;
84
85 for (d = 0; d < 3; d++)
86 fprintf(stream, "%zu\t", arr[d]);
87 }
88
89 void
90 grain_print(FILE *stream, const struct grain *g)
91 {
92 fprintf(stream, "%.*g\t", FLOATPREC, g->diameter);
93 print_padded_nd_double(stream, g->pos);
94 print_padded_nd_double(stream, g->vel);
95 print_padded_nd_double(stream, g->acc);
96 print_padded_nd_int(stream, g->acc_lock);
97 print_padded_nd_double(stream, g->force);
98 print_padded_nd_double(stream, g->angpos);
99 print_padded_nd_double(stream, g->angvel);
100 print_padded_nd_double(stream, g->angacc);
101 print_padded_nd_double(stream, g->torque);
102 print_padded_nd_double(stream, g->disp);
103 print_padded_nd_double(stream, g->forceext);
104 fprintf(stream, "%.*g\t", FLOATPREC, g->density);
105 fprintf(stream, "%d\t", g->fixed);
106 fprintf(stream, "%d\t", g->rotating);
107 fprintf(stream, "%d\t", g->enabled);
108 fprintf(stream, "%.*g\t", FLOATPREC, g->youngs_modulus);
109 fprintf(stream, "%.*g\t", FLOATPREC, g->poissons_ratio);
110 fprintf(stream, "%.*g\t", FLOATPREC, g->friction_coeff);
111 fprintf(stream, "%.*g\t", FLOATPREC, g->damping_n);
112 fprintf(stream, "%.*g\t", FLOATPREC, g->damping_t);
113 fprintf(stream, "%.*g\t", FLOATPREC, g->tensile_strength);
114 fprintf(stream, "%.*g\t", FLOATPREC, g->shear_strength);
115 fprintf(stream, "%.*g\t", FLOATPREC, g->fracture_toughness);
116 print_padded_nd_uint(stream, g->gridpos);
117 fprintf(stream, "%zu\t", g->ncontacts);
118 print_padded_nd_double(stream, g->contact_stress);
119 fprintf(stream, "%.*g\t", FLOATPREC, g->thermal_energy);
120 fprintf(stream, "%d\n", g->color);
121 }
122
123 struct grain *
124 grain_read(char *line)
125 {
126 struct grain *g;
127
128 if (!(g = malloc(sizeof(struct grain))))
129 err(1, "%s: grain malloc", __func__);
130
131 if (sscanf(line, "%lg\t" /* diameter */
132 "%lg\t%lg\t%lg\t" /* pos */
133 "%lg\t%lg\t%lg\t" /* vel */
134 "%lg\t%lg\t%lg\t" /* acc */
135 "%d\t%d\t%d\t" /* acc_lock */
136 "%lg\t%lg\t%lg\t" /* force */
137 "%lg\t%lg\t%lg\t" /* angpos */
138 "%lg\t%lg\t%lg\t" /* angvel */
139 "%lg\t%lg\t%lg\t" /* angacc */
140 "%lg\t%lg\t%lg\t" /* torque */
141 "%lg\t%lg\t%lg\t" /* disp */
142 "%lg\t%lg\t%lg\t" /* forceext */
143 "%lg\t" /* density */
144 "%d\t" /* fixed */
145 "%d\t" /* rotating */
146 "%d\t" /* enabled */
147 "%lg\t" /* youngs_modulus */
148 "%lg\t" /* poissons_ratio */
149 "%lg\t" /* friction_coeff */
150 "%lg\t" /* damping_n */
151 "%lg\t" /* damping_t */
152 "%lg\t" /* tensile_strength */
153 "%lg\t" /* shear_strength */
154 "%lg\t" /* fracture_toughness */
155 "%zu\t%zu\t%zu\t" /* gridpos */
156 "%zu\t" /* ncontacts */
157 "%lg\t%lg\t%lg\t" /* contact_stress */
158 "%lg\t" /* thermal_energy */
159 "%d\n", /* color */
160 &g->diameter,
161 &g->pos[0], &g->pos[1], &g->pos[2],
162 &g->vel[0], &g->vel[1], &g->vel[2],
163 &g->acc[0], &g->acc[1], &g->acc[2],
164 &g->acc_lock[0], &g->acc_lock[1], &g->acc_lock[2],
165 &g->force[0], &g->force[1], &g->force[2],
166 &g->angpos[0], &g->angpos[1], &g->angpos[2],
167 &g->angvel[0], &g->angvel[1], &g->angvel[2],
168 &g->angacc[0], &g->angacc[1], &g->angacc[2],
169 &g->torque[0], &g->torque[1], &g->torque[2],
170 &g->disp[0], &g->disp[1], &g->disp[2],
171 &g->forceext[0], &g->forceext[1], &g->forceext[2],
172 &g->density,
173 &g->fixed,
174 &g->rotating,
175 &g->enabled,
176 &g->youngs_modulus,
177 &g->poissons_ratio,
178 &g->friction_coeff,
179 &g->damping_n,
180 &g->damping_t,
181 &g->tensile_strength,
182 &g->shear_strength,
183 &g->fracture_toughness,
184 &g->gridpos[0], &g->gridpos[1], &g->gridpos[2],
185 &g->ncontacts,
186 &g->contact_stress[0], &g->contact_stress[1], &g->contact_stress[2],
187 &g->thermal_energy,
188 &g->color) != 55)
189 errx(1, "%s: could not read line: %s", __func__, line);
190
191 if (grain_check_values(g))
192 errx(1, "%s: invalid values in grain line: %s", __func__, line);
193
194 return g;
195 }
196
197 int
198 grain_check_values(const struct grain *g)
199 {
200 int d;
201 int status = 0;
202
203 check_float_positive("grain->diameter", g->diameter, &status);
204
205 for (d = 0; d < 3; d++) {
206 check_float("grain->pos", g->pos[d], &status);
207 check_float("grain->vel", g->vel[d], &status);
208 check_float("grain->acc", g->acc[d], &status);
209 check_int_bool("grain->acc_lock", g->acc_lock[d], &status);
210 check_float("grain->force", g->force[d], &status);
211 check_float("grain->angpos", g->angpos[d], &status);
212 check_float("grain->angvel", g->angvel[d], &status);
213 check_float("grain->angacc", g->angacc[d], &status);
214 check_float("grain->torque", g->torque[d], &status);
215 check_float("grain->disp", g->disp[d], &status);
216 check_float("grain->forceext", g->forceext[d], &status);
217 check_float("grain->contact_stress", g->contact_stress[d], &status);
218 }
219
220 check_float_non_negative("grain->density", g->density, &status);
221 check_int_bool("grain->fixed", g->fixed, &status);
222 check_int_bool("grain->rotating", g->rotating, &status);
223 check_int_bool("grain->enabled", g->enabled, &status);
224
225 check_float_non_negative("grain->youngs_modulus",
226 g->youngs_modulus,
227 &status);
228 check_float_non_negative("grain->poissons_ratio",
229 g->poissons_ratio,
230 &status);
231 check_float_non_negative("grain->friction_coeff",
232 g->friction_coeff,
233 &status);
234 check_float_non_negative("grain->damping_n",
235 g->damping_n,
236 &status);
237 check_float_non_negative("grain->damping_t",
238 g->damping_t,
239 &status);
240 check_float_non_negative("grain->tensile_strength",
241 g->tensile_strength,
242 &status);
243 check_float_non_negative("grain->shear_strength",
244 g->shear_strength,
245 &status);
246 check_float_non_negative("grain->fracture_toughness",
247 g->fracture_toughness,
248 &status);
249
250 for (d = 0; d < 3; d++)
251 if (g->gridpos[d] > 1)
252 warn_parameter_value("grain->gridpos is not 0 or 1",
253 (double)g->gridpos[d], &status);
254
255 check_float_non_negative("grain->thermal_energy", g->thermal_energy, &status);
256
257 return status;
258 }
259
260 double
261 grain_volume(const struct grain *g)
262 {
263 return 4.0 / 3.0 * PI * pow(g->diameter / 2.0, 3.0);
264 }
265
266 double
267 grain_mass(const struct grain *g)
268 {
269 return grain_volume(g) * g->density;
270 }
271
272 double
273 grain_moment_of_inertia(const struct grain *g)
274 {
275 return 2.0 / 5.0 * grain_mass(g) * pow(g->diameter / 2.0, 2.0);
276 }
277
278 void
279 grain_zero_kinematics(struct grain *g)
280 {
281 int d;
282
283 for (d = 0; d < 3; d++) {
284 g->vel[d] = 0.0;
285 g->acc[d] = 0.0;
286 g->force[d] = 0.0;
287 g->angvel[d] = 0.0;
288 g->angacc[d] = 0.0;
289 g->torque[d] = 0.0;
290 g->disp[d] = 0.0;
291 }
292 }
293
294 void
295 grain_force_reset(struct grain *g)
296 {
297 int d;
298
299 for (d = 0; d < 3; d++)
300 g->force[d] = 0.0;
301 }
302
303 double
304 grain_translational_kinetic_energy(const struct grain *g)
305 {
306 return 0.5 * grain_mass(g) * euclidean_norm(g->vel, 3.0);
307 }
308
309 double
310 grain_rotational_kinetic_energy(const struct grain *g)
311 {
312 return 0.5 * grain_moment_of_inertia(g) * euclidean_norm(g->angvel, 3.0);
313 }
314
315 double
316 grain_kinetic_energy(const struct grain *g)
317 {
318 return grain_translational_kinetic_energy(g) +
319 grain_rotational_kinetic_energy(g);
320 }
321
322 double
323 grain_thermal_energy(const struct grain *g)
324 {
325 return g->thermal_energy;
326 }
327
328 static void
329 grain_temporal_integration_two_term_taylor(struct grain *g, double dt)
330 {
331 int d;
332 double dx, mass = grain_mass(g);
333 double moment_of_inertia = grain_moment_of_inertia(g);
334
335 for (d = 0; d < 3; d++) {
336 g->acc[d] = (g->force[d] + g->forceext[d]) / mass;
337
338 if (g->rotating)
339 g->angacc[d] = g->torque[d] / moment_of_inertia;
340 }
341
342 if (g->fixed)
343 for (d = 0; d < 3; d++) {
344 g->angacc[d] = 0.0;
345 if (!g->acc_lock[d])
346 g->acc[d] = 0.0;
347 }
348
349 for (d = 0; d < 3; d++) {
350 dx = g->vel[d] * dt + 0.5 * g->acc[d] * dt * dt;
351 g->pos[d] += dx;
352 g->disp[d] += dx;
353 g->angpos[d] += g->angvel[d] * dt + 0.5 * g->angacc[d] * dt * dt;
354 g->vel[d] += g->acc[d] * dt;
355 g->angvel[d] += g->angacc[d] * dt;
356 }
357 }
358
359 void
360 grain_temporal_integration(struct grain *g, double dt)
361 {
362 grain_temporal_integration_two_term_taylor(g, dt);
363 }
364
365 void
366 grain_register_contact(struct grain *g, size_t i, size_t j,
367 double centerdist[3], double overlap)
368 {
369 int ic, d;
370
371 /* first pass: check if contact is already registered */
372 for (ic = 0; ic < MAXCONTACTS; ic++)
373 if (g->contacts[ic].active && g->contacts[ic].j == j) {
374 g->contacts[ic].overlap = overlap;
375 for (d = 0; d < 3; d++)
376 g->contacts[ic].centerdist[d] = centerdist[d];
377 return;
378 }
379
380 /* second pass: register as new contact if overlapping */
381 if (overlap > 0.0)
382 for (ic = 0; ic <= MAXCONTACTS; ic++) {
383 if (ic == MAXCONTACTS)
384 err(1, "%s: contact %zu exceeds MAXCONTACTS (%d)",
385 __func__, i, MAXCONTACTS);
386 if (!g->contacts[ic].active) {
387 contact_new(&g->contacts[ic], i, j);
388 g->ncontacts += 1;
389 g->contacts[ic].overlap = overlap;
390 for (d = 0; d < 3; d++) {
391 g->contacts[ic].centerdist[d] = centerdist[d];
392 g->contacts[ic].tandisp[d] = 0.0;
393 }
394 return;
395 }
396 }
397 }
398
399 double
400 grain_stiffness_normal(const struct grain *g)
401 {
402 return 2.0 / 3.0 * g->youngs_modulus / (1.0 - pow(g->poissons_ratio, 2.0));
403 }
404
405 double
406 grain_stiffness_tangential(const struct grain *g)
407 {
408 return 2.0 * g->youngs_modulus
409 / ((1.0 + g->poissons_ratio) * (2.0 - g->poissons_ratio));
410 }
411
412 double
413 grain_collision_time(const struct grain *g)
414 {
415 double m;
416
417 m = grain_mass(g);
418
419 return fmin(PI * pow(2.0 * grain_stiffness_normal(g) / m
420 - pow(g->damping_n, 2.0) / 4.0, -0.5),
421 PI * pow(2.0 * grain_stiffness_tangential(g) / m
422 - pow(g->damping_t, 2.0) / 4.0, -0.5));
423 }