TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
ParoiVEF_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 <ParoiVEF_TBLE.h>
18#include <Paroi_std_hyd_VEF.h>
19#include <Champ_Q1NC.h>
20#include <Dirichlet_paroi_fixe.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 <time.h>
27#include <Modele_turbulence_scal_base.h>
28#include <Probleme_base.h>
29#include <ParoiVEF_TBLE_scal.h>
30#include <Terme_Boussinesq_VEFPreP1B_Face.h>
31#include <Terme_Source_Qdm_lambdaup_VEF_Face.h>
32#include <SFichier.h>
33#include <EcrFicPartage.h>
34#include <Modele_turbulence_hyd_base.h>
35#include <Navier_Stokes_std.h>
36#include <Param.h>
37
38Implemente_instanciable_sans_constructeur(ParoiVEF_TBLE,"Paroi_TBLE_VEF",Paroi_hyd_base_VEF);
39
40ParoiVEF_TBLE::ParoiVEF_TBLE()
41{
42 alpha=1.;
43}
44
45// printOn()
46Sortie& ParoiVEF_TBLE::printOn(Sortie& os) const
47{
48 return os << que_suis_je() << " " << le_nom();
49}
50
51//// readOn
52
53Entree& ParoiVEF_TBLE::readOn(Entree& is)
54{
55 Cerr<<"Reading of data for a "<<que_suis_je()<<" wall law"<<finl;
56 Param param(que_suis_je());
57 set_param(param);
58 param.lire_avec_accolades_depuis(is);
59 return is;
60}
61
62void ParoiVEF_TBLE::set_param(Param& param) const
63{
64 Paroi_hyd_base_VEF::set_param(param);
65 Paroi_TBLE_QDM::set_param(param);
66 Paroi_log_QDM::set_param(param);
67 param.ajouter("alpha",&alpha);
68 param.ajouter_condition("(value_of_alpha_ge_0)_and_(value_of_alpha_le_1)","The following condition must be satisfied : 0 <= alpha <= 1");
69 //Ajout pour la lecture non standard
70 //Sinon probleme dans Paroi_TBLE_QDM pour (this) qui est annonce const
71 param.ajouter_non_std("nu_t_dyn",(this));
72 param.ajouter_non_std("tps_start_stat_nu_t_dyn",(this));
73 param.ajouter_non_std("tps_nu_t_dyn",(this));
74 param.ajouter_non_std("sonde_tble",(this));
75 param.ajouter_non_std("stationnaire",(this));
76 param.ajouter_non_std("mu",(this));
77 param.ajouter_non_std("lambda",(this));
78 param.ajouter_non_std("sans_source_boussinesq",(this));
79}
80
85
86/////////////////////////////////////////////////////////////////////
87//
88// Implementation des fonctions de la classe ParoiVEF_TBLE
89//
90/////////////////////////////////////////////////////////////////////
91
92// /////////////////////////////////////////////////
93// // Initialisation des tableaux
94// ////////////////////////////////////////////////
95
96int ParoiVEF_TBLE::init_lois_paroi()
97{
98 if(je_suis_maitre())
99 {
100 Cerr << "ParoiVEF_TBLE::init_lois_paroi()" << finl;
101 }
102
103 const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF, le_dom_dis_.valeur());
104 const IntTab& face_voisins = domaine_VEF.face_voisins();
105 const IntTab& elem_faces = domaine_VEF.elem_faces();
106 const Domaine& domaine = domaine_VEF.domaine();
107 const DoubleTab& face_normale = domaine_VEF.face_normales();
108 const int nfac = domaine.nb_faces_elem();
109 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
110 const Fluide_base& le_fluide = ref_cast(Fluide_base,eqn_hydr.milieu());
111 const DoubleTab& vit = eqn_hydr.inconnue().valeurs();
112
113 // OC: Je ne comprends pas pourquoi dans les autres LP VEF on resize a nb_faces_tot le cisaillement et le tab_u_star
114 // De plus aucune loi de paroi VEF n appelle la methode init_lois_paroi_ de la classe mere
115 // qui devrait elle se charger de dimensionner ce tableau
116 // Pour l'instant, ici, on appele bien l'init de la classe mere qui resize bien a nb_faces_bord.
117
118 Paroi_hyd_base_VEF::init_lois_paroi_();
119
120 Paroi_TBLE_QDM::init_lois_paroi(le_dom_dis_, le_dom_Cl_dis_.valeur());
121
122 int elem;
123
124 int compteur_faces_paroi = 0;
125
126 double surf2;
127 double dist; //distance du centre de la maille a la paroi
128
129 IntVect num(nfac);
130 DoubleVect n(dimension),t1(dimension),t2(dimension); // vecteurs orthonomes du repere local associe a la face paroi
131 DoubleVect v_tang(nb_comp);
132
133 for (int n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
134 {
135 // pour chaque condition limite on regarde son type
136 // On applique les lois de paroi uniquement
137 // aux voisinages des parois
138
139 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
140
141 if (sub_type(Dirichlet_paroi_fixe,la_cl.valeur()) )
142 {
143 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
144 int size=le_bord.nb_faces();
145
146 //Boucle sur les faces des bords parietaux
147 for (int ind_face=0; ind_face<size; ind_face++)
148 {
149 Eq_couch_lim& equation_vitesse = eq_vit[compteur_faces_paroi];
150 Diffu_totale_hyd_base& diffu_hyd = ref_cast_non_const(Diffu_totale_hyd_base, equation_vitesse.get_diffu()); //modele de viscosite turbulente
151 diffu_hyd.setKappa(Kappa_);
152
153 int num_face = le_bord.num_face(ind_face);
154 elem = face_voisins(num_face,0);
155
156 surf2=0.;
157
158 for(int i=0; i<dimension; i++)
159 {
160 surf2 += face_normale(num_face,i)*face_normale(num_face,i);
161 }
162
163 if (dimension == 2)
164 {
165 num[0]=elem_faces(elem,0);
166 num[1]=elem_faces(elem,1);
167
168 if (num[0]==num_face) num[0]=elem_faces(elem,2);
169 else if (num[1]==num_face) num[1]=elem_faces(elem,2);
170
171 n[0] = face_normale(num_face,0)/sqrt(surf2);
172 n[1] = face_normale(num_face,1)/sqrt(surf2);
173
174 t1[0] = -n[1];
175 t1[1] = n[0];
176
177 dist = distance_2D(num_face,elem,domaine_VEF)*3./2.;
178
179 v_tang[0] = ((vit(num[0],0)+vit(num[1],0))*t1[0]
180 +(vit(num[0],1)+vit(num[1],1))*t1[1])/2.;
181 }
182 else
183 {
184 num[0]=elem_faces(elem,0);
185 num[1]=elem_faces(elem,1);
186 num[2]=elem_faces(elem,2);
187
188 if (num[0]==num_face) num[0]=elem_faces(elem,3);
189 else if (num[1]==num_face) num[1]=elem_faces(elem,3);
190 else if (num[2]==num_face) num[2]=elem_faces(elem,3);
191
192 n[0] = face_normale(num_face,0)/sqrt(surf2);
193 n[1] = face_normale(num_face,1)/sqrt(surf2);
194 n[2] = face_normale(num_face,2)/sqrt(surf2);
195
196 if( est_egal(n[0],0.) && est_egal(n[1],0.) )
197 {
198 t1[0] = 0.;
199 t1[1] = 1.;
200 t1[2] = 0.;;
201 }
202 else
203 {
204 t1[0] = -n[1]/sqrt(n[0]*n[0]+n[1]*n[1]);
205 t1[1] = n[0]/sqrt(n[0]*n[0]+n[1]*n[1]);
206 t1[2] = 0.;
207 }
208
209 t2[0] = n[1]*t1[2] - n[2]*t1[1];
210 t2[1] = n[2]*t1[0] - n[0]*t1[2];
211 t2[2] = n[0]*t1[1] - n[1]*t1[0];
212
213 dist = distance_3D(num_face,elem,domaine_VEF)*4./3.;
214
215 v_tang[0] = ((vit(num[0],0)+vit(num[1],0)+vit(num[2],0))*t1[0]
216 +(vit(num[0],1)+vit(num[1],1)+vit(num[2],1))*t1[1]
217 +(vit(num[0],2)+vit(num[1],2)+vit(num[2],2))*t1[2])/3.;
218
219 v_tang[1] = ((vit(num[0],0)+vit(num[1],0)+vit(num[2],0))*t2[0]
220 +(vit(num[0],1)+vit(num[1],1)+vit(num[2],1))*t2[1]
221 +(vit(num[0],2)+vit(num[1],2)+vit(num[2],2))*t2[2])/3.;
222 }
223
224
225 equation_vitesse.set_y0(0.); //ordonnee de la paroi
226 equation_vitesse.set_yn(dist); //ordonnee du 1er centre de maille
227
228 equation_vitesse.initialiser(nb_comp, nb_pts, fac, epsilon, max_it, nu_t_dyn); //nbre de pts maillage fin
229
230
231
232 for (int i=0; i<nb_comp; i++)
233 {
234 equation_vitesse.set_u_y0(i,0.); // vitesse a la paroi
235 equation_vitesse.set_u_yn(i,v_tang[i]); // vitesse en yn
236
237 //vitesse sur le maillage fin a l'instant initial
238
239 if (reprise_ok==0)
240 equation_vitesse.set_Uinit_lin(i, 0., v_tang[i]);
241 else
242 {
243 for (int itble=0; itble<nb_pts+1; itble++)
244 equation_vitesse.set_Uinit(i,itble,valeurs_reprises(compteur_faces_paroi, i, itble)) ;
245 }
246 }
247 compteur_faces_paroi++;
248
249 }//Fin boucle sur num_face
250 }//Fin if Paroi_fixe
251 }//Fin boucle sur les bords parietaux
252
253 f_tang_moy.resize(compteur_faces_paroi, nb_comp);
254 f_tang_moy=0.;
255
256 // Si on n'a pas interdit la presence du terme de Boussi dans TBLE (mot cle "sans_source_boussinesq",
257 // Alors on verifie la presence ou non d'un terme source de Boussinesq dans le pb hydraulique :
258 if (source_boussinesq==1)
259 {
260 source_boussinesq=0;
261 const Sources& les_sources=eqn_hydr.sources();
262 int nb_sources=les_sources.size();
263
264 for (int j=0; j<nb_sources; j++)
265 {
266 const Source_base& ts = les_sources(j).valeur();
267 if (sub_type(Terme_Boussinesq_base,ts))
268 {
269 const Terme_Boussinesq_base& terme_boussi = ref_cast(Terme_Boussinesq_base,ts);
270 T0 = terme_boussi.Scalaire0(0);
271 source_boussinesq=1;
272 const Champ_Don_base& ch_beta_t = le_fluide.beta_t();
273 const DoubleTab& tab_champ_beta_t = ch_beta_t.valeurs();
274 if (sub_type(Champ_Uniforme,ch_beta_t))
275 {
276 beta_t = std::max(tab_champ_beta_t(0,0),DMINFLOAT);
277 }
278 else
279 {
280 Cerr << " On ne sait pas traiter un beta_t non uniforme dans TBLE !!!"<< finl;
281 Process::exit();
282 }
283
284 break;
285 }
286 }
287 }
288 Cerr << "fin init_lois_paroi()" << finl;
289 return 1;
290}
291
292
293/////////////////////////////////////////////////
294//Calcul du cisaillement a la paroi
295////////////////////////////////////////////////
296
297int ParoiVEF_TBLE::calculer_hyd(DoubleTab& nu_t, DoubleTab& tab_k)
298{
299 return calculer_hyd(nu_t,0,tab_k);
300}
301
302int ParoiVEF_TBLE::calculer_hyd(DoubleTab& tab_k_eps)
303{
304
305 DoubleTab bidon(0);
306 // bidon ne servira pas
307 return calculer_hyd(tab_k_eps,1,bidon);
308
309}
310
311int ParoiVEF_TBLE::calculer_hyd_BiK(DoubleTab& tab_k,DoubleTab& tab_eps)
312{
313 // on est en K-eps
314 const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF, le_dom_dis_.valeur());
315 const IntTab& face_voisins = domaine_VEF.face_voisins();
316 const IntTab& elem_faces = domaine_VEF.elem_faces();
317 const Domaine& domaine = domaine_VEF.domaine();
318 const DoubleTab& face_normale = domaine_VEF.face_normales();
319 const int nfac = domaine.nb_faces_elem();
320 const DoubleVect& volumes_entrelaces = domaine_VEF.volumes_entrelaces();
321 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
322
323 const Fluide_base& le_fluide = ref_cast(Fluide_base,eqn_hydr.milieu());
324 const Champ_Don_base& ch_visco_cin = le_fluide.viscosite_cinematique();
325 const DoubleTab& tab_visco = ch_visco_cin.valeurs();
326
327 if (source_boussinesq==1)
328 {
329 const Champ_Don_base& ch_gravite=le_fluide.gravite();
330
331 if (!sub_type(Champ_Uniforme,ch_gravite))
332 {
333 Cerr << " On ne sait pas traiter la gravite non uniforme dans TBLE !!!"<< finl;
334 Process::exit();
335 }
336 }
337
338 int methode=-1;
339 int is_champ_Q1NC=sub_type(Champ_Q1NC,eqn_hydr.inconnue());
340 remplir_face_keps_imposee( flag_face_keps_imposee_, methode, face_keps_imposee_, domaine_VEF,le_dom_Cl_dis_,!is_champ_Q1NC);
341
342 double visco=-1;
343 int l_unif;
344 if (sub_type(Champ_Uniforme,ch_visco_cin))
345 {
346 visco = std::max(tab_visco(0,0),DMINFLOAT);
347 l_unif = 1;
348 }
349 else
350 l_unif = 0;
351
352 // preparer_calcul_hyd(tab);
353 if ((!l_unif) && (tab_visco.local_min_vect()<DMINFLOAT))
354 // on ne doit pas changer tab_visco ici !
355 {
356 Cerr << "In ParoiVEF_TBLE::calculer_hyd_BiK : visco = " << tab_visco.local_min_vect() << " <= 0 ? " << finl;
357 throw;
358 }
359 //tab_visco+=DMINFLOAT;
360
361 int elem;
362
363 int itmax=0,itmax_loc;
364
365 const DoubleTab& vit = eqn_hydr.inconnue().valeurs(); //vitesse
366
367 const Navier_Stokes_std& eqnNS = ref_cast(Navier_Stokes_std, eqn_hydr);
368
369 // Calcul du gradient de pression
370 DoubleTab grad_p(vit);
371 const DoubleTab& p = eqnNS.pression().valeurs();
372 const Operateur_Grad& gradient = eqnNS.operateur_gradient();
373 gradient.calculer(p, grad_p);
374
375 const DoubleTab& visco_turb = mon_modele_turb_hyd->viscosite_turbulente().valeurs();
376
377
378 DoubleTab termes_sources(vit);
379 termes_sources=0;
380 // On calcule les termes sources, sauf celui de lambda_uprime et de Boussinesq (TBLE recalcule par lui meme ce terme s'il est demande)
381
382 const double tps = eqnNS.schema_temps().temps_courant();
383 const double dt = eqnNS.schema_temps().pas_de_temps();
384 const double dt_min = eqnNS.schema_temps().pas_temps_min();
385
386 int compteur_faces_paroi = 0;
387
388 double surf2;
389 double dist; //distance du centre de la maille a la paroi
390
391 IntVect num(nfac);
392 DoubleVect n(dimension), t1(dimension), t2(dimension); // vecteurs orthonomes du repere local associe a la face paroi
393 DoubleVect F(dimension);
394 DoubleVect v_tang(nb_comp),F_tang(nb_comp);
395 DoubleTab ts_boussi(nb_comp, nb_pts+1);
396 ts_boussi=0.;
397
398 for (int n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
399 {
400 // pour chaque condition limite on regarde son type
401 // On applique les lois de paroi uniquement
402 // aux voisinages des parois
403
404 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
405
406 if (sub_type(Dirichlet_paroi_fixe,la_cl.valeur()) )
407 {
408 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
409 int size=le_bord.nb_faces();
410
411 //Boucle sur les faces des bords parietaux
412 for (int ind_face=0; ind_face<size; ind_face++)
413 {
414 Eq_couch_lim& equation_vitesse = eq_vit[compteur_faces_paroi];
415 int num_face = le_bord.num_face(ind_face);
416 elem = face_voisins(num_face,0);
417
418 surf2=0.;
419
420 for(int i=0; i<dimension; i++)
421 {
422 surf2 += face_normale(num_face,i)*face_normale(num_face,i);
423 }
424
425 if (dimension == 2)
426 {
427 num[0]=elem_faces(elem,0);
428 num[1]=elem_faces(elem,1);
429
430 if (num[0]==num_face) num[0]=elem_faces(elem,2);
431 else if (num[1]==num_face) num[1]=elem_faces(elem,2);
432
433 n[0] = face_normale(num_face,0)/sqrt(surf2);
434 n[1] = face_normale(num_face,1)/sqrt(surf2);
435
436 t1[0] = -n[1];
437 t1[1] = n[0];
438
439 dist = distance_2D(num_face,elem,domaine_VEF)*3./2.;
440
441 v_tang[0] = ((vit(num[0],0)+vit(num[1],0))*t1[0]
442 +(vit(num[0],1)+vit(num[1],1))*t1[1])/2.;
443
444 // Terme source de Boussinesq si demande
445 if (source_boussinesq==1)
446 {
447 ts_boussi = 0.;
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 ParoiVEF_TBLE_scal& loi_tble_T = ref_cast_non_const(ParoiVEF_TBLE_scal,loi);
454 if (loi_tble_T.est_initialise())
455 {
456
457 double gt = 0.;
458 for (int idim=0; idim<dimension; idim++)
459 gt += gravite(idim)*t1[idim];
460 for(int i=0; i<nb_pts+1; i++)
461 ts_boussi(0, i) = -gt*beta_t*(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i)-T0);
462 }
463 }
464
465
466 F[0] = ((termes_sources(num[0],0)-grad_p(num[0],0))/volumes_entrelaces(num[0])
467 +(termes_sources(num[1],0)-grad_p(num[1],0))/volumes_entrelaces(num[1]))/2. ;
468
469 F[1] = ((termes_sources(num[0],1)-grad_p(num[0],1))/volumes_entrelaces(num[0])
470 +(termes_sources(num[1],1)-grad_p(num[1],1))/volumes_entrelaces(num[1]))/2. ;
471
472 F_tang[0] = F[0]*t1[0] + F[1]*t1[1] ; // projection de F suivant l'axe tangentiel_1
473 f_tang_moy(compteur_faces_paroi,0) = alpha*F_tang[0] +(1-alpha)*f_tang_moy(compteur_faces_paroi,0);
474 }
475 else
476 {
477 num[0]=elem_faces(elem,0);
478 num[1]=elem_faces(elem,1);
479 num[2]=elem_faces(elem,2);
480
481 if (num[0]==num_face) num[0]=elem_faces(elem,3);
482 else if (num[1]==num_face) num[1]=elem_faces(elem,3);
483 else if (num[2]==num_face) num[2]=elem_faces(elem,3);
484
485 n[0] = face_normale(num_face,0)/sqrt(surf2);
486 n[1] = face_normale(num_face,1)/sqrt(surf2);
487 n[2] = face_normale(num_face,2)/sqrt(surf2);
488
489
490 if( est_egal(n[0],0.) && est_egal(n[1],0.) )
491 {
492 t1[0] = 0.;
493 t1[1] = 1;
494 t1[2] = 0;
495 }
496 else
497 {
498 t1[0] = -n[1]/sqrt(n[0]*n[0]+n[1]*n[1]);
499 t1[1] = n[0]/sqrt(n[0]*n[0]+n[1]*n[1]);
500 t1[2] = 0.;
501 }
502
503 t2[0] = n[1]*t1[2] - n[2]*t1[1];
504 t2[1] = n[2]*t1[0] - n[0]*t1[2];
505 t2[2] = n[0]*t1[1] - n[1]*t1[0];
506
507 dist = distance_3D(num_face,elem,domaine_VEF)*4./3.;
508
509 v_tang[0] = ((vit(num[0],0)+vit(num[1],0)+vit(num[2],0))*t1[0]
510 +(vit(num[0],1)+vit(num[1],1)+vit(num[2],1))*t1[1]
511 +(vit(num[0],2)+vit(num[1],2)+vit(num[2],2))*t1[2])/3.;
512
513 v_tang[1] = ((vit(num[0],0)+vit(num[1],0)+vit(num[2],0))*t2[0]
514 +(vit(num[0],1)+vit(num[1],1)+vit(num[2],1))*t2[1]
515 +(vit(num[0],2)+vit(num[1],2)+vit(num[2],2))*t2[2])/3.;
516
517 F[0] = ((termes_sources(num[0],0)-grad_p(num[0],0))/volumes_entrelaces(num[0])
518 +(termes_sources(num[1],0)-grad_p(num[1],0))/volumes_entrelaces(num[1])
519 +(termes_sources(num[2],0)-grad_p(num[2],0))/volumes_entrelaces(num[2]))/3.;
520
521 F[1] = ((termes_sources(num[0],1)-grad_p(num[0],1))/volumes_entrelaces(num[0])
522 +(termes_sources(num[1],1)-grad_p(num[1],1))/volumes_entrelaces(num[1])
523 +(termes_sources(num[2],1)-grad_p(num[2],1))/volumes_entrelaces(num[2]))/3.;
524
525 F[2] = ((termes_sources(num[0],2)-grad_p(num[0],2))/volumes_entrelaces(num[0])
526 +(termes_sources(num[1],2)-grad_p(num[1],2))/volumes_entrelaces(num[1])
527 +(termes_sources(num[2],2)-grad_p(num[2],2))/volumes_entrelaces(num[2]))/3.;
528
529
530 F_tang[0] = F[0]*t1[0] + F[1]*t1[1] + F[2]*t1[2] ; // projection de F suivant l'axe tangentiel_1
531 F_tang[1] = F[0]*t2[0] + F[1]*t2[1] + F[2]*t2[2] ; // projection de F suivant l'axe tangentiel_2
532 f_tang_moy(compteur_faces_paroi,0) = alpha*F_tang[0] +(1-alpha)*f_tang_moy(compteur_faces_paroi,0);
533 f_tang_moy(compteur_faces_paroi,1) = alpha*F_tang[1] +(1-alpha)*f_tang_moy(compteur_faces_paroi,1);
534
535 // Terme source de Boussinesq si demande
536 if (source_boussinesq==1)
537 {
538 ts_boussi = 0.;
539 const Champ_Don_base& ch_gravite=le_fluide.gravite();
540 const DoubleVect& gravite = ch_gravite.valeurs();
541 const Equation_base& eqn_th = mon_modele_turb_hyd->equation().probleme().equation(1);
542 const Modele_turbulence_scal_base& le_mod_turb_th = ref_cast(Modele_turbulence_scal_base,eqn_th.get_modele(TURBULENCE).valeur());
543 const Turbulence_paroi_scal_base& loi = le_mod_turb_th.loi_paroi();
544 ParoiVEF_TBLE_scal& loi_tble_T = ref_cast_non_const(ParoiVEF_TBLE_scal,loi);
545 if (loi_tble_T.est_initialise())
546 {
547
548 double gt1 = 0.;
549 double gt2 = 0.;
550 for (int idim=0; idim<dimension; idim++)
551 {
552 gt1 += gravite(idim)*t1[idim];
553 gt2 += gravite(idim)*t2[idim];
554 }
555
556 for(int i=0; i<nb_pts+1; i++)
557 {
558 double tmp = beta_t*(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i)-T0);
559 ts_boussi(0, i) = -gt1*tmp;
560 ts_boussi(1, i) = -gt2*tmp;
561 }
562 /*Cout << "gt1 " << gt1 << finl;
563 Cout << "gt2 " << gt2 << finl;
564 Cout << "T0 " << T0 << finl;
565 Cout << "beta_t " << beta_t << finl;*/
566 }
567 }
568 }
569
570
571 if(tps<=dt_min)
572 {
573 visco_turb_moy(elem)=visco_turb(elem);
574 }
575 else if(tps>tps_start_stat_nu_t_dyn)
576 {
577 visco_turb_moy(elem) = (dt/tps)*visco_turb(elem)+(1-(dt/tps))*visco_turb_moy(elem);
578 }
579
580
581 // au premier appel, dt ne vaut pas dt_min mais 0. car un appel a mettre_a_jour dans le preparer_calcul des modeles de turb
582 // precede l'initialisation du schema en temps.
583 if(dt<dt_min) equation_vitesse.set_dt(1.e12);
584 else equation_vitesse.set_dt(dt);
585
586
587 for (int i=0; i<nb_comp; i++)
588 {
589 equation_vitesse.set_u_y0(i,0.); // vitesse a la paroi
590 equation_vitesse.set_u_yn(i,v_tang[i]); // vitesse en yn
591
592 //equation_vitesse.set_F(i, F_tang[i]);
593 for(int ipts=0 ; ipts<nb_pts+1; ipts++)
594 equation_vitesse.set_F(i, ipts, f_tang_moy(compteur_faces_paroi,i)+ts_boussi(i,ipts));
595 }
596
597
598 /////////////////////////////////////////////////////
599 // On resout les equations aux limites simplifiees //
600 // pendant un pas de temps du maillage grossier //
601 /////////////////////////////////////////////////////
602
603
604 if (statio==0)
605 equation_vitesse.aller_au_temps(tps);
606 else
607 equation_vitesse.aller_jusqu_a_convergence(max_it_statio, eps_statio);
608
609
610 itmax_loc = equation_vitesse.get_it();
611
612 if(itmax_loc>max_it) Cerr << "WARNING : TOO MANY ITERATIONS ARE NEEDED IN THE TBLE MESH" << finl;
613
614 if(itmax_loc>itmax) itmax=itmax_loc;
615
616
617 ///////////////////////////////////////////////
618 ///// Determination des composantes //////
619 ///// du cisaillement a la paroi //////
620 ///////////////////////////////////////////////
621
622
623 double d_visco;
624
625 if (l_unif==1)
626 d_visco = visco;
627 else
628 d_visco = tab_visco[elem];
629
630 if (dimension == 2)
631 {
632 double C0 = equation_vitesse.get_cis(0);
633
634 double Cx = t1[0]*C0 ;
635 double Cy = t1[1]*C0 ;
636
637 Cisaillement_paroi_(num_face,0) = Cx;
638 Cisaillement_paroi_(num_face,1) = Cy;
639
640 tab_u_star_(num_face) = pow(Cx*Cx+Cy*Cy,0.25);
641 tab_d_plus_(num_face) = dist*tab_u_star(num_face)/d_visco;
642 }
643 else
644 {
645 double C0 = equation_vitesse.get_cis(0);
646 double C1 = equation_vitesse.get_cis(1);
647
648 double Cx = t1[0]*C0 + t2[0]*C1 ;
649 double Cy = t1[1]*C0 + t2[1]*C1 ;
650 double Cz = t1[2]*C0 + t2[2]*C1 ;
651
652 Cisaillement_paroi_(num_face,0) = Cx;
653 Cisaillement_paroi_(num_face,1) = Cy;
654 Cisaillement_paroi_(num_face,2) = Cz;
655
656 tab_u_star_(num_face) = pow(Cx*Cx+Cy*Cy+Cz*Cz,0.25);
657 tab_d_plus_(num_face) = dist*tab_u_star(num_face)/d_visco;
658 }
659
660
661 //////////////////////////////////////////////////
662 ///// Mise a jour des grandeurs turbulentes //////
663 ///// nu_t - k - eps //////
664 //////////////////////////////////////////////////
665
666
667
668 traitement_keps_BiK(tab_k, tab_eps, num_face, face_voisins, elem_faces, nfac, dist, tab_d_plus(num_face), d_visco, tab_u_star(num_face));
669
670 compteur_faces_paroi++;
671
672 }//Fin boucle sur num_face
673 }//Fin if Paroi_fixe
674 }//Fin boucle sur les bords parietaux
675
676
677 if (Process::je_suis_maitre())
678 {
679 SFichier fic("iter.dat",ios::app); // ouverture du fichier iter.dat
680 fic << tps << " " << itmax << finl;
681 }
682
683 if(oui_stats==1)
684 calculer_stats();
685
686 return 1;
687}
688
689int ParoiVEF_TBLE::calculer_hyd(DoubleTab& tab1,int isKeps,DoubleTab& tab2)
690{
691 // Certains modeles de turbulence appelle cette methode avant le debut du calcul
692 // Ca coince ici alors car le tble_scal n est pas initialise.
693 // pas propre mais on laisse pour l'instant :
694 // si isKeps=1, on est dans le cas k-eps
695 const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF, le_dom_dis_.valeur());
696 const IntTab& face_voisins = domaine_VEF.face_voisins();
697 const IntTab& elem_faces = domaine_VEF.elem_faces();
698 const Domaine& domaine = domaine_VEF.domaine();
699 const DoubleTab& face_normale = domaine_VEF.face_normales();
700 const int nfac = domaine.nb_faces_elem();
701 const DoubleVect& volumes_entrelaces = domaine_VEF.volumes_entrelaces();
702 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
703
704 const Fluide_base& le_fluide = ref_cast(Fluide_base,eqn_hydr.milieu());
705 const Champ_Don_base& ch_visco_cin = le_fluide.viscosite_cinematique();
706 const DoubleTab& tab_visco = ch_visco_cin.valeurs();
707
708 if (source_boussinesq==1)
709 {
710 const Champ_Don_base& ch_gravite=le_fluide.gravite();
711
712 if (!sub_type(Champ_Uniforme,ch_gravite))
713 {
714 Cerr << " On ne sait pas traiter la gravite non uniforme dans TBLE !!!"<< finl;
715 Process::exit();
716 }
717 }
718
719 if (isKeps==1) // on est en K-eps
720 {
721 // on prend la methode par defaut
722 int methode=-1;
723 int is_champ_Q1NC=sub_type(Champ_Q1NC,eqn_hydr.inconnue());
724 remplir_face_keps_imposee( flag_face_keps_imposee_, methode, face_keps_imposee_, domaine_VEF,le_dom_Cl_dis_,!is_champ_Q1NC);
725 }
726
727 double visco=-1;
728 int l_unif;
729 if (sub_type(Champ_Uniforme,ch_visco_cin))
730 {
731 visco = std::max(tab_visco(0,0),DMINFLOAT);
732 l_unif = 1;
733 }
734 else
735 l_unif = 0;
736
737 // preparer_calcul_hyd(tab);
738 if ((!l_unif) && (tab_visco.local_min_vect()<DMINFLOAT))
739 // on ne doit pas changer tab_visco ici !
740 {
741 Cerr << "In ParoiVEF_TBLE::calculer_hyd : visco = " << tab_visco.local_min_vect() << " <= 0 ? " << finl;
742 throw;
743 }
744 //tab_visco+=DMINFLOAT;
745
746 int elem;
747
748 int itmax=0,itmax_loc;
749
750 const DoubleTab& vit = eqn_hydr.inconnue().valeurs(); //vitesse
751
752 const Navier_Stokes_std& eqnNS = ref_cast(Navier_Stokes_std, eqn_hydr);
753
754 // Calcul du gradient de pression
755 DoubleTab grad_p(vit);
756 const DoubleTab& p = eqnNS.pression().valeurs();
757 const Operateur_Grad& gradient = eqnNS.operateur_gradient();
758 gradient.calculer(p, grad_p);
759
760 const DoubleTab& visco_turb = mon_modele_turb_hyd->viscosite_turbulente().valeurs();
761
762
763 DoubleTab termes_sources(vit);
764 termes_sources=0;
765 // On calcule les termes sources, sauf celui de lambda_uprime et de Boussinesq (TBLE recalcule par lui meme ce terme s'il est demande)
766 /* const Sources& les_sources=eqn_hydr.sources();
767 int nb_sources=les_sources.size();
768 for (int j=0; j<nb_sources; j++)
769 {
770 const Source_base& ts = les_sources(j).valeur();
771 if (!sub_type(Terme_Boussinesq_scalaire_VEFPreP1B_Face,ts))
772 {
773 ts.ajouter(termes_sources);
774 }
775 }
776 */
777 // termes_sources.resize(domaine_VEF.nb_faces(),dimension);
778 // eqn_hydr.sources().calculer(termes_sources); //les termes sources : commente : Pour ne pas prendre en compte lambda.uprime dans un premier temps
779
780
781
782 const double tps = eqnNS.schema_temps().temps_courant();
783 const double dt = eqnNS.schema_temps().pas_de_temps();
784 const double dt_min = eqnNS.schema_temps().pas_temps_min();
785
786 int compteur_faces_paroi = 0;
787
788 double surf2;
789 double dist; //distance du centre de la maille a la paroi
790
791 IntVect num(nfac);
792 DoubleVect n(dimension), t1(dimension), t2(dimension); // vecteurs orthonomes du repere local associe a la face paroi
793 DoubleVect F(dimension);
794 DoubleVect v_tang(nb_comp),F_tang(nb_comp);
795 DoubleTab ts_boussi(nb_comp, nb_pts+1);
796 ts_boussi=0.;
797
798 for (int n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
799 {
800 // pour chaque condition limite on regarde son type
801 // On applique les lois de paroi uniquement
802 // aux voisinages des parois
803
804 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
805
806 if (sub_type(Dirichlet_paroi_fixe,la_cl.valeur()) )
807 {
808 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
809 int size=le_bord.nb_faces();
810
811 //Boucle sur les faces des bords parietaux
812 for (int ind_face=0; ind_face<size; ind_face++)
813 {
814 Eq_couch_lim& equation_vitesse = eq_vit[compteur_faces_paroi];
815 int num_face = le_bord.num_face(ind_face);
816 elem = face_voisins(num_face,0);
817
818 surf2=0.;
819
820 for(int i=0; i<dimension; i++)
821 {
822 surf2 += face_normale(num_face,i)*face_normale(num_face,i);
823 }
824
825 if (dimension == 2)
826 {
827 num[0]=elem_faces(elem,0);
828 num[1]=elem_faces(elem,1);
829
830 if (num[0]==num_face) num[0]=elem_faces(elem,2);
831 else if (num[1]==num_face) num[1]=elem_faces(elem,2);
832
833 n[0] = face_normale(num_face,0)/sqrt(surf2);
834 n[1] = face_normale(num_face,1)/sqrt(surf2);
835
836 t1[0] = -n[1];
837 t1[1] = n[0];
838
839 dist = distance_2D(num_face,elem,domaine_VEF)*3./2.;
840
841 v_tang[0] = ((vit(num[0],0)+vit(num[1],0))*t1[0]
842 +(vit(num[0],1)+vit(num[1],1))*t1[1])/2.;
843
844 // Terme source de Boussinesq si demande
845 if (source_boussinesq==1)
846 {
847 ts_boussi = 0.;
848 const Champ_Don_base& ch_gravite=le_fluide.gravite();
849 const DoubleVect& gravite = ch_gravite.valeurs();
850 const Equation_base& eqn_th = mon_modele_turb_hyd->equation().probleme().equation(1);
851 const Modele_turbulence_scal_base& le_mod_turb_th = ref_cast(Modele_turbulence_scal_base,eqn_th.get_modele(TURBULENCE).valeur());
852 const Turbulence_paroi_scal_base& loi = le_mod_turb_th.loi_paroi();
853 ParoiVEF_TBLE_scal& loi_tble_T = ref_cast_non_const(ParoiVEF_TBLE_scal,loi);
854 if (loi_tble_T.est_initialise())
855 {
856
857 double gt = 0.;
858 for (int idim=0; idim<dimension; idim++)
859 gt += gravite(idim)*t1[idim];
860 for(int i=0; i<nb_pts+1; i++)
861 ts_boussi(0, i) = -gt*beta_t*(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i)-T0);
862 }
863 }
864
865
866 F[0] = ((termes_sources(num[0],0)-grad_p(num[0],0))/volumes_entrelaces(num[0])
867 +(termes_sources(num[1],0)-grad_p(num[1],0))/volumes_entrelaces(num[1]))/2. ;
868
869 F[1] = ((termes_sources(num[0],1)-grad_p(num[0],1))/volumes_entrelaces(num[0])
870 +(termes_sources(num[1],1)-grad_p(num[1],1))/volumes_entrelaces(num[1]))/2. ;
871
872 F_tang[0] = F[0]*t1[0] + F[1]*t1[1] ; // projection de F suivant l'axe tangentiel_1
873 f_tang_moy(compteur_faces_paroi,0) = alpha*F_tang[0] +(1-alpha)*f_tang_moy(compteur_faces_paroi,0);
874 }
875 else
876 {
877 num[0]=elem_faces(elem,0);
878 num[1]=elem_faces(elem,1);
879 num[2]=elem_faces(elem,2);
880
881 if (num[0]==num_face) num[0]=elem_faces(elem,3);
882 else if (num[1]==num_face) num[1]=elem_faces(elem,3);
883 else if (num[2]==num_face) num[2]=elem_faces(elem,3);
884
885 n[0] = face_normale(num_face,0)/sqrt(surf2);
886 n[1] = face_normale(num_face,1)/sqrt(surf2);
887 n[2] = face_normale(num_face,2)/sqrt(surf2);
888
889
890 if( est_egal(n[0],0.) && est_egal(n[1],0.) )
891 {
892 t1[0] = 0.;
893 t1[1] = 1;
894 t1[2] = 0;
895 }
896 else
897 {
898 t1[0] = -n[1]/sqrt(n[0]*n[0]+n[1]*n[1]);
899 t1[1] = n[0]/sqrt(n[0]*n[0]+n[1]*n[1]);
900 t1[2] = 0.;
901 }
902
903 t2[0] = n[1]*t1[2] - n[2]*t1[1];
904 t2[1] = n[2]*t1[0] - n[0]*t1[2];
905 t2[2] = n[0]*t1[1] - n[1]*t1[0];
906
907 dist = distance_3D(num_face,elem,domaine_VEF)*4./3.;
908
909 v_tang[0] = ((vit(num[0],0)+vit(num[1],0)+vit(num[2],0))*t1[0]
910 +(vit(num[0],1)+vit(num[1],1)+vit(num[2],1))*t1[1]
911 +(vit(num[0],2)+vit(num[1],2)+vit(num[2],2))*t1[2])/3.;
912
913 v_tang[1] = ((vit(num[0],0)+vit(num[1],0)+vit(num[2],0))*t2[0]
914 +(vit(num[0],1)+vit(num[1],1)+vit(num[2],1))*t2[1]
915 +(vit(num[0],2)+vit(num[1],2)+vit(num[2],2))*t2[2])/3.;
916
917 F[0] = ((termes_sources(num[0],0)-grad_p(num[0],0))/volumes_entrelaces(num[0])
918 +(termes_sources(num[1],0)-grad_p(num[1],0))/volumes_entrelaces(num[1])
919 +(termes_sources(num[2],0)-grad_p(num[2],0))/volumes_entrelaces(num[2]))/3.;
920
921 F[1] = ((termes_sources(num[0],1)-grad_p(num[0],1))/volumes_entrelaces(num[0])
922 +(termes_sources(num[1],1)-grad_p(num[1],1))/volumes_entrelaces(num[1])
923 +(termes_sources(num[2],1)-grad_p(num[2],1))/volumes_entrelaces(num[2]))/3.;
924
925 F[2] = ((termes_sources(num[0],2)-grad_p(num[0],2))/volumes_entrelaces(num[0])
926 +(termes_sources(num[1],2)-grad_p(num[1],2))/volumes_entrelaces(num[1])
927 +(termes_sources(num[2],2)-grad_p(num[2],2))/volumes_entrelaces(num[2]))/3.;
928
929
930 F_tang[0] = F[0]*t1[0] + F[1]*t1[1] + F[2]*t1[2] ; // projection de F suivant l'axe tangentiel_1
931 F_tang[1] = F[0]*t2[0] + F[1]*t2[1] + F[2]*t2[2] ; // projection de F suivant l'axe tangentiel_2
932 f_tang_moy(compteur_faces_paroi,0) = alpha*F_tang[0] +(1-alpha)*f_tang_moy(compteur_faces_paroi,0);
933 f_tang_moy(compteur_faces_paroi,1) = alpha*F_tang[1] +(1-alpha)*f_tang_moy(compteur_faces_paroi,1);
934
935 // Terme source de Boussinesq si demande
936 if (source_boussinesq==1)
937 {
938 ts_boussi = 0.;
939 const Champ_Don_base& ch_gravite=le_fluide.gravite();
940 const DoubleVect& gravite = ch_gravite.valeurs();
941 const Equation_base& eqn_th = mon_modele_turb_hyd->equation().probleme().equation(1);
942 const Modele_turbulence_scal_base& le_mod_turb_th = ref_cast(Modele_turbulence_scal_base,eqn_th.get_modele(TURBULENCE).valeur());
943 const Turbulence_paroi_scal_base& loi = le_mod_turb_th.loi_paroi();
944 ParoiVEF_TBLE_scal& loi_tble_T = ref_cast_non_const(ParoiVEF_TBLE_scal,loi);
945 if (loi_tble_T.est_initialise())
946 {
947
948 double gt1 = 0.;
949 double gt2 = 0.;
950 for (int idim=0; idim<dimension; idim++)
951 {
952 gt1 += gravite(idim)*t1[idim];
953 gt2 += gravite(idim)*t2[idim];
954 }
955
956 for(int i=0; i<nb_pts+1; i++)
957 {
958 double tmp = beta_t*(loi_tble_T.get_eq_couche_lim_T(compteur_faces_paroi).get_Unp1(0,i)-T0);
959 ts_boussi(0, i) = -gt1*tmp;
960 ts_boussi(1, i) = -gt2*tmp;
961 }
962 /*Cout << "gt1 " << gt1 << finl;
963 Cout << "gt2 " << gt2 << finl;
964 Cout << "T0 " << T0 << finl;
965 Cout << "beta_t " << beta_t << finl;*/
966 }
967 }
968 }
969
970
971 if(tps<=dt_min)
972 {
973 visco_turb_moy(elem)=visco_turb(elem);
974 }
975 else if(tps>tps_start_stat_nu_t_dyn)
976 {
977 visco_turb_moy(elem) = (dt/tps)*visco_turb(elem)+(1-(dt/tps))*visco_turb_moy(elem);
978 }
979
980
981 // au premier appel, dt ne vaut pas dt_min mais 0. car un appel a mettre_a_jour dans le preparer_calcul des modeles de turb
982 // precede l'initialisation du schema en temps.
983 if(dt<dt_min) equation_vitesse.set_dt(1.e12);
984 else equation_vitesse.set_dt(dt);
985
986
987 for (int i=0; i<nb_comp; i++)
988 {
989 equation_vitesse.set_u_y0(i,0.); // vitesse a la paroi
990 equation_vitesse.set_u_yn(i,v_tang[i]); // vitesse en yn
991
992 //equation_vitesse.set_F(i, F_tang[i]);
993 for(int ipts=0 ; ipts<nb_pts+1; ipts++)
994 equation_vitesse.set_F(i, ipts, f_tang_moy(compteur_faces_paroi,i)+ts_boussi(i,ipts));
995 }
996
997
998 /////////////////////////////////////////////////////
999 // On resout les equations aux limites simplifiees //
1000 // pendant un pas de temps du maillage grossier //
1001 /////////////////////////////////////////////////////
1002
1003
1004 if (statio==0)
1005 equation_vitesse.aller_au_temps(tps);
1006 else
1007 equation_vitesse.aller_jusqu_a_convergence(max_it_statio, eps_statio);
1008
1009
1010 itmax_loc = equation_vitesse.get_it();
1011
1012 if(itmax_loc>max_it) Cerr << "WARNING : TOO MANY ITERATIONS ARE NEEDED IN THE TBLE MESH" << finl;
1013
1014 if(itmax_loc>itmax) itmax=itmax_loc;
1015
1016
1017 ///////////////////////////////////////////////
1018 ///// Determination des composantes //////
1019 ///// du cisaillement a la paroi //////
1020 ///////////////////////////////////////////////
1021
1022
1023 double d_visco;
1024
1025 if (l_unif==1)
1026 d_visco = visco;
1027 else
1028 d_visco = tab_visco[elem];
1029
1030 if (dimension == 2)
1031 {
1032 double C0 = equation_vitesse.get_cis(0);
1033
1034 double Cx = t1[0]*C0 ;
1035 double Cy = t1[1]*C0 ;
1036
1037 Cisaillement_paroi_(num_face,0) = Cx;
1038 Cisaillement_paroi_(num_face,1) = Cy;
1039
1040 tab_u_star_(num_face) = pow(Cx*Cx+Cy*Cy,0.25);
1041 tab_d_plus_(num_face) = dist*tab_u_star(num_face)/d_visco;
1042 }
1043 else
1044 {
1045 double C0 = equation_vitesse.get_cis(0);
1046 double C1 = equation_vitesse.get_cis(1);
1047
1048 double Cx = t1[0]*C0 + t2[0]*C1 ;
1049 double Cy = t1[1]*C0 + t2[1]*C1 ;
1050 double Cz = t1[2]*C0 + t2[2]*C1 ;
1051
1052 Cisaillement_paroi_(num_face,0) = Cx;
1053 Cisaillement_paroi_(num_face,1) = Cy;
1054 Cisaillement_paroi_(num_face,2) = Cz;
1055
1056 tab_u_star_(num_face) = pow(Cx*Cx+Cy*Cy+Cz*Cz,0.25);
1057 tab_d_plus_(num_face) = dist*tab_u_star(num_face)/d_visco;
1058 }
1059
1060
1061 //////////////////////////////////////////////////
1062 ///// Mise a jour des grandeurs turbulentes //////
1063 ///// nu_t - k - eps //////
1064 //////////////////////////////////////////////////
1065
1066
1067 if (isKeps==0) // on est en LES
1068 {
1069 if((nu_t_dyn==0)||(tps<tps_nu_t_dyn))
1070 {
1071 tab1(elem)=equation_vitesse.get_nu_t_yn();
1072 }
1073 else
1074 equation_vitesse.set_nu_t_yn(visco_turb_moy(elem));
1075 }
1076 else if (isKeps==1) // on est en K-eps
1077 {
1078 traitement_keps(tab1, num_face, face_voisins, elem_faces, nfac, dist, tab_d_plus(num_face), d_visco, tab_u_star(num_face));
1079 }
1080
1081 compteur_faces_paroi++;
1082
1083 }//Fin boucle sur num_face
1084 }//Fin if Paroi_fixe
1085 }//Fin boucle sur les bords parietaux
1086
1087
1088 if (Process::je_suis_maitre())
1089 {
1090 SFichier fic("iter.dat",ios::app); // ouverture du fichier iter.dat
1091 fic << tps << " " << itmax << finl;
1092 }
1093
1094 if(oui_stats==1)
1095 calculer_stats();
1096
1097 return 1;
1098
1099}
1100
1101void ParoiVEF_TBLE::traitement_keps_BiK(DoubleTab& tab_k,DoubleTab& tab_eps, int num_face, const IntTab& face_voisins,
1102 const IntTab& elem_faces, int nfac, double dist, double d_plus, double d_visco, double u_star)
1103{
1104 IntVect num(nfac);
1105 int elem=face_voisins(num_face,0);
1106 int nf2;
1107
1108 // on determine les face de l'elements qui sont autres que le num_face traite
1109 for (nf2=0; nf2<nfac; nf2++)
1110 {
1111 num[nf2] = elem_faces(elem,nf2);
1112 }
1113 // Maintenant on place le num_face en fin de tableau
1114 for (nf2=0; nf2<nfac-1; nf2++)
1115 {
1116 num[nf2] = elem_faces(elem,nf2);
1117 if (num[nf2] == num_face)
1118 {
1119 num[nf2] = num[nfac-1];
1120 num[nfac-1] = num_face;
1121 }
1122 }
1123
1124 // Boucle sur les faces :
1125 for (int nf=0; nf<nfac; nf++)
1126 {
1127 if (num[nf]==num_face)
1128 {
1129
1130
1131 // if (sub_type(Champ_P1NC,eqn_hydr.inconnue()))
1132 // Strategie pour les tetras :
1133 // On impose k et eps a la paroi :
1134 // approximation: d(k)/d(n) = 0 a la paroi
1135 // c'est faux mais ca marche
1136 tab_k(num[nf]) =0.;
1137 tab_eps(num[nf])=0.;
1138 int nk=0;
1139
1140 for (int k=0; k<nfac; k++)
1141 //if ( (num[k] >= ndebint) && (k != nf))
1142 if ( (face_keps_imposee_[num[k]]>-1) && (k != nf))
1143 {
1144
1145 tab_k(num[nf]) += tab_k(num[k]);
1146 tab_eps(num[nf])+= tab_eps(num[k]);
1147 nk++;
1148 }
1149
1150 if (nk != 0 )
1151 {
1152 tab_k(num[nf]) /=nk;
1153 tab_eps(num[nf])/=nk;
1154 }
1155
1156 }
1157
1158
1159 // On verifie si num[nf] n'est pas une face de bord :
1160
1161 else if ( (face_keps_imposee_[num[nf]]>-1))//if (num[nf]>=ndebint)
1162 {
1163 calculer_k_eps(tab_k(num[nf]),tab_eps(num[nf]),d_plus,u_star,d_visco,dist);
1164 }
1165 }
1166}
1167
1168void ParoiVEF_TBLE::traitement_keps(DoubleTab& tab_k_eps, int num_face, const IntTab& face_voisins,
1169 const IntTab& elem_faces, int nfac, double dist, double d_plus, double d_visco, double u_star)
1170{
1171 IntVect num(nfac);
1172 int elem=face_voisins(num_face,0);
1173 int nf2;
1174
1175 // on determine les face de l'elements qui sont autres que le num_face traite
1176 for (nf2=0; nf2<nfac; nf2++)
1177 {
1178 num[nf2] = elem_faces(elem,nf2);
1179 }
1180 // Maintenant on place le num_face en fin de tableau
1181 for (nf2=0; nf2<nfac-1; nf2++)
1182 {
1183 num[nf2] = elem_faces(elem,nf2);
1184 if (num[nf2] == num_face)
1185 {
1186 num[nf2] = num[nfac-1];
1187 num[nfac-1] = num_face;
1188 }
1189 }
1190
1191 // Boucle sur les faces :
1192 for (int nf=0; nf<nfac; nf++)
1193 {
1194 if (num[nf]==num_face)
1195 {
1196
1197
1198 // if (sub_type(Champ_P1NC,eqn_hydr.inconnue()))
1199 // Strategie pour les tetras :
1200 // On impose k et eps a la paroi :
1201 // approximation: d(k)/d(n) = 0 a la paroi
1202 // c'est faux mais ca marche
1203 tab_k_eps(num[nf],0)=0.;
1204 tab_k_eps(num[nf],1)=0.;
1205 int nk=0;
1206
1207 for (int k=0; k<nfac; k++)
1208 //if ( (num[k] >= ndebint) && (k != nf))
1209 if ( (face_keps_imposee_[num[k]]>-1) && (k != nf))
1210 {
1211
1212 tab_k_eps(num[nf],0)+= tab_k_eps(num[k],0);
1213 tab_k_eps(num[nf],1)+= tab_k_eps(num[k],1);
1214 nk++;
1215 }
1216
1217 if (nk != 0 )
1218 {
1219 tab_k_eps(num[nf],0)/=nk;
1220 tab_k_eps(num[nf],1)/=nk;
1221 }
1222
1223 }
1224
1225
1226 // On verifie si num[nf] n'est pas une face de bord :
1227
1228 else if ( (face_keps_imposee_[num[nf]]>-1))//if (num[nf]>=ndebint)
1229 {
1230 calculer_k_eps(tab_k_eps(num[nf],0),tab_k_eps(num[nf],1),d_plus,u_star,d_visco,dist);
1231
1232
1233 /* double y_plus = d_plus;
1234
1235 if (y_plus<=8)
1236 {
1237 tab_k_eps(elem,0)=0.;
1238 tab_k_eps(elem,1)=0.;
1239 }
1240 else if (y_plus>8 && y_plus<30)
1241 {
1242 double u_star = tab_u_star_(num_face);
1243 double lm_plus = calcul_lm_plus(y_plus);
1244 double deriv = Fdypar_direct(lm_plus);
1245 double x = lm_plus*u_star*deriv;
1246 tab_k_eps(elem,0) = x*x/sqrt(CMU_DEF);
1247 tab_k_eps(elem,1) = (tab_k_eps(elem,0)*u_star*u_star*deriv)*sqrt(CMU_DEF)/d_visco;
1248 }
1249 else if (y_plus>=30)
1250 {
1251 double us2 = tab_u_star_(num_face)*tab_u_star_(num_face);
1252 tab_k_eps(elem,0)=us2/sqrt(CMU_DEF);
1253 tab_k_eps(elem,1)=us2*u_star/K_DEF/dist;
1254 }
1255 */
1256
1257 }
1258 else
1259 {
1260 //Cerr<<"Attention, on n'a rien fait dans Paroi_std_hyd_VEF::calculer_hyd"<<finl;
1261 //Cerr<<"pour le numero de face num[nf] = "<<num[nf]<<" et num_face = "<<num_face<<finl;
1262 //Cerr<<"ndebint = "<<ndebint<<finl;
1263 }
1264 }
1265}
1266
1267int ParoiVEF_TBLE::calculer_k_eps(double& k, double& eps , double yp, double u_star,
1268 double d_visco, double dist)
1269{
1270 // PQ : 05/04/07 : formulation continue de k et epsilon
1271 // assurant le bon comportement asymptotique
1272 k = 0.07*yp*yp*(exp(-yp/9.))+(1./(sqrt(Cmu_)))*pow((1.-exp(-yp/20.)),2); // k_plus
1273 k *= u_star*u_star;
1274 // PL: 50625=15^4 on evite d'utiliser pow car lent
1275 eps = (1./(Kappa_*pow(yp*yp*yp*yp+50625,0.25))); // eps_plus
1276 eps *= pow(u_star,4)/d_visco;
1277
1278 return 1;
1279}
1280
1281int ParoiVEF_TBLE::calculer_stats()
1282{
1283 const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF, le_dom_dis_.valeur());
1284 //const Domaine& domaine = domaine_VEF.domaine();
1285 const DoubleTab& face_normale = domaine_VEF.face_normales();
1286 //const int nfac = domaine.nb_faces_elem();
1287
1288 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
1289 const double tps = eqn_hydr.inconnue().temps();
1290 const double dt = eqn_hydr.schema_temps().pas_de_temps();
1291
1292 int num_face, num_face_global;
1293 double surf;
1294 double nx,ny,nz=0.; // vecteur normal a la paroi
1295 double t1x,t1y,t1z=0.; // vecteur tangentiel_1 a la paroi
1296 double t2x,t2y,t2z=0.; // vecteur tangentiel_2 a la paroi (n^t1)
1297 double Unp1_t1,Unp1_t2=0.;
1298 double Fx, Fy, Fz;
1299 double u,v,w=0.;
1300 //IntVect num(nfac);
1301
1302
1303 //////////////////////////////////////
1304 //Moyennes Temporelles
1305 //////////////////////////////////////
1306 if(((tps>tps_deb_stats) && (tps<tps_fin_stats)) && (oui_stats!=0))
1307 {
1308 if (dimension==2)
1309 {
1310 for(int j=0; j<nb_post_pts; j++)
1311 {
1312 num_face=num_faces_post(j);
1313 num_face_global=num_global_faces_post(j);
1314
1315 surf=sqrt(face_normale(num_face_global,0)*face_normale(num_face_global,0)+
1316 face_normale(num_face_global,1)*face_normale(num_face_global,1));
1317
1318 nx = face_normale(num_face_global,0)/surf;
1319 ny = face_normale(num_face_global,1)/surf;
1320
1321 t1x = -ny;
1322 t1y = nx;
1323
1324 for(int i=0 ; i<nb_pts+1 ; i++)
1325 {
1326 Unp1_t1 = eq_vit[num_face].get_Unp1(0,i);
1327 Fx = t1x*eq_vit[num_face].get_F(0,i);
1328 Fy = t1y*eq_vit[num_face].get_F(0,i);
1329 u = t1x*Unp1_t1;
1330 v = t1y*Unp1_t1;
1331
1332 calculer_stat(j,i,Fx, Fy, 0, u,v,0,dt);
1333 }
1334 }
1335 }
1336 else if (dimension==3)
1337 {
1338 for(int j=0; j<nb_post_pts; j++)
1339 {
1340 num_face=num_faces_post(j);
1341 num_face_global=num_global_faces_post(j);
1342
1343 surf=sqrt(face_normale(num_face_global,0)*face_normale(num_face_global,0)+
1344 face_normale(num_face_global,1)*face_normale(num_face_global,1)+
1345 face_normale(num_face_global,2)*face_normale(num_face_global,2));
1346
1347 nx = face_normale(num_face_global,0)/surf;
1348 ny = face_normale(num_face_global,1)/surf;
1349 nz = face_normale(num_face_global,2)/surf;
1350
1351 if( est_egal(nx,0.) && est_egal(ny,0.) )
1352 {
1353 t1x = 0.;
1354 t1y = 1;
1355 t1z = 0;
1356 }
1357 else
1358 {
1359 t1x = -ny/sqrt(nx*nx+ny*ny);
1360 t1y = nx/sqrt(nx*nx+ny*ny);
1361 t1z = 0.;
1362 }
1363
1364 t2x = ny*t1z - nz*t1y;
1365 t2y = nz*t1x - nx*t1z;
1366 t2z = nx*t1y - ny*t1x;
1367
1368 for(int i=0 ; i<nb_pts+1 ; i++)
1369 {
1370
1371 Unp1_t1 = eq_vit[num_face].get_Unp1(0,i);
1372 Unp1_t2 = eq_vit[num_face].get_Unp1(1,i);
1373 Fx = t1x*eq_vit[num_face].get_F(0,i)+t2x*eq_vit[num_face].get_F(1,i);
1374 Fy = t1y*eq_vit[num_face].get_F(0,i)+t2y*eq_vit[num_face].get_F(1,i);
1375 Fz = t1z*eq_vit[num_face].get_F(0,i)+t2z*eq_vit[num_face].get_F(1,i);
1376
1377 u = t1x*Unp1_t1 + t2x*Unp1_t2 ;
1378 v = t1y*Unp1_t1 + t2y*Unp1_t2 ;
1379 w = t1z*Unp1_t1 + t2z*Unp1_t2 ;
1380
1381 calculer_stat(j,i,Fx, Fy, Fz, u,v,w,dt);
1382 }
1383 }
1384 }
1385 }
1386 return 1;
1387}
1388
1389
1390
1391void ParoiVEF_TBLE::imprimer_ustar(Sortie& os) const
1392{
1393 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
1394 const double tps = eqn_hydr.inconnue().temps();
1395
1396 const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF, le_dom_dis_.valeur());
1397 //const Domaine& domaine = domaine_VEF.domaine();
1398 const DoubleTab& face_normale = domaine_VEF.face_normales();
1399 //const int nfac = domaine.nb_faces_elem();
1400
1401 int num_face, num_face_global;
1402 double surf;
1403 double nx,ny,nz=0.; // vecteur normal a la paroi
1404 double t1x,t1y,t1z=0.; // vecteur tangentiel_1 a la paroi
1405 double t2x,t2y,t2z=0.; // vecteur tangentiel_2 a la paroi (n^t1)
1406 double Unp1_t1,Unp1_t2=0.;
1407 double u,v,w=0.;
1408 double norm_v=0.;
1409 double nut=0.;
1410 //IntVect num(nfac);
1411
1412 for(int j=0; j<nb_post_pts; j++)
1413 {
1414 num_face=num_faces_post(j);
1415 num_face_global=num_global_faces_post(j);
1416
1417 Nom tmp;
1418 tmp="tble_post_";
1419 tmp+=nom_pts[j];
1420 tmp+=".dat";
1421
1422 EcrFicPartage fic_post(tmp, ios::app);
1423 fic_post << "# t="<< tps << " " ;
1424 fic_post << "x= " << domaine_VEF.xv(num_face_global,0) << " " << "y= " << domaine_VEF.xv(num_face_global,1) ;
1425 if (dimension==3) fic_post << " "<< "z= " << domaine_VEF.xv(num_face_global,2);
1426 fic_post << " " << "u*= " << tab_u_star(num_face_global) ;
1427 //fic_post << " " << "dp/dt1= " << eq_vit[num_face].get_F0(0) ;
1428 //if (dimension==3) fic_post << " " << "dp/dt2= " << eq_vit[num_face].get_F0(1) ;
1429
1430 if (dimension==2)
1431 {
1432 surf=sqrt(face_normale(num_face_global,0)*face_normale(num_face_global,0)+
1433 face_normale(num_face_global,1)*face_normale(num_face_global,1));
1434
1435 nx = face_normale(num_face_global,0)/surf;
1436 ny = face_normale(num_face_global,1)/surf;
1437
1438 t1x = -ny;
1439 t1y = nx;
1440
1441 fic_post << " " << "F_x= " << eq_vit[num_face].get_F0(0)*t1x << " F_y= " << eq_vit[num_face].get_F0(0)*t1y << finl;
1442 fic_post << "# d_paroi u v norm_v nu_t" << finl;
1443
1444 for(int i=0; i<nb_pts+1; i++)
1445 {
1446 Unp1_t1 = eq_vit[num_face].get_Unp1(0,i);
1447 nut = (i<nb_pts)? eq_vit[num_face].get_visco_tot(i) : 0;
1448
1449 u = t1x*Unp1_t1;
1450 v = t1y*Unp1_t1;
1451 norm_v=sqrt(u*u+v*v);
1452
1453 fic_post << eq_vit[num_face].get_y(i) << " " << u << " " << v << " " << norm_v << " " << nut << finl;
1454 }
1455 }
1456 else if (dimension==3)
1457 {
1458 surf=sqrt(face_normale(num_face_global,0)*face_normale(num_face_global,0)+
1459 face_normale(num_face_global,1)*face_normale(num_face_global,1)+
1460 face_normale(num_face_global,2)*face_normale(num_face_global,2));
1461
1462 nx = face_normale(num_face_global,0)/surf;
1463 ny = face_normale(num_face_global,1)/surf;
1464 nz = face_normale(num_face_global,2)/surf;
1465
1466 if( est_egal(nx,0.) && est_egal(ny,0.) )
1467 {
1468 t1x = 0.;
1469 t1y = 1.;
1470 t1z = 0.;
1471 }
1472 else
1473 {
1474 t1x = -ny/sqrt(nx*nx+ny*ny);
1475 t1y = nx/sqrt(nx*nx+ny*ny);
1476 t1z = 0.;
1477 }
1478
1479
1480 t2x = ny*t1z - nz*t1y;
1481 t2y = nz*t1x - nx*t1z;
1482 t2z = nx*t1y - ny*t1x;
1483
1484 fic_post << " " << "F_x= " << eq_vit[num_face].get_F0(0)*t1x+eq_vit[num_face].get_F0(1)*t2x
1485 << " F_y= " << eq_vit[num_face].get_F0(0)*t1y+eq_vit[num_face].get_F0(1)*t2y
1486 << " F_z= " << eq_vit[num_face].get_F0(0)*t1z+eq_vit[num_face].get_F0(1)*t2z
1487 << finl;
1488 fic_post << "# dp/dt1= " << eq_vit[num_face].get_F0(0) << " dp/dt2= " << eq_vit[num_face].get_F0(1) << finl;
1489 fic_post << "# d_paroi u v w norm_v nu_t" << finl;
1490
1491 for(int i=0 ; i<nb_pts+1 ; i++)
1492 {
1493 Unp1_t1 = eq_vit[num_face].get_Unp1(0,i);
1494 Unp1_t2 = eq_vit[num_face].get_Unp1(1,i);
1495
1496 u = t1x*Unp1_t1 + t2x*Unp1_t2 ;
1497 v = t1y*Unp1_t1 + t2y*Unp1_t2 ;
1498 w = t1z*Unp1_t1 + t2z*Unp1_t2 ;
1499 norm_v=sqrt(u*u+v*v+w*w);
1500 nut = (i<nb_pts)? eq_vit[num_face].get_visco_tot(i) : 0;
1501
1502 fic_post << eq_vit[num_face].get_y(i) << " " << u << " " << v << " " << w << " " << norm_v << " " << nut << finl;
1503 }
1504 }
1505 fic_post.syncfile();
1506 fic_post.close();
1507 }
1508 Paroi_TBLE_QDM::imprimer_stat(os, tps);
1509
1510
1511}
1512
1513int ParoiVEF_TBLE::sauvegarder(Sortie& os) const
1514{
1515 const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF, le_dom_dis_.valeur());
1516 double tps = mon_modele_turb_hyd->equation().inconnue().temps();
1517 return Paroi_TBLE_QDM::sauvegarder(os, domaine_VEF, le_dom_Cl_dis_.valeur(), tps);
1518}
1519
1520
1521int ParoiVEF_TBLE::reprendre(Entree& is)
1522{
1523 if (le_dom_dis_) // test pour ne pas planter dans "avancer_fichier(...)"
1524 {
1525 const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF, le_dom_dis_.valeur());
1526 double tps_reprise = mon_modele_turb_hyd->equation().schema_temps().temps_courant();
1527 return Paroi_TBLE_QDM::reprendre(is, domaine_VEF, le_dom_Cl_dis_.valeur(), tps_reprise);
1528 }
1529 return 1;
1530}
1531
1532
1533const Probleme_base& ParoiVEF_TBLE::getPbBase() const
1534{
1535 const Probleme_base& pb_base = mon_modele_turb_hyd->equation().probleme();
1536 return pb_base;
1537}
Champ_Don_base
classe Champ_Don_base classe de base des Champs donnes (non calcules)
Definition Champ_Don_base.h:32
Champ_Don_base::valeurs
DoubleTab & valeurs() override
Surcharge Champ_base::valeurs() Renvoie le tableau des valeurs.
Definition Champ_Don_base.h:49
Champ_Inc_base::valeurs
DoubleTab & valeurs() override
Renvoie le tableau des valeurs du champ au temps courant.
Definition Champ_Inc_base.h:77
Champ_Q1NC
Definition Champ_Q1NC.h:26
Champ_Uniforme
classe Champ_Uniforme Represente un champ constant dans l'espace et dans le temps.
Definition Champ_Uniforme.h:26
Champ_base::temps
double temps() const
Renvoie le temps du champ.
Definition Champ_base.cpp:1004
Cond_lim
classe Cond_lim Classe generique servant a representer n'importe quelle classe
Definition Cond_lim.h:31
Diffu_totale_base::setKappa
void setKappa(double k)
Definition Diffu_totale_base.h:53
Diffu_totale_hyd_base
Classe Diffu_totale_hyd_base Classe abstraite calculant la diffusivite totale (somme diffusivite.
Definition Diffu_totale_hyd_base.h:31
Dirichlet_paroi_fixe
classe Dirichlet_paroi_fixe Represente une paroi immobile dans une equation de type Navier_Stokes.
Definition Dirichlet_paroi_fixe.h:26
Domaine_VEF
class Domaine_VEF
Definition Domaine_VEF.h:54
Domaine_VF::volumes_entrelaces
DoubleVect & volumes_entrelaces()
Definition Domaine_VF.h:99
Domaine_VF::face_normales
virtual double face_normales(int face, int comp) const
Definition Domaine_VF.h:47
Domaine_VF::xv
double xv(int num_face, int k) const
Definition Domaine_VF.h:76
Domaine_VF::elem_faces
int elem_faces(int i, int j) const
renvoie le numero de le ieme face de la maille num_elem la facon dont ces faces sont numerotees est
Definition Domaine_VF.h:543
Domaine_VF::face_voisins
int face_voisins(int num_face, int i) const
renvoie l'element voisin de numface dans la direction i.
Definition Domaine_VF.h:418
Domaine_dis_base::nb_front_Cl
int nb_front_Cl() const
Definition Domaine_dis_base.h:53
Domaine_dis_base::domaine
const Domaine & domaine() const
Definition Domaine_dis_base.h:46
EcrFicPartage::syncfile
Sortie & syncfile() override
Provoque l'ecriture sur disque des donnees accumulees sur les differents processeurs depuis le dernie...
Definition EcrFicPartage.cpp:173
Entree
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
Eq_couch_lim
Definition Eq_couch_lim.h:39
Eq_couch_lim::initialiser
void initialiser(int, int, double, double, int, int)
Definition Eq_couch_lim.cpp:43
Eq_couch_lim::set_u_y0
void set_u_y0(int j, double u)
Definition Eq_couch_lim.h:82
Eq_couch_lim::get_nu_t_yn
double get_nu_t_yn()
Definition Eq_couch_lim.h:160
Eq_couch_lim::set_u_yn
void set_u_yn(int j, double u)
Definition Eq_couch_lim.h:87
Eq_couch_lim::get_it
int get_it()
Definition Eq_couch_lim.h:140
Eq_couch_lim::get_cis
double get_cis(int j)
Definition Eq_couch_lim.h:197
Eq_couch_lim::aller_au_temps
void aller_au_temps(double)
Definition Eq_couch_lim.cpp:132
Eq_couch_lim::set_Uinit_lin
void set_Uinit_lin(int comp, double f1, double f2)
Definition Eq_couch_lim.h:95
Eq_couch_lim::get_diffu
Diffu_totale_base & get_diffu()
Definition Eq_couch_lim.h:247
Eq_couch_lim::set_F
void set_F(int j, double f)
Definition Eq_couch_lim.h:55
Eq_couch_lim::set_nu_t_yn
void set_nu_t_yn(double nu_t)
Definition Eq_couch_lim.h:91
Eq_couch_lim::set_yn
void set_yn(double y)
Definition Eq_couch_lim.h:78
Eq_couch_lim::aller_jusqu_a_convergence
void aller_jusqu_a_convergence(int, double)
Definition Eq_couch_lim.cpp:268
Eq_couch_lim::set_y0
void set_y0(double y)
Definition Eq_couch_lim.h:74
Eq_couch_lim::set_dt
void set_dt(double delta_t)
Definition Eq_couch_lim.h:120
Eq_couch_lim::set_Uinit
void set_Uinit(int comp, int i, double val)
Definition Eq_couch_lim.h:108
Equation_base
classe Equation_base Le role d'une equation est le calcul d'un ou plusieurs champs....
Definition Equation_base.h:76
Equation_base::milieu
virtual const Milieu_base & milieu() const =0
Equation_base::get_modele
virtual const RefObjU & get_modele(Type_modele type) const
Definition Equation_base.cpp:1894
Equation_base::sources
Sources & sources()
Renvoie les termes sources asssocies a l'equation.
Definition Equation_base.cpp:348
Equation_base::inconnue
virtual const Champ_Inc_base & inconnue() const =0
Equation_base::probleme
Probleme_base & probleme()
Renvoie le probleme associe a l'equation.
Definition Equation_base.cpp:747
Equation_base::schema_temps
Schema_Temps_base & schema_temps()
Renvoie le schema en temps associe a l'equation.
Definition Equation_base.cpp:865
Fluide_base
classe Fluide_base Cette classe represente un d'un fluide incompressible ainsi que
Definition Fluide_base.h:38
Fluide_base::viscosite_cinematique
const Champ_Don_base & viscosite_cinematique() const
Definition Fluide_base.h:58
Front_VF
class Front_VF
Definition Front_VF.h:36
Front_VF::nb_faces
int nb_faces() const
Definition Front_VF.h:53
Front_VF::num_face
int num_face(const int) const
Definition Front_VF.h:68
Milieu_base::beta_t
virtual const Champ_Don_base & beta_t() const
Renvoie beta_t du milieu.
Definition Milieu_base.cpp:887
Milieu_base::gravite
virtual const Champ_Don_base & gravite() const
Renvoie la gravite du milieu si elle a ete associe provoque une erreur sinon.
Definition Milieu_base.cpp:716
Modele_turbulence_scal_base
Classe Modele_turbulence_scal_base Cette classe represente un modele de turbulence pour une equation ...
Definition Modele_turbulence_scal_base.h:40
Modele_turbulence_scal_base::loi_paroi
const Turbulence_paroi_scal_base & loi_paroi() const
Renvoie la loi de turbulence sur la paroi (version const).
Definition Modele_turbulence_scal_base.h:90
Motcle
Une chaine de caractere (Nom) en majuscules.
Definition Motcle.h:26
Navier_Stokes_std
classe Navier_Stokes_std Cette classe porte les termes de l'equation de la dynamique
Definition Navier_Stokes_std.h:56
Navier_Stokes_std::operateur_gradient
Operateur_Grad & operateur_gradient()
Renvoie l'operateur de calcul du gradient associe a l'equation.
Definition Navier_Stokes_std.cpp:568
Navier_Stokes_std::pression
Champ_Inc_base & pression()
Definition Navier_Stokes_std.h:123
Nom
class Nom Une chaine de caractere pour nommer les objets de TRUST
Definition Nom.h:31
Objet_U::set_param
virtual void set_param(Param &) const
Definition Objet_U.h:135
Objet_U::Entree
friend class Entree
Definition Objet_U.h:76
Objet_U::dimension
static int dimension
Definition Objet_U.h:99
Objet_U::Sortie
friend class Sortie
Definition Objet_U.h:75
Objet_U::que_suis_je
const Nom & que_suis_je() const
renvoie la chaine identifiant la classe.
Definition Objet_U.cpp:104
Objet_U::readOn
virtual Entree & readOn(Entree &)
Lecture d'un Objet_U sur un flot d'entree Methode a surcharger.
Definition Objet_U.cpp:293
Objet_U::le_nom
virtual const Nom & le_nom() const
Donne le nom de l'Objet_U Methode a surcharger : renvoie "neant" dans cette implementation.
Definition Objet_U.cpp:319
Objet_U::printOn
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
Definition Objet_U.cpp:282
Operateur_Grad
Classe Operateur_Grad Classe generique de la hierarchie des operateurs calculant le gradient.
Definition Operateur_Grad.h:31
Param
Helper class to factorize the readOn method of Objet_U classes.
Definition Param.h:112
Param::ajouter_condition
void ajouter_condition(const char *condition, const char *message, const char *name=0)
Declare a post-read logical condition that must hold on the parameter values.
Definition Param.cpp:496
Param::ajouter
void ajouter(const char *keyword, const int *value, Param::Nature nat=Param::OPTIONAL)
Register an integer parameter.
Definition Param.cpp:364
Param::ajouter_non_std
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
ParoiVEF_TBLE_scal
CLASS: ParoiVEF_TBLE_scal.
Definition ParoiVEF_TBLE_scal.h:39
ParoiVEF_TBLE
CLASS: ParoiVEF_TBLE.
Definition ParoiVEF_TBLE.h:39
ParoiVEF_TBLE::set_param
void set_param(Param &param) const override
Definition ParoiVEF_TBLE.cpp:62
ParoiVEF_TBLE::ParoiVEF_TBLE
ParoiVEF_TBLE()
Definition ParoiVEF_TBLE.cpp:40
ParoiVEF_TBLE::init_lois_paroi
int init_lois_paroi() override
Definition ParoiVEF_TBLE.cpp:96
ParoiVEF_TBLE::imprimer_ustar
void imprimer_ustar(Sortie &) const override
Definition ParoiVEF_TBLE.cpp:1391
ParoiVEF_TBLE::sauvegarder
int sauvegarder(Sortie &) const override
Sauvegarde d'un Objet_U sur un flot de sortie Methode a surcharger.
Definition ParoiVEF_TBLE.cpp:1513
ParoiVEF_TBLE::reprendre
int reprendre(Entree &) override
Reprise d'un Objet_U sur un flot d'entree Methode a surcharger.
Definition ParoiVEF_TBLE.cpp:1521
ParoiVEF_TBLE::calculer_hyd
int calculer_hyd(DoubleTab &) override
Definition ParoiVEF_TBLE.cpp:302
ParoiVEF_TBLE::calculer_hyd_BiK
int calculer_hyd_BiK(DoubleTab &, DoubleTab &) override
Definition ParoiVEF_TBLE.cpp:311
ParoiVEF_TBLE::getPbBase
const Probleme_base & getPbBase() const override
Definition ParoiVEF_TBLE.cpp:1533
ParoiVEF_TBLE::lire_motcle_non_standard
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.
Definition ParoiVEF_TBLE.cpp:81
Paroi_TBLE_QDM_Scal::est_initialise
int est_initialise() const
Definition Paroi_TBLE_QDM_Scal.h:93
Paroi_TBLE_QDM_Scal::get_eq_couche_lim_T
Eq_couch_lim & get_eq_couche_lim_T(int i)
Definition Paroi_TBLE_QDM_Scal.h:84
Paroi_TBLE_QDM::tps_start_stat_nu_t_dyn
double tps_start_stat_nu_t_dyn
Definition Paroi_TBLE_QDM.h:101
Paroi_TBLE_QDM::visco_turb_moy
DoubleTab visco_turb_moy
Definition Paroi_TBLE_QDM.h:108
Paroi_TBLE_QDM::calculer_stat
void calculer_stat(int indice_post, int indice_maillage, double Fx, double Fy, double Fz, double u, double v, double w, double dt)
Definition Paroi_TBLE_QDM.cpp:505
Paroi_TBLE_QDM::nu_t_dyn
int nu_t_dyn
Definition Paroi_TBLE_QDM.h:109
Paroi_TBLE_QDM::oui_stats
int oui_stats
Definition Paroi_TBLE_QDM.h:91
Paroi_TBLE_QDM::num_faces_post
IntTab num_faces_post
Definition Paroi_TBLE_QDM.h:114
Paroi_TBLE_QDM::imprimer_stat
void imprimer_stat(Sortie &, double) const
Definition Paroi_TBLE_QDM.cpp:445
Paroi_TBLE_QDM::tps_deb_stats
double tps_deb_stats
Definition Paroi_TBLE_QDM.h:99
Paroi_TBLE_QDM::nom_pts
Noms nom_pts
Definition Paroi_TBLE_QDM.h:113
Paroi_TBLE_QDM::fac
double fac
Definition Paroi_TBLE_QDM.h:98
Paroi_TBLE_QDM::set_param
void set_param(Param &param) const
Definition Paroi_TBLE_QDM.cpp:62
Paroi_TBLE_QDM::nb_comp
int nb_comp
Definition Paroi_TBLE_QDM.h:91
Paroi_TBLE_QDM::beta_t
double beta_t
Definition Paroi_TBLE_QDM.h:147
Paroi_TBLE_QDM::max_it_statio
int max_it_statio
Definition Paroi_TBLE_QDM.h:94
Paroi_TBLE_QDM::eps_statio
double eps_statio
Definition Paroi_TBLE_QDM.h:96
Paroi_TBLE_QDM::reprendre
int reprendre(Entree &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, double tps)
Definition Paroi_TBLE_QDM.cpp:340
Paroi_TBLE_QDM::T0
double T0
Definition Paroi_TBLE_QDM.h:146
Paroi_TBLE_QDM::num_global_faces_post
IntTab num_global_faces_post
Definition Paroi_TBLE_QDM.h:114
Paroi_TBLE_QDM::source_boussinesq
int source_boussinesq
Definition Paroi_TBLE_QDM.h:145
Paroi_TBLE_QDM::epsilon
double epsilon
Definition Paroi_TBLE_QDM.h:98
Paroi_TBLE_QDM::sauvegarder
int sauvegarder(Sortie &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, double tps) const
Definition Paroi_TBLE_QDM.cpp:408
Paroi_TBLE_QDM::valeurs_reprises
DoubleTab valeurs_reprises
Definition Paroi_TBLE_QDM.h:117
Paroi_TBLE_QDM::reprise_ok
int reprise_ok
Definition Paroi_TBLE_QDM.h:116
Paroi_TBLE_QDM::max_it
int max_it
Definition Paroi_TBLE_QDM.h:91
Paroi_TBLE_QDM::lire_motcle_non_standard
int lire_motcle_non_standard(const Motcle &, Entree &)
Definition Paroi_TBLE_QDM.cpp:90
Paroi_TBLE_QDM::nb_post_pts
int nb_post_pts
Definition Paroi_TBLE_QDM.h:111
Paroi_TBLE_QDM::tps_nu_t_dyn
double tps_nu_t_dyn
Definition Paroi_TBLE_QDM.h:102
Paroi_TBLE_QDM::statio
int statio
Definition Paroi_TBLE_QDM.h:93
Paroi_TBLE_QDM::init_lois_paroi
int init_lois_paroi(const Domaine_VF &, const Domaine_Cl_dis_base &)
Definition Paroi_TBLE_QDM.cpp:192
Paroi_TBLE_QDM::tps_fin_stats
double tps_fin_stats
Definition Paroi_TBLE_QDM.h:99
Paroi_TBLE_QDM::nb_pts
int nb_pts
Definition Paroi_TBLE_QDM.h:91
Paroi_hyd_base_VEF
CLASS: Paroi_hyd_base_VEF Classe de base des lois de paroi hydraulique en VEF.
Definition Paroi_hyd_base_VEF.h:28
Paroi_hyd_base_VEF::face_keps_imposee_
IntVect face_keps_imposee_
Definition Paroi_hyd_base_VEF.h:38
Paroi_hyd_base_VEF::flag_face_keps_imposee_
int flag_face_keps_imposee_
Definition Paroi_hyd_base_VEF.h:39
Paroi_hyd_base_VEF::init_lois_paroi_
void init_lois_paroi_()
Definition Paroi_hyd_base_VEF.cpp:37
Paroi_log_QDM::Cmu_
double Cmu_
Definition Paroi_log_QDM.h:52
Paroi_log_QDM::set_param
void set_param(Param &param) const
Definition Paroi_log_QDM.cpp:22
Paroi_log_QDM::Kappa_
double Kappa_
Definition Paroi_log_QDM.h:53
Probleme_base
classe Probleme_base C'est un Probleme_U qui n'est pas un couplage.
Definition Probleme_base.h:55
Probleme_base::equation
virtual const Equation_base & equation(int) const =0
Process::exit
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
Process::je_suis_maitre
static int je_suis_maitre()
renvoie 1 si on est sur le processeur maitre du groupe courant (c'est a dire me() == 0),...
Definition Process.cpp:86
SFichier
Cette classe est a la classe C++ ofstream ce que la classe Sortie est a la classe C++ ostream Elle re...
Definition SFichier.h:27
Schema_Temps_base::temps_courant
double temps_courant() const
Renvoie le temps courant.
Definition Schema_Temps_base.h:445
Schema_Temps_base::pas_temps_min
double pas_temps_min() const
Renvoie le pas de temps minimum.
Definition Schema_Temps_base.h:407
Schema_Temps_base::pas_de_temps
double pas_de_temps() const
Renvoie le pas de temps (delta_t) courant.
Definition Schema_Temps_base.h:393
Sortie
Classe de base des flux de sortie.
Definition Sortie.h:52
Source_base
classe Source_base Un objet Source_base est un terme apparaissant au second membre d'une
Definition Source_base.h:42
Sources
class Sources Sources represente une liste de Source.
Definition Sources.h:31
TRUSTVect::local_min_vect
_TYPE_ local_min_vect(Mp_vect_options opt=VECT_REAL_ITEMS) const
Definition TRUSTVect.h:155
TRUST_Ref_Objet_U::valeur
const Objet_U & valeur() const
Definition TRUST_Ref.h:134
Terme_Boussinesq_base
Classe Terme_Boussinesq_base Cette classe represente le terme de gravite qui figure dans l'equation.
Definition Terme_Boussinesq_base.h:55
Terme_Boussinesq_base::Scalaire0
double Scalaire0(int i) const
Definition Terme_Boussinesq_base.h:63
Turbulence_paroi_base::tab_d_plus_
DoubleVect tab_d_plus_
Definition Turbulence_paroi_base.h:88
Turbulence_paroi_base::tab_d_plus
const DoubleVect & tab_d_plus() const
Definition Turbulence_paroi_base.h:141
Turbulence_paroi_base::tab_u_star
const DoubleVect & tab_u_star() const
Definition Turbulence_paroi_base.h:131
Turbulence_paroi_base::tab_u_star_
DoubleVect tab_u_star_
Definition Turbulence_paroi_base.h:87
Turbulence_paroi_base::Cisaillement_paroi_
DoubleTab Cisaillement_paroi_
Definition Turbulence_paroi_base.h:85
Turbulence_paroi_scal_base
Classe Turbulence_paroi_scal_base Classe de base pour la hierarchie des classes representant les mode...
Definition Turbulence_paroi_scal_base.h:42