1 /*============================================================================
2 * Thermodynamic laws for the compressible module
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27
28 #include "cs_defs.h"
29
30 /*----------------------------------------------------------------------------
31 * Standard C library headers
32 *----------------------------------------------------------------------------*/
33
34 #include <assert.h>
35 #include <stdio.h>
36 #include <stdlib.h>
37 #include <string.h>
38 #include <math.h>
39
40 /*----------------------------------------------------------------------------
41 * Local headers
42 *----------------------------------------------------------------------------*/
43
44 #include "bft_mem.h"
45 #include "bft_printf.h"
46
47 #include "cs_field.h"
48 #include "cs_field_pointer.h"
49 #include "cs_math.h"
50 #include "cs_mesh.h"
51 #include "cs_mesh_quantities.h"
52 #include "cs_parall.h"
53 #include "cs_parameters.h"
54 #include "cs_thermal_model.h"
55 #include "cs_cf_model.h"
56 #include "cs_hgn_thermo.h"
57
58 /*----------------------------------------------------------------------------
59 * Header for the current file
60 *----------------------------------------------------------------------------*/
61
62 #include "cs_cf_thermo.h"
63
64 /*----------------------------------------------------------------------------*/
65
66 BEGIN_C_DECLS
67
68 /*=============================================================================
69 * Additional doxygen documentation
70 *============================================================================*/
71
72 /*!
73 \file cs_cf_thermo.c
74 Define thermodynamic laws for the compressible flow module.
75 Only the perfect gas law is available for now. The molar mass has to be
76 provided either in the GUI or in cs_user_parameters.f90.
77 */
78
79 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
80
81 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
82
83 /*=============================================================================
84 * Public function definitions
85 *============================================================================*/
86
87 /*----------------------------------------------------------------------------*/
88 /*!
89 * \brief Set variability of isobaric specific heat and isochoric specific heat
90 * according to the chosen thermodynamic law.
91 */
92 /*----------------------------------------------------------------------------*/
93
94 void
cs_cf_set_thermo_options(void)95 cs_cf_set_thermo_options(void)
96 {
97 cs_fluid_properties_t *fluid_properties = cs_get_glob_fluid_properties();
98 int ieos = cs_glob_cf_model->ieos;
99 if (ieos == CS_EOS_IDEAL_GAS || ieos == CS_EOS_STIFFENED_GAS) {
100 /* Calculation options: constant Cp and Cv (perfect or stiffened gas)
101 specific heat Cv0 is calculated in a subsequent section (from Cp0) */
102 fluid_properties->icp = -1;
103 fluid_properties->icv = -1;
104 }
105 else if (ieos == CS_EOS_GAS_MIX) {
106 /* variable Cp and Cv for ideal gas mix eos. */
107 fluid_properties->icp = 0;
108 fluid_properties->icv = 0;
109 }
110 }
111
112 /*----------------------------------------------------------------------------*/
113 /*!
114 * \brief Initialize density, total energy and isochoric specific heat
115 * according to the chosen thermodynamic law using the default parameters.
116 */
117 /*----------------------------------------------------------------------------*/
118
119 void
cs_cf_thermo_default_init(void)120 cs_cf_thermo_default_init(void)
121 {
122 /* Local variables */
123 cs_real_t e0;
124
125 const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
126
127 cs_real_t r_pg = cs_physical_constants_r;
128 cs_real_t psginf= cs_glob_cf_model->psginf;
129 cs_real_t p0 = cs_glob_fluid_properties->p0;
130 cs_real_t t0 = cs_glob_fluid_properties->t0;
131 cs_real_t cp0 = cs_glob_fluid_properties->cp0;
132
133 cs_fluid_properties_t *fluid_properties = cs_get_glob_fluid_properties();
134
135 cs_real_t *cv0 = &fluid_properties->cv0;
136 cs_real_t *ro0 = &fluid_properties->ro0;
137
138 /* Default initializations
139 t0 is positive (this assumption has been checked in verini) */
140 cs_real_t *crom = CS_F_(rho)->val;
141 cs_real_t *cvar_en = CS_F_(e_tot)->val;
142
143 int ieos = cs_glob_cf_model->ieos;
144
145 /* perfect gas */
146 if (ieos == CS_EOS_IDEAL_GAS || ieos == CS_EOS_GAS_MIX) {
147 cs_real_t xmasml = cs_glob_fluid_properties->xmasmr;
148 *cv0 = cp0 - r_pg/xmasml;
149 *ro0 = p0 * xmasml/(r_pg*t0);
150 e0 = *cv0 * t0;
151 }
152 /* stiffened gas: cv0 is set by the user */
153 else if (ieos == CS_EOS_STIFFENED_GAS) {
154 cs_real_t gamma = cs_glob_cf_model->gammasg;
155 *ro0 = (p0 + psginf) / ((gamma-1.)*(*cv0)*t0);
156 e0 = *cv0*t0 + psginf / *ro0;
157 }
158 /* homogeneous two-phase TODOHGN */
159 else if (ieos == CS_EOS_HOMOGENEOUS_TWO_PHASE) {
160 *cv0 = 1.;
161 *ro0 = 1.;
162 e0 = 1.;
163 }
164 else {
165 e0 = 0;
166 }
167
168 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
169 crom[cell_id] = *ro0;
170 cvar_en[cell_id] = e0;
171 }
172 }
173
174 /*----------------------------------------------------------------------------*/
175 /*!
176 * \brief Check the positivity of the pressure.
177 *
178 * \param[in] pres array of pressure values
179 * \param[in] l_size l_size of the array
180 */
181 /*----------------------------------------------------------------------------*/
182 // FIXME: the check function should be generalized (pass the name as argument).
183
184 void
cs_cf_check_pressure(cs_real_t * pres,cs_lnum_t l_size)185 cs_cf_check_pressure(cs_real_t *pres,
186 cs_lnum_t l_size)
187 {
188 cs_real_t psginf = cs_glob_cf_model->psginf;
189
190 /* Local variables */
191 cs_gnum_t ierr;
192
193 /* If the pressure is lower or equal to zero, stop the computation.
194 Indeed, if this is the case, the thermodynamic computations will most
195 probably fail. This call is done at the end of the density calculation */
196 ierr = 0;
197 for (cs_lnum_t ii = 0; ii < l_size; ii++)
198 if (pres[ii] <= -psginf+cs_math_epzero)
199 ierr = ierr + 1;
200
201 if (cs_glob_rank_id >= 0) cs_parall_counter(&ierr, 1);
202
203 /* TODO check if message is OK in stiffened gas ("real p" = p+psginf??) */
204 /* Which pressure should be post-processed ? */
205 if (ierr > 0)
206 bft_error(__FILE__, __LINE__, 0,
207 _("Error in thermodynamics computations for compressible flows\n"
208 ":\n"
209 "Negative values of the pressure were encountered in %lu"
210 " cells.\n"), ierr);
211 }
212
213 /*----------------------------------------------------------------------------*/
214 /*!
215 * \brief Check the positivity of the internal energy.
216
217 * \param[in] ener array of total energy values
218 * \param[in] l_size l_size of the array
219 * \param[in] vel array of velocity values
220 */
221 /*----------------------------------------------------------------------------*/
222
223 void
cs_cf_check_internal_energy(cs_real_t * ener,cs_lnum_t l_size,cs_real_3_t * vel)224 cs_cf_check_internal_energy(cs_real_t *ener,
225 cs_lnum_t l_size,
226 cs_real_3_t *vel)
227 {
228 /* Local variables */
229 cs_gnum_t ierr;
230 cs_real_t enint;
231
232 /* If the internal energy <= zero: stop the computation.
233 Indeed, if this is the case, the thermodynamic computations will
234 most probably fail. */
235 ierr = 0;
236 for (cs_lnum_t ii = 0; ii < l_size; ii++) {
237 cs_real_t v2 = cs_math_3_square_norm(vel[ii]);
238 enint = ener[ii] - 0.5*v2;
239
240 if (enint <= cs_math_epzero)
241 ierr++;
242 }
243
244 if (cs_glob_rank_id >= 0)
245 cs_parall_counter(&ierr, 1);
246
247 if (ierr > 0)
248 bft_error(__FILE__, __LINE__, 0,
249 _("Error in thermodynamics computations for compressible flows\n"
250 ":\n"
251 "Negative values of the internal energy were encountered in %lu"
252 " cells.\n"), ierr);
253 }
254
255 /*----------------------------------------------------------------------------*/
256 /*!
257 * \brief Check the positivity of the density given by the user.
258 *
259 * \param[in] dens array of density values
260 * \param[in] l_size l_size of the array
261 */
262 /*----------------------------------------------------------------------------*/
263
264 void
cs_cf_check_density(cs_real_t * dens,cs_lnum_t l_size)265 cs_cf_check_density(cs_real_t *dens,
266 cs_lnum_t l_size)
267 {
268 /* Local variables */
269 cs_gnum_t ierr;
270
271 /* Verification of the values of the density
272 Stop if a negative value is detected (since the density has been
273 provided by the user, one potential cause is a wrong user
274 initialization). */
275 ierr = 0;
276 for (cs_lnum_t ii = 0; ii < l_size; ii++)
277 if (dens[ii] <= cs_math_epzero)
278 ierr = ierr + 1;
279
280 if (cs_glob_rank_id >= 0) cs_parall_counter(&ierr, 1);
281
282 if (ierr > 0)
283 bft_error(__FILE__, __LINE__, 0,
284 _("Error in thermodynamics computations for compressible flows\n"
285 ":\n"
286 "Negative values of the density were encountered in %lu"
287 " cells.\n"), ierr);
288 }
289
290 /*----------------------------------------------------------------------------*/
291 /*!
292 * \brief Check strict positivity of temperature (Celsius) given by the user.
293
294 * \param[in] temp array of temperature values
295 * \param[in] l_size l_size of the array
296 */
297 /*----------------------------------------------------------------------------*/
298
299 void
cs_cf_check_temperature(cs_real_t * temp,cs_lnum_t l_size)300 cs_cf_check_temperature(cs_real_t *temp,
301 cs_lnum_t l_size)
302 {
303 /* Local variables */
304 cs_gnum_t ierr;
305
306 /* Verification of the values of the temperature
307 Stop if a negative value is detected (since the temperature has been
308 provided by the user, one potential cause is a wrong user
309 initialization). */
310 ierr = 0;
311 for (cs_lnum_t ii = 0; ii < l_size; ii++)
312 if (temp[ii] <= cs_math_epzero)
313 ierr++;
314
315 if (cs_glob_rank_id >= 0) cs_parall_counter(&ierr, 1);
316
317 if (ierr > 0)
318 bft_error(__FILE__, __LINE__, 0,
319 _("Error in thermodynamics computations for compressible flows\n"
320 ":\n"
321 "Negative values of the temperature were encountered in %lu"
322 " cells.\n"), ierr);
323 }
324
325 /*----------------------------------------------------------------------------*/
326 /*!
327 * \brief Compute temperature and total energy from density and pressure.
328 *
329 * \param[in] cp array of isobaric specific heat values
330 * \param[in] cv array of isochoric specific heat values
331 * \param[in] pres array of pressure values
332 * \param[in] dens array of density values
333 * \param[out] temp array of temperature values
334 * \param[out] ener array of total energy values
335 * \param[in] vel array of velocity component values
336 * \param[in] l_size l_size of the array
337 */
338 /*----------------------------------------------------------------------------*/
339
340 void
cs_cf_thermo_te_from_dp(cs_real_t * cp,cs_real_t * cv,cs_real_t * pres,cs_real_t * dens,cs_real_t * temp,cs_real_t * ener,cs_real_3_t * vel,cs_lnum_t l_size)341 cs_cf_thermo_te_from_dp(cs_real_t *cp,
342 cs_real_t *cv,
343 cs_real_t *pres,
344 cs_real_t *dens,
345 cs_real_t *temp,
346 cs_real_t *ener,
347 cs_real_3_t *vel,
348 cs_lnum_t l_size)
349 {
350 /* local variables */
351 int ieos = cs_glob_cf_model->ieos;
352
353 /* calculation of temperature and energy from pressure and density */
354
355 /* single ideal gas or stiffened gas eos - constant gamma */
356 if (ieos == CS_EOS_IDEAL_GAS || ieos == CS_EOS_STIFFENED_GAS) {
357 cs_real_t gamma0;
358 cs_real_t cp0 = cs_glob_fluid_properties->cp0;
359 cs_real_t cv0 = cs_glob_fluid_properties->cv0;
360 cs_real_t psginf = cs_glob_cf_model->psginf;
361 cs_lnum_t l_size0 = 1;
362
363 cs_cf_thermo_gamma(&cp0, &cv0, &gamma0, l_size0);
364
365 for (cs_lnum_t ii = 0; ii < l_size; ii++) {
366 /* temperature */
367 temp[ii] = (pres[ii]+psginf) / ((gamma0-1.)*dens[ii]*cv0);
368 /* total energy */
369 cs_real_t v2 = cs_math_3_square_norm(vel[ii]);
370 ener[ii] = (pres[ii]+gamma0*psginf) / ((gamma0-1.)*dens[ii]) + 0.5*v2;
371 }
372 }
373 /* ideal gas mixture */
374 else if (ieos == CS_EOS_GAS_MIX) {
375 cs_real_t *gamma;
376 cs_real_t psginf = cs_glob_cf_model->psginf;
377
378 BFT_MALLOC(gamma, l_size, cs_real_t);
379
380 cs_cf_thermo_gamma(cp, cv, gamma, l_size);
381
382 for (cs_lnum_t ii = 0; ii < l_size; ii++) {
383 /* temperature */
384 temp[ii] = (pres[ii]+psginf) / ((gamma[ii]-1.)*dens[ii]*cv[ii]);
385 /* total energy */
386 cs_real_t v2 = cs_math_3_square_norm(vel[ii]);
387 ener[ii] = (pres[ii]+gamma[ii]*psginf) / ((gamma[ii]-1.)*dens[ii])
388 + 0.5*v2;
389 }
390
391 BFT_FREE(gamma);
392 }
393 }
394
395 /*----------------------------------------------------------------------------*/
396 /*!
397 * \brief Compute density and total energy from pressure and temperature
398 *
399 * \param[in] cp array of isobaric specific heat values
400 * \param[in] cv array of isochoric specific heat values
401 * \param[in] pres array of pressure values
402 * \param[in] temp array of temperature values
403 * \param[out] dens array of density values
404 * \param[out] ener array of total energy values
405 * \param[in] vel array of velocity component values
406 * \param[in] l_size l_size of the array
407 */
408 /*----------------------------------------------------------------------------*/
409
410 void
cs_cf_thermo_de_from_pt(cs_real_t * cp,cs_real_t * cv,cs_real_t * pres,cs_real_t * temp,cs_real_t * dens,cs_real_t * ener,cs_real_3_t * vel,cs_lnum_t l_size)411 cs_cf_thermo_de_from_pt(cs_real_t *cp,
412 cs_real_t *cv,
413 cs_real_t *pres,
414 cs_real_t *temp,
415 cs_real_t *dens,
416 cs_real_t *ener,
417 cs_real_3_t *vel,
418 cs_lnum_t l_size)
419 {
420 /* Local variables */
421 int ieos = cs_glob_cf_model->ieos;
422
423 /* single ideal gas or stiffened gas eos - constant gamma */
424 if (ieos == CS_EOS_IDEAL_GAS || ieos == CS_EOS_STIFFENED_GAS) {
425 cs_real_t gamma0;
426 cs_real_t cp0 = cs_glob_fluid_properties->cp0;
427 cs_real_t cv0 = cs_glob_fluid_properties->cv0;
428 cs_real_t psginf = cs_glob_cf_model->psginf;
429 cs_lnum_t l_size0 = 1;
430
431 cs_cf_thermo_gamma(&cp0, &cv0, &gamma0, l_size0);
432
433 for (cs_lnum_t ii = 0; ii < l_size; ii++) {
434 /* Density */
435 dens[ii] = (pres[ii]+psginf) / ((gamma0-1.)*temp[ii]*cv0);
436 /* Total energy */
437 cs_real_t v2 = cs_math_3_square_norm(vel[ii]);
438 ener[ii] = (pres[ii]+gamma0*psginf) / ((gamma0-1.)*dens[ii]) + 0.5*v2;
439 }
440 }
441 /* ideal gas mixture */
442 else if (ieos == CS_EOS_GAS_MIX) {
443 cs_real_t *gamma;
444 cs_real_t psginf = cs_glob_cf_model->psginf;
445
446 BFT_MALLOC(gamma, l_size, cs_real_t);
447
448 cs_cf_thermo_gamma(cp, cv, gamma, l_size);
449
450 for (cs_lnum_t ii = 0; ii < l_size; ii++) {
451 /* Density */
452 dens[ii] = (pres[ii]+psginf) / ((gamma[ii]-1.)*temp[ii]*cv[ii]);
453 /* Total energy */
454 cs_real_t v2 = cs_math_3_square_norm(vel[ii]);
455 ener[ii] = (pres[ii]+gamma[ii]*psginf) / ((gamma[ii]-1.)*dens[ii])
456 + 0.5*v2;
457 }
458
459 BFT_FREE(gamma);
460 }
461 }
462
463 /*----------------------------------------------------------------------------*/
464 /*!
465 * \brief Compute density and temperature from pressure and total energy;
466 *
467 * \param[in] cp array of isobaric specific heat values
468 * \param[in] cv array of isochoric specific heat values
469 * \param[in] pres array of pressure values
470 * \param[in] ener array of total energy values
471 * \param[out] dens array of density values
472 * \param[out] temp array of temperature values
473 * \param[in] vel array of velocity component values
474 * \param[in] l_size l_size of the array
475 */
476 /*----------------------------------------------------------------------------*/
477
478 void
cs_cf_thermo_dt_from_pe(cs_real_t * cp,cs_real_t * cv,cs_real_t * pres,cs_real_t * ener,cs_real_t * dens,cs_real_t * temp,cs_real_3_t * vel,cs_lnum_t l_size)479 cs_cf_thermo_dt_from_pe(cs_real_t *cp,
480 cs_real_t *cv,
481 cs_real_t *pres,
482 cs_real_t *ener,
483 cs_real_t *dens,
484 cs_real_t *temp,
485 cs_real_3_t *vel,
486 cs_lnum_t l_size)
487 {
488 /* Local variables */
489 cs_real_t enint;
490 int ieos = cs_glob_cf_model->ieos;
491
492 /* single ideal gas or stiffened gas eos - constant gamma */
493 if (ieos == CS_EOS_IDEAL_GAS || ieos == CS_EOS_STIFFENED_GAS) {
494 cs_real_t gamma0;
495 cs_real_t cp0 = cs_glob_fluid_properties->cp0;
496 cs_real_t cv0 = cs_glob_fluid_properties->cv0;
497 cs_real_t psginf = cs_glob_cf_model->psginf;
498 cs_lnum_t l_size0 = 1;
499
500 cs_cf_thermo_gamma(&cp0, &cv0, &gamma0, l_size0);
501
502 for (cs_lnum_t ii = 0; ii < l_size; ii++) {
503 /* Internal energy (to avoid the need to divide by the temperature
504 to compute density) */
505 cs_real_t v2 = cs_math_3_square_norm(vel[ii]);
506 enint = ener[ii] - 0.5*v2;
507
508 /* Density */
509 dens[ii] = (pres[ii]+gamma0*psginf) / ((gamma0-1.)*enint);
510 /* Temperature */
511 temp[ii] = (pres[ii]+psginf) / ((gamma0-1.)*dens[ii]*cv0);
512 }
513 }
514 /* ideal gas mixture */
515 else if (ieos == CS_EOS_GAS_MIX) {
516 cs_real_t *gamma;
517 cs_real_t psginf = cs_glob_cf_model->psginf;
518
519 BFT_MALLOC(gamma, l_size, cs_real_t);
520
521 cs_cf_thermo_gamma(cp, cv, gamma, l_size);
522
523 for (cs_lnum_t ii = 0; ii < l_size; ii++) {
524 /* Internal energy (to avoid the need to divide by the temperature
525 to compute density) */
526 cs_real_t v2 = cs_math_3_square_norm(vel[ii]);
527 enint = ener[ii] - 0.5*v2;
528
529 /* Density */
530 dens[ii] = (pres[ii]+gamma[ii]*psginf) / ((gamma[ii]-1.)*enint);
531 /* Temperature */
532 temp[ii] = (pres[ii]+psginf) / ((gamma[ii]-1.)*dens[ii]*cv[ii]);
533 }
534
535 BFT_FREE(gamma);
536 }
537 }
538
539 /*----------------------------------------------------------------------------*/
540 /*!
541 * \brief Compute pressure and total energy from density and temperature
542 *
543 * \param[in] cp array of isobaric specific heat values
544 * \param[in] cv array of isochoric specific heat values
545 * \param[in] dens array of density values
546 * \param[in] temp array of temperature values
547 * \param[out] pres array of pressure values
548 * \param[out] ener array of total energy values
549 * \param[in] vel array of velocity component values
550 * \param[in] l_size l_size of the array
551 */
552 /*----------------------------------------------------------------------------*/
553
554 void
cs_cf_thermo_pe_from_dt(cs_real_t * cp,cs_real_t * cv,cs_real_t * dens,cs_real_t * temp,cs_real_t * pres,cs_real_t * ener,cs_real_3_t * vel,cs_lnum_t l_size)555 cs_cf_thermo_pe_from_dt(cs_real_t *cp,
556 cs_real_t *cv,
557 cs_real_t *dens,
558 cs_real_t *temp,
559 cs_real_t *pres,
560 cs_real_t *ener,
561 cs_real_3_t *vel,
562 cs_lnum_t l_size)
563 {
564 /* Local variables */
565 int ieos = cs_glob_cf_model->ieos;
566
567 /* single ideal gas or stiffened gas eos - constant gamma */
568 if (ieos == CS_EOS_IDEAL_GAS || ieos == CS_EOS_STIFFENED_GAS) {
569 cs_real_t gamma0;
570 cs_real_t cp0 = cs_glob_fluid_properties->cp0;
571 cs_real_t cv0 = cs_glob_fluid_properties->cv0;
572 cs_real_t psginf = cs_glob_cf_model->psginf;
573 cs_lnum_t l_size0 = 1;
574
575 cs_cf_thermo_gamma(&cp0, &cv0, &gamma0, l_size0);
576
577 for (cs_lnum_t ii = 0; ii < l_size; ii++) {
578 /* Pressure */
579 pres[ii] = (gamma0-1.)*cv0*dens[ii]*temp[ii] - psginf;
580 /* Total energy */
581 cs_real_t v2 = cs_math_3_square_norm(vel[ii]);
582 ener[ii] = (pres[ii]+gamma0*psginf) / ((gamma0-1.)*dens[ii]) + 0.5*v2;
583 }
584 }
585 /* ideal gas mixture */
586 else if (ieos == CS_EOS_GAS_MIX) {
587 cs_real_t *gamma;
588 cs_real_t psginf = cs_glob_cf_model->psginf;
589 BFT_MALLOC(gamma, l_size, cs_real_t);
590
591 cs_cf_thermo_gamma(cp, cv, gamma, l_size);
592
593 for (cs_lnum_t ii = 0; ii < l_size; ii++) {
594 /* Pressure */
595 pres[ii] = (gamma[ii]-1.)*cv[ii]*dens[ii]*temp[ii] - psginf;
596 /* Total energy */
597 cs_real_t v2 = cs_math_3_square_norm(vel[ii]);
598 ener[ii] = (pres[ii]+gamma[ii]*psginf) / ((gamma[ii]-1.)*dens[ii]) + 0.5*v2;
599 }
600
601 BFT_FREE(gamma);
602 }
603 }
604
605 /*----------------------------------------------------------------------------*/
606 /*!
607 * \brief Compute pressure and temperature from density and total energy.
608 *
609 * \param[in] cp array of isobaric specific heat values
610 * \param[in] cv array of isochoric specific heat values
611 * \param[in] dens array of density values
612 * \param[in] ener array of total energy values
613 * \param[out] pres array of pressure values
614 * \param[out] temp array of temperature values
615 * \param[in] vel array of velocity component values
616 * \param[in,out] fracv array of volume fraction values
617 * \param[in,out] fracm array of mass fraction values
618 * \param[in,out] frace array of energy fraction values
619 * \param[in] l_size l_size of the array
620 */
621 /*----------------------------------------------------------------------------*/
622
623 void
cs_cf_thermo_pt_from_de(cs_real_t * cp,cs_real_t * cv,cs_real_t * dens,cs_real_t * ener,cs_real_t * pres,cs_real_t * temp,cs_real_3_t * vel,cs_real_t * fracv,cs_real_t * fracm,cs_real_t * frace,cs_lnum_t l_size)624 cs_cf_thermo_pt_from_de(cs_real_t *cp,
625 cs_real_t *cv,
626 cs_real_t *dens,
627 cs_real_t *ener,
628 cs_real_t *pres,
629 cs_real_t *temp,
630 cs_real_3_t *vel,
631 cs_real_t *fracv,
632 cs_real_t *fracm,
633 cs_real_t *frace,
634 cs_lnum_t l_size)
635 {
636 /* Local variables */
637 cs_real_t enint;
638 int ieos = cs_glob_cf_model->ieos;
639
640 /* single ideal gas or stiffened gas eos - constant gamma */
641 if (ieos == CS_EOS_IDEAL_GAS || ieos == CS_EOS_STIFFENED_GAS) {
642 cs_real_t gamma0;
643 cs_real_t cp0 = cs_glob_fluid_properties->cp0;
644 cs_real_t cv0 = cs_glob_fluid_properties->cv0;
645 cs_real_t psginf = cs_glob_cf_model->psginf;
646 cs_lnum_t l_size0 = 1;
647
648 cs_cf_thermo_gamma(&cp0, &cv0, &gamma0, l_size0);
649
650 for (cs_lnum_t ii = 0; ii < l_size; ii++) {
651 /* Internal energy (to avoid the need to divide by the temperature
652 to compute density) */
653 cs_real_t v2 = cs_math_3_square_norm(vel[ii]);
654 enint = ener[ii] - 0.5*v2;
655
656 /* Pressure */
657 pres[ii] = (gamma0-1.)*dens[ii]*enint - gamma0*psginf;
658 /* Temperature */
659 temp[ii] = (pres[ii]+psginf) / ((gamma0-1.)*dens[ii]*cv0);
660 }
661 }
662 /* ideal gas mixture */
663 else if (ieos == CS_EOS_GAS_MIX) {
664 cs_real_t *gamma;
665 cs_real_t psginf = cs_glob_cf_model->psginf;
666 BFT_MALLOC(gamma, l_size, cs_real_t);
667
668 cs_cf_thermo_gamma(cp, cv, gamma, l_size);
669
670 for (cs_lnum_t ii = 0; ii < l_size; ii++) {
671 /* Internal energy (to avoid the need to divide by the temperature
672 to compute density) */
673 cs_real_t v2 = cs_math_3_square_norm(vel[ii]);
674 enint = ener[ii] - 0.5*v2;
675
676 /* Pressure */
677 pres[ii] = (gamma[ii]-1.)*dens[ii]*enint - gamma[ii]*psginf;
678 /* Temperature */
679 temp[ii] = (pres[ii]+psginf) / ((gamma[ii]-1.)*dens[ii]*cv[ii]);
680 }
681
682 BFT_FREE(gamma);
683 }
684 /* homogeneous two phase */
685 else if (ieos == CS_EOS_HOMOGENEOUS_TWO_PHASE) {
686 for (cs_lnum_t ii = 0; ii < l_size; ii++) {
687 cs_real_t v2 = cs_math_3_square_norm(vel[ii]);
688
689 enint = ener[ii] - 0.5*v2;
690
691 cs_real_t tau = 1./dens[ii];
692
693 cs_hgn_thermo_pt(fracv[ii],
694 fracm[ii],
695 frace[ii],
696 enint,
697 tau,
698 &temp[ii],
699 &pres[ii]);
700 }
701 }
702 }
703
704 /*----------------------------------------------------------------------------*/
705 /*!
706 * \brief Compute square of sound velocity:
707 *
708 * \f[c^2 = \left(\frac{\partial p}{\partial \rho}\right)_s\f];
709 *
710 * for perfect gas, this expression simply writes:
711 *
712 * \f[c^2 = \gamma \frac{p}{\rho}\f]
713 *
714 * \param[in] cp array of isobaric specific heat values
715 * \param[in] cv array of isochoric specific heat values
716 * \param[in] pres array of pressure values
717 * \param[in] dens array of density values
718 * \param[in,out] fracv array of volume fraction values
719 * \param[in,out] fracm array of mass fraction values
720 * \param[in,out] frace array of energy fraction values
721 * \param[out] c2 array of the values of the square of sound velocity
722 * \param[in] l_size l_size of the array
723 */
724 /*----------------------------------------------------------------------------*/
725
726 void
cs_cf_thermo_c_square(cs_real_t * cp,cs_real_t * cv,cs_real_t * pres,cs_real_t * dens,cs_real_t * fracv,cs_real_t * fracm,cs_real_t * frace,cs_real_t * c2,cs_lnum_t l_size)727 cs_cf_thermo_c_square(cs_real_t *cp,
728 cs_real_t *cv,
729 cs_real_t *pres,
730 cs_real_t *dens,
731 cs_real_t *fracv,
732 cs_real_t *fracm,
733 cs_real_t *frace,
734 cs_real_t *c2,
735 cs_lnum_t l_size)
736 {
737 /* Local variables */
738 int ieos = cs_glob_cf_model->ieos;
739
740 /* single ideal gas or stiffened gas eos - constant gamma */
741 if (ieos == CS_EOS_IDEAL_GAS || ieos == CS_EOS_STIFFENED_GAS) {
742 cs_real_t gamma0;
743 cs_real_t cp0 = cs_glob_fluid_properties->cp0;
744 cs_real_t cv0 = cs_glob_fluid_properties->cv0;
745 cs_real_t psginf = cs_glob_cf_model->psginf;
746 cs_lnum_t l_size0 = 1;
747
748 cs_cf_thermo_gamma(&cp0, &cv0, &gamma0, l_size0);
749
750 for (cs_lnum_t ii = 0; ii < l_size; ii++)
751 c2[ii] = gamma0 * (pres[ii]+psginf) / dens[ii];
752 }
753 /* ideal gas mixture */
754 else if (ieos == CS_EOS_GAS_MIX) {
755 cs_real_t *gamma;
756 cs_real_t psginf = cs_glob_cf_model->psginf;
757 BFT_MALLOC(gamma, l_size, cs_real_t);
758
759 cs_cf_thermo_gamma(cp, cv, gamma, l_size);
760
761 for (cs_lnum_t ii = 0; ii < l_size; ii++)
762 c2[ii] = gamma[ii] * (pres[ii]+psginf) / dens[ii];
763
764 BFT_FREE(gamma);
765 }
766 else if (ieos == CS_EOS_HOMOGENEOUS_TWO_PHASE){
767 for (cs_lnum_t ii = 0; ii < l_size; ii++) {
768 cs_real_t tau = 1./dens[ii];
769
770 c2[ii] = cs_hgn_thermo_c2(fracv[ii],
771 fracm[ii],
772 frace[ii],
773 pres[ii],
774 tau);
775 }
776 }
777 }
778
779 /*----------------------------------------------------------------------------*/
780 /*!
781 * \brief Compute the thermal expansion coefficient:
782 *
783 * \f[ \beta = \left(\frac{\partial p}{\partial s}\right)_\rho \f]
784 *
785 * for a perfect gas, the explicit formula is:
786 *
787 * \f[ \beta = \rho^\gamma \f]
788 *
789 * \param[in] cp array of isobaric specific heat values
790 * \param[in] cv array of isochoric specific heat values
791 * \param[in] dens array of density values
792 * \param[out] beta array of beta values
793 * \param[in] l_size l_size of the array
794 */
795 /*----------------------------------------------------------------------------*/
796
797 void
cs_cf_thermo_beta(cs_real_t * cp,cs_real_t * cv,cs_real_t * dens,cs_real_t * beta,cs_lnum_t l_size)798 cs_cf_thermo_beta(cs_real_t *cp,
799 cs_real_t *cv,
800 cs_real_t *dens,
801 cs_real_t *beta,
802 cs_lnum_t l_size)
803 {
804 /* Local variables */
805 int ieos = cs_glob_cf_model->ieos;
806
807 /* single ideal gas or stiffened gas eos - constant gamma */
808 if (ieos == CS_EOS_IDEAL_GAS || ieos == CS_EOS_STIFFENED_GAS) {
809 cs_real_t gamma0;
810 cs_real_t cp0 = cs_glob_fluid_properties->cp0;
811 cs_real_t cv0 = cs_glob_fluid_properties->cv0;
812 cs_lnum_t l_size0 = 1;
813
814 cs_cf_thermo_gamma(&cp0, &cv0, &gamma0, l_size0);
815
816 for (cs_lnum_t ii = 0; ii < l_size; ii++)
817 beta[ii] = pow(dens[ii],gamma0);
818 }
819 /* ideal gas mixture */
820 else if (ieos == CS_EOS_GAS_MIX) {
821 cs_real_t *gamma;
822 BFT_MALLOC(gamma, l_size, cs_real_t);
823
824 cs_cf_thermo_gamma(cp, cv, gamma, l_size);
825
826 for (cs_lnum_t ii = 0; ii < l_size; ii++)
827 beta[ii] = pow(dens[ii],gamma[ii]);
828
829 BFT_FREE(gamma);
830 }
831 }
832
833 /*----------------------------------------------------------------------------*/
834 /*!
835 * \brief Compute the isochoric specific heat:
836 *
837 * \f[C_v = \left(\frac{\partial e}{\partial T}\right)_\rho\f]
838 *
839 * \param[in] cp array of isobaric specific heat values
840 * \param[in] xmasml array of molar mass values
841 * \param[out] cv array of isochoric specific heat values
842 * \param[in] l_size l_size of the array
843 */
844 /*----------------------------------------------------------------------------*/
845
846 void
cs_cf_thermo_cv(cs_real_t * cp,cs_real_t * xmasml,cs_real_t * cv,cs_lnum_t l_size)847 cs_cf_thermo_cv(cs_real_t *cp,
848 cs_real_t *xmasml,
849 cs_real_t *cv,
850 cs_lnum_t l_size)
851 {
852 int ieos = cs_glob_cf_model->ieos;
853
854 /* Cv for a single ideal gas or a mixture of ideal gas */
855 if (ieos == CS_EOS_IDEAL_GAS || ieos == CS_EOS_GAS_MIX) {
856 cs_real_t r_pg = cs_physical_constants_r;
857 for (cs_lnum_t ii = 0; ii < l_size; ii++)
858 cv[ii] = cp[ii]-r_pg/xmasml[ii];
859 }
860 /* Cv for a stiffened gas */
861 else if (ieos == CS_EOS_STIFFENED_GAS) {
862 for (cs_lnum_t ii = 0; ii < l_size; ii++)
863 cv[ii] = cs_glob_fluid_properties->cv0;
864 }
865 }
866
867 /*----------------------------------------------------------------------------*/
868 /*!
869 * \brief Compute entropy from pressure and density:
870 *
871 * \f[s = \frac{p}{\rho^\gamma}\f]
872 *
873 * \param[in] cp array of isobaric specific heat values
874 * \param[in] cv array of isochoric specific heat values
875 * \param[in] dens array of density values
876 * \param[in] pres array of pressure values
877 * \param[out] entr array of total energy values
878 * \param[in] l_size l_size of the array
879 */
880 /*----------------------------------------------------------------------------*/
881
882 void
cs_cf_thermo_s_from_dp(cs_real_t * cp,cs_real_t * cv,cs_real_t * dens,cs_real_t * pres,cs_real_t * entr,cs_lnum_t l_size)883 cs_cf_thermo_s_from_dp(cs_real_t *cp,
884 cs_real_t *cv,
885 cs_real_t *dens,
886 cs_real_t *pres,
887 cs_real_t *entr,
888 cs_lnum_t l_size)
889 {
890 /* Local variables */
891 int ieos = cs_glob_cf_model->ieos;
892
893 /* single ideal gas or stiffened gas eos - constant gamma */
894 if (ieos == CS_EOS_IDEAL_GAS || ieos == CS_EOS_STIFFENED_GAS) {
895 cs_real_t gamma0;
896 cs_real_t cp0 = cs_glob_fluid_properties->cp0;
897 cs_real_t cv0 = cs_glob_fluid_properties->cv0;
898 cs_real_t psginf = cs_glob_cf_model->psginf;
899 cs_lnum_t l_size0 = 1;
900
901 cs_cf_thermo_gamma(&cp0, &cv0, &gamma0, l_size0);
902
903 cs_cf_check_density(dens, l_size0);
904
905 for (cs_lnum_t ii = 0; ii < l_size; ii++)
906 entr[ii] = (pres[ii]+psginf) / pow(dens[ii],gamma0);
907 }
908 /* ideal gas mixture */
909 else if (ieos == CS_EOS_GAS_MIX) {
910 cs_real_t *gamma;
911 cs_real_t psginf = cs_glob_cf_model->psginf;
912
913 BFT_MALLOC(gamma, l_size, cs_real_t);
914
915 cs_cf_thermo_gamma(cp, cv, gamma, l_size);
916
917 cs_cf_check_density(dens, l_size);
918
919 for (cs_lnum_t ii = 0; ii < l_size; ii++)
920 entr[ii] = (pres[ii]+psginf) / pow(dens[ii],gamma[ii]);
921
922 BFT_FREE(gamma);
923 }
924 }
925
926 /*----------------------------------------------------------------------------*/
927 /*!
928 * \brief Compute wall boundary condition values.
929 *
930 * \param[in,out] wbfa output work array
931 * \param[in,out] wbfb output work array
932 * \param[in] face_id boundary face index
933 */
934 /*----------------------------------------------------------------------------*/
935
936 void
cs_cf_thermo_wall_bc(cs_real_t * wbfa,cs_real_t * wbfb,cs_lnum_t face_id)937 cs_cf_thermo_wall_bc(cs_real_t *wbfa,
938 cs_real_t *wbfb,
939 cs_lnum_t face_id)
940 {
941 const cs_mesh_t *m = cs_glob_mesh;
942 cs_mesh_quantities_t *fvq = cs_glob_mesh_quantities;
943
944 const cs_lnum_t *restrict b_face_cells
945 = (const cs_lnum_t *restrict)m->b_face_cells;
946 const cs_real_3_t *restrict b_face_normal
947 = (const cs_real_3_t *restrict)fvq->b_face_normal;
948 const cs_real_t *restrict b_face_surf = fvq->b_face_surf;
949
950 /* Map field arrays */
951 cs_real_3_t *vel = (cs_real_3_t *)CS_F_(vel)->val;
952 cs_real_t *cvar_pr = CS_F_(p)->val;
953 cs_real_t *crom = CS_F_(rho)->val;
954
955 cs_real_t cp, cv, gamma;
956 cs_lnum_t l_size = 1;
957 int ieos = cs_glob_cf_model->ieos;
958
959 /* single ideal gas or stiffened gas eos or ideal gas mixture */
960 if (ieos == CS_EOS_IDEAL_GAS || ieos == CS_EOS_STIFFENED_GAS || ieos == CS_EOS_GAS_MIX) {
961 cs_real_t psginf = cs_glob_cf_model->psginf;
962 cs_lnum_t cell_id = b_face_cells[face_id];
963
964 if (ieos == CS_EOS_IDEAL_GAS || ieos == CS_EOS_STIFFENED_GAS) {
965 cp = cs_glob_fluid_properties->cp0;
966 cv = cs_glob_fluid_properties->cv0;
967 }
968 else if (ieos == CS_EOS_GAS_MIX) {
969 cp = CS_F_(cp)->val[cell_id];
970 cv = CS_F_(cv)->val[cell_id];
971 }
972
973 cs_cf_thermo_gamma(&cp, &cv, &gamma, l_size);
974
975 /* Calculation of the Mach number at the boundary face, using the
976 cell center velocity projected on the vector normal to the boundary */
977 cs_real_t uni = cs_math_3_dot_product(vel[cell_id],b_face_normal[face_id])
978 / b_face_surf[face_id];
979 cs_real_t xmach = uni
980 / sqrt(gamma * (cvar_pr[cell_id]+psginf) / crom[cell_id]);
981
982 /* Pressure */
983
984 /* A Neumann boundary condition is used. This does not allow to use
985 the Rusanov scheme, but some stabilization effect is expected.
986 A test based on the value of coefb at the previous time step
987 is implemented to avoid oscillating between a rarefaction
988 situation and a shock configuration from one time step to the
989 next. */
990
991 /* Rarefaction */
992 /* Here wbfb is the ratio (pwall+psginf)/(pi+psginf)
993 at the previous time step */
994 if (xmach < 0. && wbfb[face_id] <= 1.) {
995
996 if (xmach > 2./(1.-gamma))
997 wbfb[face_id] = pow(1. + (gamma-1.)/2. * xmach, 2.*gamma/(gamma-1.));
998 else
999 /* In case the rarefaction is too strong, a zero Dirichlet value
1000 is used for pressure (the value of wbfb is used here as an
1001 indicator) */
1002 wbfb[face_id] = cs_math_infinite_r;
1003
1004 }
1005 /* Shock */
1006 else if (xmach > 0. && wbfb[face_id] >= 1.) {
1007
1008 wbfb[face_id] = 1. + gamma*xmach
1009 *((gamma+1.)/4.*xmach
1010 + sqrt(1. + pow(gamma+1.,2)/16.*xmach*xmach));
1011
1012 }
1013 /* Oscillation between rarefaction and shock or zero Mach number */
1014 else {
1015 wbfb[face_id] = 1.;
1016 }
1017
1018 /* Here wbfb is the ratio of wall pressure over cell pressure for
1019 "the perfect gas part" of the wall pressure and wbfa is the
1020 "purely stiffened gas part" of the wall pressure
1021 (wbfa=0 when in perfect gas), but wbfa is not a pressure ratio:
1022 pwall = wbfb[face_id]*cvar_pr[cell_id]+wbfa[face_id] */
1023 wbfa[face_id] = (wbfb[face_id]-1.)*psginf;
1024 }
1025
1026 }
1027
1028 /*----------------------------------------------------------------------------*/
1029 /*!
1030 * \brief Compute subsonic outlet boundary conditions values.
1031
1032 * \param[in,out] bc_en total energy values at boundary faces
1033 * \param[in,out] bc_pr pressure values at boundary faces
1034 * \param[in,out] bc_vel velocity values at boundary faces
1035 * \param[in] face_id boundary face index
1036 */
1037 /*----------------------------------------------------------------------------*/
1038
1039 void
cs_cf_thermo_subsonic_outlet_bc(cs_real_t * bc_en,cs_real_t * bc_pr,cs_real_3_t * bc_vel,cs_lnum_t face_id)1040 cs_cf_thermo_subsonic_outlet_bc(cs_real_t *bc_en,
1041 cs_real_t *bc_pr,
1042 cs_real_3_t *bc_vel,
1043 cs_lnum_t face_id)
1044 {
1045 const cs_mesh_t *m = cs_glob_mesh;
1046 cs_mesh_quantities_t *fvq = cs_glob_mesh_quantities;
1047
1048 const cs_lnum_t *restrict b_face_cells
1049 = (const cs_lnum_t *restrict)m->b_face_cells;
1050 const cs_real_3_t *restrict b_face_normal
1051 = (const cs_real_3_t *restrict)fvq->b_face_normal;
1052 const cs_real_t *restrict b_face_surf = fvq->b_face_surf;
1053
1054 cs_real_t gamma, yp;
1055 cs_real_t roi, ro1, pri, uni, un1, uns;
1056 cs_real_t ci, c1, mi, a, b, sigma1, pinf;
1057
1058 /* Map field arrays */
1059 cs_real_3_t *vel = (cs_real_3_t *)CS_F_(vel)->val;
1060 cs_real_t *cvar_pr = CS_F_(p)->val;
1061 cs_real_t *cvar_en = CS_F_(e_tot)->val;
1062 cs_real_t *crom = CS_F_(rho)->val;
1063 cs_real_t *brom = CS_F_(rho_b)->val;
1064
1065 cs_real_t cp, cv;
1066 cs_lnum_t l_size = 1;
1067 int ieos = cs_glob_cf_model->ieos;
1068
1069 /* single ideal gas or stiffened gas eos or ideal gas mixture */
1070 if (ieos == CS_EOS_IDEAL_GAS || ieos == CS_EOS_STIFFENED_GAS || ieos == CS_EOS_GAS_MIX) {
1071 cs_real_t psginf = cs_glob_cf_model->psginf;
1072 cs_lnum_t cell_id = b_face_cells[face_id];
1073
1074 if (ieos == CS_EOS_IDEAL_GAS || ieos == CS_EOS_STIFFENED_GAS) {
1075 cp = cs_glob_fluid_properties->cp0;
1076 cv = cs_glob_fluid_properties->cv0;
1077 }
1078 else if (ieos == CS_EOS_GAS_MIX) {
1079 cp = CS_F_(cp)->val[cell_id];
1080 cv = CS_F_(cv)->val[cell_id];
1081 }
1082
1083 cs_cf_thermo_gamma(&cp, &cv, &gamma, l_size);
1084
1085 pinf = bc_pr[face_id];
1086 pri = cvar_pr[cell_id];
1087 yp = (pinf+psginf) / (pri+psginf);
1088 roi = crom[cell_id];
1089
1090 ci = sqrt(gamma * pri / roi);
1091 uni = cs_math_3_dot_product(vel[cell_id],b_face_normal[face_id])
1092 /b_face_surf[face_id];
1093
1094 cs_real_t deltap = pinf-pri;
1095 cs_real_t res = CS_ABS(deltap/(pinf+psginf));
1096 /* Rarefaction case */
1097 if (deltap < 0. || res < cs_math_epzero) {
1098
1099 /* Computation of the velocity in state 1 using Riemann invariants
1100 of the 1-rarefaction */
1101 a = 2 * ci / (gamma - 1.)
1102 * (1. - pow(yp, (gamma-1.)/(2.*gamma)));
1103 un1 = uni + a;
1104
1105 /* Computation of the density in state 1 using Rieman invariants
1106 of the 1-rarefaction */
1107 ro1 = roi * pow(yp, 1./gamma);
1108
1109 /* Subsonic inlet - state 2 should be imposed but too few information
1110 is available to compute it
1111 for want of anything better, state 1 is imposed */
1112 if (un1 < 0.) {
1113
1114 /* Density */
1115 brom[face_id] = ro1;
1116 /* Velocity */
1117 bc_vel[face_id][0] = vel[cell_id][0]
1118 + a * b_face_normal[face_id][0] / b_face_surf[face_id];
1119 bc_vel[face_id][1] = vel[cell_id][1]
1120 + a * b_face_normal[face_id][1] / b_face_surf[face_id];
1121 bc_vel[face_id][2] = vel[cell_id][2]
1122 + a * b_face_normal[face_id][2] / b_face_surf[face_id];
1123 /* Total energy */
1124 bc_en[face_id] = (pinf+gamma*psginf) / ((gamma - 1.) * ro1)
1125 + 0.5 * cs_math_3_square_norm(bc_vel[face_id]);
1126
1127 }
1128 /* Outlet */
1129 else {
1130
1131 /* Computation of the sound speed in state 1 */
1132 c1 = sqrt(gamma * (pinf+psginf) / ro1);
1133
1134 /* Subsonic outlet - state 1 is imposed */
1135 if ((un1-c1) < 0.) {
1136
1137 /* Density */
1138 brom[face_id] = ro1;
1139 /* Velocity */
1140 bc_vel[face_id][0] = vel[cell_id][0]
1141 + a * b_face_normal[face_id][0] / b_face_surf[face_id];
1142 bc_vel[face_id][1] = vel[cell_id][1]
1143 + a * b_face_normal[face_id][1] / b_face_surf[face_id];
1144 bc_vel[face_id][2] = vel[cell_id][2]
1145 + a * b_face_normal[face_id][2] / b_face_surf[face_id];
1146 /* Total energy */
1147 bc_en[face_id] = (pinf+gamma*psginf) / ((gamma - 1.) * ro1)
1148 + 0.5 * cs_math_3_square_norm(bc_vel[face_id]);
1149
1150 }
1151 /* Sonic outlet */
1152 else if ((uni-ci) < 0.) {
1153
1154 /* Mach number in the domain */
1155 mi = uni / ci;
1156
1157 b = (gamma - 1.) / (gamma + 1.) * (mi + 2. / (gamma - 1));
1158
1159 /* Sonic state pressure */
1160 bc_pr[face_id] = -psginf + (pri+psginf) * pow(b, 2. * gamma / (gamma - 1.));
1161 /* Sonic state density */
1162 brom[face_id] = roi * pow(b, 2. / (gamma - 1.));
1163 /* Sonic state velocity */
1164 uns = b * ci;
1165 bc_vel[face_id][0] = uns * b_face_normal[face_id][0] / b_face_surf[face_id];
1166 bc_vel[face_id][1] = uns * b_face_normal[face_id][1] / b_face_surf[face_id];
1167 bc_vel[face_id][2] = uns * b_face_normal[face_id][2] / b_face_surf[face_id];
1168 /* Sonic state energy */
1169 bc_en[face_id] = (bc_pr[face_id]+gamma*psginf)
1170 /((gamma - 1.)*brom[face_id])
1171 + 0.5 * uns*uns;
1172
1173 }
1174 /* Supersonic outlet */
1175 else {
1176
1177 /* pb = pri */
1178 bc_pr[face_id] = pri;
1179 /* ub = uni */
1180 bc_vel[face_id][0] = vel[cell_id][0];
1181 bc_vel[face_id][1] = vel[cell_id][1];
1182 bc_vel[face_id][2] = vel[cell_id][2];
1183 /* rob = roi */
1184 brom[face_id] = roi;
1185 /* eb = ei */
1186 bc_en[face_id] = cvar_en[cell_id];
1187
1188 }
1189
1190
1191 }
1192
1193 }
1194 /* Shock case */
1195 else {
1196
1197 /* Computation of the density in state 1 with Rankine-Hugoniot relations */
1198 ro1 = roi * ((gamma - 1.) * (pri+psginf) + (gamma + 1.) * (pinf+psginf))
1199 / ((gamma - 1.) * (pinf+psginf) + (gamma + 1.) * (pri+psginf));
1200
1201 /* Computation of the velocity in state 1 with Rankine-Hugoniot relations
1202 un1 = un2 */
1203 a = sqrt( (pinf - pri) * (1./roi - 1./ro1) );
1204 un1 = uni - a;
1205
1206 /* Subsonic inlet - state 2 should be imposed but too few information
1207 is available to compute it
1208 for want of anything better, state 1 is imposed */
1209 if (un1 <= 0.) {
1210
1211 /* Density */
1212 brom[face_id] = ro1;
1213 /* Velocity */
1214 bc_vel[face_id][0] = vel[cell_id][0] - a * b_face_normal[face_id][0]
1215 / b_face_surf[face_id];
1216 bc_vel[face_id][1] = vel[cell_id][1] - a * b_face_normal[face_id][1]
1217 / b_face_surf[face_id];
1218 bc_vel[face_id][2] = vel[cell_id][2] - a * b_face_normal[face_id][2]
1219 / b_face_surf[face_id];
1220 /* Total energy */
1221 bc_en[face_id] = (pinf+gamma*psginf) / ((gamma-1.) * brom[face_id])
1222 + 0.5 * cs_math_3_square_norm(bc_vel[face_id]);
1223
1224 }
1225 /* Outlet */
1226 else {
1227
1228 /* Computation of the shock velocity */
1229 sigma1 = (roi * uni - ro1 * un1) / (roi - ro1);
1230
1231 /* Subsonic outlet - state 1 is imposed */
1232 if (sigma1 <= 0.) {
1233
1234 /* Density */
1235 brom[face_id] = ro1;
1236 /* Velocity */
1237 bc_vel[face_id][0] = vel[cell_id][0]
1238 - a * b_face_normal[face_id][0] / b_face_surf[face_id];
1239 bc_vel[face_id][1] = vel[cell_id][1]
1240 - a * b_face_normal[face_id][1] / b_face_surf[face_id];
1241 bc_vel[face_id][2] = vel[cell_id][2]
1242 - a * b_face_normal[face_id][2] / b_face_surf[face_id];
1243 /* Total energy */
1244 bc_en[face_id] = (pinf+gamma*psginf) / ((gamma-1.) * brom[face_id])
1245 + 0.5 * cs_math_3_square_norm(bc_vel[face_id]);
1246
1247 }
1248 /* Supersonic outlet */
1249 else {
1250
1251 /* pb = pri */
1252 bc_pr[face_id] = pri;
1253 /* unb = uni */
1254 bc_vel[face_id][0] = vel[cell_id][0];
1255 bc_vel[face_id][1] = vel[cell_id][1];
1256 bc_vel[face_id][2] = vel[cell_id][2];
1257 /* rob = roi */
1258 brom[face_id] = roi;
1259 /* eb = ei */
1260 bc_en[face_id] = cvar_en[cell_id];
1261
1262 } /* test on shock speed sign */
1263
1264 } /* test on state 1 velocity sign */
1265
1266 } /* test on pinf-pri sign */
1267
1268 }
1269
1270 }
1271
1272 /*----------------------------------------------------------------------------*/
1273 /*!
1274 * \brief Compute inlet boundary condition with total pressure and total
1275 * enthalpy imposed.
1276 *
1277 * \param[in,out] bc_en total energy values at boundary faces
1278 * \param[in,out] bc_pr pressure values at boundary faces
1279 * \param[in,out] bc_vel velocity values at boundary faces
1280 * \param[in] face_id boundary face index
1281 */
1282 /*----------------------------------------------------------------------------*/
1283
1284 void
cs_cf_thermo_ph_inlet_bc(cs_real_t * bc_en,cs_real_t * bc_pr,cs_real_3_t * bc_vel,cs_lnum_t face_id)1285 cs_cf_thermo_ph_inlet_bc(cs_real_t *bc_en,
1286 cs_real_t *bc_pr,
1287 cs_real_3_t *bc_vel,
1288 cs_lnum_t face_id)
1289 {
1290 const cs_mesh_t *m = cs_glob_mesh;
1291 cs_mesh_quantities_t *fvq = cs_glob_mesh_quantities;
1292
1293 const cs_lnum_t *restrict b_face_cells
1294 = (const cs_lnum_t *restrict)m->b_face_cells;
1295 const cs_real_3_t *restrict b_face_normal
1296 = (const cs_real_3_t *restrict)fvq->b_face_normal;
1297 const cs_real_t *restrict b_face_surf = fvq->b_face_surf;
1298
1299 int niter, nitermax;
1300 cs_real_t gamma, bMach, eps, pstat, old_pstat, ptot, res, rhotot;
1301 cs_real_t roi, ro1, pri, ei, uni, un1, y, uns, bc, cosalp, norm;
1302 cs_real_t ci, c1, mi, a, sigma1, utxi, utyi, utzi;
1303 cs_real_3_t dir;
1304
1305 /* Map field arrays */
1306 cs_real_3_t *vel = (cs_real_3_t *)CS_F_(vel)->val;
1307 cs_real_t *cvar_pr = CS_F_(p)->val;
1308 cs_real_t *cvar_en = CS_F_(e_tot)->val;
1309 cs_real_t *crom = CS_F_(rho)->val;
1310 cs_real_t *brom = CS_F_(rho_b)->val;
1311
1312 cs_real_t cp, cv;
1313 cs_lnum_t l_size = 1;
1314 int ieos = cs_glob_cf_model->ieos;
1315
1316 /* single ideal gas or stiffened gas eos or ideal gas mixture */
1317 if ( ieos == CS_EOS_IDEAL_GAS
1318 || ieos == CS_EOS_STIFFENED_GAS
1319 || ieos == CS_EOS_GAS_MIX) {
1320 cs_real_t psginf = cs_glob_cf_model->psginf;
1321 cs_lnum_t cell_id = b_face_cells[face_id];
1322
1323 if (ieos == CS_EOS_IDEAL_GAS || ieos == CS_EOS_STIFFENED_GAS) {
1324 cp = cs_glob_fluid_properties->cp0;
1325 cv = cs_glob_fluid_properties->cv0;
1326 }
1327 else if (ieos == CS_EOS_GAS_MIX) {
1328 cp = CS_F_(cp)->val[cell_id];
1329 cv = CS_F_(cv)->val[cell_id];
1330 }
1331
1332 cs_cf_thermo_gamma(&cp, &cv, &gamma, l_size);
1333
1334 niter = 0;
1335
1336 roi = crom[cell_id];
1337 pri = cvar_pr[cell_id];
1338
1339 /* Normalize the direction vector given by the user */
1340 norm = sqrt(cs_math_3_square_norm(bc_vel[face_id]));
1341 if (norm < cs_math_epzero)
1342 bft_error(__FILE__, __LINE__, 0,
1343 _("Error in thermodynamics computations for compressible "
1344 "flows:\n"
1345 "The computation of the subsonic inlet boundary condition\n"
1346 "with imposed total pressure and total enthalpy failed at\n"
1347 "boundary face %ld. The direction vector given by the user\n"
1348 "can't be null."),
1349 (long)face_id);
1350
1351 dir[0] = bc_vel[face_id][0] / norm;
1352 dir[1] = bc_vel[face_id][1] / norm;
1353 dir[2] = bc_vel[face_id][2] / norm;
1354
1355 /* Angle between the imposed direction and the inlet normal */
1356 cosalp = cs_math_3_dot_product(dir,b_face_normal[face_id])
1357 / b_face_surf[face_id];
1358
1359 /* If direction vector is outward, warn the user */
1360 if (cosalp > cs_math_epzero)
1361 bft_printf("Warning in thermodynamics computations for compressible"
1362 "flows:\n"
1363 "The computation of the subsonic inlet boundary condition\n"
1364 "with imposed total pressure and total enthalpy failed at\n"
1365 "boundary face %ld. The direction vector given by the user\n"
1366 "points outward the fluid domain.\n",
1367 (long)face_id);
1368
1369 /* Computation of the sound speed inside the domain */
1370 ci = sqrt(gamma * (pri+psginf) / roi);
1371
1372 uni = cs_math_3_dot_product(vel[cell_id],b_face_normal[face_id])
1373 / b_face_surf[face_id];
1374
1375 bMach = uni / ci;
1376
1377 utxi = vel[cell_id][0] - uni * b_face_normal[face_id][0]
1378 * b_face_surf[face_id];
1379 utyi = vel[cell_id][1] - uni * b_face_normal[face_id][1]
1380 * b_face_surf[face_id];
1381 utzi = vel[cell_id][2] - uni * b_face_normal[face_id][2]
1382 * b_face_surf[face_id];
1383
1384 cs_real_t v2 = cs_math_3_square_norm(vel[cell_id]);
1385 ei = cvar_en[cell_id] - 0.5 * v2;
1386
1387 ptot = bc_pr[face_id];
1388 /* bc_en holds the value of the total enthalpy given by the user */
1389 rhotot = gamma / (gamma - 1.) * (ptot+gamma*psginf) / bc_en[face_id];
1390 old_pstat = ptot;
1391
1392 int key_cal_opt_id = cs_field_key_id("var_cal_opt");
1393 cs_var_cal_opt_t var_cal_opt;
1394 cs_field_get_key_struct(CS_F_(p), key_cal_opt_id, &var_cal_opt);
1395
1396 eps = var_cal_opt.epsrsm;
1397 nitermax = 100;
1398 res = 1.;
1399
1400 while (niter <= nitermax && res > eps) {
1401
1402 pstat = -psginf + (ptot+psginf)
1403 * pow(1.+(gamma-1.)*0.5*bMach*bMach,gamma/(1.-gamma));
1404 y = pri / pstat;
1405
1406 /* 1-shock */
1407 if (y < 1.) {
1408
1409 /* Computation of the density in state 1 with Rankine-Hugoniot relations */
1410 ro1 = roi * ( (gamma-1.)*(pri+psginf )
1411 + (gamma+1.)*(pstat+psginf))
1412 / ( (gamma-1.)*(pstat+psginf)
1413 + (gamma+1.)*(pri+psginf ));
1414
1415 /* Computation of the velocity in state 1 with Rankine-Hugoniot relations
1416 un1 = un2 */
1417 un1 = uni - sqrt( (pstat - pri) * (1./roi - 1./ro1) );
1418
1419 /* Subsonic inlet */
1420 if (un1 <= 0.) {
1421
1422 /* unb = u2 */
1423 bc_vel[face_id][0] = un1 / cosalp * dir[0];
1424 bc_vel[face_id][1] = un1 / cosalp * dir[1];
1425 bc_vel[face_id][2] = un1 / cosalp * dir[2];
1426 /* rob = ro2 */
1427 brom[face_id] = pow((pstat+psginf)/(ptot+psginf),1./gamma) * rhotot;
1428 /* eb = e2 */
1429 bc_en[face_id] = (pstat+gamma*psginf) / ((gamma-1.)*brom[face_id])
1430 + 0.5 * cs_math_3_square_norm(bc_vel[face_id]);
1431 }
1432 /* Outlet */
1433 else {
1434 /* Computation of the shock velocity */
1435 sigma1 = (roi * uni - ro1 * un1) / (roi - ro1);
1436
1437 /* subsonic outlet */
1438 if (sigma1 <= 0.) {
1439
1440 /* unb = u1 */
1441 bc_vel[face_id][0] = utxi + un1 * b_face_normal[face_id][0]
1442 / b_face_surf[face_id];
1443 bc_vel[face_id][1] = utyi + un1 * b_face_normal[face_id][1]
1444 / b_face_surf[face_id];
1445 bc_vel[face_id][2] = utzi + un1 * b_face_normal[face_id][2]
1446 / b_face_surf[face_id];
1447 /* rob = ro1 */
1448 brom[face_id] = ro1;
1449 /* eb = e1 */
1450 bc_en[face_id] = ei
1451 - 0.5 * (pstat + pri) * (1. / ro1 - 1. / roi)
1452 + 0.5 * (un1*un1+utxi*utxi+utyi*utyi+utzi*utzi);
1453
1454 }
1455 /* supersonic outlet */
1456 else {
1457
1458 /* pb = pri */
1459 pstat = pri;
1460 /* unb = uni */
1461 bc_vel[face_id][0] = vel[cell_id][0];
1462 bc_vel[face_id][1] = vel[cell_id][1];
1463 bc_vel[face_id][2] = vel[cell_id][2];
1464 /* rob = roi */
1465 brom[face_id] = roi;
1466 /* eb = ei */
1467 bc_en[face_id] = cvar_en[cell_id];
1468
1469 }
1470
1471 }
1472
1473 }
1474 /* 1-rarefaction */
1475 else {
1476
1477 cs_real_t yp = (pstat+psginf) / (pri+psginf);
1478 /* Computation of the velocity in state 1 using Riemann invariants
1479 of the 1-rarefaction */
1480 un1 = uni + 2 * ci / (gamma - 1.)
1481 * (1. - pow(yp, (gamma-1.) / (2.*gamma)));
1482
1483 /* Computation of the density in state 1 using Riemann invariants
1484 of the 1-rarefaction */
1485 ro1 = pow(yp, 1. / gamma) * roi;
1486
1487 /* Subsonic inlet */
1488 if (un1 <= 0.) {
1489
1490 /* unb = u2 */
1491 bc_vel[face_id][0] = un1 / cosalp * dir[0];
1492 bc_vel[face_id][1] = un1 / cosalp * dir[1];
1493 bc_vel[face_id][2] = un1 / cosalp * dir[2];
1494 /* rob = ro2 */
1495 brom[face_id] = pow((pstat+psginf)/(ptot+psginf),1./gamma) * rhotot;
1496 /* eb = e2 */
1497 bc_en[face_id] = (pstat+gamma*psginf) / ((gamma-1.)*brom[face_id])
1498 + 0.5 * cs_math_3_square_norm(bc_vel[face_id]);
1499 }
1500 /* Outlet */
1501 else {
1502
1503 /* Computation of the sound speed in state 1 */
1504 c1 = sqrt(gamma * (pstat+psginf) / ro1);
1505
1506 /* Subsonic outlet */
1507 if ((un1 - c1) < 0.) {
1508
1509 /* unb = u1 */
1510 bc_vel[face_id][0] = utxi + un1 * b_face_normal[face_id][0]
1511 / b_face_surf[face_id];
1512 bc_vel[face_id][1] = utyi + un1 * b_face_normal[face_id][1]
1513 / b_face_surf[face_id];
1514 bc_vel[face_id][2] = utzi + un1 * b_face_normal[face_id][2]
1515 / b_face_surf[face_id];
1516 /* rob = ro1 */
1517 brom[face_id] = ro1;
1518 /* eb = e1 */
1519 bc_en[face_id] = (pstat+gamma*psginf) / (ro1 * (gamma-1.))
1520 + 0.5 * (un1*un1+utxi*utxi+utyi*utyi+utzi*utzi);
1521
1522 }
1523 /* Supersonic outlet */
1524 else if ((uni - ci) >= 0.) {
1525
1526 /* pb = pri */
1527 pstat = pri;
1528 /* ub = uni */
1529 bc_vel[face_id][0] = vel[cell_id][0];
1530 bc_vel[face_id][1] = vel[cell_id][1];
1531 bc_vel[face_id][2] = vel[cell_id][2];
1532 /* rob = roi */
1533 brom[face_id] = roi;
1534 /* eb = ei */
1535 bc_en[face_id] = cvar_en[cell_id];
1536
1537 }
1538 /* Outlet in sonic state */
1539 else {
1540
1541 /* Mach number in the domain */
1542 mi = uni / ci;
1543
1544 a = (gamma - 1.) / (gamma + 1.) * (mi + 2. / (gamma - 1));
1545
1546 /* Sonic state pressure */
1547 pstat = -psginf + (pri+psginf) * pow(a,2.*gamma/(gamma-1.));
1548 /* Sonic state density */
1549 brom[face_id] = roi * pow(a,2./(gamma-1.));
1550 /* Sonic state velocity */
1551 uns = a * ci;
1552 bc_vel[face_id][0] = uns * b_face_normal[face_id][0]
1553 / b_face_surf[face_id];
1554 bc_vel[face_id][1] = uns * b_face_normal[face_id][1]
1555 / b_face_surf[face_id];
1556 bc_vel[face_id][2] = uns * b_face_normal[face_id][2]
1557 / b_face_surf[face_id];
1558 /* Sonic state energy */
1559 bc_en[face_id] = (pstat+gamma*psginf) / ((gamma-1.)*brom[face_id])
1560 + 0.5 * uns*uns;
1561
1562 }
1563
1564 }
1565
1566 }
1567
1568
1569 bc = sqrt(gamma * (pstat+psginf) / brom[face_id]);
1570 bMach = cs_math_3_dot_product(bc_vel[face_id],b_face_normal[face_id])
1571 / b_face_surf[face_id] / bc;
1572
1573 bc_pr[face_id] = pstat;
1574
1575 /* Pressure residual */
1576 res = CS_ABS((pstat - old_pstat) / ptot);
1577
1578 /* Prepare next iteration */
1579 old_pstat = pstat;
1580 niter++;
1581 }
1582
1583 /* Warn the user if fixed point algorithm did not converge */
1584 if (niter > nitermax)
1585 bft_printf("Warning in thermodynamics computations for compressible"
1586 "flows:\n"
1587 "Fixed point algorithm did not converge when computing\n"
1588 "the subsonic inlet boundary condition with total\n"
1589 "pressure and total enthalpy imposed.\n"
1590 "At boundary face %ld,\n"
1591 "boundary Mach number residual = %12.4e,\n"
1592 "maximum number of iterations (%i) was reached.\n",
1593 (long)face_id, res, nitermax);
1594 }
1595 }
1596
1597 /*----------------------------------------------------------------------------*/
1598 /*!
1599 * \brief Compute epsilon sup:
1600 *
1601 * \f[\epsilon_{\textrm{sup}} = e - C_v T\f]
1602 *
1603 * for perfect gas: \f[\epsilon_{\textrm{sup}} = 0\f]
1604 *
1605 * \param[in] dens array of density values
1606 * \param[out] eps_sup epsilon sup array
1607 * \param[in] l_size l_size of the array
1608 */
1609 /*----------------------------------------------------------------------------*/
1610
1611 void
cs_cf_thermo_eps_sup(const cs_real_t * dens,cs_real_t * eps_sup,cs_lnum_t l_size)1612 cs_cf_thermo_eps_sup(const cs_real_t *dens,
1613 cs_real_t *eps_sup,
1614 cs_lnum_t l_size)
1615 {
1616 int ieos = cs_glob_cf_model->ieos;
1617
1618 /* single ideal gas or stiffened gas eos or ideal gas mixture
1619 (if ideal gas, infinite pressure equals 0) */
1620 if ( ieos == CS_EOS_IDEAL_GAS
1621 || ieos == CS_EOS_STIFFENED_GAS
1622 || ieos == CS_EOS_GAS_MIX) {
1623 cs_real_t psginf = cs_glob_cf_model->psginf;
1624
1625 for (cs_lnum_t ii = 0; ii < l_size; ii++)
1626 eps_sup[ii] = psginf / dens[ii];
1627 }
1628 /* TODO diffusion to be investigated for 2-phase homogeneous model */
1629 else if (ieos == CS_EOS_HOMOGENEOUS_TWO_PHASE) {
1630 for (cs_lnum_t ii = 0; ii < l_size; ii++)
1631 eps_sup[ii] = 0.;
1632 }
1633 else {
1634 for (cs_lnum_t ii = 0; ii < l_size; ii++)
1635 eps_sup[ii] = 0.;
1636 }
1637 }
1638
1639 /*----------------------------------------------------------------------------*/
1640 /*!
1641 * \brief This function is a driver allowing to call the appropriate
1642 * thermodynamical functions depending on the quantities provided by the user.
1643 * Hence it is only used during the initialization step and at the boundaries
1644 * of type supersonic inlet. It is described in the following how to select the
1645 * quantity to be returned.
1646 *
1647 * When calling the user subroutine, the integer 'iccfth' specifies which
1648 * calculation has to be performed (and which quantity has to be returned).
1649 * The values for 'iccfth' for each case are provided below.
1650 *
1651 * The variables are referred to using a different index i:
1652 *
1653 * - pressure: 2
1654 * - density: 3
1655 * - temperature: 5
1656 * - internal energy: 7
1657 * - entropy: 13
1658 *
1659 * iccfth is as follows, depending on which quantity needs to be computed:
1660 * - variables at cell centers from variable i and variable j (i<j):
1661 * iccfth = i*j*10000
1662 * - variables at boundary faces from variable i and variable j (i<j):
1663 * iccfth = i*j*10000+900
1664 *
1665 * Detailed values of iccfth and corresponding computations:
1666 *
1667 * Values at the cell centers:
1668 *
1669 * - temperature and energy from pressure and density: iccfth = 60000
1670 * - density and energy from pressure and temperature: iccfth = 100000
1671 * - density and temperature from pressure and energy: iccfth = 140000
1672 * - pressure and energy from density and temperature: iccfth = 150000
1673 * - pressure and temperature from density and energy: iccfth = 210000
1674 *
1675 * Values at the faces for boundary conditions:
1676 * - temperature and energy from pressure and density: iccfth = 60900
1677 * - density and energy from pressure and temperature: iccfth = 100900
1678 * - density and temperature from pressure and energy: iccfth = 140900
1679 * - pressure and energy from density and temperature: iccfth = 150900
1680 * - pressure and temperature from density and energy: iccfth = 210900
1681 *
1682 * \param[in] iccfth id of computation
1683 * \param[in] face_id face index if the computation is for a B.C.
1684 * -1 if the computation is located on cells
1685 * \param[in,out] bc_en total energy values at boundary faces
1686 * \param[in,out] bc_pr pressure values at boundary faces
1687 * \param[in,out] bc_tk temperature values at boundary faces
1688 * \param[in,out] bc_vel velocity values at boundary faces
1689 */
1690 /*----------------------------------------------------------------------------*/
1691
1692 void
cs_cf_thermo(const int iccfth,cs_lnum_t face_id,cs_real_t * bc_en,cs_real_t * bc_pr,cs_real_t * bc_tk,cs_real_3_t * bc_vel)1693 cs_cf_thermo(const int iccfth,
1694 cs_lnum_t face_id,
1695 cs_real_t *bc_en,
1696 cs_real_t *bc_pr,
1697 cs_real_t *bc_tk,
1698 cs_real_3_t *bc_vel)
1699 {
1700 const cs_mesh_t *m = cs_glob_mesh;
1701 const cs_lnum_t n_cells = m->n_cells;
1702 const cs_lnum_t *restrict b_face_cells
1703 = (const cs_lnum_t *restrict)m->b_face_cells;
1704
1705 /* Local variables */
1706 cs_lnum_t cell_id = -1;
1707
1708 if (face_id >= 0)
1709 cell_id = b_face_cells[face_id];
1710
1711 /* Map field arrays */
1712 cs_real_3_t *vel = (cs_real_3_t *)CS_F_(vel)->val;
1713 cs_real_t *cvar_pr = (cs_real_t *)CS_F_(p)->val;
1714 cs_real_t *crom = (cs_real_t *)CS_F_(rho)->val;
1715 cs_real_t *brom = (cs_real_t *)CS_F_(rho_b)->val;
1716 cs_real_t *cvar_tk = (cs_real_t *)CS_F_(t)->val;
1717 cs_real_t *cvar_en = (cs_real_t *)CS_F_(e_tot)->val;
1718
1719 /* Map specific heats field arrays - handle uniform cases */
1720 cs_real_t *cpro_cp = NULL;
1721 cs_real_t *cpro_cv = NULL;
1722 cs_real_t cpb = 0.;
1723 cs_real_t cvb = 0.;
1724
1725 if (CS_F_(cp) != NULL) {
1726 cpro_cp = (cs_real_t *)CS_F_(cp)->val;
1727 if (face_id >= 0) cpb = cpro_cp[cell_id];
1728 }
1729 if (CS_F_(cv) != NULL) {
1730 cpro_cv = (cs_real_t *)CS_F_(cv)->val;
1731 if (face_id >= 0) cvb = cpro_cv[cell_id];
1732 }
1733
1734 cs_real_t *cvar_fracv, *cvar_fracm, *cvar_frace;
1735 cvar_fracv = NULL;
1736 cvar_fracm = NULL;
1737 cvar_frace = NULL;
1738
1739 if (CS_F_(volume_f) != NULL){
1740 cvar_fracv = (cs_real_t *)CS_F_(volume_f)->val;
1741 cvar_fracm = (cs_real_t *)CS_F_(mass_f)->val;
1742 cvar_frace = (cs_real_t *)CS_F_(energy_f)->val;
1743 }
1744
1745 /* 0. Initialization. */
1746
1747 /* Calculation of temperature and energy from pressure and density */
1748 if (iccfth == 60000) {
1749 cs_cf_check_density(crom, n_cells);
1750 cs_cf_thermo_te_from_dp(cpro_cp, cpro_cv, cvar_pr, crom, cvar_tk, cvar_en,
1751 vel, n_cells);
1752 }
1753 /* Calculation of density and energy from pressure and temperature: */
1754 else if (iccfth == 100000) {
1755 cs_cf_check_temperature(cvar_tk, n_cells);
1756 cs_cf_thermo_de_from_pt(cpro_cp, cpro_cv, cvar_pr, cvar_tk, crom, cvar_en,
1757 vel, n_cells);
1758 }
1759 /* Calculation of density and temperature from pressure and energy */
1760 else if (iccfth == 140000) {
1761 cs_cf_thermo_dt_from_pe(cpro_cp, cpro_cv, cvar_pr, cvar_en, crom, cvar_tk,
1762 vel, n_cells);
1763 }
1764 /* Calculation of pressure and energy from density and temperature */
1765 else if (iccfth == 150000) {
1766 cs_cf_thermo_pe_from_dt(cpro_cp, cpro_cv, crom, cvar_tk, cvar_pr, cvar_en,
1767 vel, n_cells);
1768 }
1769 /* Calculation of pressure and temperature from density and energy */
1770 else if (iccfth == 210000) {
1771 cs_cf_thermo_pt_from_de(cpro_cp, cpro_cv, crom,
1772 cvar_en, cvar_pr, cvar_tk,
1773 vel,
1774 cvar_fracv,
1775 cvar_fracm,
1776 cvar_frace,
1777 n_cells);
1778 }
1779 /* Calculation of temperature and energy from pressure and density
1780 (it is postulated that the pressure and density values are strictly
1781 positive) */
1782 else if (iccfth == 60900) {
1783 cs_cf_thermo_te_from_dp(&cpb, &cvb, &bc_pr[face_id], &brom[face_id],
1784 &bc_tk[face_id], &bc_en[face_id], &bc_vel[face_id],
1785 1);
1786 }
1787 /* Calculation of density and energy from pressure and temperature */
1788 else if (iccfth == 100900) {
1789 cs_cf_thermo_de_from_pt(&cpb, &cvb, &bc_pr[face_id], &bc_tk[face_id],
1790 &brom[face_id], &bc_en[face_id], &bc_vel[face_id],
1791 1);
1792 }
1793 /* Calculation of density and temperature from pressure and total energy */
1794 else if (iccfth == 140900) {
1795 cs_cf_thermo_dt_from_pe(&cpb, &cvb, &bc_pr[face_id], &bc_en[face_id],
1796 &brom[face_id], &bc_tk[face_id], &bc_vel[face_id],
1797 1);
1798 }
1799 /* Calculation of pressure and energy from density and temperature */
1800 else if (iccfth == 150900) {
1801 cs_cf_thermo_pe_from_dt(&cpb, &cvb, &brom[face_id], &bc_tk[face_id],
1802 &bc_pr[face_id], &bc_en[face_id], &bc_vel[face_id],
1803 1);
1804 }
1805 /* Calculation of pressure and temperature from density and energy */
1806 else if (iccfth == 210900) {
1807 cs_real_t _fracv = cvar_fracv[cell_id];
1808 cs_real_t _fracm = cvar_fracm[cell_id];
1809 cs_real_t _frace = cvar_frace[cell_id];
1810 cs_cf_thermo_pt_from_de(&cpb, &cvb, &brom[face_id],
1811 &bc_en[face_id], &bc_pr[face_id],
1812 &bc_tk[face_id], &bc_vel[face_id],
1813 &_fracv,
1814 &_fracm,
1815 &_frace,
1816 1);
1817 }
1818 }
1819
1820 /*----------------------------------------------------------------------------*/
1821 /*!
1822 * \brief Compute density at boundary based on pressure and temperature.
1823 *
1824 * This is needed when imposing a mass flow rate on a given inlet.
1825 *
1826 * \param[in] face_id face id
1827 * \param[in] bc_pr pressure value at boundary face
1828 * \param[in] bc_tk temperature value at boundary face
1829 *
1830 * \return density at boundary face
1831 */
1832 /*----------------------------------------------------------------------------*/
1833
1834 cs_real_t
cs_cf_thermo_b_rho_from_pt(cs_lnum_t face_id,cs_real_t bc_pr,cs_real_t bc_tk)1835 cs_cf_thermo_b_rho_from_pt(cs_lnum_t face_id,
1836 cs_real_t bc_pr,
1837 cs_real_t bc_tk)
1838 {
1839 /* Local variables */
1840 int ieos = cs_glob_cf_model->ieos;
1841 cs_real_t psginf = cs_glob_cf_model->psginf;
1842
1843 cs_real_t b_rho = 0;
1844
1845 if (ieos == CS_EOS_IDEAL_GAS || ieos == CS_EOS_STIFFENED_GAS) {
1846 cs_real_t gamma0;
1847 cs_real_t cp0 = cs_glob_fluid_properties->cp0;
1848 cs_real_t cv0 = cs_glob_fluid_properties->cv0;
1849 cs_cf_thermo_gamma(&cp0, &cv0, &gamma0, 1);
1850
1851 b_rho = (bc_pr+psginf) / ((gamma0-1.)*bc_tk*cv0);
1852 }
1853 else if (ieos == CS_EOS_GAS_MIX) {
1854 cs_lnum_t cell_id = cs_glob_mesh->b_face_cells[face_id];
1855
1856 cs_real_t cp = CS_F_(cp)->val[cell_id];
1857 cs_real_t cv = CS_F_(cv)->val[cell_id];
1858
1859 cs_real_t gamma;
1860 cs_cf_thermo_gamma(&cp, &cv, &gamma, 1);
1861
1862 b_rho = (bc_pr+psginf) / ((gamma-1.)*bc_tk*cv);
1863 }
1864
1865 return b_rho;
1866 }
1867
1868 /*----------------------------------------------------------------------------*/
1869
1870 END_C_DECLS
1871