TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
ParoiVDF_TBLE_LRM.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
17#include <ParoiVDF_TBLE_LRM.h>
18#include <Motcle.h>
19#include <Front_VF.h>
20#include <Domaine_Cl_VDF.h>
21#include <Dirichlet_paroi_fixe.h>
22#include <Dirichlet_paroi_defilante.h>
23
24#include <Champ_Uniforme.h>
25#include <Fluide_base.h>
26#include <EFichier.h>
27#include <SFichier.h>
28#include <Diffu_k.h>
29#include <Paroi_TBLE_scal_VDF.h>
30#include <Modele_turbulence_scal_base.h>
31#include <Probleme_base.h>
32#include <Terme_Boussinesq_VDF_Face.h>
33#include <Navier_Stokes_std.h>
34#include <Modele_turbulence_hyd_base.h>
35
36Implemente_instanciable_sans_constructeur(ParoiVDF_TBLE_LRM,"Paroi_TBLE_LRM_VDF",Paroi_hyd_base_VDF);
37
38// printOn()
39/////
40
41Sortie& ParoiVDF_TBLE_LRM::printOn(Sortie& os) const
42{
43 return os << que_suis_je() << " " << le_nom();
44}
45
46//// readOn
47//
48
49Entree& ParoiVDF_TBLE_LRM::readOn(Entree& is)
50{
51
52 //Cerr << "debut readon";
53
54 Motcle mot_lu;
55
56 // Valeurs par defaut
57 nb_pts=-1; // nb de pts dans le maillage fin
58 modele_visco = "Diffu_k"; //modele de viscosite turbulente
59 nb_comp = 2; //nb de composantes dans le maillage fin
60 fac=1.; // facteur de raffinement du maillage fin
61 epsilon = 1.e-5; // seuil de resolution dans le maillage fin
62 max_it = 1000; // max d'iterations dans le maillage fin
63 // k_init=-1; //condition initiale pour k
64 avec_boussi=1;//par defaut on met les termes de Boussinesq
65 // FIN valeurs par defaut
66 nb_post_pts=0;
67
68 Motcles les_mots(13);
69 {
70 les_mots[0]="N";
71 les_mots[1]="modele_visco";
72 les_mots[2]="N_comp";
73 les_mots[3]="facteur";
74 // mot cle 4 libre
75 les_mots[5]="epsilon";
76 les_mots[6]="max_it";
77 les_mots[9]="type_regime";
78 les_mots[10]="sans_Boussinesq";
79 les_mots[11]="k_init";
80 les_mots[12]="sonde_tble";
81
82 }
83 Motcle acc_ouverte("{");
84 Motcle acc_fermee("}");
85 Motcle type_reg="vide";
86 Motcle nom_classe_mod = "Mod_echelle_";
87
88
89 is >> mot_lu;
90 if(mot_lu == acc_ouverte)
91 {
92 // on passe l'accolade ouvrante
93 is >> mot_lu;
94 }
95 while(mot_lu != acc_fermee)
96 {
97 int rang=les_mots.search(mot_lu);
98 switch(rang)
99 {
100 case 0 :
101 is >> nb_pts;
102 Cerr << "N =" << nb_pts << finl;
103 break;
104 case 1 :
105 is >> modele_visco;
106 Cerr << "Le modele de viscosite turbulente dans le maillge fin est "
107 << modele_visco << finl;
108 break;
109 case 2 :
110 is >> nb_comp;
111 Cerr << "Nombre de composantes de vitesse tangentielles a la paroi dans le maillage fin: "
112 << nb_comp << finl;
113 break;
114
115 case 3 :
116 is >> fac;
117 Cerr << "Raison geometrique du maillage fin : " << fac << finl;
118 break;
119 case 5 :
120 is >> epsilon;
121 Cerr << "Seuil de convergence pour les equations de couche limites : " << epsilon << finl;
122 break;
123 case 6 :
124 is >> max_it;
125 Cerr << "Le maximum d'iterations pour la resolution des equations TBLE est de : " << max_it << finl;
126 break;
127 case 9 :
128 is >> type_reg;
129 nom_classe_mod+=type_reg;
130 mod_ech.typer(nom_classe_mod);
131 break;
132 case 10 :
133 avec_boussi=0;
134 break;
135 case 11 :
136 is >> k_init;
137 Cerr << "k_init = " << k_init << finl;
138 break;
139 case 12 :
140 is >> nb_post_pts;
141 sonde_tble.resize(nb_post_pts, dimension);
142 for(int i=0; i<nb_post_pts; i++)
143 {
144 nom_pts.dimensionner(nb_post_pts);
145 is >> nom_pts[i];
146 Cerr << "Nom"<<i<<"=" << nom_pts[i] <<finl;
147 for(int j=0; j<dimension ; j++)
148 {
149 is >> sonde_tble(i,j);
150 Cerr << "sonde_tble( "<< i << "," << j <<" ) =" << sonde_tble(i,j) << finl;
151 }
152 }
153 break;
154
155 default :
156 {
157 Cerr << mot_lu << " n'est pas un mot compris par ParoiVDF_TBLE_LRM" << finl;
158 Cerr << "Les mots compris sont : " << les_mots << finl;
159 Process::exit();
160 }
161 }
162 is >> mot_lu;
163 }
164
165 if (type_reg=="vide")
166 {
167 Cerr << "Le modele d'echelles l_eps, l_mu n'est pas precise !" << finl;
168 Process::exit();
169 }
170
171
172 return is;
173
174 //Cerr << "fin readon";
175
176}
177
178/////////////////////////////////////////////////////////////////////
179//
180// Implementation des fonctions de la classe ParoiVDF_TBLE_LRM
181//
182/////////////////////////////////////////////////////////////////////
183
184
185// /////////////////////////////////////////////////
186// // Initialisation des tableaux
187// ////////////////////////////////////////////////
188
189int ParoiVDF_TBLE_LRM::init_lois_paroi()
190{
191 Cerr << "debut init_lois_paroi" << finl;
192
193 const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_dis_.valeur());
194 const IntVect& orientation = domaine_VDF.orientation();
195 const IntTab& face_voisins = domaine_VDF.face_voisins();
196 const IntTab& elem_faces = domaine_VDF.elem_faces();
197 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
198 const Fluide_base& le_fluide = ref_cast(Fluide_base, eqn_hydr.milieu());
199 const DoubleVect& vit = eqn_hydr.inconnue().valeurs();
200 // eq_k_U_W.dimensionner(le_dom_dis_->nb_faces_bord());
201 int compteur_faces_paroi = 0;
202 DoubleVect corresp(le_dom_dis_->nb_faces_bord());
203
204 // B.M. Nouveau codage a tester: je n'ai pas trouve ou est fait le resize()
205 // initial du tableau a nb_faces_bord. Je ne sais pas si le tableau contient
206 // deja des valeurs ou pas, s'il faut les conserver ou pas !
207 Process::exit();
208 Cisaillement_paroi_.resize(0, dimension);
209 domaine_VDF.creer_tableau_faces_bord(Cisaillement_paroi_);
210
211 int ndeb=0,nfin=0;
212 int elem, ori;
213 double vmoy = 0.;
214 double dist; //distance du premier centre de maille a la paroi
215
216
217 // SFichier fic_corresp("corresp.dat",ios::app); // impression de la correspondance
218
219 // Boucle sur les bords
220 //Cerr << "Boucle sur les bords" << finl;
221
222
223 for (int n_bord=0; n_bord<domaine_VDF.nb_front_Cl(); n_bord++)
224 {
225 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
226
227 if (sub_type(Dirichlet_paroi_fixe,la_cl.valeur()) )
228
229 {
230 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
231 ndeb = le_bord.num_premiere_face();
232 nfin = ndeb + le_bord.nb_faces();
233
234 for (int num_face=ndeb; num_face<nfin; num_face++)
235 {
236 corresp(compteur_faces_paroi)=num_face;
237 // fic_corresp << compteur_faces_paroi << " " << num_face << finl;
238 // Cerr << compteur_faces_paroi << " " << num_face << finl;
239 compteur_faces_paroi++;
240
241
242 }
243 }
244 }
245
246
247
248 compteur_faces_paroi = le_dom_dis_->nb_faces_bord();
249 // Cerr << "compteur_faces_paroi pour dimensionnement = " << compteur_faces_paroi << finl;
250 eq_k_U_W.dimensionner(compteur_faces_paroi); //Dimensionnement du vecteur eq_couch_lim
251 corresp.resize(compteur_faces_paroi); //Redimensionnement de corresp
252
253
254 num_faces_post.resize(nb_post_pts);
255 num_faces_post2.resize(nb_post_pts);
256 for(int j=0; j<nb_post_pts; j++)
257 {
258 double d0=1.e9;
259 int face=-1;
260 int face2=-1;
261 compteur_faces_paroi=0;
262 for (int n_bord=0; n_bord<domaine_VDF.nb_front_Cl(); n_bord++)
263 {
264
265 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
266
267 if (sub_type(Dirichlet_paroi_fixe,la_cl.valeur()) )
268
269 {
270 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
271 ndeb = le_bord.num_premiere_face();
272 nfin = ndeb + le_bord.nb_faces();
273
274 for (int num_face=ndeb; num_face<nfin; num_face++)
275 {
276 double d1=0.;
277 for (int d=0; d<dimension; d++)
278 {
279 double x=domaine_VDF.xv(num_face,d)-sonde_tble(j,d);
280 d1+=x*x;
281 }
282
283 if (d1<d0)
284 {
285 d0=d1;
286 face=compteur_faces_paroi;
287 face2=num_face;
288 }
289 compteur_faces_paroi++;
290 }
291 }
292 }
293 //face en local
294 num_faces_post(j)=face;
295 //face en global
296 num_faces_post2(j)=face2;
297 }
298
299
300 compteur_faces_paroi = 0; //Reinitialisation de compteur_faces_paroi
301
302
303 for (int n_bord=0; n_bord<domaine_VDF.nb_front_Cl(); n_bord++)
304 {
305
306 // pour chaque condition limite on regarde son type
307 // On applique les lois de paroi uniquement
308 // aux voisinages des parois
309
310 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
311
312 if (sub_type(Dirichlet_paroi_fixe,la_cl.valeur()) )
313
314 {
315 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
316 ndeb = le_bord.num_premiere_face();
317 nfin = ndeb + le_bord.nb_faces();
318
319
320
321 ////////////
322 /// 2D ////
323 ////////////
324
325 if (dimension == 2 )
326 {
327 //Boucle sur les faces des bords parietaux
328
329 for (int num_face=ndeb; num_face<nfin; num_face++)
330 {
331
332 //ici ori=direction perpendiculaire a la paroi
333 ori = orientation(num_face);
334
335 if ((elem = face_voisins(num_face,0)) == -1)
336 {
337 elem = face_voisins(num_face,1);
338 }
339
340
341 eq_k_U_W[compteur_faces_paroi].associer_milieu(le_fluide);
342
343 //Distance a la paroi du 1er centre de maille
344 if (axi)
345 dist = domaine_VDF.dist_norm_bord_axi(num_face);
346 else
347 dist = domaine_VDF.dist_norm_bord(num_face);
348
349 eq_k_U_W[compteur_faces_paroi].set_y0(0.); //ordonnee de la paroi
350 eq_k_U_W[compteur_faces_paroi].set_yn(dist); //ordonnee du 1er centre de maille
351
352
353 eq_k_U_W[compteur_faces_paroi].initialiser(nb_comp+1, nb_pts,
354 fac, epsilon, max_it,0); //nbre de pts maillage fin
355
356 eq_k_U_W[compteur_faces_paroi].set_u_y0(0,0.);//1ere composante de la vitesse a la paroi
357 eq_k_U_W[compteur_faces_paroi].set_u_y0(1,0.); //Composante de k a la paroi
358
359
360
361 eq_k_U_W[compteur_faces_paroi].set_diffu(modele_visco); //modele de viscosite turbulente
362
363
364 if (sub_type(Diffu_k,eq_k_U_W[compteur_faces_paroi].get_diffu()))
365 {
366 Diffu_k& diffu = ref_cast_non_const(Diffu_k, eq_k_U_W[compteur_faces_paroi].get_diffu());
367 diffu.associer_modele_echelles(mod_ech.valeur());
368 }
369
370 vmoy = 0.5*(vit(elem_faces(elem,(ori+1)%4)) + vit(elem_faces(elem,(ori+3)%4)));
371
372 //vitesse sur le maillage fin a l'instant initial
373 eq_k_U_W[compteur_faces_paroi].set_u_yn(0,vmoy); //1ere composante de la vitesse en yn
374 eq_k_U_W[compteur_faces_paroi].set_Uinit_lin(0, 0., vmoy);
375
376 vmoy = 0.5*(vit(elem_faces(elem,(ori+1)%4)) + vit(elem_faces(elem,(ori+3)%4)));
377
378 //vmoy /= 2.;
379
380
381 eq_k_U_W[compteur_faces_paroi].set_Uinit_lin(1, 0., k_init); // a ameliorer pour la condition initiale de k
382
383
384 compteur_faces_paroi++;
385
386 }//Fin boucle sur les faces de bord
387
388 }
389
390
391
392 ////////////
393 /// 3D ////
394 ////////////
395
396 if (dimension == 3)
397 {
398
399 //Boucle sur les faces des bords parietaux
400
401 for (int num_face=ndeb; num_face<nfin; num_face++)
402 {
403
404 //ici ori=direction perpendiculaire a la paroi
405 ori = orientation(num_face);
406
407 if ((elem = face_voisins(num_face,0)) == -1)
408 {
409 elem = face_voisins(num_face,1);
410 }
411
412 eq_k_U_W[compteur_faces_paroi].associer_milieu(le_fluide);
413
414 //Distance a la paroi du 1er centre de maille
415 if (axi)
416 dist = domaine_VDF.dist_norm_bord_axi(num_face);
417 else
418 dist = domaine_VDF.dist_norm_bord(num_face);
419
420 eq_k_U_W[compteur_faces_paroi].set_y0(0.); //ordonnee de la paroi
421 eq_k_U_W[compteur_faces_paroi].set_yn(dist); //ordonnee du 1er centre de maille
422
423
424 eq_k_U_W[compteur_faces_paroi].initialiser(nb_comp+1, nb_pts,
425 fac, epsilon, max_it,0);
426
427
428 eq_k_U_W[compteur_faces_paroi].set_u_y0(0,0.); //1ere composante de la vitesse a la paroi
429 eq_k_U_W[compteur_faces_paroi].set_u_y0(1,0.); //2e composante de la vitesse a la paroi
430 eq_k_U_W[compteur_faces_paroi].set_u_y0(2,0.); //k a la paroi
431
432 eq_k_U_W[compteur_faces_paroi].set_diffu(modele_visco); //modele de viscosite turbulente
433
434 if (sub_type(Diffu_k,eq_k_U_W[compteur_faces_paroi].get_diffu()))
435 {
436 Diffu_k& diffu = ref_cast_non_const(Diffu_k, eq_k_U_W[compteur_faces_paroi].get_diffu());
437 diffu.associer_modele_echelles(mod_ech.valeur());
438 }
439
440 ///a voir :!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
441
442 vmoy = 0.5*(vit(elem_faces(elem,(ori+1)%6)) + vit(elem_faces(elem,(ori+4)%6)));
443
444 eq_k_U_W[compteur_faces_paroi].set_u_yn(0,vmoy); //1ere composante de la vitesse en yn
445
446
447 //vitesse sur le maillage fin a l'instant initial
448 eq_k_U_W[compteur_faces_paroi].set_Uinit_lin(0, 0., vmoy);
449
450 vmoy = 0.5*(vit(elem_faces(elem,(ori+2)%6)) + vit(elem_faces(elem,(ori+5)%6)));
451
452 //vmoy /= 2.;
453
454 // eq_k_U_W[compteur_faces_paroi].set_u_yn(1,1.); //2e composante de la vitesse en yn
455 eq_k_U_W[compteur_faces_paroi].set_u_yn(1,vmoy); //2e composante de la vitesse en yn
456 eq_k_U_W[compteur_faces_paroi].set_Uinit_lin(1, 0., vmoy);
457
458
459 //pour initialiser le calcul, on prend le u_tau theorique
460
461 eq_k_U_W[compteur_faces_paroi].set_Uinit_lin(2, 0., k_init); // a ameliorer pour la condition initiale de k
462
463 compteur_faces_paroi++;
464
465 }//Fin boucle sur les faces de bord
466 }// Fin if dim 3
467 }//Fin if Paroi_fixe
468 }//Fin boucle sur les bords parietaux
469
470 Cerr << "fin init_lois_paroi" << finl;
471
472 return 1;
473}
474
475
476/////////////////////////////////////////////////
477//Calcul du cisaillement a la paroi
478////////////////////////////////////////////////
479
480int ParoiVDF_TBLE_LRM::calculer_hyd(DoubleTab& nu_t, DoubleTab& tab_k)
481{
482 Cerr <<"TBLE LRM non code en LES" << finl;
483 Process::exit();
484 return 0;
485}
486
487
488int ParoiVDF_TBLE_LRM::calculer_hyd_BiK(DoubleTab& tab_k, DoubleTab& tab_eps)
489{
490
491 if (Process::nproc()>1)
492 {
493 Cerr << "ParoiVDF_TBLE_LRM::calculer_hyd n'est pas parallelise." << finl;
494 Process::exit();
495 }
496 const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_dis_.valeur());
497 const IntVect& orientation = domaine_VDF.orientation();
498 const IntTab& face_voisins = domaine_VDF.face_voisins();
499 const IntTab& elem_faces = domaine_VDF.elem_faces();
500 const DoubleVect& volumes_entrelaces = domaine_VDF.volumes_entrelaces();
501 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
502 const Fluide_base& le_fluide = ref_cast(Fluide_base,eqn_hydr.milieu());
503 const Champ_Don_base& ch_visco_cin = le_fluide.viscosite_cinematique();
504 const DoubleTab& tab_visco = ch_visco_cin.valeurs();
505 Mod_echelle_LRM_base& le_mod_ech= mod_ech.valeur();
506
507 const DoubleTab& nu_t = mon_modele_turb_hyd->viscosite_turbulente().valeurs();
508
509
510
511 double visco=-1;
512 int l_unif;
513 if (sub_type(Champ_Uniforme,ch_visco_cin))
514 {
515 visco = std::max(tab_visco(0,0),DMINFLOAT);
516 l_unif = 1;
517 }
518 else
519 l_unif = 0;
520
521 // preparer_calcul_hyd(tab);
522 if ((!l_unif) && (tab_visco.local_min_vect()<DMINFLOAT))
523 // on ne doit pas changer tab_visco ici !
524 {
525 Cerr << "In ParoiVDF_TBLE_LRM::calculer_hyd_BiK : visco = " << tab_visco.local_min_vect() << " <= 0 ? " << finl;
526 throw;
527 }
528 // tab_visco+=DMINFLOAT;
529
530
531 int ndeb=0,nfin=0;
532 int ori;
533 // int elem_gauche0, elem_gauche1, elem_droite0, elem_droite1;
534 // int face_gauche0, face_gauche1, face_droite0, face_droite1;
535 int elem;
536 // int compteur=0;
537
538 int signe;
539 int itmax=0;
540
541 double vmoy = 0., ts0 =0., ts1=0.,gradient_de_pression0=0., gradient_de_pression1 = 0.;
542 const DoubleTab& vit = eqn_hydr.inconnue().valeurs(); //vitesse
543 DoubleVect grad_vit_elemx;
544 DoubleVect grad_vit_elem_moyx;
545 DoubleVect grad_vit_elemz;
546 DoubleVect grad_vit_elem_moyz;
547 DoubleVect grad_T;
548 DoubleVect grad_T_moy;
549 DoubleVect source_T_k;
550 DoubleVect source_T_U;
551
552
553 double eps0;
554 double eps;
555 double l_eps;
556 double vv_bar;
557 double u=0;
558
559 grad_vit_elemx.resize(nb_pts);
560 grad_vit_elem_moyx.resize(nb_pts);
561 grad_vit_elemz.resize(nb_pts);
562 grad_vit_elem_moyz.resize(nb_pts);
563 grad_T.resize(nb_pts);
564 grad_T_moy.resize(nb_pts);
565 source_T_k.resize(nb_pts);
566 source_T_U.resize(nb_pts);
567
568 const Navier_Stokes_std& eqnNS = ref_cast(Navier_Stokes_std, eqn_hydr);
569
570 DoubleTab grad_p(vit);
571 const DoubleTab& p = eqnNS.pression().valeurs();
572 const Operateur_Grad& gradient = eqnNS.operateur_gradient();
573 gradient.calculer(p, grad_p); // Calcul du gradient de pression
574
575 DoubleTab termes_sources;
576 termes_sources.resize(domaine_VDF.nb_faces(),1);
577 eqn_hydr.sources().calculer(termes_sources); //les termes sources
578
579 const double tps = eqnNS.schema_temps().temps_courant();
580 const double dt = eqnNS.schema_temps().pas_de_temps();
581 const double dt_min = eqnNS.schema_temps().pas_temps_min();
582
583
584 //SFichier fic_v("v.dat",ios::app); // impression de la vitesse normale pour test
585
586 int compteur_faces_paroi = 0;
587 DoubleVect corresp(le_dom_dis_->nb_faces_bord());
588
589
590 // Boucle sur les bords
591 //Cerr << "Boucle sur les bords" << finl;
592
593 for (int n_bord=0; n_bord<domaine_VDF.nb_front_Cl(); n_bord++)
594 {
595 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
596
597 if (sub_type(Dirichlet_paroi_fixe,la_cl.valeur()) )
598
599 {
600 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
601 ndeb = le_bord.num_premiere_face();
602 nfin = ndeb + le_bord.nb_faces();
603
604 for (int num_face=ndeb; num_face<nfin; num_face++)
605 {
606 corresp(compteur_faces_paroi)=num_face;
607 compteur_faces_paroi++;
608 }
609 }
610 }
611
612 corresp.resize(compteur_faces_paroi); //Redimensionnement de corresp
613 compteur_faces_paroi = 0; //Reinitialisation de compteur_faces_paroi
614
615
616 // Boucle sur les bords
617
618 //Cerr << "Boucle sur les bords " << finl;
619 for (int n_bord=0; n_bord<domaine_VDF.nb_front_Cl(); n_bord++)
620 {
621
622 // pour chaque condition limite on regarde son type
623 // On applique les lois de paroi uniquement
624 // aux voisinages des parois
625
626
627 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
628
629 if (sub_type(Dirichlet_paroi_fixe,la_cl.valeur()) )
630
631 {
632 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
633 ndeb = le_bord.num_premiere_face();
634 nfin = ndeb + le_bord.nb_faces();
635
636
637
638 ///////////////////////
639 // Si probleme 2D //
640 ///////////////////////
641
642
643 if(dimension == 2)
644 {
645 //eq_k_U_W.dimensionner(2*le_dom_dis_->nb_faces_bord());
646
647 //Boucle sur les faces des bords parietaux
648
649 for (int num_face=ndeb; num_face<nfin; num_face++)
650 {
651 //ici ori=direction perpendiculaire a la paroi
652 ori = orientation(num_face);
653 //Cerr << "ori = " << ori << finl;
654
655
656
657 //compteur++;
658 //Selection de l'element associe a une face parietale
659 //et initialisation de l'entier signe
660
661 if ((elem = face_voisins(num_face,0)) == -1)
662 {
663 elem = face_voisins(num_face,1);
664 signe = 1;
665 }
666 else signe = -1 ;
667
668
669 if(dt<dt_min)
670 eq_k_U_W[compteur_faces_paroi].set_dt(1.e12);
671 else
672 eq_k_U_W[compteur_faces_paroi].set_dt(dt);
673
674
675
676 // 1er couple de faces perpendiculaires a la paroi
677 int face1 = elem_faces(elem,(ori+1));
678 int face2 = elem_faces(elem,(ori+3)%4);
679
680
681
682 ////////////////////////////////////////////////////////////
683 ////// couple de faces perpendiculaires a la paroi /////
684 ////////////////////////////////////////////////////////////
685
686 vmoy = 0.5*(vit(face1) + vit(face2));
687
688 //vitesse sur le maillage fin a l'instant initial
689 //eq_k_U_W[compteur_faces_paroi].set_Uinit_lin(0, 0., vmoy);
690
691 ts0 = 0.5*(termes_sources(face1)/volumes_entrelaces(face1)
692 + termes_sources(face2)/volumes_entrelaces(face2));
693
694 eq_k_U_W[compteur_faces_paroi].set_u_y0(0,0.); //1ere composante de la vitesse a la paroi
695 eq_k_U_W[compteur_faces_paroi].set_u_yn(0,vmoy); //vitesse en yn
696 // eq_k_U_W[compteur_faces_paroi].set_u_yn(0,1.); //vitesse en yn
697
698
699
700 //*** Gradient de pression du 1er couple de faces
701 // perpendiculaires a la paroi ***
702
703 //gradient de pression (constant dans le maillage fin)
704 gradient_de_pression0 = 0.5*(grad_p(face1)+grad_p(face2));
705
706 //Calcul du gradient On prend U au temps n
707 for(int i=0 ; i<nb_pts ; i++)
708 {
709 grad_vit_elemx(i)=(eq_k_U_W[compteur_faces_paroi].get_Unp1(0,i+1)-eq_k_U_W[compteur_faces_paroi].get_Unp1(0,i))/(eq_k_U_W[compteur_faces_paroi].get_y(i+1)-eq_k_U_W[compteur_faces_paroi].get_y(i));
710 }
711
712 //Cacul du gradient a l'ordonnee y(i)
713 for(int i=1 ; i<nb_pts ; i++)
714 {
715 grad_vit_elem_moyx(i)=0.5*(grad_vit_elemx(i)+grad_vit_elemx(i-1));
716 }
717
718 //gradient en 0
719 grad_vit_elem_moyx(0)=(eq_k_U_W[compteur_faces_paroi].get_Unp1(0,1)-eq_k_U_W[compteur_faces_paroi].get_Unp1(0,0))/(eq_k_U_W[compteur_faces_paroi].get_y(1)-eq_k_U_W[compteur_faces_paroi].get_y(0));
720
721
722
723 double d_visco;
724
725 if (l_unif==1)
726 d_visco = visco;
727 else
728 d_visco = tab_visco[elem];
729
730 // On calcul vv_bar, l_eps, y_nu, y_star avec k au temps n
731 // a la paroi eps est proportionelle a 1
732
733 Diffu_totale_base& diffu = eq_k_U_W[compteur_faces_paroi].get_diffu();
734
735 double dist;
736 if (axi)
737 dist = domaine_VDF.dist_norm_bord_axi(num_face);
738 else
739 dist = domaine_VDF.dist_norm_bord(num_face);
740
741
742
743 // Definition de la face servant a recuperer le k de la deuxieme maille
744 int face3 = elem_faces(elem,(ori+2)%4);
745 int face4 = elem_faces(elem,ori);
746
747 int n2=face3;
748 if (face3==num_face)
749 {
750 n2=face4;
751 }
752
753 int elem1=face_voisins(n2,0);
754 if (elem1 == elem)
755 elem1=face_voisins(n2,1);
756
757
758
759 double dy = domaine_VDF.dist_norm(n2);
760
761
762 // recuperartion de l'element correspondant a la deuxieme maille
763
764 double kcl;
765
766 double a1=d_visco+0.5*(nu_t(elem)+nu_t(elem1));
767 double a2=0.5*(diffu.calculer_a_local(nb_pts)+diffu.calculer_a_local(nb_pts-1));
768
769
770 double b1=a1/dy;
771 double b2=a2/(eq_k_U_W[compteur_faces_paroi].get_y(nb_pts)-eq_k_U_W[compteur_faces_paroi].get_y(nb_pts-1));
772
773
774 kcl=(b1*tab_k(elem1)+b2*eq_k_U_W[compteur_faces_paroi].get_Unp1(1,nb_pts-1))/(b1+b2);
775
776
777 //equation en k ;
778
779 eq_k_U_W[compteur_faces_paroi].set_u_y0(1,0.); //k a la paroi
780 eq_k_U_W[compteur_faces_paroi].set_u_yn(1,kcl);//k a yn
781
782
783 // double Re, Fmu;
784
785 //initialisation de epsilon pour maillage 2D tel que nut(tble)=nut(2D)
786
787
788 if(kcl<0)
789 {
790 kcl=1.e-8;
791 tab_k(elem) =kcl;
792 tab_eps(elem)=0.09*tab_k(elem)*tab_k(elem)/(diffu.calculer_a_local(nb_pts)-d_visco);
793 }
794
795 else
796 {
797 tab_k(elem)=kcl;
798 tab_eps(elem)= 0.09*tab_k(elem)*tab_eps(elem)/(diffu.calculer_a_local(nb_pts)-d_visco);
799 }
800
801
802 double T0=0;
803 double beta_t=0.;
804
805 source_T_k=0.;
806 source_T_U=0.;
807
808 if (avec_boussi==1)
809 {
810 int boussi_ok=0;
811 const Sources& les_sources=eqnNS.sources();
812 int nb_sources=les_sources.size();
813
814 for (int j=0; j<nb_sources; j++)
815 {
816 const Source_base& ts = les_sources(j).valeur();
817 if (sub_type(Terme_Boussinesq_base,ts))
818 {
819 const Terme_Boussinesq_base& terme_boussi = ref_cast(Terme_Boussinesq_base,ts);
820 T0 = terme_boussi.Scalaire0(0);
821 boussi_ok=1;
822 }
823 }
824
825 if (boussi_ok==1)
826 {
827 const Equation_base& eqn_th = mon_modele_turb_hyd->equation().probleme().equation(1);
828 const Modele_turbulence_scal_base& le_mod_turb_th = ref_cast(Modele_turbulence_scal_base,eqn_th.get_modele(TURBULENCE).valeur());
829 const Turbulence_paroi_scal_base& loi = le_mod_turb_th.loi_paroi();
830 Paroi_TBLE_scal_VDF& loi_tble_T = ref_cast_non_const(Paroi_TBLE_scal_VDF,loi);
831
832 const Champ_Don_base& ch_beta_t = le_fluide.beta_t();
833 const DoubleTab& tab_champ_beta_t = ch_beta_t.valeurs();
834 if (sub_type(Champ_Uniforme,ch_beta_t))
835 {
836 beta_t = std::max(tab_champ_beta_t(0,0),DMINFLOAT);
837 }
838 else
839 {
840 Cerr << " On ne sait pas traiter un beta_t non uniforme dans Paroi TBLE_LRM !!!"<< finl;
841 Process::exit();
842 }
843
844
845
846 for(int i=0 ; i<nb_pts ; i++)
847 {
848 grad_T(i)=(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i+1)-loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i))/(eq_k_U_W[compteur_faces_paroi].get_y(i+1)-eq_k_U_W[compteur_faces_paroi].get_y(i));
849 }
850
851 for(int i=1 ; i<nb_pts ; i++)
852 {
853 grad_T_moy(i)=0.5*(grad_T(i)+grad_T(i-1));
854 }
855
856
857 grad_T_moy(0)=(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,1)-loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,0))/(eq_k_U_W[compteur_faces_paroi].get_y(1)-eq_k_U_W[compteur_faces_paroi].get_y(0));
858
859 const Champ_Don_base& ch_gravite=le_fluide.gravite();
860 const DoubleVect& gravite = ch_gravite.valeurs();
861
862 if (!sub_type(Champ_Uniforme,ch_gravite))
863 {
864 Cerr << " On ne sait pas traiter la gravite non uniforme dans Paroi TBLE_LRM !!!"<< finl;
865 Process::exit();
866 }
867
868
869
870 //calcul de la gravite
871
872 double norm_n=domaine_VDF.face_surfaces(num_face);
873 double gn=0.;
874
875 for (int idim=0; idim<dimension; idim++)
876 {
877 gn+=gravite(idim)*domaine_VDF.face_normales(num_face,idim)/norm_n;
878 }
879
880 DoubleVect gt_vect(dimension);
881
882 for (int idim=0; idim<dimension; idim++)
883 {
884 gt_vect(idim) = gravite(idim)-gn*domaine_VDF.face_normales(num_face,idim)/norm_n;
885 }
886
887 double g_t=gt_vect(1-ori);
888
889
890
891 for(int i=0 ; i<nb_pts; i++)
892 {
893 source_T_k(i) = g_t*beta_t*((diffu.calculer_a_local(i)-d_visco)/(0.9))*grad_T_moy(i);
894 source_T_U(i) = g_t*beta_t*(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i)-T0);
895 }
896
897 }
898 }
899
900 //terme source a la paroi
901
902 eps0=2*d_visco*eq_k_U_W[compteur_faces_paroi].get_Unp1(1,1)/(eq_k_U_W[compteur_faces_paroi].get_y(1)*eq_k_U_W[compteur_faces_paroi].get_y(1));
903 eq_k_U_W[compteur_faces_paroi].set_F(1,0,-eps0 + source_T_k(0));
904 eq_k_U_W[compteur_faces_paroi].set_F(0,0, -gradient_de_pression0 + ts0 - source_T_U(0));
905
906 for(int i=1 ; i<nb_pts; i++)
907 {
908
909 //echelles de vitesse et de longueur (u ne sert pas actuellement, il est mis a 0)
910 double k=eq_k_U_W[compteur_faces_paroi].get_Unp1(1,i);
911 double y=eq_k_U_W[compteur_faces_paroi].get_y(i);
912
913 vv_bar=le_mod_ech.calculer_vv_bar(y, k, u, d_visco);
914 l_eps=le_mod_ech.calculer_l_eps(y, k, u, d_visco);
915
916
917 //terme source en i pour k
918 eps=eq_k_U_W[compteur_faces_paroi].get_Unp1(1,i)*sqrt(vv_bar)/l_eps;
919
920 eq_k_U_W[compteur_faces_paroi].set_F(1,i,(diffu.calculer_a_local(i)-d_visco)
921 *grad_vit_elem_moyx(i)*grad_vit_elem_moyx(i)-eps+source_T_k(i));
922
923 //terme source pour U
924 eq_k_U_W[compteur_faces_paroi].set_F(0,i, -gradient_de_pression0 + ts0 - source_T_U(i) );
925
926 }
927
928 //On resoud les equations aux limites simplifiees
929 //pendant un pas de temps du maillage grossier
930 // On obtient U, V, W au temps n+1
931
932 eq_k_U_W[compteur_faces_paroi].aller_au_temps(tps);
933
934
935 //Message lorsque k est negatif. Modification de la valeur a 1e-9
936 for(int i=0 ; i<nb_pts ; i++)
937 {
938 if (eq_k_U_W[compteur_faces_paroi].get_Unp1(1,i)<0.)
939 {
940 Cerr << "Warning dans TBLE_LRM : k <0 !!! k = " << eq_k_U_W[compteur_faces_paroi].get_Unp1(1,i) << finl;
941 eq_k_U_W[compteur_faces_paroi].set_Unp1(1,i,1.e-9);
942 }
943 }
944
945 if(itmax<eq_k_U_W[compteur_faces_paroi].get_it())
946 {
947 itmax = eq_k_U_W[compteur_faces_paroi].get_it();
948 if(itmax>max_it)
949 Cerr << "WARNING : TOO MANY ITERATIONS ARE NEEDED IN THE TBLE MESH"
950 << finl;
951 }
952
953
954 ///////////////////////////////////////////////
955 ///// Determination des deux composantes //////
956 ///// du cisaillement a la paroi //////
957 ///////////////////////////////////////////////
958
959 if (ori == 0)
960 {
961 //Cerr << "ori=0" << finl;
962 Cisaillement_paroi_(num_face,0) = 0.;
963 Cisaillement_paroi_(num_face,1) = eq_k_U_W[compteur_faces_paroi].get_cis(0)*signe;
964 }
965 else if (ori == 1)
966 {
967
968 Cisaillement_paroi_(num_face,0) = eq_k_U_W[compteur_faces_paroi].get_cis(0)*signe;
969 Cisaillement_paroi_(num_face,1) = 0.;
970 }
971
972 tab_u_star_(num_face) =pow(eq_k_U_W[compteur_faces_paroi].get_cis(0)*eq_k_U_W[compteur_faces_paroi].get_cis(0),0.25);
973 tab_d_plus_(num_face) = dist*tab_u_star(num_face)/d_visco;
974
975 compteur_faces_paroi++;
976
977
978 }//fin boucle sur les faces de bords parietaux
979
980
981 }//fin if(dim2)
982
983
984 ///////////////////////
985 // Si probleme 3D //
986 ///////////////////////
987
988 else if(dimension == 3)
989 {
990 //eq_k_U_W.dimensionner(2*le_dom_dis_->nb_faces_bord());
991
992 //Boucle sur les faces des bords parietaux
993
994 for (int num_face=ndeb; num_face<nfin; num_face++)
995 {
996 //ici ori=direction perpendiculaire a la paroi
997 ori = orientation(num_face);
998 //Cerr << "ori = " << ori << finl;
999
1000 //compteur++;
1001 //Selection de l'element associe a une face parietale
1002 //et initialisation de l'entier signe
1003
1004 if ((elem = face_voisins(num_face,0)) == -1)
1005 {
1006 elem = face_voisins(num_face,1);
1007 signe = 1;
1008 }
1009 else signe = -1 ;
1010
1011
1012 if(dt<dt_min)
1013 eq_k_U_W[compteur_faces_paroi].set_dt(1.e12);
1014 else
1015 eq_k_U_W[compteur_faces_paroi].set_dt(dt);
1016
1017
1018 // 1er couple de faces perpendiculaires a la paroi
1019 int face1 = elem_faces(elem,(ori+1));
1020 int face2 = elem_faces(elem,(ori+4)%6);
1021 // 2e couple de faces perpendiculaires a la paroi
1022 int face3 = elem_faces(elem,(ori+2));
1023 int face4 = elem_faces(elem,(ori+5)%6);
1024
1025
1026
1027
1028
1029
1030 ////////////////////////////////////////////////////////////
1031 ////// 1er couple de faces perpendiculaires a la paroi /////
1032 ////////////////////////////////////////////////////////////
1033
1034 vmoy = 0.5*(vit(face1) + vit(face2));
1035
1036 //vitesse sur le maillage fin a l'instant initial
1037
1038
1039 ts0 = 0.5*(termes_sources(face1)/volumes_entrelaces(face1)
1040 + termes_sources(face2)/volumes_entrelaces(face2));
1041
1042 eq_k_U_W[compteur_faces_paroi].set_u_y0(0,0.); //1ere composante de la vitesse a la paroi
1043
1044 eq_k_U_W[compteur_faces_paroi].set_u_yn(0,vmoy); //vitesse en yn
1045
1046 //*** Gradient de pression du 1er couple de faces
1047 // perpendiculaires a la paroi ***
1048
1049 //gradient de pression (constant dans le maillage fin)
1050 gradient_de_pression0 = 0.5*(grad_p(face1)+grad_p(face2));
1051
1052 //eq_k_U_W[compteur_faces_paroi].set_F(0, gradient_de_pression0 - ts0);
1053
1054
1055 //dpzbis += (gradient_de_pression_bis);
1056
1057 ////////////////////////////////////////////////////////////
1058 ////// 2e couple de faces perpendiculaires a la paroi /////
1059 ////////////////////////////////////////////////////////////
1060
1061 vmoy = 0.5*(vit(face3) + vit(face4));
1062
1063 ts1 = 0.5*(termes_sources(face3)/volumes_entrelaces(face3)
1064 + termes_sources(face4)/volumes_entrelaces(face4));
1065
1066 eq_k_U_W[compteur_faces_paroi].set_u_y0(1,0.); //2e composante de la vitesse a la paroi
1067
1068 eq_k_U_W[compteur_faces_paroi].set_u_yn(1,vmoy); //vitesse en yn
1069
1070 //*** Gradient de pression du 2e couple de faces
1071 // perpendiculaires a la paroi ***
1072
1073
1074 //gradient de pression (constant dans le maillage fin)
1075
1076 gradient_de_pression1 = 0.5*(grad_p(face3)+grad_p(face4));
1077
1078
1079 //eq_k_U_W[compteur_faces_paroi].set_F(1, gradient_de_pression1 - ts1);
1080 //Calcul du gradient On prend U au temps n
1081
1082 //equation en k ;
1083
1084
1085 double d_visco;
1086
1087 if (l_unif==1)
1088 d_visco = visco;
1089 else
1090 d_visco = tab_visco[elem];
1091
1092
1093 Diffu_totale_base& diffu = eq_k_U_W[compteur_faces_paroi].get_diffu();
1094
1095
1096 for(int i=0 ; i<nb_pts; i++)
1097 {
1098
1099 grad_vit_elemx(i)=(eq_k_U_W[compteur_faces_paroi].get_Unp1(0,i+1)-eq_k_U_W[compteur_faces_paroi].get_Unp1(0,i))/(eq_k_U_W[compteur_faces_paroi].get_y(i+1)-eq_k_U_W[compteur_faces_paroi].get_y(i));
1100 grad_vit_elemz(i)=(eq_k_U_W[compteur_faces_paroi].get_Unp1(1,i+1)-eq_k_U_W[compteur_faces_paroi].get_Unp1(1,i))/(eq_k_U_W[compteur_faces_paroi].get_y(i+1)-eq_k_U_W[compteur_faces_paroi].get_y(i));
1101
1102 }
1103
1104 //gradient en 0
1105 grad_vit_elem_moyx(0)=(eq_k_U_W[compteur_faces_paroi].get_Unp1(0,1)-eq_k_U_W[compteur_faces_paroi].get_Unp1(0,0))/(eq_k_U_W[compteur_faces_paroi].get_y(1)-eq_k_U_W[compteur_faces_paroi].get_y(0));
1106 grad_vit_elem_moyz(0)=(eq_k_U_W[compteur_faces_paroi].get_Unp1(1,1)-eq_k_U_W[compteur_faces_paroi].get_Unp1(1,0))/(eq_k_U_W[compteur_faces_paroi].get_y(1)-eq_k_U_W[compteur_faces_paroi].get_y(0));
1107
1108 //Cacul du gradient a l'ordonnee y(i)
1109 for(int i=1 ; i<nb_pts ; i++)
1110 {
1111
1112 grad_vit_elem_moyx(i)=0.5*(grad_vit_elemx(i)+grad_vit_elemx(i-1));
1113 grad_vit_elem_moyz(i)=0.5*(grad_vit_elemz(i)+grad_vit_elemz(i-1));
1114
1115 }
1116
1117 // On calcul vv_bar, l_eps, y_nu, y_star avec k au temps n
1118
1119 double dist;
1120
1121 //Distance a la paroi du 1er centre de maille
1122 if (axi)
1123 dist = domaine_VDF.dist_norm_bord_axi(num_face);
1124 else
1125 dist = domaine_VDF.dist_norm_bord(num_face);
1126
1127
1128 // Definition de la face servant a recuperer le k de la deuxieme maille
1129 int face5 = elem_faces(elem,(ori+3)%6);
1130 int face6 = elem_faces(elem,(ori));
1131
1132
1133
1134 int n2=face5;
1135 if (face5==num_face)
1136 {
1137 n2=face6;
1138 }
1139 int elem1=face_voisins(n2,0);
1140 if (elem1 == elem)
1141 elem1=face_voisins(n2,1);
1142
1143
1144 double dy = domaine_VDF.dist_norm(n2);
1145 // recuperartion de l'element correspondant a la deuxieme maille
1146
1147 double kcl;
1148 double a1=d_visco+0.5*(nu_t(elem)+nu_t(elem1));
1149 double a2=0.5*(diffu.calculer_a_local(nb_pts)+diffu.calculer_a_local(nb_pts-1));
1150
1151 double b1=a1/dy;
1152 double b2=a2/(eq_k_U_W[compteur_faces_paroi].get_y(nb_pts)-eq_k_U_W[compteur_faces_paroi].get_y(nb_pts-1));
1153
1154
1155
1156 kcl=(b1*tab_k(elem1)+b2*eq_k_U_W[compteur_faces_paroi].get_Unp1(2,nb_pts-1))/(b1+b2);
1157
1158
1159 //initialisation de epsilon pour maillage 3D tel que nut(tble)=nut(3D)
1160
1161 if(tab_k(elem)<0)
1162 {
1163 //Cerr << "Warning dans TBLE_LRM : k3D <0 !!! k = " << tab_k_eps(elem,0) << finl;
1164 tab_k(elem) =1.e-8;
1165 tab_eps(elem)=0.09*tab_k(elem)*tab_k(elem)/(diffu.calculer_a_local(nb_pts)-d_visco);
1166 }
1167
1168 else
1169 {
1170 tab_k(elem)=kcl;
1171 tab_eps(elem)=0.09*tab_k(elem)*tab_k(elem)/(diffu.calculer_a_local(nb_pts)-d_visco);
1172 }
1173
1174 double T0=0;
1175 double beta_t=0.;
1176
1177 source_T_k=0.;
1178 source_T_U=0.;
1179
1180 if (avec_boussi==1)
1181 {
1182 int boussi_ok=0;
1183 const Sources& les_sources=eqnNS.sources();
1184 int nb_sources=les_sources.size();
1185
1186 for (int j=0; j<nb_sources; j++)
1187 {
1188 const Source_base& ts = les_sources(j).valeur();
1189 if (sub_type(Terme_Boussinesq_base,ts))
1190 {
1191 const Terme_Boussinesq_base& terme_boussi = ref_cast(Terme_Boussinesq_base,ts);
1192 T0 = terme_boussi.Scalaire0(0);
1193 boussi_ok=1;
1194 }
1195 }
1196
1197 if (boussi_ok==1)
1198 {
1199 const Equation_base& eqn_th = mon_modele_turb_hyd->equation().probleme().equation(1);
1200 const Modele_turbulence_scal_base& le_mod_turb_th = ref_cast(Modele_turbulence_scal_base,eqn_th.get_modele(TURBULENCE).valeur());
1201 const Turbulence_paroi_scal_base& loi = le_mod_turb_th.loi_paroi();
1202 Paroi_TBLE_scal_VDF& loi_tble_T = ref_cast_non_const(Paroi_TBLE_scal_VDF,loi);
1203
1204 const Champ_Don_base& ch_beta_t = le_fluide.beta_t();
1205 const DoubleTab& tab_champ_beta_t = ch_beta_t.valeurs();
1206 if (sub_type(Champ_Uniforme,ch_beta_t))
1207 {
1208 beta_t = std::max(tab_champ_beta_t(0,0),DMINFLOAT);
1209 }
1210 else
1211 {
1212 Cerr << " On ne sait pas traiter un beta_t non uniforme dans Paroi TBLE_LRM !!!"<< finl;
1213 Process::exit();
1214 }
1215
1216
1217
1218 for(int i=0 ; i<nb_pts ; i++)
1219 {
1220 grad_T(i)=(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i+1)-loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i))/(eq_k_U_W[compteur_faces_paroi].get_y(i+1)-eq_k_U_W[compteur_faces_paroi].get_y(i));
1221 }
1222
1223 for(int i=1 ; i<nb_pts ; i++)
1224 {
1225 grad_T_moy(i)=0.5*(grad_T(i)+grad_T(i-1));
1226 }
1227
1228
1229 grad_T_moy(0)=(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,1)-loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,0))/(eq_k_U_W[compteur_faces_paroi].get_y(1)-eq_k_U_W[compteur_faces_paroi].get_y(0));
1230
1231 const Champ_Don_base& ch_gravite=le_fluide.gravite();
1232 const DoubleVect& gravite = ch_gravite.valeurs();
1233
1234 if (!sub_type(Champ_Uniforme,ch_gravite))
1235 {
1236 Cerr << " On ne sait pas traiter la gravite non uniforme dans Paroi TBLE_LRM !!!"<< finl;
1237 Process::exit();
1238 }
1239
1240
1241
1242 //calcul de la gravite
1243
1244
1245 double norm_n=domaine_VDF.face_surfaces(num_face);
1246 double gn=0.;
1247
1248 for (int idim=0; idim<dimension; idim++)
1249 {
1250 gn+=gravite(idim)*domaine_VDF.face_normales(num_face,idim)/norm_n;
1251 }
1252
1253 DoubleVect gt_vect(dimension);
1254
1255 for (int idim=0; idim<dimension; idim++)
1256 {
1257 gt_vect(idim) = gravite(idim)-gn*domaine_VDF.face_normales(num_face,idim)/norm_n;
1258 }
1259
1260 double g_t=gt_vect(1-ori);
1261
1262 for(int i=1 ; i<nb_pts; i++)
1263 {
1264 source_T_k(i) =g_t*beta_t*((diffu.calculer_a_local(i)-d_visco)/(0.9))*grad_T_moy(i);
1265 source_T_U(i) =g_t*beta_t*(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i)-T0);
1266
1267 }
1268
1269 }
1270 }
1271
1272 eps0=2*d_visco*eq_k_U_W[compteur_faces_paroi].get_Unp1(2,1)/(eq_k_U_W[compteur_faces_paroi].get_y(1)*eq_k_U_W[compteur_faces_paroi].get_y(1));
1273 eq_k_U_W[compteur_faces_paroi].set_F(2,0,-eps0 + source_T_k(0));
1274 eq_k_U_W[compteur_faces_paroi].set_F(0,0, -gradient_de_pression0 + ts0 - source_T_U(0));
1275
1276 eq_k_U_W[compteur_faces_paroi].set_F(1,0, -gradient_de_pression1 + ts1 - source_T_U(0) );
1277
1278
1279 for(int i=1 ; i<nb_pts ; i++)
1280 {
1281
1282 //echelles de vitesse et de longueur (u ne sert pas actuellement, il est mis a 0)
1283 double k=eq_k_U_W[compteur_faces_paroi].get_Unp1(2,i);
1284 double y=eq_k_U_W[compteur_faces_paroi].get_y(i);
1285
1286
1287 vv_bar=le_mod_ech.calculer_vv_bar(y, k, u, d_visco);
1288 l_eps=le_mod_ech.calculer_l_eps(y, k, u, d_visco);
1289
1290 //terme source a la paroi
1291
1292
1293 eps=eq_k_U_W[compteur_faces_paroi].get_Unp1(2,i)*sqrt(vv_bar)/l_eps;
1294
1295 eq_k_U_W[compteur_faces_paroi].set_F(2,i,(diffu.calculer_a_local(i)-d_visco)*(grad_vit_elem_moyx(i)*grad_vit_elem_moyx(i)+grad_vit_elem_moyz(i)*grad_vit_elem_moyz(i))-eps
1296 + source_T_k(i));
1297
1298 // pour U et W
1299
1300 eq_k_U_W[compteur_faces_paroi].set_F(0,i, -gradient_de_pression0 + ts0 - source_T_U(i));
1301 eq_k_U_W[compteur_faces_paroi].set_F(1,i, -gradient_de_pression1 + ts1 - source_T_U(i) );
1302
1303 }
1304
1305
1306
1307 //On resoud les equations aux limites simplifiees
1308 //pendant un pas de temps du maillage grossier
1309
1310 eq_k_U_W[compteur_faces_paroi].set_u_y0(2,0.); //k a la paroi
1311 eq_k_U_W[compteur_faces_paroi].set_u_yn(2,kcl);//k a yn
1312
1313 eq_k_U_W[compteur_faces_paroi].aller_au_temps(tps);
1314
1315 //Message lorsque k est negatif. Modification de la valeur a 1e-9
1316 for(int i=0 ; i<nb_pts ; i++)
1317 {
1318 if (eq_k_U_W[compteur_faces_paroi].get_Unp1(2,i)<0.)
1319 {
1320 // Cerr << "Warning dans TBLE_LRM : k <0 !!! k = " << eq_k_U_W[compteur_faces_paroi].get_Unp1(2,i) << finl;
1321 eq_k_U_W[compteur_faces_paroi].set_Unp1(2,i,1.e-9);
1322 }
1323 }
1324
1325 if(itmax<eq_k_U_W[compteur_faces_paroi].get_it())
1326 {
1327 itmax = eq_k_U_W[compteur_faces_paroi].get_it();
1328 if(itmax>max_it)
1329 Cerr << "WARNING : TOO MANY ITERATIONS ARE NEEDED IN THE TBLE MESH"
1330 << finl;
1331 }
1332
1333
1334 ///////////////////////////////////////////////
1335 ///// Determination des deux composantes //////
1336 ///// du cisaillement a la paroi //////
1337 ///////////////////////////////////////////////
1338
1339 if (ori == 0)
1340 {
1341 //Cerr << "ori=0" << finl;
1342 Cisaillement_paroi_(num_face,0) = 0.;
1343 Cisaillement_paroi_(num_face,1) = eq_k_U_W[compteur_faces_paroi].get_cis(0)*signe;
1344 Cisaillement_paroi_(num_face,2) = eq_k_U_W[compteur_faces_paroi].get_cis(1)*signe;
1345 }
1346 else if (ori == 1)
1347 {
1348
1349 Cisaillement_paroi_(num_face,0) = eq_k_U_W[compteur_faces_paroi].get_cis(1)*signe;
1350 Cisaillement_paroi_(num_face,1) = 0.;
1351 Cisaillement_paroi_(num_face,2) = eq_k_U_W[compteur_faces_paroi].get_cis(0)*signe;
1352
1353 }
1354 else
1355 {
1356 //Cerr << "ori=2" << finl;
1357 Cisaillement_paroi_(num_face,0) = eq_k_U_W[compteur_faces_paroi].get_cis(0)*signe;
1358 Cisaillement_paroi_(num_face,1) = eq_k_U_W[compteur_faces_paroi].get_cis(1)*signe;
1359 Cisaillement_paroi_(num_face,2) = 0.;
1360 }
1361
1362 tab_u_star_(num_face) = pow((eq_k_U_W[compteur_faces_paroi].get_cis(0)*eq_k_U_W[compteur_faces_paroi].get_cis(0)
1363 +eq_k_U_W[compteur_faces_paroi].get_cis(1)*eq_k_U_W[compteur_faces_paroi].get_cis(1)),0.25);
1364 tab_d_plus_(num_face) = dist*tab_u_star(num_face)/d_visco;
1365
1366
1367 compteur_faces_paroi++;
1368
1369
1370 for(int j=0; j<nb_post_pts; j++)
1371 {
1372 Nom tmp2;
1373 tmp2="nut_post_";
1374 tmp2+=nom_pts[j];
1375 tmp2+=".dat";
1376
1377 SFichier fic_nut(tmp2); //creation fichier de post-traitement des points TBLE
1378
1379 Diffu_totale_base& diffu1 = eq_k_U_W[num_faces_post(j)].get_diffu();
1380 for(int i=0 ; i<nb_pts+1 ; i++)
1381 {
1382 fic_nut << eq_k_U_W[num_faces_post(j)].get_y(i)<< " " << diffu1.calculer_a_local(i)-d_visco << finl;
1383 }
1384 fic_nut << finl;
1385
1386 }
1387
1388 }//fin boucle sur les faces de bords parietaux
1389
1390
1391
1392 }//fin if(dim=3)
1393 }//fin if(sub_type(dirichlet))
1394
1395 }//fin boucle sur les bords
1396
1397 if (Process::je_suis_maitre())
1398 {
1399 SFichier fic("iter.dat",ios::app); // ouverture du fichier iter.dat
1400 fic << tps << " " << itmax << finl;
1401 }
1402
1403
1404 Cisaillement_paroi_.echange_espace_virtuel();
1405
1406 return 1;
1407}
1408
1409int ParoiVDF_TBLE_LRM::calculer_hyd(DoubleTab& tab_k_eps)
1410{
1411
1412
1413 if (Process::nproc()>1)
1414 {
1415 Cerr << "ParoiVDF_TBLE_LRM::calculer_hyd n'est pas parallelise." << finl;
1416 Process::exit();
1417 }
1418 const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_dis_.valeur());
1419 const IntVect& orientation = domaine_VDF.orientation();
1420 const IntTab& face_voisins = domaine_VDF.face_voisins();
1421 const IntTab& elem_faces = domaine_VDF.elem_faces();
1422 const DoubleVect& volumes_entrelaces = domaine_VDF.volumes_entrelaces();
1423 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
1424 const Fluide_base& le_fluide = ref_cast(Fluide_base, eqn_hydr.milieu());
1425 const Champ_Don_base& ch_visco_cin = le_fluide.viscosite_cinematique();
1426 const DoubleTab& tab_visco = ch_visco_cin.valeurs();
1427 Mod_echelle_LRM_base& le_mod_ech= mod_ech.valeur();
1428
1429 const DoubleTab& nu_t = mon_modele_turb_hyd->viscosite_turbulente().valeurs();
1430
1431
1432
1433 double visco=-1;
1434 int l_unif;
1435 if (sub_type(Champ_Uniforme,ch_visco_cin))
1436 {
1437 visco = std::max(tab_visco(0,0),DMINFLOAT);
1438 l_unif = 1;
1439 }
1440 else
1441 l_unif = 0;
1442
1443 // preparer_calcul_hyd(tab);
1444 if ((!l_unif) && (tab_visco.local_min_vect()<DMINFLOAT))
1445 // on ne doit pas changer tab_visco ici !
1446 {
1447 Cerr << "In ParoiVDF_TBLE_LRM::calculer_hyd : visco = " << tab_visco.local_min_vect() << " <= 0 ? " << finl;
1448 throw;
1449 }
1450 // tab_visco+=DMINFLOAT;
1451
1452
1453 int ndeb=0,nfin=0;
1454 int ori;
1455 // int elem_gauche0, elem_gauche1, elem_droite0, elem_droite1;
1456 // int face_gauche0, face_gauche1, face_droite0, face_droite1;
1457 int elem;
1458 // int compteur=0;
1459
1460 int signe;
1461 int itmax=0;
1462
1463 double vmoy = 0., ts0 =0., ts1=0.,gradient_de_pression0=0., gradient_de_pression1 = 0.;
1464 const DoubleTab& vit = eqn_hydr.inconnue().valeurs(); //vitesse
1465 DoubleVect grad_vit_elemx;
1466 DoubleVect grad_vit_elem_moyx;
1467 DoubleVect grad_vit_elemz;
1468 DoubleVect grad_vit_elem_moyz;
1469 DoubleVect grad_T;
1470 DoubleVect grad_T_moy;
1471 DoubleVect source_T_k;
1472 DoubleVect source_T_U;
1473
1474
1475 double eps0;
1476 double eps;
1477 double l_eps;
1478 double vv_bar;
1479 double u=0;
1480
1481 grad_vit_elemx.resize(nb_pts);
1482 grad_vit_elem_moyx.resize(nb_pts);
1483 grad_vit_elemz.resize(nb_pts);
1484 grad_vit_elem_moyz.resize(nb_pts);
1485 grad_T.resize(nb_pts);
1486 grad_T_moy.resize(nb_pts);
1487 source_T_k.resize(nb_pts);
1488 source_T_U.resize(nb_pts);
1489
1490 const Navier_Stokes_std& eqnNS = ref_cast(Navier_Stokes_std, eqn_hydr);
1491
1492 DoubleTab grad_p(vit);
1493 const DoubleTab& p = eqnNS.pression().valeurs();
1494 const Operateur_Grad& gradient = eqnNS.operateur_gradient();
1495 gradient.calculer(p, grad_p); // Calcul du gradient de pression
1496
1497 DoubleTab termes_sources;
1498 termes_sources.resize(domaine_VDF.nb_faces(),1);
1499 eqn_hydr.sources().calculer(termes_sources); //les termes sources
1500
1501 const double tps = eqnNS.schema_temps().temps_courant();
1502 const double dt = eqnNS.schema_temps().pas_de_temps();
1503 const double dt_min = eqnNS.schema_temps().pas_temps_min();
1504
1505
1506 //SFichier fic_v("v.dat",ios::app); // impression de la vitesse normale pour test
1507
1508 int compteur_faces_paroi = 0;
1509 DoubleVect corresp(le_dom_dis_->nb_faces_bord());
1510
1511
1512 // Boucle sur les bords
1513 //Cerr << "Boucle sur les bords" << finl;
1514
1515 for (int n_bord=0; n_bord<domaine_VDF.nb_front_Cl(); n_bord++)
1516 {
1517 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
1518
1519 if (sub_type(Dirichlet_paroi_fixe,la_cl.valeur()) )
1520
1521 {
1522 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1523 ndeb = le_bord.num_premiere_face();
1524 nfin = ndeb + le_bord.nb_faces();
1525
1526 for (int num_face=ndeb; num_face<nfin; num_face++)
1527 {
1528 corresp(compteur_faces_paroi)=num_face;
1529 compteur_faces_paroi++;
1530 }
1531 }
1532 }
1533
1534 corresp.resize(compteur_faces_paroi); //Redimensionnement de corresp
1535 compteur_faces_paroi = 0; //Reinitialisation de compteur_faces_paroi
1536
1537
1538 // Boucle sur les bords
1539
1540 //Cerr << "Boucle sur les bords " << finl;
1541 for (int n_bord=0; n_bord<domaine_VDF.nb_front_Cl(); n_bord++)
1542 {
1543
1544 // pour chaque condition limite on regarde son type
1545 // On applique les lois de paroi uniquement
1546 // aux voisinages des parois
1547
1548
1549 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
1550
1551 if (sub_type(Dirichlet_paroi_fixe,la_cl.valeur()) )
1552
1553 {
1554 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1555 ndeb = le_bord.num_premiere_face();
1556 nfin = ndeb + le_bord.nb_faces();
1557
1558
1559
1560 ///////////////////////
1561 // Si probleme 2D //
1562 ///////////////////////
1563
1564
1565 if(dimension == 2)
1566 {
1567 //eq_k_U_W.dimensionner(2*le_dom_dis_->nb_faces_bord());
1568
1569 //Boucle sur les faces des bords parietaux
1570
1571 for (int num_face=ndeb; num_face<nfin; num_face++)
1572 {
1573 //ici ori=direction perpendiculaire a la paroi
1574 ori = orientation(num_face);
1575 //Cerr << "ori = " << ori << finl;
1576
1577
1578
1579 //compteur++;
1580 //Selection de l'element associe a une face parietale
1581 //et initialisation de l'entier signe
1582
1583 if ((elem = face_voisins(num_face,0)) == -1)
1584 {
1585 elem = face_voisins(num_face,1);
1586 signe = 1;
1587 }
1588 else signe = -1 ;
1589
1590
1591 if(dt<dt_min)
1592 eq_k_U_W[compteur_faces_paroi].set_dt(1.e12);
1593 else
1594 eq_k_U_W[compteur_faces_paroi].set_dt(dt);
1595
1596
1597
1598 // 1er couple de faces perpendiculaires a la paroi
1599 int face1 = elem_faces(elem,(ori+1));
1600 int face2 = elem_faces(elem,(ori+3)%4);
1601
1602
1603
1604 ////////////////////////////////////////////////////////////
1605 ////// couple de faces perpendiculaires a la paroi /////
1606 ////////////////////////////////////////////////////////////
1607
1608 vmoy = 0.5*(vit(face1) + vit(face2));
1609
1610 //vitesse sur le maillage fin a l'instant initial
1611 //eq_k_U_W[compteur_faces_paroi].set_Uinit_lin(0, 0., vmoy);
1612
1613 ts0 = 0.5*(termes_sources(face1)/volumes_entrelaces(face1)
1614 + termes_sources(face2)/volumes_entrelaces(face2));
1615
1616 eq_k_U_W[compteur_faces_paroi].set_u_y0(0,0.); //1ere composante de la vitesse a la paroi
1617 eq_k_U_W[compteur_faces_paroi].set_u_yn(0,vmoy); //vitesse en yn
1618 // eq_k_U_W[compteur_faces_paroi].set_u_yn(0,1.); //vitesse en yn
1619
1620
1621
1622 //*** Gradient de pression du 1er couple de faces
1623 // perpendiculaires a la paroi ***
1624
1625 //gradient de pression (constant dans le maillage fin)
1626 gradient_de_pression0 = 0.5*(grad_p(face1)+grad_p(face2));
1627
1628 //Calcul du gradient On prend U au temps n
1629 for(int i=0 ; i<nb_pts ; i++)
1630 {
1631 grad_vit_elemx(i)=(eq_k_U_W[compteur_faces_paroi].get_Unp1(0,i+1)-eq_k_U_W[compteur_faces_paroi].get_Unp1(0,i))/(eq_k_U_W[compteur_faces_paroi].get_y(i+1)-eq_k_U_W[compteur_faces_paroi].get_y(i));
1632 }
1633
1634 //Cacul du gradient a l'ordonnee y(i)
1635 for(int i=1 ; i<nb_pts ; i++)
1636 {
1637 grad_vit_elem_moyx(i)=0.5*(grad_vit_elemx(i)+grad_vit_elemx(i-1));
1638 }
1639
1640 //gradient en 0
1641 grad_vit_elem_moyx(0)=(eq_k_U_W[compteur_faces_paroi].get_Unp1(0,1)-eq_k_U_W[compteur_faces_paroi].get_Unp1(0,0))/(eq_k_U_W[compteur_faces_paroi].get_y(1)-eq_k_U_W[compteur_faces_paroi].get_y(0));
1642
1643
1644
1645 double d_visco;
1646
1647 if (l_unif==1)
1648 d_visco = visco;
1649 else
1650 d_visco = tab_visco[elem];
1651
1652 // On calcul vv_bar, l_eps, y_nu, y_star avec k au temps n
1653 // a la paroi eps est proportionelle a 1
1654
1655 Diffu_totale_base& diffu = eq_k_U_W[compteur_faces_paroi].get_diffu();
1656
1657 double dist;
1658 if (axi)
1659 dist = domaine_VDF.dist_norm_bord_axi(num_face);
1660 else
1661 dist = domaine_VDF.dist_norm_bord(num_face);
1662
1663
1664
1665 // Definition de la face servant a recuperer le k de la deuxieme maille
1666 int face3 = elem_faces(elem,(ori+2)%4);
1667 int face4 = elem_faces(elem,ori);
1668
1669 int n2=face3;
1670 if (face3==num_face)
1671 {
1672 n2=face4;
1673 }
1674
1675 int elem1=face_voisins(n2,0);
1676 if (elem1 == elem)
1677 elem1=face_voisins(n2,1);
1678
1679
1680
1681 double dy = domaine_VDF.dist_norm(n2);
1682
1683
1684 // recuperartion de l'element correspondant a la deuxieme maille
1685
1686 double kcl;
1687
1688 double a1=d_visco+0.5*(nu_t(elem)+nu_t(elem1));
1689 double a2=0.5*(diffu.calculer_a_local(nb_pts)+diffu.calculer_a_local(nb_pts-1));
1690
1691
1692 double b1=a1/dy;
1693 double b2=a2/(eq_k_U_W[compteur_faces_paroi].get_y(nb_pts)-eq_k_U_W[compteur_faces_paroi].get_y(nb_pts-1));
1694
1695
1696 kcl=(b1*tab_k_eps(elem1,0)+b2*eq_k_U_W[compteur_faces_paroi].get_Unp1(1,nb_pts-1))/(b1+b2);
1697
1698
1699 //equation en k ;
1700
1701 eq_k_U_W[compteur_faces_paroi].set_u_y0(1,0.); //k a la paroi
1702 eq_k_U_W[compteur_faces_paroi].set_u_yn(1,kcl);//k a yn
1703
1704
1705 // double Re, Fmu;
1706
1707 //initialisation de epsilon pour maillage 2D tel que nut(tble)=nut(2D)
1708
1709
1710 /* if(tab_k_eps(elem,0)<0)
1711 {
1712 // Cerr << "Warning dans TBLE_LRM : k3D <0 !!! k = " << tab_k_eps(elem,0) << finl;
1713 tab_k_eps(elem,0)=1.e-8;
1714 tab_k_eps(elem,1)=0.09*tab_k_eps(elem,0)*tab_k_eps(elem,0)/(diffu.calculer_a_local(nb_pts)-d_visco);
1715 } */
1716
1717 if(kcl<0)
1718 {
1719 // Cerr << "Warning dans TBLE_LRM : k3D <0 !!! k = " << tab_k_eps(elem,0) << finl;
1720 kcl=1.e-8;
1721 tab_k_eps(elem,0)=kcl;
1722 tab_k_eps(elem,1)=0.09*tab_k_eps(elem,0)*tab_k_eps(elem,0)/(diffu.calculer_a_local(nb_pts)-d_visco);
1723 }
1724
1725 else
1726 {
1727 tab_k_eps(elem,0)=kcl;
1728 tab_k_eps(elem,1)= 0.09*tab_k_eps(elem,0)*tab_k_eps(elem,0)/(diffu.calculer_a_local(nb_pts)-d_visco);
1729 }
1730
1731
1732 double T0=0;
1733 double beta_t=0.;
1734
1735 source_T_k=0.;
1736 source_T_U=0.;
1737
1738 if (avec_boussi==1)
1739 {
1740 int boussi_ok=0;
1741 const Sources& les_sources=eqnNS.sources();
1742 int nb_sources=les_sources.size();
1743
1744 for (int j=0; j<nb_sources; j++)
1745 {
1746 const Source_base& ts = les_sources(j).valeur();
1747 if (sub_type(Terme_Boussinesq_base,ts))
1748 {
1749 const Terme_Boussinesq_base& terme_boussi = ref_cast(Terme_Boussinesq_base,ts);
1750 T0 = terme_boussi.Scalaire0(0);
1751 boussi_ok=1;
1752 }
1753 }
1754
1755 if (boussi_ok==1)
1756 {
1757 const Equation_base& eqn_th = mon_modele_turb_hyd->equation().probleme().equation(1);
1758 const Modele_turbulence_scal_base& le_mod_turb_th = ref_cast(Modele_turbulence_scal_base,eqn_th.get_modele(TURBULENCE).valeur());
1759 const Turbulence_paroi_scal_base& loi = le_mod_turb_th.loi_paroi();
1760 Paroi_TBLE_scal_VDF& loi_tble_T = ref_cast_non_const(Paroi_TBLE_scal_VDF,loi);
1761
1762 const Champ_Don_base& ch_beta_t = le_fluide.beta_t();
1763 const DoubleTab& tab_champ_beta_t = ch_beta_t.valeurs();
1764 if (sub_type(Champ_Uniforme,ch_beta_t))
1765 {
1766 beta_t = std::max(tab_champ_beta_t(0,0),DMINFLOAT);
1767 }
1768 else
1769 {
1770 Cerr << " On ne sait pas traiter un beta_t non uniforme dans Paroi TBLE_LRM !!!"<< finl;
1771 Process::exit();
1772 }
1773
1774
1775
1776 for(int i=0 ; i<nb_pts ; i++)
1777 {
1778 grad_T(i)=(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i+1)-loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i))/(eq_k_U_W[compteur_faces_paroi].get_y(i+1)-eq_k_U_W[compteur_faces_paroi].get_y(i));
1779 }
1780
1781 for(int i=1 ; i<nb_pts ; i++)
1782 {
1783 grad_T_moy(i)=0.5*(grad_T(i)+grad_T(i-1));
1784 }
1785
1786
1787 grad_T_moy(0)=(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,1)-loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,0))/(eq_k_U_W[compteur_faces_paroi].get_y(1)-eq_k_U_W[compteur_faces_paroi].get_y(0));
1788
1789 const Champ_Don_base& ch_gravite=le_fluide.gravite();
1790 const DoubleVect& gravite = ch_gravite.valeurs();
1791
1792 if (!sub_type(Champ_Uniforme,ch_gravite))
1793 {
1794 Cerr << " On ne sait pas traiter la gravite non uniforme dans Paroi TBLE_LRM !!!"<< finl;
1795 Process::exit();
1796 }
1797
1798
1799
1800 //calcul de la gravite
1801
1802 double norm_n=domaine_VDF.face_surfaces(num_face);
1803 double gn=0.;
1804
1805 for (int idim=0; idim<dimension; idim++)
1806 {
1807 gn+=gravite(idim)*domaine_VDF.face_normales(num_face,idim)/norm_n;
1808 }
1809
1810 DoubleVect gt_vect(dimension);
1811
1812 for (int idim=0; idim<dimension; idim++)
1813 {
1814 gt_vect(idim) = gravite(idim)-gn*domaine_VDF.face_normales(num_face,idim)/norm_n;
1815 }
1816
1817 double g_t=gt_vect(1-ori);
1818
1819
1820
1821 for(int i=0 ; i<nb_pts; i++)
1822 {
1823 source_T_k(i) = g_t*beta_t*((diffu.calculer_a_local(i)-d_visco)/(0.9))*grad_T_moy(i);
1824 source_T_U(i) = g_t*beta_t*(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i)-T0);
1825 }
1826
1827 }
1828 }
1829
1830 //terme source a la paroi
1831
1832 eps0=2*d_visco*eq_k_U_W[compteur_faces_paroi].get_Unp1(1,1)/(eq_k_U_W[compteur_faces_paroi].get_y(1)*eq_k_U_W[compteur_faces_paroi].get_y(1));
1833 eq_k_U_W[compteur_faces_paroi].set_F(1,0,-eps0 + source_T_k(0));
1834 eq_k_U_W[compteur_faces_paroi].set_F(0,0, -gradient_de_pression0 + ts0 - source_T_U(0));
1835
1836 for(int i=1 ; i<nb_pts; i++)
1837 {
1838
1839 //echelles de vitesse et de longueur (u ne sert pas actuellement, il est mis a 0)
1840 double k=eq_k_U_W[compteur_faces_paroi].get_Unp1(1,i);
1841 double y=eq_k_U_W[compteur_faces_paroi].get_y(i);
1842
1843 vv_bar=le_mod_ech.calculer_vv_bar(y, k, u, d_visco);
1844 l_eps=le_mod_ech.calculer_l_eps(y, k, u, d_visco);
1845
1846
1847 //terme source en i pour k
1848 eps=eq_k_U_W[compteur_faces_paroi].get_Unp1(1,i)*sqrt(vv_bar)/l_eps;
1849
1850 eq_k_U_W[compteur_faces_paroi].set_F(1,i,(diffu.calculer_a_local(i)-d_visco)
1851 *grad_vit_elem_moyx(i)*grad_vit_elem_moyx(i)-eps+source_T_k(i));
1852
1853 //terme source pour U
1854 eq_k_U_W[compteur_faces_paroi].set_F(0,i, -gradient_de_pression0 + ts0 - source_T_U(i) );
1855
1856 }
1857
1858 //On resoud les equations aux limites simplifiees
1859 //pendant un pas de temps du maillage grossier
1860 // On obtient U, V, W au temps n+1
1861
1862 eq_k_U_W[compteur_faces_paroi].aller_au_temps(tps);
1863
1864
1865 //Message lorsque k est negatif. Modification de la valeur a 1e-9
1866 for(int i=0 ; i<nb_pts ; i++)
1867 {
1868 if (eq_k_U_W[compteur_faces_paroi].get_Unp1(1,i)<0.)
1869 {
1870 Cerr << "Warning dans TBLE_LRM : k <0 !!! k = " << eq_k_U_W[compteur_faces_paroi].get_Unp1(1,i) << finl;
1871 eq_k_U_W[compteur_faces_paroi].set_Unp1(1,i,1.e-9);
1872 }
1873 }
1874
1875 if(itmax<eq_k_U_W[compteur_faces_paroi].get_it())
1876 {
1877 itmax = eq_k_U_W[compteur_faces_paroi].get_it();
1878 if(itmax>max_it)
1879 Cerr << "WARNING : TOO MANY ITERATIONS ARE NEEDED IN THE TBLE MESH"
1880 << finl;
1881 }
1882
1883
1884 ///////////////////////////////////////////////
1885 ///// Determination des deux composantes //////
1886 ///// du cisaillement a la paroi //////
1887 ///////////////////////////////////////////////
1888
1889 if (ori == 0)
1890 {
1891 //Cerr << "ori=0" << finl;
1892 Cisaillement_paroi_(num_face,0) = 0.;
1893 Cisaillement_paroi_(num_face,1) = eq_k_U_W[compteur_faces_paroi].get_cis(0)*signe;
1894 }
1895 else if (ori == 1)
1896 {
1897
1898 Cisaillement_paroi_(num_face,0) = eq_k_U_W[compteur_faces_paroi].get_cis(0)*signe;
1899 Cisaillement_paroi_(num_face,1) = 0.;
1900 }
1901
1902 tab_u_star_(num_face) =pow(eq_k_U_W[compteur_faces_paroi].get_cis(0)*eq_k_U_W[compteur_faces_paroi].get_cis(0),0.25);
1903 tab_d_plus_(num_face) = dist*tab_u_star(num_face)/d_visco;
1904
1905 compteur_faces_paroi++;
1906
1907
1908 }//fin boucle sur les faces de bords parietaux
1909
1910
1911 }//fin if(dim2)
1912
1913
1914 ///////////////////////
1915 // Si probleme 3D //
1916 ///////////////////////
1917
1918 else if(dimension == 3)
1919 {
1920 //eq_k_U_W.dimensionner(2*le_dom_dis_->nb_faces_bord());
1921
1922 //Boucle sur les faces des bords parietaux
1923
1924 for (int num_face=ndeb; num_face<nfin; num_face++)
1925 {
1926 //ici ori=direction perpendiculaire a la paroi
1927 ori = orientation(num_face);
1928 //Cerr << "ori = " << ori << finl;
1929
1930 //compteur++;
1931 //Selection de l'element associe a une face parietale
1932 //et initialisation de l'entier signe
1933
1934 if ((elem = face_voisins(num_face,0)) == -1)
1935 {
1936 elem = face_voisins(num_face,1);
1937 signe = 1;
1938 }
1939 else signe = -1 ;
1940
1941
1942 if(dt<dt_min)
1943 eq_k_U_W[compteur_faces_paroi].set_dt(1.e12);
1944 else
1945 eq_k_U_W[compteur_faces_paroi].set_dt(dt);
1946
1947
1948 // 1er couple de faces perpendiculaires a la paroi
1949 int face1 = elem_faces(elem,(ori+1));
1950 int face2 = elem_faces(elem,(ori+4)%6);
1951 // 2e couple de faces perpendiculaires a la paroi
1952 int face3 = elem_faces(elem,(ori+2));
1953 int face4 = elem_faces(elem,(ori+5)%6);
1954
1955
1956
1957
1958
1959
1960 ////////////////////////////////////////////////////////////
1961 ////// 1er couple de faces perpendiculaires a la paroi /////
1962 ////////////////////////////////////////////////////////////
1963
1964 vmoy = 0.5*(vit(face1) + vit(face2));
1965
1966 //vitesse sur le maillage fin a l'instant initial
1967
1968
1969 ts0 = 0.5*(termes_sources(face1)/volumes_entrelaces(face1)
1970 + termes_sources(face2)/volumes_entrelaces(face2));
1971
1972 eq_k_U_W[compteur_faces_paroi].set_u_y0(0,0.); //1ere composante de la vitesse a la paroi
1973
1974 eq_k_U_W[compteur_faces_paroi].set_u_yn(0,vmoy); //vitesse en yn
1975
1976 //*** Gradient de pression du 1er couple de faces
1977 // perpendiculaires a la paroi ***
1978
1979 //gradient de pression (constant dans le maillage fin)
1980 gradient_de_pression0 = 0.5*(grad_p(face1)+grad_p(face2));
1981
1982 //eq_k_U_W[compteur_faces_paroi].set_F(0, gradient_de_pression0 - ts0);
1983
1984
1985 //dpzbis += (gradient_de_pression_bis);
1986
1987 ////////////////////////////////////////////////////////////
1988 ////// 2e couple de faces perpendiculaires a la paroi /////
1989 ////////////////////////////////////////////////////////////
1990
1991 vmoy = 0.5*(vit(face3) + vit(face4));
1992
1993 ts1 = 0.5*(termes_sources(face3)/volumes_entrelaces(face3)
1994 + termes_sources(face4)/volumes_entrelaces(face4));
1995
1996 eq_k_U_W[compteur_faces_paroi].set_u_y0(1,0.); //2e composante de la vitesse a la paroi
1997
1998 eq_k_U_W[compteur_faces_paroi].set_u_yn(1,vmoy); //vitesse en yn
1999
2000 //*** Gradient de pression du 2e couple de faces
2001 // perpendiculaires a la paroi ***
2002
2003
2004 //gradient de pression (constant dans le maillage fin)
2005
2006 gradient_de_pression1 = 0.5*(grad_p(face3)+grad_p(face4));
2007
2008
2009 //eq_k_U_W[compteur_faces_paroi].set_F(1, gradient_de_pression1 - ts1);
2010 //Calcul du gradient On prend U au temps n
2011
2012 //equation en k ;
2013
2014
2015 double d_visco;
2016
2017 if (l_unif==1)
2018 d_visco = visco;
2019 else
2020 d_visco = tab_visco[elem];
2021
2022
2023 Diffu_totale_base& diffu = eq_k_U_W[compteur_faces_paroi].get_diffu();
2024
2025
2026 for(int i=0 ; i<nb_pts; i++)
2027 {
2028
2029 grad_vit_elemx(i)=(eq_k_U_W[compteur_faces_paroi].get_Unp1(0,i+1)-eq_k_U_W[compteur_faces_paroi].get_Unp1(0,i))/(eq_k_U_W[compteur_faces_paroi].get_y(i+1)-eq_k_U_W[compteur_faces_paroi].get_y(i));
2030 grad_vit_elemz(i)=(eq_k_U_W[compteur_faces_paroi].get_Unp1(1,i+1)-eq_k_U_W[compteur_faces_paroi].get_Unp1(1,i))/(eq_k_U_W[compteur_faces_paroi].get_y(i+1)-eq_k_U_W[compteur_faces_paroi].get_y(i));
2031
2032 }
2033
2034 //gradient en 0
2035 grad_vit_elem_moyx(0)=(eq_k_U_W[compteur_faces_paroi].get_Unp1(0,1)-eq_k_U_W[compteur_faces_paroi].get_Unp1(0,0))/(eq_k_U_W[compteur_faces_paroi].get_y(1)-eq_k_U_W[compteur_faces_paroi].get_y(0));
2036 grad_vit_elem_moyz(0)=(eq_k_U_W[compteur_faces_paroi].get_Unp1(1,1)-eq_k_U_W[compteur_faces_paroi].get_Unp1(1,0))/(eq_k_U_W[compteur_faces_paroi].get_y(1)-eq_k_U_W[compteur_faces_paroi].get_y(0));
2037
2038 //Cacul du gradient a l'ordonnee y(i)
2039 for(int i=1 ; i<nb_pts ; i++)
2040 {
2041
2042 grad_vit_elem_moyx(i)=0.5*(grad_vit_elemx(i)+grad_vit_elemx(i-1));
2043 grad_vit_elem_moyz(i)=0.5*(grad_vit_elemz(i)+grad_vit_elemz(i-1));
2044
2045 }
2046
2047 // On calcul vv_bar, l_eps, y_nu, y_star avec k au temps n
2048
2049 double dist;
2050
2051 //Distance a la paroi du 1er centre de maille
2052 if (axi)
2053 dist = domaine_VDF.dist_norm_bord_axi(num_face);
2054 else
2055 dist = domaine_VDF.dist_norm_bord(num_face);
2056
2057
2058 // Definition de la face servant a recuperer le k de la deuxieme maille
2059 int face5 = elem_faces(elem,(ori+3)%6);
2060 int face6 = elem_faces(elem,(ori));
2061
2062
2063
2064 int n2=face5;
2065 if (face5==num_face)
2066 {
2067 n2=face6;
2068 }
2069 int elem1=face_voisins(n2,0);
2070 if (elem1 == elem)
2071 elem1=face_voisins(n2,1);
2072
2073
2074 double dy = domaine_VDF.dist_norm(n2);
2075 // recuperartion de l'element correspondant a la deuxieme maille
2076
2077 double kcl;
2078 double a1=d_visco+0.5*(nu_t(elem)+nu_t(elem1));
2079 double a2=0.5*(diffu.calculer_a_local(nb_pts)+diffu.calculer_a_local(nb_pts-1));
2080
2081 double b1=a1/dy;
2082 double b2=a2/(eq_k_U_W[compteur_faces_paroi].get_y(nb_pts)-eq_k_U_W[compteur_faces_paroi].get_y(nb_pts-1));
2083
2084
2085
2086 kcl=(b1*tab_k_eps(elem1,0)+b2*eq_k_U_W[compteur_faces_paroi].get_Unp1(2,nb_pts-1))/(b1+b2);
2087
2088
2089 //initialisation de epsilon pour maillage 3D tel que nut(tble)=nut(3D)
2090
2091 if(tab_k_eps(elem,0)<0)
2092 {
2093 //Cerr << "Warning dans TBLE_LRM : k3D <0 !!! k = " << tab_k_eps(elem,0) << finl;
2094 tab_k_eps(elem,0)=1.e-8;
2095 tab_k_eps(elem,1)=0.09*tab_k_eps(elem,0)*tab_k_eps(elem,0)/(diffu.calculer_a_local(nb_pts)-d_visco);
2096 }
2097
2098 else
2099 {
2100 tab_k_eps(elem,0)=kcl;
2101 tab_k_eps(elem,1)=0.09*tab_k_eps(elem,0)*tab_k_eps(elem,0)/(diffu.calculer_a_local(nb_pts)-d_visco);
2102 }
2103
2104 double T0=0;
2105 double beta_t=0.;
2106
2107 source_T_k=0.;
2108 source_T_U=0.;
2109
2110 if (avec_boussi==1)
2111 {
2112 int boussi_ok=0;
2113 const Sources& les_sources=eqnNS.sources();
2114 int nb_sources=les_sources.size();
2115
2116 for (int j=0; j<nb_sources; j++)
2117 {
2118 const Source_base& ts = les_sources(j).valeur();
2119 if (sub_type(Terme_Boussinesq_base,ts))
2120 {
2121 const Terme_Boussinesq_base& terme_boussi = ref_cast(Terme_Boussinesq_base,ts);
2122 T0 = terme_boussi.Scalaire0(0);
2123 boussi_ok=1;
2124 }
2125 }
2126
2127 if (boussi_ok==1)
2128 {
2129 const Equation_base& eqn_th = mon_modele_turb_hyd->equation().probleme().equation(1);
2130 const Modele_turbulence_scal_base& le_mod_turb_th = ref_cast(Modele_turbulence_scal_base,eqn_th.get_modele(TURBULENCE).valeur());
2131 const Turbulence_paroi_scal_base& loi = le_mod_turb_th.loi_paroi();
2132 Paroi_TBLE_scal_VDF& loi_tble_T = ref_cast_non_const(Paroi_TBLE_scal_VDF,loi);
2133
2134 const Champ_Don_base& ch_beta_t = le_fluide.beta_t();
2135 const DoubleTab& tab_champ_beta_t = ch_beta_t.valeurs();
2136 if (sub_type(Champ_Uniforme,ch_beta_t))
2137 {
2138 beta_t = std::max(tab_champ_beta_t(0,0),DMINFLOAT);
2139 }
2140 else
2141 {
2142 Cerr << " On ne sait pas traiter un beta_t non uniforme dans Paroi TBLE_LRM !!!"<< finl;
2143 Process::exit();
2144 }
2145
2146
2147
2148 for(int i=0 ; i<nb_pts ; i++)
2149 {
2150 grad_T(i)=(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i+1)-loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i))/(eq_k_U_W[compteur_faces_paroi].get_y(i+1)-eq_k_U_W[compteur_faces_paroi].get_y(i));
2151 }
2152
2153 for(int i=1 ; i<nb_pts ; i++)
2154 {
2155 grad_T_moy(i)=0.5*(grad_T(i)+grad_T(i-1));
2156 }
2157
2158
2159 grad_T_moy(0)=(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,1)-loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,0))/(eq_k_U_W[compteur_faces_paroi].get_y(1)-eq_k_U_W[compteur_faces_paroi].get_y(0));
2160
2161 const Champ_Don_base& ch_gravite=le_fluide.gravite();
2162 const DoubleVect& gravite = ch_gravite.valeurs();
2163
2164 if (!sub_type(Champ_Uniforme,ch_gravite))
2165 {
2166 Cerr << " On ne sait pas traiter la gravite non uniforme dans Paroi TBLE_LRM !!!"<< finl;
2167 Process::exit();
2168 }
2169
2170
2171
2172 //calcul de la gravite
2173
2174
2175 double norm_n=domaine_VDF.face_surfaces(num_face);
2176 double gn=0.;
2177
2178 for (int idim=0; idim<dimension; idim++)
2179 {
2180 gn+=gravite(idim)*domaine_VDF.face_normales(num_face,idim)/norm_n;
2181 }
2182
2183 DoubleVect gt_vect(dimension);
2184
2185 for (int idim=0; idim<dimension; idim++)
2186 {
2187 gt_vect(idim) = gravite(idim)-gn*domaine_VDF.face_normales(num_face,idim)/norm_n;
2188 }
2189
2190 double g_t=gt_vect(1-ori);
2191
2192 for(int i=1 ; i<nb_pts; i++)
2193 {
2194 source_T_k(i) =g_t*beta_t*((diffu.calculer_a_local(i)-d_visco)/(0.9))*grad_T_moy(i);
2195 source_T_U(i) =g_t*beta_t*(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i)-T0);
2196
2197 }
2198
2199 }
2200 }
2201
2202 eps0=2*d_visco*eq_k_U_W[compteur_faces_paroi].get_Unp1(2,1)/(eq_k_U_W[compteur_faces_paroi].get_y(1)*eq_k_U_W[compteur_faces_paroi].get_y(1));
2203 eq_k_U_W[compteur_faces_paroi].set_F(2,0,-eps0 + source_T_k(0));
2204 eq_k_U_W[compteur_faces_paroi].set_F(0,0, -gradient_de_pression0 + ts0 - source_T_U(0));
2205
2206 eq_k_U_W[compteur_faces_paroi].set_F(1,0, -gradient_de_pression1 + ts1 - source_T_U(0) );
2207
2208
2209 for(int i=1 ; i<nb_pts ; i++)
2210 {
2211
2212 //echelles de vitesse et de longueur (u ne sert pas actuellement, il est mis a 0)
2213 double k=eq_k_U_W[compteur_faces_paroi].get_Unp1(2,i);
2214 double y=eq_k_U_W[compteur_faces_paroi].get_y(i);
2215
2216
2217 vv_bar=le_mod_ech.calculer_vv_bar(y, k, u, d_visco);
2218 l_eps=le_mod_ech.calculer_l_eps(y, k, u, d_visco);
2219
2220 //terme source a la paroi
2221
2222
2223 eps=eq_k_U_W[compteur_faces_paroi].get_Unp1(2,i)*sqrt(vv_bar)/l_eps;
2224
2225 eq_k_U_W[compteur_faces_paroi].set_F(2,i,(diffu.calculer_a_local(i)-d_visco)*(grad_vit_elem_moyx(i)*grad_vit_elem_moyx(i)+grad_vit_elem_moyz(i)*grad_vit_elem_moyz(i))-eps
2226 + source_T_k(i));
2227
2228 // pour U et W
2229
2230 eq_k_U_W[compteur_faces_paroi].set_F(0,i, -gradient_de_pression0 + ts0 - source_T_U(i));
2231 eq_k_U_W[compteur_faces_paroi].set_F(1,i, -gradient_de_pression1 + ts1 - source_T_U(i) );
2232
2233 }
2234
2235
2236
2237 //On resoud les equations aux limites simplifiees
2238 //pendant un pas de temps du maillage grossier
2239
2240 eq_k_U_W[compteur_faces_paroi].set_u_y0(2,0.); //k a la paroi
2241 eq_k_U_W[compteur_faces_paroi].set_u_yn(2,kcl);//k a yn
2242
2243 eq_k_U_W[compteur_faces_paroi].aller_au_temps(tps);
2244
2245 //Message lorsque k est negatif. Modification de la valeur a 1e-9
2246 for(int i=0 ; i<nb_pts ; i++)
2247 {
2248 if (eq_k_U_W[compteur_faces_paroi].get_Unp1(2,i)<0.)
2249 {
2250 // Cerr << "Warning dans TBLE_LRM : k <0 !!! k = " << eq_k_U_W[compteur_faces_paroi].get_Unp1(2,i) << finl;
2251 eq_k_U_W[compteur_faces_paroi].set_Unp1(2,i,1.e-9);
2252 }
2253 }
2254
2255 if(itmax<eq_k_U_W[compteur_faces_paroi].get_it())
2256 {
2257 itmax = eq_k_U_W[compteur_faces_paroi].get_it();
2258 if(itmax>max_it)
2259 Cerr << "WARNING : TOO MANY ITERATIONS ARE NEEDED IN THE TBLE MESH"
2260 << finl;
2261 }
2262
2263
2264 ///////////////////////////////////////////////
2265 ///// Determination des deux composantes //////
2266 ///// du cisaillement a la paroi //////
2267 ///////////////////////////////////////////////
2268
2269 if (ori == 0)
2270 {
2271 //Cerr << "ori=0" << finl;
2272 Cisaillement_paroi_(num_face,0) = 0.;
2273 Cisaillement_paroi_(num_face,1) = eq_k_U_W[compteur_faces_paroi].get_cis(0)*signe;
2274 Cisaillement_paroi_(num_face,2) = eq_k_U_W[compteur_faces_paroi].get_cis(1)*signe;
2275 }
2276 else if (ori == 1)
2277 {
2278
2279 Cisaillement_paroi_(num_face,0) = eq_k_U_W[compteur_faces_paroi].get_cis(1)*signe;
2280 Cisaillement_paroi_(num_face,1) = 0.;
2281 Cisaillement_paroi_(num_face,2) = eq_k_U_W[compteur_faces_paroi].get_cis(0)*signe;
2282
2283 }
2284 else
2285 {
2286 //Cerr << "ori=2" << finl;
2287 Cisaillement_paroi_(num_face,0) = eq_k_U_W[compteur_faces_paroi].get_cis(0)*signe;
2288 Cisaillement_paroi_(num_face,1) = eq_k_U_W[compteur_faces_paroi].get_cis(1)*signe;
2289 Cisaillement_paroi_(num_face,2) = 0.;
2290 }
2291
2292 tab_u_star_(num_face) = pow((eq_k_U_W[compteur_faces_paroi].get_cis(0)*eq_k_U_W[compteur_faces_paroi].get_cis(0)
2293 +eq_k_U_W[compteur_faces_paroi].get_cis(1)*eq_k_U_W[compteur_faces_paroi].get_cis(1)),0.25);
2294 tab_d_plus_(num_face) = dist*tab_u_star(num_face)/d_visco;
2295
2296
2297 compteur_faces_paroi++;
2298
2299
2300 for(int j=0; j<nb_post_pts; j++)
2301 {
2302 Nom tmp2;
2303 tmp2="nut_post_";
2304 tmp2+=nom_pts[j];
2305 tmp2+=".dat";
2306
2307 SFichier fic_nut(tmp2); //creation fichier de post-traitement des points TBLE
2308
2309 Diffu_totale_base& diffu1 = eq_k_U_W[num_faces_post(j)].get_diffu();
2310 for(int i=0 ; i<nb_pts+1 ; i++)
2311 {
2312 fic_nut << eq_k_U_W[num_faces_post(j)].get_y(i)<< " " << diffu1.calculer_a_local(i)-d_visco << finl;
2313 }
2314 fic_nut << finl;
2315
2316 }
2317
2318 }//fin boucle sur les faces de bords parietaux
2319
2320
2321
2322 }//fin if(dim=3)
2323 }//fin if(sub_type(dirichlet))
2324
2325 }//fin boucle sur les bords
2326
2327 if (Process::je_suis_maitre())
2328 {
2329 SFichier fic("iter.dat",ios::app); // ouverture du fichier iter.dat
2330 fic << tps << " " << itmax << finl;
2331 }
2332
2333
2334 Cisaillement_paroi_.echange_espace_virtuel();
2335
2336
2337 return 1;
2338}
2339
2340void ParoiVDF_TBLE_LRM::imprimer_ustar(Sortie& os) const
2341{
2342 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
2343 const double tps = eqn_hydr.inconnue().temps();
2344 const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_dis_.valeur());
2345 const IntVect& orientation = domaine_VDF.orientation();
2346
2347
2348 // const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_dis_.valeur());
2349 // int ndeb,nfin;
2350 //double upmoy,dpmoy,utaumoy;
2351 //upmoy=0.;dpmoy=0.;utaumoy=0.;
2352 //int compt=0;
2353
2354 for(int j=0; j<nb_post_pts; j++)
2355 {
2356 Nom tmp;
2357 tmp="tble_post_";
2358 tmp+=nom_pts[j];
2359 tmp+=".dat";
2360
2361 SFichier fic_post(tmp /*,ios::app*/);
2362
2363 int ori=orientation(num_faces_post2(j));
2364
2365 if (dimension == 2)
2366 {
2367 double x=domaine_VDF.xv(num_faces_post2(j),0);
2368 double y=domaine_VDF.xv(num_faces_post2(j),1);
2369 fic_post << "#" << "x= " << x << " " << "y= " << y << finl;
2370 }
2371
2372 if(dimension == 3)
2373 {
2374 double x=domaine_VDF.xv(num_faces_post2(j),0);
2375 double y=domaine_VDF.xv(num_faces_post2(j),1);
2376 double z=domaine_VDF.xv(num_faces_post2(j),2);
2377 fic_post << "#" << "x= " << x << " " << "y= " << y << " " << "z= " << z << finl;
2378 }
2379
2380
2381 fic_post << "#"<< tps << " " << "u*= " << tab_u_star(num_faces_post2(j)) << finl;
2382
2383
2384 if(dimension == 2)
2385 {
2386
2387 for(int i=0; i<nb_pts+1; i++)
2388 {
2389 fic_post << eq_k_U_W[num_faces_post(j)].get_y(i) << " " << eq_k_U_W[num_faces_post(j)].get_Unp1(0,i) <<" " <<
2390 eq_k_U_W[num_faces_post(j)].get_Unp1(1,i) << finl;
2391 }
2392 }
2393
2394
2395 if(dimension == 3)
2396 {
2397 if(ori == 0)
2398 {
2399 for(int i=0; i<nb_pts+1; i++)
2400 {
2401 fic_post << eq_k_U_W[num_faces_post(j)].get_y(i) << " " << eq_k_U_W[num_faces_post(j)].get_Unp1(0,i) <<
2402 " " << eq_k_U_W[num_faces_post(j)].get_Unp1(1,i) << " " << eq_k_U_W[num_faces_post(j)].get_Unp1(2,i) << finl;
2403 }
2404 }
2405 else
2406 {
2407 for(int i=0; i<nb_pts+1; i++)
2408 {
2409 fic_post << eq_k_U_W[num_faces_post(j)].get_y(i) << " " <<
2410 eq_k_U_W[num_faces_post(j)].get_Unp1(1,i) <<
2411 " " << eq_k_U_W[num_faces_post(j)].get_Unp1(0,i) << " " << eq_k_U_W[num_faces_post(j)].get_Unp1(2,i) << finl;
2412 }
2413 }
2414 }
2415
2416
2417 }
2418
2419}
2420
2421
2422
2423
2424
2425
2426
2427
2428
2429
2430
2431
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
Champ_base::temps
double temps() const
Renvoie le temps du champ.
Definition Champ_base.cpp:1004
Cond_lim
classe Cond_lim Classe generique servant a representer n'importe quelle classe
Definition Cond_lim.h:31
Diffu_k
Classe Diffu_k Calsse derivant de la classe Diffu_totale_hyd__base et specifiant la valeur.
Definition Diffu_k.h:33
Diffu_k::associer_modele_echelles
void associer_modele_echelles(Mod_echelle_LRM_base &mod)
Definition Diffu_k.h:39
Diffu_totale_base
Classe Diffu_totale_base Classe abstraite calculant la diffusivite totale (somme diffusivite.
Definition Diffu_totale_base.h:33
Diffu_totale_base::calculer_a_local
virtual double calculer_a_local(int ind)=0
Dirichlet_paroi_fixe
classe Dirichlet_paroi_fixe Represente une paroi immobile dans une equation de type Navier_Stokes.
Definition Dirichlet_paroi_fixe.h:26
Domaine_VDF
class Domaine_VDF
Definition Domaine_VDF.h:64
Domaine_VDF::dist_norm
double dist_norm(int num_face) const override
Definition Domaine_VDF.h:342
Domaine_VDF::dist_norm_bord_axi
double dist_norm_bord_axi(int num_face) const
Definition Domaine_VDF.h:400
Domaine_VDF::orientation
int orientation(int) const override
inline DoubleVect& Domaine_VDF::porosite_face() {
Definition Domaine_VDF.h:297
Domaine_VDF::face_normales
double face_normales(int, int) const override
Definition Domaine_VDF.h:187
Domaine_VDF::dist_norm_bord
double dist_norm_bord(int num_face) const override
Definition Domaine_VDF.h:382
Domaine_VF::face_surfaces
virtual const DoubleVect & face_surfaces() const
Definition Domaine_VF.h:51
Domaine_VF::nb_faces
int nb_faces() const
renvoie le nombre global de faces.
Definition Domaine_VF.h:471
Domaine_VF::volumes_entrelaces
DoubleVect & volumes_entrelaces()
Definition Domaine_VF.h:99
Domaine_VF::xv
double xv(int num_face, int k) const
Definition Domaine_VF.h:76
Domaine_VF::creer_tableau_faces_bord
void creer_tableau_faces_bord(Array_base &, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT) const
Definition Domaine_VF.cpp:892
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
Entree
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
Eq_couch_lim::get_Unp1
double get_Unp1(int j, int i)
Definition Eq_couch_lim.h:185
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::get_modele
virtual const RefObjU & get_modele(Type_modele type) const
Definition Equation_base.cpp:1894
Equation_base::sources
Sources & sources()
Renvoie les termes sources asssocies a l'equation.
Definition Equation_base.cpp:348
Equation_base::inconnue
virtual const Champ_Inc_base & inconnue() const =0
Equation_base::probleme
Probleme_base & probleme()
Renvoie le probleme associe a l'equation.
Definition Equation_base.cpp:747
Equation_base::schema_temps
Schema_Temps_base & schema_temps()
Renvoie le schema en temps associe a l'equation.
Definition Equation_base.cpp:865
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
Milieu_base::beta_t
virtual const Champ_Don_base & beta_t() const
Renvoie beta_t du milieu.
Definition Milieu_base.cpp:887
Milieu_base::gravite
virtual const Champ_Don_base & gravite() const
Renvoie la gravite du milieu si elle a ete associe provoque une erreur sinon.
Definition Milieu_base.cpp:716
Mod_echelle_LRM_base
Classe Mod_echelle_LRM_base.
Definition Mod_echelle_LRM_base.h:31
Mod_echelle_LRM_base::calculer_vv_bar
virtual double calculer_vv_bar(double y, double k, double u, double visco_cin)=0
Mod_echelle_LRM_base::calculer_l_eps
virtual double calculer_l_eps(double y, double k, double u, double visco_cin)=0
Modele_turbulence_scal_base
Classe Modele_turbulence_scal_base Cette classe represente un modele de turbulence pour une equation ...
Definition Modele_turbulence_scal_base.h:40
Modele_turbulence_scal_base::loi_paroi
const Turbulence_paroi_scal_base & loi_paroi() const
Renvoie la loi de turbulence sur la paroi (version const).
Definition Modele_turbulence_scal_base.h:90
Navier_Stokes_std
classe Navier_Stokes_std Cette classe porte les termes de l'equation de la dynamique
Definition Navier_Stokes_std.h:56
Navier_Stokes_std::operateur_gradient
Operateur_Grad & operateur_gradient()
Renvoie l'operateur de calcul du gradient associe a l'equation.
Definition Navier_Stokes_std.cpp:568
Navier_Stokes_std::pression
Champ_Inc_base & pression()
Definition Navier_Stokes_std.h:123
Nom
class Nom Une chaine de caractere pour nommer les objets de TRUST
Definition Nom.h:31
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::axi
static int axi
Definition Objet_U.h:101
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
Operateur_Grad
Classe Operateur_Grad Classe generique de la hierarchie des operateurs calculant le gradient.
Definition Operateur_Grad.h:31
ParoiVDF_TBLE_LRM
CLASS: ParoiVDF_TBLE_LRM.
Definition ParoiVDF_TBLE_LRM.h:39
ParoiVDF_TBLE_LRM::max_it
int max_it
Definition ParoiVDF_TBLE_LRM.h:57
ParoiVDF_TBLE_LRM::modele_visco
Motcle modele_visco
Definition ParoiVDF_TBLE_LRM.h:59
ParoiVDF_TBLE_LRM::init_lois_paroi
int init_lois_paroi() override
Definition ParoiVDF_TBLE_LRM.cpp:189
ParoiVDF_TBLE_LRM::nb_post_pts
int nb_post_pts
Definition ParoiVDF_TBLE_LRM.h:64
ParoiVDF_TBLE_LRM::avec_boussi
int avec_boussi
Definition ParoiVDF_TBLE_LRM.h:62
ParoiVDF_TBLE_LRM::epsilon
double epsilon
Definition ParoiVDF_TBLE_LRM.h:58
ParoiVDF_TBLE_LRM::num_faces_post2
IntTab num_faces_post2
Definition ParoiVDF_TBLE_LRM.h:68
ParoiVDF_TBLE_LRM::fac
double fac
Definition ParoiVDF_TBLE_LRM.h:58
ParoiVDF_TBLE_LRM::calculer_hyd
int calculer_hyd(DoubleTab &) override
Definition ParoiVDF_TBLE_LRM.cpp:1409
ParoiVDF_TBLE_LRM::nb_pts
int nb_pts
Definition ParoiVDF_TBLE_LRM.h:57
ParoiVDF_TBLE_LRM::imprimer_ustar
void imprimer_ustar(Sortie &os) const override
Definition ParoiVDF_TBLE_LRM.cpp:2340
ParoiVDF_TBLE_LRM::nom_pts
Noms nom_pts
Definition ParoiVDF_TBLE_LRM.h:66
ParoiVDF_TBLE_LRM::sonde_tble
DoubleTab sonde_tble
Definition ParoiVDF_TBLE_LRM.h:65
ParoiVDF_TBLE_LRM::k_init
double k_init
Definition ParoiVDF_TBLE_LRM.h:63
ParoiVDF_TBLE_LRM::nb_comp
int nb_comp
Definition ParoiVDF_TBLE_LRM.h:57
ParoiVDF_TBLE_LRM::calculer_hyd_BiK
int calculer_hyd_BiK(DoubleTab &, DoubleTab &) override
Definition ParoiVDF_TBLE_LRM.cpp:488
ParoiVDF_TBLE_LRM::num_faces_post
IntTab num_faces_post
Definition ParoiVDF_TBLE_LRM.h:67
Paroi_TBLE_QDM_Scal::get_eq_couche_lim_T
Eq_couch_lim & get_eq_couche_lim_T(int i)
Definition Paroi_TBLE_QDM_Scal.h:84
Paroi_TBLE_scal_VDF
CLASS: Paroi_ODVM_scal_VDF.
Definition Paroi_TBLE_scal_VDF.h:39
Paroi_hyd_base_VDF
Definition Paroi_hyd_base_VDF.h:29
Probleme_base::equation
virtual const Equation_base & equation(int) const =0
Process::nproc
static int nproc()
renvoie le nombre de processeurs dans le groupe courant Voir Comm_Group::nproc() et PE_Groups::curren...
Definition Process.cpp:104
Process::exit
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
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
SFichier
Cette classe est a la classe C++ ofstream ce que la classe Sortie est a la classe C++ ostream Elle re...
Definition SFichier.h:27
Schema_Temps_base::temps_courant
double temps_courant() const
Renvoie le temps courant.
Definition Schema_Temps_base.h:445
Schema_Temps_base::pas_temps_min
double pas_temps_min() const
Renvoie le pas de temps minimum.
Definition Schema_Temps_base.h:407
Schema_Temps_base::pas_de_temps
double pas_de_temps() const
Renvoie le pas de temps (delta_t) courant.
Definition Schema_Temps_base.h:393
Sortie
Classe de base des flux de sortie.
Definition Sortie.h:52
Source_base
classe Source_base Un objet Source_base est un terme apparaissant au second membre d'une
Definition Source_base.h:42
Sources
class Sources Sources represente une liste de Source.
Definition Sources.h:31
Sources::calculer
DoubleTab & calculer(DoubleTab &) const
Calcule la contribution de toutes les sources de la liste stocke le resultat dans le tableau passe en...
Definition Sources.cpp:98
TRUSTTab::resize
void resize(_SIZE_ n, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
Definition TRUSTTab.tpp:469
TRUSTVect::local_min_vect
_TYPE_ local_min_vect(Mp_vect_options opt=VECT_REAL_ITEMS) const
Definition TRUSTVect.h:155
TRUSTVect::resize
void resize(_SIZE_, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
Definition TRUSTVect.tpp:91
TRUST_Ref_Objet_U::valeur
const Objet_U & valeur() const
Definition TRUST_Ref.h:134
Terme_Boussinesq_base
Classe Terme_Boussinesq_base Cette classe represente le terme de gravite qui figure dans l'equation.
Definition Terme_Boussinesq_base.h:55
Terme_Boussinesq_base::Scalaire0
double Scalaire0(int i) const
Definition Terme_Boussinesq_base.h:63
Turbulence_paroi_base::tab_d_plus_
DoubleVect tab_d_plus_
Definition Turbulence_paroi_base.h:88
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::Cisaillement_paroi_
DoubleTab Cisaillement_paroi_
Definition Turbulence_paroi_base.h:85
Turbulence_paroi_scal_base
Classe Turbulence_paroi_scal_base Classe de base pour la hierarchie des classes representant les mode...
Definition Turbulence_paroi_scal_base.h:42