1 /*============================================================================
2 * This function is called at the end of each time step, and has a very
3 * general purpose
4 * (i.e. anything that does not have another dedicated user function)
5 *============================================================================*/
6
7 /* VERS */
8
9 /*
10 This file is part of Code_Saturne, a general-purpose CFD tool.
11
12 Copyright (C) 1998-2021 EDF S.A.
13
14 This program is free software; you can redistribute it and/or modify it under
15 the terms of the GNU General Public License as published by the Free Software
16 Foundation; either version 2 of the License, or (at your option) any later
17 version.
18
19 This program is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
21 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
22 details.
23
24 You should have received a copy of the GNU General Public License along with
25 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
26 Street, Fifth Floor, Boston, MA 02110-1301, USA.
27 */
28
29 /*----------------------------------------------------------------------------*/
30
31 #include "cs_defs.h"
32
33 /*----------------------------------------------------------------------------
34 * Standard C library headers
35 *----------------------------------------------------------------------------*/
36
37 #include <assert.h>
38 #include <math.h>
39 #include <stdio.h>
40
41 #if defined(HAVE_MPI)
42 #include <mpi.h>
43 #endif
44
45 /*----------------------------------------------------------------------------
46 * PLE library headers
47 *----------------------------------------------------------------------------*/
48
49 #include <ple_coupling.h>
50
51 /*----------------------------------------------------------------------------
52 * Local headers
53 *----------------------------------------------------------------------------*/
54
55 #include "cs_headers.h"
56
57 /*----------------------------------------------------------------------------*/
58
59 BEGIN_C_DECLS
60
61 /*----------------------------------------------------------------------------*/
62 /*!
63 * \file cs_user_extra_operations-turbomachinery.c
64 *
65 * \brief This function is called at the end of each time step, and has a very
66 * general purpose (i.e. anything that does not have another dedicated
67 * user function).
68 *
69 * This example is a part of the \ref turbomachinery example.
70 */
71 /*----------------------------------------------------------------------------*/
72
73 /*============================================================================
74 * User function definitions
75 *============================================================================*/
76
77 /*----------------------------------------------------------------------------
78 * User function to find the closest cell and the associated rank of a point
79 * in a given rotation zone.
80 *
81 * parameters:
82 * domain <-- pointer to a cs_domain_t structure
83 * r <-- id of the rotation zone
84 * coords <-- point coordinates
85 * node --> cell id
86 * rank --> rank id
87 *----------------------------------------------------------------------------*/
88
89 static void
_findpt_r(cs_domain_t * domain,const cs_rotation_t * r,const cs_real_3_t coords,cs_lnum_t * node,int * rank)90 _findpt_r(cs_domain_t *domain,
91 const cs_rotation_t *r,
92 const cs_real_3_t coords,
93 cs_lnum_t *node,
94 int *rank)
95 {
96 cs_real_t d[3];
97
98 cs_lnum_t n_cells = domain->mesh-> n_cells;
99 cs_real_3_t *cell_cen = (cs_real_3_t *)domain->mesh_quantities->cell_cen;
100
101 *node = (int)(n_cells + 1)/2 - 1;
102
103 for (int i = 0; i < 3; i++)
104 d[i] = coords[i] - cell_cen[*node][i];
105 cs_real_t dis2mn = cs_math_3_square_norm(d);
106
107 /*! [extra_tbm_get_rotor] */
108
109 const int *rotor_num = NULL;
110 rotor_num = cs_turbomachinery_get_cell_rotor_num();
111
112 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
113 cs_rotation_t *rot = cs_glob_rotation + rotor_num[cell_id];
114 /* ... */
115
116 /*! [extra_tbm_get_rotor] */
117
118 if (r == rot) {
119 for (int i = 0; i < 3; i++)
120 d[i] = coords[i] - cell_cen[cell_id][i];
121 cs_real_t dis2 = cs_math_3_square_norm(d);
122 if (dis2 < dis2mn) {
123 *node = cell_id;
124 dis2mn = dis2;
125 }
126 }
127 }
128
129 if (cs_glob_rank_id >= 0)
130 cs_parall_min_id_rank_r(node, rank, dis2mn);
131 else
132 *rank = -1;
133 }
134
135 /*----------------------------------------------------------------------------*/
136 /*!
137 * \brief Example of extra operations for turbomachinery studies.
138 *
139 * \param[in, out] domain pointer to a cs_domain_t structure
140 */
141 /*----------------------------------------------------------------------------*/
142
143 void
cs_user_extra_operations(cs_domain_t * domain)144 cs_user_extra_operations(cs_domain_t *domain)
145 {
146
147 /*! [loc_def_init] */
148
149 /* Mesh-related variables */
150
151 const cs_lnum_t n_b_faces = domain->mesh->n_b_faces;
152
153 const cs_real_3_t *restrict cell_cen
154 = (const cs_real_3_t *restrict)domain->mesh_quantities->cell_cen;
155
156 /* 0. Initialization
157 ================= */
158
159 /* Local constants defintion */
160
161 const cs_real_t _PI = acos(-1.);
162 const cs_real_t _G = 9.81;
163
164 /* Flow properties */
165
166 double ro0 = cs_glob_fluid_properties->ro0;
167 cs_real_t flowrate = 0.238;
168
169 /*! [loc_def_init] */
170
171 /* Assert there is a rotation zone */
172
173 if (cs_glob_rotation == NULL) return;
174
175 /*! [extra_tbm_get_rotor_info] */
176
177 /* Access to a specific rotation zone (the first one) */
178
179 cs_lnum_t rotor_id = 1;
180 const cs_rotation_t *ref_rot = cs_glob_rotation + rotor_id;
181
182 /* Rotation zone parameters */
183
184 double omega = ref_rot->omega; /* angular velocity [rad/s] */
185 cs_real_3_t axis = {ref_rot->axis[0], /* rotation axis (normalized vector) */
186 ref_rot->axis[1], ref_rot->axis[2]};
187
188 /*! [extra_tbm_get_rotor_info] */
189
190 /* Example 1: calculation of the machinery characteristics
191 ======================================================= */
192
193 /*! [extra_tbm_post_util] */
194
195 /* Postprocessing of the couple: axial moment resulting from the stresses */
196 /* on the rotor blades */
197
198 cs_lnum_t n_elts;
199 cs_lnum_t *elt_list;
200
201 BFT_MALLOC(elt_list, n_b_faces, cs_lnum_t);
202 cs_selector_get_b_face_list("rotor_blades", &n_elts, elt_list);
203
204 cs_real_t c = cs_post_moment_of_force(n_elts, elt_list, axis);
205
206 BFT_FREE(elt_list);
207
208 /* Postprocessing of the head: total pressure increase through the machinery */
209
210 cs_real_t turbomachinery_head
211 = cs_post_turbomachinery_head
212 ("inlet", /* selection criteria at suction */
213 CS_MESH_LOCATION_BOUNDARY_FACES, /* associated mesh location */
214 "z > 2.725 and z < 2.775", /* selection criteria at discharge */
215 CS_MESH_LOCATION_CELLS); /* associated mesh location */
216
217 /*! [extra_tbm_post_util] */
218
219 cs_real_t manometric_head = turbomachinery_head/(ro0*_G);
220
221 cs_real_t power = c*omega;
222
223 cs_real_t efficiency = turbomachinery_head*flowrate/power;
224
225 /* Print in the log */
226
227 bft_printf("Turbomachinery characteristics:\n\n"
228 " %17s%17s%17s%17s\n"
229 " %17.9e%17.9e%17.9e%17.9e\n",
230 "Flowrate [m3/s]","Head [m]", "Power [W]", "Efficiency [1]",
231 flowrate, manometric_head, power, efficiency);
232
233 /* Example 2: extraction of a velocity profile in cylindrical coordinates
234 ====================================================================== */
235
236 /*! [extra_tbm_velocity_cylinder] */
237
238 if (domain->time_step->nt_cur == domain->time_step->nt_max){
239
240 cs_real_3_t *vel = (cs_real_3_t *)CS_F_(vel)->val;
241
242 cs_field_t *f_mom_wr = NULL, *f_mom_wt = NULL;
243
244 /* In case of transient rotor/stator computation, the moments for
245 the velocity in cylindrical coordinates are assumed to be defined
246 in the dedicated user file (cs_user_parameters.c, see the user example
247 for time moments definition) */
248
249 cs_turbomachinery_model_t tbm_model = cs_turbomachinery_get_model();
250
251 if (tbm_model == CS_TURBOMACHINERY_TRANSIENT) {
252 f_mom_wr = cs_time_moment_get_field(1);
253 f_mom_wt = cs_time_moment_get_field(2);
254 }
255
256 FILE *f1 = NULL, *f2 = NULL;
257
258 /* Only process of rank 0 (parallel) or -1 (scalar) writes to this file. */
259
260 if (cs_glob_rank_id <= 0) {
261
262 f1 = fopen("Wr_mean.dat","w");
263 fprintf(f1, "# %17s%17s\n", "angle", "Wr");
264
265 f2 = fopen("Wt_mean.dat","w");
266 fprintf(f2, "# %17s%17s\n", "angle", "Wt");
267 }
268
269 cs_real_t radius = 0.214;
270 cs_real_t axicoo = 2.e-2;
271
272 cs_lnum_t npoint = 360;
273 cs_lnum_t cell_id1 = -999;
274 int rank_id1 = -999;
275
276 for (cs_lnum_t point_id = 0; point_id < npoint; point_id++) {
277
278 cs_real_t xth0 = -_PI/npoint;
279 cs_real_3_t xyz = {radius*cos(xth0 - (double)point_id/npoint*2.*_PI),
280 radius*sin(xth0 - (double)point_id/npoint*2.*_PI),
281 axicoo};
282
283 cs_lnum_t cell_id, rank_id;
284
285 /* Find the closest cell in this rotor */
286 _findpt_r(domain, ref_rot, xyz, &cell_id, &rank_id);
287
288 if ((cell_id != cell_id1) || (rank_id != rank_id1)) {
289 cell_id1 = cell_id;
290 rank_id1 = rank_id;
291
292 cs_real_t xr, xtheta, xvr, xvt;
293
294 /* Set temporary variables for the process containing
295 * the point and then send it to other processes. */
296 if (cs_glob_rank_id == rank_id) {
297
298 /* Radius (first component of the x vector in cylindrical coords) */
299 cs_real_t xc[3];
300 cs_rotation_cyl_v(ref_rot, cell_cen[cell_id], cell_cen[cell_id], xc);
301 xr = xc[0];
302
303 /* Angle in [0, 2pi] */
304 double xthet1, xthet2;
305 xthet1 = acos((cell_cen[cell_id][0] - ref_rot->invariant[0])/xr);
306 xthet2 = asin((cell_cen[cell_id][1] - ref_rot->invariant[1])/xr);
307 if (xthet2 > 0)
308 xtheta = xthet1;
309 else if (xthet1 < _PI/2.)
310 xtheta = 2.*_PI + xthet2;
311 else
312 xtheta = _PI - xthet2;
313
314 if (tbm_model == CS_TURBOMACHINERY_FROZEN) {
315 cs_real_3_t xvc;
316 cs_rotation_cyl_v(ref_rot, cell_cen[cell_id], vel[cell_id], xvc);
317 xvr = xvc[0];
318 /* Relative tangential velocity */
319 xvt = xvc[1] - omega*xr;
320 }
321 else {
322 xvr = f_mom_wr->val[cell_id];
323 xvt = f_mom_wt->val[cell_id];
324 }
325
326 } else {
327 xtheta = 0.;
328 xvr = 0.;
329 xvt = 0.;
330 }
331
332 /* Broadcast to other ranks in parallel */
333 cs_parall_bcast(rank_id, 1, CS_DOUBLE, &xtheta);
334 cs_parall_bcast(rank_id, 1, CS_DOUBLE, &xvr);
335 cs_parall_bcast(rank_id, 1, CS_DOUBLE, &xvt);
336
337 if (cs_glob_rank_id <= 0) {
338 fprintf(f1," %17.9e%17.9e\n", xtheta, xvr);
339 fprintf(f2," %17.9e%17.9e\n", xtheta, xvt);
340 }
341 }
342
343 } /* end of loop on point_id */
344 if (cs_glob_rank_id <= 0) {
345 fclose(f1);
346 fclose(f2);
347 }
348
349 } /* end of if statement on current time step */
350
351 /*! [extra_tbm_velocity_cylinder] */
352 }
353
354 END_C_DECLS
355