tadd acceleration as constant body force - granular - granular dynamics simulation (HTM) git clone git://src.adamsgaard.dk/granular (DIR) Log (DIR) Files (DIR) Refs (DIR) README (DIR) LICENSE --- (DIR) commit cce359bcaf94938276bdd3e28515f088b2ceef5f (DIR) parent 0dc426c8b11d38d85d984ece057b53911b74330a (HTM) Author: Anders Damsgaard <anders@adamsgaard.dk> Date: Thu, 22 Apr 2021 16:13:17 +0200 add acceleration as constant body force Diffstat: M grain.c | 136 ++++++++++++++++---------------- M granular.c | 4 +++- M simulation.c | 21 +++++++++++++++------ M simulation.h | 4 +--- 4 files changed, 87 insertions(+), 78 deletions(-) --- (DIR) diff --git a/grain.c b/grain.c t@@ -21,24 +21,24 @@ grain_new(void) void grain_defaults(struct grain *g) { - int i; + int d; g->diameter = 1.0; - for (i = 0; i < 3; i++) { - g->pos[i] - = g->vel[i] - = g->acc[i] - = g->force[i] - = g->angpos[i] - = g->angvel[i] - = g->angacc[i] - = g->torque[i] - = g->disp[i] - = g->forceext[i] - = g->contact_stress[i] + for (d = 0; d < 3; d++) { + g->pos[d] + = g->vel[d] + = g->acc[d] + = g->force[d] + = g->angpos[d] + = g->angvel[d] + = g->angacc[d] + = g->torque[d] + = g->disp[d] + = g->forceext[d] + = g->contact_stress[d] = 0.0; - g->gridpos[i] - = g->acc_lock[i] = 0; + g->gridpos[d] + = g->acc_lock[d] = 0; } g->density = 2.5e3; g->fixed = 0; t@@ -53,8 +53,8 @@ grain_defaults(struct grain *g) g->shear_strength = 1e8; g->fracture_toughness = 1e8; g->ncontacts = 0; - for (i = 0; i < MAXCONTACTS; i++) - contact_defaults(&g->contacts[i]); + for (d = 0; d < MAXCONTACTS; d++) + contact_defaults(&g->contacts[d]); g->thermal_energy = 0.0; g->color = 0; } t@@ -62,28 +62,28 @@ grain_defaults(struct grain *g) static void print_padded_nd_double(FILE *stream, const double *arr) { - int i; + int d; - for (i = 0; i < 3; i++) - fprintf(stream, "%.*g\t", FLOATPREC, arr[i]); + for (d = 0; d < 3; d++) + fprintf(stream, "%.*g\t", FLOATPREC, arr[d]); } static void print_padded_nd_int(FILE *stream, const int *arr) { - int i; + int d; - for (i = 0; i < 3; i++) - fprintf(stream, "%d\t", arr[i]); + for (d = 0; d < 3; d++) + fprintf(stream, "%d\t", arr[d]); } static void print_padded_nd_uint(FILE *stream, const size_t *arr) { - int i; + int d; - for (i = 0; i < 3; i++) - fprintf(stream, "%zu\t", arr[i]); + for (d = 0; d < 3; d++) + fprintf(stream, "%zu\t", arr[d]); } void t@@ -197,24 +197,24 @@ grain_read(char *line) int grain_check_values(const struct grain *g) { - int i; + int d; int status = 0; check_float_positive("grain->diameter", g->diameter, &status); - for (i = 0; i < 3; i++) { - check_float("grain->pos", g->pos[i], &status); - check_float("grain->vel", g->vel[i], &status); - check_float("grain->acc", g->acc[i], &status); - check_int_bool("grain->acc_lock", g->acc_lock[i], &status); - check_float("grain->force", g->force[i], &status); - check_float("grain->angpos", g->angpos[i], &status); - check_float("grain->angvel", g->angvel[i], &status); - check_float("grain->angacc", g->angacc[i], &status); - check_float("grain->torque", g->torque[i], &status); - check_float("grain->disp", g->disp[i], &status); - check_float("grain->forceext", g->forceext[i], &status); - check_float("grain->contact_stress", g->contact_stress[i], &status); + for (d = 0; d < 3; d++) { + check_float("grain->pos", g->pos[d], &status); + check_float("grain->vel", g->vel[d], &status); + check_float("grain->acc", g->acc[d], &status); + check_int_bool("grain->acc_lock", g->acc_lock[d], &status); + check_float("grain->force", g->force[d], &status); + check_float("grain->angpos", g->angpos[d], &status); + check_float("grain->angvel", g->angvel[d], &status); + check_float("grain->angacc", g->angacc[d], &status); + check_float("grain->torque", g->torque[d], &status); + check_float("grain->disp", g->disp[d], &status); + check_float("grain->forceext", g->forceext[d], &status); + check_float("grain->contact_stress", g->contact_stress[d], &status); } check_float_non_negative("grain->density", g->density, &status); t@@ -247,10 +247,10 @@ grain_check_values(const struct grain *g) g->fracture_toughness, &status); - for (i = 0; i < 3; i++) - if (g->gridpos[i] > 1) + for (d = 0; d < 3; d++) + if (g->gridpos[d] > 1) warn_parameter_value("grain->gridpos is not 0 or 1", - (double)g->gridpos[i], &status); + (double)g->gridpos[d], &status); check_float_non_negative("grain->thermal_energy", g->thermal_energy, &status); t@@ -278,16 +278,16 @@ grain_moment_of_inertia(const struct grain *g) void grain_zero_kinematics(struct grain *g) { - int i; - - for (i = 0; i < 3; i++) { - g->vel[i] = 0.0; - g->acc[i] = 0.0; - g->force[i] = 0.0; - g->angvel[i] = 0.0; - g->angacc[i] = 0.0; - g->torque[i] = 0.0; - g->disp[i] = 0.0; + int d; + + for (d = 0; d < 3; d++) { + g->vel[d] = 0.0; + g->acc[d] = 0.0; + g->force[d] = 0.0; + g->angvel[d] = 0.0; + g->angacc[d] = 0.0; + g->torque[d] = 0.0; + g->disp[d] = 0.0; } } t@@ -313,29 +313,29 @@ grain_kinetic_energy(const struct grain *g) static void grain_temporal_integration_two_term_taylor(struct grain *g, double dt) { - int i; + int d; double dx, mass = grain_mass(g); double moment_of_inertia = grain_moment_of_inertia(g); - for (i = 0; i < 3; i++) { - if (!g->acc_lock[i]) - g->acc[i] = g->force[i] / mass; + for (d = 0; d < 3; d++) { + if (!g->acc_lock[d]) + g->acc[d] = (g->force[d] + g->forceext[d]) / mass; if (g->rotating) - g->angacc[i] = g->torque[i] / moment_of_inertia; + g->angacc[d] = g->torque[d] / moment_of_inertia; } if (g->fixed) - for (i = 0; i < 3; i++) - g->acc[i] = 0.0; - - for (i = 0; i < 3; i++) { - dx = g->vel[i] * dt + 0.5 * g->acc[i] * dt * dt; - g->pos[i] += dx; - g->disp[i] += dx; - g->angpos[i] += g->angvel[i] * dt + 0.5 * g->angacc[i] * dt * dt; - g->vel[i] += g->acc[i] * dt; - g->angvel[i] += g->angacc[i] * dt; + for (d = 0; d < 3; d++) + g->acc[d] = 0.0; + + for (d = 0; d < 3; d++) { + dx = g->vel[d] * dt + 0.5 * g->acc[d] * dt * dt; + g->pos[d] += dx; + g->disp[d] += dx; + g->angpos[d] += g->angvel[d] * dt + 0.5 * g->angacc[d] * dt * dt; + g->vel[d] += g->acc[d] * dt; + g->angvel[d] += g->angacc[d] * dt; } } (DIR) diff --git a/granular.c b/granular.c t@@ -24,6 +24,7 @@ main(int argc, char *argv[]) { int ret; struct simulation sim = sim_new(); + double grav = 0.0; #ifdef __OpenBSD__ if (pledge("stdio wpath cpath", NULL) == -1) t@@ -37,7 +38,7 @@ main(int argc, char *argv[]) sim.t_end = atof(EARGF(usage())); break; case 'g': - sim.constacc[1] = atof(EARGF(usage())); + grav = atof(EARGF(usage())); break; case 'h': usage(); t@@ -63,6 +64,7 @@ main(int argc, char *argv[]) usage(); sim_read_grains(&sim, stdin); + sim_add_acceleration_scalar(&sim, grav, 1); sim_set_timestep(&sim); if (sim.t < sim.t_end) sim_run_time_loop(&sim); (DIR) diff --git a/simulation.c b/simulation.c t@@ -24,17 +24,16 @@ sim_new(void) void sim_defaults(struct simulation *sim) { - int i, ret; + int d, ret; ret = snprintf(sim->name, sizeof(sim->name), DEFAULT_SIMULATION_NAME); if (ret < 0 || (size_t)ret >= sizeof(sim->name)) errx(1, "%s: could not write simulation name", __func__); sim->ng = 0; - for (i = 0; i < 3; i++) { - sim->constacc[i] = 0.0; - sim->nd[i] = 1; - sim->origo[i] = 0.0; - sim->L[i] = 1.0; + for (d = 0; d < 3; d++) { + sim->nd[d] = 1; + sim->origo[d] = 0.0; + sim->L[d] = 1.0; } sim->t = 0.0; sim->t_end = 0.0; t@@ -144,6 +143,16 @@ sim_check_add_contact(struct simulation *sim, size_t i, size_t j) grain_register_contact(&sim->grains[i], j, i, centerdist, overlap); } +void +sim_add_acceleration_scalar(struct simulation *sim, double acc, int dimension) +{ + size_t i; + + for (i = 0; i < sim->ng; i++) + sim->grains[i].forceext[dimension] + += acc * grain_mass(&sim->grains[i]); +} + /* Silbert et al. 2001 */ void sim_set_timestep(struct simulation *sim) (DIR) diff --git a/simulation.h b/simulation.h t@@ -12,9 +12,6 @@ struct simulation { /* simulation name to use for output files */ char name[255]; - /* constant acceleration [m/s^2] */ - double constacc[3]; - /* grain sorting grid */ size_t nd[3]; double origo[3]; t@@ -45,6 +42,7 @@ void sim_print_grains(FILE *stream, const struct simulation *sim); void sim_print_grains_vtk(FILE *stream, const struct simulation *sim); void sim_write_output_files(struct simulation *sim); +void sim_add_acceleration_scalar(struct simulation *sim, double acc, int dimension); void sim_set_timestep(struct simulation *sim); void sim_detect_contacts(struct simulation *sim); void sim_step_time(struct simulation *sim);