1 /*============================================================================
2 * Enthaly to and from temperature conversion.
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 <stdarg.h>
36 #include <stdlib.h>
37 #include <stdio.h>
38 #include <string.h>
39
40 /*----------------------------------------------------------------------------
41 * Local headers
42 *----------------------------------------------------------------------------*/
43
44 #include "bft_mem.h"
45 #include "bft_error.h"
46
47 #include "cs_elec_model.h"
48 #include "cs_field.h"
49 #include "cs_field_pointer.h"
50 #include "cs_log.h"
51 #include "cs_mesh_location.h"
52 #include "cs_physical_constants.h"
53 #include "cs_physical_model.h"
54 #include "cs_prototypes.h"
55 #include "cs_volume_zone.h"
56
57 /*----------------------------------------------------------------------------
58 * Header for the current file
59 *----------------------------------------------------------------------------*/
60
61 #include "cs_ht_convert.h"
62
63 /*----------------------------------------------------------------------------*/
64
65 BEGIN_C_DECLS
66
67 /*=============================================================================
68 * Additional doxygen documentation
69 *============================================================================*/
70
71 /*!
72 \file cs_ht_convert.c
73 Enthalpy to and from temperature conversion.
74
75 Other fields may be involved in the conversion.
76
77 TODO: when possible (based on calling functions's conversion to C)
78 a function pointer-based logic would allow migrating this
79 functionnality to cs_physical_properties, without adding
80 high level physical model dependencies (i.e. cs_physical_model.h)
81 to that lower-level API.
82 */
83
84 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
85
86 /*============================================================================
87 * Prototypes for Fortran subroutines
88 *============================================================================*/
89
90 extern void CS_PROCF(coh2tb, COH2TB)
91 (
92 const cs_real_t *h,
93 cs_real_t *t
94 );
95
96 extern void CS_PROCF(cot2hb, COT2HB)
97 (
98 const cs_lnum_t *n_faces,
99 const cs_lnum_t *face_ids,
100 const cs_real_t *t,
101 cs_real_t *h
102 );
103
104 /*=============================================================================
105 * Local macro definitions
106 *============================================================================*/
107
108 /*============================================================================
109 * Private function definitions
110 *============================================================================*/
111
112 /*============================================================================
113 * Type definitions
114 *============================================================================*/
115
116 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
117
118 /*=============================================================================
119 * Public function definitions
120 *============================================================================*/
121
122 /*----------------------------------------------------------------------------*/
123 /*!
124 * \brief Convert enthalpy to temperature at all cells.
125 *
126 * This handles both user and model enthalpy conversions, so can be used
127 * safely whenever conversion is needed.
128 *
129 * \param[in] h enthalpy values
130 * \param[out] t temperature values
131 */
132 /*----------------------------------------------------------------------------*/
133
134 void
cs_ht_convert_h_to_t_cells(const cs_real_t h[],cs_real_t t[])135 cs_ht_convert_h_to_t_cells(const cs_real_t h[],
136 cs_real_t t[])
137 {
138 const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
139
140 const int *pm_flag = cs_glob_physical_model_flag;
141 const cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
142 bool need_solid_default = (mq->has_disable_flag) ? true : false;
143
144 cs_real_t *cpro_t = NULL;
145
146 /* Gas combustion: premixed or diffusion flame */
147 if ( pm_flag[CS_COMBUSTION_EBU] >= 0
148 || pm_flag[CS_COMBUSTION_3PT] >= 0)
149 cpro_t = CS_F_(t)->val;
150
151 /* Pulverized coal or fuel combustion */
152 else if ( pm_flag[CS_COMBUSTION_COAL] >= 0
153 || pm_flag[CS_COMBUSTION_FUEL] >= 0)
154 cpro_t = CS_F_(t)->val;
155
156 /* Electric arcs */
157 else if ( pm_flag[CS_JOULE_EFFECT] >= 1
158 || pm_flag[CS_ELECTRIC_ARCS] >= 1)
159 cpro_t = CS_F_(t)->val;
160
161 /* When temperature maintained by model is available
162 ------------------------------------------------- */
163
164 if (cpro_t != NULL) {
165
166 for (cs_lnum_t i = 0; i < n_cells; i++)
167 t[i] = cpro_t[i];
168
169 }
170
171 /* Default for other cases
172 ----------------------- */
173
174 else {
175
176 const cs_field_t *f_cp = cs_field_by_name_try("specific_heat");
177 if (f_cp != NULL) {
178 const cs_real_t *cpro_cp = f_cp->val;
179 for (cs_lnum_t i = 0; i < n_cells; i++)
180 t[i] = h[i] / cpro_cp[i];
181 }
182 else {
183 const double cp0 = cs_glob_fluid_properties->cp0;
184 for (cs_lnum_t i = 0; i < n_cells; i++)
185 t[i] = h[i] / cp0;
186 }
187
188 need_solid_default = false;
189
190 }
191
192 /* Handle solid zones */
193
194 if (need_solid_default) {
195
196 const int *disable_flag = mq->c_disable_flag;
197
198 const cs_field_t *f_cp = cs_field_by_name_try("specific_heat");
199 if (f_cp != NULL) {
200 const cs_real_t *cpro_cp = f_cp->val;
201 for (cs_lnum_t i = 0; i < n_cells; i++) {
202 if (disable_flag[i])
203 t[i] = h[i] / cpro_cp[i];
204 }
205 }
206 else {
207 const double cp0 = cs_glob_fluid_properties->cp0;
208 for (cs_lnum_t i = 0; i < n_cells; i++) {
209 if (disable_flag[i])
210 t[i] = h[i] / cp0;
211 }
212 }
213
214 }
215
216 /* Allow user functions */
217
218 int n_zones = cs_volume_zone_n_zones();
219 for (int z_id = 0; z_id < n_zones; z_id++) {
220 const cs_zone_t *z = cs_volume_zone_by_id(z_id);
221
222 /* Note: we could restrict this to
223 z->type & CS_VOLUME_ZONE_PHYSICAL_PROPERTIES
224 but the user can also easily handle this */
225
226 cs_user_physical_properties_h_to_t(cs_glob_domain,
227 z,
228 false, /* z_local */
229 h,
230 t);
231
232 }
233 }
234
235 /*----------------------------------------------------------------------------*/
236 /*!
237 * \brief Convert enthalpy to temperature at solid cells only.
238 *
239 * This handles both user and model enthalpy conversions, so can be used
240 * safely whenever conversion is needed.
241 */
242 /*----------------------------------------------------------------------------*/
243
244 void
cs_ht_convert_h_to_t_cells_solid(void)245 cs_ht_convert_h_to_t_cells_solid(void)
246 {
247 const cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
248 if (mq->has_disable_flag == 0 || CS_F_(h) == NULL || CS_F_(t) == NULL)
249 return;
250
251 const cs_real_t *h = CS_F_(h)->val;
252 cs_real_t *t = CS_F_(t)->val;
253
254 int n_zones = cs_volume_zone_n_zones();
255 for (int z_id = 0; z_id < n_zones; z_id++) {
256 const cs_zone_t *z = cs_volume_zone_by_id(z_id);
257
258 if ( z->type & CS_VOLUME_ZONE_SOLID
259 && z->type & CS_VOLUME_ZONE_PHYSICAL_PROPERTIES) {
260
261 const cs_field_t *f_cp = cs_field_by_name_try("specific_heat");
262 if (f_cp != NULL) {
263 const cs_real_t *cpro_cp = f_cp->val;
264 for (cs_lnum_t i = 0; i < z->n_elts; i++) {
265 cs_lnum_t c_id = z->elt_ids[i];
266 t[c_id] = h[c_id] / cpro_cp[c_id];
267 }
268 }
269 else {
270 const double cp0 = cs_glob_fluid_properties->cp0;
271 for (cs_lnum_t i = 0; i < z->n_elts; i++) {
272 cs_lnum_t c_id = z->elt_ids[i];
273 t[c_id] = h[c_id] / cp0;
274 }
275 }
276
277 cs_user_physical_properties_h_to_t(cs_glob_domain,
278 z,
279 false, /* z_local */
280 h,
281 t);
282
283 }
284
285 }
286 }
287
288 /*----------------------------------------------------------------------------*/
289 /*!
290 * \brief Convert enthalpy to temperature at all boundary faces.
291 *
292 * This handles both user and model enthalpy conversions, so can be used
293 * safely whenever conversion is needed.
294 *
295 * \param[in] h enthalpy values
296 * \param[out] t temperature values
297 */
298 /*----------------------------------------------------------------------------*/
299
300 void
cs_ht_convert_h_to_t_faces(const cs_real_t h[],cs_real_t t[])301 cs_ht_convert_h_to_t_faces(const cs_real_t h[],
302 cs_real_t t[])
303 {
304 const cs_mesh_t *m = cs_glob_mesh;
305 const cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
306 const cs_lnum_t *b_face_cells = m->b_face_cells;
307 const cs_lnum_t n_b_faces = m->n_b_faces;
308 const cs_lnum_t has_dc = mq->has_disable_flag;
309
310 bool need_solid_default = (has_dc) ? true : false;
311
312 const int *pm_flag = cs_glob_physical_model_flag;
313
314 /* Gas combustion: premixed or diffusion flame */
315 if ( pm_flag[CS_COMBUSTION_EBU] >= 0
316 || pm_flag[CS_COMBUSTION_3PT] >= 0)
317 CS_PROCF(coh2tb, COH2TB)(h, t);
318
319 /* Pulverized coal combustion */
320 else if (pm_flag[CS_COMBUSTION_COAL] >= 0)
321 cs_coal_thfieldconv1(CS_MESH_LOCATION_BOUNDARY_FACES, h, t);
322
323 /* Pulverized fuel combustion */
324 else if (pm_flag[CS_COMBUSTION_FUEL] >= 0)
325 cs_fuel_thfieldconv1(CS_MESH_LOCATION_BOUNDARY_FACES, h, t);
326
327 /* Electric arcs */
328 else if (pm_flag[CS_JOULE_EFFECT] < 1 && pm_flag[CS_ELECTRIC_ARCS] >= 1)
329 cs_elec_convert_h_to_t_faces(h, t);
330
331 /* Default for other cases
332 ----------------------- */
333
334 else {
335
336 const cs_field_t *f_cp = cs_field_by_name_try("specific_heat");
337 if (f_cp != NULL) {
338 const cs_real_t *cpro_cp = f_cp->val;
339 for (cs_lnum_t i = 0; i < n_b_faces; i++) {
340 cs_lnum_t c_id = b_face_cells[i];
341 t[i] = h[i] / cpro_cp[c_id];
342 }
343 }
344 else {
345 const double cp0 = cs_glob_fluid_properties->cp0;
346 for (cs_lnum_t i = 0; i < n_b_faces; i++)
347 t[i] = h[i] / cp0;
348 }
349
350 need_solid_default = false;
351
352 }
353
354 /* Default for solid zones
355 ----------------------- */
356
357 if (need_solid_default) {
358
359 assert(has_dc == 1);
360 const int *disable_flag = mq->c_disable_flag;
361
362 const cs_field_t *f_cp = cs_field_by_name_try("specific_heat");
363 if (f_cp != NULL) {
364 const cs_real_t *cpro_cp = f_cp->val;
365 for (cs_lnum_t i = 0; i < n_b_faces; i++) {
366 cs_lnum_t c_id = b_face_cells[i];
367 if (disable_flag[c_id])
368 t[i] = h[i] / cpro_cp[c_id];
369 }
370 }
371 else {
372 const double cp0 = cs_glob_fluid_properties->cp0;
373 for (cs_lnum_t i = 0; i < n_b_faces; i++) {
374 cs_lnum_t c_id = b_face_cells[i];
375 if (disable_flag[c_id])
376 t[i] = h[i] / cp0;
377 }
378 }
379
380 }
381
382 /* Allow user functions */
383
384 int n_zones = cs_boundary_zone_n_zones();
385 for (int z_id = 0; z_id < n_zones; z_id++) {
386 const cs_zone_t *z = cs_boundary_zone_by_id(z_id);
387
388 cs_user_physical_properties_h_to_t(cs_glob_domain,
389 z,
390 false, /* z_local */
391 h,
392 t);
393
394 }
395 }
396
397 /*----------------------------------------------------------------------------*/
398 /*!
399 * \brief Convert temperature to enthalpy at selected boundary faces.
400 *
401 * This handles both user and model enthalpy conversions, so can be used
402 * safely whenever conversion is needed.
403 *
404 * \param[in] n_faces number of selected boundary faces
405 * \param[in] face_ids list of associated face ids
406 * \param[in] t temperature values
407 * \param[out] h enthalpy values
408 */
409 /*----------------------------------------------------------------------------*/
410
411 void
cs_ht_convert_t_to_h_faces_l(cs_lnum_t n_faces,const cs_lnum_t face_ids[],const cs_real_t t[],cs_real_t h[])412 cs_ht_convert_t_to_h_faces_l(cs_lnum_t n_faces,
413 const cs_lnum_t face_ids[],
414 const cs_real_t t[],
415 cs_real_t h[])
416 {
417 const cs_mesh_t *m = cs_glob_mesh;
418 const cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
419 const cs_lnum_t *b_face_cells = m->b_face_cells;
420 const cs_lnum_t n_b_faces = m->n_b_faces;
421 const cs_lnum_t has_dc = mq->has_disable_flag;
422
423 const int *pm_flag = cs_glob_physical_model_flag;
424
425 bool need_solid_default = (has_dc) ? true : false;
426
427 /* Gas combustion: premixed or diffusion flame */
428 if ( pm_flag[CS_COMBUSTION_EBU] >= 0
429 || pm_flag[CS_COMBUSTION_3PT] >= 0)
430 CS_PROCF(cot2hb, COT2HB)(&n_faces, face_ids, t, h);
431
432 /* Pulverized coal combustion */
433 else if (pm_flag[CS_COMBUSTION_COAL] >= 0)
434 cs_coal_bt2h(n_faces, face_ids, t, h);
435
436 /* Pulverized fuel combustion */
437 else if (pm_flag[CS_COMBUSTION_FUEL] >= 0)
438 cs_fuel_bt2h(n_faces, face_ids, t, h);
439
440 /* Electric arcs */
441 else if (pm_flag[CS_JOULE_EFFECT] < 1 && pm_flag[CS_ELECTRIC_ARCS] >= 1)
442 cs_elec_convert_t_to_h_faces(n_faces, face_ids, t, h);
443
444 /* Default for other cases
445 ----------------------- */
446
447 else {
448
449 const cs_field_t *f_cp = cs_field_by_name_try("specific_heat");
450 if (f_cp != NULL) {
451 const cs_real_t *cpro_cp = f_cp->val;
452 for (cs_lnum_t i = 0; i < n_faces; i++) {
453 cs_lnum_t f_id = face_ids[i];
454 cs_lnum_t c_id = b_face_cells[f_id];
455 h[f_id] = t[f_id] * cpro_cp[c_id];
456 }
457 }
458 else {
459 const double cp0 = cs_glob_fluid_properties->cp0;
460 for (cs_lnum_t i = 0; i < n_faces; i++) {
461 cs_lnum_t f_id = face_ids[i];
462 h[f_id] = t[f_id] * cp0;
463 }
464 }
465
466 need_solid_default = false;
467
468 }
469
470 /* Default for solid zones
471 ----------------------- */
472
473 if (need_solid_default) {
474
475 assert(has_dc == 1);
476 const int *disable_flag = mq->c_disable_flag;
477
478 const cs_field_t *f_cp = cs_field_by_name_try("specific_heat");
479 if (f_cp != NULL) {
480 const cs_real_t *cpro_cp = f_cp->val;
481 for (cs_lnum_t i = 0; i < n_faces; i++) {
482 cs_lnum_t f_id = face_ids[i];
483 cs_lnum_t c_id = b_face_cells[f_id];
484 if (disable_flag[c_id])
485 h[f_id] = t[f_id] * cpro_cp[c_id];
486 }
487 }
488 else {
489 const double cp0 = cs_glob_fluid_properties->cp0;
490 for (cs_lnum_t i = 0; i < n_faces; i++) {
491 cs_lnum_t f_id = face_ids[i];
492 cs_lnum_t c_id = b_face_cells[f_id];
493 if (disable_flag[c_id])
494 h[f_id] = t[f_id] * cp0;
495 }
496 }
497
498 }
499
500 /* Allow user functions */
501
502 cs_real_t *h_f;
503 BFT_MALLOC(h_f, n_b_faces, cs_real_t);
504 for (cs_lnum_t i = 0; i < n_b_faces; i++)
505 h_f[i] = h[i];
506
507 int n_zones = cs_boundary_zone_n_zones();
508 for (int z_id = 0; z_id < n_zones; z_id++) {
509 const cs_zone_t *z = cs_boundary_zone_by_id(z_id);
510
511 cs_user_physical_properties_t_to_h(cs_glob_domain,
512 z,
513 false, /* z_local */
514 t,
515 h_f);
516 }
517
518 for (cs_lnum_t i = 0; i < n_faces; i++) {
519 cs_lnum_t f_id = face_ids[i];
520 h[f_id] = h_f[f_id];
521 }
522
523 BFT_FREE(h_f);
524 }
525
526 /*----------------------------------------------------------------------------*/
527
528 END_C_DECLS
529