1 /*============================================================================
2 * Compute rotation/curvature correction.
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 <stdio.h>
35 #include <stdlib.h>
36 #include <string.h>
37 #include <math.h>
38
39 /*----------------------------------------------------------------------------
40 * Local headers
41 *----------------------------------------------------------------------------*/
42
43 #include "bft_mem.h"
44
45 #include "cs_base.h"
46 #include "cs_field.h"
47 #include "cs_field_pointer.h"
48 #include "cs_field_operator.h"
49 #include "cs_gradient.h"
50 #include "cs_math.h"
51 #include "cs_physical_constants.h"
52 #include "cs_rotation.h"
53 #include "cs_time_step.h"
54 #include "cs_turbulence_model.h"
55
56 /*----------------------------------------------------------------------------
57 * Header for the current file
58 *----------------------------------------------------------------------------*/
59
60 #include "cs_turbulence_rotation.h"
61
62 /*----------------------------------------------------------------------------*/
63
64 BEGIN_C_DECLS
65
66 /*=============================================================================
67 * Additional doxygen documentation
68 *============================================================================*/
69
70 /*!
71 * \file cs_turbulence_rotation.c
72 *
73 * Compute rotation/curvature correction for eddy-viscosity models.
74 *
75 * Two types of rotation/curvature correction can be computed, depending on
76 * the specific eddy-viscosity model:
77 *
78 * - itycor = 1: - Cazalbou correction (variable Ce2 coefficient in the
79 * destruction term of dissipation equation)
80 * - default correction for \f$ k - \epsilon \f$ type models,
81 * including elliptic relaxation/blending models
82 * (iturb = 20, 21, 50 or 51)
83 *
84 * - itycor = 2: - Spalart-Shur correction (production terms are multiplied
85 * by a rotation function)
86 * - default correction for \f$ k - \omega \f$ SST or
87 * Spalart-Allmaras (iturb = 60 or 70)
88 */
89
90 /*----------------------------------------------------------------------------*/
91
92 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
93
94 /*============================================================================
95 * Private function definitions
96 *============================================================================*/
97
98 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
99
100 /*=============================================================================
101 * Public function definitions
102 *============================================================================*/
103
104 /*----------------------------------------------------------------------------*/
105 /*!
106 * \brief Compute rotation/curvature correction for eddy-viscosity models.
107 *
108 * This function is called for the linear eddy viscosity RANS models,
109 * when irccor = 1 is verified.
110 *
111 * \param[in] dt time step (per cell)
112 * \param[out] rotfct rotation function of Spalart-Shur correction
113 * at cell center
114 * \param[out] ce2rc modified ce2 coeficient of Cazalbou correction
115 * at cell center
116 */
117 /*----------------------------------------------------------------------------*/
118
119 void
cs_turbulence_rotation_correction(const cs_real_t dt[],cs_real_t rotfct[],cs_real_t ce2rc[])120 cs_turbulence_rotation_correction(const cs_real_t dt[],
121 cs_real_t rotfct[],
122 cs_real_t ce2rc[])
123 {
124 const cs_mesh_t *m = cs_glob_mesh;
125 const cs_lnum_t n_cells = m->n_cells;
126 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
127
128 cs_real_t d1s2 =1./2.;
129
130 /* Initialization
131 ============== */
132
133 /* Empty field at the first iteration */
134
135 cs_real_6_t *cpro_straio
136 = (cs_real_6_t *)cs_field_by_name("strain_rate_tensor")->val;
137
138 /* Map field arrays */
139
140 const cs_real_3_t *vela = (const cs_real_3_t *)CS_F_(vel)->val_pre;
141 const cs_real_t *cvara_k = NULL;
142 const cs_real_t *cvara_ep = NULL;
143 const cs_real_t *cvara_omg = NULL;
144
145 if (cs_glob_turb_rans_model->itycor == 1){
146 cvara_k = (const cs_real_t *)CS_F_(k)->val_pre;
147 cvara_ep = (const cs_real_t *)CS_F_(eps)->val_pre;
148 }
149 else if (cs_glob_turb_rans_model->itycor == 2){
150 if (cs_glob_turb_model->iturb == 60){
151 cvara_omg = (const cs_real_t *)CS_F_(omg)->val_pre;
152 }
153 }
154
155 cs_real_t matrot[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
156
157 if (cs_glob_physical_constants->icorio == 1) {
158 /* In case of a rotating frame, all cells belong to the same "rotor" */
159 int r_num = 1;
160 const cs_rotation_t *r = cs_glob_rotation + r_num;
161 cs_rotation_add_coriolis_t(r, 1., matrot);
162 }
163
164 /* Preliminary calculations
165 ======================== */
166
167 /* Compute the strain rate and absolute vorticity tensor
168 ----------------------------------------------------- */
169
170 cs_real_6_t *strain = NULL;
171 cs_real_3_t *vortab = NULL;
172 cs_real_33_t *gradv = NULL;
173
174 BFT_MALLOC(strain, n_cells_ext, cs_real_6_t);
175 BFT_MALLOC(vortab, n_cells_ext, cs_real_3_t);
176 BFT_MALLOC(gradv, n_cells_ext, cs_real_33_t);
177
178 cs_field_gradient_vector(CS_F_(vel),
179 true, /* use_previous_t */
180 1 , /* inc */
181 gradv);
182
183 /* Compute the strain rate tensor (symmetric)
184 * S_ij = 0.5(dU_i/dx_j+dU_j/dx_i)
185 *
186 * and the absolute vorticity tensor (anti-symmetric)
187 * W_ij = 0.5(dU_i/dx_j-dU_j/dx_i) + e_imj*Omega_m
188 *
189 * Only the non zero components in the upper triangle are stored */
190
191 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
192 /* S11 */
193 strain[c_id][0] = gradv[c_id][0][0];
194 /* S22 */
195 strain[c_id][1] = gradv[c_id][1][1];
196 /* S33 */
197 strain[c_id][2] = gradv[c_id][2][2];
198 /* S12 */
199 strain[c_id][3] = d1s2*(gradv[c_id][0][1] + gradv[c_id][1][0]);
200 /* S13 */
201 strain[c_id][4] = d1s2*(gradv[c_id][0][2] + gradv[c_id][2][0]);
202 /* S23 */
203 strain[c_id][5] = d1s2*(gradv[c_id][1][2] + gradv[c_id][2][1]);
204 /* W12 */
205 vortab[c_id][0] = d1s2*(gradv[c_id][0][1] - gradv[c_id][1][0]) + matrot[1][0];
206 /* W13 */
207 vortab[c_id][1] = d1s2*(gradv[c_id][0][2] - gradv[c_id][2][0]) + matrot[2][0];
208 /* W23 */
209 vortab[c_id][2] = d1s2*(gradv[c_id][1][2] - gradv[c_id][2][1]) + matrot[2][1];
210 }
211
212 /* Partially free memory (strain and vortab arrays are deallocated later) */
213
214 BFT_FREE(gradv);
215
216 /* Computation of:
217 * --------------
218 *
219 * brtild = 2.W_ik.S_jk(DS_ij/Dt + (e_imn.S_jn + e_jmn.S_in)*Omega_m)
220 * eta1 = S_ij.S_ij
221 * eta2 = W_ij.W_ij
222 */
223
224 cs_real_63_t *grdsij = NULL;
225 cs_real_t *brtild = NULL;
226 cs_real_t *eta1 = NULL;
227 cs_real_t *eta2 = NULL;
228
229 BFT_MALLOC(grdsij, n_cells_ext, cs_real_63_t);
230 BFT_MALLOC(brtild, n_cells, cs_real_t);
231 BFT_MALLOC(eta1, n_cells, cs_real_t);
232 BFT_MALLOC(eta2, n_cells, cs_real_t);
233
234 /* Index connectivity */
235
236 /* istrai(i,j): position of the (i,j) component of the tensor */
237 /* in the strain and straio arrays */
238 /* ivorab(i,j): position of the (i,j) component of the tensor */
239 /* in the vortab array */
240 /* sigvor(i,j): sign of the (i,j) component of the absolute vorticity tensor */
241 /* = 1 if i > j */
242 /* = -1 if i < j */
243 /* = 0 if i = j */
244
245 int istrai[3][3];
246 int ivorab[3][3];
247 cs_real_t sigvor[3][3];
248
249 istrai[0][0] = 0;
250 istrai[1][1] = 1;
251 istrai[2][2] = 2;
252 istrai[1][0] = 3;
253 istrai[2][0] = 4;
254 istrai[2][1] = 5;
255 istrai[0][1] = istrai[1][0];
256 istrai[0][2] = istrai[2][0];
257 istrai[1][2] = istrai[2][1];
258
259 ivorab[0][0] = 0;
260 ivorab[1][1] = 0;
261 ivorab[2][2] = 0;
262 ivorab[1][0] = 0;
263 ivorab[2][0] = 1;
264 ivorab[2][1] = 2;
265 ivorab[0][1] = ivorab[1][0];
266 ivorab[0][2] = ivorab[2][0];
267 ivorab[1][2] = ivorab[2][1];
268
269 for (int i = 0; i < 3; i++) {
270 for (int j = 0; j < 3; j++) {
271 if (i < j)
272 sigvor[i][j] = -1.;
273 else if (i == j)
274 sigvor[i][j] = 0.;
275 else
276 sigvor[i][j] = 1.;
277 }
278 }
279
280 for (cs_lnum_t i = 0; i < n_cells; i++) {
281 brtild[i] = 0.;
282 eta1[i] = 0.;
283 eta2[i] = 0.;
284 }
285
286 const int key_cal_opt_id = cs_field_key_id("var_cal_opt");
287 cs_var_cal_opt_t var_cal_opt;
288
289 if ( cs_glob_turb_model->itytur == 2
290 || cs_glob_turb_model->itytur == 5
291 || cs_glob_turb_model->iturb == CS_TURB_K_OMEGA) {
292 cs_field_get_key_struct(CS_F_(k), key_cal_opt_id, &var_cal_opt);
293 }
294 else if (cs_glob_turb_model->iturb == CS_TURB_SPALART_ALLMARAS){
295 cs_field_get_key_struct(CS_F_(nusa), key_cal_opt_id, &var_cal_opt);
296 }
297
298 int nswrgp = var_cal_opt.nswrgr;
299 int imligp = var_cal_opt.imligr;
300 int iwarnp = var_cal_opt.verbosity;
301 cs_real_t epsrgp = var_cal_opt.epsrgr;
302 cs_real_t climgp = var_cal_opt.climgr;
303
304 int inc = 1;
305
306 cs_halo_type_t halo_type = CS_HALO_STANDARD;
307 cs_gradient_type_t gradient_type = CS_GRADIENT_GREEN_ITER;
308
309 cs_gradient_type_by_imrgra(var_cal_opt.imrgra,
310 &gradient_type,
311 &halo_type);
312
313 cs_gradient_tensor("Sij_gradient",
314 gradient_type,
315 halo_type,
316 inc,
317 nswrgp,
318 iwarnp,
319 imligp,
320 epsrgp,
321 climgp,
322 NULL, /* Use default Neumann BC */
323 NULL,
324 strain,
325 grdsij);
326
327 cs_real_t dsijdt, trrota, wiksjk;
328
329 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
330 for (cs_lnum_t i = 0; i < 3; i++) {
331 for (cs_lnum_t j = 0; j < 3; j++) {
332
333 /* material derivative of S_ij */
334 if (cs_glob_time_step_options->idtvar < 0)
335 dsijdt = 0.;
336 else
337 dsijdt = (strain[c_id][istrai[i][j]]
338 - cpro_straio[c_id][istrai[i][j]])/dt[c_id];
339
340 dsijdt += cs_math_3_dot_product(vela[c_id],
341 grdsij[c_id][istrai[i][j]]);
342
343 /* (e_imn.S_jn+e_jmn.S_in)*Omega_m term */
344
345 trrota = 0.;
346 for (cs_lnum_t k = 0; k < 3; k++) {
347 trrota += matrot[k][i] * strain[c_id][istrai[j][k]]
348 + matrot[k][j] * strain[c_id][istrai[i][k]];
349 }
350
351 /* W_ik.S_jk term */
352
353 wiksjk = 0.;
354 for (cs_lnum_t k = 0; k < 3; k++) {
355 wiksjk += sigvor[i][k]*vortab[c_id][ivorab[i][k]]
356 * strain[c_id][istrai[j][k]];
357 }
358
359 /* brtild, eta1, eta2 (see the definitions above) */
360
361 brtild[c_id] += 2. * wiksjk * (dsijdt + trrota);
362 eta1[c_id] += cs_math_pow2(strain[c_id][istrai[i][j]]);
363 eta2[c_id] += cs_math_pow2(sigvor[i][j] * vortab[c_id][ivorab[i][j]]);
364 }
365 }
366 }
367
368 /* Effective computation of the rotation correction
369 ================================================ */
370
371 cs_real_t stilde;
372 cs_real_t wtilde;
373
374 if (cs_glob_turb_rans_model->itycor == 1) {
375
376 /* 2.1 Cazalbou correction
377 ----------------------- */
378
379 cs_real_t xk;
380 cs_real_t xe;
381 cs_real_t rotild;
382 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
383
384 /* Computation of stilde = sqrt(2.S_ij.S_ij)
385 * and wtilde = sqrt(W_ij.W_ij/2) */
386
387 stilde = CS_MAX(sqrt(eta1[c_id] * 2.), 1.e-15);
388 wtilde = CS_MAX(sqrt(eta2[c_id] * 0.5), 1.e-15);
389
390 xk = CS_MAX(cvara_k[c_id], 1.e-15);
391 xe = CS_MAX(cvara_ep[c_id], 1.e-15);
392 rotild = xe/wtilde/xk;
393 brtild[c_id] = -brtild[c_id]*xk/xe/pow(stilde,3.);
394
395 /* Variable C_eps_2 coefficient of Cazalbou */
396
397 ce2rc[c_id] = cs_turb_ccaze2 + (cs_turb_ccaze2 - 1.)
398 / (1. + cs_turb_ccaza*pow(rotild,1.5))
399 + cs_turb_ccaze2*cs_turb_ccazsc*stilde*xk/xe
400 * ( tanh(cs_turb_ccazb*brtild[c_id] + cs_turb_ccazc)
401 - cs_turb_ccazd);
402
403 ce2rc[c_id] = CS_MAX(ce2rc[c_id], 0.);
404 }
405
406 }
407 else if (cs_glob_turb_rans_model->itycor == 2) {
408
409 /* Spalart-Shur correction
410 *------------------------
411 * (including modifications of Smirnov & Menter, ASME, 2009) */
412
413 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
414
415 /* Computation of stilde = 2.S_ij.S_ij and wtilde = 2.W_ij.W_ij */
416
417 stilde = CS_MAX(sqrt(eta1[c_id] * 2.), 1.e-15);
418 wtilde = CS_MAX(sqrt(eta2[c_id] * 2.), 1.e-15);
419
420 cs_real_t echtm2 = stilde;
421
422 if (cs_glob_turb_model->iturb == CS_TURB_K_OMEGA)
423 echtm2 = CS_MAX(echtm2, cs_turb_cmu*cs_math_pow2(cvara_omg[c_id]));
424
425 brtild[c_id] = brtild[c_id] / sqrt(wtilde*cs_math_pow3(echtm2));
426
427 cs_real_t rstar = sqrt(stilde)/sqrt(wtilde);
428
429 /* Rotation function of Spalart & Shur */
430
431 rotfct[c_id] = (1. + cs_turb_cssr1)*2.*rstar/(1. + rstar)
432 * (1. - cs_turb_cssr3*atan(cs_turb_cssr2*brtild[c_id]))
433 - cs_turb_cssr1;
434
435 rotfct[c_id] = CS_MIN(CS_MAX((double)rotfct[c_id],0.),1.25);
436
437 }
438 }
439
440 /* Finalization
441 ============ */
442
443 /* Save the strain rate tensor for the next time step */
444
445 if (cs_glob_time_step_options->idtvar >= 0) {
446 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
447 for (cs_lnum_t i = 0; i < 6; i++)
448 cpro_straio[c_id][i] = strain[c_id][i];
449 }
450 }
451
452 /* Free memory */
453
454 BFT_FREE(strain);
455 BFT_FREE(vortab);
456 BFT_FREE(grdsij);
457 BFT_FREE(brtild);
458 BFT_FREE(eta1);
459 BFT_FREE(eta2);
460 }
461
462 /*----------------------------------------------------------------------------*/
463
464 END_C_DECLS
465