TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Paroi_2couches_VDF.cpp
1/****************************************************************************
2* Copyright (c) 2017, CEA
3* All rights reserved.
4*
5* Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
6* 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
7* 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
8* 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
9*
10* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
11* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
12* OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
13*
14*****************************************************************************/
15
16#include <Paroi_2couches_VDF.h>
17#include <Modele_turbulence_hyd_K_Eps_2_Couches.h>
18#include <Fluide_base.h>
19#include <Champ_Uniforme.h>
20#include <Domaine_Cl_VDF.h>
21#include <Dirichlet_paroi_fixe.h>
22#include <Dirichlet_paroi_defilante.h>
23#include <Paroi_std_hyd_VDF.h>
24#include <Debog.h>
25#include <EcrFicPartage.h>
26#include <Param.h>
27
28Implemente_instanciable_sans_constructeur(Paroi_2couches_VDF,"loi_paroi_2_couches_VDF",Paroi_hyd_base_VDF);
29
30// printOn()
31/////
32
33Sortie& Paroi_2couches_VDF::printOn(Sortie& s) const
34{
35 return s << que_suis_je() << " " << le_nom();
36}
37
38//// readOn
39//
40
41Entree& Paroi_2couches_VDF::readOn(Entree& s)
42{
43 Paroi_hyd_base_VDF::readOn(s);
44 return s ;
45}
46
52
53
54/////////////////////////////////////////////////////////////////////
55//
56// Implementation des fonctions de la classe Paroi_2couches_VDF
57//
58/////////////////////////////////////////////////////////////////////
59
60int Paroi_2couches_VDF::init_lois_paroi()
61{
62 uplus_.resize(le_dom_dis_->nb_faces_bord());
63 init_lois_paroi_();
64
65 return init_lois_paroi_hydraulique();
66}
67
68// Remplissage de la table
69
70int Paroi_2couches_VDF::init_lois_paroi_hydraulique()
71{
72 Cmu_ = mon_modele_turb_hyd->get_Cmu();
73 init_lois_paroi_hydraulique_();
74 return 1;
75}
76
77int Paroi_2couches_VDF::calculer_hyd(DoubleTab& tab_k_eps)
78{
79 Cerr << "Paroi_2couches_VDF::calculer_hyd(DoubleTab& tab_k_eps) : Erreur on ne devrait pas etre ici!!!" << finl;
80 return 1;
81}
82
83int Paroi_2couches_VDF::calculer_hyd_BiK(DoubleTab& tab_k,DoubleTab& tab_eps)
84{
85 Cerr << "Paroi_2couches_VDF::calculer_hyd_BiK(DoubleTab& tab_k,DoubleTab& tab_eps) : Erreur on ne devrait pas etre ici!!!" << finl;
86 return 1;
87}
88
89int Paroi_2couches_VDF::calculer_hyd(DoubleTab& tab_nu_t,DoubleTab& tab_k_eps)
90{
91 //Cerr << "Dans Paroi_2couches_VDF::calculer_hyd(DoubleTab& tab_nu_t,DoubleTab& tab_k) : Ok ! " << finl;
92 const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_dis_.valeur());
93 const IntVect& orientation = domaine_VDF.orientation();
94 const IntTab& face_voisins = domaine_VDF.face_voisins();
95 const IntTab& elem_faces = domaine_VDF.elem_faces();
96 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
97 const Fluide_base& le_fluide = ref_cast(Fluide_base, eqn_hydr.milieu());
98 const Champ_Don_base& ch_visco_cin = le_fluide.viscosite_cinematique();
99 const DoubleVect& vit = eqn_hydr.inconnue().valeurs();
100 const DoubleTab& tab_visco = ch_visco_cin.valeurs();
101 double visco=-1.;
102 Modele_turbulence_hyd_K_Eps_2_Couches& mod_2couches = ref_cast(Modele_turbulence_hyd_K_Eps_2_Couches,mon_modele_turb_hyd.valeur());
103 const int nbcouches = mod_2couches.get_nbcouches();
104
105 Loi_2couches_base& loi2couches = ref_cast(Transport_K_KEps,mod_2couches.get_set_eq_transport()).loi2couches();
106
107 int valswitch;
108
109 const int typeswitch = mod_2couches.get_switch();
110 if ( typeswitch == 0 )
111 valswitch = mod_2couches.get_yswitch();
112 else
113 valswitch = mod_2couches.get_nutswitch();
114
115 int l_unif;
116 if (sub_type(Champ_Uniforme,ch_visco_cin))
117 {
118 visco = std::max(tab_visco(0,0),DMINFLOAT);
119 l_unif = 1;
120 }
121 else
122 l_unif = 0;
123
124 // preparer_calcul_hyd(tab);
125 if ((!l_unif) && (tab_visco.local_min_vect()<DMINFLOAT))
126 // on ne doit pas changer tab_visco ici !
127 {
128 Cerr << "In Paroi_2couches_VDF::calculer_hyd : visco = " << tab_visco.local_min_vect() << " <= 0 ? " << finl;
129 throw;
130 }
131 // tab_visco+=DMINFLOAT;
132
133 int ndeb,nfin,face_courante,elem_courant;
134 int elem,ori;
135 double norm_v;
136 double dist;
137 double u_plus_d_plus,d_visco,y_etoile, critere_switch;
138 double val,val1,val2;
139
140 // Boucle sur les bords
141 //Cerr << " nb frontieres = " << domaine_VDF.nb_front_Cl() << finl;
142 for (int n_bord=0; n_bord<domaine_VDF.nb_front_Cl(); n_bord++)
143 {
144
145 // pour chaque condition limite on regarde son type
146 // On applique les lois de paroi uniquement
147 // aux voisinages des parois
148
149 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
150 // Cerr << "n = " << n_bord << " = " << la_cl.valeur() << finl;
151 if (sub_type(Dirichlet_paroi_fixe,la_cl.valeur()) )
152 {
153 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
154 ndeb = le_bord.num_premiere_face();
155 nfin = ndeb + le_bord.nb_faces();
156
157 if (dimension == 2 )
158 for (int num_face=ndeb; num_face<nfin; num_face++)
159 {
160 ori = orientation(num_face);
161 if ( (elem =face_voisins(num_face,0)) != -1)
162 norm_v=norm_2D_vit(vit,elem,ori,domaine_VDF,val);
163 else
164 {
165 elem = face_voisins(num_face,1);
166 norm_v=norm_2D_vit(vit,elem,ori,domaine_VDF,val);
167 }
168 dist=domaine_VDF.distance_normale(num_face);
169 if (l_unif)
170 d_visco = visco;
171 else
172 d_visco = tab_visco[elem];
173 u_plus_d_plus = norm_v*dist/d_visco;
174
175
176 face_courante = num_face;
177 elem_courant = elem;
178
179 // Loi pour Epsilon et visco_turb
180 for(int icouche=0; (icouche<nbcouches) && (elem_courant != -1); icouche++)
181 {
182 //Cerr << "***********Couche " << icouche << "**************" << finl;
183 double kSqRt = sqrt(tab_k_eps(elem_courant,0));
184 y_etoile = kSqRt*dist/d_visco;
185 if (typeswitch == 0) //selon le critere de switch choisit:y* ou nu_t
186 critere_switch = y_etoile;
187 else
188 critere_switch = tab_nu_t(elem_courant)/d_visco;
189 if (critere_switch < valswitch) // si on est en proche paroi
190 {
191 double vvSqRt,Leps,Lmu;
192 loi2couches.LepsLmu(tab_k_eps(elem_courant,0),d_visco,dist,y_etoile,Leps,Lmu,vvSqRt);
193 if ( (Leps != 0) && (Lmu!=0.) )
194 {
195 tab_k_eps(elem_courant,1) = tab_k_eps(elem_courant,0)*vvSqRt/Leps;
196 tab_nu_t(elem_courant) = vvSqRt*Lmu;
197 }
198 else
199 {
200 tab_k_eps(elem_courant,1) = 0;
201 tab_nu_t(elem_courant) = 0;
202 }
203 }
204 if ( elem_faces(elem_courant,0) == face_courante)
205 face_courante = elem_faces(elem_courant,2);
206 else if ( elem_faces(elem_courant,1) == face_courante)
207 face_courante = elem_faces(elem_courant,3);
208 else if ( elem_faces(elem_courant,2) == face_courante)
209 face_courante = elem_faces(elem_courant,0);
210 else if ( elem_faces(elem_courant,3) == face_courante)
211 face_courante = elem_faces(elem_courant,1);
212 //Cerr << "elem courant " << elem_courant << finl;
213 //Cerr << "nve elem = " << face_voisins(face_courante,0) << finl;
214 if ( face_voisins(face_courante,0) != elem_courant)
215 elem_courant = face_voisins(face_courante,0);
216 else
217 elem_courant = face_voisins(face_courante,1);
218 //Cerr << "dist = " << dist << finl;
219 dist+=domaine_VDF.distance_normale(face_courante);
220
221 }
222
223 // Calcul de u* et des grandeurs turbulentes:
224
225 double valmin = table_hyd.val(5);
226 double valmax = table_hyd.val(40);
227 if (u_plus_d_plus <= valmin)
228 tab_u_star_(num_face) = sqrt(norm_v*d_visco/dist);
229 else if ((u_plus_d_plus > valmin) && (u_plus_d_plus < valmax))
230 {
231 double d_plus;
232 calculer_u_star_sous_couche_tampon(d_plus,u_plus_d_plus,d_visco,dist,num_face);
233 }
234 else if (u_plus_d_plus >= valmax)
235 calculer_u_star_sous_couche_log(norm_v,d_visco,dist,num_face);
236
237
238 // Calcul de la contrainte tangentielle
239
240 Cisaillement_paroi_(num_face,1-ori) = tab_u_star(num_face)*tab_u_star(num_face)*val;
241
242 // Calcul de u+ d+
243 calculer_uplus_dplus(uplus_, tab_d_plus_, tab_u_star_, num_face, dist, d_visco, norm_v) ;
244 }
245
246 else if (dimension == 3)
247 for (int num_face=ndeb; num_face<nfin; num_face++)
248 {
249 ori = orientation(num_face);
250 if ( (elem = face_voisins(num_face,0)) != -1)
251 norm_v=norm_3D_vit(vit,elem,ori,domaine_VDF,val1,val2);
252 else
253 {
254 elem = face_voisins(num_face,1);
255 norm_v=norm_3D_vit(vit,elem,ori,domaine_VDF,val1,val2);
256 }
257 dist = domaine_VDF.distance_normale(num_face);
258 if (l_unif)
259 d_visco = visco;
260 else
261 d_visco = tab_visco[elem];
262 u_plus_d_plus = norm_v*dist/d_visco;
263
264 face_courante = num_face;
265 elem_courant = elem;
266
267 // Loi pour Epsilon et visco_turb
268 for(int icouche=0; (icouche<nbcouches) && (elem_courant != -1); icouche++)
269 {
270 //Cerr << "***********Couche " << icouche << "**************" << finl;
271 double kSqRt = sqrt(tab_k_eps(elem_courant,0));
272 y_etoile = kSqRt*dist/d_visco;
273 if (typeswitch == 0) //selon le critere de switch choisit:y* ou nu_t
274 critere_switch = y_etoile;
275 else
276 critere_switch = tab_nu_t(elem_courant)/d_visco;
277 if (critere_switch < valswitch) // si on est en proche paroi
278 {
279 double vvSqRt,Leps,Lmu;
280 loi2couches.LepsLmu(tab_k_eps(elem_courant,0),d_visco,dist,y_etoile,Leps,Lmu,vvSqRt);
281 if ( (Leps != 0) && (Lmu!=0.) )
282 {
283 tab_k_eps(elem_courant,1) = tab_k_eps(elem_courant,0)*vvSqRt/Leps;
284 tab_nu_t(elem_courant) = vvSqRt*Lmu;
285 }
286 else
287 {
288 tab_k_eps(elem_courant,1) = 0;
289 tab_nu_t(elem_courant) = 0;
290 }
291 }
292
293 //Cerr << "face cournate " << face_courante << finl;
294 //Cerr << "face 0 " << elem_faces(elem_courant,0) << finl;
295 if ( elem_faces(elem_courant,0) == face_courante)
296 face_courante = elem_faces(elem_courant,3);
297 else if ( elem_faces(elem_courant,1) == face_courante)
298 face_courante = elem_faces(elem_courant,4);
299 else if ( elem_faces(elem_courant,2) == face_courante)
300 face_courante = elem_faces(elem_courant,5);
301 else if ( elem_faces(elem_courant,3) == face_courante)
302 face_courante = elem_faces(elem_courant,0);
303 else if ( elem_faces(elem_courant,4) == face_courante)
304 face_courante = elem_faces(elem_courant,1);
305 else if ( elem_faces(elem_courant,5) == face_courante)
306 face_courante = elem_faces(elem_courant,2);
307
308 //Cerr << "elem courant " << elem_courant << finl;
309 //Cerr << "nve elem = " << face_voisins(face_courante,0) << finl;
310 if ( face_voisins(face_courante,0) != elem_courant)
311 elem_courant = face_voisins(face_courante,0);
312 else
313 elem_courant = face_voisins(face_courante,1);
314 //Cerr << "dist = " << dist << finl;
315 dist+=domaine_VDF.distance_normale(face_courante);
316
317 }
318
319 // Calcul de u* et des grandeurs turbulentes:
320
321 double valmin = table_hyd.val(5);
322 double valmax = table_hyd.val(40);
323
324 if (u_plus_d_plus <= valmin)
325 tab_u_star_(num_face) = sqrt(norm_v*d_visco/dist);
326 else if ((u_plus_d_plus > valmin) && (u_plus_d_plus < valmax))
327 {
328 double d_plus;
329 calculer_u_star_sous_couche_tampon(d_plus,u_plus_d_plus,d_visco,dist,num_face);
330 }
331 else if (u_plus_d_plus >= valmax)
332 calculer_u_star_sous_couche_log(norm_v,d_visco,dist,num_face);
333
334 // Calcul des deux composantes de la contrainte tangentielle:
335
336 double vit_frot = tab_u_star(num_face)*tab_u_star(num_face);
337
338 if (ori == 0)
339 {
340 Cisaillement_paroi_(num_face,1) = vit_frot*val1;
341 Cisaillement_paroi_(num_face,2) = vit_frot*val2;
342 }
343 else if (ori == 1)
344 {
345 Cisaillement_paroi_(num_face,0) = vit_frot*val1;
346 Cisaillement_paroi_(num_face,2) = vit_frot*val2;
347 }
348 else
349 {
350 Cisaillement_paroi_(num_face,0) = vit_frot*val1;
351 Cisaillement_paroi_(num_face,1) = vit_frot*val2;
352 }
353 // Calcul de u+ d+
354 calculer_uplus_dplus(uplus_, tab_d_plus_, tab_u_star_, num_face, dist, d_visco, norm_v) ;
355 }
356 }
357
358 else if (sub_type(Dirichlet_paroi_defilante,la_cl.valeur()) )
359 {
360 const Dirichlet_paroi_defilante& cl_diri = ref_cast(Dirichlet_paroi_defilante,la_cl.valeur());
361 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
362 ndeb = le_bord.num_premiere_face();
363 nfin = ndeb + le_bord.nb_faces();
364 DoubleTab vitesse_imposee_face_bord(le_bord.nb_faces(),dimension);
365 for (int face=ndeb; face<nfin; face++)
366 for (int k=0; k<dimension; k++)
367 vitesse_imposee_face_bord(face-ndeb,k) = cl_diri.val_imp(face-ndeb,k);
368
369 if (dimension == 2 )
370
371 // On calcule au centre des mailles de bord
372 // la vitesse parallele a l'axe de la face de bord
373 //
374 // On calcule la norme de cette vitesse
375
376 for (int num_face=ndeb; num_face<nfin; num_face++)
377 {
378 ori = orientation(num_face);
379 int rang = num_face-ndeb;
380 if ( (elem = face_voisins(num_face,0)) != -1)
381 norm_v=norm_2D_vit(vit,elem,ori,domaine_VDF,vitesse_imposee_face_bord(rang,0),
382 vitesse_imposee_face_bord(rang,1),val);
383 else
384 {
385 elem = face_voisins(num_face,1);
386 norm_v=norm_2D_vit(vit,elem,ori,domaine_VDF,vitesse_imposee_face_bord(rang,0),
387 vitesse_imposee_face_bord(rang,1),val);
388 }
389
390 dist = domaine_VDF.distance_normale(num_face);
391 if (l_unif)
392 d_visco = visco;
393 else
394 d_visco = tab_visco[elem];
395
396 u_plus_d_plus = norm_v*dist/d_visco;
397
398 face_courante = num_face;
399 elem_courant = elem;
400
401 // Loi pour Epsilon et visco_turb
402 for(int icouche=0; (icouche<nbcouches) && (elem_courant != -1); icouche++)
403 {
404 //Cerr << "***********Couche " << icouche << "**************" << finl;
405 double kSqRt = sqrt(tab_k_eps(elem_courant,0));
406 y_etoile = kSqRt*dist/d_visco;
407 if (typeswitch == 0) //selon le critere de switch choisit:y* ou nu_t
408 critere_switch = y_etoile;
409 else
410 critere_switch = tab_nu_t(elem_courant)/d_visco;
411 if (critere_switch < valswitch) // si on est en proche paroi
412 {
413 double vvSqRt,Leps,Lmu;
414 loi2couches.LepsLmu(tab_k_eps(elem_courant,0),d_visco,dist,y_etoile,Leps,Lmu,vvSqRt);
415 if ( (Leps != 0) && (Lmu!=0.) )
416 {
417 tab_k_eps(elem_courant,1) = tab_k_eps(elem_courant,0)*vvSqRt/Leps;
418 tab_nu_t(elem_courant) = vvSqRt*Lmu;
419 }
420 else
421 {
422 tab_k_eps(elem_courant,1) = 0;
423 tab_nu_t(elem_courant) = 0;
424 }
425 }
426 //Cerr << "face courante " << face_courante << finl;
427 //Cerr << "face 0 " << elem_faces(elem_courant,0) << finl;
428 if ( elem_faces(elem_courant,0) == face_courante)
429 face_courante = elem_faces(elem_courant,2);
430 else if ( elem_faces(elem_courant,1) == face_courante)
431 face_courante = elem_faces(elem_courant,3);
432 else if ( elem_faces(elem_courant,2) == face_courante)
433 face_courante = elem_faces(elem_courant,0);
434 else if ( elem_faces(elem_courant,3) == face_courante)
435 face_courante = elem_faces(elem_courant,1);
436 //Cerr << "elem courant " << elem_courant << finl;
437 //Cerr << "nve elem = " << face_voisins(face_courante,0) << finl;
438 if ( face_voisins(face_courante,0) != elem_courant)
439 elem_courant = face_voisins(face_courante,0);
440 else
441 elem_courant = face_voisins(face_courante,1);
442 //Cerr << "dist = " << dist << finl;
443 dist+=domaine_VDF.distance_normale(face_courante);
444
445 }
446 // Calcul de u* et des grandeurs turbulentes:
447
448 // calculer_local(u_plus_d_plus,d_visco,tab_nu_t,tab_k,norm_v,dist,elem,num_face);
449 double valmin = table_hyd.val(5);
450 double valmax = table_hyd.val(40);
451
452 if (u_plus_d_plus <= valmin)
453 tab_u_star_(num_face) = sqrt(norm_v*d_visco/dist);
454 else if ((u_plus_d_plus > valmin) && (u_plus_d_plus < valmax))
455 {
456 double d_plus;
457 calculer_u_star_sous_couche_tampon(d_plus,u_plus_d_plus,d_visco,dist,num_face);
458 }
459 else if (u_plus_d_plus >= valmax)
460 calculer_u_star_sous_couche_log(norm_v,d_visco,dist,num_face);
461
462 // Calcul de la contrainte tangentielle a la paroi:
463
464 Cisaillement_paroi_(num_face,1-ori) = tab_u_star(num_face)*tab_u_star(num_face)*val;
465 // Calcul de u+ d+
466 calculer_uplus_dplus(uplus_, tab_d_plus_, tab_u_star_, num_face, dist, d_visco, norm_v) ;
467 }
468
469 else if (dimension == 3)
470 // On calcule au centre des mailles de bord
471 // la vitesse dans le plan parallele a celui de la face de bord
472 // On calcule la norme de cette vitesse
473
474 for (int num_face=ndeb; num_face<nfin; num_face++)
475 {
476 ori = orientation(num_face);
477 int rang = num_face-ndeb;
478 if ( (elem = face_voisins(num_face,0)) != -1)
479 norm_v=norm_3D_vit(vit,elem,ori,domaine_VDF,vitesse_imposee_face_bord(rang,0),
480 vitesse_imposee_face_bord(rang,1),
481 vitesse_imposee_face_bord(rang,2),val1,val2);
482 else
483 {
484 elem = face_voisins(num_face,1);
485 norm_v=norm_3D_vit(vit,elem,ori,domaine_VDF,vitesse_imposee_face_bord(rang,0),
486 vitesse_imposee_face_bord(rang,1),
487 vitesse_imposee_face_bord(rang,2),val1,val2);
488 }
489 dist = domaine_VDF.distance_normale(num_face);
490 if (l_unif)
491 d_visco = visco;
492 else
493 d_visco = tab_visco[elem];
494
495 u_plus_d_plus = norm_v*dist/d_visco;
496
497 face_courante = num_face;
498 elem_courant = elem;
499
500 // Loi pour Epsilon et visco_turb
501 for(int icouche=0; (icouche<nbcouches) && (elem_courant != -1); icouche++)
502 {
503 //Cerr << "***********Couche " << icouche << "**************" << finl;
504 double kSqRt = sqrt(tab_k_eps(elem_courant,0));
505 y_etoile = kSqRt*dist/d_visco;
506 if (typeswitch == 0) //selon le critere de switch choisit:y* ou nu_t
507 critere_switch = y_etoile;
508 else
509 critere_switch = tab_nu_t(elem_courant)/d_visco;
510 if (critere_switch < valswitch) // si on est en proche paroi
511 {
512 double vvSqRt,Leps,Lmu;
513 loi2couches.LepsLmu(tab_k_eps(elem_courant,0),d_visco,dist,y_etoile,Leps,Lmu,vvSqRt);
514 if ( (Leps != 0) && (Lmu!=0.) )
515 {
516 tab_k_eps(elem_courant,1) = tab_k_eps(elem_courant,0)*vvSqRt/Leps;
517 tab_nu_t(elem_courant) = vvSqRt*Lmu;
518 }
519 else
520 {
521 tab_k_eps(elem_courant,1) = 0;
522 tab_nu_t(elem_courant) = 0;
523 }
524 }
525 //Cerr << "face cournate " << face_courante << finl;
526 //Cerr << "face 0 " << elem_faces(elem_courant,0) << finl;
527 if ( elem_faces(elem_courant,0) == face_courante)
528 face_courante = elem_faces(elem_courant,3);
529 else if ( elem_faces(elem_courant,1) == face_courante)
530 face_courante = elem_faces(elem_courant,4);
531 else if ( elem_faces(elem_courant,2) == face_courante)
532 face_courante = elem_faces(elem_courant,5);
533 else if ( elem_faces(elem_courant,3) == face_courante)
534 face_courante = elem_faces(elem_courant,0);
535 else if ( elem_faces(elem_courant,4) == face_courante)
536 face_courante = elem_faces(elem_courant,1);
537 else if ( elem_faces(elem_courant,5) == face_courante)
538 face_courante = elem_faces(elem_courant,2);
539
540 //Cerr << "elem courant " << elem_courant << finl;
541 //Cerr << "nve elem = " << face_voisins(face_courante,0) << finl;
542 if ( face_voisins(face_courante,0) != elem_courant)
543 elem_courant = face_voisins(face_courante,0);
544 else
545 elem_courant = face_voisins(face_courante,1);
546 //Cerr << "dist = " << dist << finl;
547 dist+=domaine_VDF.distance_normale(face_courante);
548
549 }
550 // Calcul de u* et des grandeurs turbulentes:
551
552 // calculer_local(u_plus_d_plus,d_visco,tab_nu_t,tab_k,norm_v,dist,elem,num_face);
553 double valmin = table_hyd.val(5);
554 double valmax = table_hyd.val(40);
555
556 if (u_plus_d_plus <= valmin)
557 tab_u_star_(num_face) = sqrt(norm_v*d_visco/dist);
558 else if ((u_plus_d_plus > valmin) && (u_plus_d_plus < valmax))
559 {
560 double d_plus;
561 calculer_u_star_sous_couche_tampon(d_plus,u_plus_d_plus,d_visco,dist,num_face);
562 }
563 else if (u_plus_d_plus >= valmax)
564 calculer_u_star_sous_couche_log(norm_v,d_visco,dist,num_face);
565
566 // Calcul des deux composantes de la contrainte tangentielle:
567
568 double vit_frot = tab_u_star(num_face)*tab_u_star(num_face);
569
570 if (ori == 0)
571 {
572 Cisaillement_paroi_(num_face,1) = vit_frot*val1;
573 Cisaillement_paroi_(num_face,2) = vit_frot*val2;
574 }
575 else if (ori == 1)
576 {
577 Cisaillement_paroi_(num_face,0) = vit_frot*val1;
578 Cisaillement_paroi_(num_face,2) = vit_frot*val2;
579 }
580 else
581 {
582 Cisaillement_paroi_(num_face,0) = vit_frot*val1;
583 Cisaillement_paroi_(num_face,1) = vit_frot*val2;
584 }
585 // Calcul de u+ d+
586 calculer_uplus_dplus(uplus_, tab_d_plus_, tab_u_star_, num_face, dist, d_visco, norm_v) ;
587 }
588
589 }
590 }
591 Cisaillement_paroi_.echange_espace_virtuel();
592 tab_nu_t.echange_espace_virtuel();
593 tab_k_eps.echange_espace_virtuel();
594
595 Debog::verifier("Cisaillement_paroi : ", Cisaillement_paroi_);
596 Debog::verifier("tab_nu_t : ", tab_nu_t);
597 Debog::verifier("tab_k_eps : ", tab_k_eps);
598
599 return 1;
600}
601
602// Calcul de u+ d+
603void Paroi_2couches_VDF::calculer_uplus_dplus(DoubleVect& uplus, DoubleVect& dplus, DoubleVect& tab_ustar,
604 int num_face, double dist, double d_visco, double norm_v)
605{
606 double ustar=tab_ustar(num_face);
607 dplus(num_face)=ustar*dist/d_visco;
608 if (ustar != 0)
609 uplus(num_face)=norm_v/ustar;
610}
611
612
613int Paroi_2couches_VDF::calculer_u_star_sous_couche_tampon(double& d_plus,double u_plus_d_plus,
614 double d_visco,double dist,int face)
615{
616
617 // Calcul de d_plus solution de table_hyd(inconnue) = u_plus_d_plus
618 // puis calcul de u* dans la sous-couche tampon : u* = nu*d+/dist
619 double d_plus_min = 5;
620 double d_plus_max = 40;
621 double valeur,u_star;
622 double epsilon = 1.e-12;
623 double gauche = table_hyd.val(d_plus_min);
624 double droite = table_hyd.val(d_plus_max);
625 double deriv;
626 if (gauche == droite)
627 d_plus = d_plus_min;
628 else
629 {
630 while(1)
631 {
632 deriv = (droite-gauche)/(d_plus_max-d_plus_min);
633 d_plus = d_plus_min + (u_plus_d_plus - gauche)/deriv;
634 valeur = table_hyd.val(d_plus);
635 if(std::fabs(valeur-u_plus_d_plus) < epsilon)
636 {
637 u_star = (d_visco*d_plus)/dist;
638 tab_u_star_(face) = u_star;
639 return 1;
640 }
641 if(valeur>u_plus_d_plus)
642 {
643 droite=valeur;
644 d_plus_max=d_plus;
645 }
646 else
647 {
648 gauche=valeur;
649 d_plus_min=d_plus;
650 }
651 }
652 }
653 return -1;
654}
655
656int Paroi_2couches_VDF::calculer_u_star_sous_couche_log(double norm_vit,double d_visco,
657 double dist, int face)
658{
659 // Dans la sous couche inertielle u* est solution d'une equation
660 // differentielle resolue de maniere iterative
661
662 const double Xpini = 40.;
663 const int itmax = 25;
664 const double seuil = 0.01;
665 const double c1 = Kappa*norm_vit;
666 const double c2 = log(Erugu*dist/d_visco); // log = logarithme neperien
667 double u_star,u_star1;
668
669 u_star = Xpini*d_visco/dist;
670
671 int iter = 0;
672 u_star1 = 1;
673 while ((iter++<itmax) && (std::fabs(u_star1) > seuil))
674 {
675 u_star1 = (c1-u_star*(log(u_star) + c2))/(c1 + u_star);
676 u_star = (1+u_star1)*u_star;
677 }
678 if (std::fabs(u_star1) >= seuil) erreur_non_convergence();
679
680 tab_u_star_(face)= u_star;
681
682 return 1;
683}
684
685void Paroi_2couches_VDF::imprimer_ustar(Sortie& os) const
686{
687 const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_dis_.valeur());
688 int ndeb,nfin;
689 double upmoy,dpmoy,utaumoy;
690 upmoy=0.;
691 dpmoy=0.;
692 utaumoy=0.;
693 int compt=0;
694
695 EcrFicPartage Ustar;
696 ouvrir_fichier_partage(Ustar,"Ustar");
697
698 for (int n_bord=0; n_bord<domaine_VDF.nb_front_Cl(); n_bord++)
699 {
700 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
701 if ( (sub_type(Dirichlet_paroi_fixe,la_cl.valeur())) ||
702 (sub_type(Dirichlet_paroi_defilante,la_cl.valeur()) ))
703 {
704 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
705 if(je_suis_maitre())
706 {
707 Ustar << finl;
708 Ustar << "Bord " << le_bord.le_nom() << finl;
709 if (dimension == 2)
710 {
711 Ustar << "----------------------------------------------------------------------------------------------------------------------------------------------------------------" << finl;
712 Ustar << "\tFace a\t\t\t\t|\t\t\t\t\t\t\t\t\t| TAU=Nu.Grad(Ut) [m2/s2]" << finl;
713 Ustar << "----------------------------------------|-----------------------------------------------------------------------|-----------------------------------------------" << finl;
714 Ustar << "X\t\t| Y\t\t\t| u+\t\t\t| d+\t\t\t| u*\t\t\t| |TAUx|\t\t| |TAUy|" << finl;
715 Ustar << "----------------|-----------------------|-----------------------|-----------------------|-----------------------|-----------------------|-----------------------" << finl;
716 }
717 if (dimension == 3)
718 {
719 Ustar << "---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------" << finl;
720 Ustar << "\tFace a\t\t\t\t\t\t\t|\t\t\t\t\t\t\t\t\t| TAU=Nu.Grad(Ut) [m2/s2]" << finl;
721 Ustar << "----------------------------------------------------------------|-----------------------------------------------------------------------|----------------------------------------------------------------" << finl;
722 Ustar << "X\t\t| Y\t\t\t| Z\t\t\t| u+\t\t\t| d+\t\t\t| u*\t\t\t| |TAUx|\t\t| |TAUy|\t\t| |TAUz|" << finl;
723 Ustar << "----------------|-----------------------|-----------------------|-----------------------|-----------------------|-----------------------|-----------------------|-----------------------|----------------" << finl;
724 }
725 }
726 ndeb = le_bord.num_premiere_face();
727 nfin = ndeb + le_bord.nb_faces();
728 for (int num_face=ndeb; num_face<nfin; num_face++)
729 {
730 double x=domaine_VDF.xv(num_face,0);
731 double y=domaine_VDF.xv(num_face,1);
732 if (dimension == 2)
733 Ustar << x << "\t| " << y;
734 if (dimension == 3)
735 {
736 double z=domaine_VDF.xv(num_face,2);
737 Ustar << x << "\t| " << y << "\t| " << z;
738 }
739 Ustar << "\t| " << uplus_(num_face) << "\t| " << tab_d_plus(num_face) << "\t| " << tab_u_star(num_face);
740 Ustar << "\t| " << Cisaillement_paroi_(num_face,1) << "\t| " << Cisaillement_paroi_(num_face,0) ;
741 if (dimension == 3)
742 Ustar << "\t| " << Cisaillement_paroi_(num_face,2) << finl;
743 else
744 Ustar << finl;
745
746 upmoy +=uplus_(num_face);
747 dpmoy +=tab_d_plus(num_face);
748 utaumoy +=tab_u_star(num_face);
749 compt +=1;
750 }
751 Ustar.syncfile();
752 }
753 }
754 upmoy = mp_sum(upmoy);
755 dpmoy = mp_sum(dpmoy);
756 utaumoy = mp_sum(utaumoy);
757 compt = check_int_overflow(mp_sum(compt));
758 if (compt && je_suis_maitre())
759 {
760 Ustar << finl;
761 Ustar << "-------------------------------------------------------------" << finl;
762 Ustar << "Calcul des valeurs moyennes (en supposant maillage regulier):" << finl;
763 Ustar << "<u+>= " << upmoy/compt << " <d+>= " << dpmoy/compt << " <u*>= " << utaumoy/compt << finl;
764 }
765 if(je_suis_maitre())
766 Ustar << finl << finl;
767 Ustar.syncfile();
768}
769
770
771
Champ_Don_base
classe Champ_Don_base classe de base des Champs donnes (non calcules)
Definition Champ_Don_base.h:32
Champ_Don_base::valeurs
DoubleTab & valeurs() override
Surcharge Champ_base::valeurs() Renvoie le tableau des valeurs.
Definition Champ_Don_base.h:49
Champ_Inc_base::valeurs
DoubleTab & valeurs() override
Renvoie le tableau des valeurs du champ au temps courant.
Definition Champ_Inc_base.h:77
Champ_Uniforme
classe Champ_Uniforme Represente un champ constant dans l'espace et dans le temps.
Definition Champ_Uniforme.h:26
Cond_lim
classe Cond_lim Classe generique servant a representer n'importe quelle classe
Definition Cond_lim.h:31
Debog::verifier
static void verifier(const char *const msg, double)
Definition Debog.cpp:21
Dirichlet_paroi_defilante
classe Dirichlet_paroi_defilante Impose la vitesse de paroi dnas une equation de type Navier_Stokes.
Definition Dirichlet_paroi_defilante.h:26
Dirichlet_paroi_fixe
classe Dirichlet_paroi_fixe Represente une paroi immobile dans une equation de type Navier_Stokes.
Definition Dirichlet_paroi_fixe.h:26
Dirichlet::val_imp
virtual double val_imp(int i) const
Renvoie la valeur imposee sur la i-eme composante du champ a la frontiere au temps par defaut du cham...
Definition Dirichlet.cpp:35
Domaine_VDF
class Domaine_VDF
Definition Domaine_VDF.h:64
Domaine_VDF::orientation
int orientation(int) const override
inline DoubleVect& Domaine_VDF::porosite_face() {
Definition Domaine_VDF.h:297
Domaine_VDF::distance_normale
double distance_normale(int num_face) const
Definition Domaine_VDF.h:457
Domaine_VF::xv
double xv(int num_face, int k) const
Definition Domaine_VF.h:76
Domaine_VF::elem_faces
int elem_faces(int i, int j) const
renvoie le numero de le ieme face de la maille num_elem la facon dont ces faces sont numerotees est
Definition Domaine_VF.h:543
Domaine_VF::face_voisins
int face_voisins(int num_face, int i) const
renvoie l'element voisin de numface dans la direction i.
Definition Domaine_VF.h:418
Domaine_dis_base::nb_front_Cl
int nb_front_Cl() const
Definition Domaine_dis_base.h:53
EcrFicPartage::syncfile
Sortie & syncfile() override
Provoque l'ecriture sur disque des donnees accumulees sur les differents processeurs depuis le dernie...
Definition EcrFicPartage.cpp:173
Entree
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
Equation_base
classe Equation_base Le role d'une equation est le calcul d'un ou plusieurs champs....
Definition Equation_base.h:76
Equation_base::milieu
virtual const Milieu_base & milieu() const =0
Equation_base::inconnue
virtual const Champ_Inc_base & inconnue() const =0
Fluide_base
classe Fluide_base Cette classe represente un d'un fluide incompressible ainsi que
Definition Fluide_base.h:38
Fluide_base::viscosite_cinematique
const Champ_Don_base & viscosite_cinematique() const
Definition Fluide_base.h:58
Front_VF
class Front_VF
Definition Front_VF.h:36
Front_VF::nb_faces
int nb_faces() const
Definition Front_VF.h:53
Front_VF::num_premiere_face
int num_premiere_face() const
Definition Front_VF.h:63
Frontiere_dis_base::le_nom
const Nom & le_nom() const override
Renvoie le nom de la frontiere geometrique.
Definition Frontiere_dis_base.cpp:81
Loi_2couches_base
classe Loi_2couches_base Cette classe de base represente les modeles 1 equation pour le modele a deux...
Definition Loi_2couches_base.h:26
Loi_2couches_base::LepsLmu
virtual void LepsLmu(double k, double nu, double dist, double y_etoile, double &Leps, double &Lmu, double &vvSqRt)=0
Modele_turbulence_hyd_K_Eps_2_Couches
Classe Modele_turbulence_hyd_K_Eps_2_Couches Cette classe represente le modele de turbulence (k,...
Definition Modele_turbulence_hyd_K_Eps_2_Couches.h:28
Modele_turbulence_hyd_K_Eps_2_Couches::get_nutswitch
int get_nutswitch() const
Renvoie la valeur de nut/nu qui delimite les deux couches.
Definition Modele_turbulence_hyd_K_Eps_2_Couches.h:100
Modele_turbulence_hyd_K_Eps_2_Couches::get_set_eq_transport
Transport_K_KEps & get_set_eq_transport() override
Definition Modele_turbulence_hyd_K_Eps_2_Couches.h:51
Modele_turbulence_hyd_K_Eps_2_Couches::get_nbcouches
int get_nbcouches() const
Renvoie le nombre couches utilises pour les lois de paroi Ce nombre est porte par l'equation de trans...
Definition Modele_turbulence_hyd_K_Eps_2_Couches.h:67
Modele_turbulence_hyd_K_Eps_2_Couches::get_switch
int get_switch() const
Renvoie 0 si on choisit le switch par y*, 1 par nu_t.
Definition Modele_turbulence_hyd_K_Eps_2_Couches.h:89
Modele_turbulence_hyd_K_Eps_2_Couches::get_yswitch
int get_yswitch() const
Renvoie le y* de switch entre les deux couches pour le modele a deux couches.
Definition Modele_turbulence_hyd_K_Eps_2_Couches.h:78
Objet_U::set_param
virtual void set_param(Param &) const
Definition Objet_U.h:135
Objet_U::dimension
static int dimension
Definition Objet_U.h:99
Objet_U::Sortie
friend class Sortie
Definition Objet_U.h:75
Objet_U::que_suis_je
const Nom & que_suis_je() const
renvoie la chaine identifiant la classe.
Definition Objet_U.cpp:104
Objet_U::readOn
virtual Entree & readOn(Entree &)
Lecture d'un Objet_U sur un flot d'entree Methode a surcharger.
Definition Objet_U.cpp:293
Objet_U::le_nom
virtual const Nom & le_nom() const
Donne le nom de l'Objet_U Methode a surcharger : renvoie "neant" dans cette implementation.
Definition Objet_U.cpp:319
Objet_U::printOn
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
Definition Objet_U.cpp:282
Param
Helper class to factorize the readOn method of Objet_U classes.
Definition Param.h:112
Paroi_2couches_VDF
CLASS: Paroi_2couches_VDF.
Definition Paroi_2couches_VDF.h:35
Paroi_2couches_VDF::init_lois_paroi
int init_lois_paroi() override
Definition Paroi_2couches_VDF.cpp:60
Paroi_2couches_VDF::calculer_hyd
int calculer_hyd(DoubleTab &) override
Definition Paroi_2couches_VDF.cpp:77
Paroi_2couches_VDF::uplus_
DoubleVect uplus_
Definition Paroi_2couches_VDF.h:65
Paroi_2couches_VDF::imprimer_ustar
void imprimer_ustar(Sortie &os) const override
Definition Paroi_2couches_VDF.cpp:685
Paroi_2couches_VDF::set_param
void set_param(Param &param) const override
Definition Paroi_2couches_VDF.cpp:47
Paroi_2couches_VDF::calculer_uplus_dplus
void calculer_uplus_dplus(DoubleVect &, DoubleVect &, DoubleVect &, int, double, double, double)
Definition Paroi_2couches_VDF.cpp:603
Paroi_2couches_VDF::calculer_hyd_BiK
int calculer_hyd_BiK(DoubleTab &, DoubleTab &) override
Definition Paroi_2couches_VDF.cpp:83
Paroi_2couches_VDF::Erugu
double Erugu
Definition Paroi_2couches_VDF.h:61
Paroi_2couches_VDF::calculer_u_star_sous_couche_log
int calculer_u_star_sous_couche_log(double, double, double, int)
Definition Paroi_2couches_VDF.cpp:656
Paroi_2couches_VDF::init_lois_paroi_hydraulique
int init_lois_paroi_hydraulique()
Definition Paroi_2couches_VDF.cpp:70
Paroi_2couches_VDF::Kappa
double Kappa
Definition Paroi_2couches_VDF.h:60
Paroi_2couches_VDF::calculer_u_star_sous_couche_tampon
int calculer_u_star_sous_couche_tampon(double &, double, double, double, int)
Definition Paroi_2couches_VDF.cpp:613
Paroi_hyd_base_VDF
Definition Paroi_hyd_base_VDF.h:29
Paroi_hyd_base_VDF::init_lois_paroi_
void init_lois_paroi_()
Definition Paroi_hyd_base_VDF.cpp:37
Paroi_log_QDM::init_lois_paroi_hydraulique_
void init_lois_paroi_hydraulique_()
Definition Paroi_log_QDM.cpp:32
Paroi_log_QDM::Cmu_
double Cmu_
Definition Paroi_log_QDM.h:52
Paroi_log_QDM::set_param
void set_param(Param &param) const
Definition Paroi_log_QDM.cpp:22
Paroi_log_QDM::table_hyd
Table table_hyd
Definition Paroi_log_QDM.h:56
Process::check_int_overflow
static int check_int_overflow(trustIdType)
Definition Process.cpp:428
Process::mp_sum
static double mp_sum(double)
Calcule la somme de x sur tous les processeurs du groupe courant.
Definition Process.cpp:146
Process::je_suis_maitre
static int je_suis_maitre()
renvoie 1 si on est sur le processeur maitre du groupe courant (c'est a dire me() == 0),...
Definition Process.cpp:86
Sortie
Classe de base des flux de sortie.
Definition Sortie.h:52
TRUSTVect::local_min_vect
_TYPE_ local_min_vect(Mp_vect_options opt=VECT_REAL_ITEMS) const
Definition TRUSTVect.h:155
TRUSTVect::echange_espace_virtuel
virtual void echange_espace_virtuel(IsExchangeBlocking exchange_type=IsExchangeBlocking::DefaultBlocking, const std::string kernel_name="noname")
Definition TRUSTVect.tpp:282
Transport_K_KEps
classe Transport_K_KEps Cette classe represente l'equation de transport de l'energie cinetique
Definition Transport_K_KEps.h:40
Turbulence_paroi_base::tab_d_plus_
DoubleVect tab_d_plus_
Definition Turbulence_paroi_base.h:88
Turbulence_paroi_base::tab_d_plus
const DoubleVect & tab_d_plus() const
Definition Turbulence_paroi_base.h:141
Turbulence_paroi_base::tab_u_star
const DoubleVect & tab_u_star() const
Definition Turbulence_paroi_base.h:131
Turbulence_paroi_base::tab_u_star_
DoubleVect tab_u_star_
Definition Turbulence_paroi_base.h:87
Turbulence_paroi_base::ouvrir_fichier_partage
void ouvrir_fichier_partage(EcrFicPartage &, const Nom &) const
Ouverture/creation d'un fichier d'impression de Face, uplus_, dplus_, tab_u_star, Cisaillement_paroi_...
Definition Turbulence_paroi_base.cpp:172
Turbulence_paroi_base::Cisaillement_paroi_
DoubleTab Cisaillement_paroi_
Definition Turbulence_paroi_base.h:85