tgrains.c - granular - granular dynamics simulation
(HTM) git clone git://src.adamsgaard.dk/granular
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
tgrains.c (7803B)
---
1 #include <stdio.h>
2 #include <math.h>
3 #include "grain.h"
4 #include "arrays.h"
5
6 #define VTK_FLOAT_FMT "%.17g "
7
8 #define VTK_XML_SCALAR(M, N, T, F) \
9 fprintf(stream,\
10 "\t\t\t\t<DataArray type=\"" T "\" Name=\"" N "\" "\
11 "NumberOfComponents=\"1\" format=\"ascii\">\n");\
12 for (i = 0; i < ng; i++)\
13 fprintf(stream, F, grains[i].M);\
14 fprintf(stream, "\n\t\t\t\t</DataArray>\n");
15
16 #define VTK_XML_VECTOR(M, N, T, F) \
17 fprintf(stream,\
18 "\t\t\t\t<DataArray type=\"" T "\" Name=\"" N "\" "\
19 "NumberOfComponents=\"3\" format=\"ascii\">\n");\
20 for (i = 0; i < ng; i++)\
21 for (d = 0; d < 3; d++)\
22 fprintf(stream, F, grains[i].M[d]);\
23 fprintf(stream, "\n\t\t\t\t</DataArray>\n");
24
25 void
26 grains_print(FILE *stream, const struct grain *grains, size_t ng)
27 {
28 size_t i;
29
30 for (i = 0; i < ng; i++)
31 grain_print(stream, &grains[i]);
32 }
33
34 void
35 grains_print_vtk(FILE *stream, const struct grain *grains, size_t ng)
36 {
37 size_t i, d;
38
39 fprintf(stream,
40 "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n"
41 "<VTKFile type=\"UnstructuredGrid\" version=\"1.0\" "
42 "byte_order=\"LittleEndian\">\n"
43 "\t<UnstructuredGrid>\n"
44 "\t\t<Piece NumberOfPoints=\"%zu\" NumberOfCells=\"0\">\n", ng);
45 fprintf(stream, "\t\t\t<Points>\n");
46 VTK_XML_VECTOR(pos, "Position [m]", "Float64", VTK_FLOAT_FMT);
47 fprintf(stream, "\t\t\t</Points>\n");
48
49 fprintf(stream,
50 "\t\t\t<Cells>\n"
51 "\t\t\t\t<DataArray type=\"Int32\" Name=\"connectivity\" "
52 "NumberOfComponents=\"1\" format=\"ascii\"/>\n"
53 "\t\t\t\t<DataArray type=\"Int32\" Name=\"offsets\" "
54 "NumberOfComponents=\"1\" format=\"ascii\"/>\n"
55 "\t\t\t\t<DataArray type=\"UInt8\" Name=\"types\" "
56 "NumberOfComponents=\"1\" format=\"ascii\"/>\n"
57 "\t\t\t</Cells>\n");
58
59 fprintf(stream,
60 "\t\t\t<PointData Scalars=\"Diameter [m]\" "
61 "Vectors=\"Angular position [-]\">\n");
62
63 VTK_XML_SCALAR(diameter, "Diameter [m]", "Float64", VTK_FLOAT_FMT);
64 VTK_XML_VECTOR(vel, "Velocity [m/s]", "Float64", VTK_FLOAT_FMT);
65 VTK_XML_VECTOR(acc, "Acceleration [m/s^2]", "Float64", VTK_FLOAT_FMT);
66 VTK_XML_VECTOR(acc_lock, "Acceleration lock [-]", "Int64", "%d ");
67 VTK_XML_VECTOR(force, "Force [N]", "Float64", VTK_FLOAT_FMT);
68 VTK_XML_VECTOR(angpos, "Angular position [rad]", "Float64", VTK_FLOAT_FMT);
69 VTK_XML_VECTOR(angvel, "Angular velocity [rad/s]", "Float64", VTK_FLOAT_FMT);
70 VTK_XML_VECTOR(angacc, "Angular acceleration [rad/s^2]", "Float64", VTK_FLOAT_FMT);
71 VTK_XML_VECTOR(torque, "Torque [N/m]", "Float64", VTK_FLOAT_FMT);
72 VTK_XML_VECTOR(disp, "Displacement [m]", "Float64", VTK_FLOAT_FMT);
73 VTK_XML_VECTOR(forceext, "External body force [N]", "Float64", VTK_FLOAT_FMT);
74 VTK_XML_SCALAR(density, "Density [kg/m^3]", "Float64", VTK_FLOAT_FMT);
75 VTK_XML_SCALAR(fixed, "Fixed [-]", "Int64", "%d ");
76 VTK_XML_SCALAR(rotating, "Rotating [-]", "Int64", "%d ");
77 VTK_XML_SCALAR(enabled, "Enabled [-]", "Int64", "%d ");
78 VTK_XML_SCALAR(youngs_modulus, "Young's modulus [Pa]", "Float64", VTK_FLOAT_FMT);
79 VTK_XML_SCALAR(poissons_ratio, "Poisson's ratio [-]", "Float64", VTK_FLOAT_FMT);
80 VTK_XML_SCALAR(friction_coeff, "Friction coefficient [-]", "Float64", VTK_FLOAT_FMT);
81 VTK_XML_SCALAR(damping_n, "Contact-normal damping constant [1/s]", "Float64", VTK_FLOAT_FMT);
82 VTK_XML_SCALAR(damping_t, "Contact-tangential damping constant [1/s]", "Float64", VTK_FLOAT_FMT);
83 VTK_XML_SCALAR(tensile_strength, "Tensile strength [Pa]", "Float64", VTK_FLOAT_FMT);
84 VTK_XML_SCALAR(shear_strength, "Shear strength [Pa]", "Float64", VTK_FLOAT_FMT);
85 VTK_XML_SCALAR(fracture_toughness, "Fracture toughness [Pa]", "Float64", VTK_FLOAT_FMT);
86 VTK_XML_VECTOR(gridpos, "Grid position [-]", "UInt64", "%zu ");
87 VTK_XML_SCALAR(ncontacts, "Number of contacts [-]", "UInt64", "%zu ");
88 VTK_XML_VECTOR(contact_stress, "Contact stress [Pa]", "Float64", VTK_FLOAT_FMT);
89 VTK_XML_SCALAR(thermal_energy, "Thermal energy [J]", "Float64", VTK_FLOAT_FMT);
90 VTK_XML_SCALAR(color, "Color [-]", "Int64", "%d ");
91
92 fprintf(stream, "\t\t\t</PointData>\n");
93
94 fprintf(stream,
95 "\t\t</Piece>\n"
96 "\t</UnstructuredGrid>\n"
97 "</VTKFile>\n");
98 }
99
100 void
101 grains_print_energy(FILE *stream, const struct grain *grains, size_t ng)
102 {
103 size_t i;
104 double E_kin = 0.0,
105 E_rot = 0.0,
106 E_thermal = 0.0;
107
108 for (i = 0; i < ng; i++) {
109 E_kin += grain_translational_kinetic_energy(&grains[i]);
110 E_rot += grain_rotational_kinetic_energy(&grains[i]);
111 E_thermal += grain_thermal_energy(&grains[i]);
112 }
113
114 fprintf(stream, "%.17g\t%.17g\t%.17g\t%.17g\n",
115 E_kin + E_rot + E_thermal,
116 E_kin, E_rot, E_thermal);
117 }
118
119 double
120 grains_minimum_grain_edge(const struct grain *grains, size_t ng, int d)
121 {
122 size_t i;
123 double res = INFINITY;
124
125 for (i = 0; i < ng; i++)
126 if (res > grains[i].pos[d] - grains[i].diameter / 2.0)
127 res = grains[i].pos[d] - grains[i].diameter / 2.0;
128
129 return res;
130 }
131
132 double
133 grains_maximum_grain_edge(const struct grain *grains, size_t ng, int d)
134 {
135 size_t i;
136 double res = -INFINITY;
137
138 for (i = 0; i < ng; i++)
139 if (res < grains[i].pos[d] + grains[i].diameter / 2.0)
140 res = grains[i].pos[d] + grains[i].diameter / 2.0;
141
142 return res;
143 }
144
145 /* Hertz-Mindlin contact model for compressible spheres */
146 void
147 grains_interact(struct grain *g_i, struct grain *g_j, int ic, double dt)
148 {
149 int d;
150 double f_n_ij[3], f_t_ij[3],
151 v_ij[3], v_n_ij[3], v_t_ij[3], v_n_dot_n_ij,
152 delta_ij,
153 n_ij[3],
154 d_i, d_j,
155 r_ij[3], r_ij_norm,
156 m_i, m_j, m_ij,
157 mu_ij,
158 gamma_n_ij, gamma_t_ij,
159 k_n_ij, k_t_ij, A_ij,
160 angvel_ij[3], u_t_dot_v_ij,
161 *angvel_cross_r_ij, *f_t_cross_r_ij;
162
163 delta_ij = g_i->contacts[ic].overlap;
164
165 if (delta_ij < 0.0) { /* TODO: implement tensile strength */
166 g_i->contacts[ic].active = 0;
167 g_i->ncontacts--;
168 return;
169 }
170
171 d_i = g_i->diameter;
172 d_j = g_j->diameter;
173 m_i = grain_mass(g_i);
174 m_j = grain_mass(g_j);
175 m_ij = m_i * m_j / (m_i + m_j);
176 r_ij_norm = euclidean_norm(g_i->contacts[ic].centerdist, 3);
177
178 for (d = 0; d < 3; d++) {
179 r_ij[d] = g_i->contacts[ic].centerdist[d];
180 n_ij[d] = r_ij[d]/r_ij_norm;
181 v_ij[d] = g_i->vel[d] - g_j->vel[d];
182 angvel_ij[d] = g_i->angvel[d] + g_j->angvel[d];
183 }
184
185 v_n_dot_n_ij = dot(v_ij, n_ij, 3);
186 angvel_cross_r_ij = cross(angvel_ij, r_ij);
187
188 for (d = 0; d < 3; d++) {
189 v_n_ij[d] = v_n_dot_n_ij * n_ij[d];
190 v_t_ij[d] = v_ij[d] - v_n_ij[d] - 0.5 * angvel_cross_r_ij[d];
191 }
192
193 mu_ij = fmin(g_i->friction_coeff, g_j->friction_coeff);
194 gamma_n_ij = fmin(g_i->damping_n, g_j->damping_n);
195 gamma_t_ij = fmin(g_i->damping_t, g_j->damping_t);
196
197 k_n_ij = fmin(grain_stiffness_normal(g_i),
198 grain_stiffness_normal(g_j));
199 k_t_ij = fmin(grain_stiffness_tangential(g_i),
200 grain_stiffness_tangential(g_j));
201
202 A_ij = sqrt(delta_ij) * sqrt(d_i * d_j / (2.0 * (d_i + d_j)));
203
204 for (d = 0; d < 3; d++) {
205 f_n_ij[d] = A_ij * (k_n_ij * delta_ij * n_ij[d]
206 - m_ij * gamma_n_ij * v_n_ij[d]);
207 f_t_ij[d] = A_ij * (-k_t_ij * g_i->contacts[ic].tandisp[d]
208 - m_ij * gamma_t_ij * v_t_ij[d]);
209 }
210
211 if (euclidean_norm(f_t_ij, 3) >= euclidean_norm(f_n_ij, 3) * mu_ij)
212 return; /* TODO: limit f_t_ij according to Coulomb criterion */
213
214 f_t_cross_r_ij = cross(f_t_ij, r_ij);
215
216 for (d = 0; d < 3; d++) {
217 g_i->force[d] += f_n_ij[d] + f_t_ij[d];
218 g_j->force[d] -= f_n_ij[d] + f_t_ij[d];
219 g_i->torque[d] += -0.5 * f_t_cross_r_ij[d];
220 g_j->torque[d] -= -0.5 * f_t_cross_r_ij[d];
221 }
222
223 g_i->contacts[ic].age += dt;
224 u_t_dot_v_ij = dot(g_i->contacts[ic].tandisp, v_ij, 3);
225 for (d = 0; d < 3; d++)
226 g_i->contacts[ic].tandisp[d] += (v_t_ij[d]
227 - u_t_dot_v_ij * r_ij[d]
228 / (r_ij_norm * r_ij_norm))
229 * dt;
230
231 free(angvel_cross_r_ij);
232 free(f_t_cross_r_ij);
233 }