1 /*============================================================================
2 * Calculation of DLVO 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 the DLVO forces
29 *============================================================================*/
30
31 #include "cs_defs.h"
32
33 /*----------------------------------------------------------------------------
34 * Standard C library headers
35 *----------------------------------------------------------------------------*/
36
37 #include <stdio.h>
38 #include <stdlib.h>
39 #include <math.h>
40 #include <float.h>
41 #include <assert.h>
42
43 /*----------------------------------------------------------------------------
44 * Local headers
45 *----------------------------------------------------------------------------*/
46
47 #include "bft_printf.h"
48 #include "bft_error.h"
49 #include "bft_mem.h"
50
51 #include "fvm_periodicity.h"
52
53 #include "cs_base.h"
54 #include "cs_interface.h"
55 #include "cs_mesh.h"
56 #include "cs_mesh_quantities.h"
57 #include "cs_parall.h"
58 #include "cs_search.h"
59 #include "cs_halo.h"
60
61 #include "cs_lagr.h"
62 #include "cs_lagr_roughness.h"
63
64 /*----------------------------------------------------------------------------
65 * Header for the current file
66 *----------------------------------------------------------------------------*/
67
68 #include "cs_lagr_dlvo.h"
69
70 /*----------------------------------------------------------------------------*/
71
72 BEGIN_C_DECLS
73
74 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
75
76 /*============================================================================
77 * Local macro declarations
78 *============================================================================*/
79
80 #define PG_CST 8.314 /* Ideal gas constant */
81
82 /*============================================================================
83 * Local structure declarations
84 *============================================================================*/
85
86 static cs_lagr_dlvo_param_t cs_lagr_dlvo_param;
87
88 /*============================================================================
89 * Static global variables
90 *============================================================================*/
91
92 static const double _pi = 3.14159265358979323846;
93
94 /* Cut-off distance for adhesion forces (assumed to be the Born distance) */
95 static const cs_real_t _d_cut_off = 1.65e-10;
96
97 /* Boltzmann constant */
98 static const double _k_boltzmann = 1.38e-23;
99
100 /* Free space permittivity */
101 static const cs_real_t _free_space_permit = 8.854e-12;
102
103 /* Faraday constant */
104 static const cs_real_t _faraday_cst = 9.648e4;
105
106 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
107
108 /*============================================================================
109 * Public function prototypes for Fortran API
110 *============================================================================*/
111
112 /*----------------------------------------------------------------------------
113 * DLVO initialization:
114 * - Retrieve various parameters for storing in global structure
115 * - Compute and store the Debye screening length
116 *----------------------------------------------------------------------------*/
117
118 void
cs_lagr_dlvo_init(const cs_real_t water_permit,const cs_real_t ionic_strength,const cs_real_t temperature[],const cs_real_t valen,const cs_real_t phi_p,const cs_real_t phi_s,const cs_real_t cstham,const cs_real_t csthpp,const cs_real_t lambda_vdw)119 cs_lagr_dlvo_init(const cs_real_t water_permit,
120 const cs_real_t ionic_strength,
121 const cs_real_t temperature[],
122 const cs_real_t valen,
123 const cs_real_t phi_p,
124 const cs_real_t phi_s,
125 const cs_real_t cstham,
126 const cs_real_t csthpp,
127 const cs_real_t lambda_vdw)
128 {
129 cs_lnum_t iel;
130
131 const cs_mesh_t *mesh = cs_glob_mesh;
132
133 /* Retrieve physical parameters related to clogging modeling */
134 /* and fill the global structure cs_lagr_clog_param */
135
136 cs_lagr_dlvo_param.water_permit = water_permit;
137 cs_lagr_dlvo_param.ionic_strength = ionic_strength;
138 cs_lagr_dlvo_param.valen = valen;
139 cs_lagr_dlvo_param.phi_p = phi_p;
140 cs_lagr_dlvo_param.phi_s = phi_s;
141 cs_lagr_dlvo_param.cstham = cstham;
142 cs_lagr_dlvo_param.cstham = csthpp;
143 cs_lagr_dlvo_param.lambda_vdw = lambda_vdw;
144
145 /* Allocate memory for the temperature and Debye length arrays */
146
147 if (cs_lagr_dlvo_param.temperature == NULL)
148 BFT_MALLOC(cs_lagr_dlvo_param.temperature, mesh->n_cells, cs_real_t);
149
150 if (cs_lagr_dlvo_param.debye_length == NULL)
151 BFT_MALLOC(cs_lagr_dlvo_param.debye_length, mesh->n_cells, cs_real_t);
152
153 /* Store the temperature */
154
155 for (iel = 0; iel < mesh->n_cells; iel++)
156 cs_lagr_dlvo_param.temperature[iel] = temperature[iel];
157
158 /* Computation and storage of the Debye length */
159
160 for (iel = 0; iel < mesh->n_cells ; iel++)
161
162 cs_lagr_dlvo_param.debye_length[iel]
163 = pow(2e3 * pow(_faraday_cst,2)
164 * cs_lagr_dlvo_param.ionic_strength /
165 (cs_lagr_dlvo_param.water_permit
166 * _free_space_permit * PG_CST
167 * cs_lagr_dlvo_param.temperature[iel]), -0.5);
168
169 #if 0 && defined(DEBUG) && !defined(NDEBUG)
170 bft_printf(" epseau = %g\n", cs_lagr_dlvo_param.water_permit);
171 bft_printf(" fion = %g\n", cs_lagr_dlvo_param.ionic_strength);
172 bft_printf(" temp[0] = %g\n", cs_lagr_dlvo_param.temperature[0]);
173 bft_printf(" valen = %g\n", cs_lagr_dlvo_param.valen);
174 bft_printf(" debye[0] = %g\n", cs_lagr_dlvo_param.debye_length[0]);
175 bft_printf(" phi_p = %g\n", cs_lagr_dlvo_param.phi_p);
176 bft_printf(" phi_s = %g\n", cs_lagr_dlvo_param.phi_s);
177 #endif
178
179 }
180
181 /*============================================================================
182 * Public function prototypes
183 *============================================================================*/
184
185 /*----------------------------------------------------------------------------
186 * Deallocate the arrays storing temperature and Debye length.
187 *----------------------------------------------------------------------------*/
188
189 void
cs_lagr_dlvo_finalize()190 cs_lagr_dlvo_finalize()
191 {
192 BFT_FREE(cs_lagr_dlvo_param.temperature);
193 BFT_FREE(cs_lagr_dlvo_param.debye_length);
194 }
195
196 /*----------------------------------------------------------------------------
197 * Compute the energy barrier for a smooth wall.
198 *----------------------------------------------------------------------------*/
199
200 void
cs_lagr_barrier(const void * particle,const cs_lagr_attribute_map_t * attr_map,cs_lnum_t iel,cs_real_t * energy_barrier)201 cs_lagr_barrier(const void *particle,
202 const cs_lagr_attribute_map_t *attr_map,
203 cs_lnum_t iel,
204 cs_real_t *energy_barrier)
205 {
206 cs_lnum_t i;
207 cs_real_t rpart = cs_lagr_particle_get_real(particle, attr_map,
208 CS_LAGR_DIAMETER) * 0.5;
209
210 *energy_barrier = 0.;
211
212 cs_real_t barr = 0.;
213 cs_real_t distp = 0.;
214
215 /* Computation of the energy barrier */
216
217 for (i = 0; i < 1001; i++) {
218
219 cs_real_t step = cs_lagr_dlvo_param.debye_length[iel]/30.0;
220
221 /* Interaction between the sphere and the plate */
222
223 distp = _d_cut_off + i * step;
224
225 cs_real_t var1
226 = cs_lagr_van_der_waals_sphere_plane(distp,
227 rpart,
228 cs_lagr_dlvo_param.lambda_vdw,
229 cs_lagr_dlvo_param.cstham);
230
231 cs_real_t var2
232 = cs_lagr_edl_sphere_plane(distp,
233 rpart,
234 cs_lagr_dlvo_param.valen,
235 cs_lagr_dlvo_param.phi_p,
236 cs_lagr_dlvo_param.phi_s,
237 cs_lagr_dlvo_param.temperature[iel],
238 cs_lagr_dlvo_param.debye_length[iel],
239 cs_lagr_dlvo_param.water_permit);
240
241 barr = (var1 + var2);
242
243 if (barr > *energy_barrier)
244 *energy_barrier = barr;
245 if ( *energy_barrier < 0)
246 *energy_barrier = 0;
247 }
248
249 *energy_barrier = *energy_barrier / rpart;
250 }
251
252 /*----------------------------------------------------------------------------
253 * Compute the energy barrier for two particles.
254 *----------------------------------------------------------------------------*/
255
256 void
cs_lagr_barrier_pp(cs_real_t dpart,cs_lnum_t iel,cs_real_t * energy_barrier)257 cs_lagr_barrier_pp(cs_real_t dpart,
258 cs_lnum_t iel,
259 cs_real_t *energy_barrier)
260 {
261 cs_real_t rpart = dpart * 0.5;
262
263 *energy_barrier = 0.;
264
265 /* Computation of the energy barrier */
266
267 for (int i = 0; i < 1001; i++) {
268
269 cs_real_t step = cs_lagr_dlvo_param.debye_length[iel] / 30.0;
270
271 /* Interaction between two spheres */
272
273 cs_real_t distcc = _d_cut_off + i * step + 2.0 * rpart;
274
275 cs_real_t var1
276 = cs_lagr_van_der_waals_sphere_sphere(distcc,
277 rpart,
278 rpart,
279 cs_lagr_dlvo_param.lambda_vdw,
280 cs_lagr_dlvo_param.csthpp);
281
282 cs_real_t var2
283 = cs_lagr_edl_sphere_sphere(distcc,
284 rpart,
285 rpart,
286 cs_lagr_dlvo_param.valen,
287 cs_lagr_dlvo_param.phi_p,
288 cs_lagr_dlvo_param.phi_p,
289 cs_lagr_dlvo_param.temperature[iel],
290 cs_lagr_dlvo_param.debye_length[iel],
291 cs_lagr_dlvo_param.water_permit);
292
293 cs_real_t barr = (var1 + var2);
294
295 if (barr > *energy_barrier)
296 *energy_barrier = barr;
297 if (*energy_barrier < 0)
298 *energy_barrier = 0;
299 }
300
301 *energy_barrier = *energy_barrier / rpart;
302 }
303
304 /*----------------------------------------------------------------------------
305 * Van der Waals interaction between a sphere and a plane
306 * using formulas from Czarnecki (large distances)
307 * and Gregory (small distances)
308 *----------------------------------------------------------------------------*/
309
310 cs_real_t
cs_lagr_van_der_waals_sphere_plane(cs_real_t distp,cs_real_t rpart,cs_real_t lambda_vdw,cs_real_t cstham)311 cs_lagr_van_der_waals_sphere_plane (cs_real_t distp,
312 cs_real_t rpart,
313 cs_real_t lambda_vdw,
314 cs_real_t cstham)
315 {
316 cs_real_t var;
317
318 if (distp < (lambda_vdw / 2 / _pi)) {
319 var = -cstham * rpart / (6 * distp)
320 * (1 / (1 + 14 * distp / lambda_vdw + 5 * _pi/4.9
321 * pow(distp,3) / lambda_vdw / pow(rpart,2)));
322 }
323 else {
324 var = cstham
325 * ((2.45 * lambda_vdw ) /(60. * _pi)
326 * ( (distp - rpart) / pow(distp,2)
327 - (distp + 3. * rpart) / pow(distp + 2. * rpart,2))
328 - 2.17 / 720. / pow(_pi,2) * pow(lambda_vdw,2) * ((distp - 2. * rpart)
329 / pow(distp,3) - (distp + 4. * rpart) / pow(distp + 2. * rpart,3))
330 + 0.59 / 5040. / pow(_pi,3) * pow(lambda_vdw,3)
331 * ((distp - 3. * rpart) / pow(distp,4) - (distp + 5. * rpart)
332 / pow(distp + 2. * rpart,4)));
333 }
334
335 return var;
336 }
337
338 /*----------------------------------------------------------------------------
339 * Calculation of the Van der Waals interaction between two spheres
340 * following the formula from Gregory (1981a)
341 *----------------------------------------------------------------------------*/
342
343 cs_real_t
cs_lagr_van_der_waals_sphere_sphere(cs_real_t distcc,cs_real_t rpart1,cs_real_t rpart2,cs_real_t lambda_vdw,cs_real_t cstham)344 cs_lagr_van_der_waals_sphere_sphere(cs_real_t distcc,
345 cs_real_t rpart1,
346 cs_real_t rpart2,
347 cs_real_t lambda_vdw,
348 cs_real_t cstham)
349 {
350 cs_real_t var = - cstham * rpart1 * rpart2 / (6 * (distcc - rpart1 - rpart2)
351 * (rpart1 + rpart2)) * (1 - 5.32 * (distcc - rpart1 - rpart2)
352 / lambda_vdw * log(1 + lambda_vdw / (distcc - rpart1 - rpart2) / 5.32));
353
354 return var;
355 }
356
357 /*----------------------------------------------------------------------------
358 * Electric Double Layer (EDL) interaction between a sphere and a plane
359 * using the formula from Bell & al (1970)
360 * based on the McCartney & Levine method
361 *----------------------------------------------------------------------------*/
362
363 cs_real_t
cs_lagr_edl_sphere_plane(cs_real_t distp,cs_real_t rpart,cs_real_t valen,cs_real_t phi1,cs_real_t phi2,cs_real_t temp,cs_real_t debye_length,cs_real_t water_permit)364 cs_lagr_edl_sphere_plane(cs_real_t distp,
365 cs_real_t rpart,
366 cs_real_t valen,
367 cs_real_t phi1,
368 cs_real_t phi2,
369 cs_real_t temp,
370 cs_real_t debye_length,
371 cs_real_t water_permit)
372 {
373
374 cs_real_t charge = 1.6e-19;
375 /* Reduced zeta potential */
376 cs_real_t lphi1 = valen * charge * phi1 / _k_boltzmann / temp;
377 cs_real_t lphi2 = valen * charge * phi2 / _k_boltzmann / temp;
378
379 cs_real_t tau = rpart / debye_length;
380
381 /* Extended reduced zeta potential */
382 /* (following the work from Ohshima et al, 1982, JCIS, 90, 17-26) */
383
384 lphi1 = 8. * tanh(lphi1 / 4.) /
385 ( 1. + pow(1. - (2. * tau + 1.) / (pow(tau + 1,2))
386 * pow(tanh(lphi1 / 4.),2),0.5));
387
388 lphi2 = 4. * tanh(lphi2 / 4.) ;
389
390 cs_real_t alpha = sqrt((distp + rpart) / rpart)
391 + sqrt(rpart / (distp + rpart));
392 cs_real_t omega1 = pow(lphi1,2) + pow(lphi2,2) + alpha * lphi1 * lphi2;
393 cs_real_t omega2 = pow(lphi1,2) + pow(lphi2,2) - alpha * lphi1 * lphi2;
394 cs_real_t gamma = sqrt(rpart / (distp + rpart)) * exp(-1./debye_length * distp);
395
396 cs_real_t var = 2 * _pi * _free_space_permit * water_permit
397 * pow((_k_boltzmann * temp / (1. * valen) / charge),2)
398 * rpart * (distp + rpart) / (distp + 2 * rpart)
399 * (omega1 * log(1 + gamma) + omega2 * log(1 - gamma));
400
401 return var;
402 }
403
404 /*----------------------------------------------------------------------------
405 * Calculation of the EDL interaction between two spheres
406 * using the formula from Bell & al (1970)
407 * based on the McCartney & Levine method
408 *----------------------------------------------------------------------------*/
409
410 cs_real_t
cs_lagr_edl_sphere_sphere(cs_real_t distcc,cs_real_t rpart1,cs_real_t rpart2,cs_real_t valen,cs_real_t phi1,cs_real_t phi2,cs_real_t temp,cs_real_t debye_length,cs_real_t water_permit)411 cs_lagr_edl_sphere_sphere(cs_real_t distcc,
412 cs_real_t rpart1,
413 cs_real_t rpart2,
414 cs_real_t valen,
415 cs_real_t phi1,
416 cs_real_t phi2,
417 cs_real_t temp,
418 cs_real_t debye_length,
419 cs_real_t water_permit)
420 {
421 cs_real_t charge = 1.6e-19;
422
423 /* Reduced zeta potential */
424 cs_real_t lphi1 = valen * charge * phi1 / _k_boltzmann / temp;
425 cs_real_t lphi2 = valen * charge * phi2 / _k_boltzmann / temp;
426
427
428 /* Extended reduced zeta potential */
429 /* (following the work from Ohshima et al, 1982, JCIS, 90, 17-26) */
430
431 cs_real_t tau1 = rpart1 / debye_length;
432 lphi1 = 8. * tanh(lphi1 / 4.) /
433 ( 1. + pow(1. - (2. * tau1 + 1.) / (pow(tau1 + 1,2))
434 * pow(tanh(lphi1 / 4.),2),0.5));
435
436 cs_real_t tau2 = rpart2 / debye_length;
437 lphi2 = 8. * tanh(lphi2 / 4.) /
438 ( 1. + pow(1. - (2. * tau2 + 1.) / (pow(tau2 + 1,2))
439 * pow(tanh(lphi2 / 4.),2),0.5));
440
441 cs_real_t alpha = sqrt(rpart2 * (distcc - rpart2)
442 / (rpart1 * (distcc - rpart1)))
443 + sqrt(rpart1 * (distcc - rpart1)
444 / (rpart2 * (distcc - rpart2)));
445
446 cs_real_t omega1 = pow(lphi1,2) + pow(lphi2,2) + alpha * lphi1 * lphi2;
447
448 cs_real_t omega2 = pow(lphi1,2) + pow(lphi2,2) - alpha * lphi1 * lphi2;
449
450 cs_real_t gamma = sqrt(rpart1 * rpart2 / (distcc-rpart1) / (distcc-rpart2))
451 *exp(1. / debye_length * (rpart1 + rpart2 - distcc));
452
453 cs_real_t var = 2 * _pi * _free_space_permit * water_permit
454 * pow((_k_boltzmann * temp / charge),2)
455 * rpart1 * rpart2 * (distcc - rpart1) * (distcc - rpart2)
456 / (distcc * ( distcc * (rpart1 + rpart2)
457 - pow(rpart1,2) - pow(rpart2,2)))
458 * (omega1 * log(1 + gamma) + omega2 * log(1 - gamma));
459
460 return var;
461 }
462
463 /*----------------------------------------------------------------------------*/
464
465 END_C_DECLS
466