1 /*============================================================================
2 * Base electrical model data.
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
31 /*----------------------------------------------------------------------------
32 * Standard C library headers
33 *----------------------------------------------------------------------------*/
34
35 #include <assert.h>
36 #include <stdio.h>
37 #include <stdlib.h>
38 #include <string.h>
39 #include <math.h>
40
41 /*----------------------------------------------------------------------------
42 * Local headers
43 *----------------------------------------------------------------------------*/
44
45 #include "bft_mem.h"
46 #include "bft_error.h"
47 #include "bft_printf.h"
48
49 #include "cs_field.h"
50 #include "cs_field_default.h"
51 #include "cs_file.h"
52 #include "cs_log.h"
53 #include "cs_parall.h"
54 #include "cs_math.h"
55 #include "cs_mesh_quantities.h"
56 #include "cs_mesh_location.h"
57 #include "cs_time_step.h"
58 #include "cs_parameters.h"
59 #include "cs_field_pointer.h"
60 #include "cs_gradient.h"
61 #include "cs_field_operator.h"
62 #include "cs_physical_constants.h"
63 #include "cs_physical_model.h"
64 #include "cs_thermal_model.h"
65 #include "cs_turbulence_model.h"
66 #include "cs_gui_specific_physics.h"
67 #include "cs_gui_util.h"
68 #include "cs_post.h"
69 #include "cs_prototypes.h"
70
71 /*----------------------------------------------------------------------------
72 * Header for the current file
73 *----------------------------------------------------------------------------*/
74
75 #include "cs_elec_model.h"
76
77 /*----------------------------------------------------------------------------*/
78
79 BEGIN_C_DECLS
80
81 /*=============================================================================
82 * Additional doxygen documentation
83 *============================================================================*/
84
85 /*!
86 * \file cs_elec_model.c
87 *
88 * \brief Base electrical model data.
89 *
90 * Please refer to the
91 * <a href="../../theory.pdf#electric"><b>electric arcs</b></a>
92 * section of the theory guide for more informations.
93 */
94
95 /*----------------------------------------------------------------------------*/
96
97 /*! \struct cs_elec_option_t
98
99 \brief option for electric model
100
101 \var cs_elec_option_t::ixkabe
102 Model for radiative properties
103 - 0: last column read but not use
104 - 1: last column : absorption coefficient
105 - 2: last column : radiative ST
106 \var cs_elec_option_t::ntdcla
107 First iteration to take into account restrike model
108 \var cs_elec_option_t::irestrike
109 Indicate if restrike or not
110 \var cs_elec_option_t::restrike_point
111 Coordinates for restrike point
112 \var cs_elec_option_t::crit_reca
113 Defines plane coordinates component used to calculate
114 current in a plane.\n
115 Useful if \ref modrec = 2.
116 \var cs_elec_option_t::ielcor
117 Indicate if scaling or not.\n
118 When \ref ielcor = 1, the boundary conditions for the potential
119 will be tuned at each time step in order to reach a user-specified
120 target dissipated power \ref puisim (Joule effect) or a user-specified
121 target current intensity \ref couimp (electric arcs).\n The boundary
122 condition tuning is controlled by subroutines \ref elreca or
123 \ref cs_user_electric_scaling.
124 \var cs_elec_option_t::modrec
125 Model for scaling
126 - 1: volumic power for boundary conditions tuning,
127 - 2: by plane for boundary conditions tuning,
128 - 3: user function for boundary conditions tuning.
129 In this case, we need to define the plane and the
130 current density component used.
131 \var cs_elec_option_t::idreca
132 Defines the current density component used to calculate
133 current in a plane.\n
134 Useful if \ref modrec = 2.
135 \var cs_elec_option_t::izreca
136 Indicator for faces for scaling
137 \var cs_elec_option_t::couimp
138 Imposed current.\n
139 With the electric arcs module, \ref couimp is the target current
140 intensity (\f$A\f$) for the calculations with boundary condition
141 tuning for the potential.\n The target intensity will be reached
142 if the boundary conditions are expressed using the variable
143 \ref pot_diff or if the initial boundary conditions are multiplied by
144 the variable \ref coejou.\n
145 Useful with the electric arcs module if \ref ielcor = 1.
146 \var cs_elec_option_t::pot_diff
147 Potential difference.\n
148 \ref pot_diff is the potential difference (\f$V\f$) which generates
149 the current (and the Joule effect) for the calculations with boundary
150 conditions tuning for the potential. This value is initialised set by
151 the user (\ref cs_user_parameters). It is then automatically tuned
152 depending on the value of dissipated power (Joule effect module) or the
153 intensity of current (electric arcs module). In order for the correct
154 power or intensity to be reached, the boundary conditions for the
155 potential must be expressed with \ref pot_diff . The tuning can be
156 controlled in \ref cs_user_electric_scaling.\n
157 Useful if \ref ielcor = 1.
158 \var cs_elec_option_t::puisim
159 Imposed power.\n
160 With the Joule effect module, \ref puisim is the target dissipated power ($W$)
161 for the calculations with boundary condition tuning for the potential.\n
162 The target power will be reached if the boundary conditions are expressed
163 using the variable \ref pot_diff or if the initial boundary conditions are
164 multiplied by the variable \ref coejou .
165 Useful with the Joule effect module if \ref ielcor = 1.
166 \var cs_elec_option_t::coejou
167 coefficient for scaling
168 \var cs_elec_option_t::elcou
169 current in scaling plane
170 */
171
172 /*! \struct cs_data_joule_effect_t
173 ² \brief Structure to read transformer parameters in dp_ELE
174
175 \var cs_data_joule_effect_t::nbelec
176 transformer number
177 \var cs_data_joule_effect_t::ielecc
178 \var cs_data_joule_effect_t::ielect
179 \var cs_data_joule_effect_t::ielecb
180 \var cs_data_joule_effect_t::nbtrf
181 \var cs_data_joule_effect_t::ntfref
182 \var cs_data_joule_effect_t::ibrpr
183 \var cs_data_joule_effect_t::ibrsec
184 \var cs_data_joule_effect_t::tenspr
185 \var cs_data_joule_effect_t::rnbs
186 \var cs_data_joule_effect_t::zr
187 \var cs_data_joule_effect_t::zi
188 \var cs_data_joule_effect_t::uroff
189 \var cs_data_joule_effect_t::uioff
190 */
191
192 /*! \struct cs_data_elec_t
193
194 \brief physical properties for electric model descriptor.
195
196 \var cs_data_elec_t::ngaz
197 number of gaz in electrical data file
198 \var cs_data_elec_t::npoint
199 number of point in electrical data file for each gas
200 \var cs_data_elec_t::th
201 temperature values
202 \var cs_data_elec_t::ehgaz
203 enthalpy values
204 \var cs_data_elec_t::rhoel
205 density values
206 \var cs_data_elec_t::cpel
207 specific heat values
208 \var cs_data_elec_t::sigel
209 electric conductivity values
210 \var cs_data_elec_t::visel
211 dynamic viscosity
212 \var cs_data_elec_t::xlabel
213 thermal conductivity
214 \var cs_data_elec_t::xkabel
215 absorption coefficent
216 */
217
218 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
219
220 /*=============================================================================
221 * Macro definitions
222 *============================================================================*/
223
224 #define LG_MAX 1000
225
226 /*============================================================================
227 * Type definitions
228 *============================================================================*/
229
230 static cs_elec_option_t _elec_option = {.ixkabe = -1,
231 .ntdcla = -1,
232 .irestrike = -1,
233 .restrike_point = {0., 0., 0.},
234 .crit_reca = {0, 0, 0, 0, 0},
235 .ielcor = -1,
236 .modrec = -1,
237 .idreca = -1,
238 .izreca = NULL,
239 .couimp = 0.,
240 .pot_diff = 0.,
241 .puisim = 0.,
242 .coejou = 0.,
243 .elcou = 0.,
244 .srrom = 0.};
245
246 static cs_data_elec_t _elec_properties = {.ngaz = 0,
247 .npoint = 0,
248 .th = NULL,
249 .ehgaz = NULL,
250 .rhoel = NULL,
251 .cpel = NULL,
252 .sigel = NULL,
253 .visel = NULL,
254 .xlabel = NULL,
255 .xkabel = NULL};
256
257 static cs_data_joule_effect_t *_transformer = NULL;
258
259 const cs_elec_option_t *cs_glob_elec_option = NULL;
260 const cs_data_elec_t *cs_glob_elec_properties = NULL;
261 const cs_data_joule_effect_t *cs_glob_transformer = NULL;
262
263 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
264
265 /*!
266 * vacuum magnetic permeability constant (H/m). (= 1.2566e-6)
267 *
268 */
269 const cs_real_t cs_elec_permvi = 1.2566e-6;
270
271 /*!
272 * vacuum permittivity constant (F/m). (= 8.854e-12)
273 *
274 */
275 const cs_real_t cs_elec_epszer = 8.854e-12;
276
277 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
278
279 void
280 cs_f_elec_model_get_pointers(int **ngazge,
281 int **ielcor,
282 double **pot_diff,
283 double **coejou,
284 double **elcou,
285 double **couimp,
286 int **irestrike,
287 int **ntdcla,
288 double **restrike_point_x,
289 double **restrike_point_y,
290 double **restrike_point_z);
291
292 /*----------------------------------------------------------------------------
293 * Get pointers to members of the global electric model structure.
294 *
295 * This function is intended for use by Fortran wrappers, and
296 * enables mapping to Fortran global pointers.
297 *
298 * parameters:
299 * ngazge --> pointer to cs_glob_elec_properties->ngaz
300 * ielcor --> pointer to cs_glob_elec_option->ielcor
301 * pot_diff --> pointer to cs_glob_elec_option->pot_diff
302 * coejou --> pointer to cs_glob_elec_option->coejou
303 * elcou --> pointer to cs_glob_elec_option->elcou
304 * irestrike --> pointer to cs_glob_elec_option->irestrike
305 * ntdcla --> pointer to cs_glob_elec_option->ntdcla
306 * restrike_point --> pointer to cs_glob_elec_option->restrike_point
307 *----------------------------------------------------------------------------*/
308
309 void
cs_f_elec_model_get_pointers(int ** ngazge,int ** ielcor,double ** pot_diff,double ** coejou,double ** elcou,double ** couimp,int ** irestrike,int ** ntdcla,double ** restrike_point_x,double ** restrike_point_y,double ** restrike_point_z)310 cs_f_elec_model_get_pointers(int **ngazge,
311 int **ielcor,
312 double **pot_diff,
313 double **coejou,
314 double **elcou,
315 double **couimp,
316 int **irestrike,
317 int **ntdcla,
318 double **restrike_point_x,
319 double **restrike_point_y,
320 double **restrike_point_z)
321 {
322 *ngazge = &(_elec_properties.ngaz);
323 *ielcor = &(_elec_option.ielcor);
324 *pot_diff = &(_elec_option.pot_diff);
325 *coejou = &(_elec_option.coejou);
326 *elcou = &(_elec_option.elcou);
327 *couimp = &(_elec_option.couimp);
328 *irestrike = &(_elec_option.irestrike);
329 *ntdcla = &(_elec_option.ntdcla);
330 *restrike_point_x = &(_elec_option.restrike_point[0]);
331 *restrike_point_y = &(_elec_option.restrike_point[1]);
332 *restrike_point_z = &(_elec_option.restrike_point[2]);
333 }
334
335 /*============================================================================
336 * Private function definitions
337 *============================================================================*/
338
339 static void
_cs_electrical_model_verify(void)340 _cs_electrical_model_verify(void)
341 {
342 bool verif = true;
343
344 int ieljou = cs_glob_physical_model_flag[CS_JOULE_EFFECT];
345 int ielarc = cs_glob_physical_model_flag[CS_ELECTRIC_ARCS];
346
347 if (ielarc != -1 && ielarc != 2)
348 bft_error(__FILE__, __LINE__, 0,
349 _("Error for electric arc model\n"
350 "only choice -1 or 2 are permitted yet\n"
351 "model selected : \"%i\";\n"),
352 ielarc);
353
354 if ( ieljou != -1 && ieljou != 1 && ieljou != 2 && ieljou != 3
355 && ieljou != 4)
356 bft_error(__FILE__, __LINE__, 0,
357 _("Error for joule model\n"
358 "only choice -1, 1, 2, 3 or 4 are permitted yet\n"
359 "model selected : \"%i\";\n"),
360 ieljou);
361
362 /* options */
363 if (cs_glob_elec_option->ielcor != 0 && cs_glob_elec_option->ielcor != 1)
364 bft_error(__FILE__, __LINE__, 0,
365 _("Error for scaling model\n"
366 "only choice -1 or 2 are permitted yet\n"
367 "model selected : \"%i\";\n"),
368 cs_glob_elec_option->ielcor);
369
370 if (cs_glob_elec_option->ielcor == 1) {
371 if (ielarc > 0) {
372 if (cs_glob_elec_option->couimp < 0.) {
373 bft_printf("value for COUIMP must be strictly positive\n");
374 verif = false;
375 }
376 if (cs_glob_elec_option->pot_diff < 0.) {
377 bft_printf("value for DPOT must be strictly positive\n");
378 verif = false;
379 }
380 }
381 if (ieljou > 0) {
382 if (cs_glob_elec_option->puisim < 0.) {
383 bft_printf("value for PUISIM must be strictly positive\n");
384 verif = false;
385 }
386 if (cs_glob_elec_option->coejou < 0.) {
387 bft_printf("value for COEJOU must be strictly positive\n");
388 verif = false;
389 }
390 if (cs_glob_elec_option->pot_diff < 0.) {
391 bft_printf("value for DPOT must be strictly positive\n");
392 verif = false;
393 }
394 }
395 }
396
397 if (!verif) {
398 bft_error(__FILE__, __LINE__, 0,
399 _("Invalid or incomplete calculation parameter\n"
400 "Verify parameters\n"));
401 }
402 }
403
404 /*----------------------------------------------------------------------------*/
405 /*!
406 * \brief Map base fields to enumerated pointers for electric arcs
407 *
408 * \param[in] n_gasses number of gasses
409 */
410 /*----------------------------------------------------------------------------*/
411
412 static void
_field_pointer_map_electric_arcs(int n_gasses)413 _field_pointer_map_electric_arcs(int n_gasses)
414 {
415 char s[64];
416
417 cs_field_pointer_map(CS_ENUMF_(h),
418 cs_field_by_name_try("enthalpy"));
419
420 cs_field_pointer_map(CS_ENUMF_(potr), cs_field_by_name_try("elec_pot_r"));
421 cs_field_pointer_map(CS_ENUMF_(poti), cs_field_by_name_try("elec_pot_i"));
422
423 cs_field_pointer_map(CS_ENUMF_(potva), cs_field_by_name_try("vec_potential"));
424
425 for (int i = 0; i < n_gasses - 1; i++) {
426 snprintf(s, 63, "esl_fraction_%02d", i+1); s[63] = '\0';
427 cs_field_pointer_map_indexed(CS_ENUMF_(ycoel), i, cs_field_by_name_try(s));
428 }
429 }
430
431 /*----------------------------------------------------------------------------*/
432 /*!
433 * \brief Map base fields to enumerated pointers properties for electric arcs
434 */
435 /*----------------------------------------------------------------------------*/
436
437 static void
_field_pointer_properties_map_electric_arcs(void)438 _field_pointer_properties_map_electric_arcs(void)
439 {
440 cs_field_pointer_map(CS_ENUMF_(t),
441 cs_field_by_name_try("temperature"));
442
443 cs_field_pointer_map(CS_ENUMF_(joulp),
444 cs_field_by_name_try("joule_power"));
445 cs_field_pointer_map(CS_ENUMF_(radsc),
446 cs_field_by_name_try("radiation_source"));
447 cs_field_pointer_map(CS_ENUMF_(elech),
448 cs_field_by_name_try("elec_charge"));
449
450 cs_field_pointer_map(CS_ENUMF_(curre),
451 cs_field_by_name_try("current_re"));
452 cs_field_pointer_map(CS_ENUMF_(curim),
453 cs_field_by_name_try("current_im"));
454 cs_field_pointer_map(CS_ENUMF_(laplf),
455 cs_field_by_name_try("laplace_force"));
456 cs_field_pointer_map(CS_ENUMF_(magfl),
457 cs_field_by_name_try("magnetic_field"));
458 cs_field_pointer_map(CS_ENUMF_(elefl),
459 cs_field_by_name_try("electric_field"));
460 }
461
462 /*----------------------------------------------------------------------------
463 * convert temperature to enthalpy
464 *----------------------------------------------------------------------------*/
465
466 static cs_real_t
_cs_elec_convert_t_to_h(const cs_real_t ym[],cs_real_t temp)467 _cs_elec_convert_t_to_h(const cs_real_t ym[],
468 cs_real_t temp)
469 {
470 int ngaz = cs_glob_elec_properties->ngaz;
471 int it = cs_glob_elec_properties->npoint;
472
473 cs_real_t enthal = 0.;
474
475 if (temp >= cs_glob_elec_properties->th[it - 1]) {
476 for (int iesp = 0; iesp < ngaz; iesp++)
477 enthal += ym[iesp] * cs_glob_elec_properties->ehgaz[iesp * (it-1) + it-1];
478 }
479 else if (temp <= cs_glob_elec_properties->th[0]) {
480 for (int iesp = 0; iesp < ngaz; iesp++)
481 enthal += ym[iesp] * cs_glob_elec_properties->ehgaz[iesp * (it-1) + 0];
482 }
483 else {
484 for (int itt = 0; itt < cs_glob_elec_properties->npoint - 1; itt++) {
485 if ( temp > cs_glob_elec_properties->th[itt]
486 && temp <= cs_glob_elec_properties->th[itt + 1]) {
487 cs_real_t eh0 = 0.;
488 cs_real_t eh1 = 0.;
489
490 for (int iesp = 0; iesp < ngaz; iesp++) {
491 eh0 += ym[iesp] * cs_glob_elec_properties->ehgaz[iesp * (it-1) + itt];
492 eh1 += ym[iesp] * cs_glob_elec_properties->ehgaz[iesp * (it-1) + itt+1];
493 }
494
495 enthal = eh0 + (eh1 - eh0) * (temp - cs_glob_elec_properties->th[itt]) /
496 ( cs_glob_elec_properties->th[itt + 1]
497 - cs_glob_elec_properties->th[itt]);
498
499 break;
500 }
501 }
502 }
503
504 return enthal;
505 }
506
507 /*----------------------------------------------------------------------------
508 * Convert enthalpy to temperature
509 *----------------------------------------------------------------------------*/
510
511 static cs_real_t
_cs_elec_convert_h_to_t(const cs_real_t ym[],cs_real_t enthal)512 _cs_elec_convert_h_to_t(const cs_real_t ym[],
513 cs_real_t enthal)
514 {
515 int ngaz = cs_glob_elec_properties->ngaz;
516 int it = cs_glob_elec_properties->npoint;
517
518 cs_real_t eh1 = 0.;
519
520 for (int iesp = 0; iesp < ngaz; iesp++)
521 eh1 += ym[iesp] * cs_glob_elec_properties->ehgaz[iesp * (it-1) + it - 1];
522
523 if (enthal >= eh1) {
524 return cs_glob_elec_properties->th[it - 1];
525 }
526
527 eh1 = 0.;
528
529 for (int iesp = 0; iesp < ngaz; iesp++)
530 eh1 += ym[iesp] * cs_glob_elec_properties->ehgaz[iesp * (it-1) + 0];
531
532 if (enthal <= eh1) {
533 return cs_glob_elec_properties->th[0];
534 }
535
536 for (int itt = 0; itt < cs_glob_elec_properties->npoint - 1; itt++) {
537 cs_real_t eh0 = 0.;
538 eh1 = 0.;
539
540 for (int iesp = 0; iesp < ngaz; iesp++) {
541 eh0 += ym[iesp] * cs_glob_elec_properties->ehgaz[iesp * (it-1) + itt];
542 eh1 += ym[iesp] * cs_glob_elec_properties->ehgaz[iesp * (it-1) + itt+1];
543 }
544
545 if (enthal > eh0 && enthal <= eh1) {
546 cs_real_t temp = cs_glob_elec_properties->th[itt]
547 + (enthal - eh0) * ( cs_glob_elec_properties->th[itt+1]
548 - cs_glob_elec_properties->th[itt])
549 / (eh1 - eh0);
550 return temp;
551 }
552 }
553
554 assert(0); /* Should not arrive here */
555 return 0;
556 }
557
558 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
559
560 /*=============================================================================
561 * Public function definitions for Fortran API
562 *============================================================================*/
563
564 /*----------------------------------------------------------------------------*/
565
566 void
CS_PROCF(elini1,ELINI1)567 CS_PROCF (elini1, ELINI1) (cs_real_t *diftl0)
568 {
569 cs_electrical_model_specific_initialization(diftl0);
570 }
571
572 void
CS_PROCF(elflux,ELFLUX)573 CS_PROCF (elflux, ELFLUX) (int *iappel)
574 {
575 cs_elec_compute_fields(cs_glob_mesh,
576 *iappel);
577 }
578
579 void
CS_PROCF(elthht,ELTHHT)580 CS_PROCF (elthht, ELTHHT) (int *mode,
581 cs_real_t *ym,
582 cs_real_t *enthal,
583 cs_real_t *temp)
584 {
585 if (*mode == -1)
586 *enthal = _cs_elec_convert_t_to_h(ym, *temp);
587 else if (*mode == 1)
588 *temp = _cs_elec_convert_h_to_t(ym, *enthal);
589 else
590 bft_error(__FILE__, __LINE__, 0,
591 _("electric module:\n"
592 "bad value for mode (integer equal to -1 or 1 : %i here.\n"),
593 *mode);
594 }
595
596 void
CS_PROCF(ellecd,ELLECD)597 CS_PROCF (ellecd, ELLECD) (void)
598 {
599 cs_electrical_model_initialize();
600 cs_electrical_properties_read();
601 }
602
603 void
CS_PROCF(elphyv,ELPHYV)604 CS_PROCF (elphyv, ELPHYV) (void)
605 {
606 cs_elec_physical_properties(cs_glob_domain);
607 }
608
609 void
CS_PROCF(eltssc,ELTSSC)610 CS_PROCF (eltssc, ELTSSC) (const int *isca,
611 cs_real_t *smbrs)
612 {
613 const cs_mesh_t *mesh = cs_glob_mesh;
614 const cs_mesh_quantities_t *mesh_quantities = cs_glob_mesh_quantities;
615 const int keysca = cs_field_key_id("scalar_id");
616
617 for (int f_id = 0; f_id < cs_field_n_fields(); f_id++) {
618 cs_field_t *f = cs_field_by_id(f_id);
619 if (cs_field_get_key_int(f, keysca) == *isca)
620 cs_elec_source_terms(mesh, mesh_quantities, f->id, smbrs);
621 }
622 }
623
624 void
CS_PROCF(eltsvv,ELTSVV)625 CS_PROCF (eltsvv, ELTSVV) (const int *f_id,
626 cs_real_t *smbrv)
627 {
628 const cs_mesh_t *mesh = cs_glob_mesh;
629 const cs_mesh_quantities_t *mesh_quantities = cs_glob_mesh_quantities;
630
631 cs_elec_source_terms_v(mesh, mesh_quantities, *f_id, (cs_real_3_t *)smbrv);
632 }
633
634 void
CS_PROCF(eliniv,ELINIV)635 CS_PROCF (eliniv, ELINIV) (int *isuite)
636 {
637 cs_elec_fields_initialize(cs_glob_mesh, *isuite);
638 }
639
640 void
CS_PROCF(elreca,ELRECA)641 CS_PROCF (elreca, ELRECA) (cs_real_t *dt)
642 {
643 cs_elec_scaling_function(cs_glob_mesh, cs_glob_mesh_quantities, dt);
644 }
645
646 /*=============================================================================
647 * Public function definitions
648 *============================================================================*/
649
650 /*----------------------------------------------------------------------------
651 * Provide access to cs_elec_option
652 *----------------------------------------------------------------------------*/
653
654 cs_elec_option_t *
cs_get_glob_elec_option(void)655 cs_get_glob_elec_option(void)
656 {
657 return &_elec_option;
658 }
659
660 /*----------------------------------------------------------------------------
661 * Provide access to cs_glob_transformer
662 *----------------------------------------------------------------------------*/
663
664 cs_data_joule_effect_t *
cs_get_glob_transformer(void)665 cs_get_glob_transformer(void)
666 {
667 return _transformer;
668 }
669
670 /*----------------------------------------------------------------------------
671 * Initialize structures for electrical model
672 *----------------------------------------------------------------------------*/
673
674 void
cs_electrical_model_initialize(void)675 cs_electrical_model_initialize(void)
676 {
677 int ieljou = cs_glob_physical_model_flag[CS_JOULE_EFFECT];
678
679 if (ieljou >= 3)
680 BFT_MALLOC(_transformer, 1, cs_data_joule_effect_t);
681
682 _elec_option.ixkabe = 0;
683 _elec_option.ntdcla = 1;
684 _elec_option.irestrike = 0;
685 for (int i = 0; i < 3; i++)
686 _elec_option.restrike_point[i] = 0.;
687 _elec_option.izreca = NULL;
688 _elec_option.elcou = 0.;
689 _elec_option.ielcor = 0;
690 _elec_option.couimp = 0.;
691 _elec_option.puisim = 0.;
692 _elec_option.pot_diff = 0.;
693 _elec_option.coejou = 1.;
694 _elec_option.modrec = 1; /* standard model */
695 _elec_option.idreca = 3;
696 _elec_option.srrom = 0.;
697
698 for (int i = 0; i < 3; i++)
699 _elec_option.crit_reca[i] = 0.;
700 _elec_option.crit_reca[4] = 0.0002;
701
702 cs_glob_elec_option = &_elec_option;
703 cs_glob_elec_properties = &_elec_properties;
704 cs_glob_transformer = _transformer;
705
706 cs_fluid_properties_t *fluid_properties = cs_get_glob_fluid_properties();
707 fluid_properties->icp = 0;
708 fluid_properties->irovar = 1;
709 fluid_properties->ivivar = 1;
710 }
711
712 /*----------------------------------------------------------------------------
713 * Destroy structures for electrical model
714 *----------------------------------------------------------------------------*/
715
716 void
cs_electrical_model_finalize(void)717 cs_electrical_model_finalize(void)
718 {
719 int ieljou = cs_glob_physical_model_flag[CS_JOULE_EFFECT];
720 int ielarc = cs_glob_physical_model_flag[CS_ELECTRIC_ARCS];
721
722 if (ielarc > 0) {
723 BFT_FREE(_elec_properties.th);
724 BFT_FREE(_elec_properties.ehgaz);
725 BFT_FREE(_elec_properties.rhoel);
726 BFT_FREE(_elec_properties.cpel);
727 BFT_FREE(_elec_properties.sigel);
728 BFT_FREE(_elec_properties.visel);
729 BFT_FREE(_elec_properties.xlabel);
730 BFT_FREE(_elec_properties.xkabel);
731 }
732
733 if (ieljou >= 3) {
734 BFT_FREE(_transformer->tenspr);
735 BFT_FREE(_transformer->rnbs);
736 BFT_FREE(_transformer->zr);
737 BFT_FREE(_transformer->zi);
738 BFT_FREE(_transformer->ibrpr);
739 BFT_FREE(_transformer->ibrsec);
740 BFT_FREE(_transformer->tenspr);
741 BFT_FREE(_transformer->uroff);
742 BFT_FREE(_transformer->uioff);
743 BFT_FREE(_transformer);
744 }
745
746 BFT_FREE(_elec_option.izreca);
747 }
748
749 /*----------------------------------------------------------------------------
750 * Specific initialization for electric arc
751 *----------------------------------------------------------------------------*/
752
753 void
cs_electrical_model_specific_initialization(cs_real_t * diftl0)754 cs_electrical_model_specific_initialization(cs_real_t *diftl0)
755 {
756 cs_field_t *f = NULL;
757 int key_cal_opt_id = cs_field_key_id("var_cal_opt");
758 const int kvisls0 = cs_field_key_id("diffusivity_ref");
759 const int ksigmas = cs_field_key_id("turbulent_schmidt");
760 cs_var_cal_opt_t var_cal_opt;
761
762 /* specific initialization for field */
763 f = CS_F_(potr);
764 cs_field_get_key_struct(f, key_cal_opt_id, &var_cal_opt);
765 var_cal_opt.iconv = 0;
766 var_cal_opt.istat = 0;
767 var_cal_opt.idiff = 1;
768 var_cal_opt.idifft = 0;
769
770 cs_field_set_key_struct(f, key_cal_opt_id, &var_cal_opt);
771
772 int ieljou = cs_glob_physical_model_flag[CS_JOULE_EFFECT];
773 int ielarc = cs_glob_physical_model_flag[CS_ELECTRIC_ARCS];
774
775 if (ieljou == 2 || ieljou == 4) {
776 f = CS_F_(poti);
777 cs_field_get_key_struct(f, key_cal_opt_id, &var_cal_opt);
778 var_cal_opt.iconv = 0;
779 var_cal_opt.istat = 0;
780 var_cal_opt.idiff = 1;
781 var_cal_opt.idifft = 0;
782 cs_field_set_key_struct(f, key_cal_opt_id, &var_cal_opt);
783 }
784
785 if (ielarc > 1) {
786 cs_field_t *fp = cs_field_by_name_try("vec_potential");
787 cs_field_get_key_struct(fp, key_cal_opt_id, &var_cal_opt);
788 var_cal_opt.iconv = 0;
789 var_cal_opt.istat = 0;
790 var_cal_opt.idiff = 1;
791 var_cal_opt.idifft = 0;
792 cs_field_set_key_struct(fp, key_cal_opt_id, &var_cal_opt);
793 cs_field_set_key_double(fp, kvisls0, 1.0);
794 }
795
796 /* for all specific field */
797 f = CS_F_(h);
798 cs_field_get_key_struct(f, key_cal_opt_id, &var_cal_opt);
799 var_cal_opt.blencv = 1.;
800 cs_field_set_key_double(f, ksigmas, 0.7);
801 cs_field_set_key_struct(f, key_cal_opt_id, &var_cal_opt);
802
803 f = CS_F_(potr);
804 cs_field_get_key_struct(f, key_cal_opt_id, &var_cal_opt);
805 var_cal_opt.blencv = 1.;
806 cs_field_set_key_double(f, ksigmas, 0.7);
807 cs_field_set_key_struct(f, key_cal_opt_id, &var_cal_opt);
808
809 if (ieljou == 2 || ieljou == 4) {
810 f = CS_F_(poti);
811 cs_field_get_key_struct(f, key_cal_opt_id, &var_cal_opt);
812 var_cal_opt.blencv = 1.;
813 cs_field_set_key_double(f, ksigmas, 0.7);
814 cs_field_set_key_struct(f, key_cal_opt_id, &var_cal_opt);
815 }
816
817 if (ielarc > 1) {
818 cs_field_t *fp = cs_field_by_name_try("vec_potential");
819 cs_field_get_key_struct(fp, key_cal_opt_id, &var_cal_opt);
820 var_cal_opt.blencv = 1.;
821 cs_field_set_key_double(f, ksigmas, 0.7);
822 cs_field_set_key_struct(fp, key_cal_opt_id, &var_cal_opt);
823 }
824
825 if (cs_glob_elec_properties->ngaz > 1) {
826 for (int igaz = 0; igaz < cs_glob_elec_properties->ngaz - 1; igaz++) {
827 f = CS_FI_(ycoel, igaz);
828 cs_field_get_key_struct(f, key_cal_opt_id, &var_cal_opt);
829 var_cal_opt.blencv = 1.;
830 cs_field_set_key_double(f, ksigmas, 0.7);
831 cs_field_set_key_struct(f, key_cal_opt_id, &var_cal_opt);
832 }
833 }
834
835 CS_PROCF(uicpi1,UICPI1) (&(_elec_option.srrom), diftl0);
836 cs_gui_elec_model();
837 _elec_option.pot_diff = 1000.;//FIXME
838
839 _cs_electrical_model_verify();
840 }
841
842 /*----------------------------------------------------------------------------
843 * Read properties file
844 *----------------------------------------------------------------------------*/
845
846 void
cs_electrical_properties_read(void)847 cs_electrical_properties_read(void)
848 {
849 int ieljou = cs_glob_physical_model_flag[CS_JOULE_EFFECT];
850 int ielarc = cs_glob_physical_model_flag[CS_ELECTRIC_ARCS];
851
852 if (ielarc <= 0 && ieljou < 3)
853 return;
854
855 char str[LG_MAX];
856
857 if (ielarc > 0) {
858
859 /* read local file for electric properties if present,
860 default otherwise */
861
862 FILE *file = cs_base_open_properties_data_file("dp_ELE");
863
864 /* Position at the beginning of the file */
865 fseek(file, 0, SEEK_SET);
866
867 int nb_line_tot = 0;
868
869 int iesp = 0;
870 int it = 0;
871
872 while (fgets(str, LG_MAX, file) != NULL) {
873 nb_line_tot++;
874 if (nb_line_tot < 8)
875 continue;
876
877 /* read number of fluids and number of points */
878 if (nb_line_tot == 8)
879 sscanf(str, "%d %d",
880 &(_elec_properties.ngaz),
881 &(_elec_properties.npoint));
882
883 if (_elec_properties.ngaz <= 0)
884 bft_error(__FILE__, __LINE__, 0,
885 _("incorrect number of species \"%i\";\n"),
886 _elec_properties.ngaz);
887
888 cs_lnum_t size = cs_glob_elec_properties->ngaz
889 * cs_glob_elec_properties->npoint;
890
891 if (nb_line_tot == 8) {
892 BFT_MALLOC(_elec_properties.th,
893 cs_glob_elec_properties->npoint,
894 cs_real_t);
895 BFT_MALLOC(_elec_properties.ehgaz, size, cs_real_t);
896 BFT_MALLOC(_elec_properties.rhoel, size, cs_real_t);
897 BFT_MALLOC(_elec_properties.cpel, size, cs_real_t);
898 BFT_MALLOC(_elec_properties.sigel, size, cs_real_t);
899 BFT_MALLOC(_elec_properties.visel, size, cs_real_t);
900 BFT_MALLOC(_elec_properties.xlabel, size, cs_real_t);
901 BFT_MALLOC(_elec_properties.xkabel, size, cs_real_t);
902 }
903
904 if (nb_line_tot < 14)
905 continue;
906
907 if (nb_line_tot == 14)
908 sscanf(str, "%i", &(_elec_option.ixkabe));
909
910 if (cs_glob_elec_option->ixkabe < 0 || cs_glob_elec_option->ixkabe >= 3)
911 bft_error(__FILE__, __LINE__, 0,
912 _("incorrect choice for radiative model \"%i\";\n"),
913 cs_glob_elec_option->ixkabe < 0);
914
915 if (nb_line_tot < 22)
916 continue;
917
918 if (nb_line_tot >= 22) {
919 sscanf(str, "%lf %lf %lf %lf %lf %lf %lf %lf",
920 &(_elec_properties.th[it]),
921 &(_elec_properties.ehgaz[iesp * (cs_glob_elec_properties->npoint - 1) + it]),
922 &(_elec_properties.rhoel[iesp * (cs_glob_elec_properties->npoint - 1) + it]),
923 &(_elec_properties.cpel[iesp * (cs_glob_elec_properties->npoint - 1) + it]),
924 &(_elec_properties.sigel[iesp * (cs_glob_elec_properties->npoint - 1) + it]),
925 &(_elec_properties.visel[iesp * (cs_glob_elec_properties->npoint - 1) + it]),
926 &(_elec_properties.xlabel[iesp * (cs_glob_elec_properties->npoint - 1) + it]),
927 &(_elec_properties.xkabel[iesp * (cs_glob_elec_properties->npoint - 1) + it]));
928 it++;
929 if (it == cs_glob_elec_properties->npoint) {
930 iesp++;
931 it = 0;
932 }
933 }
934 }
935
936 fclose(file);
937 }
938
939 #if 0
940 for (int it = 0; it < cs_glob_elec_properties->npoint; it++)
941 bft_printf("read dp_ELE "
942 "%15.8E %15.8E %15.8E %15.8E "
943 "%15.8E %15.8E %15.8E %15.8E\n",
944 _elec_properties.th[it],
945 _elec_properties.ehgaz[it],
946 _elec_properties.rhoel[it],
947 _elec_properties.cpel[it],
948 _elec_properties.sigel[it],
949 _elec_properties.visel[it],
950 _elec_properties.xlabel[it],
951 _elec_properties.xkabel[it]);
952 #endif
953
954 if (ieljou >= 3) {
955
956 /* read local file for Joule effect if present,
957 default otherwise */
958
959 FILE *file = cs_base_open_properties_data_file("dp_transformers");
960
961 /* Position at the beginning of the file */
962 fseek(file, 0, SEEK_SET);
963
964 int nb_line_tot = 0;
965
966 int iesp = 0;
967 int it = 0;
968 while (fgets(str, LG_MAX, file) != NULL) {
969 nb_line_tot++;
970 if (nb_line_tot == 1)
971 sscanf(str, "%i", &(_transformer->ntfref));
972
973 if (nb_line_tot < 4)
974 continue;
975
976 if (nb_line_tot == 4) {
977 sscanf(str, "%i", &(_transformer->nbtrf));
978
979 BFT_MALLOC(_transformer->tenspr, cs_glob_transformer->nbtrf, cs_real_t);
980 BFT_MALLOC(_transformer->rnbs, cs_glob_transformer->nbtrf, cs_real_t);
981 BFT_MALLOC(_transformer->zr, cs_glob_transformer->nbtrf, cs_real_t);
982 BFT_MALLOC(_transformer->zi, cs_glob_transformer->nbtrf, cs_real_t);
983 BFT_MALLOC(_transformer->ibrpr, cs_glob_transformer->nbtrf, int);
984 BFT_MALLOC(_transformer->ibrsec, cs_glob_transformer->nbtrf, int);
985
986 // alloc for boundary conditions
987 BFT_MALLOC(_transformer->uroff, cs_glob_transformer->nbtrf, cs_real_t);
988 BFT_MALLOC(_transformer->uioff, cs_glob_transformer->nbtrf, cs_real_t);
989 }
990
991 if (nb_line_tot > 4 && nb_line_tot <= 4 + cs_glob_transformer->nbtrf * 6) {
992 it++;
993 if (it == 1)
994 continue;
995 if (it == 2)
996 sscanf(str, "%lf", &(_transformer->tenspr[iesp]));
997 if (it == 3)
998 sscanf(str, "%lf", &(_transformer->rnbs[iesp]));
999 if (it == 4)
1000 sscanf(str, "%lf %lf", &(_transformer->zr[iesp]), &(_transformer->zi[iesp]));
1001 if (it == 5)
1002 sscanf(str, "%i", &(_transformer->ibrpr[iesp]));
1003 if (it == 6) {
1004 sscanf(str, "%i", &(_transformer->ibrsec[iesp]));
1005 it = 0;
1006 iesp++;
1007 }
1008 }
1009
1010 if (nb_line_tot < 7 + cs_glob_transformer->nbtrf * 6)
1011 continue;
1012
1013 if (nb_line_tot == 7 + cs_glob_transformer->nbtrf * 6) {
1014 sscanf(str, "%i", &(_transformer->nbelec));
1015 BFT_MALLOC(_transformer->ielecc, cs_glob_transformer->nbelec, int);
1016 BFT_MALLOC(_transformer->ielect, cs_glob_transformer->nbelec, int);
1017 BFT_MALLOC(_transformer->ielecb, cs_glob_transformer->nbelec, int);
1018 iesp = 0;
1019 }
1020
1021 if (nb_line_tot > 7 + cs_glob_transformer->nbelec * 6) {
1022 sscanf(str, "%i %i %i",
1023 &(_transformer->ielecc[iesp]),
1024 &(_transformer->ielect[iesp]),
1025 &(_transformer->ielecb[iesp]));
1026 iesp++;
1027 }
1028 }
1029
1030 fclose(file);
1031 }
1032 }
1033
1034 /*----------------------------------------------------------------------------
1035 * compute physical properties
1036 *----------------------------------------------------------------------------*/
1037
1038 void
cs_elec_physical_properties(cs_domain_t * domain)1039 cs_elec_physical_properties(cs_domain_t *domain)
1040 {
1041 static long ipass = 0;
1042 int nt_cur = cs_glob_time_step->nt_cur;
1043 int isrrom = 0;
1044 const cs_lnum_t n_cells = domain->mesh->n_cells;
1045 const int keysca = cs_field_key_id("diffusivity_id");
1046 int diff_id = cs_field_get_key_int(CS_F_(potr), keysca);
1047 cs_field_t *c_prop = NULL;
1048 if (diff_id > -1)
1049 c_prop = cs_field_by_id(diff_id);
1050 ipass++;
1051
1052 const cs_data_elec_t *e_props = cs_glob_elec_properties; /* local name */
1053
1054 if (nt_cur > 1 && cs_glob_elec_option->srrom > 0.)
1055 isrrom = 1;
1056
1057 /* Joule effect (law must be specified by user) */
1058
1059 int ifcvsl = cs_field_get_key_int(CS_F_(h), keysca);
1060 cs_field_t *diff_th = NULL;
1061 if (ifcvsl >= 0)
1062 diff_th = cs_field_by_id(ifcvsl);
1063
1064 int ielarc = cs_glob_physical_model_flag[CS_ELECTRIC_ARCS];
1065
1066 /* Electric arc */
1067
1068 if (ielarc > 0) {
1069 if (ipass == 1)
1070 bft_printf("electric arc module: properties read on file.\n");
1071
1072 /* compute temperature from enthalpy */
1073 int ngaz = e_props->ngaz;
1074 int npt = e_props->npoint;
1075
1076 cs_real_t *ym, *yvol, *roesp, *visesp, *cpesp,
1077 *sigesp, *xlabes, *xkabes, *coef;
1078 BFT_MALLOC(ym, ngaz, cs_real_t);
1079 BFT_MALLOC(yvol, ngaz, cs_real_t);
1080 BFT_MALLOC(roesp, ngaz, cs_real_t);
1081 BFT_MALLOC(visesp, ngaz, cs_real_t);
1082 BFT_MALLOC(cpesp, ngaz, cs_real_t);
1083 BFT_MALLOC(sigesp, ngaz, cs_real_t);
1084 BFT_MALLOC(xlabes, ngaz, cs_real_t);
1085 BFT_MALLOC(xkabes, ngaz, cs_real_t);
1086 BFT_MALLOC(coef, ngaz * ngaz, cs_real_t);
1087
1088 int ifcsig = cs_field_get_key_int(CS_F_(potr), keysca);
1089
1090 if (ngaz == 1) {
1091 ym[0] = 1.;
1092
1093 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
1094 CS_F_(t)->val[iel] = _cs_elec_convert_h_to_t(ym, CS_F_(h)->val[iel]);
1095 }
1096 else {
1097
1098 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1099 ym[ngaz - 1] = 1.;
1100
1101 for (int ii = 0; ii < ngaz - 1; ii++) {
1102 ym[ii] = CS_FI_(ycoel, ii)->val[iel];
1103 ym[ngaz - 1] -= ym[ii];
1104 }
1105
1106 CS_F_(t)->val[iel] = _cs_elec_convert_h_to_t(ym, CS_F_(h)->val[iel]);
1107 }
1108 }
1109
1110 /* Map some fields */
1111
1112 cs_real_t *cpro_absco = NULL;
1113
1114 if (cs_glob_elec_option->ixkabe == 1) {
1115 if (CS_FI_(rad_cak, 0) != NULL)
1116 cpro_absco = CS_FI_(rad_cak, 0)->val;
1117 }
1118
1119 /* Interpolate properties */
1120
1121 # pragma omp parallel for
1122 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1123 // temperature
1124 cs_real_t tp = CS_F_(t)->val[iel];
1125
1126 // determine point
1127 int it = -1;
1128 if (tp <= e_props->th[0])
1129 it = 0;
1130 else if (tp >= e_props->th[npt - 1])
1131 it = npt - 1;
1132 else {
1133 for (int iiii = 0; iiii < npt - 1; iiii++)
1134 if (tp > e_props->th[iiii] &&
1135 tp <= e_props->th[iiii + 1]) {
1136 it = iiii;
1137 break;
1138 }
1139 }
1140 if (it == -1)
1141 bft_error(__FILE__, __LINE__, 0,
1142 _("electric module: properties read on file\n"
1143 "Warning: error in cs_elec_physical_properties\n"
1144 "Invalid reading with temperature : %f.\n"),
1145 tp);
1146
1147 /* mass frction */
1148 ym[ngaz - 1] = 1.;
1149
1150 for (int ii = 0; ii < ngaz - 1; ii++) {
1151 ym[ii] = CS_FI_(ycoel, ii)->val[iel];
1152 ym[ngaz - 1] -= ym[ii];
1153 }
1154
1155 /* density, viscosity, ... for each species */
1156 if (it == 0) {
1157 for (int ii = 0; ii < ngaz; ii++) {
1158 roesp[ii] = e_props->rhoel[ii * (npt - 1)];
1159 visesp[ii] = e_props->visel[ii * (npt - 1)];
1160 cpesp[ii] = e_props->cpel[ii * (npt - 1)];
1161 sigesp[ii] = e_props->sigel[ii * (npt - 1)];
1162 xlabes[ii] = e_props->xlabel[ii * (npt - 1)];
1163
1164 if (cs_glob_elec_option->ixkabe > 0)
1165 xkabes[ii] = e_props->xkabel[ii * (npt - 1)];
1166 }
1167 }
1168 else if (it == npt - 1) {
1169 for (int ii = 0; ii < ngaz; ii++) {
1170 roesp[ii] = e_props->rhoel[ii * (npt - 1) + npt - 1];
1171 visesp[ii] = e_props->visel[ii * (npt - 1) + npt - 1];
1172 cpesp[ii] = e_props->cpel[ii * (npt - 1) + npt - 1];
1173 sigesp[ii] = e_props->sigel[ii * (npt - 1) + npt - 1];
1174 xlabes[ii] = e_props->xlabel[ii * (npt - 1) + npt - 1];
1175
1176 if (cs_glob_elec_option->ixkabe > 0)
1177 xkabes[ii] = e_props->xkabel[ii * (npt - 1) + npt - 1];
1178 }
1179 }
1180 else {
1181 double delt = e_props->th[it + 1] - e_props->th[it];
1182
1183 for (int ii = 0; ii < ngaz; ii++) {
1184 double alpro = (e_props->rhoel[ii * (npt - 1) + it + 1] -
1185 e_props->rhoel[ii * (npt - 1) + it]) / delt;
1186 roesp[ii] = e_props->rhoel[ii * (npt - 1) + it]
1187 + alpro * (tp -e_props->th[it]);
1188
1189 double alpvis = (e_props->visel[ii * (npt - 1) + it + 1]
1190 - e_props->visel[ii * (npt - 1) + it]) / delt;
1191 visesp[ii] = e_props->visel[ii * (npt - 1) + it]
1192 + alpvis * (tp -e_props->th[it]);
1193
1194 double alpcp = (e_props->cpel[ii * (npt - 1) + it + 1] -
1195 e_props->cpel[ii * (npt - 1) + it]) / delt;
1196 cpesp[ii] = e_props->cpel[ii * (npt - 1) + it] +
1197 alpcp * (tp -e_props->th[it]);
1198
1199 double alpsig = (e_props->sigel[ii * (npt - 1) + it + 1] -
1200 e_props->sigel[ii * (npt - 1) + it]) / delt;
1201 sigesp[ii] = e_props->sigel[ii * (npt - 1) + it] +
1202 alpsig * (tp -e_props->th[it]);
1203
1204 double alplab = (e_props->xlabel[ii * (npt - 1) + it + 1] -
1205 e_props->xlabel[ii * (npt - 1) + it]) / delt;
1206 xlabes[ii] = e_props->xlabel[ii * (npt - 1) + it] +
1207 alplab * (tp -e_props->th[it]);
1208
1209 if (cs_glob_elec_option->ixkabe > 0) {
1210 double alpkab = (e_props->xkabel[ii * (npt - 1) + it + 1] -
1211 e_props->xkabel[ii * (npt - 1) + it]) / delt;
1212 xkabes[ii] = e_props->xkabel[ii * (npt - 1) + it] +
1213 alpkab * (tp -e_props->th[it]);
1214 }
1215 }
1216 }
1217
1218 /* compute density */
1219 double rhonp1 = 0.;
1220
1221 for (int ii = 0; ii < ngaz; ii++)
1222 rhonp1 += ym[ii] / roesp[ii];
1223
1224 rhonp1 = 1. / rhonp1;
1225
1226 if (isrrom == 1)
1227 CS_F_(rho)->val[iel] = CS_F_(rho)->val[iel] * cs_glob_elec_option->srrom
1228 + (1. - cs_glob_elec_option->srrom) * rhonp1;
1229 else
1230 CS_F_(rho)->val[iel] = rhonp1;
1231
1232 for (int ii = 0; ii < ngaz; ii++) {
1233 yvol[ii] = ym[ii] * roesp[ii] / CS_F_(rho)->val[iel];
1234 if (yvol[ii] <= 0.)
1235 yvol[ii] = cs_math_epzero * cs_math_epzero;
1236 }
1237
1238 /* compute molecular viscosity : kg/(m s) */
1239 for (int iesp1 = 0; iesp1 < ngaz; iesp1++) {
1240 for (int iesp2 = 0; iesp2 < ngaz; iesp2++) {
1241 coef[iesp1 * (ngaz - 1) + iesp2]
1242 = 1. + sqrt(visesp[iesp1] / visesp[iesp2])
1243 * sqrt(sqrt(roesp[iesp2] / roesp[iesp1]));
1244 coef[iesp1 * (ngaz - 1) + iesp2] *= coef[iesp1 * (ngaz - 1) + iesp2];
1245 coef[iesp1 * (ngaz - 1) + iesp2] /= (sqrt(1. + roesp[iesp1]
1246 / roesp[iesp2]) * sqrt(8.));
1247 }
1248 }
1249
1250 CS_F_(mu)->val[iel] = 0.;
1251
1252 for (int iesp1 = 0; iesp1 < ngaz; iesp1++) {
1253 if (yvol[iesp1] > 1e-30) {
1254 double somphi = 0.;
1255 for (int iesp2 = 0; iesp2 < ngaz; iesp2++) {
1256 if (iesp1 != iesp2)
1257 somphi += coef[iesp1 * (ngaz - 1) + iesp2]
1258 * yvol[iesp2] / yvol[iesp1];
1259 }
1260
1261 CS_F_(mu)->val[iel] += visesp[iesp1] / (1. + somphi);
1262 }
1263 }
1264
1265 /* compute specific heat : J/(kg degres) */
1266 if (cs_glob_fluid_properties->icp > 0) {
1267 CS_F_(cp)->val[iel] = 0.;
1268 for (int iesp1 = 0; iesp1 < ngaz; iesp1++)
1269 CS_F_(cp)->val[iel] += ym[iesp1] * cpesp[iesp1];
1270 }
1271
1272 /* compute Lambda/Cp : kg/(m s) */
1273 if (diff_th != NULL) {
1274
1275 for (int iesp1 = 0; iesp1 < ngaz; iesp1++) {
1276 for (int iesp2 = 0; iesp2 < ngaz; iesp2++) {
1277 coef[iesp1 * (ngaz - 1) + iesp2]
1278 = 1. + sqrt(xlabes[iesp1] / xlabes[iesp2])
1279 * sqrt(sqrt(roesp[iesp2] / roesp[iesp1]));
1280 coef[iesp1 * (ngaz - 1) + iesp2]
1281 *= coef[iesp1 * (ngaz - 1) + iesp2];
1282 coef[iesp1 * (ngaz - 1) + iesp2]
1283 /= (sqrt(1. + roesp[iesp1] / roesp[iesp2]) * sqrt(8.));
1284 }
1285 }
1286 /* Lambda */
1287 diff_th->val[iel] = 0.;
1288
1289 for (int iesp1 = 0; iesp1 < ngaz; iesp1++) {
1290 if (yvol[iesp1] > 1e-30) {
1291 double somphi = 0.;
1292 for (int iesp2 = 0; iesp2 < ngaz; iesp2++) {
1293 if (iesp1 != iesp2)
1294 somphi += coef[iesp1 * (ngaz - 1) + iesp2]
1295 * yvol[iesp2] / yvol[iesp1];
1296 }
1297
1298 diff_th->val[iel] += xlabes[iesp1] / (1. + 1.065 * somphi);
1299 }
1300 }
1301
1302 /* Lambda/Cp */
1303 if (cs_glob_fluid_properties->icp <= 0)
1304 diff_th->val[iel] /= cs_glob_fluid_properties->cp0;
1305 else
1306 diff_th->val[iel] /= CS_F_(cp)->val[iel];
1307 }
1308
1309 /* compute electric conductivity: S/m */
1310 if (ifcsig >= 0) {
1311 c_prop->val[iel] = 0.;
1312 double val = 0.;
1313
1314 for (int iesp1 = 0; iesp1 < ngaz; iesp1++)
1315 val += yvol[iesp1] / sigesp[iesp1];
1316
1317 c_prop->val[iel] = 1. / val;
1318 }
1319
1320 /* compute radiative transfer : W/m3 */
1321 if (cs_glob_elec_option->ixkabe == 1) {
1322 if (cpro_absco != NULL) { /* May be NULL if no active radiation model */
1323 double val = 0.;
1324 for (int iesp1 = 0; iesp1 < ngaz; iesp1++)
1325 val += yvol[iesp1] * xkabes[iesp1];
1326 cpro_absco[iel] = val;
1327 }
1328 }
1329 else if (cs_glob_elec_option->ixkabe == 2) {
1330 CS_F_(radsc)->val[iel] = 0.;
1331 double val = 0.;
1332
1333 for (int iesp1 = 0; iesp1 < ngaz; iesp1++)
1334 val += yvol[iesp1] * xkabes[iesp1];
1335
1336 CS_F_(radsc)->val[iel] = val;
1337 }
1338
1339 /* diffusivity for other properties
1340 * nothing to do
1341 * no other properties in this case */
1342
1343 } /* End of loop on cells */
1344
1345 BFT_FREE(ym);
1346 BFT_FREE(yvol);
1347 BFT_FREE(roesp);
1348 BFT_FREE(visesp);
1349 BFT_FREE(cpesp);
1350 BFT_FREE(sigesp);
1351 BFT_FREE(xlabes);
1352 BFT_FREE(xkabes);
1353 BFT_FREE(coef);
1354 }
1355
1356 /* now user properties (for joule effect particulary) */
1357 cs_user_physical_properties(domain);
1358 }
1359
1360 /*----------------------------------------------------------------------------
1361 * compute specific electric arc fields
1362 *----------------------------------------------------------------------------*/
1363
1364 void
cs_elec_compute_fields(const cs_mesh_t * mesh,int call_id)1365 cs_elec_compute_fields(const cs_mesh_t *mesh,
1366 int call_id)
1367 {
1368 cs_lnum_t n_cells = mesh->n_cells;
1369 cs_lnum_t n_cells_ext = mesh->n_cells_with_ghosts;
1370 const int keysca = cs_field_key_id("diffusivity_id");
1371
1372 int ieljou = cs_glob_physical_model_flag[CS_JOULE_EFFECT];
1373 int ielarc = cs_glob_physical_model_flag[CS_ELECTRIC_ARCS];
1374
1375 bool log_active = cs_log_default_is_active();
1376
1377 /* Reconstructed value */
1378 cs_real_3_t *grad;
1379 BFT_MALLOC(grad, n_cells_ext, cs_real_3_t);
1380
1381 /* ----------------------------------------------------- */
1382 /* first call : J, E => J.E */
1383 /* ----------------------------------------------------- */
1384
1385 if (call_id == 1) {
1386 /* compute grad(potR) */
1387
1388 /* Get the calculation option from the field */
1389 cs_real_3_t *cpro_elefl = (cs_real_3_t *)(CS_F_(elefl)->val);
1390
1391 cs_field_gradient_scalar(CS_F_(potr),
1392 false, /* use_previous_t */
1393 1, /* inc */
1394 true, /* recompute_cocg */
1395 grad);
1396
1397 /* compute electric field E = - grad (potR) */
1398 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1399 cpro_elefl[iel][0] = grad[iel][0];
1400 cpro_elefl[iel][1] = grad[iel][1];
1401 cpro_elefl[iel][2] = grad[iel][2];
1402 }
1403
1404 /* compute current density j = sig E */
1405 int diff_id = cs_field_get_key_int(CS_F_(potr), keysca);
1406 cs_field_t *c_prop = NULL;
1407 if (diff_id > -1)
1408 c_prop = cs_field_by_id(diff_id);
1409
1410 if (ieljou > 0 || ielarc > 0) {
1411 cs_real_3_t *cpro_curre = (cs_real_3_t *)(CS_F_(curre)->val);
1412 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1413 cpro_curre[iel][0] = -c_prop->val[iel] * grad[iel][0];
1414 cpro_curre[iel][1] = -c_prop->val[iel] * grad[iel][1];
1415 cpro_curre[iel][2] = -c_prop->val[iel] * grad[iel][2];
1416 }
1417 }
1418
1419 /* compute joule effect : j . E */
1420 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1421 CS_F_(joulp)->val[iel] = c_prop->val[iel] *
1422 (grad[iel][0] * grad[iel][0] +
1423 grad[iel][1] * grad[iel][1] +
1424 grad[iel][2] * grad[iel][2]);
1425 }
1426
1427 /* compute min max for E and J */
1428 if (log_active) {
1429 bft_printf("-----------------------------------------\n"
1430 " Variable Minimum Maximum\n"
1431 "-----------------------------------------\n");
1432
1433 /* Grad PotR = -E */
1434 double vrmin[3], vrmax[3];
1435
1436 for (int i = 0; i < 3; i++) {
1437 vrmin[i] = HUGE_VAL;
1438 vrmax[i] = -HUGE_VAL;
1439 }
1440
1441 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1442 for (int i = 0; i < 3; i++) {
1443 vrmin[i] = CS_MIN(vrmin[i], grad[iel][i]);
1444 vrmax[i] = CS_MAX(vrmax[i], grad[iel][i]);
1445 }
1446 }
1447
1448 cs_parall_min(3, CS_DOUBLE, vrmin);
1449 cs_parall_max(3, CS_DOUBLE, vrmax);
1450
1451 for (int i = 0; i < 3; i++) {
1452 bft_printf("v Gr_PotR%s %12.5e %12.5e\n",
1453 cs_glob_field_comp_name_3[i],
1454 vrmin[i], vrmax[i]);
1455 }
1456
1457 /* current real */
1458 for (int i = 0; i < 3; i++) {
1459 vrmin[i] = HUGE_VAL;
1460 vrmax[i] = -HUGE_VAL;
1461 }
1462
1463 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1464 for (int i = 0; i < 3; i++) {
1465 vrmin[i] = CS_MIN(vrmin[i], -c_prop->val[iel] * grad[iel][i]);
1466 vrmax[i] = CS_MAX(vrmax[i], -c_prop->val[iel] * grad[iel][i]);
1467 }
1468 }
1469
1470 cs_parall_min(3, CS_DOUBLE, vrmin);
1471 cs_parall_max(3, CS_DOUBLE, vrmax);
1472
1473 for (int i = 0; i < 3; i++) {
1474 bft_printf("v Cour_Re%s %12.5E %12.5E\n",
1475 cs_glob_field_comp_name_3[i],
1476 vrmin[i], vrmax[i]);
1477 }
1478 bft_printf("-----------------------------------------\n");
1479 }
1480
1481 if (ieljou == 2 || ieljou == 4) {
1482 /* compute grad(potI) */
1483
1484 cs_field_gradient_scalar(CS_F_(poti),
1485 false, /* use_previous_t */
1486 1, /* inc */
1487 true, /* recompute_cocg */
1488 grad);
1489
1490 /* compute electric field E = - grad (potI) */
1491
1492 /* compute current density j = sig E */
1493 int diff_id_i = cs_field_get_key_int(CS_F_(poti), keysca);
1494 cs_field_t *c_propi = NULL;
1495 if (diff_id_i > -1)
1496 c_propi = cs_field_by_id(diff_id_i);
1497
1498 if (ieljou == 4) {
1499 cs_real_3_t *cpro_curim = (cs_real_3_t *)(CS_F_(curim)->val);
1500 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1501 cpro_curim[iel][0] = -c_propi->val[iel] * grad[iel][0];
1502 cpro_curim[iel][1] = -c_propi->val[iel] * grad[iel][1];
1503 cpro_curim[iel][2] = -c_propi->val[iel] * grad[iel][2];
1504 }
1505 }
1506
1507 /* compute joule effect : j . E */
1508 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1509 CS_F_(joulp)->val[iel] += c_propi->val[iel]
1510 * cs_math_3_square_norm(grad[iel]);
1511 }
1512
1513 /* compute min max for E and J */
1514 if (log_active) {
1515
1516 double vrmin[3], vrmax[3];
1517
1518 /* Grad PotR = -Ei */
1519
1520 for (int i = 0; i < 3; i++) {
1521 vrmin[i] = HUGE_VAL;
1522 vrmax[i] = -HUGE_VAL;
1523 }
1524
1525 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1526 for (int i = 0; i < 3; i++) {
1527 vrmin[i] = CS_MIN(vrmin[i], grad[iel][0]);
1528 vrmax[i] = CS_MAX(vrmax[i], grad[iel][0]);
1529 }
1530 }
1531
1532 cs_parall_min(3, CS_DOUBLE, vrmin);
1533 cs_parall_max(3, CS_DOUBLE, vrmax);
1534
1535 for (int i = 0; i < 3; i++) {
1536 bft_printf("v Gr_PotI%s %12.5E %12.5E\n",
1537 cs_glob_field_comp_name_3[i],
1538 vrmin[i], vrmax[i]);
1539 }
1540
1541 /* Imaginary current */
1542
1543 for (int i = 0; i < 3; i++) {
1544 vrmin[i] = HUGE_VAL;
1545 vrmax[i] = -HUGE_VAL;
1546 }
1547
1548 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1549 for (int i = 0; i < 3; i++) {
1550 vrmin[i] = CS_MIN(vrmin[i], -c_propi->val[iel] * grad[iel][i]);
1551 vrmax[i] = CS_MAX(vrmax[i], -c_propi->val[iel] * grad[iel][i]);
1552 }
1553 }
1554
1555 cs_parall_min(3, CS_DOUBLE, vrmin);
1556 cs_parall_max(3, CS_DOUBLE, vrmax);
1557
1558 for (int i = 0; i < 3; i++) {
1559 bft_printf("v Cour_Im%s %12.5E %12.5E\n",
1560 cs_glob_field_comp_name_3[i],
1561 vrmin[i], vrmax[i]);
1562 }
1563 }
1564 }
1565 }
1566
1567 /* ----------------------------------------------------- */
1568 /* second call : A, B, JXB */
1569 /* ----------------------------------------------------- */
1570
1571 else if (call_id == 2) {
1572
1573 cs_real_3_t *cpro_magfl = (cs_real_3_t *)(CS_F_(magfl)->val);
1574
1575 if (ielarc == 2) {
1576 /* compute magnetic field component B */
1577 cs_field_t *fp = cs_field_by_name_try("vec_potential");
1578
1579 cs_real_33_t *gradv = NULL;
1580 BFT_MALLOC(gradv, n_cells_ext, cs_real_33_t);
1581
1582 cs_field_gradient_vector(fp,
1583 false, /* use_previous_t */
1584 1, /* inc */
1585 gradv);
1586
1587 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1588 cpro_magfl[iel][0] = -gradv[iel][1][2]+gradv[iel][2][1];
1589 cpro_magfl[iel][1] = gradv[iel][0][2]-gradv[iel][2][0];
1590 cpro_magfl[iel][2] = -gradv[iel][0][1]+gradv[iel][1][0];
1591 }
1592
1593 BFT_FREE(gradv);
1594 }
1595 else if (ielarc == 1)
1596 bft_error(__FILE__, __LINE__, 0,
1597 _("Error electric arc with ampere theorem not available\n"));
1598
1599 /* compute laplace effect j x B */
1600 cs_real_3_t *cpro_laplf = (cs_real_3_t *)(CS_F_(laplf)->val);
1601 cs_real_3_t *cpro_curre = (cs_real_3_t *)(CS_F_(curre)->val);
1602 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1603 cpro_laplf[iel][0] = cpro_curre[iel][1] * cpro_magfl[iel][2]
1604 - cpro_curre[iel][2] * cpro_magfl[iel][1];
1605 cpro_laplf[iel][1] = cpro_curre[iel][2] * cpro_magfl[iel][0]
1606 - cpro_curre[iel][0] * cpro_magfl[iel][2];
1607 cpro_laplf[iel][2] = cpro_curre[iel][0] * cpro_magfl[iel][1]
1608 - cpro_curre[iel][1] * cpro_magfl[iel][0];
1609 }
1610
1611 /* compute min max for B */
1612 if (ielarc > 1 && log_active) {
1613 /* Grad PotR = -E */
1614 double vrmin[3], vrmax[3];
1615
1616 for (int i = 0; i < 3; i++) {
1617 vrmin[i] = HUGE_VAL;
1618 vrmax[i] = -HUGE_VAL;
1619 }
1620
1621 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1622 for (int i = 0; i < 3; i++) {
1623 vrmin[i] = CS_MIN(vrmin[i], cpro_magfl[iel][i]);
1624 vrmax[i] = CS_MAX(vrmax[i], cpro_magfl[iel][i]);
1625 }
1626 }
1627
1628 cs_parall_min(3, CS_DOUBLE, vrmin);
1629 cs_parall_max(3, CS_DOUBLE, vrmax);
1630
1631 for (int i = 0; i < 3; i++) {
1632 bft_printf("v Magnetic_field%s %12.5E %12.5E\n",
1633 cs_glob_field_comp_name_3[i], vrmin[i], vrmax[i]);
1634 }
1635 }
1636 }
1637
1638 /* Free memory */
1639 BFT_FREE(grad);
1640 }
1641
1642 /*----------------------------------------------------------------------------
1643 * compute source terms for energy
1644 *----------------------------------------------------------------------------*/
1645
1646 void
cs_elec_source_terms(const cs_mesh_t * mesh,const cs_mesh_quantities_t * mesh_quantities,int f_id,cs_real_t * smbrs)1647 cs_elec_source_terms(const cs_mesh_t *mesh,
1648 const cs_mesh_quantities_t *mesh_quantities,
1649 int f_id,
1650 cs_real_t *smbrs)
1651 {
1652 const cs_field_t *f = cs_field_by_id(f_id);
1653 const char *name = f->name;
1654 cs_lnum_t n_cells = mesh->n_cells;
1655 cs_lnum_t n_cells_ext = mesh->n_cells_with_ghosts;
1656 const cs_real_t *volume = mesh_quantities->cell_vol;
1657
1658 int key_cal_opt_id = cs_field_key_id("var_cal_opt");
1659 cs_var_cal_opt_t var_cal_opt;
1660 cs_field_get_key_struct(f, key_cal_opt_id, &var_cal_opt);
1661
1662 int ielarc = cs_glob_physical_model_flag[CS_ELECTRIC_ARCS];
1663
1664 cs_real_t *w1;
1665 BFT_MALLOC(w1, n_cells_ext, cs_real_t);
1666
1667 /* enthalpy source term */
1668 if (strcmp(name, "enthalpy") == 0) {
1669 if (var_cal_opt.verbosity > 0)
1670 bft_printf("compute source terms for variable : %s\n", name);
1671
1672 if (cs_glob_time_step->nt_cur > 2) {
1673 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
1674 w1[iel] = CS_F_(joulp)->val[iel] * volume[iel];
1675
1676 if (ielarc >= 1)
1677 if (cs_glob_elec_option->ixkabe == 2)
1678 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
1679 w1[iel] -= CS_F_(radsc)->val[iel] * volume[iel];
1680
1681 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
1682 smbrs[iel] += w1[iel];
1683
1684 if (var_cal_opt.verbosity > 0) {
1685 double valmin = w1[0];
1686 double valmax = w1[0];
1687
1688 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1689 valmin = CS_MIN(valmin, w1[iel]);
1690 valmax = CS_MAX(valmax, w1[iel]);
1691 }
1692
1693 cs_parall_min(1, CS_DOUBLE, &valmin);
1694 cs_parall_max(1, CS_DOUBLE, &valmax);
1695
1696 bft_printf(" source terms for H min= %14.5E, max= %14.5E\n",
1697 valmin, valmax);
1698 }
1699 }
1700 }
1701
1702 BFT_FREE(w1);
1703 }
1704
1705 /*----------------------------------------------------------------------------
1706 * compute source terms for vector potential
1707 *----------------------------------------------------------------------------*/
1708
1709 void
cs_elec_source_terms_v(const cs_mesh_t * mesh,const cs_mesh_quantities_t * mesh_quantities,int f_id,cs_real_3_t * smbrv)1710 cs_elec_source_terms_v(const cs_mesh_t *mesh,
1711 const cs_mesh_quantities_t *mesh_quantities,
1712 int f_id,
1713 cs_real_3_t *smbrv)
1714 {
1715 const cs_field_t *f = cs_field_by_id(f_id);
1716 cs_lnum_t n_cells = mesh->n_cells;
1717 const cs_real_t *volume = mesh_quantities->cell_vol;
1718
1719 int key_cal_opt_id = cs_field_key_id("var_cal_opt");
1720 cs_var_cal_opt_t var_cal_opt;
1721 cs_field_get_key_struct(f, key_cal_opt_id, &var_cal_opt);
1722
1723 int ielarc = cs_glob_physical_model_flag[CS_ELECTRIC_ARCS];
1724
1725 /* source term for potential vector */
1726
1727 if (ielarc >= 2 && f_id == (CS_F_(potva)->id)) {
1728 cs_real_3_t *cpro_curre = (cs_real_3_t *)(CS_F_(curre)->val);
1729
1730 if (var_cal_opt.verbosity > 0)
1731 bft_printf("compute source terms for variable: %s\n", f->name);
1732
1733 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
1734 for (int isou = 0; isou < 3; isou++)
1735 smbrv[iel][isou] += cs_elec_permvi * cpro_curre[iel][isou] * volume[iel];
1736 }
1737 }
1738
1739 /*----------------------------------------------------------------------------
1740 * add variables fields
1741 *----------------------------------------------------------------------------*/
1742
1743 void
cs_elec_add_variable_fields(void)1744 cs_elec_add_variable_fields(void)
1745 {
1746 cs_field_t *f;
1747 int dim1 = 1;
1748 int dim3 = 3;
1749 double grand = 1.e12;
1750
1751 const int kscmin = cs_field_key_id("min_scalar_clipping");
1752 const int kscmax = cs_field_key_id("max_scalar_clipping");
1753 const int kivisl = cs_field_key_id("diffusivity_id");
1754
1755 const cs_data_elec_t *e_props = cs_glob_elec_properties; /* local name */
1756
1757 int ielarc = cs_glob_physical_model_flag[CS_ELECTRIC_ARCS];
1758 int ieljou = cs_glob_physical_model_flag[CS_JOULE_EFFECT];
1759
1760 {
1761 int f_id = cs_variable_field_create("enthalpy", "Enthalpy",
1762 CS_MESH_LOCATION_CELLS, dim1);
1763 f = cs_field_by_id(f_id);
1764 cs_field_set_key_double(f, kscmin, -grand);
1765 cs_field_set_key_int(f, kivisl, 0);
1766 int isca = cs_add_model_field_indexes(f->id);
1767
1768 /* set thermal model */
1769 cs_thermal_model_t *thermal = cs_get_glob_thermal_model();
1770 thermal->itherm = CS_THERMAL_MODEL_ENTHALPY;
1771 thermal->iscalt = isca;
1772 }
1773
1774 {
1775 int f_id = cs_variable_field_create("elec_pot_r", "POT_EL_R",
1776 CS_MESH_LOCATION_CELLS, dim1);
1777 f = cs_field_by_id(f_id);
1778 cs_field_set_key_double(f, kscmin, -grand);
1779 cs_field_set_key_double(f, kscmax, grand);
1780 cs_field_set_key_int(f, kivisl, 0);
1781 cs_add_model_field_indexes(f->id);
1782 }
1783
1784 if (ieljou == 2 || ieljou == 4) {
1785 int f_id = cs_variable_field_create("elec_pot_i", "POT_EL_I",
1786 CS_MESH_LOCATION_CELLS, dim1);
1787 f = cs_field_by_id(f_id);
1788 cs_field_set_key_double(f, kscmin, -grand);
1789 cs_field_set_key_double(f, kscmax, grand);
1790 cs_field_set_key_int(f, kivisl, 0);
1791 cs_add_model_field_indexes(f->id);
1792 }
1793
1794 if (ielarc > 1) {
1795 int f_id = cs_variable_field_create("vec_potential", "POT_VEC",
1796 CS_MESH_LOCATION_CELLS, dim3);
1797 f = cs_field_by_id(f_id);
1798 //cs_field_set_key_double(f, kscmin, -grand);
1799 //cs_field_set_key_double(f, kscmax, grand);
1800 cs_field_set_key_int(f, kivisl, -1);
1801 cs_add_model_field_indexes(f->id);
1802 }
1803
1804 if (e_props->ngaz > 1) {
1805 for (int igaz = 0; igaz < e_props->ngaz - 1; igaz++) {
1806 char *name = NULL;
1807 char *label = NULL;
1808 char *suf = NULL;
1809 BFT_MALLOC(name, strlen("esl_fraction_") + 2 + 1, char);
1810 BFT_MALLOC(label, strlen("YM_ESL") + 2 + 1, char);
1811 BFT_MALLOC(suf, 3, char);
1812 strcpy(name, "esl_fraction_");
1813 strcpy(label, "YM_ESL");
1814 sprintf(suf, "%02d", igaz + 1);
1815 strcat(name, suf);
1816 strcat(label, suf);
1817
1818 int f_id = cs_variable_field_create(name, label,
1819 CS_MESH_LOCATION_CELLS, dim1);
1820 f = cs_field_by_id(f_id);
1821
1822 cs_field_set_key_double(f, kscmin, 0.);
1823 cs_field_set_key_double(f, kscmax, 1.);
1824 cs_field_set_key_int(f, kivisl, 0);
1825 cs_add_model_field_indexes(f->id);
1826 BFT_FREE(name);
1827 BFT_FREE(label);
1828 BFT_FREE(suf);
1829 }
1830 }
1831
1832 _field_pointer_map_electric_arcs(e_props->ngaz);
1833 }
1834
1835 /*----------------------------------------------------------------------------
1836 * add properties fields
1837 *----------------------------------------------------------------------------*/
1838
1839 void
cs_elec_add_property_fields(void)1840 cs_elec_add_property_fields(void)
1841 {
1842 cs_field_t *f;
1843 int field_type = CS_FIELD_INTENSIVE | CS_FIELD_PROPERTY;
1844 bool has_previous = false;
1845 const int klbl = cs_field_key_id("label");
1846 const int keyvis = cs_field_key_id("post_vis");
1847 const int keylog = cs_field_key_id("log");
1848 const int post_flag = CS_POST_ON_LOCATION | CS_POST_MONITOR;
1849
1850 int ieljou = cs_glob_physical_model_flag[CS_JOULE_EFFECT];
1851
1852 {
1853 f = cs_field_create("temperature",
1854 field_type,
1855 CS_MESH_LOCATION_CELLS,
1856 1, /* dim */
1857 has_previous);
1858 cs_field_set_key_int(f, keyvis, post_flag);
1859 cs_field_set_key_int(f, keylog, 1);
1860 cs_field_set_key_str(f, klbl, "Temperature");
1861 }
1862
1863 {
1864 f = cs_field_create("joule_power",
1865 field_type,
1866 CS_MESH_LOCATION_CELLS,
1867 1, /* dim */
1868 has_previous);
1869 cs_field_set_key_int(f, keyvis, post_flag);
1870 cs_field_set_key_int(f, keylog, 1);
1871 cs_field_set_key_str(f, klbl, "PowJoul");
1872 }
1873
1874 {
1875 f = cs_field_create("current_re",
1876 field_type,
1877 CS_MESH_LOCATION_CELLS,
1878 3, /* dim */
1879 has_previous);
1880 cs_field_set_key_int(f, keyvis, post_flag);
1881 cs_field_set_key_int(f, keylog, 1);
1882 cs_field_set_key_str(f, klbl, "Current_Real");
1883 }
1884
1885 {
1886 f = cs_field_create("electric_field",
1887 field_type,
1888 CS_MESH_LOCATION_CELLS,
1889 3, /* dim */
1890 has_previous);
1891 cs_field_set_key_int(f, keyvis, post_flag);
1892 cs_field_set_key_int(f, keylog, 1);
1893 cs_field_set_key_str(f, klbl, "Elec_Field");
1894 }
1895
1896 /* specific for joule effect */
1897 if (ieljou == 2 || ieljou == 4) {
1898 f = cs_field_create("current_im",
1899 field_type,
1900 CS_MESH_LOCATION_CELLS,
1901 3, /* dim */
1902 has_previous);
1903 cs_field_set_key_int(f, keyvis, post_flag);
1904 cs_field_set_key_int(f, keylog, 1);
1905 cs_field_set_key_str(f, klbl, "Current_Imag");
1906 }
1907
1908 /* specific for electric arcs */
1909 {
1910 f = cs_field_create("laplace_force",
1911 field_type,
1912 CS_MESH_LOCATION_CELLS,
1913 3, /* dim */
1914 has_previous);
1915 cs_field_set_key_int(f, keyvis, post_flag);
1916 cs_field_set_key_int(f, keylog, 1);
1917 cs_field_set_key_str(f, klbl, "For_Lap");
1918
1919 f = cs_field_create("magnetic_field",
1920 field_type,
1921 CS_MESH_LOCATION_CELLS,
1922 3, /* dim */
1923 has_previous);
1924 cs_field_set_key_int(f, keyvis, post_flag);
1925 cs_field_set_key_int(f, keylog, 1);
1926 cs_field_set_key_str(f, klbl, "Mag_Field");
1927 }
1928
1929 if (cs_glob_elec_option->ixkabe == 1) {
1930 f = cs_field_create("absorption_coeff",
1931 field_type,
1932 CS_MESH_LOCATION_CELLS,
1933 1, /* dim */
1934 has_previous);
1935 cs_field_set_key_int(f, keyvis, post_flag);
1936 cs_field_set_key_int(f, keylog, 1);
1937 cs_field_set_key_str(f, klbl, "Coef_Abso");
1938 }
1939 else if (cs_glob_elec_option->ixkabe == 2) {
1940 f = cs_field_create("radiation_source",
1941 field_type,
1942 CS_MESH_LOCATION_CELLS,
1943 1, /* dim */
1944 has_previous);
1945 cs_field_set_key_int(f, keyvis, post_flag);
1946 cs_field_set_key_int(f, keylog, 1);
1947 cs_field_set_key_str(f, klbl, "ST_radia");
1948 }
1949
1950
1951 _field_pointer_properties_map_electric_arcs();
1952 }
1953
1954 /*----------------------------------------------------------------------------
1955 * initialize electric fields
1956 *----------------------------------------------------------------------------*/
1957
1958 void
cs_elec_fields_initialize(const cs_mesh_t * mesh,int isuite)1959 cs_elec_fields_initialize(const cs_mesh_t *mesh,
1960 int isuite)
1961 {
1962 BFT_MALLOC(_elec_option.izreca, mesh->n_i_faces, int);
1963 for (cs_lnum_t i = 0; i < mesh->n_i_faces; i++)
1964 _elec_option.izreca[i] = 0;
1965
1966 cs_lnum_t n_cells = mesh->n_cells;
1967
1968 static int ipass = 0;
1969 ipass += 1;
1970
1971 int ielarc = cs_glob_physical_model_flag[CS_ELECTRIC_ARCS];
1972
1973 if (isuite == 0 && ipass == 1) {
1974 /* enthalpy */
1975 cs_real_t hinit = 0.;
1976 if (ielarc > 0) {
1977 cs_real_t *ym;
1978 BFT_MALLOC(ym, cs_glob_elec_properties->ngaz, cs_real_t);
1979 ym[0] = 1.;
1980 if (cs_glob_elec_properties->ngaz > 1)
1981 for (int i = 1; i < cs_glob_elec_properties->ngaz; i++)
1982 ym[i] = 0.;
1983
1984 cs_real_t tinit = cs_glob_fluid_properties->t0;
1985 hinit = _cs_elec_convert_t_to_h(ym, tinit);
1986 BFT_FREE(ym);
1987 }
1988
1989 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1990 CS_F_(h)->val[iel] = hinit;
1991 }
1992
1993 /* mass fraction of first gas */
1994 if (cs_glob_elec_properties->ngaz > 1) {
1995 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
1996 CS_FI_(ycoel, 0)->val[iel] = 1.;
1997 }
1998 }
1999 }
2000
2001 /*----------------------------------------------------------------------------
2002 * scaling electric quantities
2003 *----------------------------------------------------------------------------*/
2004
2005 void
cs_elec_scaling_function(const cs_mesh_t * mesh,const cs_mesh_quantities_t * mesh_quantities,cs_real_t * dt)2006 cs_elec_scaling_function(const cs_mesh_t *mesh,
2007 const cs_mesh_quantities_t *mesh_quantities,
2008 cs_real_t *dt)
2009 {
2010 cs_real_t *volume = mesh_quantities->cell_vol;
2011 cs_real_t *surfac = mesh_quantities->i_face_normal;
2012 cs_lnum_t n_cells = mesh->n_cells;
2013 cs_lnum_t nfac = mesh->n_i_faces;
2014
2015 double coepot = 0.;
2016 double coepoa = 1.;
2017
2018 int ieljou = cs_glob_physical_model_flag[CS_JOULE_EFFECT];
2019 int ielarc = cs_glob_physical_model_flag[CS_ELECTRIC_ARCS];
2020
2021 if (ielarc >= 1) {
2022 if (cs_glob_elec_option->modrec == 1) {
2023 /* standard model */
2024 double somje = 0.;
2025 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
2026 somje += CS_F_(joulp)->val[iel] * volume[iel];
2027
2028 cs_parall_sum(1, CS_DOUBLE, &somje);
2029
2030 coepot = cs_glob_elec_option->couimp * cs_glob_elec_option->pot_diff
2031 / CS_MAX(somje, cs_math_epzero);
2032 coepoa = coepot;
2033
2034 if (coepot > 1.5)
2035 coepot = 1.5;
2036 if (coepot < 0.75)
2037 coepot = 0.75;
2038
2039 bft_printf("imposed current / current %14.5e, scaling coef. %14.5e\n",
2040 coepoa, coepot);
2041 }
2042 else if (cs_glob_elec_option->modrec == 2) {
2043 /* restrike model */
2044 cs_gui_elec_model_rec();
2045 double elcou = 0.;
2046 cs_real_3_t *cpro_curre = (cs_real_3_t *)(CS_F_(curre)->val);
2047 if (mesh->halo != NULL)
2048 cs_halo_sync_var_strided(mesh->halo, CS_HALO_STANDARD,
2049 (cs_real_t *)cpro_curre, 3);
2050 for (cs_lnum_t ifac = 0; ifac < nfac; ifac++) {
2051 if (cs_glob_elec_option->izreca[ifac] > 0) {
2052 bool ok = true;
2053 for (int idir = 0; idir < 3; idir++)
2054 if ( fabs(surfac[3 * ifac + idir]) > 0.
2055 && idir != (cs_glob_elec_option->idreca - 1))
2056 ok = false;
2057 if (ok) {
2058 cs_lnum_t iel = mesh->i_face_cells[ifac][0];
2059 elcou += cpro_curre[iel][cs_glob_elec_option->idreca - 1]
2060 * surfac[3 * ifac + cs_glob_elec_option->idreca - 1];
2061 }
2062 }
2063 }
2064 cs_parall_sum(1, CS_DOUBLE, &elcou);
2065 if (fabs(elcou) > 1.e-6)
2066 elcou = fabs(elcou);
2067 else
2068 elcou = 0.;
2069
2070 if (fabs(elcou) > 1.e-20)
2071 coepoa = cs_glob_elec_option->couimp / elcou;
2072
2073 bft_printf("ELCOU %15.8E\n", elcou);
2074 _elec_option.elcou = elcou;
2075 }
2076
2077 if ( cs_glob_elec_option->modrec == 1
2078 || cs_glob_elec_option->modrec == 2) {
2079 double dtj = 1.e15;
2080 double dtjm = dtj;
2081 double delhsh = 0.;
2082 double cdtj = 20.;
2083
2084 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
2085 if (CS_F_(rho)->val[iel] > 0.)
2086 delhsh = CS_F_(joulp)->val[iel] * dt[iel]
2087 / CS_F_(rho)->val[iel];
2088
2089 if (fabs(delhsh) > 1.e-20)
2090 dtjm = CS_F_(h)->val[iel] / delhsh;
2091 else
2092 dtjm = dtj;
2093 dtjm = fabs(dtjm);
2094 dtj = CS_MIN(dtj, dtjm);
2095 }
2096 cs_parall_min(1, CS_DOUBLE, &dtj);
2097 bft_printf("DTJ %15.8E\n", dtj);
2098
2099 double cpmx = pow(cdtj * dtj, 0.5);
2100 coepot = cpmx;
2101
2102 if (cs_glob_time_step->nt_cur > 2) {
2103 if (coepoa > 1.05)
2104 coepot = cpmx;
2105 else
2106 coepot = coepoa;
2107 }
2108
2109 bft_printf(" Cpmx = %14.5E\n", cpmx);
2110 bft_printf(" COEPOA = %14.5E\n", coepoa);
2111 bft_printf(" COEPOT = %14.5E\n", coepot);
2112 bft_printf(" Dpot recale = %14.5E\n", _elec_option.pot_diff * coepot);
2113
2114 /* scaling electric fields */
2115 _elec_option.pot_diff *= coepot;
2116
2117 /* electric potential (for post treatment) */
2118 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
2119 CS_F_(potr)->val[iel] *= coepot;
2120
2121 /* current density */
2122 if (ielarc > 0) {
2123 cs_real_3_t *cpro_curre = (cs_real_3_t *)(CS_F_(curre)->val);
2124 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
2125 for (cs_lnum_t i = 0; i < 3; i++)
2126 cpro_curre[iel][i] *= coepot;
2127 }
2128 }
2129
2130 /* joule effect */
2131 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
2132 CS_F_(joulp)->val[iel] *= coepot * coepot;
2133 }
2134 }
2135
2136 /* joule effect */
2137 if (ieljou > 0) {
2138 /* standard model */
2139 double somje = 0.;
2140 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
2141 somje += CS_F_(joulp)->val[iel] * volume[iel];
2142
2143 cs_parall_sum(1, CS_DOUBLE, &somje);
2144
2145 coepot = cs_glob_elec_option->puisim / CS_MAX(somje, cs_math_epzero);
2146 double coefav = coepot;
2147
2148 if (coepot > 1.5)
2149 coepot = 1.5;
2150 if (coepot < 0.75)
2151 coepot = 0.75;
2152
2153 bft_printf("imposed power / sum(jE) %14.5E, scaling coef. %14.5E\n",
2154 coefav, coepot);
2155
2156 /* scaling electric fields */
2157 _elec_option.pot_diff *= coepot;
2158 _elec_option.coejou *= coepot;
2159
2160 /* electric potential (for post treatment) */
2161 if (ieljou != 3 && ieljou != 4)
2162 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
2163 CS_F_(potr)->val[iel] *= coepot;
2164
2165 /* current density */
2166 if (ieljou == 2)
2167 for (int i = 0; i < 3; i++)
2168 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
2169 CS_F_(poti)->val[iel] *= coepot;
2170
2171 /* joule effect */
2172 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
2173 CS_F_(joulp)->val[iel] *= coepot * coepot;
2174 }
2175
2176 cs_user_scaling_elec(mesh, mesh_quantities, dt);
2177 }
2178
2179 /*----------------------------------------------------------------------------*/
2180 /*!
2181 * \brief Convert enthalpy to temperature at all boundary faces.
2182 *
2183 * This handles both user and model enthalpy conversions, so can be used
2184 * safely whenever conversion is needed.
2185 *
2186 * \param[in] h enthalpy values
2187 * \param[out] t temperature values
2188 */
2189 /*----------------------------------------------------------------------------*/
2190
2191 void
cs_elec_convert_h_to_t_faces(const cs_real_t h[],cs_real_t t[])2192 cs_elec_convert_h_to_t_faces(const cs_real_t h[],
2193 cs_real_t t[])
2194 {
2195 const cs_mesh_t *m = cs_glob_mesh;
2196 const cs_lnum_t n_b_faces = m->n_b_faces;
2197
2198 const cs_data_elec_t *el_p = cs_glob_elec_properties;
2199 const int n_gasses = el_p->ngaz;
2200
2201 if (n_gasses == 1) {
2202
2203 cs_real_t ym[1] = {1.};
2204
2205 for (cs_lnum_t i = 0; i < n_b_faces; i++)
2206 t[i] = _cs_elec_convert_h_to_t(ym, h[i]);
2207
2208 }
2209 else {
2210
2211 const cs_lnum_t *b_face_cells = m->b_face_cells;
2212
2213 cs_real_t *ym;
2214 BFT_MALLOC(ym, n_gasses, cs_real_t);
2215
2216 for (cs_lnum_t f_id = 0; f_id < n_b_faces; f_id++) {
2217
2218 cs_lnum_t c_id = b_face_cells[f_id];
2219
2220 ym[n_gasses - 1] = 1.;
2221 for (int gas_id = 0; gas_id < n_gasses - 1; gas_id++) {
2222 ym[gas_id] = CS_FI_(ycoel, gas_id)->val[c_id];
2223 ym[n_gasses - 1] -= ym[gas_id];
2224 }
2225
2226 t[f_id] = _cs_elec_convert_h_to_t(ym, h[f_id]);
2227
2228 }
2229
2230 BFT_FREE(ym);
2231
2232 }
2233 }
2234
2235 /*----------------------------------------------------------------------------*/
2236 /*!
2237 * \brief Convert temperature to enthalpy at all cells.
2238 *
2239 * This handles both user and model temperature conversions, so can be used
2240 * safely whenever conversion is needed.
2241 *
2242 * \param[in] t temperature values
2243 * \param[out] h enthalpy values
2244 */
2245 /*----------------------------------------------------------------------------*/
2246
2247 void
cs_elec_convert_t_to_h_cells(const cs_real_t t[],cs_real_t h[])2248 cs_elec_convert_t_to_h_cells(const cs_real_t t[],
2249 cs_real_t h[])
2250 {
2251 const cs_mesh_t *m = cs_glob_mesh;
2252 const cs_lnum_t n_cells = m->n_cells;
2253
2254 const cs_data_elec_t *el_p = cs_glob_elec_properties;
2255 const int n_gasses = el_p->ngaz;
2256
2257 if (n_gasses == 1) {
2258
2259 cs_real_t ym[1] = {1.};
2260
2261 for (cs_lnum_t i = 0; i < n_cells; i++)
2262 h[i] = _cs_elec_convert_t_to_h(ym, t[i]);
2263
2264 }
2265 else {
2266
2267 cs_real_t *ym;
2268 BFT_MALLOC(ym, n_gasses, cs_real_t);
2269
2270 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
2271
2272 ym[n_gasses - 1] = 1.;
2273 for (int gas_id = 0; gas_id < n_gasses - 1; gas_id++) {
2274 ym[gas_id] = CS_FI_(ycoel, gas_id)->val[c_id];
2275 ym[n_gasses - 1] -= ym[gas_id];
2276 }
2277
2278 h[c_id] = _cs_elec_convert_t_to_h(ym, t[c_id]);
2279
2280 }
2281
2282 BFT_FREE(ym);
2283
2284 }
2285 }
2286
2287 /*----------------------------------------------------------------------------*/
2288 /*!
2289 * \brief Convert temperature to enthalpy at selected boundary faces.
2290 *
2291 * This handles both user and model temperature conversions, so can be used
2292 * safely whenever conversion is needed.
2293 *
2294 * \param[in] n_faces number of selected faces
2295 * \param[in] face_ids ids of selected faces
2296 * \param[in] t temperature values (defined on all boundary faces)
2297 * \param[out] h enthalpy values (defined on all boundary faces)
2298 */
2299 /*----------------------------------------------------------------------------*/
2300
2301 void
cs_elec_convert_t_to_h_faces(const cs_lnum_t n_faces,const cs_lnum_t face_ids[],const cs_real_t t[],cs_real_t h[])2302 cs_elec_convert_t_to_h_faces(const cs_lnum_t n_faces,
2303 const cs_lnum_t face_ids[],
2304 const cs_real_t t[],
2305 cs_real_t h[])
2306 {
2307 const cs_mesh_t *m = cs_glob_mesh;
2308
2309 const cs_data_elec_t *el_p = cs_glob_elec_properties;
2310 const int n_gasses = el_p->ngaz;
2311
2312 if (n_gasses == 1) {
2313
2314 cs_real_t ym[1] = {1.};
2315
2316 for (cs_lnum_t i = 0; i < n_faces; i++) {
2317 cs_lnum_t f_id = face_ids[i];
2318 h[f_id] = _cs_elec_convert_t_to_h(ym, t[f_id]);
2319 }
2320
2321 }
2322 else {
2323
2324 const cs_lnum_t *b_face_cells = m->b_face_cells;
2325
2326 cs_real_t *ym;
2327 BFT_MALLOC(ym, n_gasses, cs_real_t);
2328
2329 for (cs_lnum_t i = 0; i < n_faces; i++) {
2330
2331 cs_lnum_t f_id = face_ids[i];
2332 cs_lnum_t c_id = b_face_cells[f_id];
2333 for (int gas_id = 0; gas_id < n_gasses - 1; gas_id++) {
2334 ym[gas_id] = CS_FI_(ycoel, gas_id)->val[c_id];
2335 ym[n_gasses - 1] -= ym[gas_id];
2336 }
2337
2338 h[f_id] = _cs_elec_convert_t_to_h(ym, t[f_id]);
2339
2340 }
2341
2342 BFT_FREE(ym);
2343
2344 }
2345 }
2346
2347 /*----------------------------------------------------------------------------*/
2348
2349 END_C_DECLS
2350