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