TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
ParoiVDF_TBLE.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.h>
18#include <Domaine_Cl_VDF.h>
19#include <Dirichlet_paroi_fixe.h>
20#include <Dirichlet_paroi_defilante.h>
21#include <Champ_Uniforme.h>
22#include <Fluide_base.h>
23#include <EFichier.h>
24#include <Diffu_lm.h>
25#include <Schema_Temps_base.h>
26#include <Operateur_Grad.h>
27#include <Probleme_base.h>
28#include <Modele_turbulence_scal_base.h>
29#include <Paroi_TBLE_scal_VDF.h>
30#include <Terme_Boussinesq_VDF_Face.h>
31#include <SFichier.h>
32#include <EcrFicPartage.h>
33#include <Navier_Stokes_std.h>
34#include <Modele_turbulence_hyd_base.h>
35#include <Param.h>
36
37Implemente_instanciable(ParoiVDF_TBLE, "Paroi_TBLE_VDF", Paroi_hyd_base_VDF);
38// XD Paroi_TBLE_VDF turbulence_paroi_base Paroi_TBLE_VDF BRACE Wall function using the thin boundary layer formulation.
39
41{
42 return os << que_suis_je() << " " << le_nom();
43}
44
46{
47 Cerr << "Reading of data for a " << que_suis_je() << " wall law" << finl;
48 Param param(que_suis_je());
49 set_param(param);
50 param.lire_avec_accolades_depuis(is);
51 return is;
52}
53
55{
59 param.ajouter_flag("alpha_cv",&alpha_cv); // XD_ADD_P rien
60 // XD_CONT Activate if TBLE and convective terms (not sure) (default = false)
61 param.ajouter_non_std("nu_t_dyn", (this));
62 param.ajouter_non_std("tps_start_stat_nu_t_dyn", (this));
63 param.ajouter_non_std("tps_nu_t_dyn", (this));
64 param.ajouter_non_std("sonde_tble", (this));
65 param.ajouter_non_std("stationnaire", (this));
66 param.ajouter_non_std("mu", (this));
67 param.ajouter_non_std("lambda", (this));
68 param.ajouter_non_std("sans_source_boussinesq", (this));
69}
70
75
76/////////////////////////////////////////////////////////////////////
77//
78// Implementation des fonctions de la classe ParoiVDF_TBLE
79//
80/////////////////////////////////////////////////////////////////////
81
82// /////////////////////////////////////////////////
83// // Initialisation des tableaux
84// ////////////////////////////////////////////////
85
87{
88 Cerr << "debut init_lois_paroi" << finl;
89
90 const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_dis_.valeur());
91 const IntVect& orientation = domaine_VDF.orientation();
92 const IntTab& face_voisins = domaine_VDF.face_voisins();
93 const IntTab& elem_faces = domaine_VDF.elem_faces();
94 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
95 const DoubleVect& vit = eqn_hydr.inconnue().valeurs();
96 const Fluide_base& le_fluide = ref_cast(Fluide_base,eqn_hydr.milieu());
97 int compteur_faces_paroi = 0;
98 int elem {0};
99 double vmoy = 0.;
100 double dist {0.0}; //distance du premier centre de maille a la paroi
101
102
104 Paroi_TBLE_QDM::init_lois_paroi(domaine_VDF, le_dom_Cl_dis_.valeur());
105
106 corresp.resize(le_dom_dis_->nb_faces_bord()); // Ce tableau ete introduit a l origine pour les termes convectifs
107 // On le garde ?
108 compteur_faces_paroi = 0; //Reinitialisation de compteur_faces_paroi
109
110 for (int n_bord = 0; n_bord < domaine_VDF.nb_front_Cl(); n_bord++)
111 {
112 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
113
114 if (sub_type(Dirichlet_paroi_fixe, la_cl.valeur()))
115 {
116 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
117 const int ndeb = le_bord.num_premiere_face();
118 const int nfin = ndeb + le_bord.nb_faces();
119
120 for (int num_face = ndeb; num_face < nfin; num_face++)
121 {
122 corresp[num_face]=compteur_faces_paroi;
123 compteur_faces_paroi++;
124 }
125 }
126 }
127
128 compteur_faces_paroi = 0; //Reinitialisation de compteur_faces_paroi
129
130 for (int n_bord = 0; n_bord < domaine_VDF.nb_front_Cl(); n_bord++)
131 {
132 // pour chaque condition limite on regarde son type
133 // On applique les lois de paroi uniquement
134 // aux voisinages des parois
135
136 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
137
138 if (sub_type(Dirichlet_paroi_fixe,la_cl.valeur()))
139 {
140 const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
141 const int ndeb = le_bord.num_premiere_face();
142 const int nfin = ndeb + le_bord.nb_faces();
143
144 if (dimension == 2)
145 {
146 // Boucle sur les faces des bords parietaux
147 for (int num_face = ndeb; num_face < nfin; num_face++)
148 {
149 //ici ori=direction perpendiculaire a la paroi
150 int ori = orientation(num_face);
151
152 if ((elem = face_voisins(num_face, 0)) == -1)
153 elem = face_voisins(num_face, 1);
154
155 // Distance a la paroi du 1er centre de maille
156 if (axi)
157 dist = domaine_VDF.dist_norm_bord_axi(num_face);
158 else
159 dist = domaine_VDF.dist_norm_bord(num_face);
160
161 Diffu_totale_hyd_base& diffu_hyd = ref_cast_non_const(Diffu_totale_hyd_base, eq_vit[corresp[num_face]].get_diffu()); //modele de viscosite turbulente
162 diffu_hyd.setKappa(Kappa_);
163
164 eq_vit[corresp[num_face]].set_y0(0.); //ordonnee de la paroi
165 eq_vit[corresp[num_face]].set_yn(dist); //ordonnee du 1er centre de maille
166 eq_vit[corresp[num_face]].initialiser(nb_comp, nb_pts,
167 fac, epsilon, max_it, nu_t_dyn); //nbre de pts maillage fin
168 eq_vit[corresp[num_face]].set_u_y0(0,0.); //1ere composante de la vitesse a la paroi
169
170 vmoy = 0.5*(vit(elem_faces(elem, (ori + 1) % 4)) + vit(elem_faces(elem, (ori + 3) % 4)));
171
172 eq_vit[corresp[num_face]].set_u_yn(0, vmoy); //1ere composante de la vitesse en yn
173
174 //vitesse sur le maillage fin a l'instant initial
175 if (reprise_ok == 0)
176 eq_vit[corresp[num_face]].set_Uinit_lin(0, 0., vmoy);
177 else
178 for (int itble = 0; itble < nb_pts + 1; itble++)
179 eq_vit[corresp[num_face]].set_Uinit(0, itble, valeurs_reprises(corresp[num_face], 0, itble));
180
181 compteur_faces_paroi++;
182 }//Fin boucle sur les faces de bord
183 }
184 else if (dimension == 3)
185 {
186 // Boucle sur les faces des bords parietaux
187 for (int num_face = ndeb; num_face < nfin; num_face++)
188 {
189 //ici ori=direction perpendiculaire a la paroi
190 const int ori = orientation(num_face);
191
192 if ((elem = face_voisins(num_face, 0)) == -1)
193 elem = face_voisins(num_face, 1);
194
195 // Distance a la paroi du 1er centre de maille
196 if (axi)
197 dist = domaine_VDF.dist_norm_bord_axi(num_face);
198 else
199 dist = domaine_VDF.dist_norm_bord(num_face);
200
201 eq_vit[corresp[num_face]].set_y0(0.); //ordonnee de la paroi
202 eq_vit[corresp[num_face]].set_yn(dist); //ordonnee du 1er centre de maille
203 eq_vit[corresp[num_face]].initialiser(nb_comp, nb_pts,
205 nu_t_dyn);
206 eq_vit[corresp[num_face]].set_u_y0(0, 0.); //1ere composante de la vitesse a la paroi
207 eq_vit[corresp[num_face]].set_u_y0(1, 0.); //2e composante de la vitesse a la paroi
208
209 vmoy = 0.5*(vit(elem_faces(elem, (ori + 1) % 6)) + vit(elem_faces(elem, (ori + 4) % 6)));
210
211 eq_vit[corresp[num_face]].set_u_yn(0, vmoy); //1ere composante de la vitesse en yn
212
213 // vitesse sur le maillage fin a l'instant initial
214 if (reprise_ok == 0)
215 eq_vit[corresp[num_face]].set_Uinit_lin(0, 0., vmoy);
216 else
217 for (int itble = 0; itble < nb_pts + 1; itble++)
218 eq_vit[corresp[num_face]].set_Uinit(0, itble, valeurs_reprises(corresp[num_face], 0, itble));
219
220 vmoy = 0.5*(vit(elem_faces(elem, (ori + 2) % 6)) + vit(elem_faces(elem, (ori + 5) % 6)));
221 eq_vit[corresp[num_face]].set_u_yn(1, vmoy); //2e composante de la vitesse en yn
222 if (reprise_ok == 0)
223 eq_vit[corresp[num_face]].set_Uinit_lin(1, 0., vmoy);
224 else
225 for (int itble = 0; itble < nb_pts + 1; itble++)
226 eq_vit[corresp[num_face]].set_Uinit(1, itble, valeurs_reprises(corresp[num_face], 1, itble));
227
228 compteur_faces_paroi++;
229
230 }//Fin boucle sur les faces de bord
231 }// Fin if dim 3
232 }//Fin if Paroi_fixe
233 }//Fin boucle sur les bords parietaux
234
235 // Si on n'a pas interdit la presence du terme de Boussi dans TBLE (mot cle "sans_source_boussinesq",
236 // Alors on verifie la presence ou non d'un terme source de Boussinesq dans le pb hydraulique :
237 if (source_boussinesq==1)
238 {
240 const Sources& les_sources=eqn_hydr.sources();
241 int nb_sources=les_sources.size();
242
243 for (int j=0; j<nb_sources; j++)
244 {
245 const Source_base& ts = les_sources(j).valeur();
246 if (sub_type(Terme_Boussinesq_VDF_Face,ts))
247 {
248 const Terme_Boussinesq_base& terme_boussi = ref_cast(Terme_Boussinesq_base,ts);
249 T0 = terme_boussi.Scalaire0(0);
251 const Champ_Don_base& ch_beta_t = le_fluide.beta_t();
252 const DoubleTab& tab_champ_beta_t = ch_beta_t.valeurs();
253 if (sub_type(Champ_Uniforme,ch_beta_t))
254 {
255 beta_t = std::max(tab_champ_beta_t(0,0),DMINFLOAT);
256 }
257 else
258 {
259 Cerr << " On ne sait pas traiter un beta_t non uniforme dans TBLE !!!"<< finl;
261 }
262
263 break;
264 }
265 }
266 }
267
268
269
270 Cerr << "fin init_lois_paroi" << finl;
271 return 1;
272}
273
274
275/////////////////////////////////////////////////
276//Calcul du cisaillement a la paroi
277////////////////////////////////////////////////
278
279int ParoiVDF_TBLE::calculer_hyd(DoubleTab& nu_t, DoubleTab& tab_k)
280{
281 return calculer_hyd(nu_t, 0, tab_k);
282}
283
284int ParoiVDF_TBLE::calculer_hyd(DoubleTab& tab_k_eps)
285{
286
287 DoubleTab bidon(0);
288 // bidon ne servira pas
289 return calculer_hyd(tab_k_eps, 1, bidon);
290
291}
292
293int ParoiVDF_TBLE::calculer_hyd_BiK(DoubleTab& tab_k, DoubleTab& tab_eps)
294{
295// on est dans le cas k-eps
296
297 const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_dis_.valeur());
298 const IntVect& orientation = domaine_VDF.orientation();
299 const IntTab& face_voisins = domaine_VDF.face_voisins();
300 const IntTab& elem_faces = domaine_VDF.elem_faces();
301 const DoubleVect& volumes_entrelaces = domaine_VDF.volumes_entrelaces();
302 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
303
304 const Fluide_base& le_fluide = ref_cast(Fluide_base,eqn_hydr.milieu());
305 const Champ_Don_base& ch_visco_cin = le_fluide.viscosite_cinematique();
306 const DoubleTab& tab_visco = ch_visco_cin.valeurs();
307
308 if (source_boussinesq==1)
309 {
310 const Champ_Don_base& ch_gravite=le_fluide.gravite();
311
312 if (!sub_type(Champ_Uniforme,ch_gravite))
313 {
314 Cerr << " On ne sait pas traiter la gravite non uniforme dans TBLE !!!"<< finl;
316 }
317 }
318
319 double visco=-1;
320 int l_unif;
321 if (sub_type(Champ_Uniforme,ch_visco_cin))
322 {
323 visco = std::max(tab_visco(0,0),DMINFLOAT);
324 l_unif = 1;
325 }
326 else
327 l_unif = 0;
328
329 // preparer_calcul_hyd(tab);
330 if ((!l_unif) && (tab_visco.local_min_vect()<DMINFLOAT))
331 // on ne doit pas changer tab_visco ici !
332 {
333 Cerr << "In ParoiVDF_TBLE::calculer_hyd_BiK : visco = " << tab_visco.local_min_vect() << " <= 0 ? " << finl;
334 throw;
335 }
336 // tab_visco+=DMINFLOAT;
337
338 int ndeb=0,nfin=0;
339 int ori;
340 int elem;
341
342 int signe;
343 int itmax=0;
344
345 double gradient_de_pression0 = 0., gradient_de_pression1 = 0., vmoy = 0., ts0 =0., ts1=0.;
346 const DoubleTab& vit = eqn_hydr.inconnue().valeurs(); //vitesse
347 Navier_Stokes_std& eqnNS = ref_cast_non_const(Navier_Stokes_std, eqn_hydr);
348
349 // Calcul du gradient de pression
350 DoubleTab grad_p(vit);
351 const DoubleTab& p = eqnNS.pression().valeurs();
352 const Operateur_Grad& gradient = eqnNS.operateur_gradient();
353 gradient.calculer(p, grad_p);
354 const DoubleTab& visco_turb = mon_modele_turb_hyd->viscosite_turbulente().valeurs();
355
356
357 DoubleTab termes_sources;
358 termes_sources.resize(domaine_VDF.nb_faces(),1);
359 termes_sources = 0.;
360 // On calcule les termes sources, sauf celui de Boussinesq (TBLE recalcule par lui meme ce terme s'il est demande)
361 const Sources& les_sources=eqn_hydr.sources();
362 int nb_sources=les_sources.size();
363 for (int j=0; j<nb_sources; j++)
364 {
365 const Source_base& ts = les_sources(j).valeur();
366 if (!sub_type(Terme_Boussinesq_base,ts))
367 {
368 ts.ajouter(termes_sources);
369 }
370 }
371
372
373 const double tps = eqnNS.schema_temps().temps_courant();
374 const double dt = eqnNS.schema_temps().pas_de_temps();
375 const double dt_min = eqnNS.schema_temps().pas_temps_min();
376
377 int compteur_faces_paroi = 0;
378
379
380
381
382 compteur_faces_paroi = 0; //Reinitialisation de compteur_faces_paroi
383
384 for (int n_bord=0; n_bord<domaine_VDF.nb_front_Cl(); n_bord++)
385 {
386
387 // pour chaque condition limite on regarde son type
388 // On applique les lois de paroi uniquement
389 // aux voisinages des parois
390
391 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
392
393 if (sub_type(Dirichlet_paroi_fixe,la_cl.valeur()) )
394
395 {
396 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
397 ndeb = le_bord.num_premiere_face();
398 nfin = ndeb + le_bord.nb_faces();
399
400 ///////////////////////
401 // Si probleme 2D //
402 ///////////////////////
403
404
405 if(dimension == 2)
406 {
407 //eq_vit.dimensionner(2*le_dom_dis_->nb_faces_bord());
408
409 //Boucle sur les faces des bords parietaux
410
411 for (int num_face=ndeb; num_face<nfin; num_face++)
412 {
413 ori = orientation(num_face);
414
415 //Selection de l'element associe a une face parietale
416 //et initialisation de l'entier signe
417
418 if ((elem = face_voisins(num_face,0)) == -1)
419 {
420 elem = face_voisins(num_face,1);
421 signe = 1;
422 }
423 else signe = -1 ;
424
425
426 if(dt<dt_min)
427 eq_vit[corresp[num_face]].set_dt(1.e12); // Ca ne devrait pas servir ???
428 else
429 eq_vit[corresp[num_face]].set_dt(dt);
430
431
432
433 // 1er couple de faces perpendiculaires a la paroi
434 int face1 = elem_faces(elem,(ori+1));
435 int face2 = elem_faces(elem,(ori+3)%4);
436
437
438
439
440
441 // Terme source de Boussinesq si demande
442 DoubleVect ts_boussi(nb_pts+1);
443 ts_boussi = 0.;
444 if (source_boussinesq==1)
445 {
446 const Champ_Don_base& ch_gravite=le_fluide.gravite();
447 const DoubleVect& gravite = ch_gravite.valeurs();
448 const Equation_base& eqn_th = mon_modele_turb_hyd->equation().probleme().equation(1);
449 const Modele_turbulence_scal_base& le_mod_turb_th = ref_cast(Modele_turbulence_scal_base,eqn_th.get_modele(TURBULENCE).valeur());
450 const Turbulence_paroi_scal_base& loi = le_mod_turb_th.loi_paroi();
451 Paroi_TBLE_scal_VDF& loi_tble_T = ref_cast_non_const(Paroi_TBLE_scal_VDF,loi);
452
453 if (loi_tble_T.est_initialise())
454 {
455
456 double norm_n=domaine_VDF.face_surfaces(num_face);
457 double gn=0.;
458
459 for (int idim=0; idim<dimension; idim++)
460 {
461 gn+=gravite(idim)*domaine_VDF.face_normales(num_face,idim)/norm_n;
462 }
463
464 DoubleVect gt_vect(dimension);
465
466 for (int idim=0; idim<dimension; idim++)
467 {
468 gt_vect(idim) = gravite(idim)-gn*domaine_VDF.face_normales(num_face,idim)/norm_n;
469 }
470
471 double g_t=gt_vect(1-ori);
472
473 for(int i=0; i<nb_pts+1; i++)
474 {
475 ts_boussi(i) = -g_t*beta_t*(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i)-T0);
476 }
477 }
478 }
479
480 ////////////////////////////////////////////////////////////
481 ////// couple de faces perpendiculaires a la paroi /////
482 ////////////////////////////////////////////////////////////
483
484 vmoy = 0.5*(vit(face1) + vit(face2));
485
486 //vitesse sur le maillage fin a l'instant initial
487
488 ts0 = 0.5*(termes_sources(face1)/volumes_entrelaces(face1)
489 + termes_sources(face2)/volumes_entrelaces(face2));
490
491 eq_vit[corresp[num_face]].set_u_y0(0,0.); //1ere composante de la vitesse a la paroi
492
493 eq_vit[corresp[num_face]].set_u_yn(0,vmoy); //vitesse en yn
494
495 //*** Gradient de pression du 1er couple de faces
496 // perpendiculaires a la paroi ***
497
498 //gradient de pression (constant dans le maillage fin)
499 gradient_de_pression0 = 0.5*(grad_p(face1)/volumes_entrelaces(face1)+grad_p(face2)/volumes_entrelaces(face2));
500
501 for(int i=0 ; i<nb_pts+1; i++)
502 eq_vit[compteur_faces_paroi].set_F(0, i, -gradient_de_pression0+ts0+ts_boussi(i));
503
504
505 //On resoud les equations aux limites simplifiees
506 //pendant un pas de temps du maillage grossier
507
508 if (statio==0)
509 eq_vit[compteur_faces_paroi].aller_au_temps(tps);
510 else
511 eq_vit[compteur_faces_paroi].aller_jusqu_a_convergence(max_it_statio, eps_statio);
512
513
514 if(itmax<eq_vit[corresp[num_face]].get_it())
515 {
516 itmax = eq_vit[compteur_faces_paroi].get_it();
517 if(itmax>max_it)
518 Cerr << "WARNING : TOO MANY ITERATIONS ARE NEEDED IN THE TBLE MESH"
519 << finl;
520 }
521
522
523 ///////////////////////////////////////////////
524 ///// Determination des deux composantes //////
525 ///// du cisaillement a la paroi //////
526 ///////////////////////////////////////////////
527
528 if (ori == 0)
529 {
530 //Cerr << "ori=0" << finl;
531 Cisaillement_paroi_(num_face,0) = 0.;
532 Cisaillement_paroi_(num_face,1) = eq_vit[corresp[num_face]].get_cis(0)*signe;
533 }
534 else if (ori == 1)
535 {
536
537 Cisaillement_paroi_(num_face,0) = eq_vit[corresp[num_face]].get_cis(0)*signe;
538 Cisaillement_paroi_(num_face,1) = 0.;
539 }
540
541 double dist;
542
543 if (axi)
544 dist = domaine_VDF.dist_norm_bord_axi(num_face);
545 else
546 dist = domaine_VDF.dist_norm_bord(num_face);
547
548 double d_visco;
549
550 if (l_unif==1)
551 d_visco = visco;
552 else
553 d_visco = tab_visco[elem];
554
555 tab_u_star_(num_face) = pow(eq_vit[corresp[num_face]].get_cis(0)*eq_vit[corresp[num_face]].get_cis(0),0.25);
556 tab_d_plus_(num_face) = dist*tab_u_star(num_face)/d_visco;
557
558
559 double y_plus = tab_d_plus_(num_face);
560
561 if (y_plus<8)
562 {
563 tab_k(elem) =0.;
564 tab_eps(elem)=0.;
565 }
566 else if (y_plus>8 && y_plus<30)
567 {
568 double u_star = tab_u_star(num_face);
569 double lm_plus = calcul_lm_plus(y_plus);
570 double deriv = Fdypar_direct(lm_plus);
571 double x = lm_plus*u_star*deriv;
572 tab_k(elem) = x*x/sqrt(Cmu_) ;
573 tab_eps(elem) = (tab_k(elem)*u_star*u_star*deriv)*sqrt(Cmu_)/d_visco;
574 }
575 else if (y_plus>30)
576 {
577 double us2 = tab_u_star(num_face)*tab_u_star(num_face);
578 tab_k(elem) =us2/sqrt(Cmu_);
579 tab_eps(elem)=us2*tab_u_star(num_face)/Kappa_/dist;
580 }
581
582
583 compteur_faces_paroi++;
584
585 }//fin boucle sur les faces de bords parietaux
586
587
588 }//fin if(dim2)
589
590
591 ///////////////////////
592 // Si probleme 3D //
593 ///////////////////////
594
595 else if(dimension == 3)
596 {
597 //Boucle sur les faces des bords parietaux
598
599 for (int num_face=ndeb; num_face<nfin; num_face++)
600 {
601 //ici ori=direction perpendiculaire a la paroi
602 ori = orientation(num_face);
603
604
605 //compteur++;
606 //Selection de l'element associe a une face parietale
607 //et initialisation de l'entier signe
608
609 if ((elem = face_voisins(num_face,0)) == -1)
610 {
611 elem = face_voisins(num_face,1);
612 signe = 1;
613 }
614 else signe = -1 ;
615
616
617 if(dt<dt_min)
618 eq_vit[corresp[num_face]].set_dt(1.e12); // Ca ne devrait pas servir ???
619 else
620 eq_vit[corresp[num_face]].set_dt(dt);
621
622
623 // 1er couple de faces perpendiculaires a la paroi
624 int face1 = elem_faces(elem,(ori+1));
625 int face2 = elem_faces(elem,(ori+4)%6);
626 // 2e couple de faces perpendiculaires a la paroi
627 int face3 = elem_faces(elem,(ori+2));
628 int face4 = elem_faces(elem,(ori+5)%6);
629
630
631 ////////////////////////////////////////////////////////////
632 ////// 1er couple de faces perpendiculaires a la paroi /////
633 ////////////////////////////////////////////////////////////
634
635 vmoy = 0.5*(vit(face1) + vit(face2));
636
637 //vitesse sur le maillage fin a l'instant initial
638
639 ts0 = 0.5*(termes_sources(face1)/volumes_entrelaces(face1)
640 + termes_sources(face2)/volumes_entrelaces(face2));
641
642 eq_vit[corresp[num_face]].set_u_y0(0,0.); //1ere composante de la vitesse a la paroi
643
644 eq_vit[corresp[num_face]].set_u_yn(0,vmoy); //vitesse en yn
645
646 //*** Gradient de pression du 1er couple de faces
647 // perpendiculaires a la paroi ***
648
649 //gradient de pression (constant dans le maillage fin)
650 gradient_de_pression0 = 0.5*(grad_p(face1)/volumes_entrelaces(face1)+grad_p(face2)/volumes_entrelaces(face2));
651
652 ////////////////////////////////////////////////////////////
653 ////// 2e couple de faces perpendiculaires a la paroi /////
654 ////////////////////////////////////////////////////////////
655
656 vmoy = 0.5*(vit(face3) + vit(face4));
657
658 ts1 = 0.5*(termes_sources(face3)/volumes_entrelaces(face3)
659 + termes_sources(face4)/volumes_entrelaces(face4));
660
661 eq_vit[corresp[num_face]].set_u_y0(1,0.); //2e composante de la vitesse a la paroi
662
663 eq_vit[corresp[num_face]].set_u_yn(1,vmoy); //vitesse en yn
664
665 //*** Gradient de pression du 2e couple de faces
666 // perpendiculaires a la paroi ***
667
668
669 //gradient de pression (constant dans le maillage fin)
670
671 gradient_de_pression1 = 0.5*(grad_p(face3)/volumes_entrelaces(face3)+grad_p(face4)/volumes_entrelaces(face4));
672
673 // Terme source de Boussinesq si demande
674 DoubleTab ts_boussi(nb_pts+1,2);
675 ts_boussi = 0.;
676 if (source_boussinesq==1)
677 {
678 const Champ_Don_base& ch_gravite=le_fluide.gravite();
679 const DoubleVect& gravite = ch_gravite.valeurs();
680 const Equation_base& eqn_th = mon_modele_turb_hyd->equation().probleme().equation(1);
681 const Modele_turbulence_scal_base& le_mod_turb_th = ref_cast(Modele_turbulence_scal_base,eqn_th.get_modele(TURBULENCE).valeur());
682 const Turbulence_paroi_scal_base& loi = le_mod_turb_th.loi_paroi();
683 Paroi_TBLE_scal_VDF& loi_tble_T = ref_cast_non_const(Paroi_TBLE_scal_VDF,loi);
684 if (loi_tble_T.est_initialise())
685 {
686 double norm_n=domaine_VDF.face_surfaces(num_face);
687 double gn=0.;
688
689 for (int idim=0; idim<dimension; idim++)
690 {
691 gn+=gravite(idim)*domaine_VDF.face_normales(num_face,idim)/norm_n;
692 }
693
694 DoubleVect gt_vect(dimension);
695
696 for (int idim=0; idim<dimension; idim++)
697 {
698 gt_vect(idim) = gravite(idim)-gn*domaine_VDF.face_normales(num_face,idim)/norm_n;
699 }
700
701 for(int i=0; i<nb_pts+1; i++)
702 {
703 ts_boussi(i,0) = -gt_vect((ori+1)%dimension)*beta_t*(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i)-T0);
704 ts_boussi(i,1) = -gt_vect((ori+2)%dimension)*beta_t*(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i)-T0);
705 }
706 }
707 }
708
709
710
711 if(alpha_cv == 1) //Si TBLE + termes convectifs
712 calculer_convection(num_face,face1, face2, face3, face4, elem, ndeb, nfin, ori, gradient_de_pression0, ts0, gradient_de_pression1, ts1);
713 else
714 {
715 for(int i=0 ; i<nb_pts+1 ; i++)
716 {
717 eq_vit[corresp[num_face]].set_F(0, i, - gradient_de_pression0 + ts0 + ts_boussi(i,0));
718 eq_vit[corresp[num_face]].set_F(1, i, - gradient_de_pression1 + ts1 + ts_boussi(i,1));
719 }
720 }
721
722 if(tps<=dt_min)
723 {
724 visco_turb_moy(elem)=visco_turb(elem);
725 }
726 else if(tps>tps_start_stat_nu_t_dyn)
727 {
728 visco_turb_moy(elem) = (dt/tps)*visco_turb(elem)+(1-(dt/tps))*visco_turb_moy(elem);
729 }
730
731 //On resoud les equations aux limites simplifiees
732 //pendant un pas de temps du maillage grossier
733 if (statio==0)
734 eq_vit[corresp[num_face]].aller_au_temps(tps);
735 else
736 eq_vit[corresp[num_face]].aller_jusqu_a_convergence(max_it_statio, eps_statio);
737
738 if(itmax<eq_vit[corresp[num_face]].get_it())
739 {
740 itmax = eq_vit[corresp[num_face]].get_it();
741 if(itmax>max_it)
742 Cerr << "WARNING : TOO MANY ITERATIONS ARE NEEDED IN THE TBLE MESH"
743 << finl;
744 }
745
746 ///////////////////////////////////////////////
747 ///// Determination des deux composantes //////
748 ///// du cisaillement a la paroi //////
749 ///////////////////////////////////////////////
750
751 if (ori == 0)
752 {
753 Cisaillement_paroi_(num_face,0) = 0.;
754 Cisaillement_paroi_(num_face,1) = eq_vit[corresp[num_face]].get_cis(0)*signe;
755 Cisaillement_paroi_(num_face,2) = eq_vit[corresp[num_face]].get_cis(1)*signe;
756 }
757 else if (ori == 1)
758 {
759 Cisaillement_paroi_(num_face,0) = eq_vit[corresp[num_face]].get_cis(1)*signe;
760 Cisaillement_paroi_(num_face,1) = 0.;
761 Cisaillement_paroi_(num_face,2) = eq_vit[corresp[num_face]].get_cis(0)*signe;
762 }
763 else
764 {
765 Cisaillement_paroi_(num_face,0) = eq_vit[corresp[num_face]].get_cis(0)*signe;
766 Cisaillement_paroi_(num_face,1) = eq_vit[corresp[num_face]].get_cis(1)*signe;
767 Cisaillement_paroi_(num_face,2) = 0.;
768 }
769
770 double dist;
771
772 if (axi)
773 dist = domaine_VDF.dist_norm_bord_axi(num_face);
774 else
775 dist = domaine_VDF.dist_norm_bord(num_face);
776 double d_visco;
777
778 if (l_unif==1)
779 d_visco = visco;
780 else
781 d_visco = tab_visco[elem];
782
783 tab_u_star_(num_face) = sqrt(sqrt(eq_vit[corresp[num_face]].get_cis(0)*eq_vit[corresp[num_face]].get_cis(0)
784 +eq_vit[corresp[num_face]].get_cis(1)*eq_vit[corresp[num_face]].get_cis(1)));
785 tab_d_plus_(num_face) = dist*tab_u_star(num_face)/d_visco;
786
787
788
789 double y_plus = tab_d_plus(num_face);
790
791 if (y_plus<8)
792 {
793 tab_k(elem) =0.;
794 tab_eps(elem)=0.;
795 }
796 else if (y_plus>8 && y_plus<30)
797 {
798 double u_star = tab_u_star(num_face);
799 double lm_plus = calcul_lm_plus(y_plus);
800 double deriv = Fdypar_direct(lm_plus);
801 double x = lm_plus*u_star*deriv;
802 tab_k(elem) = x*x/sqrt(Cmu_);
803 tab_eps(elem) = (tab_k(elem)*u_star*u_star*deriv)*sqrt(Cmu_)/d_visco;
804 }
805 else if (y_plus>30)
806 {
807 double us2 = tab_u_star(num_face)*tab_u_star(num_face);
808 tab_k(elem)=us2/sqrt(Cmu_);
809 tab_eps(elem)=us2*tab_u_star(num_face)/Kappa_/dist;
810 }
811
812 compteur_faces_paroi++;
813
814 }//fin boucle sur les faces de bords parietaux
815 }//fin if(dim=3)
816 }//fin if(sub_type(dirichlet))
817
818 }//fin boucle sur les bords
820 {
821 SFichier fic("iter.dat",ios::app); // ouverture du fichier iter.dat
822 fic << tps << " " << itmax << finl;
823 }
824 if(oui_stats==1)
825 calculer_stats();
826
827 return 1;
828}
829
830int ParoiVDF_TBLE::calculer_hyd(DoubleTab& tab1,int isKeps,DoubleTab& tab2)
831{
832 // si isKeps=1, on est dans le cas k-eps
833
834 // clock_t clock0,clock1;
835
836 // clock0=clock();
837
838 /* if (Process::nproc()>1)
839 {
840 Cerr << "ParoiVDF_TBLE::calculer_hyd n'est pas parallelise." << finl;
841 Process::exit();
842 }
843 */
844
845 const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_dis_.valeur());
846 const IntVect& orientation = domaine_VDF.orientation();
847 const IntTab& face_voisins = domaine_VDF.face_voisins();
848 const IntTab& elem_faces = domaine_VDF.elem_faces();
849 const DoubleVect& volumes_entrelaces = domaine_VDF.volumes_entrelaces();
850 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
851
852 const Fluide_base& le_fluide = ref_cast(Fluide_base,eqn_hydr.milieu());
853 const Champ_Don_base& ch_visco_cin = le_fluide.viscosite_cinematique();
854 const DoubleTab& tab_visco = ch_visco_cin.valeurs();
855
856 if (source_boussinesq==1)
857 {
858 const Champ_Don_base& ch_gravite=le_fluide.gravite();
859
860 if (!sub_type(Champ_Uniforme,ch_gravite))
861 {
862 Cerr << " On ne sait pas traiter la gravite non uniforme dans TBLE !!!"<< finl;
864 }
865 }
866
867 double visco=-1;
868 int l_unif;
869 if (sub_type(Champ_Uniforme,ch_visco_cin))
870 {
871 visco = std::max(tab_visco(0,0),DMINFLOAT);
872 l_unif = 1;
873 }
874 else
875 l_unif = 0;
876
877 // preparer_calcul_hyd(tab);
878 if ((!l_unif) && (tab_visco.local_min_vect()<DMINFLOAT))
879 // on ne doit pas changer tab_visco ici !
880 {
881 Cerr << "In ParoiVDF_TBLE::calculer_hyd : visco = " << tab_visco.local_min_vect() << " <= 0 ? " << finl;
882 throw;
883 }
884 // tab_visco+=DMINFLOAT;
885
886 int ndeb=0,nfin=0;
887 int ori;
888 int elem;
889
890 int signe;
891 int itmax=0;
892
893 double gradient_de_pression0 = 0., gradient_de_pression1 = 0., vmoy = 0., ts0 =0., ts1=0.;
894 const DoubleTab& vit = eqn_hydr.inconnue().valeurs(); //vitesse
895 Navier_Stokes_std& eqnNS = ref_cast_non_const(Navier_Stokes_std, eqn_hydr);
896
897 // Calcul du gradient de pression
898 DoubleTab grad_p(vit);
899 const DoubleTab& p = eqnNS.pression().valeurs();
900 const Operateur_Grad& gradient = eqnNS.operateur_gradient();
901 gradient.calculer(p, grad_p);
902 const DoubleTab& visco_turb = mon_modele_turb_hyd->viscosite_turbulente().valeurs();
903
904
905 DoubleTab termes_sources;
906 termes_sources.resize(domaine_VDF.nb_faces(),1);
907 termes_sources = 0.;
908 // On calcule les termes sources, sauf celui de Boussinesq (TBLE recalcule par lui meme ce terme s'il est demande)
909 const Sources& les_sources=eqn_hydr.sources();
910 int nb_sources=les_sources.size();
911 for (int j=0; j<nb_sources; j++)
912 {
913 const Source_base& ts = les_sources(j).valeur();
914 if (!sub_type(Terme_Boussinesq_base,ts))
915 {
916 ts.ajouter(termes_sources);
917 }
918 }
919
920
921 const double tps = eqnNS.schema_temps().temps_courant();
922 const double dt = eqnNS.schema_temps().pas_de_temps();
923 const double dt_min = eqnNS.schema_temps().pas_temps_min();
924
925 int compteur_faces_paroi = 0;
926
927
928
929
930 compteur_faces_paroi = 0; //Reinitialisation de compteur_faces_paroi
931
932 for (int n_bord=0; n_bord<domaine_VDF.nb_front_Cl(); n_bord++)
933 {
934
935 // pour chaque condition limite on regarde son type
936 // On applique les lois de paroi uniquement
937 // aux voisinages des parois
938
939 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
940
941 if (sub_type(Dirichlet_paroi_fixe,la_cl.valeur()) )
942
943 {
944 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
945 ndeb = le_bord.num_premiere_face();
946 nfin = ndeb + le_bord.nb_faces();
947
948 ///////////////////////
949 // Si probleme 2D //
950 ///////////////////////
951
952
953 if(dimension == 2)
954 {
955 //eq_vit.dimensionner(2*le_dom_dis_->nb_faces_bord());
956
957 //Boucle sur les faces des bords parietaux
958
959 for (int num_face=ndeb; num_face<nfin; num_face++)
960 {
961 ori = orientation(num_face);
962
963 //Selection de l'element associe a une face parietale
964 //et initialisation de l'entier signe
965
966 if ((elem = face_voisins(num_face,0)) == -1)
967 {
968 elem = face_voisins(num_face,1);
969 signe = 1;
970 }
971 else signe = -1 ;
972
973
974 if(dt<dt_min)
975 eq_vit[corresp[num_face]].set_dt(1.e12); // Ca ne devrait pas servir ???
976 else
977 eq_vit[corresp[num_face]].set_dt(dt);
978
979
980
981 // 1er couple de faces perpendiculaires a la paroi
982 int face1 = elem_faces(elem,(ori+1));
983 int face2 = elem_faces(elem,(ori+3)%4);
984
985
986
987
988
989 // Terme source de Boussinesq si demande
990 DoubleVect ts_boussi(nb_pts+1);
991 ts_boussi = 0.;
992 if (source_boussinesq==1)
993 {
994 const Champ_Don_base& ch_gravite=le_fluide.gravite();
995 const DoubleVect& gravite = ch_gravite.valeurs();
996 const Equation_base& eqn_th = mon_modele_turb_hyd->equation().probleme().equation(1);
997 const Modele_turbulence_scal_base& le_mod_turb_th = ref_cast(Modele_turbulence_scal_base,eqn_th.get_modele(TURBULENCE).valeur());
998 const Turbulence_paroi_scal_base& loi = le_mod_turb_th.loi_paroi();
999 Paroi_TBLE_scal_VDF& loi_tble_T = ref_cast_non_const(Paroi_TBLE_scal_VDF,loi);
1000
1001 if (loi_tble_T.est_initialise())
1002 {
1003
1004 double norm_n=domaine_VDF.face_surfaces(num_face);
1005 double gn=0.;
1006
1007 for (int idim=0; idim<dimension; idim++)
1008 {
1009 gn+=gravite(idim)*domaine_VDF.face_normales(num_face,idim)/norm_n;
1010 }
1011
1012 DoubleVect gt_vect(dimension);
1013
1014 for (int idim=0; idim<dimension; idim++)
1015 {
1016 gt_vect(idim) = gravite(idim)-gn*domaine_VDF.face_normales(num_face,idim)/norm_n;
1017 }
1018
1019 double g_t=gt_vect(1-ori);
1020
1021 for(int i=0; i<nb_pts+1; i++)
1022 {
1023 ts_boussi(i) = -g_t*beta_t*(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i)-T0);
1024 }
1025 /*Cout << "gt " << g_t << finl;
1026 Cout << "T0 " << T0 << finl;
1027 Cout << "beta_t " << beta_t << finl;*/
1028 }
1029 }
1030
1031 ////////////////////////////////////////////////////////////
1032 ////// couple de faces perpendiculaires a la paroi /////
1033 ////////////////////////////////////////////////////////////
1034
1035 vmoy = 0.5*(vit(face1) + vit(face2));
1036
1037 //vitesse sur le maillage fin a l'instant initial
1038
1039 ts0 = 0.5*(termes_sources(face1)/volumes_entrelaces(face1)
1040 + termes_sources(face2)/volumes_entrelaces(face2));
1041
1042 eq_vit[corresp[num_face]].set_u_y0(0,0.); //1ere composante de la vitesse a la paroi
1043
1044 eq_vit[corresp[num_face]].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)/volumes_entrelaces(face1)+grad_p(face2)/volumes_entrelaces(face2));
1051
1052 for(int i=0 ; i<nb_pts+1; i++)
1053 eq_vit[compteur_faces_paroi].set_F(0, i, -gradient_de_pression0+ts0+ts_boussi(i));
1054
1055
1056 //On resoud les equations aux limites simplifiees
1057 //pendant un pas de temps du maillage grossier
1058
1059 if (statio==0)
1060 eq_vit[compteur_faces_paroi].aller_au_temps(tps);
1061 else
1062 eq_vit[compteur_faces_paroi].aller_jusqu_a_convergence(max_it_statio, eps_statio);
1063
1064
1065 if(itmax<eq_vit[corresp[num_face]].get_it())
1066 {
1067 itmax = eq_vit[compteur_faces_paroi].get_it();
1068 if(itmax>max_it)
1069 Cerr << "WARNING : TOO MANY ITERATIONS ARE NEEDED IN THE TBLE MESH"
1070 << finl;
1071 }
1072
1073
1074 ///////////////////////////////////////////////
1075 ///// Determination des deux composantes //////
1076 ///// du cisaillement a la paroi //////
1077 ///////////////////////////////////////////////
1078
1079 if (ori == 0)
1080 {
1081 //Cerr << "ori=0" << finl;
1082 Cisaillement_paroi_(num_face,0) = 0.;
1083 Cisaillement_paroi_(num_face,1) = eq_vit[corresp[num_face]].get_cis(0)*signe;
1084 }
1085 else if (ori == 1)
1086 {
1087
1088 Cisaillement_paroi_(num_face,0) = eq_vit[corresp[num_face]].get_cis(0)*signe;
1089 Cisaillement_paroi_(num_face,1) = 0.;
1090 }
1091
1092 double dist;
1093
1094 if (axi)
1095 dist = domaine_VDF.dist_norm_bord_axi(num_face);
1096 else
1097 dist = domaine_VDF.dist_norm_bord(num_face);
1098
1099 double d_visco;
1100
1101 if (l_unif==1)
1102 d_visco = visco;
1103 else
1104 d_visco = tab_visco[elem];
1105
1106 tab_u_star_(num_face) = pow(eq_vit[corresp[num_face]].get_cis(0)*eq_vit[corresp[num_face]].get_cis(0),0.25);
1107 tab_d_plus_(num_face) = dist*tab_u_star(num_face)/d_visco;
1108
1109 if (isKeps==0) // on est en LES
1110 {
1111
1112 ////////////////////////////////////////////////////////////////////
1113 //nu_t longueur de melange en premiere maille calcule dans Diffu_lm
1114 ////////////////////////////////////////////////////////////////////
1115
1116 tab1(elem)=eq_vit[corresp[num_face]].get_nu_t_yn();
1117 }
1118 else if (isKeps==1) // on est en K-eps
1119 {
1120 double y_plus = tab_d_plus_(num_face);
1121
1122 if (y_plus<8)
1123 {
1124 tab1(elem,0)=0.;
1125 tab1(elem,1)=0.;
1126 }
1127 else if (y_plus>8 && y_plus<30)
1128 {
1129 double u_star = tab_u_star(num_face);
1130 double lm_plus = calcul_lm_plus(y_plus);
1131 double deriv = Fdypar_direct(lm_plus);
1132 double x = lm_plus*u_star*deriv;
1133 tab1(elem,0) = x*x/sqrt(Cmu_) ;
1134 tab1(elem,1) = (tab1(elem,0)*u_star*u_star*deriv)*sqrt(Cmu_)/d_visco;
1135 }
1136 else if (y_plus>30)
1137 {
1138 double us2 = tab_u_star(num_face)*tab_u_star(num_face);
1139 tab1(elem,0)=us2/sqrt(Cmu_);
1140 tab1(elem,1)=us2*tab_u_star(num_face)/Kappa_/dist;
1141 }
1142 }
1143
1144 compteur_faces_paroi++;
1145
1146 }//fin boucle sur les faces de bords parietaux
1147
1148
1149 }//fin if(dim2)
1150
1151
1152 ///////////////////////
1153 // Si probleme 3D //
1154 ///////////////////////
1155
1156 else if(dimension == 3)
1157 {
1158 //Boucle sur les faces des bords parietaux
1159
1160 for (int num_face=ndeb; num_face<nfin; num_face++)
1161 {
1162 //ici ori=direction perpendiculaire a la paroi
1163 ori = orientation(num_face);
1164
1165
1166 //compteur++;
1167 //Selection de l'element associe a une face parietale
1168 //et initialisation de l'entier signe
1169
1170 if ((elem = face_voisins(num_face,0)) == -1)
1171 {
1172 elem = face_voisins(num_face,1);
1173 signe = 1;
1174 }
1175 else signe = -1 ;
1176
1177
1178 if(dt<dt_min)
1179 eq_vit[corresp[num_face]].set_dt(1.e12); // Ca ne devrait pas servir ???
1180 else
1181 eq_vit[corresp[num_face]].set_dt(dt);
1182
1183
1184 // 1er couple de faces perpendiculaires a la paroi
1185 int face1 = elem_faces(elem,(ori+1));
1186 int face2 = elem_faces(elem,(ori+4)%6);
1187 // 2e couple de faces perpendiculaires a la paroi
1188 int face3 = elem_faces(elem,(ori+2));
1189 int face4 = elem_faces(elem,(ori+5)%6);
1190
1191
1192 ////////////////////////////////////////////////////////////
1193 ////// 1er couple de faces perpendiculaires a la paroi /////
1194 ////////////////////////////////////////////////////////////
1195
1196 vmoy = 0.5*(vit(face1) + vit(face2));
1197
1198 //vitesse sur le maillage fin a l'instant initial
1199
1200 ts0 = 0.5*(termes_sources(face1)/volumes_entrelaces(face1)
1201 + termes_sources(face2)/volumes_entrelaces(face2));
1202
1203 eq_vit[corresp[num_face]].set_u_y0(0,0.); //1ere composante de la vitesse a la paroi
1204
1205 eq_vit[corresp[num_face]].set_u_yn(0,vmoy); //vitesse en yn
1206
1207 //*** Gradient de pression du 1er couple de faces
1208 // perpendiculaires a la paroi ***
1209
1210 //gradient de pression (constant dans le maillage fin)
1211 gradient_de_pression0 = 0.5*(grad_p(face1)/volumes_entrelaces(face1)+grad_p(face2)/volumes_entrelaces(face2));
1212
1213 ////////////////////////////////////////////////////////////
1214 ////// 2e couple de faces perpendiculaires a la paroi /////
1215 ////////////////////////////////////////////////////////////
1216
1217 vmoy = 0.5*(vit(face3) + vit(face4));
1218
1219 ts1 = 0.5*(termes_sources(face3)/volumes_entrelaces(face3)
1220 + termes_sources(face4)/volumes_entrelaces(face4));
1221
1222 eq_vit[corresp[num_face]].set_u_y0(1,0.); //2e composante de la vitesse a la paroi
1223
1224 eq_vit[corresp[num_face]].set_u_yn(1,vmoy); //vitesse en yn
1225
1226 //*** Gradient de pression du 2e couple de faces
1227 // perpendiculaires a la paroi ***
1228
1229
1230 //gradient de pression (constant dans le maillage fin)
1231
1232 gradient_de_pression1 = 0.5*(grad_p(face3)/volumes_entrelaces(face3)+grad_p(face4)/volumes_entrelaces(face4));
1233
1234 // Terme source de Boussinesq si demande
1235 DoubleTab ts_boussi(nb_pts+1,2);
1236 ts_boussi = 0.;
1237 if (source_boussinesq==1)
1238 {
1239 const Champ_Don_base& ch_gravite=le_fluide.gravite();
1240 const DoubleVect& gravite = ch_gravite.valeurs();
1241 const Equation_base& eqn_th = mon_modele_turb_hyd->equation().probleme().equation(1);
1242 const Modele_turbulence_scal_base& le_mod_turb_th = ref_cast(Modele_turbulence_scal_base,eqn_th.get_modele(TURBULENCE).valeur());
1243 const Turbulence_paroi_scal_base& loi = le_mod_turb_th.loi_paroi();
1244 Paroi_TBLE_scal_VDF& loi_tble_T = ref_cast_non_const(Paroi_TBLE_scal_VDF,loi);
1245 if (loi_tble_T.est_initialise())
1246 {
1247 double norm_n=domaine_VDF.face_surfaces(num_face);
1248 double gn=0.;
1249
1250 for (int idim=0; idim<dimension; idim++)
1251 {
1252 gn+=gravite(idim)*domaine_VDF.face_normales(num_face,idim)/norm_n;
1253 }
1254
1255 DoubleVect gt_vect(dimension);
1256
1257 for (int idim=0; idim<dimension; idim++)
1258 {
1259 gt_vect(idim) = gravite(idim)-gn*domaine_VDF.face_normales(num_face,idim)/norm_n;
1260 }
1261
1262 for(int i=0; i<nb_pts+1; i++)
1263 {
1264 ts_boussi(i,0) = -gt_vect((ori+1)%dimension)*beta_t*(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i)-T0);
1265 ts_boussi(i,1) = -gt_vect((ori+2)%dimension)*beta_t*(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i)-T0);
1266 /*Cout << "gt_vect((ori+1)%dimension) " << gt_vect((ori+1)%dimension) << finl;
1267 Cout << "gt_vect((ori+2)%dimension) " << gt_vect((ori+2)%dimension) << finl;
1268 Cout << "vit(face1) " << vit(face1) << finl;
1269 Cout << "vit(face3) " << vit(face3) << finl;
1270 Cout << "T0 " << T0 << finl;
1271 Cout << "beta_t " << beta_t << finl;*/
1272 }
1273 }
1274 }
1275
1276
1277
1278 if(alpha_cv == 1) //Si TBLE + termes convectifs
1279 calculer_convection(num_face,face1, face2, face3, face4, elem, ndeb, nfin, ori, gradient_de_pression0, ts0, gradient_de_pression1, ts1);
1280 else
1281 {
1282 for(int i=0 ; i<nb_pts+1 ; i++)
1283 {
1284 eq_vit[corresp[num_face]].set_F(0, i, - gradient_de_pression0 + ts0 + ts_boussi(i,0));
1285 eq_vit[corresp[num_face]].set_F(1, i, - gradient_de_pression1 + ts1 + ts_boussi(i,1));
1286 }
1287 }
1288
1289 if(tps<=dt_min)
1290 {
1291 visco_turb_moy(elem)=visco_turb(elem);
1292 }
1293 else if(tps>tps_start_stat_nu_t_dyn)
1294 {
1295 visco_turb_moy(elem) = (dt/tps)*visco_turb(elem)+(1-(dt/tps))*visco_turb_moy(elem);
1296 }
1297
1298 //On resoud les equations aux limites simplifiees
1299 //pendant un pas de temps du maillage grossier
1300 if (statio==0)
1301 eq_vit[corresp[num_face]].aller_au_temps(tps);
1302 else
1303 eq_vit[corresp[num_face]].aller_jusqu_a_convergence(max_it_statio, eps_statio);
1304
1305 if(itmax<eq_vit[corresp[num_face]].get_it())
1306 {
1307 itmax = eq_vit[corresp[num_face]].get_it();
1308 if(itmax>max_it)
1309 Cerr << "WARNING : TOO MANY ITERATIONS ARE NEEDED IN THE TBLE MESH"
1310 << finl;
1311 }
1312
1313 ///////////////////////////////////////////////
1314 ///// Determination des deux composantes //////
1315 ///// du cisaillement a la paroi //////
1316 ///////////////////////////////////////////////
1317
1318 if (ori == 0)
1319 {
1320 Cisaillement_paroi_(num_face,0) = 0.;
1321 Cisaillement_paroi_(num_face,1) = eq_vit[corresp[num_face]].get_cis(0)*signe;
1322 Cisaillement_paroi_(num_face,2) = eq_vit[corresp[num_face]].get_cis(1)*signe;
1323 }
1324 else if (ori == 1)
1325 {
1326 Cisaillement_paroi_(num_face,0) = eq_vit[corresp[num_face]].get_cis(1)*signe;
1327 Cisaillement_paroi_(num_face,1) = 0.;
1328 Cisaillement_paroi_(num_face,2) = eq_vit[corresp[num_face]].get_cis(0)*signe;
1329 }
1330 else
1331 {
1332 Cisaillement_paroi_(num_face,0) = eq_vit[corresp[num_face]].get_cis(0)*signe;
1333 Cisaillement_paroi_(num_face,1) = eq_vit[corresp[num_face]].get_cis(1)*signe;
1334 Cisaillement_paroi_(num_face,2) = 0.;
1335 }
1336
1337 double dist;
1338
1339 if (axi)
1340 dist = domaine_VDF.dist_norm_bord_axi(num_face);
1341 else
1342 dist = domaine_VDF.dist_norm_bord(num_face);
1343 double d_visco;
1344
1345 if (l_unif==1)
1346 d_visco = visco;
1347 else
1348 d_visco = tab_visco[elem];
1349
1350 tab_u_star_(num_face) = sqrt(sqrt(eq_vit[corresp[num_face]].get_cis(0)*eq_vit[corresp[num_face]].get_cis(0)
1351 +eq_vit[corresp[num_face]].get_cis(1)*eq_vit[corresp[num_face]].get_cis(1)));
1352 tab_d_plus_(num_face) = dist*tab_u_star(num_face)/d_visco;
1353
1354
1355 if (isKeps==0) // on est en LES
1356 {
1357 if((nu_t_dyn==0)||(tps<tps_nu_t_dyn))
1358 {
1359 tab1(elem)=eq_vit[corresp[num_face]].get_nu_t_yn();
1360 }
1361 else
1362 eq_vit[corresp[num_face]].set_nu_t_yn(visco_turb_moy(elem));
1363 }
1364 else if (isKeps==1) // on est en K-eps
1365 {
1366 double y_plus = tab_d_plus(num_face);
1367
1368 if (y_plus<8)
1369 {
1370 tab1(elem,0)=0.;
1371 tab1(elem,1)=0.;
1372 }
1373 else if (y_plus>8 && y_plus<30)
1374 {
1375 double u_star = tab_u_star(num_face);
1376 double lm_plus = calcul_lm_plus(y_plus);
1377 double deriv = Fdypar_direct(lm_plus);
1378 double x = lm_plus*u_star*deriv;
1379 tab1(elem,0) = x*x/sqrt(Cmu_);
1380 tab1(elem,1) = (tab1(elem,0)*u_star*u_star*deriv)*sqrt(Cmu_)/d_visco;
1381 }
1382 else if (y_plus>30)
1383 {
1384 double us2 = tab_u_star(num_face)*tab_u_star(num_face);
1385 tab1(elem,0)=us2/sqrt(Cmu_);
1386 tab1(elem,1)=us2*tab_u_star(num_face)/Kappa_/dist;
1387 }
1388 }
1389
1390 compteur_faces_paroi++;
1391
1392 }//fin boucle sur les faces de bords parietaux
1393 }//fin if(dim=3)
1394 }//fin if(sub_type(dirichlet))
1395
1396 }//fin boucle sur les bords
1398 {
1399 SFichier fic("iter.dat", ios::app); // ouverture du fichier iter.dat
1400 fic << tps << " " << itmax << finl;
1401 }
1402
1403 if(oui_stats==1)
1404 calculer_stats();
1405
1406
1407 // clock1=clock();
1408 // double time = (clock1-clock0)*1.e-6;
1409 // if(je_suis_maitre()) { Cerr << " CPU time TBLE : " << Process::mp_max(time) << finl; }
1410
1411 return 1;
1412}
1413
1414
1415int ParoiVDF_TBLE::calculer_stats()
1416{
1417 const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_dis_.valeur());
1418 const IntVect& orientation = domaine_VDF.orientation();
1419
1420 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
1421 const double tps = eqn_hydr.inconnue().temps();
1422 const double dt = eqn_hydr.schema_temps().pas_de_temps();
1423
1424
1425 //////////////////////////////////////
1426 //Moyennes Temporelles
1427 //////////////////////////////////////
1428 if(((tps>tps_deb_stats) && (tps<tps_fin_stats)) && (oui_stats!=0))
1429 {
1430 if (dimension==2)
1431 {
1432 for(int j=0; j<nb_post_pts; j++)
1433 {
1434 for(int i=0 ; i<nb_pts+1 ; i++)
1435 {
1436 double u = eq_vit[num_faces_post(j)].get_Unp1(0,i);
1437 double F=eq_vit[num_faces_post(j)].get_F(0,i);
1438 calculer_stat(j,i,F,0,0,u,0.,0.,dt);
1439 }
1440 }
1441 }
1442 else if (dimension==3)
1443 {
1444 for(int j=0; j<nb_post_pts; j++)
1445 {
1446 int ori=orientation(num_global_faces_post(j));
1447 if (ori == 0)
1448 {
1449 for(int i=0 ; i<nb_pts+1 ; i++)
1450 {
1451 double u=eq_vit[num_faces_post(j)].get_v(i);
1452 double v=eq_vit[num_faces_post(j)].get_Unp1(0,i);
1453 double w=eq_vit[num_faces_post(j)].get_Unp1(1,i);
1454 double Fy=eq_vit[num_faces_post(j)].get_F(0,i);
1455 double Fz=eq_vit[num_faces_post(j)].get_F(1,i);
1456 calculer_stat(j,i,0, Fy, Fz, u,v,w,dt);
1457 }
1458 }
1459 else if (ori == 1)
1460 {
1461 for(int i=0 ; i<nb_pts+1 ; i++)
1462 {
1463 double u=eq_vit[num_faces_post(j)].get_Unp1(1,i);
1464 double v=eq_vit[num_faces_post(j)].get_v(i);
1465 double w=eq_vit[num_faces_post(j)].get_Unp1(0,i);
1466 double Fx=eq_vit[num_faces_post(j)].get_F(1,i);
1467 double Fz=eq_vit[num_faces_post(j)].get_F(0,i);
1468 calculer_stat(j,i,Fx, 0, Fz, u,v,w,dt);
1469 }
1470 }
1471 else if (ori == 2)
1472 {
1473 for(int i=0 ; i<nb_pts+1 ; i++)
1474 {
1475 double u=eq_vit[num_faces_post(j)].get_Unp1(0,i);
1476 double v=eq_vit[num_faces_post(j)].get_Unp1(1,i);
1477 double w=eq_vit[num_faces_post(j)].get_v(i);
1478 double Fx=eq_vit[num_faces_post(j)].get_F(0,i);
1479 double Fy=eq_vit[num_faces_post(j)].get_F(1,i);
1480 calculer_stat(j,i,Fx,Fy,0,u,v,w,dt);
1481 }
1482 }
1483 }
1484 }
1485 }
1486
1487 return 1;
1488}
1489
1491{
1492 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
1493 const double tps = eqn_hydr.inconnue().temps();
1494
1495 const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_dis_.valeur());
1496 const IntVect& orientation = domaine_VDF.orientation();
1497
1498 for(int j=0; j<nb_post_pts; j++)
1499 {
1500 Nom tmp;
1501 tmp="tble_post_";
1502 tmp+=nom_pts[j];
1503 tmp+=".dat";
1504 EcrFicPartage fic_post(tmp, ios::app);
1505
1506 int ori=orientation(num_global_faces_post(j));
1507
1508 fic_post << "# "<< tps << " " << "u*= " << tab_u_star(num_global_faces_post(j)) << finl;
1509 if (dimension==2)
1510 {
1511 for(int i=0; i<nb_pts+1; i++)
1512 fic_post << eq_vit[num_faces_post(j)].get_y(i) << " " << eq_vit[num_faces_post(j)].get_Unp1(0,i) << finl;
1513 }
1514 else if (dimension==3)
1515 {
1516 if (ori == 0)
1517 for(int i=0; i<nb_pts+1; i++)
1518 fic_post << eq_vit[num_faces_post(j)].get_y(i) << " 0. " << eq_vit[num_faces_post(j)].get_Unp1(0,i) <<
1519 " " << eq_vit[num_faces_post(j)].get_Unp1(1,i) << finl;
1520 else if (ori == 1)
1521 for(int i=0; i<nb_pts+1; i++)
1522 fic_post << eq_vit[num_faces_post(j)].get_y(i) << " " << eq_vit[num_faces_post(j)].get_Unp1(1,i) <<
1523 " 0. " << eq_vit[num_faces_post(j)].get_Unp1(0,i) << finl;
1524 else if (ori == 2)
1525 for(int i=0; i<nb_pts+1; i++)
1526 fic_post << eq_vit[num_faces_post(j)].get_y(i) << " " << eq_vit[num_faces_post(j)].get_Unp1(0,i) <<
1527 " " << eq_vit[num_faces_post(j)].get_Unp1(1,i) << " 0." << finl;
1528 }
1529 fic_post.syncfile();
1530 fic_post.close();
1531 }
1532
1534
1535}
1536
1538{
1539 const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_dis_.valeur());
1540 double tps = mon_modele_turb_hyd->equation().inconnue().temps();
1541 return Paroi_TBLE_QDM::sauvegarder(os, domaine_VDF, le_dom_Cl_dis_.valeur(), tps);
1542}
1543
1544
1546{
1547 if (le_dom_dis_) // test pour ne pas planter dans "avancer_fichier(...)"
1548 {
1549 const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_dis_.valeur());
1550 double tps_reprise = mon_modele_turb_hyd->equation().schema_temps().temps_courant();
1551 return Paroi_TBLE_QDM::reprendre(is, domaine_VDF, le_dom_Cl_dis_.valeur(), tps_reprise);
1552 }
1553 else return 1;
1554}
1555
1557{
1558 const Probleme_base& pb_base = mon_modele_turb_hyd->equation().probleme();
1559 return pb_base;
1560}
1561
1562void ParoiVDF_TBLE::calculer_convection(int num_face, int face1, int face2, int face3, int face4, int elem, int ndeb, int nfin, int ori, double gradient_de_pression0, double ts0, double gradient_de_pression1, double ts1)
1563{
1564 const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_dis_.valeur());
1565 const IntTab& face_voisins = domaine_VDF.face_voisins();
1566 const IntTab& elem_faces = domaine_VDF.elem_faces();
1567 double d0d0_1=0;
1568 double d1d1_1=0;
1569 double d00d0 = 0., d0vdy= 0., d01d1= 0., d10d0= 0., d1vdy= 0., d11d1=0.;
1570 int elem_gauche0, elem_gauche1, elem_droite0, elem_droite1;
1571 int face_gauche0, face_gauche1, face_droite0, face_droite1;
1572
1573 ////////////////////////////////////////////////////////////////////////////////////
1574 //Calcul de la vitesse normale a la paroi a partir de la conservation de la masse
1575 ////////////////////////////////////////////////////////////////////////////////////
1576 double v = 0.;
1577
1578 //On cherche les elements voisins de l'element considere
1579
1580 elem_gauche0 = face_voisins(face1,0);
1581 elem_droite0 = face_voisins(face2,1);
1582 if (elem_gauche0 == elem)
1583 {
1584 elem_droite0 = face_voisins(face1,1);
1585 elem_gauche0 = face_voisins(face2,0);
1586 }
1587
1588 elem_gauche1 = face_voisins(face3,0);
1589 elem_droite1 = face_voisins(face4,1);
1590 if (elem_gauche1 == elem)
1591 {
1592 elem_droite1 = face_voisins(face3,1);
1593 elem_gauche1 = face_voisins(face4,0);
1594 }
1595
1596 //On evite les coins interieurs
1597 if((elem_droite1!=-1)&&(elem_droite0!=-1)&&(elem_gauche1!=-1)&&(elem_gauche0!=-1))
1598 {
1599 //On cherche les faces parietales voisines de la face
1600 face_gauche0 = elem_faces(elem_gauche0, ori);
1601 if ((face_voisins(face_gauche0,0) != -1)&&(face_voisins(face_gauche0,1) != -1))
1602 {
1603 face_gauche0 = elem_faces(elem_gauche0, ori+3);
1604 }
1605 face_droite0 = elem_faces(elem_droite0, ori);
1606 if ((face_voisins(face_droite0,0) != -1)&&(face_voisins(face_droite0,1) != -1))
1607 {
1608 face_droite0 = elem_faces(elem_droite0, ori+3);
1609 }
1610
1611 face_gauche1 = elem_faces(elem_gauche1, ori);
1612 if ((face_voisins(face_gauche1,0) != -1)&&(face_voisins(face_gauche1,1) != -1))
1613 {
1614 face_gauche1 = elem_faces(elem_gauche1, ori+3);
1615 }
1616 face_droite1 = elem_faces(elem_droite1, ori);
1617 if ((face_voisins(face_droite1,0) != -1)&&(face_voisins(face_droite1,1) != -1))
1618 {
1619 face_droite1 = elem_faces(elem_droite1, ori+3);
1620 }
1621
1622 //On evite les coins exterieurs (on evite de calculer des gradients sur des faces non parietales)
1623 if((face_droite1<nfin)&&(face_droite0<nfin)&&(face_gauche1<nfin)&&(face_gauche0<nfin)&&
1624 (ndeb<=face_droite1)&&(ndeb<=face_droite0)&&(ndeb<=face_gauche1)&&(ndeb<=face_gauche0))
1625 {
1626 int i=1;
1627
1628 ////////////////
1629 //Calvul de v(1)
1630 ////////////////
1631
1632
1633 if(eq_vit[corresp[num_face]].get_Unp1(0,1)>0)
1634 d0d0_1 = (eq_vit[corresp[face_droite0]].get_Unp1(0,1)-eq_vit[corresp[num_face]].get_Unp1(0,1))
1635 /domaine_VDF.dist_elem_period(elem,elem_droite0,(ori+1)%3);
1636 else
1637 d0d0_1 = (eq_vit[corresp[num_face]].get_Unp1(0,1)-eq_vit[corresp[face_gauche0]].get_Unp1(0,1))
1638 /domaine_VDF.dist_elem_period(elem,elem_droite0,(ori+1)%3);
1639
1640 if(eq_vit[corresp[num_face]].get_Unp1(0,1)>0)
1641 d1d1_1 = (eq_vit[corresp[face_droite1]].get_Unp1(1,1)-eq_vit[corresp[num_face]].get_Unp1(1,1))
1642 /domaine_VDF.dist_elem_period(elem,elem_droite0,(ori+1)%3);
1643 else
1644 d1d1_1 = (eq_vit[corresp[num_face]].get_Unp1(1,1)-eq_vit[corresp[face_gauche1]].get_Unp1(1,1))
1645 /domaine_VDF.dist_elem_period(elem,elem_droite0,(ori+1)%3);
1646
1647 v += -(d0d0_1+d1d1_1)
1648 *0.5*(eq_vit[corresp[num_face]].get_yc(1)-eq_vit[corresp[num_face]].get_yc(0));
1649
1650 eq_vit[corresp[num_face]].set_v(i,v);
1651 //fic_v << eq_vit[corresp[num_face]].get_y(i) << " " << -eq_vit[corresp[num_face]].get_v(i) << finl;
1652
1653 ////////////////
1654 //Calvul de v(i)
1655 ////////////////
1656 for(i=2 ; i<nb_pts ; i++)
1657 {
1658 v += -(
1659 (eq_vit[corresp[face_droite0]].get_Unp1(0,i)-eq_vit[corresp[face_gauche0]].get_Unp1(0,i))
1660 /(domaine_VDF.dist_elem_period(elem_gauche0,elem,(ori+1)%3)+domaine_VDF.dist_elem_period(elem,elem_droite0,(ori+1)%3))
1661 +
1662 (eq_vit[corresp[face_droite1]].get_Unp1(1,i)-eq_vit[corresp[face_gauche1]].get_Unp1(1,i))
1663 /(domaine_VDF.dist_elem_period(elem_gauche1,elem,(ori+2)%3)+domaine_VDF.dist_elem_period(elem,elem_droite1,(ori+2)%3))
1664 )
1665 *0.5*(eq_vit[corresp[num_face]].get_yc(i)-eq_vit[corresp[num_face]].get_yc(i-1))
1666 -(
1667 (eq_vit[corresp[face_droite0]].get_Unp1(0,i-1)-eq_vit[corresp[face_gauche0]].get_Unp1(0,i-1))
1668 /(domaine_VDF.dist_elem_period(elem_gauche0,elem,(ori+1)%3)+domaine_VDF.dist_elem_period(elem,elem_droite0,(ori+1)%3))
1669 +
1670 (eq_vit[corresp[face_droite1]].get_Unp1(1,i-1)-eq_vit[corresp[face_gauche1]].get_Unp1(1,i-1))
1671 /(domaine_VDF.dist_elem_period(elem_gauche1,elem,(ori+2)%3)+domaine_VDF.dist_elem_period(elem,elem_droite1,(ori+2)%3))
1672 )
1673 *0.5*(eq_vit[corresp[num_face]].get_yc(i-1)-eq_vit[corresp[num_face]].get_yc(i-2));
1674
1675 eq_vit[corresp[num_face]].set_v(i,v);
1676 //fic_v << eq_vit[corresp[num_face]].get_y(i) << " " << -eq_vit[corresp[num_face]].get_v(i) << finl;
1677 }
1678
1679 ////////////////
1680 //Calvul de v(N)
1681 ////////////////
1682 //Cerr << "nb_pts = " << nb_pts << finl;
1683 i=nb_pts;
1684
1685 v += -(
1686 (eq_vit[corresp[face_droite0]].get_Unp1(0,i)-eq_vit[corresp[face_gauche0]].get_Unp1(0,i))
1687 /(domaine_VDF.dist_elem_period(elem_gauche0,elem,(ori+1)%3)+domaine_VDF.dist_elem_period(elem,elem_droite0,(ori+1)%3))
1688 +
1689 (eq_vit[corresp[face_droite1]].get_Unp1(1,i)-eq_vit[corresp[face_gauche1]].get_Unp1(1,i))
1690 /(domaine_VDF.dist_elem_period(elem_gauche1,elem,(ori+2)%3)+domaine_VDF.dist_elem_period(elem,elem_droite1,(ori+2)%3))
1691 )
1692 *0.5*(eq_vit[corresp[num_face]].get_yc(nb_pts-1)-eq_vit[corresp[num_face]].get_yc(nb_pts-2));
1693
1694 eq_vit[corresp[num_face]].set_v(nb_pts,v);
1695
1696
1697 // //////////////////////////////////////////////////////////////////////////////////////
1698 // //Fin calcul de la vitesse normale a la paroi a partir de la conservation de la masse
1699 // //////////////////////////////////////////////////////////////////////////////////////
1700
1701 // ///////////////////////////////
1702 // //Calcul des termes convectifs
1703 // ///////////////////////////////
1704
1705 eq_vit[corresp[num_face]].set_F(0, 0, - gradient_de_pression0 + ts0);
1706 eq_vit[corresp[num_face]].set_F(1, 0, - gradient_de_pression1 + ts1);
1707
1708 for(i=1 ; i<nb_pts ; i++)
1709 {
1710 d00d0 = (eq_vit[corresp[face_droite0]].get_Unp1(0,i)*eq_vit[corresp[face_droite0]].get_Unp1(0,i)
1711 -eq_vit[corresp[face_gauche0]].get_Unp1(0,i)*eq_vit[corresp[face_gauche0]].get_Unp1(0,i))
1712 /(domaine_VDF.dist_elem_period(elem_gauche0,elem,(ori+1)%3)
1713 +domaine_VDF.dist_elem_period(elem,elem_droite0,(ori+1)%3));
1714
1715 d01d1 = (eq_vit[corresp[face_droite1]].get_Unp1(0,i)*eq_vit[corresp[face_droite1]].get_Unp1(1,i)
1716 -eq_vit[corresp[face_gauche1]].get_Unp1(0,i)*eq_vit[corresp[face_gauche1]].get_Unp1(1,i))
1717 /(domaine_VDF.dist_elem_period(elem_gauche1,elem,(ori+2)%3)
1718 +domaine_VDF.dist_elem_period(elem,elem_droite1,(ori+2)%3));
1719
1720 d10d0 = (eq_vit[corresp[face_droite0]].get_Unp1(0,i)*eq_vit[corresp[face_droite0]].get_Unp1(1,i)
1721 -eq_vit[corresp[face_gauche0]].get_Unp1(0,i)*eq_vit[corresp[face_gauche0]].get_Unp1(1,i))
1722 /(domaine_VDF.dist_elem_period(elem_gauche0,elem,(ori+1)%3)
1723 +domaine_VDF.dist_elem_period(elem,elem_droite0,(ori+1)%3));
1724
1725 d11d1 = (eq_vit[corresp[face_droite1]].get_Unp1(1,i)*eq_vit[corresp[face_droite1]].get_Unp1(1,i)
1726 -eq_vit[corresp[face_gauche1]].get_Unp1(1,i)*eq_vit[corresp[face_gauche1]].get_Unp1(1,i))
1727 /(domaine_VDF.dist_elem_period(elem_gauche1,elem,(ori+2)%3)
1728 +domaine_VDF.dist_elem_period(elem,elem_droite1,(ori+2)%3));
1729
1730 d0vdy = (eq_vit[corresp[num_face]].get_Unp1(0,i+1)*eq_vit[corresp[num_face]].get_v(i+1)
1731 -eq_vit[corresp[num_face]].get_Unp1(0,i-1)*eq_vit[corresp[num_face]].get_v(i-1))
1732 /(eq_vit[corresp[num_face]].get_y(i+1)-eq_vit[corresp[num_face]].get_y(i-1));
1733
1734 d1vdy = (eq_vit[corresp[num_face]].get_Unp1(1,i+1)*eq_vit[corresp[num_face]].get_v(i+1)
1735 -eq_vit[corresp[num_face]].get_Unp1(1,i-1)*eq_vit[corresp[num_face]].get_v(i-1))
1736 /(eq_vit[corresp[num_face]].get_y(i+1)-eq_vit[corresp[num_face]].get_y(i-1));
1737
1738 eq_vit[corresp[num_face]].set_F(0, i, - gradient_de_pression0 + ts0 - (d00d0+d01d1+d0vdy));
1739 eq_vit[corresp[num_face]].set_F(1, i, - gradient_de_pression1 + ts1 - (d11d1+d10d0+d1vdy));
1740
1741 }
1742
1743 i=nb_pts;
1744
1745 d00d0 = (eq_vit[corresp[face_droite0]].get_Unp1(0,i)*eq_vit[corresp[face_droite0]].get_Unp1(0,i)
1746 -eq_vit[corresp[face_gauche0]].get_Unp1(0,i)*eq_vit[corresp[face_gauche0]].get_Unp1(0,i))
1747 /(domaine_VDF.dist_elem_period(elem_gauche0,elem,(ori+1)%3)
1748 +domaine_VDF.dist_elem_period(elem,elem_droite0,(ori+1)%3));
1749
1750 d01d1 = (eq_vit[corresp[face_droite1]].get_Unp1(0,i)*eq_vit[corresp[face_droite1]].get_Unp1(1,i)
1751 -eq_vit[corresp[face_gauche1]].get_Unp1(0,i)*eq_vit[corresp[face_gauche1]].get_Unp1(1,i))
1752 /(domaine_VDF.dist_elem_period(elem_gauche1,elem,(ori+2)%3)
1753 +domaine_VDF.dist_elem_period(elem,elem_droite1,(ori+2)%3));
1754
1755 d10d0 = (eq_vit[corresp[face_droite0]].get_Unp1(0,i)*eq_vit[corresp[face_droite0]].get_Unp1(1,i)
1756 -eq_vit[corresp[face_gauche0]].get_Unp1(0,i)*eq_vit[corresp[face_gauche0]].get_Unp1(1,i))
1757 /(domaine_VDF.dist_elem_period(elem_gauche0,elem,(ori+1)%3)
1758 +domaine_VDF.dist_elem_period(elem,elem_droite0,(ori+1)%3));
1759
1760 d11d1 = (eq_vit[corresp[face_droite1]].get_Unp1(1,i)*eq_vit[corresp[face_droite1]].get_Unp1(1,i)
1761 -eq_vit[corresp[face_gauche1]].get_Unp1(1,i)*eq_vit[corresp[face_gauche1]].get_Unp1(1,i))
1762 /(domaine_VDF.dist_elem_period(elem_gauche1,elem,(ori+2)%3)
1763 +domaine_VDF.dist_elem_period(elem,elem_droite1,(ori+2)%3));
1764
1765 d0vdy = (eq_vit[corresp[num_face]].get_Unp1(0,i)*eq_vit[corresp[num_face]].get_v(i)
1766 -eq_vit[corresp[num_face]].get_Unp1(0,i-1)*eq_vit[corresp[num_face]].get_v(i-1))
1767 /(eq_vit[corresp[num_face]].get_y(i)-eq_vit[corresp[num_face]].get_y(i-1));
1768
1769 d1vdy = (eq_vit[corresp[num_face]].get_Unp1(1,i)*eq_vit[corresp[num_face]].get_v(i)
1770 -eq_vit[corresp[num_face]].get_Unp1(1,i-1)*eq_vit[corresp[num_face]].get_v(i-1))
1771 /(eq_vit[corresp[num_face]].get_y(i)-eq_vit[corresp[num_face]].get_y(i-1));
1772
1773 // //////////////////////////////////////
1774 // //Fin du calcul des termes convectifs
1775 // //////////////////////////////////////
1776
1777 eq_vit[corresp[num_face]].set_F(0, i, - gradient_de_pression0 + ts0 - (d00d0+d01d1+d0vdy));
1778 eq_vit[corresp[num_face]].set_F(1, i, - gradient_de_pression1 + ts1 - (d11d1+d10d0+d1vdy));
1779
1780 //Cerr << "je suis dans les termes convectifs vitesse" << finl;
1781 }//fin if pas de coins exterieurs
1782 }//fin if pas de coins interieurs
1783
1784 else
1785 {
1786 for(int i=0 ; i<nb_pts+1 ; i++)
1787 {
1788 eq_vit[corresp[num_face]].set_F(0, i, - gradient_de_pression0 + ts0);
1789 eq_vit[corresp[num_face]].set_F(1, i, - gradient_de_pression1 + ts1);
1790 }
1791 }
1792}
classe Champ_Don_base classe de base des Champs donnes (non calcules)
DoubleTab & valeurs() override
Surcharge Champ_base::valeurs() Renvoie le tableau des valeurs.
DoubleTab & valeurs() override
Renvoie le tableau des valeurs du champ au temps courant.
classe Champ_Uniforme Represente un champ constant dans l'espace et dans le temps.
double temps() const
Renvoie le temps du champ.
classe Cond_lim Classe generique servant a representer n'importe quelle classe
Definition Cond_lim.h:31
void setKappa(double k)
Classe Diffu_totale_hyd_base Classe abstraite calculant la diffusivite totale (somme diffusivite.
classe Dirichlet_paroi_fixe Represente une paroi immobile dans une equation de type Navier_Stokes.
class Domaine_VDF
Definition Domaine_VDF.h:64
double dist_norm_bord_axi(int num_face) const
int orientation(int) const override
inline DoubleVect& Domaine_VDF::porosite_face() {
double dist_elem_period(int, int, int) const
double face_normales(int, int) const override
double dist_norm_bord(int num_face) const override
virtual const DoubleVect & face_surfaces() const
Definition Domaine_VF.h:51
int nb_faces() const
renvoie le nombre global de faces.
Definition Domaine_VF.h:471
DoubleVect & volumes_entrelaces()
Definition Domaine_VF.h:99
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
int face_voisins(int num_face, int i) const
renvoie l'element voisin de numface dans la direction i.
Definition Domaine_VF.h:418
int nb_front_Cl() const
Sortie & syncfile() override
Provoque l'ecriture sur disque des donnees accumulees sur les differents processeurs depuis le dernie...
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
classe Equation_base Le role d'une equation est le calcul d'un ou plusieurs champs....
virtual const Milieu_base & milieu() const =0
virtual const RefObjU & get_modele(Type_modele type) const
Sources & sources()
Renvoie les termes sources asssocies a l'equation.
virtual const Champ_Inc_base & inconnue() const =0
Probleme_base & probleme()
Renvoie le probleme associe a l'equation.
Schema_Temps_base & schema_temps()
Renvoie le schema en temps associe a l'equation.
classe Fluide_base Cette classe represente un d'un fluide incompressible ainsi que
Definition Fluide_base.h:38
const Champ_Don_base & viscosite_cinematique() const
Definition Fluide_base.h:58
class Front_VF
Definition Front_VF.h:36
int nb_faces() const
Definition Front_VF.h:53
int num_premiere_face() const
Definition Front_VF.h:63
virtual const Champ_Don_base & beta_t() const
Renvoie beta_t du milieu.
virtual const Champ_Don_base & gravite() const
Renvoie la gravite du milieu si elle a ete associe provoque une erreur sinon.
Classe Modele_turbulence_scal_base Cette classe represente un modele de turbulence pour une equation ...
const Turbulence_paroi_scal_base & loi_paroi() const
Renvoie la loi de turbulence sur la paroi (version const).
Une chaine de caractere (Nom) en majuscules.
Definition Motcle.h:26
classe Navier_Stokes_std Cette classe porte les termes de l'equation de la dynamique
Operateur_Grad & operateur_gradient()
Renvoie l'operateur de calcul du gradient associe a l'equation.
Champ_Inc_base & pression()
class Nom Une chaine de caractere pour nommer les objets de TRUST
Definition Nom.h:31
virtual void set_param(Param &) const
Definition Objet_U.h:135
friend class Entree
Definition Objet_U.h:76
static int dimension
Definition Objet_U.h:99
friend class Sortie
Definition Objet_U.h:75
const Nom & que_suis_je() const
renvoie la chaine identifiant la classe.
Definition Objet_U.cpp:104
virtual Entree & readOn(Entree &)
Lecture d'un Objet_U sur un flot d'entree Methode a surcharger.
Definition Objet_U.cpp:293
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
static int axi
Definition Objet_U.h:101
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
Definition Objet_U.cpp:282
Classe Operateur_Grad Classe generique de la hierarchie des operateurs calculant le gradient.
Helper class to factorize the readOn method of Objet_U classes.
Definition Param.h:112
void ajouter_flag(const char *keyword, const bool *value)
Register a boolean flag whose mere presence switches it to true.
Definition Param.cpp:474
void ajouter_non_std(const char *keyword, const Objet_U *value, Param::Nature nat=Param::OPTIONAL)
Register a keyword handled by Objet_U::lire_motcle_non_standard.
Definition Param.cpp:489
CLASS: ParoiVDF_TBLE.
void set_param(Param &param) const override
int reprendre(Entree &) override
Reprise d'un Objet_U sur un flot d'entree Methode a surcharger.
int calculer_hyd_BiK(DoubleTab &, DoubleTab &) override
void imprimer_ustar(Sortie &) const override
int lire_motcle_non_standard(const Motcle &, Entree &) override
Lecture des parametres de type non simple d'un objet_U a partir d'un flot d'entree.
int init_lois_paroi() override
int sauvegarder(Sortie &) const override
Sauvegarde d'un Objet_U sur un flot de sortie Methode a surcharger.
int calculer_hyd(DoubleTab &) override
const Probleme_base & getPbBase() const override
Eq_couch_lim & get_eq_couche_lim_T(int i)
double tps_start_stat_nu_t_dyn
DoubleTab visco_turb_moy
void calculer_stat(int indice_post, int indice_maillage, double Fx, double Fy, double Fz, double u, double v, double w, double dt)
void imprimer_stat(Sortie &, double) const
void set_param(Param &param) const
int reprendre(Entree &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, double tps)
IntTab num_global_faces_post
int sauvegarder(Sortie &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, double tps) const
DoubleTab valeurs_reprises
int lire_motcle_non_standard(const Motcle &, Entree &)
int init_lois_paroi(const Domaine_VF &, const Domaine_Cl_dis_base &)
CLASS: Paroi_ODVM_scal_VDF.
void set_param(Param &param) const
double calcul_lm_plus(double d_plus)
classe Probleme_base C'est un Probleme_U qui n'est pas un couplage.
virtual const Equation_base & equation(int) const =0
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
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
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
double temps_courant() const
Renvoie le temps courant.
double pas_temps_min() const
Renvoie le pas de temps minimum.
double pas_de_temps() const
Renvoie le pas de temps (delta_t) courant.
Classe de base des flux de sortie.
Definition Sortie.h:52
classe Source_base Un objet Source_base est un terme apparaissant au second membre d'une
Definition Source_base.h:42
virtual DoubleTab & ajouter(DoubleTab &) const
class Sources Sources represente une liste de Source.
Definition Sources.h:31
void resize(_SIZE_ n, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
Definition TRUSTTab.tpp:469
_TYPE_ local_min_vect(Mp_vect_options opt=VECT_REAL_ITEMS) const
Definition TRUSTVect.h:155
const Objet_U & valeur() const
Definition TRUST_Ref.h:134
class Terme_Boussinesq_scalaire_VDF_Face
Classe Terme_Boussinesq_base Cette classe represente le terme de gravite qui figure dans l'equation.
double Scalaire0(int i) const
const DoubleVect & tab_d_plus() const
const DoubleVect & tab_u_star() const
Classe Turbulence_paroi_scal_base Classe de base pour la hierarchie des classes representant les mode...