1 /*============================================================================
2 * Functions to handle the resolution of the turbulence modelling
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 <stdlib.h>
35 #include <string.h>
36
37 #if defined(HAVE_MPI)
38 #include <mpi.h>
39 #endif
40
41 /*----------------------------------------------------------------------------
42 * Local headers
43 *----------------------------------------------------------------------------*/
44
45 #include <bft_mem.h>
46
47 #include "cs_post.h"
48 #include "cs_turbulence_model.h"
49 #include "cs_reco.h"
50
51 /*----------------------------------------------------------------------------
52 * Header for the current file
53 *----------------------------------------------------------------------------*/
54
55 #include "cs_cdo_turbulence.h"
56
57 /*----------------------------------------------------------------------------*/
58
59 BEGIN_C_DECLS
60
61 /*!
62 * \file cs_cdo_turbulence.c
63 *
64 * \brief Functions to handle the resolution of the turbulence modelling
65 * within the CDO framework
66 */
67
68 /*=============================================================================
69 * Local Macro definitions
70 *============================================================================*/
71
72 #define CS_CDO_TURBULENCE_DBG 0
73
74 /*============================================================================
75 * Type definitions
76 *============================================================================*/
77
78 /* --------------------------------------------------------------------------
79 * Context stucture for the k-epsilon turbulence modelling
80 * -------------------------------------------------------------------------- */
81
82 typedef struct {
83
84 /* Equations */
85
86 cs_equation_t *tke;
87 cs_equation_t *eps;
88
89 /* Properties associated to the two equations */
90
91 cs_property_t *tke_diffusivity; /* mu_t/sigma_k */
92 cs_property_t *eps_diffusivity; /* mu_t/sigma_e */
93 cs_property_t *sigma_k; /* TKE Schmidt */
94 cs_property_t *sigma_eps; /* epsilon Schmidt */
95
96 cs_xdef_t *tke_reaction; /* eps/tke by default + ... if needed */
97 cs_xdef_t *eps_reaction; /* by default + ... if needed */
98
99 cs_xdef_t *tke_source_term; /* Production + buoyancy if needed for k */
100 cs_xdef_t *eps_source_term; /* Same for epsilon */
101
102 } cs_turb_context_k_eps_t;
103
104 /* --------------------------------------------------------------------------
105 * Context structure for the k-epsilon turbulence modelling
106 * -------------------------------------------------------------------------- */
107
108 typedef struct {
109
110 /* High level structures */
111 const cs_mesh_t *mesh;
112 const cs_cdo_connect_t *connect;
113 const cs_cdo_quantities_t *quant;
114 const cs_time_step_t *time_step;
115
116 /* Turbulence structure */
117 cs_turbulence_t *tbs;
118
119 /* Velocities arrays */
120 cs_real_t *u_cell;
121 cs_real_t *u_face;
122
123 } cs_turb_source_term_t;
124
125
126 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
127
128 /*============================================================================
129 * Private variables
130 *============================================================================*/
131
132 /*============================================================================
133 * Private function prototypes
134 *============================================================================*/
135
136 /*----------------------------------------------------------------------------*/
137 /*!
138 * \brief Generic function pointer for defining a quantity at known locations
139 * Here at cells with a function. elt_ids is optional. If not NULL,
140 * the function works on a sub-list of elements. Moreover, it enables
141 * to fill retval with an indirection if dense_output is set to false
142 * Case of k-epsilon model
143 *
144 */
145 /*----------------------------------------------------------------------------*/
146
147 static void
_prepare_ke(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * time_step,const cs_turbulence_t * tbs,cs_real_t * tke_reaction,cs_real_t * eps_reaction,cs_real_t * tke_source_term,cs_real_t * eps_source_term)148 _prepare_ke(const cs_mesh_t *mesh,
149 const cs_cdo_connect_t *connect,
150 const cs_cdo_quantities_t *quant,
151 const cs_time_step_t *time_step,
152 const cs_turbulence_t *tbs,
153 cs_real_t *tke_reaction,
154 cs_real_t *eps_reaction,
155 cs_real_t *tke_source_term,
156 cs_real_t *eps_source_term)
157 {
158 const cs_turbulence_param_t *tbp = tbs->param;
159 cs_turb_context_k_eps_t *kec = (cs_turb_context_k_eps_t *)tbs->context;
160
161 const cs_real_t *u_cell = cs_equation_get_cell_values(tbs->mom_eq, false);
162 const cs_real_t *u_face = cs_equation_get_face_values(tbs->mom_eq, false);
163
164 cs_real_t *mu_t = tbs->mu_t_field->val;
165
166 /* Get the evaluation of rho */
167 int rho_stride = 0;
168 cs_real_t *rho = NULL;
169
170 const cs_real_t *tke_cell = cs_equation_get_cell_values(kec->tke, false);
171 const cs_real_t *eps_cell = cs_equation_get_cell_values(kec->eps, false);
172
173 /* Get mass density values in each cell */
174 cs_property_iso_get_cell_values(time_step->t_cur, tbs->rho,
175 &rho_stride, &rho);
176
177 const cs_real_t d1s3 = 1./3.;
178 const cs_real_t d2s3 = 2./3.;
179
180 /* Production term */
181 if (tbp->model->iturb == CS_TURB_K_EPSILON) {
182 # pragma omp parallel for if (mesh->n_cells > CS_THR_MIN)
183 for (cs_lnum_t c_id = 0; c_id < mesh->n_cells; c_id++) {
184
185 cs_real_t grd_uc[3][3];
186 /* Compute the velocity gradient */
187 cs_reco_grad_33_cell_from_fb_dofs(c_id, connect, quant,
188 u_cell, u_face, (cs_real_t *)grd_uc);
189
190 cs_real_t strain_sq = 2.
191 * ( cs_math_pow2( d2s3*grd_uc[0][0]
192 - d1s3*grd_uc[1][1]
193 - d1s3*grd_uc[2][2])
194 + cs_math_pow2(- d1s3*grd_uc[0][0]
195 + d2s3*grd_uc[1][1]
196 - d1s3*grd_uc[2][2])
197 + cs_math_pow2(- d1s3*grd_uc[0][0]
198 - d1s3*grd_uc[1][1]
199 + d2s3*grd_uc[2][2]))
200 + cs_math_pow2(grd_uc[0][1] + grd_uc[1][0])
201 + cs_math_pow2(grd_uc[0][2] + grd_uc[2][0])
202 + cs_math_pow2(grd_uc[1][2] + grd_uc[2][1]);
203
204 tke_source_term[c_id] = mu_t[c_id] * strain_sq;
205 }
206 }
207 else if (tbp->model->iturb == CS_TURB_K_EPSILON_LIN_PROD) {
208 # pragma omp parallel for if (mesh->n_cells > CS_THR_MIN)
209 for (cs_lnum_t c_id = 0; c_id < mesh->n_cells; c_id++) {
210
211 cs_real_t grd_uc[3][3];
212 /* Compute the velocity gradient */
213 cs_reco_grad_33_cell_from_fb_dofs(c_id, connect, quant,
214 u_cell, u_face, (cs_real_t *)grd_uc);
215
216 cs_real_t strain_sq = 2.
217 * ( cs_math_pow2( d2s3*grd_uc[0][0]
218 - d1s3*grd_uc[1][1]
219 - d1s3*grd_uc[2][2])
220 + cs_math_pow2(- d1s3*grd_uc[0][0]
221 + d2s3*grd_uc[1][1]
222 - d1s3*grd_uc[2][2])
223 + cs_math_pow2(- d1s3*grd_uc[0][0]
224 - d1s3*grd_uc[1][1]
225 + d2s3*grd_uc[2][2]))
226 + cs_math_pow2(grd_uc[0][1] + grd_uc[1][0])
227 + cs_math_pow2(grd_uc[0][2] + grd_uc[2][0])
228 + cs_math_pow2(grd_uc[1][2] + grd_uc[2][1]);
229
230 cs_real_t strain = sqrt(strain_sq);
231 const cs_real_t sqrcmu = sqrt(cs_turb_cmu);
232 cs_real_t cmueta =
233 fmin(cs_turb_cmu*tke_cell[c_id]/eps_cell[c_id] * strain, sqrcmu);
234
235 tke_source_term[c_id] = rho[c_id*rho_stride] * cmueta * strain * tke_cell[c_id];
236
237 }
238 }
239
240 /* Implicit dissipation term and epsilon source term */
241 # pragma omp parallel for if (mesh->n_cells > CS_THR_MIN)
242 for (cs_lnum_t c_id = 0; c_id < mesh->n_cells; c_id++) {
243 /* Inverse integral time scale */
244 cs_real_t d_ttke = eps_cell[c_id] / tke_cell[c_id];
245
246 /* Ce1 * eps/k * P */
247 eps_source_term[c_id] = cs_turb_ce1 * d_ttke * tke_source_term[c_id];
248
249 tke_reaction[c_id] = rho[c_id*rho_stride] * d_ttke;
250
251 /* TODO Get Ce2 from curvature correction, to be merged with legacy */
252 eps_reaction[c_id] = cs_turb_ce2 * rho[c_id*rho_stride] * d_ttke;
253
254 }
255
256 }
257
258 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
259
260 /*============================================================================
261 * Public function prototypes
262 *============================================================================*/
263
264 /*----------------------------------------------------------------------------*/
265 /*!
266 * \brief Allocate the structure storing the set of parameters for the
267 * turbulence modelling
268 *
269 * \return a pointer to a new allocated cs_turbulence_param_t structure
270 */
271 /*----------------------------------------------------------------------------*/
272
273 cs_turbulence_param_t *
cs_turbulence_param_create(void)274 cs_turbulence_param_create(void)
275 {
276 cs_turbulence_param_t *tbp = NULL;
277
278 BFT_MALLOC(tbp, 1, cs_turbulence_param_t);
279
280 /* The following structures are shared with the Legacy part */
281
282 tbp->model = cs_get_glob_turb_model();
283 tbp->rans_param = cs_get_glob_turb_rans_model();
284 tbp->les_param = cs_get_glob_turb_les_model();
285 tbp->reference_values = cs_get_glob_turb_ref_values();
286
287 return tbp;
288 }
289
290 /*----------------------------------------------------------------------------*/
291 /*!
292 * \brief Allocate the structure managing the turbulence modelling
293 *
294 * \param[in] tbp pointer to a \ref cs_turbulence_param_t structure
295 *
296 * \return a pointer to a new allocated cs_turbulence_t structure
297 */
298 /*----------------------------------------------------------------------------*/
299
300 cs_turbulence_t *
cs_turbulence_create(cs_turbulence_param_t * tbp)301 cs_turbulence_create(cs_turbulence_param_t *tbp)
302 {
303 cs_turbulence_t *tbs = NULL;
304
305 BFT_MALLOC(tbs, 1, cs_turbulence_t);
306
307 /* All the members of the following structures are shared with the Legacy
308 * part. This structure is owned by cs_navsto_param_t
309 */
310 tbs->param = tbp;
311 tbs->mom_eq = NULL;
312
313 /* Properties */
314 tbs->rho = NULL; /* Mass density */
315 tbs->mu_tot = NULL; /* Total viscosity */
316 tbs->mu_l = NULL; /* Laminar viscosity */
317 tbs->mu_t = NULL; /* Turbulent viscosity */
318
319 tbs->mu_tot_array = NULL;
320
321 /* Fields */
322 tbs->mu_t_field = NULL; /* Turbulent viscosity */
323 tbs->rij = NULL; /* Reynolds stress tensor */
324
325 /* Main structure (cast on-the-fly according to the turbulence model) */
326 tbs->context = NULL;
327
328 /* Function pointers */
329 tbs->init_context = NULL;
330 tbs->free_context = NULL;
331 tbs->compute = NULL;
332 tbs->compute_steady = NULL;
333 tbs->update = NULL;
334
335 return tbs;
336 }
337
338 /*----------------------------------------------------------------------------*/
339 /*!
340 * \brief Free the structure managing the turbulence modelling
341 *
342 * \param[in, out] p_turb_struct pointer to the structure to free
343 */
344 /*----------------------------------------------------------------------------*/
345
346 void
cs_turbulence_free(cs_turbulence_t ** p_turb_struct)347 cs_turbulence_free(cs_turbulence_t **p_turb_struct)
348 {
349 cs_turbulence_t *tbs = *p_turb_struct;
350
351 /* Set of parameters (members are shared and freed elsewhere).
352 * Properties, equations and fields are freed in an other part of the code
353 */
354
355 BFT_FREE(tbs->mu_tot_array);
356
357 if (tbs->free_context != NULL)
358 tbs->context = tbs->free_context(tbs->context);
359
360 assert(tbs->context == NULL);
361 BFT_FREE(tbs);
362 *p_turb_struct = NULL;
363 }
364
365 /*----------------------------------------------------------------------------*/
366 /*!
367 * \brief Initialize the structure managing the turbulence modelling
368 *
369 * \param[in, out] tbs pointer to the structure to initialize
370 * \param[in] mom_eq pointer to the structure mom_eq
371 */
372 /*----------------------------------------------------------------------------*/
373
374 void
cs_turbulence_init_setup(cs_turbulence_t * tbs,cs_equation_t * mom_eq)375 cs_turbulence_init_setup(cs_turbulence_t *tbs,
376 cs_equation_t *mom_eq)
377 {
378 if (tbs == NULL)
379 return;
380
381 const cs_turbulence_param_t *tbp = tbs->param;
382 const cs_turb_model_t *model = tbp->model;
383
384 if (model->type == CS_TURB_NONE)
385 return; /* Nothing to do if there is a laminar flow */
386
387 tbs->mom_eq = mom_eq;
388
389 /* Set field metadata */
390 const int log_key = cs_field_key_id("log");
391 const int post_key = cs_field_key_id("post_vis");
392 const bool has_previous = false;
393 const int field_post_flag = CS_POST_ON_LOCATION | CS_POST_MONITOR;
394
395 int field_mask = CS_FIELD_INTENSIVE | CS_FIELD_PROPERTY | CS_FIELD_CDO;
396 int location_id = cs_mesh_location_get_id_by_name("cells");
397
398 tbs->mu_t_field = cs_field_find_or_create(CS_NAVSTO_TURB_VISCOSITY,
399 field_mask,
400 location_id,
401 1, /* dimension */
402 has_previous);
403
404 /* Set default value for keys related to log and post-processing */
405 cs_field_set_key_int(tbs->mu_t_field, log_key, 1);
406 cs_field_set_key_int(tbs->mu_t_field, post_key, field_post_flag);
407
408 /* Properties (shared) */
409 tbs->rho = cs_property_by_name(CS_PROPERTY_MASS_DENSITY);
410 tbs->mu_tot = cs_property_by_name(CS_NAVSTO_TOTAL_VISCOSITY);
411 tbs->mu_l = cs_property_by_name(CS_NAVSTO_LAM_VISCOSITY);
412
413 assert(tbs->rho != NULL && tbs->mu_l != NULL && tbs->mu_tot != NULL);
414
415 /* Add a mu_t property and define it with the associated field */
416 tbs->mu_t = cs_property_add(CS_NAVSTO_TURB_VISCOSITY, CS_PROPERTY_ISO);
417
418 cs_property_def_by_field(tbs->mu_t, tbs->mu_t_field);
419
420 /* Set function pointers and initialize the context structure */
421 switch (model->iturb) {
422
423 case CS_TURB_K_EPSILON:
424 case CS_TURB_K_EPSILON_LIN_PROD:
425 tbs->init_context = cs_turb_init_k_eps_context;
426 tbs->free_context = cs_turb_free_k_eps_context;
427 tbs->compute = cs_turb_compute_k_eps;
428 tbs->update = cs_turb_update_k_eps;
429
430 tbs->context = tbs->init_context(model);
431 break;
432
433 case CS_TURB_NONE:
434 break;
435
436 default:
437 bft_error(__FILE__, __LINE__, 0,
438 " %s: Invalid turbulence model with CDO schemes.\n"
439 " Please check your settings.", __func__);
440 break;
441 }
442
443 }
444
445 /*----------------------------------------------------------------------------*/
446 /*!
447 * \brief Finalize the setup of the turbulence modelling and especially the
448 * equations/properties and other related quantities
449 *
450 * \param[in] mesh pointer to a \ref cs_mesh_t structure
451 * \param[in] connect pointer to a cs_cdo_connect_t structure
452 * \param[in] quant pointer to a cs_cdo_quantities_t structure
453 * \param[in] time_step pointer to a cs_time_step_t structure
454 * \param[in, out] tbs pointer to the turbulence main structure
455 */
456 /*----------------------------------------------------------------------------*/
457
458 void
cs_turbulence_finalize_setup(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * time_step,cs_turbulence_t * tbs)459 cs_turbulence_finalize_setup(const cs_mesh_t *mesh,
460 const cs_cdo_connect_t *connect,
461 const cs_cdo_quantities_t *quant,
462 const cs_time_step_t *time_step,
463 cs_turbulence_t *tbs)
464 {
465 CS_UNUSED(mesh);
466 CS_UNUSED(connect);
467 CS_UNUSED(time_step);
468
469 if (tbs == NULL)
470 return;
471
472 const cs_turbulence_param_t *tbp = tbs->param;
473 const cs_turb_model_t *model = tbp->model;
474
475 if (model->type == CS_TURB_NONE)
476 return; /* Nothing to do */
477
478 /* Define the property related to the total viscosity */
479 BFT_MALLOC(tbs->mu_tot_array, quant->n_cells, cs_real_t);
480 memset(tbs->mu_tot_array, 0, quant->n_cells*sizeof(cs_real_t));
481
482 cs_property_def_by_array(tbs->mu_tot,
483 cs_flag_primal_cell,
484 tbs->mu_tot_array,
485 false, /* definition is owner ? */
486 NULL); /* no index */
487
488 /* Last setup for each turbulence model */
489 switch (model->iturb) {
490
491 case CS_TURB_K_EPSILON:
492 case CS_TURB_K_EPSILON_LIN_PROD:
493 {
494 /* Add a source term after having retrieved the equation param related to
495 the turbulent kinetic energy equation */
496 cs_turb_context_k_eps_t *kec = (cs_turb_context_k_eps_t *)tbs->context;
497 cs_equation_param_t *tke_eqp = cs_equation_get_param(kec->tke);
498 kec->tke_source_term =
499 cs_equation_add_source_term_by_array(tke_eqp,
500 NULL, /* all cells */
501 cs_flag_primal_cell,
502 NULL,
503 false, /*is owner */
504 NULL); /* index */
505
506 kec->tke_reaction =
507 cs_property_def_by_array(cs_property_by_name("k_reaction"),
508 cs_flag_primal_cell,
509 NULL,
510 false, /* definition is owner ? */
511 NULL); /* no index */
512
513 cs_equation_param_t *eps_eqp = cs_equation_get_param(kec->eps);
514 kec->eps_source_term =
515 cs_equation_add_source_term_by_array(eps_eqp,
516 NULL, /* all cells */
517 cs_flag_primal_cell,
518 NULL,
519 false, /*is owner */
520 NULL); /* index */
521
522 kec->eps_reaction =
523 cs_property_def_by_array(cs_property_by_name("eps_reaction"),
524 cs_flag_primal_cell,
525 NULL,
526 false, /* definition is owner ? */
527 NULL); /* no index */
528
529 cs_property_def_by_array(tbs->mu_tot,
530 cs_flag_primal_cell,
531 tbs->mu_tot_array,
532 false, /* definition is owner ? */
533 NULL); /* no index */
534
535 /* Initialize TKE */
536 cs_turb_ref_values_t *t_ref= cs_get_glob_turb_ref_values();
537 cs_real_t tke_ref = 1.5 * cs_math_pow2(0.02 * t_ref->uref);
538 cs_equation_add_ic_by_value(tke_eqp,
539 NULL,
540 &tke_ref);
541
542 /* Initialize epsilon */
543 cs_real_t eps_ref = powf(tke_ref, 1.5) * cs_turb_cmu / t_ref->almax;
544 cs_equation_add_ic_by_value(eps_eqp,
545 NULL,
546 &eps_ref);
547
548 }
549 break;
550
551 default:
552 bft_error(__FILE__, __LINE__, 0,
553 " %s: Invalid turbulence model with CDO schemes.\n"
554 " Please check your settings.", __func__);
555 break;
556 }
557
558 }
559
560 /*----------------------------------------------------------------------------*/
561 /*!
562 * \brief Initialize quantities related to a turbulence model.
563 *
564 * \param[in] mesh pointer to a \ref cs_mesh_t structure
565 * \param[in] connect pointer to a cs_cdo_connect_t structure
566 * \param[in] quant pointer to a cs_cdo_quantities_t structure
567 * \param[in] time_step pointer to a cs_time_step_t structure
568 * \param[in, out] tbs pointer to the turbulence main structure
569 */
570 /*----------------------------------------------------------------------------*/
571
572 void
cs_turbulence_initialize(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * time_step,cs_turbulence_t * tbs)573 cs_turbulence_initialize(const cs_mesh_t *mesh,
574 const cs_cdo_connect_t *connect,
575 const cs_cdo_quantities_t *quant,
576 const cs_time_step_t *time_step,
577 cs_turbulence_t *tbs)
578 {
579 CS_UNUSED(mesh);
580 CS_UNUSED(connect);
581
582 if (tbs == NULL)
583 return;
584
585 const cs_turbulence_param_t *tbp = tbs->param;
586 const cs_turb_model_t *model = tbp->model;
587
588 if (model->iturb == CS_TURB_NONE)
589 return; /* Nothing to do */
590
591 tbs->update(mesh, connect, quant, time_step, tbs);
592
593 }
594
595 /*----------------------------------------------------------------------------*/
596 /*!
597 * \brief Allocate and initialize the context structure related to the
598 * k-epsilon turbulence model
599 *
600 * \param[in] tbm structure which defines the turbulence model
601 *
602 * \return a pointer to a new allocated context structure
603 */
604 /*----------------------------------------------------------------------------*/
605
606 void *
cs_turb_init_k_eps_context(const cs_turb_model_t * tbm)607 cs_turb_init_k_eps_context(const cs_turb_model_t *tbm)
608 {
609 if (tbm == NULL)
610 return NULL;
611
612 /* Sanity checks */
613 assert((tbm->iturb == CS_TURB_K_EPSILON) ||
614 (tbm->iturb == CS_TURB_K_EPSILON_LIN_PROD));
615 assert(tbm->type == CS_TURB_RANS);
616 assert(tbm->order == CS_TURB_FIRST_ORDER);
617
618 cs_turb_context_k_eps_t *kec = NULL;
619
620 BFT_MALLOC(kec, 1, cs_turb_context_k_eps_t);
621
622 /* Add new equations for the turbulent kinetic energy (tke) and the
623 dissipation (epsilon) */
624
625 kec->tke = cs_equation_add("k", /* equation name */
626 "k", /* variable name */
627 CS_EQUATION_TYPE_NAVSTO, /* related to NS */
628 1,
629 CS_PARAM_BC_HMG_NEUMANN);
630
631 kec->eps = cs_equation_add("epsilon", /* equation name */
632 "epsilon", /* variable name */
633 CS_EQUATION_TYPE_NAVSTO,
634 1,
635 CS_PARAM_BC_HMG_NEUMANN);
636
637 /* Add new related properties which will be associated to discretization
638 terms in tke and epsilon */
639
640 kec->tke_diffusivity = cs_property_add("k_diffusivity", CS_PROPERTY_ISO);
641
642 kec->eps_diffusivity = cs_property_add("epsilon_diffusivity",
643 CS_PROPERTY_ISO);
644
645 /* Turbulent Schmidt coefficients : creation and set the reference value */
646
647 kec->sigma_k = cs_property_add("k_turb_schmidt", CS_PROPERTY_ISO);
648 cs_property_set_reference_value(kec->sigma_k, 1.0);
649
650 kec->sigma_eps = cs_property_add("epsilon_turb_schmidt", CS_PROPERTY_ISO);
651 cs_property_set_reference_value(kec->sigma_eps, 1.3);
652
653 /* Reaction (implicit source terms) coefficients */
654
655 cs_property_t *k_reaction
656 = cs_property_add("k_reaction", CS_PROPERTY_ISO);
657 cs_property_t *eps_reaction
658 = cs_property_add("epsilon_reaction", CS_PROPERTY_ISO);
659
660 /* Retrieve the mass density */
661
662 cs_property_t *mass_density = cs_property_by_name(CS_PROPERTY_MASS_DENSITY);
663
664 /* Retrieve the advection field from Navier--Stokes (the mass flux) */
665
666 cs_adv_field_t *adv = cs_advection_field_by_name("mass_flux");
667 if (adv == NULL)
668 bft_error(__FILE__, __LINE__, 0,
669 " %s: Advection field not found.", __func__);
670
671 /* Add terms to the TKE equation */
672
673 cs_equation_param_t *tke_eqp = cs_equation_get_param(kec->tke);
674
675 cs_equation_add_time(tke_eqp, mass_density);
676 cs_equation_add_diffusion(tke_eqp, kec->tke_diffusivity);
677 cs_equation_add_reaction(tke_eqp, k_reaction);
678 cs_equation_add_advection(tke_eqp, adv);
679
680 /* Source term is defined elsewhere since it depends on the choice of the
681 * sub-model */
682
683 /* Add terms to the epsilon equation */
684
685 cs_equation_param_t *eps_eqp = cs_equation_get_param(kec->eps);
686
687 cs_equation_add_time(eps_eqp, mass_density);
688 cs_equation_add_diffusion(eps_eqp, kec->eps_diffusivity);
689 cs_equation_add_reaction(eps_eqp, eps_reaction);
690 cs_equation_add_advection(eps_eqp, adv);
691
692 return kec;
693 }
694
695 /*----------------------------------------------------------------------------*/
696 /*!
697 * \brief Free the context structure related to the k-epsilon turbulence model
698 *
699 * \param[in, out] tbc pointer to a structure context cast on-the-fly
700 *
701 * \return a NULL pointer
702 */
703 /*----------------------------------------------------------------------------*/
704
705 void *
cs_turb_free_k_eps_context(void * tbc)706 cs_turb_free_k_eps_context(void *tbc)
707 {
708 cs_turb_context_k_eps_t *kec = (cs_turb_context_k_eps_t *)tbc;
709
710 if (kec == NULL)
711 return kec;
712
713 BFT_FREE(kec);
714
715 return kec;
716 }
717
718 /*----------------------------------------------------------------------------*/
719 /*!
720 * \brief Update for the current time step the new state for the turbulence
721 * model. This is used to update the turbulent viscosity.
722 *
723 * \param[in] mesh pointer to a \ref cs_mesh_t structure
724 * \param[in] connect pointer to a cs_cdo_connect_t structure
725 * \param[in] quant pointer to a cs_cdo_quantities_t structure
726 * \param[in] time_step structure managing the time stepping
727 * \param[in] tbs pointer to a \ref cs_turbulence_t structure
728 */
729 /*----------------------------------------------------------------------------*/
730
731 void
cs_turb_update_k_eps(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * time_step,const cs_turbulence_t * tbs)732 cs_turb_update_k_eps(const cs_mesh_t *mesh,
733 const cs_cdo_connect_t *connect,
734 const cs_cdo_quantities_t *quant,
735 const cs_time_step_t *time_step,
736 const cs_turbulence_t *tbs)
737 {
738 CS_UNUSED(connect);
739 CS_UNUSED(quant);
740
741 if (tbs == NULL)
742 return;
743
744 cs_lnum_t n_cells = mesh->n_cells;
745 CS_UNUSED(n_cells); /* avoid a compiler warning without OpenMP */
746
747 cs_turb_context_k_eps_t *kec =
748 (cs_turb_context_k_eps_t *)tbs->context;
749
750 /* Update turbulent viscosity field */
751 cs_real_t *mu_t = tbs->mu_t_field->val;
752 cs_real_t *k = cs_equation_get_cell_values(kec->tke, false);
753 cs_real_t *eps = cs_equation_get_cell_values(kec->eps, false);
754
755 /* Get the evaluation of rho */
756 int rho_stride = 0;
757 cs_real_t *rho = NULL;
758
759 /* Get mass density values in each cell */
760 cs_property_iso_get_cell_values(time_step->t_cur, tbs->rho,
761 &rho_stride, &rho);
762
763 /* Get laminar viscosity values in each cell */
764 int mu_stride = 0;
765 cs_real_t *mu_l = NULL;
766 cs_property_iso_get_cell_values(time_step->t_cur, tbs->mu_l,
767 &mu_stride, &mu_l);
768
769
770 /* Compute mu_t in each cell and define mu_tot = mu_t + mu_l */
771 # pragma omp parallel for if (n_cells > CS_THR_MIN)
772 for (cs_lnum_t cell_id = 0; cell_id < mesh->n_cells; cell_id++) {
773
774 mu_t[cell_id] = cs_turb_cmu * rho[cell_id*rho_stride] *
775 cs_math_pow2(k[cell_id]) / eps[cell_id];
776
777 tbs->mu_tot_array[cell_id] = mu_t[cell_id] + mu_l[cell_id*mu_stride];
778
779 }
780
781 /* Free memory */
782 BFT_FREE(rho);
783 BFT_FREE(mu_l);
784 }
785
786 /*----------------------------------------------------------------------------*/
787 /*!
788 * \brief Compute for the current time step the new state for the turbulence
789 * model. This means that all related equations are built and then
790 * solved.
791 *
792 * \param[in] mesh pointer to a \ref cs_mesh_t structure
793 * \param[in] connect pointer to a cs_cdo_connect_t structure
794 * \param[in] quant pointer to a cs_cdo_quantities_t structure
795 * \param[in] time_step structure managing the time stepping
796 * \param[in, out] tbs pointer to turbulence structure
797 */
798 /*----------------------------------------------------------------------------*/
799
800 void
cs_turb_compute_k_eps(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * time_step,cs_turbulence_t * tbs)801 cs_turb_compute_k_eps(const cs_mesh_t *mesh,
802 const cs_cdo_connect_t *connect,
803 const cs_cdo_quantities_t *quant,
804 const cs_time_step_t *time_step,
805 cs_turbulence_t *tbs)
806 {
807 if (tbs == NULL)
808 return;
809
810 /* Get k epsilon context */
811 cs_turb_context_k_eps_t *kec = (cs_turb_context_k_eps_t *)tbs->context;
812 cs_equation_t *tke_eq = kec->tke;
813 cs_equation_t *eps_eq = kec->eps;
814 assert(kec != NULL);
815
816 /* Prepare source term and reaction term */
817 cs_real_t *tke_source_term = NULL, *eps_source_term = NULL;
818 cs_real_t *tke_reaction = NULL, *eps_reaction = NULL;
819 BFT_MALLOC(tke_source_term, mesh->n_cells, cs_real_t);
820 BFT_MALLOC(eps_source_term, mesh->n_cells, cs_real_t);
821 BFT_MALLOC(tke_reaction, mesh->n_cells, cs_real_t);
822 BFT_MALLOC(eps_reaction, mesh->n_cells, cs_real_t);
823
824 /* Set xdefs */
825 cs_xdef_set_array(kec->tke_reaction,
826 false, /* is_owner */
827 tke_reaction);
828
829 cs_xdef_set_array(kec->eps_reaction,
830 false, /* is_owner */
831 eps_reaction);
832
833 cs_xdef_set_array(kec->tke_source_term,
834 false, /* is_owner */
835 tke_source_term);
836
837 cs_xdef_set_array(kec->eps_source_term,
838 false, /* is_owner */
839 eps_source_term);
840
841 _prepare_ke(mesh,
842 connect,
843 quant,
844 time_step,
845 tbs,
846 tke_reaction,
847 eps_reaction,
848 tke_source_term,
849 eps_source_term);
850
851 /* Solve k */
852 cs_equation_solve(true, /* cur2prev */
853 mesh,
854 tke_eq);
855
856 /* Solve epsilon */
857 cs_equation_solve(true, /* cur2prev */
858 mesh,
859 eps_eq);
860
861 /* Free memory */
862 BFT_FREE(tke_source_term);
863 BFT_FREE(eps_source_term);
864 BFT_FREE(tke_reaction);
865 BFT_FREE(eps_reaction);
866
867 }
868
869 /*----------------------------------------------------------------------------*/
870
871 END_C_DECLS
872