1
2 /*============================================================================
3 * Methods for head losses due to particle deposit
4 *============================================================================*/
5
6 /*
7 This file is part of Code_Saturne, a general-purpose CFD tool.
8
9 Copyright (C) 1998-2021 EDF S.A.
10
11 This program is free software; you can redistribute it and/or modify it under
12 the terms of the GNU General Public License as published by the Free Software
13 Foundation; either version 2 of the License, or (at your option) any later
14 version.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
19 details.
20
21 You should have received a copy of the GNU General Public License along with
22 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
23 Street, Fifth Floor, Boston, MA 02110-1301, USA.
24 */
25
26 /*----------------------------------------------------------------------------*/
27
28 /*============================================================================
29 * Functions dealing with head losses due to particle deposit
30 *============================================================================*/
31
32 #include "cs_defs.h"
33
34 /*----------------------------------------------------------------------------
35 * Standard C library headers
36 *----------------------------------------------------------------------------*/
37
38 #include <limits.h>
39 #include <stdio.h>
40 #include <stddef.h>
41 #include <stdlib.h>
42 #include <string.h>
43 #include <math.h>
44 #include <ctype.h>
45 #include <float.h>
46 #include <assert.h>
47
48 /*----------------------------------------------------------------------------
49 * Local headers
50 *----------------------------------------------------------------------------*/
51
52 #include "bft_mem.h"
53 #include "bft_printf.h"
54
55 #include "cs_field.h"
56 #include "cs_field_operator.h"
57 #include "cs_gradient.h"
58 #include "cs_halo.h"
59 #include "cs_halo_perio.h"
60 #include "cs_math.h"
61 #include "cs_mesh.h"
62 #include "cs_mesh_quantities.h"
63 #include "cs_parameters.h"
64 #include "cs_parameters_check.h"
65 #include "cs_parall.h"
66 #include "cs_prototypes.h"
67
68 #include "cs_lagr.h"
69 #include "cs_lagr_tracking.h"
70 #include "cs_lagr_stat.h"
71
72 /*----------------------------------------------------------------------------
73 * Header for the current file
74 *----------------------------------------------------------------------------*/
75
76 #include "cs_lagr_head_losses.h"
77
78 /*----------------------------------------------------------------------------*/
79
80 BEGIN_C_DECLS
81
82 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
83
84 /*============================================================================
85 * Private function definitions
86 *============================================================================*/
87
88 /*----------------------------------------------------------------------------
89 * Compute the porosity in wall-normal cells based on the mean deposit height
90 * (which is only known at the boundary faces)
91 *----------------------------------------------------------------------------*/
92
93 static void
_porcel(cs_real_t mdiam[],cs_real_t porosi[],const cs_lnum_t bc_type[])94 _porcel(cs_real_t mdiam[],
95 cs_real_t porosi[],
96 const cs_lnum_t bc_type[])
97 {
98 /* Initialization */
99
100 const cs_mesh_t *m = cs_glob_mesh;
101
102 cs_lnum_t nfabor = m->n_b_faces;
103 cs_lnum_t ncel = m->n_cells;
104 cs_lnum_t ncelet = m->n_cells_with_ghosts;
105
106 cs_real_t *surfbn = cs_glob_mesh_quantities->b_face_surf;
107 cs_real_t *volume = cs_glob_mesh_quantities->cell_vol;
108
109 const cs_real_3_t *restrict i_face_normal
110 = (const cs_real_3_t *restrict)cs_glob_mesh_quantities->i_face_normal;
111 cs_lnum_t *ifabor = cs_glob_mesh->b_face_cells;
112
113 cs_lagr_boundary_interactions_t
114 *lag_bdy_i = cs_glob_lagr_boundary_interactions;
115
116 /*Compute distance to deposition wall, and check deposition exists on
117 at least one face.*/
118
119 cs_lnum_t indic = 0;
120
121 for (cs_lnum_t ifac = 0; ifac < nfabor; ifac++) {
122 if ( (bc_type[ifac] == CS_SMOOTHWALL || bc_type[ifac] == CS_ROUGHWALL)
123 && bound_stat[ifac + lag_bdy_i->ihdepm * nfabor] > 0.0)
124 indic = 1;
125 }
126
127 cs_parall_counter_max(&indic, 1);
128
129 /* mean particle diameter and porosity value due to
130 the deposit in each cell*/
131 for (cs_lnum_t iel = 0; iel < ncelet; iel++) {
132 porosi[iel] = 1.0;
133 mdiam[iel] = 0.0;
134 }
135
136 /* Nothing more to do if no deposition at this stage */
137 if (indic == 0) {
138 return;
139 }
140
141 /* Allocate temporary arrays for the distance resolution */
142 cs_real_3_t *q;
143 cs_real_t *masflu, *depvol;
144
145 cs_field_t *f_wall_dist = cs_field_by_name("wall_distance");
146 BFT_MALLOC(q, ncelet, cs_real_3_t);
147 BFT_MALLOC(masflu, ncelet, cs_real_t);
148 BFT_MALLOC(depvol, ncelet, cs_real_t);
149
150 /* Compute n = Grad(DISTPW)/|Grad(DISTPW)|
151 * ======================================== */
152
153 /* Compute the gradient of the distance to the wall */
154
155 cs_field_gradient_scalar(f_wall_dist,
156 false, /* use_previous_t */
157 1, /* inc */
158 true, /* recompute_cocg */
159 q);
160
161 /* Normalisation (caution, gradient can be zero sometimes) */
162
163 for (cs_lnum_t iel = 0; iel < ncel; iel++) {
164
165 cs_real_t xnorme = CS_MAX(cs_math_3_norm(q[iel]),
166 cs_math_epzero);
167
168 for (cs_lnum_t isou = 0; isou < 3; isou++)
169 q[iel][isou] /= xnorme;
170
171 }
172
173 /* Paralellism and periodicity */
174
175 cs_halo_sync_var_strided(m->halo, CS_HALO_STANDARD, (cs_real_t *)q, 3);
176
177 if (m->n_init_perio > 0)
178 cs_halo_perio_sync_var_vect(m->halo, CS_HALO_STANDARD,
179 (cs_real_t *)q, 3);
180
181 /* Compute porosity
182 * ================= */
183
184 cs_real_t epsi = 1e-06;
185
186 /* volume of the deposited particles
187 (calculated from mean deposit height)*/
188 for (cs_lnum_t iel = 0; iel < ncelet; iel++)
189 depvol[iel] = 0.0;
190
191 indic = 0;
192 for (cs_lnum_t ifac = 0; ifac < nfabor; ifac++) {
193
194 cs_lnum_t iel = ifabor[ifac];
195 depvol[iel] += bound_stat[lag_bdy_i->ihdepm * nfabor + ifac] * surfbn[ifac];
196 mdiam[iel] += bound_stat[lag_bdy_i->ihdiam * nfabor + ifac];
197
198 }
199
200 const cs_real_t fmult = (1.0 - cs_glob_lagr_clogging_model->mporos);
201
202 for (cs_lnum_t ifac = 0; ifac < nfabor; ifac++) {
203
204 cs_lnum_t iel = ifabor[ifac];
205 porosi[iel] = (volume[iel] - fmult * depvol[iel]) / volume[iel];
206
207 if (porosi[iel] < cs_glob_lagr_clogging_model->mporos)
208 indic = 1;
209
210 }
211
212 /* Paralellism and periodicity */
213
214 cs_mesh_sync_var_scal(porosi);
215 cs_parall_counter_max(&indic, 1);
216
217 cs_lnum_t nn = 0;
218 while (indic > 0) {
219
220 for (cs_lnum_t iel = 0; iel < ncelet; iel++)
221 masflu[iel] = 0.0;
222
223 for (cs_lnum_t ifac = 0; ifac < cs_glob_mesh->n_i_faces; ifac++) {
224
225 cs_lnum_t iel1 = cs_glob_mesh->i_face_cells[ifac][0];
226 cs_lnum_t iel2 = cs_glob_mesh->i_face_cells[ifac][1];
227
228 cs_real_t prod1 = cs_math_3_dot_product(q[iel1], i_face_normal[ifac]);
229 cs_real_t prod2 =-cs_math_3_dot_product(q[iel2], i_face_normal[ifac]);
230
231 if (porosi[iel1] < cs_glob_lagr_clogging_model->mporos && prod1 > epsi) {
232 masflu[iel1] -= (porosi[iel1] - cs_glob_lagr_clogging_model->mporos)
233 * volume[iel1];
234 masflu[iel2] += (porosi[iel1] - cs_glob_lagr_clogging_model->mporos)
235 * volume[iel1];
236 mdiam[iel2] = mdiam[iel1];
237 }
238
239 if (porosi[iel2] < cs_glob_lagr_clogging_model->mporos && prod2 > epsi) {
240 masflu[iel2] -= (porosi[iel2] - cs_glob_lagr_clogging_model->mporos)
241 * volume[iel2];
242 masflu[iel1] += (porosi[iel2] - cs_glob_lagr_clogging_model->mporos)
243 * volume[iel2];
244 mdiam[iel1] = mdiam[iel2];
245 }
246
247 }
248 indic = 0;
249
250 for (cs_lnum_t cell_id = 0; cell_id < ncel; cell_id++)
251 porosi[cell_id] = porosi[cell_id] + masflu[cell_id] / volume[cell_id];
252
253
254 /* Paralellism and periodicity */
255 if (cs_glob_rank_id >= 0 || cs_glob_mesh->n_init_perio > 0) {
256 cs_mesh_sync_var_scal(porosi);
257 cs_mesh_sync_var_scal(mdiam);
258 }
259
260 for (cs_lnum_t cell_id = 0; cell_id < ncel; cell_id++) {
261
262 if (porosi[cell_id] < cs_glob_lagr_clogging_model->mporos)
263 indic = 1;
264
265 }
266
267 cs_parall_counter_max(&indic, 1);
268
269 nn++;
270
271 if (nn >= 100)
272 bft_error(__FILE__, __LINE__, 0,
273 "==========================================\n"
274 " Error nn > 100\n"
275 "\n"
276 " Stop inside the porcel subroutine\n"
277 "==========================================\n");
278
279 }
280
281 /* Free memory */
282
283 BFT_FREE(masflu);
284 BFT_FREE(depvol);
285 BFT_FREE(q);
286 }
287
288 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
289
290 /*============================================================================
291 * Public function definitions
292 *============================================================================*/
293
294 /*----------------------------------------------------------------------------*/
295 /*!
296 * \brief Define Head losses to take into account deposit in the flow
297 *
298 * \param[in] n_hl_cells number of cells on which to apply head losses
299 * \param[in] cell_ids ids of cells on which to apply head losses
300 * \param[in] bc_type boundary face type
301 * \param[out] cku head loss coefficients at matchin cells
302 */
303 /*----------------------------------------------------------------------------*/
304
305 void
cs_lagr_head_losses(cs_lnum_t n_hl_cells,const cs_lnum_t cell_ids[],const cs_lnum_t bc_type[],cs_real_t cku[][6])306 cs_lagr_head_losses(cs_lnum_t n_hl_cells,
307 const cs_lnum_t cell_ids[],
308 const cs_lnum_t bc_type[],
309 cs_real_t cku[][6])
310 {
311 cs_lnum_t ncel = cs_glob_mesh->n_cells;
312 cs_lnum_t ncelet = cs_glob_mesh->n_cells_with_ghosts;
313
314 cs_real_t *volume = cs_glob_mesh_quantities->cell_vol;
315
316 cs_lagr_extra_module_t *extra = cs_glob_lagr_extra_module;
317
318 /* Check that head loss zone definitions are consistent */
319 if (n_hl_cells != ncel)
320 cs_parameters_error
321 (CS_ABORT_IMMEDIATE,
322 _("in Lagrangian module"),
323 _("The number of cells in the head loss zones must cover\n"
324 "the whole mesh (though the local head loss may be zero).\n"));
325
326 /* ==================================================================
327 * cku: compute head loss coefficients in the calculation coordinates,
328 * organized in order k11, k22, k33, k12, k13, k23
329 * Note:
330 * - make sure diagonal coefficients are positive. The calculation
331 * may crash if this is not the case, and no further check will
332 * be done
333 * ==================================================================*/
334
335 /* ===================================================================
336 * Porosity calculation for the influence of the deposit on the flow
337 * by head losses
338 * ====================================================================*/
339
340 cs_real_t *mdiam;
341 BFT_MALLOC(mdiam, ncelet, cs_real_t);
342
343 /* cs_lnum_t poro_id; */
344 /* field_get_id_try ("clogging_porosity", &poro_id); */
345
346 /* TODO this field is never built */
347 cs_field_t *f_poro = cs_field_by_name_try("clogging_porosity");
348
349 cs_real_t *lporo;
350 if (f_poro == NULL)
351 BFT_MALLOC(lporo, ncelet, cs_real_t);
352 else
353 lporo = f_poro->val;
354
355 _porcel(mdiam, lporo, bc_type);
356
357 /* Calculation of the head loss term with the Ergun law
358 * mdiam : mean diameter of deposited particles
359 * lcell : characteristic length in the flow direction */
360
361 for (cs_lnum_t ielpdc = 0; ielpdc < n_hl_cells; ielpdc++) {
362
363 cs_lnum_t iel = cell_ids[ielpdc];
364
365 if (mdiam[iel] > 0.0) {
366
367 cs_real_t lcell = pow(volume[iel], 1.0 / 3.0);
368 cs_real_t romf = extra->cromf->val[iel];
369 cs_real_t visccf = extra->viscl->val[iel] / romf;
370 cs_real_t v = sqrt ( pow(extra->vel->vals[1][iel * 3 + 0], 2)
371 + pow(extra->vel->vals[1][iel * 3 + 1], 2)
372 + pow(extra->vel->vals[1][iel * 3 + 2], 2));
373
374 cs_real_t ck = v * 1.75 * (1 - lporo[iel]) / pow (lporo[iel], 3.0)
375 * lcell / mdiam[iel]
376 + (lcell * 150.0 * visccf) / (romf * pow(mdiam[iel], 2))
377 * pow((1 - lporo[iel]), 2) / lporo[iel] * 3;
378 cku[iel][0] = ck;
379 cku[iel][1] = ck;
380 cku[iel][2] = ck;
381 cku[iel][3] = 0.0;
382 cku[iel][4] = 0.0;
383 cku[iel][5] = 0.0;
384
385 }
386
387 }
388
389 if (f_poro == NULL)
390 BFT_FREE(lporo);
391
392 BFT_FREE(mdiam);
393 }
394
395 /*----------------------------------------------------------------------------*/
396
397 END_C_DECLS
398