granular

granular dynamics simulation
git clone git://src.adamsgaard.dk/granular # fast
git clone https://src.adamsgaard.dk/granular.git # slow
Log | Files | Refs | README | LICENSE Back to index

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 }