1 /*============================================================================
2 * Fan modeling through velocity source terms.
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 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *----------------------------------------------------------------------------*/
32
33 #include <assert.h>
34 #include <math.h>
35 #include <stdio.h>
36 #include <stdlib.h>
37 #include <string.h>
38
39 /*----------------------------------------------------------------------------
40 * Local headers
41 *----------------------------------------------------------------------------*/
42
43 #include "bft_mem.h"
44
45 #include "cs_base.h"
46 #include "cs_field.h"
47 #include "cs_log.h"
48 #include "cs_math.h"
49 #include "cs_parall.h"
50
51 /*----------------------------------------------------------------------------
52 * Header for the current file
53 *----------------------------------------------------------------------------*/
54
55 #include "cs_fan.h"
56
57 /*----------------------------------------------------------------------------*/
58
59 BEGIN_C_DECLS
60
61 /*=============================================================================
62 * Additional doxygen documentation
63 *============================================================================*/
64
65 /*!
66 \file cs_fan.c
67 Fan modeling through velocity source terms.
68 */
69
70 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
71
72 /*============================================================================
73 * Local Type Definitions
74 *============================================================================*/
75
76 /* Structure associated to a fan */
77
78 struct _cs_fan_t {
79
80 int id; /* Fan id */
81 int dim; /* 2D or 3D geometry */
82 int mode; /* Use of the fan as a wind turbine (1)
83 or as a fan (0)(default) */
84 double inlet_axis_coords[3]; /* Axis point coordinates of the
85 inlet face */
86 double outlet_axis_coords[3]; /* Axis point coordinates of the
87 outlet face */
88 double axis_dir[3]; /* Unit vector of the axis
89 (inlet to outlet) */
90 double thickness; /* Fan thickness */
91 double surface; /* Fan total surface */
92 double volume; /* Fan total volume */
93 double volume_expected; /* Fan theoretical volume */
94
95 double fan_radius; /* Fan radius */
96 double blades_radius; /* Blades radius */
97 double hub_radius; /* Hub radius */
98 double curve_coeffs[3]; /* Coefficients of the terms of
99 degree 0, 1 and 2 of the
100 pressure drop/flow rate
101 characteristic curve */
102 double axial_torque; /* Fan axial torque */
103
104 cs_lnum_t n_cells; /* Number of cells */
105
106 cs_lnum_t *cell_list; /* List of the cells belonging
107 to the fan */
108
109 double in_flow; /* Current inlet flow */
110 double out_flow; /* Current outlet flow */
111 double delta_p; /* Pressure drop */
112
113 };
114
115 /*============================================================================
116 * Global variables
117 *============================================================================*/
118
119 /* Fans array */
120
121 static cs_lnum_t _cs_glob_n_fans_max = 0;
122
123 static cs_lnum_t _cs_glob_n_fans = 0;
124 static cs_fan_t ** _cs_glob_fans = NULL;
125
126 /*============================================================================
127 * Macro definitions
128 *============================================================================*/
129
130 enum {X, Y, Z};
131
132 /*============================================================================
133 * Private function definitions
134 *============================================================================*/
135
136 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
137
138 /*============================================================================
139 * Public function definitions for Fortran API
140 *============================================================================*/
141
142 /*----------------------------------------------------------------------------
143 * Compute the flows through the fans
144 *
145 * Fortran interface:
146 *
147 * subroutine debvtl
148 * *****************
149 *
150 * double precision flumas(*) : <-- : interior faces mass flux
151 * double precision flumab(*) : <-- : boundary faces mass flux
152 * double precision rhofac(*) : <-- : density at cells
153 * double precision rhofab(*) : <-- : density at boundary faces
154 *----------------------------------------------------------------------------*/
155
CS_PROCF(debvtl,DEBVTL)156 void CS_PROCF (debvtl, DEBVTL)
157 (
158 cs_real_t flumas[],
159 cs_real_t flumab[],
160 cs_real_t rho[],
161 cs_real_t rhofab[]
162 )
163 {
164 cs_fan_compute_flows(cs_glob_mesh,
165 cs_glob_mesh_quantities,
166 flumas,
167 flumab,
168 rho,
169 rhofab);
170 }
171
172 /*----------------------------------------------------------------------------
173 * Compute the force induced by the fans (needs a previous calculation
174 * of the flows through each fan).
175 *
176 * The induced force is added to the array crvxep (which can have other
177 * contributions).
178 *
179 * Fortran interface:
180 *
181 * subroutine tsvvtl
182 * *****************
183 *
184 * parameters:
185 * crvexp <-> Explicit source term (velocity)
186 *----------------------------------------------------------------------------*/
187
CS_PROCF(tsvvtl,TSVVTL)188 void CS_PROCF (tsvvtl, TSVVTL)
189 (
190 cs_real_3_t crvexp[]
191 )
192 {
193 cs_fan_compute_force(cs_glob_mesh_quantities,
194 crvexp);
195 }
196
197 /*============================================================================
198 * Public function definitions
199 *============================================================================*/
200
201 /*----------------------------------------------------------------------------*/
202 /*!
203 * \brief Fan definition (added to the ones previously defined)
204 *
205 * Fans are handled as explicit momentum source terms at the given location,
206 * based on the fan's axis and diameter.
207 * The fan's pressure characteristic curve is defined by 3 coefficients,
208 * such that:
209 * \f$\delta P = C_0 + C_1.flow + C_2.flow^2\f$.
210 * An axial torque may also be defined for the 3D model.
211 *
212 * \param[in] fan_dim fan dimension:
213 * 2: pseudo-2D (extruded mesh)
214 * 3: 3D (standard)
215 * \param[in] mode mode:
216 * 0: fan
217 * 1: wind turbine
218 * \param[in] inlet_axis_coords intersection of axis and inlet face
219 * \param[in] outlet_axis_coords intersection of axis and outlet face
220 * \param[in] fan_radius fan radius
221 * \param[in] blades_radius blades radius
222 * \param[in] hub_radius hub radius
223 * \param[in] curve_coeffs coefficients of degre 0, 1 and 2 of
224 * the pressure drop/flow rate
225 characteristic curve
226 * \param[in] axial_torque fan axial torque
227 */
228 /*----------------------------------------------------------------------------*/
229
230 void
cs_fan_define(int fan_dim,int mode,const cs_real_t inlet_axis_coords[3],const cs_real_t outlet_axis_coords[3],cs_real_t fan_radius,cs_real_t blades_radius,cs_real_t hub_radius,const cs_real_t curve_coeffs[3],cs_real_t axial_torque)231 cs_fan_define(int fan_dim,
232 int mode,
233 const cs_real_t inlet_axis_coords[3],
234 const cs_real_t outlet_axis_coords[3],
235 cs_real_t fan_radius,
236 cs_real_t blades_radius,
237 cs_real_t hub_radius,
238 const cs_real_t curve_coeffs[3],
239 cs_real_t axial_torque)
240 {
241 cs_fan_t *fan = NULL;
242
243 /* Define a new fan */
244
245 BFT_MALLOC(fan, 1, cs_fan_t);
246
247 fan->id = _cs_glob_n_fans;
248
249 fan->dim = fan_dim;
250
251 fan->mode = mode;
252
253 for (int i = 0; i < 3; i++) {
254 fan->inlet_axis_coords[i] = inlet_axis_coords[i];
255 fan->outlet_axis_coords[i] = outlet_axis_coords[i];
256 }
257
258 fan->fan_radius = fan_radius;
259 fan->blades_radius = blades_radius;
260 fan->hub_radius = hub_radius;
261
262 for (int i = 0; i < 3; i++)
263 fan->curve_coeffs[i] = curve_coeffs[i];
264 fan->axial_torque = axial_torque;
265
266 fan->n_cells = 0;
267 fan->cell_list = NULL;
268
269 /* Compute the axis vector */
270
271 fan->thickness = 0.0;
272
273 for (int i = 0; i < 3; i++) {
274 fan->axis_dir[i] = outlet_axis_coords[i] - inlet_axis_coords[i];
275 }
276
277 fan->thickness = cs_math_3_norm(fan->axis_dir);
278 cs_math_3_normalize(fan->axis_dir, fan->axis_dir);
279
280 /* Compute fan theoretical volume */
281
282 fan->volume_expected = cs_math_pi*cs_math_pow2(fan->fan_radius)*fan->thickness;
283
284 /* Surface/volume initialized to 0, will be set by cs_fan_build_all */
285
286 fan->surface = 0.0;
287 fan->volume = 0.0;
288
289 /* Flows initialized to 0 */
290
291 fan->in_flow = 0.0;
292 fan->out_flow = 0.0;
293
294 /* Increase the fans array if necessary */
295
296 if (_cs_glob_n_fans == _cs_glob_n_fans_max) {
297 _cs_glob_n_fans_max = (_cs_glob_n_fans_max + 1) * 2;
298 BFT_REALLOC(_cs_glob_fans, _cs_glob_n_fans_max, cs_fan_t *);
299 }
300
301 /* Adds in the fans array */
302
303 _cs_glob_fans[_cs_glob_n_fans] = fan;
304 _cs_glob_n_fans += 1;
305 }
306
307 /*----------------------------------------------------------------------------*/
308 /*!
309 * \brief Destroy the structures associated with fans.
310 */
311 /*----------------------------------------------------------------------------*/
312
313 void
cs_fan_destroy_all(void)314 cs_fan_destroy_all(void)
315 {
316 for (int i = 0; i < _cs_glob_n_fans; i++) {
317 cs_fan_t *fan = _cs_glob_fans[i];
318 BFT_FREE(fan->cell_list);
319 BFT_FREE(fan);
320 }
321
322 _cs_glob_n_fans_max = 0;
323 _cs_glob_n_fans = 0;
324 BFT_FREE(_cs_glob_fans);
325 }
326
327 /*----------------------------------------------------------------------------*/
328 /*!
329 * \brief Return number of fans.
330 *
331 * \return number of defined fans
332 */
333 /*----------------------------------------------------------------------------*/
334
335 int
cs_fan_n_fans(void)336 cs_fan_n_fans(void)
337 {
338 return _cs_glob_n_fans;
339 }
340
341 /*----------------------------------------------------------------------------*/
342 /*!
343 * \brief Log fans definition setup information.
344 */
345 /*----------------------------------------------------------------------------*/
346
347 void
cs_fan_log_setup(void)348 cs_fan_log_setup(void)
349 {
350 if (_cs_glob_n_fans < 1)
351 return;
352
353 cs_log_printf(CS_LOG_SETUP,
354 _("\n"
355 "Fans\n"
356 "----\n"));
357
358 for (int i = 0; i < _cs_glob_n_fans; i++) {
359 cs_fan_t *fan = _cs_glob_fans[i];
360 cs_log_printf
361 (CS_LOG_SETUP,
362 _(" Fan id: %d\n"
363 " Fan mesh dimension: %d\n"
364 " Wind turbine: %d\n"
365 " Axis coordinates: [%11.4e, %11.4e, %11.4e,\n"
366 " %11.4e, %11.4e, %11.4e]\n"
367 " Fan radius: %11.4e\n"
368 " Blades radius: %11.4e\n"
369 " Hub radius: %11.4e\n"
370 " Curve coefficients: C0: %10.3e, C1: %10.3e, C2: %10.3e\n"
371 " Axial torque: %10.3e\n"),
372 fan->id, fan->dim, fan->mode,
373 fan->inlet_axis_coords[0],
374 fan->inlet_axis_coords[1],
375 fan->inlet_axis_coords[2],
376 fan->outlet_axis_coords[0],
377 fan->outlet_axis_coords[1],
378 fan->outlet_axis_coords[2],
379 fan->fan_radius, fan->blades_radius, fan->hub_radius,
380 fan->curve_coeffs[0],
381 fan->curve_coeffs[1],
382 fan->curve_coeffs[2],
383 fan->axial_torque);
384 }
385 }
386
387 /*----------------------------------------------------------------------------*/
388 /*!
389 * \brief Log fan information for a given iteration.
390 */
391 /*----------------------------------------------------------------------------*/
392
393 void
cs_fan_log_iteration(void)394 cs_fan_log_iteration(void)
395 {
396 if (_cs_glob_n_fans < 1)
397 return;
398
399 cs_log_printf(CS_LOG_DEFAULT,
400 _("\n"
401 "Fans\n"
402 "----\n"));
403
404 cs_log_printf(CS_LOG_DEFAULT,
405 _(" id surface volume(real) volume(th.)"
406 " flow deltaP\n"
407 " ---- ----------- ----------- -----------"
408 " --------- -----------\n" ));
409
410 for (int i = 0; i < _cs_glob_n_fans; i++) {
411 cs_fan_t *fan = _cs_glob_fans[i];
412 cs_log_printf(CS_LOG_DEFAULT,
413 " %5d %11.4e %11.4e %11.4e %11.4e %11.4e\n",
414 fan->id, fan->surface, fan->volume, fan->volume_expected,
415 0.5*(fan->out_flow - fan->in_flow),
416 fan->delta_p);
417 }
418 }
419
420 /*----------------------------------------------------------------------------*/
421 /*!
422 * \brief Define the cells belonging to the different fans.
423 *
424 * \param[in] mesh associated mesh structure
425 * \param[in] mesh_quantities mesh quantities
426 */
427 /*----------------------------------------------------------------------------*/
428
429 void
cs_fan_build_all(const cs_mesh_t * mesh,const cs_mesh_quantities_t * mesh_quantities)430 cs_fan_build_all(const cs_mesh_t *mesh,
431 const cs_mesh_quantities_t *mesh_quantities)
432 {
433 cs_real_t coo_axis;
434 cs_real_t d_cel_axis[3];
435 cs_real_t l_surf;
436
437 cs_fan_t *fan = NULL;
438 cs_lnum_t *cpt_cel_vtl = NULL;
439 int *cell_fan_id = NULL;
440
441 const cs_lnum_t n_cells = mesh->n_cells;
442 const cs_lnum_t n_cells_ext = mesh->n_cells_with_ghosts;
443 const cs_lnum_2_t *i_face_cells = (const cs_lnum_2_t *)(mesh->i_face_cells);
444 const cs_lnum_t *b_face_cells = mesh->b_face_cells;
445 const cs_real_3_t *restrict cell_cen
446 = (const cs_real_3_t *restrict)mesh_quantities->cell_cen;
447 const cs_real_3_t *restrict i_face_normal
448 = (const cs_real_3_t *restrict)mesh_quantities->i_face_normal;
449 const cs_real_3_t *restrict b_face_normal
450 = (const cs_real_3_t *restrict)mesh_quantities->b_face_normal;
451
452 /* Reset fans in case already built */
453
454 for (int fan_id = 0; fan_id < _cs_glob_n_fans; fan_id++) {
455 fan = _cs_glob_fans[fan_id];
456 fan->n_cells = 0;
457 fan->surface = 0;
458 fan->volume = 0;
459 }
460
461 /* Create an array for cells flaging */
462 /*-----------------------------------*/
463
464 BFT_MALLOC(cell_fan_id, n_cells_ext, int);
465
466 for (cs_lnum_t cell_id = 0; cell_id < n_cells_ext; cell_id++)
467 cell_fan_id[cell_id] = -1;
468
469 /* Main loop on cells */
470
471 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
472
473 /* Loop on fans */
474
475 for (int fan_id = 0; fan_id < _cs_glob_n_fans; fan_id++) {
476
477 fan = _cs_glob_fans[fan_id];
478
479 /* Vector from the outlet face axis point to the cell center */
480
481 for (int coo_id = 0; coo_id < 3; coo_id++) {
482 d_cel_axis[coo_id] = (cell_cen[cell_id][coo_id])
483 - fan->inlet_axis_coords[coo_id];
484 }
485
486 /* Dot product with the axis vector */
487
488 coo_axis = cs_math_3_dot_product(d_cel_axis, fan->axis_dir);
489
490 /* Cell potentially in the fan if its center projection on the axis
491 is within the thickness */
492
493 if (coo_axis >= 0. && coo_axis <= fan->thickness) {
494
495 /* Projection of the vector from the outlet face axis point
496 to the cell center in the fan plane */
497
498 for (int coo_id = 0; coo_id < 3; coo_id++)
499 d_cel_axis[coo_id] -= coo_axis * fan->axis_dir[coo_id];
500
501 /* Distance to the axis */
502 cs_real_t d_axis = cs_math_3_norm(d_cel_axis);
503
504 /* If the cell is in the fan */
505
506 if (d_axis <= fan->fan_radius) {
507
508 cell_fan_id[cell_id] = fan_id;
509 fan->n_cells += 1;
510 fan->volume += mesh_quantities->cell_vol[cell_id];
511 }
512
513 }
514
515 } /* End of loop on fans */
516
517 } /* End of main loop on cells */
518
519 for (int fan_id = 0; fan_id < _cs_glob_n_fans; fan_id++) {
520 cs_parall_sum(1, CS_DOUBLE, &((_cs_glob_fans[fan_id])->volume));
521 }
522
523 /* Synchronize cell_fan_id */
524 if (mesh->halo != NULL)
525 cs_halo_sync_untyped(mesh->halo,
526 CS_HALO_EXTENDED,
527 sizeof(int),
528 cell_fan_id);
529
530 /* Create the lists of cells belonging to each fan */
531 /*-------------------------------------------------*/
532
533 BFT_MALLOC(cpt_cel_vtl, _cs_glob_n_fans, cs_lnum_t);
534
535 for (int fan_id = 0; fan_id < _cs_glob_n_fans; fan_id++) {
536
537 fan = _cs_glob_fans[fan_id];
538
539 BFT_REALLOC(fan->cell_list, fan->n_cells, cs_lnum_t);
540
541 cpt_cel_vtl[fan_id] = 0;
542 }
543
544 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
545
546 if (cell_fan_id[cell_id] > -1) {
547 int fan_id = cell_fan_id[cell_id];
548 fan = _cs_glob_fans[fan_id];
549 fan->cell_list[cpt_cel_vtl[fan_id]] = cell_id;
550 cpt_cel_vtl[fan_id] += 1;
551 }
552
553 }
554
555 #if defined(DEBUG) && !defined(NDEBUG)
556 for (int fan_id = 0; fan_id < _cs_glob_n_fans; fan_id++) {
557 fan = _cs_glob_fans[fan_id];
558 assert(cpt_cel_vtl[fan_id] == fan->n_cells);
559 }
560 #endif
561
562 /* Compute each fan surface */
563 /*--------------------------*/
564
565 /* Contribution to the domain interior */
566
567 for (cs_lnum_t face_id = 0; face_id < mesh->n_i_faces; face_id++) {
568
569 cs_lnum_t cell_id_1 = i_face_cells[face_id][0];
570 cs_lnum_t cell_id_2 = i_face_cells[face_id][1];
571
572 if ( cell_id_1 < mesh->n_cells /* ensure the contrib is from one domain */
573 && cell_fan_id[cell_id_1] != cell_fan_id[cell_id_2]) {
574
575 l_surf = cs_math_3_norm((i_face_normal[face_id]));
576 if (cell_fan_id[cell_id_1] > -1) {
577 int fan_id = cell_fan_id[cell_id_1];
578 fan = _cs_glob_fans[fan_id];
579 fan->surface += l_surf;
580 }
581 if (cell_fan_id[cell_id_2] > -1) {
582 int fan_id = cell_fan_id[cell_id_2];
583 fan = _cs_glob_fans[fan_id];
584 fan->surface += l_surf;
585 }
586 }
587 }
588
589 /* Contribution to the domain boundary */
590
591 for (cs_lnum_t face_id = 0; face_id < mesh->n_b_faces; face_id++) {
592
593 if (cell_fan_id[b_face_cells[face_id]] > -1) {
594 l_surf = cs_math_3_norm((b_face_normal[face_id]));
595 int fan_id = cell_fan_id[b_face_cells[face_id]];
596 fan = _cs_glob_fans[fan_id];
597 fan->surface += l_surf;
598 }
599 }
600
601 for (int fan_id = 0; fan_id < _cs_glob_n_fans; fan_id++)
602 cs_parall_sum(1, CS_DOUBLE, &((_cs_glob_fans[fan_id])->surface));
603
604 /* Free memory */
605
606
607 BFT_FREE(cpt_cel_vtl);
608 BFT_FREE(cell_fan_id);
609 }
610
611 /*----------------------------------------------------------------------------*/
612 /*!
613 * \brief Compute the flows through the fans.
614 *
615 * \param[in] mesh mesh structure
616 * \param[in] mesh_quantities mesh quantities
617 * \param[in] i_mass_flux interior faces mass flux
618 * \param[in] b_mass_flux boundary faces mass flux
619 * \param[in] c_rho density at cells
620 * \param[in] b_rho density at boundary faces
621 */
622 /*----------------------------------------------------------------------------*/
623
624 void
cs_fan_compute_flows(const cs_mesh_t * mesh,const cs_mesh_quantities_t * mesh_quantities,const cs_real_t i_mass_flux[],const cs_real_t b_mass_flux[],const cs_real_t c_rho[],const cs_real_t b_rho[])625 cs_fan_compute_flows(const cs_mesh_t *mesh,
626 const cs_mesh_quantities_t *mesh_quantities,
627 const cs_real_t i_mass_flux[],
628 const cs_real_t b_mass_flux[],
629 const cs_real_t c_rho[],
630 const cs_real_t b_rho[])
631 {
632 cs_lnum_t direction;
633
634 cs_real_t flow;
635
636 cs_fan_t *fan = NULL;
637 int *cell_fan_id = NULL;
638
639 const cs_lnum_t n_cells_ext = mesh->n_cells_with_ghosts;
640 const cs_lnum_t nbr_fac = mesh->n_i_faces;
641 const cs_lnum_t nbr_fbr = mesh->n_b_faces;
642 const cs_lnum_2_t *i_face_cells = (const cs_lnum_2_t *)(mesh->i_face_cells);
643 const cs_lnum_t *b_face_cells = mesh->b_face_cells;
644 const cs_real_3_t *restrict i_face_normal
645 = (const cs_real_3_t *restrict)mesh_quantities->i_face_normal;
646 const cs_real_3_t *restrict b_face_normal
647 = (const cs_real_3_t *restrict)mesh_quantities->b_face_normal;
648
649 /* Flag the cells */
650
651 BFT_MALLOC(cell_fan_id, n_cells_ext, int);
652
653 cs_fan_flag_cells(mesh, cell_fan_id);
654
655 /* Set the fans flows to zero */
656
657 for (int fan_id = 0; fan_id < _cs_glob_n_fans; fan_id++) {
658 fan = _cs_glob_fans[fan_id];
659 fan->in_flow = 0.0;
660 fan->out_flow = 0.0;
661 }
662
663 /* Contribution to the domain interior */
664
665 for (cs_lnum_t face_id = 0; face_id < nbr_fac; face_id++) {
666
667 cs_lnum_t cell_id_1 = i_face_cells[face_id][0];
668 cs_lnum_t cell_id_2 = i_face_cells[face_id][1];
669
670 if ( cell_id_1 < mesh->n_cells /* Make sure the contrib is from one domain */
671 && cell_fan_id[cell_id_1] != cell_fan_id[cell_id_2]) {
672
673 for (int i = 0; i < 2; i++) {
674
675 cs_lnum_t cell_id = i_face_cells[face_id][i];
676 int fan_id = cell_fan_id[cell_id];
677
678 if (fan_id > -1) {
679 fan = _cs_glob_fans[fan_id];
680 direction = (i == 0 ? 1 : - 1);
681 flow = i_mass_flux[face_id]/c_rho[cell_id] * direction;
682 if ( cs_math_3_dot_product(fan->axis_dir, i_face_normal[face_id])
683 * direction > 0.0)
684 fan->out_flow += flow;
685 else
686 fan->in_flow += flow;
687 }
688 }
689
690 }
691
692 }
693
694 /* Contribution to the domain boundary */
695
696 for (cs_lnum_t face_id = 0; face_id < nbr_fbr; face_id++) {
697
698 int fan_id = cell_fan_id[b_face_cells[face_id]];
699
700 if (fan_id > -1) {
701
702 fan = _cs_glob_fans[fan_id];
703
704 flow = b_mass_flux[face_id]/b_rho[face_id];
705 if (cs_math_3_dot_product(fan->axis_dir, b_face_normal[face_id]) > 0.0)
706 fan->out_flow += flow;
707 else
708 fan->in_flow += flow;
709
710 }
711 }
712
713 #if defined(HAVE_MPI)
714 if (cs_glob_n_ranks > 1) {
715
716 for (int fan_id = 0; fan_id < _cs_glob_n_fans; fan_id++) {
717
718 fan = _cs_glob_fans[fan_id];
719
720 cs_real_t flow_glob[2];
721 cs_real_t flow_loc[2] = {fan->out_flow, fan->in_flow};
722
723 MPI_Allreduce(flow_loc, flow_glob, 2, CS_MPI_REAL, MPI_SUM,
724 cs_glob_mpi_comm);
725
726 fan->out_flow = flow_glob[0];
727 fan->in_flow = flow_glob[1];
728
729 }
730 }
731 #endif
732
733 /* In 2D, the flow is normalized by the surface */
734
735 if (fan->dim == 2) {
736 cs_real_t surf_2d;
737 surf_2d = (0.5*fan->surface - 2*fan->fan_radius*fan->thickness)
738 / (2*fan->fan_radius+fan->thickness);
739 fan->out_flow = fan->out_flow / surf_2d;
740 fan->in_flow = fan->in_flow / surf_2d;
741 }
742
743 /* Free memory */
744
745 BFT_FREE(cell_fan_id);
746 }
747
748 /*----------------------------------------------------------------------------*/
749 /*!
750 * \brief Compute the force induced by the fans
751 * (needs a previous calculation of the flows through each fan).
752 *
753 * \param[in] mesh_quantities mesh quantities
754 * \param[in] source_t explicit source term for the velocity
755 */
756 /*----------------------------------------------------------------------------*/
757
758 void
cs_fan_compute_force(const cs_mesh_quantities_t * mesh_quantities,cs_real_3_t source_t[])759 cs_fan_compute_force(const cs_mesh_quantities_t *mesh_quantities,
760 cs_real_3_t source_t[])
761 {
762 cs_real_t f_z, f_theta;
763 cs_real_t f_rot[3];
764
765 const cs_real_3_t *restrict cell_cen
766 = (const cs_real_3_t *restrict)mesh_quantities->cell_cen;
767 const cs_real_t *cell_f_vol = mesh_quantities->cell_f_vol;
768 const cs_real_t pi = 4.*atan(1.);
769
770 /* Compute the force induced by fans */
771
772 /* Loop on fans */
773 /*--------------*/
774
775 for (int fan_id = 0; fan_id < _cs_glob_n_fans; fan_id++) {
776
777 cs_fan_t *fan = _cs_glob_fans[fan_id];
778
779 const cs_real_t r_hub = fan->hub_radius;
780 const cs_real_t r_blades = fan->blades_radius;
781 const cs_real_t r_fan = fan->fan_radius;
782
783 const cs_real_t mean_flow = 0.5 * (fan->out_flow - fan->in_flow);
784
785 fan->delta_p = (fan->curve_coeffs[2] * mean_flow*mean_flow)
786 + (fan->curve_coeffs[1] * mean_flow)
787 + (fan->curve_coeffs[0]);
788
789 /* Loop on fan cells */
790 /*-------------------*/
791
792 for (cs_lnum_t i = 0; i < fan->n_cells; i++) {
793
794 cs_lnum_t cell_id = fan->cell_list[i];
795
796 f_z = 0.0;
797 f_theta = 0.0;
798 f_rot[0] = 0.0, f_rot[1] = 0.0, f_rot[2] = 0.0;
799
800 if (r_blades < 1.0e-12 && r_hub < 1.0e-12) {
801
802 f_z = fan->delta_p / fan->thickness;
803 f_theta = 0.0;
804
805 }
806 else if (r_hub < r_blades) {
807
808 cs_real_t r_1, r_2, aux_1, aux_2, coo_axis, d_axis, d_cel_axis[3];
809
810 r_1 = 0.7 * fan->blades_radius;
811 r_2 = 0.85 * fan->blades_radius;
812
813 if (fan->dim == 2) {
814 if (fan->mode == 1) {
815 aux_1 = - (fan->delta_p * 2.0 * r_fan)
816 / (fan->thickness * (1.15*r_blades - r_hub));
817 }
818 else{
819 aux_1 = (fan->delta_p * 2.0 * r_fan)
820 / (fan->thickness * (1.15*r_blades - r_hub));
821 }
822 aux_2 = 0.0;
823 }
824 else {
825 const cs_real_t r_hub4 = r_hub * r_hub * r_hub * r_hub;
826 const cs_real_t r_hub3 = r_hub * r_hub * r_hub;
827 const cs_real_t r_blades4 = r_blades * r_blades * r_blades * r_blades;
828 const cs_real_t r_blades3 = r_blades * r_blades * r_blades;
829 const cs_real_t r_blades2 = r_blades * r_blades;
830 const cs_real_t r_fan2 = r_fan * r_fan;
831 cs_real_t f_base = (0.7*r_blades - r_hub)
832 / ( 1.0470*fan->thickness
833 * ( r_hub3
834 + 1.4560*r_blades3
835 - 2.570*r_blades2*r_hub));
836 cs_real_t f_orth = (0.7*r_blades - r_hub)
837 / ( fan->thickness
838 * ( 1.042*r_blades4
839 + 0.523*r_hub4
840 - 1.667*r_blades3*r_hub));
841 if (fan->mode == 1)
842 aux_1 = - f_base * fan->delta_p * pi * r_fan2;
843 else
844 aux_1 = f_base * fan->delta_p * pi * r_fan2;
845 aux_2 = f_orth * fan->axial_torque;
846 }
847
848 /* Vector from the outlet face axis point to the cell center */
849
850 for (int coo_id = 0; coo_id < 3; coo_id++) {
851 d_cel_axis[coo_id] = (cell_cen[cell_id][coo_id])
852 - fan->inlet_axis_coords[coo_id];
853 }
854
855 /* Projection of the cell center on the fan axis */
856
857 coo_axis = cs_math_3_dot_product(d_cel_axis, fan->axis_dir);
858
859 /* Projection of the vector from the outlet face axis point
860 to the cell center in the fan plane */
861
862 for (int coo_id = 0; coo_id < 3; coo_id++)
863 d_cel_axis[coo_id] -= coo_axis * fan->axis_dir[coo_id];
864
865 d_axis = cs_math_3_norm(d_cel_axis); /* Distance to the axis */
866
867 cs_math_3_cross_product(fan->axis_dir, d_cel_axis, f_rot);
868
869 cs_math_3_normalize(f_rot, f_rot);
870
871 if (d_axis < r_hub) {
872 f_z = 0.0;
873 f_theta = 0.0;
874 }
875 else if (d_axis < r_1) {
876 f_z = aux_1 * (d_axis - r_hub) / (r_1 - r_hub);
877 f_theta = aux_2 * (d_axis - r_hub) / (r_1 - r_hub);
878 }
879 else if (d_axis < r_2) {
880 f_z = aux_1;
881 f_theta = aux_2;
882 }
883 else if (d_axis < r_blades) {
884 f_z = aux_1 * (r_blades - d_axis) / (r_blades - r_2);
885 f_theta = aux_2 * (r_blades - d_axis) / (r_blades - r_2);
886 }
887 else {
888 f_z = 0.0;
889 f_theta = 0.0;
890 }
891
892 }
893
894 for (int coo_id = 0; coo_id < 3; coo_id++)
895 source_t[cell_id][coo_id]
896 += ( (f_z * fan->axis_dir[coo_id])
897 + (f_theta * f_rot[coo_id]))
898 * fan->volume_expected / fan->volume /* Correction factor */
899 * cell_f_vol[cell_id];
900
901 } /* End of loop on fan cells */
902
903 } /* End of loop on fans */
904
905 }
906
907 /*----------------------------------------------------------------------------*/
908 /*!
909 * \brief Flag the cells belonging to the different fans
910 * (by the fan id, -1 otherwise)
911 *
912 * \param[in] mesh assosiated mesh structure
913 * \param[out] cell_fan_id fan id (or -1) for each cell
914 */
915 /*----------------------------------------------------------------------------*/
916
917 void
cs_fan_flag_cells(const cs_mesh_t * mesh,int cell_fan_id[])918 cs_fan_flag_cells(const cs_mesh_t *mesh,
919 int cell_fan_id[])
920 {
921 cs_fan_t *fan;
922
923 const cs_lnum_t n_cells_ext = mesh->n_cells_with_ghosts;
924
925 /* Flag the cells */
926 for (cs_lnum_t cell_id = 0; cell_id < n_cells_ext; cell_id++)
927 cell_fan_id[cell_id] = -1;
928
929 for (int fan_id = 0; fan_id < _cs_glob_n_fans; fan_id++) {
930
931 fan = _cs_glob_fans[fan_id];
932
933 for (cs_lnum_t i = 0; i < fan->n_cells; i++) {
934 cs_lnum_t cell_id = fan->cell_list[i];
935 cell_fan_id[cell_id] = fan_id;
936 }
937
938 }
939 /* Synchronize cell_fan_id */
940 if (mesh->halo != NULL)
941 cs_halo_sync_untyped(mesh->halo,
942 CS_HALO_EXTENDED,
943 sizeof(int),
944 cell_fan_id);
945
946 /* Store the cell_fan_id in the postprocessing field */
947
948 cs_field_t *c_fan_id = cs_field_by_name("fan_id");
949 for (cs_lnum_t cell_id = 0; cell_id < n_cells_ext; cell_id++)
950 c_fan_id->val[cell_id] = (cs_real_t)cell_fan_id[cell_id];
951 }
952
953 /*----------------------------------------------------------------------------*/
954 /*!
955 * \brief Selection function for cells belonging to fans.
956 *
957 * This function may be used for the definition of postprocessing meshes.
958 *
959 * \param[in, out] input pointer to input (unused here)
960 * \param[out] n_cells number of selected cells
961 * \param[out] cell_ids array of selected cell ids (0 to n-1 numbering)
962 */
963 /*----------------------------------------------------------------------------*/
964
965 void
cs_fan_cells_select(void * input,cs_lnum_t * n_cells,cs_lnum_t ** cell_ids)966 cs_fan_cells_select(void *input,
967 cs_lnum_t *n_cells,
968 cs_lnum_t **cell_ids)
969 {
970 CS_UNUSED(input);
971
972 cs_lnum_t _n_cells;
973
974 int *cell_fan_id = NULL;
975 cs_lnum_t *_cell_ids = NULL;
976
977 const cs_mesh_t *m = cs_glob_mesh;
978
979 /* Preallocate selection list */
980
981 BFT_MALLOC(_cell_ids, m->n_cells, cs_lnum_t);
982
983 /* Allocate working array */
984
985 BFT_MALLOC(cell_fan_id, m->n_cells_with_ghosts, int);
986
987 /* Now flag cells and build list */
988
989 cs_fan_build_all(cs_glob_mesh, cs_glob_mesh_quantities);
990 cs_fan_flag_cells(m, cell_fan_id);
991
992 _n_cells = 0;
993
994 for (cs_lnum_t i = 0; i < m->n_cells; i++) {
995 if (cell_fan_id[i] > -1)
996 _cell_ids[_n_cells++] = i;
997 }
998
999 /* Free memory */
1000 BFT_FREE(cell_fan_id);
1001 BFT_REALLOC(_cell_ids, _n_cells, cs_lnum_t);
1002
1003 /* Set return values */
1004
1005 *n_cells = _n_cells;
1006 *cell_ids = _cell_ids;
1007 }
1008
1009 /*----------------------------------------------------------------------------*/
1010
1011 END_C_DECLS
1012