tcohesion.cuh - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
HTML git clone git://src.adamsgaard.dk/sphere
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tcohesion.cuh (11619B)
---
1 #ifndef COHESION_CUH_
2 #define COHESION_CUH_
3
4 // cohesion.cuh
5 // Functions governing attractive forces between contacts
6
7 // Check bond pair list, apply linear contact model to pairs
8 __global__ void bondsLinear(
9 uint2* __restrict__ dev_bonds,
10 Float4* __restrict__ dev_bonds_delta, // Contact displacement
11 Float4* __restrict__ dev_bonds_omega, // Contact rotational displacement
12 const Float4* __restrict__ dev_x,
13 const Float4* __restrict__ dev_vel,
14 const Float4* __restrict__ dev_angvel,
15 Float4* __restrict__ dev_force,
16 Float4* __restrict__ dev_torque)
17 {
18 // Find thread index
19 unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x;
20
21 if (idx < devC_params.nb0) {
22
23 //// Read values
24
25 // Read bond data
26 __syncthreads();
27 const uint2 bond = dev_bonds[idx]; // Particle indexes in bond pair
28
29 // Check if the bond has been erased
30 if (bond.x < devC_np) {
31
32 const Float4 delta0_4 = dev_bonds_delta[idx];
33 const Float4 omega0_4 = dev_bonds_omega[idx];
34
35 // Convert tangential vectors to Float3's
36 // Uncorrected tangential component of displacement
37 Float3 delta0_t = MAKE_FLOAT3(
38 delta0_4.x,
39 delta0_4.y,
40 delta0_4.z);
41 const Float delta0_n = delta0_4.w;
42
43 // Uncorrected tangential component of rotation
44 Float3 omega0_t = MAKE_FLOAT3(
45 omega0_4.x,
46 omega0_4.y,
47 omega0_4.z);
48 const Float omega0_n = omega0_4.w;
49
50 // Read particle data
51 const Float4 x_i = dev_x[bond.x];
52 const Float4 x_j = dev_x[bond.y];
53 const Float4 vel_i = dev_vel[bond.x];
54 const Float4 vel_j = dev_vel[bond.y];
55 const Float4 angvel4_i = dev_angvel[bond.x];
56 const Float4 angvel4_j = dev_angvel[bond.y];
57
58 const Float3 angvel_i = MAKE_FLOAT3(angvel4_i.x, angvel4_i.y, angvel4_i.z);
59 const Float3 angvel_j = MAKE_FLOAT3(angvel4_j.x, angvel4_j.y, angvel4_j.z);
60
61
62 //// Bond geometry and inertia
63
64 // Parallel-bond radius (Potyondy and Cundall 2004, eq. 12)
65 const Float R_bar = devC_params.lambda_bar * fmin(x_i.w, x_j.w);
66
67 // Bond cross section area (Potyondy and Cundall 2004, eq. 15)
68 const Float A = PI * R_bar*R_bar;
69
70 // Bond moment of inertia (Potyondy and Cundall 2004, eq. 15)
71 const Float I = 0.25 * PI * R_bar*R_bar*R_bar*R_bar;
72
73 // Bond polar moment of inertia (Potyondy and Cundall 2004, eq. 15)
74 //const Float J = 0.50 * PI * R_bar*R_bar*R_bar*R_bar;
75 const Float J = I*2.0;
76
77 // Inter-particle vector
78 const Float3 x = MAKE_FLOAT3(
79 x_i.x - x_j.x,
80 x_i.y - x_j.y,
81 x_i.z - x_j.z);
82 const Float x_length = length(x);
83
84 // Find overlap (negative value if overlapping)
85 const Float overlap = fmin(0.0, x_length - (x_i.w + x_j.w));
86
87 // Normal vector of contact (points from i to j)
88 Float3 n = x/x_length;
89
90
91 //// Force
92
93 // Correct tangential displacement vector for rotation of the contact plane
94 //const Float3 delta_t0 = delta_t0_uncor - dot(delta_t0_uncor, n);
95 delta0_t = delta0_t - (n * dot(n, delta0_t));
96
97 // Contact displacement, Luding 2008 eq. 10
98 const Float3 ddelta = (
99 MAKE_FLOAT3(
100 vel_i.x - vel_j.x,
101 vel_i.y - vel_j.y,
102 vel_i.z - vel_j.z)
103 + (x_i.w + overlap/2.0) * cross(n, angvel_i)
104 + (x_j.w + overlap/2.0) * cross(n, angvel_j)
105 ) * devC_dt;
106
107 // Normal component of the displacement increment
108 //const Float ddelta_n = dot(ddelta, n);
109 const Float ddelta_n = -dot(ddelta, n);
110
111 // Normal component of the total displacement
112 const Float delta_n = delta0_n + ddelta_n;
113
114 // Tangential component of the displacement increment
115 // Luding 2008, eq. 9
116 const Float3 ddelta_t = ddelta - n * dot(n, ddelta);
117
118 // Tangential component of the total displacement
119 const Float3 delta_t = delta0_t + ddelta_t;
120
121 // Normal force: Visco-elastic contact model
122 // The elastic component caused by the overlap is subtracted.
123 //f_n = devC_params.k_n * A * delta_n * n;
124 const Float3 f_n = (devC_params.k_n * A * delta_n + devC_params.gamma_n * ddelta_n/devC_dt) * n;
125 //f_n += devC_params.k_n * overlap * n;
126
127
128 // Tangential force: Visco-elastic contact model
129 //f_t = -devC_params.k_t * A * delta_t;
130 const Float3 f_t = -devC_params.k_t * A * delta_t - devC_params.gamma_t * ddelta_t/devC_dt;
131
132 // Force vector
133 const Float3 f = f_n + f_t;
134
135
136 //// Torque
137
138 // Correct tangential rotational vector for rotation of the contact plane
139 omega0_t = omega0_t - (-n) * dot(omega0_t, -n);
140 //omega0_t = omega0_t - (n * dot(n, omega0_t));
141
142 // Contact rotational velocity
143 Float3 domega = MAKE_FLOAT3(
144 angvel_j.x - angvel_i.x,
145 angvel_j.y - angvel_i.y,
146 angvel_j.z - angvel_i.z) * devC_dt;
147 /*const Float3 domega = MAKE_FLOAT3(
148 angvel_i.x - angvel_j.x,
149 angvel_i.y - angvel_j.y,
150 angvel_i.z - angvel_j.z) * devC_dt;*/
151
152 // Normal component of the rotational increment
153 const Float domega_n = dot(domega, -n);
154 //const Float domega_n = dot(-n, domega);
155 //const Float domega_n = -dot(n, domega);
156
157 // Normal component of the total rotation
158 const Float omega_n = omega0_n + domega_n;
159
160 // Tangential component of the rotational increment
161 //const Float3 domega_t = domega - (-n) * dot(domega, -n);
162 const Float3 domega_t = domega - domega_n * (-n);
163 //const Float3 domega_t = domega - n * dot(n, domega);
164
165 // Tangential component of the total rotation
166 const Float3 omega_t = omega0_t + domega_t;
167
168 // Twisting torque: Visco-elastic contact model
169 //const Float3 t_n = -devC_params.k_t * J * omega_n * n;
170 const Float3 t_n = -devC_params.k_t * J * omega_n * -n;
171 //t_n = devC_params.k_t * J * omega_n * n;
172 //t_n = (devC_params.k_t * J * omega_n + devC_params.gamma_t * domega_n/devC_dt) * n;
173
174 // Bending torque: Visco-elastic contact model
175 //t_t = -devC_params.k_n * I * omega_t;
176 //const Float3 t_t = devC_params.k_n * I * omega_t;
177 const Float3 t_t = -devC_params.k_n * I * omega_t;
178 //t_t = -devC_params.k_n * I * omega_t - devC_params.gamma_n * domega_t/devC_dt;
179 //t_t = devC_params.k_n * I * omega_t - devC_params.gamma_n * domega_t/devC_dt;
180
181 // Torque vector
182 //t = t_n + t_t;
183 //Float3 t_i = t_n + cross(-(x_i.w + overlap*0.5) * n, t_t);
184 //Float3 t_j = t_n + cross(-(x_j.w + overlap*0.5) * n, t_t);
185 //const Float3 t_i = t_n + cross(-(x_i.w + overlap*0.5) * n, f_t + t_t);
186 //const Float3 t_j = t_n + cross(-(x_j.w + overlap*0.5) * n, f_t - t_t);
187 const Float3 t_j = t_n + t_t;
188 //const Float3 t_i = t_n - t_t; //t_n - t_t;
189 //const Float3 t_j = t_n + t_t;
190
191
192 //// Bond strength (Potyondy & Cundall 2004)
193 // Extensions of Euler-Bernoulli beam bending theory
194 // Max. stresses in bond periphery
195
196 // Tensile stress
197 const Float sigma_max = length(f_n) / A + length(t_t) * R_bar / I;
198
199 // Shear stress
200 const Float tau_max = length(f_t) / A + length(t_n) * R_bar / J;
201
202 // Break bond if tensile and shear stresses exceed strengths
203 if (sigma_max >= devC_params.sigma_b || tau_max >= devC_params.tau_b) {
204 __syncthreads();
205 dev_bonds[idx].x = devC_params.nb0;
206 return;
207 }
208
209
210
211 //// Save values
212 __syncthreads();
213
214 // Save updated displacements in global memory
215 dev_bonds_delta[idx] = MAKE_FLOAT4(delta_t.x, delta_t.y, delta_t.z, delta_n);
216 dev_bonds_omega[idx] = MAKE_FLOAT4(omega_t.x, omega_t.y, omega_t.z, omega_n);
217
218 // Save forces and torques to the particle pairs
219 // !!! This is probably wrong, see Obermayer et al. 2013, C & GT (49)
220 dev_force[bond.x] += MAKE_FLOAT4(f.x, f.y, f.z, 0.0);
221 dev_force[bond.y] -= MAKE_FLOAT4(f.x, f.y, f.z, 0.0);
222 //dev_torque[bond.x] += MAKE_FLOAT4(t.x, t.y, t.z, 0.0);
223 //dev_torque[bond.y] += MAKE_FLOAT4(t.x, t.y, t.z, 0.0);
224 //dev_torque[bond.x] += MAKE_FLOAT4(t_i.x, t_i.y, t_i.z, 0.0);
225 dev_torque[bond.x] -= MAKE_FLOAT4(t_j.x, t_j.y, t_j.z, 0.0);
226 dev_torque[bond.y] += MAKE_FLOAT4(t_j.x, t_j.y, t_j.z, 0.0);
227 //dev_torque[bond.y] += MAKE_FLOAT4(t_j.x, t_j.y, t_j.z, 0.0);
228 //dev_torque[bond.y] -= MAKE_FLOAT4(t_j.x, t_j.y, t_j.z, 0.0);
229 //dev_torque[bond.y] -= MAKE_FLOAT4(t.x, t.y, t.z, 0.0);
230 // make sure to remove write conflicts
231 }
232 }
233 }
234
235
236 // Capillary cohesion after Richefeu et al. (2006)
237 __device__ void capillaryCohesion_exp(
238 Float3* N,
239 const Float radius_a,
240 const Float radius_b,
241 const Float delta_ab,
242 const Float3 x_ab,
243 const Float x_ab_length,
244 const Float kappa)
245 {
246
247 // Normal vector
248 Float3 n_ab = x_ab/x_ab_length;
249
250 Float3 f_c;
251 Float lambda, R_geo, R_har, r, h;
252
253 // Determine the ratio; r = max{Ri/Rj;Rj/Ri}
254 if ((radius_a/radius_b) > (radius_b/radius_a))
255 r = radius_a/radius_b;
256 else
257 r = radius_b/radius_a;
258
259 // Exponential decay function
260 h = -sqrtf(r);
261
262 // The harmonic mean
263 R_har = (2.0f * radius_a * radius_b) / (radius_a + radius_b);
264
265 // The geometrical mean
266 R_geo = sqrtf(radius_a * radius_b);
267
268 // The exponential falloff of the capillary force with distance
269 lambda = 0.9f * h * sqrtf(devC_params.V_b/R_har);
270
271 // Calculate cohesional force
272 f_c = -kappa * R_geo * expf(-delta_ab/lambda) * n_ab;
273
274 // Add force components from this collision to total force for particle
275 *N += f_c;
276
277 } // End of capillaryCohesion_exp
278
279
280 // Capillary cohesion after Richefeu et al. (2008)
281 __device__ void capillaryCohesion2_exp(
282 Float3* N,
283 const Float radius_a,
284 const Float radius_b,
285 const Float delta_ab,
286 const Float3 x_ab,
287 const Float x_ab_length,
288 const Float kappa)
289 {
290
291 // Normal vector
292 const Float3 n_ab = x_ab/x_ab_length;
293
294 // Determine the ratio; r = max{Ri/Rj;Rj/Ri}
295 Float r;
296 if ((radius_a/radius_b) > (radius_b/radius_a))
297 r = radius_a/radius_b;
298 else
299 r = radius_b/radius_a;
300
301 const Float lambda = 0.9/1.4142135623730951 * pow(devC_params.V_b, 2.0)
302 * pow(r, -0.5) * pow(1.0/radius_a + 1.0/radius_b, 0.5);
303
304 // Calculate cohesional force
305 const Float3 f_c =
306 -kappa * sqrtf(radius_a*radius_b) * expf(-delta_ab/lambda) * n_ab;
307
308 // Add force components from this collision to total force for particle
309 *N += f_c;
310 }
311
312 #endif
313 // vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4