grains.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 }