TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Modele_Lam_Bremhorst_VEF.cpp
1/****************************************************************************
2* Copyright (c) 2015 - 2016, 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#include <Modele_Lam_Bremhorst_VEF.h>
16#include <Domaine_VEF.h>
17#include <Domaine_Cl_VEF.h>
18#include <Periodique.h>
19#include <Champ_Uniforme.h>
20#include <Scatter.h>
21#include <Champ_P1NC.h>
22#include <Champ_P0_VEF.h>
23#include <Discretisation_base.h>
24#include <Check_espace_virtuel.h>
25#include <LecFicDiffuse.h>
26#include <EcritureLectureSpecial.h>
27#include <Debog.h>
28
29Implemente_instanciable(Modele_Lam_Bremhorst_VEF,"Modele_Lam_Bremhorst_VEF",Modele_Fonc_Bas_Reynolds_Base);
30// XD Lam_Bremhorst modele_fonction_bas_reynolds_base Lam_Bremhorst INHERITS_BRACE Model described in ' C.K.G.Lam and
31// XD_CONT K.Bremhorst, A modified form of the k- epsilon model for predicting wall turbulence, ASME J. Fluids Engng.,
32// XD_CONT Vol.103, p456, (1981)'. Only in VEF.
33// printOn et readOn
34
35Sortie& Modele_Lam_Bremhorst_VEF::printOn(Sortie& s ) const
36{
37 return s;
38}
39
40Entree& Modele_Lam_Bremhorst_VEF::readOn(Entree& is )
41{
42 Param param(que_suis_je());
43 set_param(param);
44 param.lire_avec_accolades_depuis(is);
45 if (is_Reynolds_stress_isotrope_ == 0)
46 is_Cmu_constant_ = 0;
47 return is;
48}
49
50void Modele_Lam_Bremhorst_VEF::set_param(Param& param) const
51{
52 param.ajouter("fichier_distance_paroi",&nom_fic); // XD_ADD_P chaine
53 // XD_CONT refer to distance_paroi keyword
54 param.ajouter("Reynolds_stress_isotrope",&is_Reynolds_stress_isotrope_); // XD_ADD_P int
55 // XD_CONT keyword for isotropic Reynolds stress
56}
57
58int Modele_Lam_Bremhorst_VEF::Calcul_is_Reynolds_stress_isotrope() const
59{
60 return is_Reynolds_stress_isotrope_;
61}
62
63int Modele_Lam_Bremhorst_VEF::Calcul_is_Cmu_constant() const
64{
65 return is_Cmu_constant_;
66}
67
68Entree& Modele_Lam_Bremhorst_VEF::lire(const Motcle& , Entree& is)
69{
70 return is;
71}
72///////////////////////////////////////////////////////////////
73// Implementation des fonctions de la classe
74///////////////////////////////////////////////////////////////
75
76void Modele_Lam_Bremhorst_VEF::associer(const Domaine_dis_base& domaine_dis,
77 const Domaine_Cl_dis_base& domaine_Cl_dis)
78{
79 le_dom_VEF = ref_cast(Domaine_VEF,domaine_dis);
80 le_dom_Cl_VEF = ref_cast(Domaine_Cl_VEF,domaine_Cl_dis);
81}
82
83DoubleTab& Modele_Lam_Bremhorst_VEF::Calcul_D(DoubleTab& D,const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis,
84 const DoubleTab& vitesse,const DoubleTab& K_eps_Bas_Re, const Champ_Don_base& ch_visco ) const
85{
86 D = 0;
87 return D;
88}
89
90DoubleTab& Modele_Lam_Bremhorst_VEF::Calcul_E(DoubleTab& E,const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis, const DoubleTab& transporte,const DoubleTab& K_eps_Bas_Re,const Champ_Don_base& ch_visco, const DoubleTab& visco_turb ) const
91{
92 E = 0;
93 return E;
94}
95
96/* DoubleTab& Modele_Lam_Bremhorst_VEF::Calcul_F1( DoubleTab& F1, const Domaine_dis_base& domaine_dis) const
97{
98 const Domaine_VEF& le_dom = ref_cast(Domaine_VEF,domaine_dis);
99 int nb_faces = le_dom.nb_faces();
100 for (int num_face=0; num_face <nb_faces; num_face ++ )
101 F1[num_face] = 1.;
102 return F1;
103}
104*/
105DoubleTab& Modele_Lam_Bremhorst_VEF::Calcul_F1( DoubleTab& F1, const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis, const DoubleTab& P, const DoubleTab& K_eps_Bas_Re,const Champ_base& ch_visco) const
106{
107
108 double visco=-1;
109 const DoubleTab& tab_visco=ch_visco.valeurs();
110 int is_visco_const=sub_type(Champ_Uniforme,ch_visco);
111 if (is_visco_const)
112 visco=tab_visco(0,0);
113 const Domaine_VEF& le_dom = ref_cast(Domaine_VEF,domaine_dis);
114 const DoubleTab& wall_length = BR_wall_length_->valeurs();
115 DoubleTab wall_length_face(0);
116 le_dom.creer_tableau_faces(wall_length_face);
117 DoubleTab Pderive(0);
118 le_dom.creer_tableau_faces(Pderive);
119 DoubleTab Fmu_loc(0);
120 le_dom.creer_tableau_faces(Fmu_loc);
121 int nb_faces = le_dom.nb_faces();
122 const IntTab& face_voisins = le_dom.face_voisins();
123 int num_face;
124 double Rey,Re;
125
126 // Calcul de la distance a la paroi aux faces
127 for (num_face=0; num_face<nb_faces ; num_face++)
128 {
129 int elem0 = face_voisins(num_face,0);
130 int elem1 = face_voisins(num_face,1);
131 if (elem1!=-1)
132 wall_length_face(num_face) = 0.5*wall_length(elem0)+0.5*wall_length(elem1);
133 else
134 wall_length_face(num_face) = wall_length(elem0);
135 }
136
137 for (num_face=0; num_face <nb_faces; num_face ++ )
138 {
139 if (visco>BR_EPS && K_eps_Bas_Re(num_face,1)>BR_EPS && K_eps_Bas_Re(num_face,0))
140 {
141 Re = (K_eps_Bas_Re(num_face,0)*K_eps_Bas_Re(num_face,0))/(visco*K_eps_Bas_Re(num_face,1));
142 Rey = wall_length_face(num_face)*sqrt(K_eps_Bas_Re(num_face,0))/visco;
143 Fmu_loc(num_face) = (1. - exp(-0.0165*Rey-BR_EPS))*(1. - exp(-0.0165*Rey-BR_EPS))*(1.+20.5/(Re+BR_EPS));
144 }
145 else
146 Fmu_loc(num_face) = 1.;
147
148 if (Fmu_loc(num_face) > 1. + 1.e-8)
149 {
150 Fmu_loc(num_face) = 1.;
151 //Cerr << " On force Fmu_loc a 1 " << finl;
152 //exit();
153 }
154 }
155// Cerr<< Fmu_loc.mp_min_vect()<< " Fmu_loc "<<Fmu_loc.mp_max_vect()<<finl;
156
157 // Calcul de F1
158 for (num_face=0; num_face <nb_faces; num_face ++ )
159 {
160 if (Fmu_loc(num_face) > BR_EPS && (visco>BR_EPS && K_eps_Bas_Re(num_face,1)>BR_EPS && K_eps_Bas_Re(num_face,0)))
161 F1[num_face] = 1.+ (0.05/(Fmu_loc(num_face)+BR_EPS))*(0.05/(Fmu_loc(num_face)+BR_EPS))*(0.05/(Fmu_loc(num_face)+BR_EPS));
162 else
163 F1[num_face] = 1.;
164 }
165// Cerr<< F1.mp_min_vect()<< " F1 "<<F1.mp_max_vect()<<finl;
166 return F1;
167}
168
169DoubleTab& Modele_Lam_Bremhorst_VEF::Calcul_F2( DoubleTab& F2, DoubleTab& Deb, const Domaine_dis_base& domaine_dis,const DoubleTab& K_eps_Bas_Re,const Champ_base& ch_visco ) const
170{
171 double visco=-1;
172 const DoubleTab& tab_visco=ch_visco.valeurs();
173 int is_visco_const=sub_type(Champ_Uniforme,ch_visco);
174 if (is_visco_const)
175 visco=tab_visco(0,0);
176 const Domaine_VEF& le_dom = ref_cast(Domaine_VEF,domaine_dis);
177 int nb_faces = le_dom.nb_faces();
178 int num_face;
179 double Re;
180
181 for (num_face=0; num_face<nb_faces ; num_face++)
182 {
183 if (!is_visco_const)
184 {
185 int elem0 = le_dom.face_voisins(num_face,0);
186 int elem1 = le_dom.face_voisins(num_face,1);
187 if (elem1!=-1)
188 {
189 visco = tab_visco(elem0)*le_dom.volumes(elem0)+tab_visco(elem1)*le_dom.volumes(elem1);
190 visco /= le_dom.volumes(elem0) + le_dom.volumes(elem1);
191 }
192 else
193 visco = tab_visco(elem0);
194
195 }
196 if (visco>BR_EPS && K_eps_Bas_Re(num_face,1)>BR_EPS)
197 {
198 Re = (K_eps_Bas_Re(num_face,0)*K_eps_Bas_Re(num_face,0))/(visco*K_eps_Bas_Re(num_face,1));
199 F2[num_face] = 1. - (exp(-Re*Re));
200 }
201 else
202 F2[num_face] = 1.;
203
204 }
205// Cerr<<F2.mp_min_vect()<<" F2 "<<F2.mp_max_vect()<<finl;
206 return F2;
207}
208
209DoubleTab& Modele_Lam_Bremhorst_VEF::Calcul_Fmu( DoubleTab& Fmu,const Domaine_dis_base& domaine_dis,const Domaine_Cl_dis_base& domaine_Cl_dis,const DoubleTab& K_eps_Bas_Re,const Champ_Don_base& ch_visco ) const
210{
211 double visco=-1;
212 const DoubleTab& tab_visco=ch_visco.valeurs();
213 int is_visco_const=sub_type(Champ_Uniforme,ch_visco);
214 if (is_visco_const)
215 visco=tab_visco(0,0);
216 const Domaine_VEF& le_dom = ref_cast(Domaine_VEF,domaine_dis);
217 const DoubleTab& wall_length = BR_wall_length_->valeurs();
218 DoubleTab wall_length_face(0);
219 le_dom.creer_tableau_faces(wall_length_face);
220 DoubleTab Pderive(0);
221 le_dom.creer_tableau_faces(Pderive);
222 int nb_faces = le_dom.nb_faces();
223 const IntTab& face_voisins = le_dom.face_voisins();
224 int num_face;
225 double Rey,Re;
226
227 // Calcul de la distance a la paroi aux faces
228 for (num_face=0; num_face<nb_faces ; num_face++)
229 {
230 int elem0 = face_voisins(num_face,0);
231 int elem1 = face_voisins(num_face,1);
232 if (elem1!=-1)
233 wall_length_face(num_face) = 0.5*wall_length(elem0)+0.5*wall_length(elem1);
234 else
235 wall_length_face(num_face) = wall_length(elem0);
236 }
237
238 for (num_face=0; num_face <nb_faces; num_face ++ )
239 {
240 if (visco>BR_EPS && K_eps_Bas_Re(num_face,1)>BR_EPS && K_eps_Bas_Re(num_face,0))
241 {
242 Re = (K_eps_Bas_Re(num_face,0)*K_eps_Bas_Re(num_face,0))/(visco*K_eps_Bas_Re(num_face,1));
243 Rey = wall_length_face(num_face)*sqrt(K_eps_Bas_Re(num_face,0))/visco;
244 Fmu(num_face) = (1. - exp(-0.0165*Rey-BR_EPS))*(1. - exp(-0.0165*Rey-BR_EPS))*(1.+20.5/(Re+BR_EPS));
245
246 }
247 else
248 {
249 Fmu(num_face) = 1.;
250 }
251 if (Fmu(num_face) > 1. + 1.e-8)
252 {
253 Fmu(num_face) = 1.;
254 //Cerr << " On force Fmu a 1 " << finl;
255 //exit();
256 }
257 }
258// Cerr<< Fmu.mp_min_vect()<< " Fmu "<<Fmu.mp_max_vect()<<finl;
259
260 return Fmu;
261}
262
263void Modele_Lam_Bremhorst_VEF::mettre_a_jour(double temps)
264{
265 ;
266}
267
268// Initialisation d'une matrice aux elements
269void Modele_Lam_Bremhorst_VEF::init_tenseur_elem(DoubleTab& Tenseur, const Domaine_VEF& domaine_VEF, const int ndim) const
270{
271 if(!Tenseur.get_md_vector())
272 {
273 if (ndim==1)
274 Tenseur.resize(0, Objet_U::dimension);
275 else if (ndim==2)
276 Tenseur.resize(0, Objet_U::dimension, Objet_U::dimension);
277 domaine_VEF.domaine().creer_tableau_elements(Tenseur);
278 }
279 Tenseur = 0.;
280}
281
282
283// Initialisation d'une matrice aux faces
284void Modele_Lam_Bremhorst_VEF::init_tenseur_face(DoubleTab& Tenseur,const Domaine_VEF& domaine_VEF, const int ndim) const
285{
286 if(!Tenseur.get_md_vector())
287 {
288 if (ndim==1)
289 Tenseur.resize(0, Objet_U::dimension);
290 else if (ndim==2)
291 Tenseur.resize(0, Objet_U::dimension, Objet_U::dimension);
292 domaine_VEF.creer_tableau_faces(Tenseur);
293 }
294 Tenseur = 0.;
295}
296
297// Calcul d'un tenseur aux faces a partir d'un tenseur aux elements
298DoubleTab& Modele_Lam_Bremhorst_VEF::calcul_tenseur_face(DoubleTab& Tenseur_face, const DoubleTab& Tenseur_elem,
299 const Domaine_VEF& domaine_VEF, const Domaine_Cl_VEF& domaine_Cl_VEF) const
300{
301 assert_espace_virtuel_vect(Tenseur_elem);
302 const IntTab& face_voisins = domaine_VEF.face_voisins();
303 int nb_faces = domaine_VEF.nb_faces();
304
305 const Conds_lim& les_cl = domaine_Cl_VEF.les_conditions_limites();
306 int nb_cl=les_cl.size();
307 const DoubleVect& volumes = domaine_VEF.volumes();
308
309 for (int n_bord=0; n_bord<nb_cl; n_bord++)
310 {
311 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
312 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
313 int ndeb = le_bord.num_premiere_face();
314 int nfin = ndeb + le_bord.nb_faces();
315
316 if (sub_type(Periodique,la_cl.valeur()))
317 {
318 for (int fac=ndeb; fac<nfin; fac++)
319 {
320 int poly1 = face_voisins(fac,0);
321 int poly2 = face_voisins(fac,1);
322 double a=volumes(poly1)/(volumes(poly1)+volumes(poly2));
323 double b=volumes(poly2)/(volumes(poly1)+volumes(poly2));
324 for (int i=0; i<dimension; i++)
325 for (int j=0; j<dimension; j++)
326 Tenseur_face(fac,i,j) = a*Tenseur_elem(poly1,i,j) + b*Tenseur_elem(poly2,i,j);
327 }
328 }
329 else
330 {
331 for (int fac=ndeb; fac<nfin; fac++)
332 {
333 int poly1 = face_voisins(fac,0);
334 for (int i=0; i<dimension; i++)
335 for (int j=0; j<dimension; j++)
336 Tenseur_face(fac,i,j) = Tenseur_elem(poly1,i,j);
337 }
338 }
339 }
340 int n0 = domaine_VEF.premiere_face_int();
341 for (int fac = n0; fac<nb_faces; fac++)
342 {
343 int poly1 = face_voisins(fac,0);
344 int poly2 = face_voisins(fac,1);
345 double a=volumes(poly1)/(volumes(poly1)+volumes(poly2));
346 double b=volumes(poly2)/(volumes(poly1)+volumes(poly2));
347 for (int i=0; i<dimension; i++)
348 for (int j=0; j<dimension; j++)
349 Tenseur_face(fac,i,j) = a*Tenseur_elem(poly1,i,j) + b*Tenseur_elem(poly2,i,j);
350 }
351
352 return Tenseur_face;
353}
354
355bool Modele_Lam_Bremhorst_VEF::calcul_tenseur_Re(const DoubleTab& nu_turb, const DoubleTab& G, DoubleTab& Re) const
356{
357 const Domaine_dis_base& domaine_dis = mon_equation->domaine_dis();
358 const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF,domaine_dis);
359 int nelem = G.dimension(0);
360
361 DoubleTab S, R;
362 init_tenseur_elem(S,domaine_VEF,2);
363 init_tenseur_elem(R,domaine_VEF,2);
364
365 for (int elem=0; elem<nelem; elem++)
366 for (int i=0; i<Objet_U::dimension; i++)
367 for (int j=0; j<Objet_U::dimension; j++)
368 {
369 S(elem,i,j) = G(elem,i,j) + G(elem,j,i);
370 R(elem,i,j) = G(elem,i,j) - G(elem,j,i);
371 }
372
373 DoubleTab Sn = calcul_norme_elem(domaine_VEF,S);
374
375 Re = calcul_tenseur_Re_elem(mon_equation->discretisation(),domaine_dis,G,S,R,Sn,mon_equation->inconnue());
376
377 return true;
378}
379
380bool Modele_Lam_Bremhorst_VEF::calcul_tenseur_Re_BiK(const DoubleTab& nu_turb, const DoubleTab& G, DoubleTab& Re) const
381{
382 const Domaine_dis_base& domaine_dis = mon_equation->domaine_dis();
383 const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF,domaine_dis);
384 int nelem = G.dimension(0);
385
386 DoubleTab S, R;
387 init_tenseur_elem(S,domaine_VEF,2);
388 init_tenseur_elem(R,domaine_VEF,2);
389
390 for (int elem=0; elem<nelem; elem++)
391 for (int i=0; i<Objet_U::dimension; i++)
392 for (int j=0; j<Objet_U::dimension; j++)
393 {
394 S(elem,i,j) = G(elem,i,j) + G(elem,j,i);
395 R(elem,i,j) = G(elem,i,j) - G(elem,j,i);
396 }
397
398 DoubleTab Sn = calcul_norme_elem(domaine_VEF,S);
399
400 Re = calcul_tenseur_Re_elem_BiK(mon_equation->discretisation(),domaine_dis,G,S,R,Sn,mon_equation->inconnue(),ma_seconde_equation->inconnue());
401
402 return true;
403}
404
405DoubleTab Modele_Lam_Bremhorst_VEF::calcul_tenseur_Re_elem(const Discretisation_base& dis,
406 const Domaine_dis_base& domaine_dis,
407 const DoubleTab& G, const DoubleTab& S,
408 const DoubleTab& R, const DoubleTab& Sn,
409 const Champ_base& K_Eps) const
410{
411 int nelem = domaine_dis.domaine().nb_elem();
412
413 OWN_PTR(Champ_Fonc_base) K_Eps_elem;
414 Noms noms(2), unites(2);
415 noms[0]="K";
416 noms[1]="eps";
417 unites[0]="m2/s2";
418 unites[1]="m2/s3";
419 const Motcle wtype= "champ_elem";
420 dis.discretiser_champ(wtype,domaine_dis,multi_scalaire,noms,unites,2,0.,K_Eps_elem);
421 K_Eps_elem->affecter(K_Eps);
422 DoubleTab tab_K_Eps = K_Eps_elem->valeurs();
423
424 DoubleTab ReNL;
425 const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF,domaine_dis);
426 init_tenseur_elem(ReNL,domaine_VEF,2);
427
428 for (int elem=0; elem<nelem; elem++)
429 {
430 double kseps;
431 if (tab_K_Eps(elem,1) <= BR_EPS)
432 {
433 kseps = tab_K_Eps(elem,0)/BR_EPS;
434 }
435 else
436 {
437 kseps = tab_K_Eps(elem,0)/tab_K_Eps(elem,1);
438 }
439 double Cmu = (2./3.)/(A1+sqrt(0.5)*Sn(elem)*kseps);
440 double C1 = CNL1/((CNL4+CNL5*pow(sqrt(0.5)*Sn(elem)*kseps,3.))*Cmu)* kseps;
441 double C2 = CNL2/((CNL4+CNL5*pow(sqrt(0.5)*Sn(elem)*kseps,3.))*Cmu)* kseps;
442 double C3 = CNL3/((CNL4+CNL5*pow(sqrt(0.5)*Sn(elem)*kseps,3.))*Cmu)* kseps;
443 for (int i=0; i<Objet_U::dimension; i++)
444 for (int j=0; j<Objet_U::dimension; j++)
445 {
446 ReNL(elem,i,j) = S(elem,i,j);
447 for (int k=0; k<Objet_U::dimension; k++)
448 ReNL(elem,i,j) -= C1* S(elem,i,k)*S(elem,k,j) +
449 C2*(R(elem,i,k)*S(elem,k,j) + R(elem,j,k)*S(elem,k,i)) +
450 C3 *R(elem,i,k)*R(elem,j,k);
451 if (i==j)
452 {
453 for (int k=0; k<Objet_U::dimension; k++)
454 for (int l=0; l<Objet_U::dimension; l++)
455 ReNL(elem,i,j) += 1./3.*(C1*S(elem,k,l)*S(elem,k,l)+C3*R(elem,k,l)*R(elem,k,l));
456 }
457 }
458 }
459
460 Debog::verifier("calcul_tenseur_Re_elem ReNL", ReNL);
461
462 return ReNL;
463}
464
465DoubleTab Modele_Lam_Bremhorst_VEF::calcul_tenseur_Re_elem_BiK(const Discretisation_base& dis,
466 const Domaine_dis_base& domaine_dis,
467 const DoubleTab& G, const DoubleTab& S,
468 const DoubleTab& R, const DoubleTab& Sn,
469 const Champ_base& K, const Champ_base& Eps) const
470{
471 int nelem = G.dimension(0);
472
473 OWN_PTR(Champ_Fonc_base) K_elem;
474 OWN_PTR(Champ_Fonc_base) Eps_elem;
475 Noms noms(1), unites(1), noms2(1), unites2(1);
476 noms[0]="K";
477 noms2[0]="eps";
478 unites[0]="m2/s2";
479 unites2[0]="m2/s3";
480 const Motcle wtype= "champ_elem";
481 dis.discretiser_champ(wtype,domaine_dis,multi_scalaire,noms,unites,1,0.,K_elem);
482 dis.discretiser_champ(wtype,domaine_dis,multi_scalaire,noms2,unites2,1,0.,Eps_elem);
483
484 K_elem->affecter(K);
485 Eps_elem->affecter(Eps);
486 DoubleTab tab_K = K_elem->valeurs();
487 DoubleTab tab_Eps = Eps_elem->valeurs();
488
489 DoubleTab ReNL;
490 const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF,domaine_dis);
491 init_tenseur_elem(ReNL,domaine_VEF,2);
492
493 ReNL = 0.;
494 for (int elem=0; elem<nelem; elem++)
495 {
496 double kseps;
497 if (tab_Eps(elem,0) <= BR_EPS)
498 {
499 kseps = tab_K(elem,0)/BR_EPS;
500 }
501 else
502 {
503 kseps = tab_K(elem,0)/tab_Eps(elem,0);
504 }
505 double Cmu = (2./3.)/(A1+sqrt(0.5)*Sn(elem)*kseps);
506 double C1 = CNL1/((CNL4+CNL5*pow(sqrt(0.5)*Sn(elem)*kseps,3.))*Cmu)* kseps;
507 double C2 = CNL2/((CNL4+CNL5*pow(sqrt(0.5)*Sn(elem)*kseps,3.))*Cmu)* kseps;
508 double C3 = CNL3/((CNL4+CNL5*pow(sqrt(0.5)*Sn(elem)*kseps,3.))*Cmu)* kseps;
509 for (int i=0; i<Objet_U::dimension; i++)
510 for (int j=0; j<Objet_U::dimension; j++)
511 {
512 ReNL(elem,i,j) = S(elem,i,j);
513 for (int k=0; k<Objet_U::dimension; k++)
514 ReNL(elem,i,j) -= C1* S(elem,i,k)*S(elem,k,j) +
515 C2*(R(elem,i,k)*S(elem,k,j) + R(elem,j,k)*S(elem,k,i)) +
516 C3 *R(elem,i,k)*R(elem,j,k);
517 if (i==j)
518 {
519 for (int k=0; k<Objet_U::dimension; k++)
520 for (int l=0; l<Objet_U::dimension; l++)
521 ReNL(elem,i,j) += 1./3.*(C1*S(elem,k,l)*S(elem,k,l)+C3*R(elem,k,l)*R(elem,k,l));
522 }
523 }
524 }
525
526 return ReNL;
527}
528
529
530DoubleTab Modele_Lam_Bremhorst_VEF::calcul_tenseur_Re_shih(const Discretisation_base& dis,
531 const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis,
532 const DoubleTab& G,
533 const Champ_base& K_Eps) const
534{
535
536 const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF,domaine_dis);
537 int nelem = G.dimension(0);
538
539 DoubleTab S, R;
540 init_tenseur_elem(S,domaine_VEF,2);
541 init_tenseur_elem(R,domaine_VEF,2);
542
543 for (int elem=0; elem<nelem; elem++)
544 for (int i=0; i<Objet_U::dimension; i++)
545 for (int j=0; j<Objet_U::dimension; j++)
546 {
547 S(elem,i,j) = 0.5*(G(elem,i,j) + G(elem,j,i));
548 R(elem,i,j) = 0.5*(G(elem,i,j) - G(elem,j,i));
549 if (i==j)
550 for (int k=0; k<Objet_U::dimension; k++)
551 {
552 S(elem,i,j) -= 1/3 * G(elem,k,k);
553 }
554 }
555
556 DoubleTab Sn = calcul_norme_elem(domaine_VEF,S);
557
558 DoubleTab Rn = calcul_norme_elem(domaine_VEF,R);
559
560 DoubleTab ReNL = calcul_tenseur_Re_elem_shih(dis,domaine_dis,G,S,R,Sn,Rn,K_Eps);
561
562 return ReNL;
563}
564
565DoubleTab Modele_Lam_Bremhorst_VEF::calcul_tenseur_Re_elem_shih(const Discretisation_base& dis,
566 const Domaine_dis_base& domaine_dis,
567 const DoubleTab& G, const DoubleTab& S,
568 const DoubleTab& R, const DoubleTab& Sn,const DoubleTab& Rn,
569 const Champ_base& K_Eps) const
570{
571 int nelem = G.dimension(0);
572
573 OWN_PTR(Champ_Fonc_base) K_Eps_elem;
574 Noms noms(2), unites(2);
575 noms[0]="K";
576 noms[1]="eps";
577 unites[0]="m2/s2";
578 unites[1]="m2/s3";
579 const Motcle wtype= "champ_elem";
580 dis.discretiser_champ(wtype,domaine_dis,multi_scalaire,noms,unites,2,0.,K_Eps_elem);
581 K_Eps_elem->affecter(K_Eps);
582 DoubleTab tab_K_Eps = K_Eps_elem->valeurs();
583
584 DoubleTab ReNL;
585 const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF,domaine_dis);
586 init_tenseur_elem(ReNL,domaine_VEF,2);
587
588 ReNL = 0.;
589 for (int elem=0; elem<nelem; elem++)
590 {
591 double kseps;
592 if (tab_K_Eps(elem,1) <= BR_EPS)
593 {
594 kseps = tab_K_Eps(elem,0)/BR_EPS;
595 }
596 else
597 {
598 kseps = tab_K_Eps(elem,0)/tab_K_Eps(elem,1);
599 }
600 double Cmu = (1.)/(6.5+sqrt(pow(Sn(elem),2)+pow(Rn(elem),2))*1.8*kseps);
601 double C2 = sqrt(1-pow(3*Cmu*Sn(elem)*kseps,2))/(1+6.*Sn(elem)*kseps*Rn(elem)*kseps);
602 for (int i=0; i<Objet_U::dimension; i++)
603 for (int j=0; j<Objet_U::dimension; j++)
604 {
605 ReNL(elem,i,j) = 2*S(elem,i,j);
606 for (int k=0; k<Objet_U::dimension; k++)
607 ReNL(elem,i,j) -= 2*C2*kseps/Cmu*(-S(elem,i,k)*R(elem,k,j)+S(elem,k,j)*R(elem,i,k));
608 }
609 }
610
611 return ReNL;
612}
613
614// Calcul de la norme d'un tenseur aux elements
615DoubleTab Modele_Lam_Bremhorst_VEF::calcul_norme_elem(const Domaine_VEF& domaine_VEF,const DoubleTab Tenseur) const
616{
617 int dim_tens = 0;
618 DoubleTab Tnorme;
619 init_tenseur_elem(Tnorme,domaine_VEF,dim_tens);
620 int nb_elems = Tnorme.dimension(0);
621 for (int num_elem=0; num_elem<nb_elems; num_elem++)
622 {
623 for (int i=0; i<Objet_U::dimension; i++)
624 for (int j=0; j<Objet_U::dimension; j++)
625 Tnorme(num_elem) += Tenseur(num_elem,i,j)*Tenseur(num_elem,i,j);
626 Tnorme(num_elem) = sqrt(Tnorme(num_elem));
627 }
628 return Tnorme;
629}
630
631void Modele_Lam_Bremhorst_VEF::lire_distance_paroi( )
632{
633
634 // PQ : 25/02/04 recuperation de la distance a la paroi dans Wall_length.xyz
635
636 const Domaine_VEF& domaine_VEF = le_dom_VEF.valeur();
637 DoubleTab& wall_length = BR_wall_length_->valeurs();
638 wall_length=-1.;
639// return;
640
641 LecFicDiffuse fic;
642 fic.set_bin(1);
643 if(!fic.ouvrir(nom_fic))
644 {
645 Cerr << " File " <<nom_fic<< " doesn't exist. To generate it, please, refer to html.doc (Distance_paroi) " << finl;
646 Process::exit();
647 }
648
649 Noms nom_paroi;
650
651 fic >> nom_paroi;
652
653 if(je_suis_maitre())
654 {
655
656 Cerr << "Recall : " <<nom_fic<< " obtained from boundaries nammed : "<< finl;
657 for (int b=0; b<nom_paroi.size(); b++)
658 {
659 Cerr << nom_paroi[b]<< finl;
660 //test pour s'assurer de la coherence de Wall_length.xyz avec le jeu de donnees :
661 domaine_VEF.rang_frontiere(nom_paroi[b]);
662 }
663 }
664 EcritureLectureSpecial::lecture_special(domaine_VEF, fic, wall_length);
665
666}
667
668DoubleTab& Modele_Lam_Bremhorst_VEF::Calcul_Cmu(DoubleTab& Cmu,
669 const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis,
670 const DoubleTab& vitesse, const DoubleTab& K_Eps, const double EPS_MIN) const
671{
672
673 const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF,domaine_dis);
674 const Domaine_Cl_VEF& domaine_Cl_VEF = ref_cast(Domaine_Cl_VEF,domaine_Cl_dis);
675
676 DoubleTab gradient_elem;
677 init_tenseur_elem(gradient_elem,domaine_VEF,2);
678 Champ_P1NC::calcul_gradient(vitesse,gradient_elem,domaine_Cl_VEF);
679
680 DoubleTab S_elem;
681 init_tenseur_elem(S_elem,domaine_VEF,2);
682 int nelem = S_elem.dimension(0);
683 for (int elem=0; elem<nelem; elem++)
684 for (int i=0; i<dimension; i++)
685 for (int j=0; j<dimension; j++)
686 S_elem(elem,i,j) = gradient_elem(elem,i,j) + gradient_elem(elem,j,i);
687 S_elem.echange_espace_virtuel();
688
689 DoubleTab S_face;
690 init_tenseur_face(S_face,domaine_VEF,2);
691
692 calcul_tenseur_face(S_face,S_elem,domaine_VEF,domaine_Cl_VEF);
693
694 int nfaces = S_face.dimension(0);
695 DoubleTab Snorme_face;
696 domaine_VEF.creer_tableau_faces(Snorme_face);
697
698 for (int face=0; face<nfaces; face++)
699 {
700 double somme = 0.;
701 for (int i=0; i<Objet_U::dimension; i++)
702 for (int j=0; j<Objet_U::dimension; j++)
703 somme += S_face(face,i,j)*S_face(face,i,j);
704 Snorme_face(face) = sqrt(somme);
705 }
706
707 for (int face=0; face<nfaces; face++)
708 {
709 // Definition d'un Cmu minimum base sur EPS_MIN = 1e-10
710 if (K_Eps(face,1) <= EPS_MIN)
711 Cmu[face] = (2./3.)/(A1+sqrt(0.5)*Snorme_face[face]*K_Eps(face,0)/BR_EPS);
712 else
713 Cmu[face] = (2./3.)/(A1+sqrt(0.5)*Snorme_face[face]*K_Eps(face,0)/K_Eps(face,1));
714 }
715
716 return Cmu;
717}
718
719DoubleTab& Modele_Lam_Bremhorst_VEF::Calcul_Cmu_Paroi(DoubleTab& Cmu,
720 const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis,
721 const DoubleTab& visco, const DoubleTab& visco_turb,
722 const DoubleTab& loi_paroi,const int idt,
723 const DoubleTab& vitesse, const DoubleTab& K_Eps, const double EPS_MIN) const
724{
725
726 const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF,domaine_dis);
727 const Domaine_Cl_VEF& domaine_Cl_VEF = ref_cast(Domaine_Cl_VEF,domaine_Cl_dis);
728
729 DoubleTab gradient_elem;
730 init_tenseur_elem(gradient_elem,domaine_VEF,2);
731 Champ_P1NC::calcul_gradient(vitesse,gradient_elem,domaine_Cl_VEF);
732
733 if (idt>0)
734 Champ_P1NC::calcul_duidxj_paroi(gradient_elem,visco,visco_turb,loi_paroi,domaine_Cl_VEF);
735 DoubleTab S_elem;
736 init_tenseur_elem(S_elem,domaine_VEF,2);
737 int nelem = S_elem.dimension(0);
738 for (int elem=0; elem<nelem; elem++)
739 for (int i=0; i<dimension; i++)
740 for (int j=0; j<dimension; j++)
741 S_elem(elem,i,j) = gradient_elem(elem,i,j) + gradient_elem(elem,j,i);
742 S_elem.echange_espace_virtuel();
743
744 DoubleTab S_face;
745 init_tenseur_face(S_face,domaine_VEF,2);
746
747 calcul_tenseur_face(S_face,S_elem,domaine_VEF,domaine_Cl_VEF);
748
749 int nfaces = S_face.dimension(0);
750 DoubleTab Snorme_face;
751 domaine_VEF.creer_tableau_faces(Snorme_face);
752
753 for (int face=0; face<nfaces; face++)
754 {
755 double somme = 0.;
756 for (int i=0; i<Objet_U::dimension; i++)
757 for (int j=0; j<Objet_U::dimension; j++)
758 somme += S_face(face,i,j)*S_face(face,i,j);
759 Snorme_face(face) = sqrt(somme);
760 }
761
762 for (int face=0; face<nfaces; face++)
763 {
764 // Definition d'un Cmu minimum base sur EPS_MIN = 1e-10
765 if (K_Eps(face,1) <= EPS_MIN)
766 Cmu[face] = (2./3.)/(A1+sqrt(0.5)*Snorme_face[face]*K_Eps(face,0)/BR_EPS);
767 else
768 Cmu[face] = (2./3.)/(A1+sqrt(0.5)*Snorme_face[face]*K_Eps(face,0)/K_Eps(face,1));
769 }
770
771 // TMA nettoyage de restes de debug
772 // Cerr<<Cmu.mp_min_vect()<<" Cmuuuuuuuuuuuuuuuuuuuu " <<Cmu.mp_max_vect()<<finl;
773 return Cmu;
774}
775
776
777DoubleTab& Modele_Lam_Bremhorst_VEF::Calcul_Cmu_BiK(DoubleTab& Cmu,
778 const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis,
779 const DoubleTab& vitesse, const DoubleTab& K, const DoubleTab& Eps, const double EPS_MIN) const
780{
781
782 const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF,domaine_dis);
783 const Domaine_Cl_VEF& domaine_Cl_VEF = ref_cast(Domaine_Cl_VEF,domaine_Cl_dis);
784
785 DoubleTab gradient_elem;
786 init_tenseur_elem(gradient_elem,domaine_VEF,2);
787 Champ_P1NC::calcul_gradient(vitesse,gradient_elem,domaine_Cl_VEF);
788
789 DoubleTab S_elem;
790 init_tenseur_elem(S_elem,domaine_VEF,2);
791 int nelem = S_elem.dimension(0);
792 for (int elem=0; elem<nelem; elem++)
793 for (int i=0; i<dimension; i++)
794 for (int j=0; j<dimension; j++)
795 S_elem(elem,i,j) = gradient_elem(elem,i,j) + gradient_elem(elem,j,i);
796 S_elem.echange_espace_virtuel();
797
798 DoubleTab S_face;
799 init_tenseur_face(S_face,domaine_VEF,2);
800
801 calcul_tenseur_face(S_face,S_elem,domaine_VEF,domaine_Cl_VEF);
802
803 int nfaces = S_face.dimension(0);
804 DoubleTab Snorme_face;
805 domaine_VEF.creer_tableau_faces(Snorme_face);
806
807 for (int face=0; face<nfaces; face++)
808 {
809 double somme = 0.;
810 for (int i=0; i<Objet_U::dimension; i++)
811 for (int j=0; j<Objet_U::dimension; j++)
812 somme += S_face(face,i,j)*S_face(face,i,j);
813 Snorme_face(face) = sqrt(somme);
814 }
815
816 for (int face=0; face<nfaces; face++)
817 {
818 // Definition d'un Cmu minimum base sur EPS_MIN = 1e-10
819 if (Eps(face) <= EPS_MIN)
820 Cmu[face] = (2./3.)/(A1+sqrt(0.5)*Snorme_face[face]*K(face)/BR_EPS);
821 else
822 Cmu[face] = (2./3.)/(A1+sqrt(0.5)*Snorme_face[face]*K(face)/Eps(face));
823 }
824
825 return Cmu;
826}
827
828DoubleTab& Modele_Lam_Bremhorst_VEF::Calcul_Cmu_Paroi_BiK(DoubleTab& Cmu,
829 const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis,
830 const DoubleTab& visco, const DoubleTab& visco_turb,
831 const DoubleTab& loi_paroi,const int idt,
832 const DoubleTab& vitesse, const DoubleTab& K, const DoubleTab& Eps, const double EPS_MIN) const
833{
834
835 const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF,domaine_dis);
836 const Domaine_Cl_VEF& domaine_Cl_VEF = ref_cast(Domaine_Cl_VEF,domaine_Cl_dis);
837
838 DoubleTab gradient_elem;
839 init_tenseur_elem(gradient_elem,domaine_VEF,2);
840 Champ_P1NC::calcul_gradient(vitesse,gradient_elem,domaine_Cl_VEF);
841
842 if (idt>0)
843 Champ_P1NC::calcul_duidxj_paroi(gradient_elem,visco,visco_turb,loi_paroi,domaine_Cl_VEF);
844 DoubleTab S_elem;
845 init_tenseur_elem(S_elem,domaine_VEF,2);
846 int nelem = S_elem.dimension(0);
847 for (int elem=0; elem<nelem; elem++)
848 for (int i=0; i<dimension; i++)
849 for (int j=0; j<dimension; j++)
850 S_elem(elem,i,j) = gradient_elem(elem,i,j) + gradient_elem(elem,j,i);
851 S_elem.echange_espace_virtuel();
852
853 DoubleTab S_face;
854 init_tenseur_face(S_face,domaine_VEF,2);
855
856 calcul_tenseur_face(S_face,S_elem,domaine_VEF,domaine_Cl_VEF);
857
858 int nfaces = S_face.dimension(0);
859 DoubleTab Snorme_face;
860 domaine_VEF.creer_tableau_faces(Snorme_face);
861
862 for (int face=0; face<nfaces; face++)
863 {
864 double somme = 0.;
865 for (int i=0; i<Objet_U::dimension; i++)
866 for (int j=0; j<Objet_U::dimension; j++)
867 somme += S_face(face,i,j)*S_face(face,i,j);
868 Snorme_face(face) = sqrt(somme);
869 }
870
871 for (int face=0; face<nfaces; face++)
872 {
873 // Definition d'un Cmu minimum base sur EPS_MIN = 1e-10
874 if (Eps(face) <= EPS_MIN)
875 Cmu[face] = (2./3.)/(A1+sqrt(0.5)*Snorme_face[face]*K(face)/BR_EPS);
876 else
877 Cmu[face] = (2./3.)/(A1+sqrt(0.5)*Snorme_face[face]*K(face)/Eps(face));
878 }
879
880 // TMA nettoyage de restes de debug
881 // Cerr<<Cmu.mp_min_vect()<<" Cmuuuuuuuuuuuuuuuuuuuu " <<Cmu.mp_max_vect()<<finl;
882 return Cmu;
883}
884
885
886DoubleTab& Modele_Lam_Bremhorst_VEF::Calcul_Fmu_BiK( DoubleTab& Fmu,const Domaine_dis_base& domaine_dis,const Domaine_Cl_dis_base& domaine_Cl_dis,const DoubleTab& K_Bas_Re,const DoubleTab& eps_Bas_Re,const Champ_Don_base& ch_visco ) const
887{
888 double visco=-1;
889 const DoubleTab& tab_visco=ch_visco.valeurs();
890 int is_visco_const=sub_type(Champ_Uniforme,ch_visco);
891 if (is_visco_const)
892 visco=tab_visco(0,0);
893 const Domaine_VEF& le_dom = ref_cast(Domaine_VEF,domaine_dis);
894 const DoubleTab& wall_length = BR_wall_length_->valeurs();
895 DoubleTab wall_length_face(0);
896 le_dom.creer_tableau_faces(wall_length_face);
897 int nb_faces = le_dom.nb_faces();
898 const IntTab& face_voisins = le_dom.face_voisins();
899 int num_face;
900 double Rey,Re;
901
902 // Calcul de la distance a la paroi aux faces
903 for (num_face=0; num_face<nb_faces ; num_face++)
904 {
905 int elem0 = face_voisins(num_face,0);
906 int elem1 = face_voisins(num_face,1);
907 if (elem1!=-1)
908 wall_length_face(num_face) = 0.5*wall_length(elem0)+0.5*wall_length(elem1);
909 else
910 wall_length_face(num_face) = wall_length(elem0);
911 }
912
913 for (num_face=0; num_face <nb_faces; num_face ++ )
914 {
915 if (visco>BR_EPS && eps_Bas_Re(num_face)>BR_EPS && K_Bas_Re(num_face))
916 {
917 Re = (K_Bas_Re(num_face)*K_Bas_Re(num_face))/(visco*eps_Bas_Re(num_face));
918 Rey = wall_length_face(num_face)*sqrt(K_Bas_Re(num_face))/visco;
919 Fmu(num_face) = (1. - exp(-0.0165*Rey-BR_EPS))*(1. - exp(-0.0165*Rey-BR_EPS))*(1.+20.5/(Re+BR_EPS));
920
921 }
922 else
923 {
924 //Cerr << " visco " << visco << " Eps " << K_eps_Bas_Re(num_face,1) << " K " << K_eps_Bas_Re(num_face,0) <<finl;
925 Fmu(num_face) = 1.;
926 }
927 if (Fmu(num_face) > 1. + 1.e-8)
928 {
929 //Cerr << " Fmu=" << Fmu(num_face) << " Eps=" << K_eps_Bas_Re(num_face,1) << " K=" << K_eps_Bas_Re(num_face,0) << " Re=" << Re << " Rey=" << Rey <<finl;
930 Fmu(num_face) = 1.;
931 //Cerr << " On force Fmu a 1 " << finl;
932 //exit();
933 }
934 }
935 Cerr<< Fmu.mp_min_vect()<< " Fmu "<<Fmu.mp_max_vect()<<finl;
936
937 return Fmu;
938}
939
940
941DoubleTab& Modele_Lam_Bremhorst_VEF::Calcul_F2_BiK( DoubleTab& F2, DoubleTab& Deb, const Domaine_dis_base& domaine_dis,const DoubleTab& K_Bas_Re,const DoubleTab& eps_Bas_Re,const Champ_base& ch_visco ) const
942{
943 double visco=-1;
944 const DoubleTab& tab_visco=ch_visco.valeurs();
945 int is_visco_const=sub_type(Champ_Uniforme,ch_visco);
946 if (is_visco_const)
947 visco=tab_visco(0,0);
948 const Domaine_VEF& le_dom = ref_cast(Domaine_VEF,domaine_dis);
949 int nb_faces = le_dom.nb_faces();
950 int num_face;
951 double Re;
952
953 for (num_face=0; num_face<nb_faces ; num_face++)
954 {
955 if (!is_visco_const)
956 {
957 int elem0 = le_dom.face_voisins(num_face,0);
958 int elem1 = le_dom.face_voisins(num_face,1);
959 if (elem1!=-1)
960 {
961 visco = tab_visco(elem0)*le_dom.volumes(elem0)+tab_visco(elem1)*le_dom.volumes(elem1);
962 visco /= le_dom.volumes(elem0) + le_dom.volumes(elem1);
963 }
964 else
965 visco = tab_visco(elem0);
966
967 }
968 if (visco>BR_EPS && eps_Bas_Re(num_face)>BR_EPS)
969 {
970 Re = (K_Bas_Re(num_face)*K_Bas_Re(num_face))/(visco*eps_Bas_Re(num_face));
971 F2[num_face] = 1. - (exp(-Re*Re));
972 }
973 else
974 F2[num_face] = 1.;
975
976 }
977 Cerr<<F2.mp_min_vect()<<" F2 "<<F2.mp_max_vect()<<finl;
978 return F2;
979}
980
981DoubleTab& Modele_Lam_Bremhorst_VEF::Calcul_F1_BiK( DoubleTab& F1, const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis, const DoubleTab& P, const DoubleTab& K_Bas_Re, const DoubleTab& eps_Bas_Re,const Champ_base& ch_visco) const
982{
983
984 double visco=-1;
985 const DoubleTab& tab_visco=ch_visco.valeurs();
986 int is_visco_const=sub_type(Champ_Uniforme,ch_visco);
987 if (is_visco_const)
988 visco=tab_visco(0,0);
989 const Domaine_VEF& le_dom = ref_cast(Domaine_VEF,domaine_dis);
990 const DoubleTab& wall_length = BR_wall_length_->valeurs();
991 DoubleTab wall_length_face(0);
992 le_dom.creer_tableau_faces(wall_length_face);
993 DoubleTab Fmu_loc(0);
994 le_dom.creer_tableau_faces(Fmu_loc);
995 int nb_faces = le_dom.nb_faces();
996 const IntTab& face_voisins = le_dom.face_voisins();
997 int num_face;
998 double Rey,Re;
999
1000 // Calcul de la distance a la paroi aux faces
1001 for (num_face=0; num_face<nb_faces ; num_face++)
1002 {
1003 int elem0 = face_voisins(num_face,0);
1004 int elem1 = face_voisins(num_face,1);
1005 if (elem1!=-1)
1006 wall_length_face(num_face) = 0.5*wall_length(elem0)+0.5*wall_length(elem1);
1007 else
1008 wall_length_face(num_face) = wall_length(elem0);
1009 }
1010
1011 for (num_face=0; num_face <nb_faces; num_face ++ )
1012 {
1013 if (visco>BR_EPS && eps_Bas_Re(num_face)>BR_EPS && K_Bas_Re(num_face))
1014 {
1015 Re = (K_Bas_Re(num_face)*K_Bas_Re(num_face))/(visco*eps_Bas_Re(num_face));
1016 Rey = wall_length_face(num_face)*sqrt(K_Bas_Re(num_face))/visco;
1017 Fmu_loc(num_face) = (1. - exp(-0.0165*Rey-BR_EPS))*(1. - exp(-0.0165*Rey-BR_EPS))*(1.+20.5/(Re+BR_EPS));
1018 }
1019 else
1020 Fmu_loc(num_face) = 1.;
1021
1022 if (Fmu_loc(num_face) > 1. + 1.e-8)
1023 {
1024 Fmu_loc(num_face) = 1.;
1025 //Cerr << " On force Fmu_loc a 1 " << finl;
1026 //exit();
1027 }
1028 }
1029 Cerr<< Fmu_loc.mp_min_vect()<< " Fmu_loc "<<Fmu_loc.mp_max_vect()<<finl;
1030
1031 // Calcul de F1
1032 for (num_face=0; num_face <nb_faces; num_face ++ )
1033 {
1034 if (Fmu_loc(num_face) > BR_EPS && (visco>BR_EPS && eps_Bas_Re(num_face)>BR_EPS && K_Bas_Re(num_face)))
1035 F1[num_face] = 1.+ (0.05/(Fmu_loc(num_face)+BR_EPS))*(0.05/(Fmu_loc(num_face)+BR_EPS))*(0.05/(Fmu_loc(num_face)+BR_EPS));
1036 else
1037 F1[num_face] = 1.;
1038 }
1039 Cerr<< F1.mp_min_vect()<< " F1 "<<F1.mp_max_vect()<<finl;
1040 return F1;
1041}
1042
1043
1044
1045DoubleTab& Modele_Lam_Bremhorst_VEF::Calcul_E_BiK(DoubleTab& E,const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis, const DoubleTab& transporte,const DoubleTab& K_Bas_Re,const DoubleTab& eps_Bas_Re,const Champ_Don_base& ch_visco, const DoubleTab& visco_turb ) const
1046{
1047 E = 0;
1048 return E;
1049}
1050
1051DoubleTab& Modele_Lam_Bremhorst_VEF::Calcul_D_BiK(DoubleTab& D,const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis,
1052 const DoubleTab& vitesse,const DoubleTab& K_Bas_Re,const DoubleTab& eps_Bas_Re, const Champ_Don_base& ch_visco ) const
1053{
1054 D = 0;
1055 return D;
1056}
1057
1058
1059
1060
1061
1062
1063
1064
1065
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_Fonc_base
classe Champ_Fonc_base Classe de base des champs qui sont fonction d'une grandeur calculee
Definition Champ_Fonc_base.h:37
Champ_P1NC::calcul_gradient
static DoubleTab & calcul_gradient(const DoubleTab &, DoubleTab &, const Domaine_Cl_VEF &)
Definition Champ_P1NC.cpp:915
Champ_P1NC::calcul_duidxj_paroi
static DoubleTab & calcul_duidxj_paroi(DoubleTab &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const Domaine_Cl_VEF &)
Definition Champ_P1NC.cpp:923
Champ_Proto::valeurs
virtual DoubleTab & valeurs()=0
Champ_Uniforme
classe Champ_Uniforme Represente un champ constant dans l'espace et dans le temps.
Definition Champ_Uniforme.h:26
Champ_base
classe Champ_base Cette classe est la base de la hierarchie des champs.
Definition Champ_base.h:43
Cond_lim
classe Cond_lim Classe generique servant a representer n'importe quelle classe
Definition Cond_lim.h:31
Conds_lim
classe Conds_lim Cette classe represente un vecteur de conditions aux limites.
Definition Conds_lim.h:32
Debog::verifier
static void verifier(const char *const msg, double)
Definition Debog.cpp:21
Discretisation_base
classe Discretisation_base Cette classe represente un schema de discretisation en espace,...
Definition Discretisation_base.h:45
Discretisation_base::discretiser_champ
void discretiser_champ(const Motcle &directive, const Domaine_dis_base &z, const Nom &nom, const Nom &unite, int nb_comp, int nb_pas_dt, double temps, OWN_PTR(Champ_Inc_base)&champ, const Nom &sous_type=NOM_VIDE) const
Definition Discretisation_base.cpp:59
Domaine_32_64::creer_tableau_elements
virtual void creer_tableau_elements(Array_base &, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT) const
creation d'un tableau parallele de valeurs aux elements.
Definition Domaine.cpp:851
Domaine_32_64::nb_elem
int_t nb_elem() const
Definition Domaine.h:131
Domaine_Cl_VEF
Definition Domaine_Cl_VEF.h:40
Domaine_Cl_dis_base
classe Domaine_Cl_dis_base Les objets Domaine_Cl_dis_base representent les conditions aux limites
Definition Domaine_Cl_dis_base.h:37
Domaine_Cl_dis_base::les_conditions_limites
const Cond_lim & les_conditions_limites(int) const
Renvoie la i-ieme condition aux limites.
Definition Domaine_Cl_dis_base.cpp:386
Domaine_VEF
class Domaine_VEF
Definition Domaine_VEF.h:54
Domaine_VF::nb_faces
int nb_faces() const
renvoie le nombre global de faces.
Definition Domaine_VF.h:471
Domaine_VF::creer_tableau_faces
void creer_tableau_faces(Array_base &, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT) const
Definition Domaine_VF.cpp:881
Domaine_VF::volumes
double volumes(int i) const
Definition Domaine_VF.h:113
Domaine_VF::premiere_face_int
int premiere_face_int() const
une face est interne ssi elle separe deux elements.
Definition Domaine_VF.h:463
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
classe Domaine_dis_base Cette classe est la base de la hierarchie des domaines discretisees.
Definition Domaine_dis_base.h:39
Domaine_dis_base::rang_frontiere
int rang_frontiere(const Nom &)
Definition Domaine_dis_base.cpp:95
Domaine_dis_base::domaine
const Domaine & domaine() const
Definition Domaine_dis_base.h:46
EcritureLectureSpecial::lecture_special
static void lecture_special(Champ_base &ch, Entree &fich)
simple appel a EcritureLectureSpecial::lecture_special (const Domaine_VF& zvf,Entree& fich,...
Definition EcritureLectureSpecial.cpp:311
Entree
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
Front_VF
class Front_VF
Definition Front_VF.h:36
Front_VF::nb_faces
int nb_faces() const
Definition Front_VF.h:53
Front_VF::num_premiere_face
int num_premiere_face() const
Definition Front_VF.h:63
LecFicDiffuse
Cette classe implemente les operateurs et les methodes virtuelles de la classe EFichier de la facon s...
Definition LecFicDiffuse.h:32
LecFicDiffuse::ouvrir
int ouvrir(const char *name, IOS_OPEN_MODE mode=ios::in) override
Ouverture du fichier.
Definition LecFicDiffuse.cpp:61
Lec_Diffuse_base::set_bin
void set_bin(bool bin) override
appelle get_entree_master().
Definition Lec_Diffuse_base.cpp:129
Modele_Fonc_Bas_Reynolds_Base
Definition Modele_Fonc_Bas_Reynolds_Base.h:38
Modele_Fonc_Bas_Reynolds_Base::nom_fic
Nom nom_fic
Definition Modele_Fonc_Bas_Reynolds_Base.h:114
Modele_Fonc_Bas_Reynolds_Base::OWN_PTR
OWN_PTR(Champ_Fonc_base) BR_wall_length_
Modele_Lam_Bremhorst_VEF
Definition Modele_Lam_Bremhorst_VEF.h:45
Modele_Lam_Bremhorst_VEF::calcul_tenseur_Re_BiK
bool calcul_tenseur_Re_BiK(const DoubleTab &, const DoubleTab &, DoubleTab &) const override
Definition Modele_Lam_Bremhorst_VEF.cpp:380
Modele_Lam_Bremhorst_VEF::Calcul_D
DoubleTab & Calcul_D(DoubleTab &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, const DoubleTab &, const DoubleTab &, const Champ_Don_base &) const override
Definition Modele_Lam_Bremhorst_VEF.cpp:83
Modele_Lam_Bremhorst_VEF::calcul_tenseur_face
virtual DoubleTab & calcul_tenseur_face(DoubleTab &, const DoubleTab &, const Domaine_VEF &, const Domaine_Cl_VEF &) const
Definition Modele_Lam_Bremhorst_VEF.cpp:298
Modele_Lam_Bremhorst_VEF::calcul_norme_elem
virtual DoubleTab calcul_norme_elem(const Domaine_VEF &, const DoubleTab) const
Definition Modele_Lam_Bremhorst_VEF.cpp:615
Modele_Lam_Bremhorst_VEF::Calcul_F1
DoubleTab & Calcul_F1(DoubleTab &F1, const Domaine_dis_base &domaine_dis, const Domaine_Cl_dis_base &domaine_Cl_dis, const DoubleTab &P, const DoubleTab &K_eps_Bas_Re, const Champ_base &ch_visco) const override
Definition Modele_Lam_Bremhorst_VEF.cpp:105
Modele_Lam_Bremhorst_VEF::Calcul_Fmu_BiK
DoubleTab & Calcul_Fmu_BiK(DoubleTab &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, const DoubleTab &, const DoubleTab &, const Champ_Don_base &) const override
Definition Modele_Lam_Bremhorst_VEF.cpp:886
Modele_Lam_Bremhorst_VEF::Calcul_E_BiK
DoubleTab & Calcul_E_BiK(DoubleTab &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const Champ_Don_base &, const DoubleTab &) const override
Definition Modele_Lam_Bremhorst_VEF.cpp:1045
Modele_Lam_Bremhorst_VEF::Calcul_is_Reynolds_stress_isotrope
int Calcul_is_Reynolds_stress_isotrope() const override
Definition Modele_Lam_Bremhorst_VEF.cpp:58
Modele_Lam_Bremhorst_VEF::init_tenseur_face
void init_tenseur_face(DoubleTab &, const Domaine_VEF &, const int) const
Definition Modele_Lam_Bremhorst_VEF.cpp:284
Modele_Lam_Bremhorst_VEF::Calcul_F2_BiK
DoubleTab & Calcul_F2_BiK(DoubleTab &, DoubleTab &, const Domaine_dis_base &, const DoubleTab &, const DoubleTab &, const Champ_base &) const override
Definition Modele_Lam_Bremhorst_VEF.cpp:941
Modele_Lam_Bremhorst_VEF::calcul_tenseur_Re_elem
virtual DoubleTab calcul_tenseur_Re_elem(const Discretisation_base &dis, const Domaine_dis_base &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const Champ_base &K_Eps) const
Definition Modele_Lam_Bremhorst_VEF.cpp:405
Modele_Lam_Bremhorst_VEF::set_param
virtual void set_param(Param &param) const override
Definition Modele_Lam_Bremhorst_VEF.cpp:50
Modele_Lam_Bremhorst_VEF::calcul_tenseur_Re
bool calcul_tenseur_Re(const DoubleTab &, const DoubleTab &, DoubleTab &) const override
Definition Modele_Lam_Bremhorst_VEF.cpp:355
Modele_Lam_Bremhorst_VEF::mettre_a_jour
void mettre_a_jour(double) override
Definition Modele_Lam_Bremhorst_VEF.cpp:263
Modele_Lam_Bremhorst_VEF::lire_distance_paroi
void lire_distance_paroi() override
Definition Modele_Lam_Bremhorst_VEF.cpp:631
Modele_Lam_Bremhorst_VEF::Calcul_F2
DoubleTab & Calcul_F2(DoubleTab &, DoubleTab &, const Domaine_dis_base &, const DoubleTab &, const Champ_base &) const override
Definition Modele_Lam_Bremhorst_VEF.cpp:169
Modele_Lam_Bremhorst_VEF::Calcul_D_BiK
DoubleTab & Calcul_D_BiK(DoubleTab &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const Champ_Don_base &) const override
Definition Modele_Lam_Bremhorst_VEF.cpp:1051
Modele_Lam_Bremhorst_VEF::Calcul_is_Cmu_constant
int Calcul_is_Cmu_constant() const override
Definition Modele_Lam_Bremhorst_VEF.cpp:63
Modele_Lam_Bremhorst_VEF::calcul_tenseur_Re_elem_BiK
virtual DoubleTab calcul_tenseur_Re_elem_BiK(const Discretisation_base &dis, const Domaine_dis_base &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const Champ_base &K, const Champ_base &Eps) const
Definition Modele_Lam_Bremhorst_VEF.cpp:465
Modele_Lam_Bremhorst_VEF::associer
void associer(const Domaine_dis_base &, const Domaine_Cl_dis_base &) override
Definition Modele_Lam_Bremhorst_VEF.cpp:76
Modele_Lam_Bremhorst_VEF::lire
Entree & lire(const Motcle &, Entree &)
Definition Modele_Lam_Bremhorst_VEF.cpp:68
Modele_Lam_Bremhorst_VEF::Calcul_F1_BiK
DoubleTab & Calcul_F1_BiK(DoubleTab &F1, const Domaine_dis_base &domaine_dis, const Domaine_Cl_dis_base &domaine_Cl_dis, const DoubleTab &P, const DoubleTab &K_Bas_Re, const DoubleTab &eps_Bas_Re, const Champ_base &ch_visco) const override
Definition Modele_Lam_Bremhorst_VEF.cpp:981
Modele_Lam_Bremhorst_VEF::Calcul_Cmu_Paroi_BiK
DoubleTab & Calcul_Cmu_Paroi_BiK(DoubleTab &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const int, const DoubleTab &, const DoubleTab &, const DoubleTab &, const double) const override
Definition Modele_Lam_Bremhorst_VEF.cpp:828
Modele_Lam_Bremhorst_VEF::Calcul_Cmu
DoubleTab & Calcul_Cmu(DoubleTab &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, const DoubleTab &, const DoubleTab &, const double) const override
Definition Modele_Lam_Bremhorst_VEF.cpp:668
Modele_Lam_Bremhorst_VEF::Calcul_Cmu_Paroi
DoubleTab & Calcul_Cmu_Paroi(DoubleTab &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const int, const DoubleTab &, const DoubleTab &, const double) const override
Definition Modele_Lam_Bremhorst_VEF.cpp:719
Modele_Lam_Bremhorst_VEF::calcul_tenseur_Re_elem_shih
DoubleTab calcul_tenseur_Re_elem_shih(const Discretisation_base &, const Domaine_dis_base &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const Champ_base &K_Eps) const
Definition Modele_Lam_Bremhorst_VEF.cpp:565
Modele_Lam_Bremhorst_VEF::Calcul_Fmu
DoubleTab & Calcul_Fmu(DoubleTab &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, const DoubleTab &, const Champ_Don_base &) const override
Definition Modele_Lam_Bremhorst_VEF.cpp:209
Modele_Lam_Bremhorst_VEF::calcul_tenseur_Re_shih
DoubleTab calcul_tenseur_Re_shih(const Discretisation_base &dis, const Domaine_dis_base &domaine_dis, const Domaine_Cl_dis_base &domaine_Cl_dis, const DoubleTab &G, const Champ_base &K_Eps) const
Definition Modele_Lam_Bremhorst_VEF.cpp:530
Modele_Lam_Bremhorst_VEF::init_tenseur_elem
void init_tenseur_elem(DoubleTab &, const Domaine_VEF &, const int) const
Definition Modele_Lam_Bremhorst_VEF.cpp:269
Modele_Lam_Bremhorst_VEF::Calcul_Cmu_BiK
DoubleTab & Calcul_Cmu_BiK(DoubleTab &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const double) const override
Definition Modele_Lam_Bremhorst_VEF.cpp:777
Modele_Lam_Bremhorst_VEF::Calcul_E
DoubleTab & Calcul_E(DoubleTab &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, const DoubleTab &, const DoubleTab &, const Champ_Don_base &, const DoubleTab &) const override
Definition Modele_Lam_Bremhorst_VEF.cpp:90
Motcle
Une chaine de caractere (Nom) en majuscules.
Definition Motcle.h:26
Noms
Un tableau de chaine de caracteres (VECT(Nom)).
Definition Noms.h:26
Objet_U::Entree
friend class Entree
Definition Objet_U.h:76
Objet_U::dimension
static int dimension
Definition Objet_U.h:99
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::printOn
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
Definition Objet_U.cpp:282
Param
Helper class to factorize the readOn method of Objet_U classes.
Definition Param.h:112
Param::ajouter
void ajouter(const char *keyword, const int *value, Param::Nature nat=Param::OPTIONAL)
Register an integer parameter.
Definition Param.cpp:364
Periodique
classe Periodique Cette classe represente une condition aux limites periodique.
Definition Periodique.h:31
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
Sortie
Classe de base des flux de sortie.
Definition Sortie.h:52
TRUSTTab::resize
void resize(_SIZE_ n, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
Definition TRUSTTab.tpp:469
TRUSTTab::dimension
_SIZE_ dimension(int d) const
Definition TRUSTTab.tpp:133
TRUSTVect::mp_max_vect
_TYPE_ mp_max_vect(Mp_vect_options opt=VECT_REAL_ITEMS) const
Definition TRUSTVect.h:158
TRUSTVect::mp_min_vect
_TYPE_ mp_min_vect(Mp_vect_options opt=VECT_REAL_ITEMS) const
Definition TRUSTVect.h:159
TRUSTVect::get_md_vector
virtual const MD_Vector & get_md_vector() const
Definition TRUSTVect.h:123
TRUSTVect::echange_espace_virtuel
virtual void echange_espace_virtuel(IsExchangeBlocking exchange_type=IsExchangeBlocking::DefaultBlocking, const std::string kernel_name="noname")
Definition TRUSTVect.tpp:282