1 /*============================================================================
2  * Code couplings definition with SYRTHES and Code_Saturne.
3  *
4  * 1) Define conjuguate heat transfer couplings with the SYRTHES code
5  * 2) Define couplings with other instances of Code_Saturne
6  *============================================================================*/
7 
8 /* VERS */
9 
10 /*
11   This file is part of Code_Saturne, a general-purpose CFD tool.
12 
13   Copyright (C) 1998-2021 EDF S.A.
14 
15   This program is free software; you can redistribute it and/or modify it under
16   the terms of the GNU General Public License as published by the Free Software
17   Foundation; either version 2 of the License, or (at your option) any later
18   version.
19 
20   This program is distributed in the hope that it will be useful, but WITHOUT
21   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
22   FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
23   details.
24 
25   You should have received a copy of the GNU General Public License along with
26   this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
27   Street, Fifth Floor, Boston, MA 02110-1301, USA.
28 */
29 
30 /*----------------------------------------------------------------------------*/
31 
32 #include "cs_defs.h"
33 
34 /*----------------------------------------------------------------------------
35  * Standard C library headers
36  *----------------------------------------------------------------------------*/
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_coupling-syrthes.c
57  *
58  * \brief Code couplings definition with SYRTHES and Code_Saturne.
59  *
60  * See \ref user_coupling for examples.
61  */
62 /*----------------------------------------------------------------------------*/
63 
64 /*============================================================================
65  * User function definitions
66  *============================================================================*/
67 
68 /*----------------------------------------------------------------------------*/
69 /*!
70  * \brief Define couplings with SYRTHES code.
71  *
72  * This is done by calling the \ref cs_syr_coupling_define function for each
73  * coupling to add.
74  */
75 /*----------------------------------------------------------------------------*/
76 
77 void
cs_user_syrthes_coupling(void)78 cs_user_syrthes_coupling(void)
79 {
80   /*! [coupling_syrthes_1] */
81   {
82     int  verbosity = 1, plot = 1;
83     float tolerance = 0.1;
84     bool allow_nonmatching = false;
85 
86     /*-------------------------------------------------------------------------
87      * Example 1:
88      *
89      * Boundary faces of group '3' coupled with instance named 'SYRTHES_01'.
90      *-------------------------------------------------------------------------*/
91 
92     cs_syr_coupling_define("SYRTHES_01",
93                            "3",               /* boundary criteria */
94                            NULL,              /* volume_criteria */
95                            ' ',               /* projection_axis */
96                            allow_nonmatching,
97                            tolerance,
98                            verbosity,
99                            plot);
100 
101   }
102   /*! [coupling_syrthes_1] */
103 
104   /*! [coupling_syrthes_2] */
105   {
106     int  verbosity = 1, plot = 1;
107     float tolerance = 0.1;
108     bool allow_nonmatching = false;
109 
110     /*-------------------------------------------------------------------------
111      * Example 2:
112      *
113      * Boundary faces of group 'Wall' coupled with 2D SYRTHES instance
114      * named 'SYRTHES_02'.
115      *-------------------------------------------------------------------------*/
116 
117     cs_syr_coupling_define("SYRTHES_02",
118                            "Wall",            /* boundary criteria */
119                            NULL,              /* volume_criteria */
120                            'z',               /* projection_axis */
121                            allow_nonmatching,
122                            tolerance,
123                            verbosity,
124                            plot);
125 
126   }
127   /*! [coupling_syrthes_2] */
128 
129   /*! [coupling_syrthes_3] */
130   {
131     int  verbosity = 1, plot = 1;
132     float tolerance = 0.1;
133     bool allow_nonmatching = false;
134 
135     /*-------------------------------------------------------------------------
136      * Example 3:
137      *
138      * Cells in box with corners (0, 0, 0) and (1, 1, 1) coupled with
139      * SYRTHES instance named 'Solid' (volume coupling).
140      *-------------------------------------------------------------------------*/
141 
142     cs_syr_coupling_define("Solid",
143                            NULL,                          /* boundary */
144                            "box[0., 0., 0., 1., 1., 1.]", /* volume */
145                            ' ',                           /* projection */
146                            allow_nonmatching,
147                            tolerance,
148                            verbosity,
149                            plot);
150   }
151   /*! [coupling_syrthes_3] */
152 
153   /* By default, conservativity forcing flag is switched off (value 0)
154      If one wants to switch on the conservativity forcing flag:
155 
156      cs_syr_coupling_set_conservativity(1);
157   */
158 
159   /* Only for a volume coupling:
160       By default, implicit treatment is done. You can switch to
161       an explicit treatment by using the following function:
162 
163      cs_syr_coupling_set_explicit_treatment();
164   */
165 
166 }
167 
168 /*----------------------------------------------------------------------------*/
169 /*!
170  * \brief Compute a volume exchange coefficient for SYRTHES couplings.
171  *
172  * \param[in]   coupling_id   Syrthes coupling id
173  * \param[in]   syrthes_name  name of associated Syrthes instance
174  * \param[in]   n_elts        number of associated cells
175  * \param[in]   elt_ids       associated cell ids
176  * \param[out]  h_vol         associated exchange coefficient (size: n_elts)
177  */
178 /*----------------------------------------------------------------------------*/
179 
180 void
cs_user_syrthes_coupling_volume_h(int coupling_id,const char * syrthes_name,cs_lnum_t n_elts,const cs_lnum_t elt_ids[],cs_real_t h_vol[])181 cs_user_syrthes_coupling_volume_h(int               coupling_id,
182                                   const char       *syrthes_name,
183                                   cs_lnum_t         n_elts,
184                                   const cs_lnum_t   elt_ids[],
185                                   cs_real_t         h_vol[])
186 {
187   CS_NO_WARN_IF_UNUSED(coupling_id);
188   CS_NO_WARN_IF_UNUSED(syrthes_name);
189 
190   /* Example 1 of the computation of a constant volumic exchange coefficient
191      ----------------------------------------------------------------------- */
192 
193   /*! [example_1] */
194   cs_real_t hvol_cst = 1.0e6;
195 
196   for (cs_lnum_t i = 0; i < n_elts; i++) {
197     h_vol[i] = hvol_cst;
198   }
199   /*! [example_1] */
200 
201   /* Example 2 of the computation of a volumic exchange coefficient
202    * --------------------------------------------------------------
203    *
204    *  hvol(iel) =  hsurf(iel) * exchange_surface_by_unit_vol
205    *
206    *  with: hsurf = Nusselt * lambda / L
207    *
208    *  lambda is the thermal conductivity coefficient
209    *   L is a characteristic length
210    *
211    *  Nusselt is computed by means of the Colburn correlation
212    *
213    *  Nu = 0.023 * Re^(0.8) * Pr^(1/3)
214    *
215    *  Re is the Reynolds number and Pr is the Prandtl number
216    */
217 
218   /*! [example_2_init] */
219   const cs_real_3_t  *cvar_vel = (const cs_real_3_t *)CS_F_(vel)->val;
220   const cs_real_t  *cpro_rom = (const cs_real_t *)CS_F_(rho)->val;
221   const cs_real_t  *cpro_mu = (const cs_real_t *)CS_F_(mu)->val;
222 
223   const cs_real_t  cp0 = cs_glob_fluid_properties->cp0;
224   cs_real_t  visls_0 = -1;
225 
226   const cs_real_t  *cpro_cp = NULL, *cpro_viscls = NULL;
227   cs_lnum_t  cp_step = 0, viscls_step = 0;
228 
229   if (CS_F_(cp) != NULL) {
230     cpro_cp = (const cs_real_t *)CS_F_(cp)->val;
231     cp_step = 1;
232   }
233   else {
234     cpro_cp = &cp0;
235   }
236 
237   /* Get thermal field and associated diffusivity
238      (temperature only handled here) */
239 
240   const cs_field_t *fth = cs_thermal_model_field();
241 
242   const int viscl_id = cs_field_get_key_int(fth,
243                                             cs_field_key_id("diffusivity_id"));
244 
245   if (viscl_id > -1) {
246     cpro_viscls = (const cs_real_t *)cs_field_by_id(viscl_id)->val;
247     viscls_step = 1;
248   }
249   else {
250     visls_0 = cs_field_get_key_double(fth, cs_field_key_id("diffusivity_ref"));
251     cpro_viscls = &visls_0;
252   }
253 
254   const int is_temperature
255     = cs_field_get_key_int(fth,
256                            cs_field_key_id("is_temperature"));
257   /*! [example_2_init] */
258 
259   /*! [example_2] */
260   cs_real_t sexcvo = 36.18;  /* Surface area where exchanges take
261                                 place by unit of volume */
262   cs_real_t l0 = 0.03;       /* Characteristic length */
263 
264   for (cs_lnum_t i = 0; i < n_elts; i++) {
265 
266     cs_lnum_t c_id = elt_ids[i];
267 
268     cs_real_t rho = cpro_rom[c_id];
269     cs_real_t mu = cpro_mu[c_id];
270     cs_real_t cp = cpro_cp[c_id*cp_step];
271 
272     cs_real_t lambda, lambda_over_cp;
273 
274     if (is_temperature) {
275       lambda = cpro_viscls[c_id*viscls_step];
276       lambda_over_cp = lambda / cp;
277     }
278     else {
279       lambda_over_cp = cpro_viscls[c_id*viscls_step];
280       lambda =  lambda_over_cp * cp;
281     }
282 
283     /* Compute a local molecular Prandtl **(1/3) */
284 
285     cs_real_t  pr = mu / lambda_over_cp;
286 
287     /* Compute a local Reynolds number */
288 
289     cs_real_t uloc = cs_math_3_norm(cvar_vel[c_id]);
290     cs_real_t re = fmax(uloc*rho*l0/mu, 1.);  /* To avoid division by zero */
291 
292     /* Compute Nusselt number using Colburn correlation */
293 
294     cs_real_t nu = 0.023 * pow(re, 0.8) * pow(pr, 1./3.);
295     cs_real_t h_corr = nu * lambda / l0;
296 
297     /* Compute hvol */
298     h_vol[i] = h_corr * sexcvo;
299   }
300   /*! [example_2] */
301 }
302 
303 /*----------------------------------------------------------------------------*/
304 
305 END_C_DECLS
306