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