1 /*============================================================================
2 * Methods for particle adhesion forces
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 /*============================================================================
28 * Functions dealing with lagrangian adhesion forces
29 *============================================================================*/
30
31 #include "cs_defs.h"
32
33 /*----------------------------------------------------------------------------
34 * Standard C library headers
35 *----------------------------------------------------------------------------*/
36
37 #include <limits.h>
38 #include <stdio.h>
39 #include <stddef.h>
40 #include <stdlib.h>
41 #include <string.h>
42 #include <math.h>
43 #include <ctype.h>
44 #include <float.h>
45 #include <assert.h>
46
47 /*----------------------------------------------------------------------------
48 * Local headers
49 *----------------------------------------------------------------------------*/
50
51 #include "cs_base.h"
52 #include "cs_math.h"
53
54 #include "bft_mem.h"
55 #include "bft_error.h"
56
57 #include "cs_physical_constants.h"
58 #include "cs_random.h"
59
60 #include "cs_lagr.h"
61 #include "cs_lagr_tracking.h"
62 #include "cs_lagr_roughness.h"
63 #include "cs_lagr_clogging.h"
64
65 /*----------------------------------------------------------------------------
66 * Header for the current file
67 *----------------------------------------------------------------------------*/
68
69 #include "cs_lagr_adh.h"
70
71 /*----------------------------------------------------------------------------*/
72
73 BEGIN_C_DECLS
74
75 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
76
77 /*============================================================================
78 * Local macro declarations
79 *============================================================================*/
80
81 /*============================================================================
82 * Local structure declarations
83 *============================================================================*/
84
85 /*============================================================================
86 * Static global variables
87 *============================================================================*/
88
89 /* Cut-off distance for adhesion forces (assumed to be the Born distance) */
90 static const cs_real_t _d_cut_off = 1.65e-10;
91
92 /* Boltzmann constant */
93 static const double _k_boltz = 1.38e-23;
94
95 /* Free space permittivity */
96 static const cs_real_t _free_space_permit = 8.854e-12;
97
98 /* Faraday constant */
99 static const cs_real_t _faraday_cst = 9.648e4;
100
101 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
102
103 /*============================================================================
104 * Private function definitions
105 *============================================================================*/
106
107 /*----------------------------------------------------------------------------*/
108 /* VDW interaction between a sphere and a plate */
109 /*----------------------------------------------------------------------------*/
110
111 static void
_vdwsp(cs_real_t cstham,cs_real_t lambda_vdw,cs_real_t * distp,cs_real_t * rpart,cs_real_t * var)112 _vdwsp(cs_real_t cstham,
113 cs_real_t lambda_vdw,
114 cs_real_t *distp,
115 cs_real_t *rpart,
116 cs_real_t *var)
117 {
118 if (*distp < lambda_vdw / 2.0 / cs_math_pi)
119 *var = - cstham * *rpart / (6.0 * *distp)
120 * (1.0 / (1.0 + 14.0 * *distp / lambda_vdw + 5.0 * cs_math_pi / 4.9 * pow (*distp, 3.0)
121 / lambda_vdw / pow (*rpart, 2.0)));
122 else
123 *var = cstham * (2.45 / 60.0 / cs_math_pi * lambda_vdw
124 * ( (*distp - *rpart) / pow (*distp, 2)
125 - (*distp + 3.0 * *rpart) / pow ((*distp + 2.0 * *rpart), 2.0))
126 - 2.17 / 720.0 / pow (cs_math_pi, 2) * pow (lambda_vdw, 2.0)
127 * ((*distp - 2.0* *rpart) / pow (*distp, 3.0) - (*distp + 4.0 * *rpart) / pow ((*distp + 2.0 * *rpart), 3.0))
128 + 0.59 / 5040.0 / pow (cs_math_pi, 3.0) * pow (lambda_vdw, 3.0)
129 * ((*distp - 3.0 * *rpart) / pow (*distp, 4.0) - (*distp + 5.0 * *rpart) / pow ((*distp + 2.0 * *rpart), 4.0)));
130
131 }
132
133 /*----------------------------------------------------------------------------*/
134 /* VDW interaction between two spheres */
135 /*----------------------------------------------------------------------------*/
136
137 static void
_vdwss(cs_real_t cstham,cs_real_t lambda_vdw,cs_real_t * distcc,cs_real_t * rpart1,cs_real_t * rpart2,cs_real_t * var)138 _vdwss(cs_real_t cstham,
139 cs_real_t lambda_vdw,
140 cs_real_t *distcc,
141 cs_real_t *rpart1,
142 cs_real_t *rpart2,
143 cs_real_t *var)
144 {
145
146 *var = - cstham * *rpart1 * *rpart2
147 / (6.0 * (*distcc - *rpart1 - *rpart2) * (*rpart1 + *rpart2))
148 * ( 1.0
149 - 5.32 * (*distcc - *rpart1 - *rpart2) / lambda_vdw
150 * log (1.0 + lambda_vdw / (*distcc - *rpart1 - *rpart2) / 5.32));
151
152 }
153
154 /*----------------------------------------------------------------------------*/
155 /*----------------------------------------------------------------------------*/
156
157 static void
_edlsp(cs_real_t * distp,cs_real_t * rpart,cs_real_t phi_p,cs_real_t phi_s,cs_real_t tempf,cs_real_t * var)158 _edlsp(cs_real_t *distp,
159 cs_real_t *rpart,
160 cs_real_t phi_p,
161 cs_real_t phi_s,
162 cs_real_t tempf,
163 cs_real_t *var)
164 {
165 cs_lagr_physico_chemical_t *lag_pc = cs_glob_lagr_physico_chemical;
166 cs_real_t charge = 1.6e-19;
167
168 cs_real_t ldebye = pow( 2000.0 * pow (_faraday_cst,2) * lag_pc->fion
169 / (lag_pc->epseau * _free_space_permit * cs_physical_constants_r * tempf), -0.5);
170
171 /* Reduced zeta potential */
172 cs_real_t lphi1 = lag_pc->valen * charge * phi_p / _k_boltz / tempf;
173 cs_real_t lphi2 = lag_pc->valen * charge * phi_s / _k_boltz / tempf;
174
175 /* xtended reduced zeta potential */
176 /* (following the work from Ohshima et al, 1982, JCIS, 90, 17-26) */
177 /* or the particle */
178
179 cs_real_t tau = *rpart / ldebye;
180 lphi1 = 8.0 * tanh (lphi1 / 4.0)
181 / (1.0 + sqrt (1.0 - (2.0 * tau + 1.0) / pow(tau + 1.0, 2.0) * pow(tanh(lphi1 / 4.0), 2.0)));
182
183 /* For the plate */
184
185 lphi2 = 4.0 * tanh (lphi2 / 4.0);
186
187 /* alculation for the EDL force */
188
189 cs_real_t alpha = sqrt ((*distp + *rpart) / *rpart) + sqrt (*rpart / (*distp + *rpart));
190
191 cs_real_t omega1 = pow (lphi1, 2.0) + pow (lphi2, 2.0) + alpha * lphi1 * lphi2;
192
193 cs_real_t omega2 = pow (lphi1, 2.0) + pow (lphi2, 2.0) - alpha * lphi1 * lphi2;
194
195 cs_real_t gamma = sqrt (*rpart / (*distp + *rpart)) * exp ( -1.0 / ldebye * *distp);
196
197 *var = 2.0 * cs_math_pi * lag_pc->epseau * _free_space_permit
198 * pow(_k_boltz * tempf / charge, 2.0) * *rpart
199 * (*distp + *rpart) / (*distp + 2.0 * *rpart)
200 * (omega1 * log (1.0 + gamma) + omega2 * log (1.0 - gamma));
201
202 }
203
204 /*----------------------------------------------------------------------------*/
205 /*----------------------------------------------------------------------------*/
206
207 static void
_edlss(cs_real_t * distcc,cs_real_t * rpart1,cs_real_t * rpart2,cs_real_t phi1,cs_real_t phi2,cs_real_t tempf,cs_real_t * var)208 _edlss(cs_real_t *distcc,
209 cs_real_t *rpart1,
210 cs_real_t *rpart2,
211 cs_real_t phi1,
212 cs_real_t phi2,
213 cs_real_t tempf,
214 cs_real_t *var)
215 {
216 cs_lagr_physico_chemical_t *lag_pc = cs_glob_lagr_physico_chemical;
217 cs_real_t charge = 1.6e-19;
218
219 cs_real_t ldebye = pow( 2000.0 * pow (_faraday_cst,2) * lag_pc->fion
220 / (lag_pc->epseau * _free_space_permit * cs_physical_constants_r * tempf), -0.5);
221
222 /* Reduced zeta potential */
223
224 cs_real_t lphi1 = lag_pc->valen * charge * phi1 / _k_boltz / tempf;
225 cs_real_t lphi2 = lag_pc->valen * charge * phi2 / _k_boltz / tempf;
226
227 /* xtended reduced zeta potential */
228 /* (following the work from Ohshima et al, 1982, JCIS, 90, 17-26) */
229 /* For the first particle */
230
231 cs_real_t tau = *rpart1 / ldebye;
232 lphi1 = 8.0 * tanh (lphi1 / 4.0)
233 / (1.0 + sqrt (1.0 - (2.0 * tau + 1.0) / pow(tau + 1, 2.0) * pow(tanh(lphi1 / 4.0), 2.0)));
234
235 /* For the second particle */
236
237 tau = *rpart2 / ldebye;
238 lphi2 = 8.0 * tanh (lphi2 / 4.0)
239 / (1.0 + sqrt (1.0 - (2.0 * tau + 1.0) / pow ((tau + 1.0), 2.0) * pow (tanh (lphi2 / 4.0), 2.0)));
240
241 /* Calculation of the EDL force */
242
243 cs_real_t alpha = sqrt(*rpart2 * (*distcc - *rpart2) / (*rpart1 * (*distcc - *rpart1)))
244 + sqrt(*rpart1 * (*distcc - *rpart1) / (*rpart2 * (*distcc - *rpart2)));
245
246 cs_real_t omega1 = pow (lphi1, 2.0) + pow (lphi2, 2.0) + alpha * lphi1 * lphi2;
247
248 cs_real_t omega2 = pow (lphi1, 2.0) + pow (lphi2, 2.0) - alpha * lphi1 * lphi2;
249
250 cs_real_t gamma = sqrt (*rpart1 * *rpart2 / (*distcc - *rpart1) / (*distcc - *rpart2))
251 * exp (1.0 / ldebye * (*rpart1 + *rpart2 - *distcc));
252
253 *var = 2.0 * cs_math_pi * lag_pc->epseau * _free_space_permit
254 * pow ((_k_boltz * tempf / charge), 2.0) * *rpart1 * *rpart2
255 * (*distcc - *rpart1) * (*distcc - *rpart2)
256 / (*distcc * (*distcc * (*rpart1 + *rpart2) - pow (*rpart1, 2.0) - pow (*rpart2, 2.0)))
257 * (omega1 * log (1.0 + gamma) + omega2 * log (1.0 - gamma));
258 }
259
260 /*============================================================================
261 * Public function definitions
262 *============================================================================*/
263
264 /*----------------------------------------------------------------------------*/
265 /*!
266 * \brief Calculation of the adhesion force and adhesion energy
267 *
268 * \param[in] ip particle number
269 * \param[in] tempf thermal scalar value at current time step
270 * \param[out] adhesion_energ particle adhesion energy
271 */
272 /*----------------------------------------------------------------------------*/
273
274 void
cs_lagr_adh(cs_lnum_t ip,cs_real_t tempf,cs_real_t * adhesion_energ)275 cs_lagr_adh(cs_lnum_t ip,
276 cs_real_t tempf,
277 cs_real_t *adhesion_energ)
278 {
279 /* ================================================================ */
280 /* 0. initialization */
281 /* ================================================================ */
282
283 cs_lagr_particle_set_t *p_set = cs_glob_lagr_particle_set;
284 const cs_lagr_attribute_map_t *p_am = p_set->p_am;
285 unsigned char *part = p_set->p_buffer + p_am->extents * ip;
286 cs_lagr_physico_chemical_t *lag_pc = cs_glob_lagr_physico_chemical;
287
288 /* step = step used to calculate the adhesion force following */
289 /* F = U(_d_cut_off+step)-U(_d_cut_off-step)/(2*step) */
290
291 cs_real_t denasp = cs_glob_lagr_reentrained_model->denasp;
292 cs_real_t rayasg = cs_glob_lagr_reentrained_model->rayasg;
293 cs_real_t rayasp = cs_glob_lagr_reentrained_model->rayasp;
294 cs_real_t espasg = cs_glob_lagr_reentrained_model->espasg;
295 cs_real_t modyeq = cs_glob_lagr_reentrained_model->modyeq;
296
297 cs_real_t step = 1e-11;
298 cs_real_t scovap = denasp * cs_math_pi * pow (rayasp, 2);
299 cs_real_t scovag = cs_math_pi * pow (rayasg, 2) / pow (espasg, 2);
300
301 cs_real_t dismin;
302 /* ================================================================ */
303 /* 3. calculation of the adhesion force */
304 /* ================================================================ */
305
306 /* determination of the number of contacts with asperities */
307 /* ======================================================= */
308
309 /* *************************************/
310 /* Number of large-scale asperities */
311 /* *************************************/
312
313 cs_real_t rpart = 0.5 * cs_lagr_particle_get_real(part, p_am, CS_LAGR_DIAMETER);
314
315 cs_real_t nmoyag = (2.0 * rpart + rayasg) / rayasg * scovag;
316
317 if (nmoyag > 600.0) {
318
319 cs_lnum_t tmp = 1;
320 cs_real_t rtmp;
321
322 cs_random_normal(1, &rtmp);
323
324 tmp = (int)nmoyag + sqrt(nmoyag) * rtmp;
325 tmp = CS_MAX(0, tmp);
326
327 cs_lagr_particle_set_lnum(part, p_am, CS_LAGR_N_LARGE_ASPERITIES, tmp);
328
329 }
330 else {
331
332 int ntmp;
333
334 cs_random_poisson(1, nmoyag, &ntmp);
335
336 cs_lagr_particle_set_lnum(part, p_am, CS_LAGR_N_LARGE_ASPERITIES, ntmp);
337
338 }
339
340 int nbasg = 0;
341 if (cs_lagr_particle_get_lnum(part, p_am, CS_LAGR_N_LARGE_ASPERITIES) > 1) {
342
343 nmoyag = 1.0
344 + 2.0 * _d_cut_off * (2.0 * rpart + 2.0 * rayasg + 4.0 * _d_cut_off)
345 / pow (rayasg, 2) * scovag;
346
347 if (nmoyag > 600.0) {
348
349 cs_real_t rtmp;
350
351 cs_random_normal(1, &rtmp);
352
353 nbasg = (int)nmoyag + sqrt (nmoyag) * rtmp;
354 nbasg = CS_MAX(0, nbasg);
355
356 }
357 else {
358
359 int ntmp;
360
361 cs_random_poisson(1, nmoyag, &ntmp);
362
363 nbasg = ntmp;
364
365 }
366
367 nbasg = CS_MAX(1, nbasg);
368
369 }
370 else {
371
372 nbasg = cs_lagr_particle_get_lnum(part, p_am, CS_LAGR_N_LARGE_ASPERITIES);
373
374 }
375
376 /* *******************************/
377 /* Nb of small-scale asperities */
378 /* *******************************/
379
380 /* 1st case: no large-scale asperities*/
381 cs_lnum_t nbasp = 0;
382 if (nbasg == 0) {
383
384 cs_real_t nmoyap = (2.0 * rpart + rayasp) / rayasp * scovap;
385
386 if (nmoyap > 600.0) {
387
388 cs_lnum_t tmp = 1;
389 cs_real_t rtmp;
390
391 cs_random_normal(1, &rtmp);
392
393 tmp = (int)nmoyap + sqrt(nmoyap) * rtmp;
394 tmp = CS_MAX(0, tmp);
395
396 cs_lagr_particle_set_lnum(part, p_am, CS_LAGR_N_SMALL_ASPERITIES, tmp);
397
398 }
399 else {
400
401 int ntmp;
402
403 cs_random_poisson(1, nmoyap, &ntmp);
404
405 cs_lagr_particle_set_lnum(part, p_am, CS_LAGR_N_SMALL_ASPERITIES, ntmp);
406
407 }
408
409 if (cs_lagr_particle_get_lnum(part, p_am, CS_LAGR_N_SMALL_ASPERITIES) > 1) {
410
411 nmoyap = 1
412 + 2.0 * _d_cut_off * (2.0 * rpart + 2.0 * rayasp + 4.0 * _d_cut_off)
413 / pow (rayasp, 2) * scovap;
414
415 if (nmoyap > 600.0) {
416
417 cs_real_t rtmp;
418
419 cs_random_normal(1, &rtmp);
420
421 nbasp = (int)nmoyap + sqrt (nmoyap) * rtmp;
422 nbasp = CS_MAX(0, nbasp);
423
424 cs_lagr_particle_set_lnum(part, p_am, CS_LAGR_N_SMALL_ASPERITIES, nbasp);
425
426 }
427 else {
428
429 int ntmp;
430
431 cs_random_poisson(1, nmoyap, &ntmp);
432 nbasp = ntmp;
433
434 }
435
436 nbasp = CS_MAX(1, nbasp);
437
438 }
439 else {
440
441 nbasp = cs_lagr_particle_get_lnum(part, p_am, CS_LAGR_N_SMALL_ASPERITIES);
442
443 }
444
445 /* Determination of the minimal distance between the particle and the plate */
446
447 dismin = rayasp * CS_MIN (1.0, cs_lagr_particle_get_lnum(part, p_am, CS_LAGR_N_SMALL_ASPERITIES));
448
449 }
450 /* 2nd case: contact with large-scale asperities */
451 else {
452
453 cs_real_t paramh = 0.5 * (2.0 * rpart + rayasp) * rayasp / (rpart + rayasg);
454 cs_real_t nmoyap = paramh * (2 * rayasg - paramh) / pow (rayasp, 2) * scovap;
455
456 if (nmoyap > 600.0) {
457
458 cs_lnum_t tmp = 1;
459 cs_real_t rtmp;
460
461 cs_random_normal(1, &rtmp);
462
463 tmp = (int)nmoyap + sqrt(nmoyap) * rtmp;
464 tmp = CS_MAX(0, tmp);
465
466 cs_lagr_particle_set_lnum(part, p_am, CS_LAGR_N_SMALL_ASPERITIES, tmp);
467
468 }
469 else {
470
471 int ntmp;
472
473 cs_random_poisson(1, nmoyap, &ntmp);
474
475 cs_lagr_particle_set_lnum(part, p_am, CS_LAGR_N_SMALL_ASPERITIES, ntmp);
476
477 }
478
479 if (cs_lagr_particle_get_lnum(part, p_am, CS_LAGR_N_SMALL_ASPERITIES) > 1) {
480
481 paramh = 0.5 * (2.0 * rpart + 2 * rayasp + 4.0 * _d_cut_off)
482 * 2.0 * _d_cut_off / (rpart + rayasg + rayasp + _d_cut_off);
483 nmoyap = 1 + paramh * (2 * rayasg - paramh) / pow (rayasp, 2) * scovap;
484
485 if (nmoyap > 600.0) {
486
487 cs_real_t rtmp;
488
489 cs_random_normal(1, &rtmp);
490
491 nbasp = (int)nmoyap + sqrt (nmoyap) * rtmp;
492 nbasp = CS_MAX(0, nbasp);
493
494 }
495 else {
496
497 int ntmp;
498
499 cs_random_poisson(1, nmoyap, &ntmp);
500
501 nbasp = ntmp;
502
503 }
504
505 nbasp = CS_MAX(1, nbasp);
506
507 }
508 else
509 nbasp = cs_lagr_particle_get_lnum(part, p_am, CS_LAGR_N_SMALL_ASPERITIES);
510
511 /* Mutliple contacts with large scale asperities? */
512 nbasp = nbasp * nbasg;
513 cs_lagr_particle_set_lnum(part, p_am, CS_LAGR_N_SMALL_ASPERITIES,
514 cs_lagr_particle_get_lnum(part, p_am, CS_LAGR_N_SMALL_ASPERITIES)
515 *cs_lagr_particle_get_lnum(part, p_am, CS_LAGR_N_LARGE_ASPERITIES));
516
517 /* Determination of the minimal distance between the particle and the plate */
518 dismin = rayasp * CS_MIN(1.0, nbasp * 1.0) + rayasg * CS_MIN (1.0, nbasg * 1.0);
519
520 }
521
522 /* Sum of {particle-plane} and {particle-asperity} energies */
523
524 cs_real_t udlvor[2];
525 /* Interaction between the sphere and the plate */
526 for (cs_lnum_t np = 0; np < 2; np++) {
527
528 udlvor[np] = 0.0;
529 cs_real_t distp = dismin + _d_cut_off + step * (3 - 2 * (np+1));
530
531 cs_real_t uvdwsp, uedlsp;
532 _vdwsp(lag_pc->cstham, lag_pc->lambda_vdw, &distp, &rpart, &uvdwsp);
533 _edlsp(&distp, &rpart, lag_pc->phi_p, lag_pc->phi_s, tempf, &uedlsp);
534
535 udlvor[np] = (uvdwsp + uedlsp) * (1 - scovag - scovap);
536
537 }
538
539 cs_real_t fadhes = (udlvor[1] - udlvor[0]) / (2.0 * step);
540
541 *adhesion_energ = udlvor[0];
542
543 /* Interaction between the sphere and small-scale asperities */
544 for (cs_lnum_t np = 0; np < 2; np++) {
545
546 udlvor[np] = 0.0;
547
548 cs_real_t distcc = _d_cut_off + step * (3 - 2 * (np + 1)) + rpart + rayasp;
549
550 cs_real_t uvdwss, uedlss;
551
552 _vdwss(lag_pc->cstham, lag_pc->lambda_vdw, &distcc, &rpart, &rayasp, &uvdwss);
553 _edlss(&distcc, &rpart, &rayasp, lag_pc->phi_p, lag_pc->phi_s, tempf, &uedlss);
554
555 udlvor[np] = uvdwss + uedlss;
556
557 }
558
559 fadhes = fadhes + (udlvor[1] - udlvor[0]) / (2.0 * step) * nbasp;
560
561 *adhesion_energ = *adhesion_energ + udlvor[0] * nbasp;
562
563 /* Interaction between the sphere and large-scale asperities */
564 for (cs_lnum_t np = 0; np < 2; np++) {
565
566 udlvor[np] = 0.0;
567 cs_real_t distcc;
568
569 if (nbasp == 0)
570 distcc = _d_cut_off + step * (3 - 2 * (np+1)) + rpart + rayasg;
571
572 else if (nbasp > 0)
573 distcc = _d_cut_off + rayasp + step * (3 - 2 * (np+1)) + rpart + rayasg;
574
575 cs_real_t uvdwss, uedlss;
576
577 _vdwss(lag_pc->cstham, lag_pc->lambda_vdw, &distcc, &rpart, &rayasg, &uvdwss);
578 _edlss(&distcc, &rpart, &rayasg, lag_pc->phi_p, lag_pc->phi_s, tempf, &uedlss);
579
580 udlvor[np] = uvdwss + uedlss;
581
582 }
583
584 fadhes = fadhes + (udlvor[1] - udlvor[0]) / (2.0 * step) * nbasg;
585
586 *adhesion_energ = *adhesion_energ + udlvor[0] * nbasg;
587
588 /* The force is negative when it is attractive */
589
590 if (fadhes >= 0.0)
591 cs_lagr_particle_set_real(part, p_am, CS_LAGR_ADHESION_FORCE, 0.0);
592
593 else
594 cs_lagr_particle_set_real(part, p_am, CS_LAGR_ADHESION_FORCE, -fadhes);
595
596 /* The interaction should be negative to prevent reentrainment (attraction) */
597 if (*adhesion_energ >= 0.0)
598 *adhesion_energ = 0.0;
599
600 else
601 *adhesion_energ = CS_ABS(*adhesion_energ);
602
603 /* Calculation of adhesion torques exerted on the particle */
604 cs_real_t rtmp;
605 cs_random_uniform(1, &rtmp);
606
607 cs_real_t dismom = rtmp;
608 if (nbasp > 0)
609 dismom = (pow(rtmp, 1.0 / cs_lagr_particle_get_lnum(part, p_am, CS_LAGR_N_SMALL_ASPERITIES)) * 2.0 - 1.0)
610 * sqrt ((2.0 * rpart + rayasp) * rayasp);
611
612 else {
613
614 /* In the sphere-plate case, we use the deformation given by the DMT theory,
615 * which is close to our approach */
616
617 cs_real_t omsurf = lag_pc->cstham / (24.0 * cs_math_pi * pow (_d_cut_off, 2));
618 dismom = pow ((4.0 * cs_math_pi * omsurf * (pow (rpart, 2.0)) / modyeq), (1.0 / 3.0));
619 /* dismom = (12.0d0 * cs_math_pi * omsurf * (rpart**2)/modyeq)**(1.0d0/3.0d0) */
620
621 }
622
623 dismom *= cs_lagr_particle_get_real(part, p_am, CS_LAGR_ADHESION_FORCE);
624
625 cs_lagr_particle_set_real(part, p_am, CS_LAGR_ADHESION_TORQUE, dismom);
626
627 }
628
629 /*============================================================================
630 * Public function definitions
631 *============================================================================*/
632
633 /*----------------------------------------------------------------------------*/
634 /*!
635 * \brief Calculation of the adhesion force and adhesion energy
636 * between two particles
637 *
638 * \param[in] dpart particle diameter
639 * \param[in] tempf thermal scalar value at current time step
640 * \param[out] adhesion_energ particle adhesion energy
641 * \param[out] adhesion_force particle adhesion force
642 */
643 /*----------------------------------------------------------------------------*/
644
645 void
cs_lagr_adh_pp(cs_real_t dpart,cs_real_t tempf,cs_real_t * adhesion_energ,cs_real_t * adhesion_force)646 cs_lagr_adh_pp(cs_real_t dpart,
647 cs_real_t tempf,
648 cs_real_t *adhesion_energ,
649 cs_real_t *adhesion_force)
650 {
651 /* ================================================================ */
652 /* 0. initialization */
653 /* ================================================================ */
654
655 cs_lagr_physico_chemical_t *lag_pc = cs_glob_lagr_physico_chemical;
656 cs_lagr_clogging_model_t *lag_cm = cs_glob_lagr_clogging_model;
657
658 cs_real_t step = 1e-11;
659
660 /* ================================================================ */
661 /* 1. calculation of the adhesion force */
662 /* ================================================================ */
663
664 cs_real_t rpart = 0.5 * dpart;
665
666 /* Sum of {particle-plane} and {particle-asperity} energies */
667
668 cs_real_t udlvor[2];
669
670 /* Interaction between the spheres */
671
672 for (cs_lnum_t np = 0; np < 2; np++) {
673
674 udlvor[np] = 0.0;
675
676 cs_real_t distcc = _d_cut_off + step * (3 - 2 * (np + 1)) + 2.0 * rpart;
677
678 cs_real_t uvdwss, uedlss;
679
680 _vdwss(lag_cm->csthpp, lag_pc->lambda_vdw, &distcc, &rpart, &rpart, &uvdwss);
681 _edlss(&distcc, &rpart, &rpart, lag_pc->phi_p, lag_pc->phi_p, tempf, &uedlss);
682
683 udlvor[np] = uvdwss + uedlss;
684
685 }
686
687 *adhesion_force = CS_MAX( -(udlvor[1] - udlvor[0]) / (2.0 * step), 0.0);
688
689 *adhesion_energ = CS_MAX(-udlvor[0], 0.0);
690
691 }
692
693 END_C_DECLS
694