1
2 /*
3 This file is part of Code_Saturne, a general-purpose CFD tool.
4
5 Copyright (C) 1998-2021 EDF S.A.
6
7 This program is free software; you can redistribute it and/or modify it under
8 the terms of the GNU General Public License as published by the Free Software
9 Foundation; either version 2 of the License, or (at your option) any later
10 version.
11
12 This program is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
15 details.
16
17 You should have received a copy of the GNU General Public License along with
18 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
19 Street, Fifth Floor, Boston, MA 02110-1301, USA.
20 */
21
22 /*----------------------------------------------------------------------------*/
23
24 #include "cs_defs.h"
25
26 /*----------------------------------------------------------------------------
27 * Standard C library headers
28 *----------------------------------------------------------------------------*/
29
30 #include <assert.h>
31 #include <math.h>
32 #include <string.h>
33
34 #if defined(HAVE_MPI)
35 #include <mpi.h>
36 #endif
37
38 /*----------------------------------------------------------------------------
39 * PLE library headers
40 *----------------------------------------------------------------------------*/
41
42 #include <ple_coupling.h>
43
44 /*----------------------------------------------------------------------------
45 * Local headers
46 *----------------------------------------------------------------------------*/
47
48 #include "cs_headers.h"
49
50 /*----------------------------------------------------------------------------*/
51
52 BEGIN_C_DECLS
53
54 /*----------------------------------------------------------------------------*/
55 /*!
56 * \file cs_user_les_inflow-base.c
57 *
58 * \brief Generation of synthetic turbulence at LES inlets initialization.
59 *
60 * See \ref les_inflow for examples.
61 */
62 /*----------------------------------------------------------------------------*/
63
64 /*============================================================================
65 * Local (user defined) function definitions
66 *============================================================================*/
67
68 /*----------------------------------------------------------------------------*/
69 /*!
70 * \brief Define parameters of synthetic turbulence at LES inflow.
71 */
72 /*----------------------------------------------------------------------------*/
73
74 void
cs_user_les_inflow_define(void)75 cs_user_les_inflow_define(void)
76 {
77 /*! [set_restart] */
78 cs_les_inflow_set_restart(false, /* allow_read */
79 false); /* allow_write */
80 /*! [set_restart] */
81
82 /* First synthetic turbulence inlet: the Batten Method is used
83 for boundary faces of zone "INLET_1" */
84
85 /*! [init_1] */
86 {
87 /* Number of turbulence structures */
88 int n_entities = 50;
89
90 /* Velocity, turbulent kinetic energy and dissipation rate */
91 cs_real_t vel_r[3] = {18.0, 0, 0};
92 cs_real_t k_r = 4.0;
93 cs_real_t eps_r = 4.0;
94
95 cs_les_inflow_add_inlet(CS_INFLOW_BATTEN,
96 false,
97 cs_boundary_zone_by_name("INLET_1"),
98 n_entities,
99 0, /* verbosity */
100 vel_r,
101 k_r,
102 eps_r);
103 }
104 /*! [init_1] */
105
106 /* Second synthetic turbulence inlet: the Synthetic Eddy Method
107 Batten Method is used for boundary faces of zone "INLET_2" */
108
109 /*! [init_2] */
110 {
111 /* Number of turbulence structures */
112 int n_entities = 50;
113
114 /* Velocity, turbulent kinetic energy and dissipation rate */
115 cs_real_t vel_r[3] = {12.0, 0, 0};
116 cs_real_t k_r = 3.0;
117 cs_real_t eps_r = 3.0;
118
119 cs_les_inflow_add_inlet(CS_INFLOW_SEM,
120 false, /* volume mode */
121 cs_boundary_zone_by_name("INLET_2"),
122 n_entities,
123 0, /* verbosity */
124 vel_r,
125 k_r,
126 eps_r);
127 }
128 /*! [init_2] */
129 }
130
131 /*----------------------------------------------------------------------------*/
132 /*!
133 * \brief Update of the characteristics of a given synthetic turbulence inlet.
134 *
135 * \param[in] zone pointer to associated boundary zone
136 * \param[out] vel_r reference mean velocity
137 * \param[out] k_r reference turbulent kinetic energy
138 * \param[out] eps_r reference turbulent dissipation
139 */
140 /*----------------------------------------------------------------------------*/
141
142 void
cs_user_les_inflow_update(const cs_zone_t * zone,cs_real_t vel_r[3],cs_real_t * k_r,cs_real_t * eps_r)143 cs_user_les_inflow_update(const cs_zone_t *zone,
144 cs_real_t vel_r[3],
145 cs_real_t *k_r,
146 cs_real_t *eps_r)
147 {
148 /*! [update_1] */
149 if (strcmp(zone->name, "INLET_1") == 0) {
150 /* Velocity, turbulent kinetic energy and dissipation rate */
151 vel_r[0] = 19.0;
152 vel_r[1] = 0.0;
153 vel_r[2] = 0.0;
154 *k_r = 5.0;
155 *eps_r = 5.0;
156 }
157 /*! [update_1] */
158 }
159
160 /*----------------------------------------------------------------------------*/
161 /*!
162 * \brief Definition of mean velocity, Reynolds stresses and dissipation rate
163 * for each boundary face of the given synthetic turbulence inlet.
164 *
165 * Rij components are ordered as usual: XX, YY, ZZ, XY, YZ, XZ
166 *
167 * Arrays are pre-initialized before this function is called
168 * (see \ref cs_user_les_inflow_define).
169 *
170 * \param[in] zone pointer to associated boundary zone
171 * \param[in, out] vel_l velocity a zone faces
172 * \param[in, out] rij_l reynods stresses at zone faces
173 * \param[in, out] eps_l reference turbulent dissipation
174 */
175 /*----------------------------------------------------------------------------*/
176
177 void
cs_user_les_inflow_advanced(const cs_zone_t * zone,cs_real_3_t vel_l[],cs_real_6_t rij_l[],cs_real_t eps_l[])178 cs_user_les_inflow_advanced(const cs_zone_t *zone,
179 cs_real_3_t vel_l[],
180 cs_real_6_t rij_l[],
181 cs_real_t eps_l[])
182 {
183 /* Example 1:
184 * - mean velocity, Reynolds stresses an dissipation are deduced
185 * from a wall law for the first synthetic turbulence inlet,
186 * - no modification of the statistics of the flow is provided for the
187 * second synthetic turbulence inlet */
188
189 /*! [example_1] */
190 if (strcmp(zone->name, "INLET_1") == 0) {
191
192 const cs_lnum_t n_faces = zone->n_elts;
193 const cs_lnum_t *face_ids = zone->elt_ids;
194
195 const cs_real_t d2_s3 = 2.0/3.0;
196
197 const cs_mesh_t *m = cs_glob_mesh;
198 const cs_real_3_t *b_face_cog
199 = (const cs_real_3_t *)cs_glob_mesh_quantities->b_face_cog;
200
201 const cs_real_t *cpro_mu = CS_F_(mu)->val;
202
203 const cs_real_t uref = cs_glob_turb_ref_values->uref;
204
205 /* Approximation of the friction velocity */
206 const cs_real_t utau = uref/20.0;
207
208 /* Reference length scale */
209 const cs_real_t h_ref = 1.0;
210
211 const cs_real_t kappa = cs_turb_xkappa;
212
213 for (cs_lnum_t i = 0; i < n_faces; i++) {
214
215 const cs_lnum_t f_id = face_ids[i];
216 const cs_lnum_t c_id = m->b_face_cells[f_id];
217
218 const cs_real_t x_mu = cpro_mu[c_id];
219 const cs_real_t r_fro = utau * h_ref / x_mu;
220
221 /* Dimensionless wall distance */
222 const cs_real_t yy = 1.0 - b_face_cog[f_id][1];
223 const cs_real_t yplus = yy / h_ref * r_fro;
224
225 /* Reichart laws (dimensionless) */
226 const cs_real_t uplus = log(1. + 0.4*yplus)/kappa
227 + 7.8*(1. - exp(-yplus/11.0)
228 - yplus/11.0*exp(-0.33*yplus));
229 const cs_real_t kplus = 0.07*yplus*yplus*exp(-yplus/8.0)
230 + (1.0 - exp(-yplus/20.0)) *4.5
231 / (1. + 4.0*yplus/r_fro);
232 const cs_real_t eplus = (1.0/kappa)
233 / pow(( cs_math_pow4(yplus)
234 + cs_math_pow4(15.)), 0.25);
235
236 /* Arrays are filled with dimensionnal stats */
237 vel_l[i][0] = uplus * utau;
238 vel_l[i][1] = 0.0;
239 vel_l[i][2] = 0.0;
240
241 rij_l[i][0] = d2_s3 * kplus * cs_math_pow2(utau);
242 rij_l[i][1] = d2_s3 * kplus * cs_math_pow2(utau);
243 rij_l[i][2] = d2_s3 * kplus * cs_math_pow2(utau);
244 rij_l[i][3] = 0;
245 rij_l[i][4] = 0;
246 rij_l[i][5] = 0;
247
248 eps_l[i] = eplus * cs_math_pow4(utau) / x_mu;
249 }
250 }
251 /*! [example_1] */
252
253 /* Example 2:
254 * - Reynolds stresses and dissipation at the inlet are computed
255 * using the turbulence intensity and standard laws for
256 * a circular pipe for the first synthetic turbulence inlet,
257 * - no modification of the statistics of the flow is provided for the
258 * second synthetic turbulence inlet */
259
260
261 /*! [example_2] */
262 if (strcmp(zone->name, "INLET_1") == 0) {
263
264 const cs_lnum_t n_faces = zone->n_elts;
265
266 const cs_real_t d2_s3 = 2.0/3.0;
267
268 for (cs_lnum_t i = 0; i < n_faces; i++) {
269
270 vel_l[i][0] = 1.1;
271 vel_l[i][1] = 1.1;
272 vel_l[i][2] = 1.1;
273
274 cs_real_t uref2 = cs_math_3_square_norm(vel_l[i]);
275 uref2 = cs_math_fmax(uref2, 1e-12);
276
277 /* Hydraulic diameter */
278 const cs_real_t x_hd = 0.075;
279
280 /* Turbulence intensity */
281 const cs_real_t x_ti = 0.02;
282
283 cs_real_t x_k = 1e-12;
284 cs_real_t x_eps = 1e-12;
285
286 cs_turbulence_bc_ke_turb_intensity(uref2,
287 x_ti,
288 x_hd,
289 &x_k,
290 &x_eps);
291
292 rij_l[i][0] = d2_s3 * x_k;
293 rij_l[i][1] = d2_s3 * x_k;
294 rij_l[i][2] = d2_s3 * x_k;
295 rij_l[i][3] = 0;
296 rij_l[i][4] = 0;
297 rij_l[i][5] = 0;
298
299 eps_l[i] = x_eps;
300 }
301 }
302 /*! [example_2] */
303 }
304
305 /*----------------------------------------------------------------------------*/
306
307 END_C_DECLS
308