TRUST 1.9.8
HPC thermohydraulic platform
Loading...
Searching...
No Matches
Traitement_particulier_NS_Brech_VEF.cpp
1/****************************************************************************
2* Copyright (c) 2024, 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#include <Traitement_particulier_NS_Brech_VEF.h>
17#include <Navier_Stokes_std.h>
18#include <Fluide_Incompressible.h>
19#include <Probleme_base.h>
20#include <TRUSTTrav.h>
21#include <Periodique.h>
22#include <Champ_P1NC.h>
23
24Implemente_instanciable_sans_constructeur(Traitement_particulier_NS_Brech_VEF,"Traitement_particulier_NS_Brech_VEF",Traitement_particulier_NS_VEF);
25
29
30Sortie& Traitement_particulier_NS_Brech_VEF::printOn(Sortie& is) const
31{
32 return is;
33}
34
35Entree& Traitement_particulier_NS_Brech_VEF::readOn(Entree& is)
36{
37 return is;
38}
39
40/*
41 inline int Traitement_particulier_NS_Brech_VEF::comprend_champ(const Motcle& mot) const
42 {
43 if ( (mot == "Richardson") || (mot == "Pression_porosite") ) return 1 ;
44 return 0 ;
45 }
46
47 int Traitement_particulier_NS_Brech_VEF::a_pour_Champ_Fonc(const Motcle& mot,
48 OBS_PTR(Champ_base)& ch_ref) const
49 {
50 if (mot == "Richardson")
51 {
52 ch_ref = ch_ri;
53 return 1 ;
54 } else if (mot == "Pression_porosite")
55 {
56 ch_ref = ch_p;
57 return 1 ;
58 }
59 return 0 ;
60 }
61*/
62Entree& Traitement_particulier_NS_Brech_VEF::lire(Entree& is)
63{
64 Motcle accouverte = "{" , accfermee = "}" ;
65 Motcle motbidon, motlu;
66
67 is >> motbidon ;
68
69 Motcles les_mots(3);
70 {
71 les_mots[0] = "calcul_flux";
72 les_mots[1] = "Richardson";
73 les_mots[2] = "Pression_porosite";
74 }
75
76
77 if (motbidon == accouverte)
78 {
79 // is >> motlu;
80
81 while(motlu != accfermee)
82 {
83 is >> motlu;
84 int rang=les_mots.search(motlu);
85 switch(rang)
86 {
87 case 0 :
88 {
89 is >> motbidon ;
90 if (motbidon != accouverte)
91 {
92 Cerr << " on attent { et pas " << motbidon << finl ;
93 }
94
95 is >> nb ;
96
97 C_trans.resize(nb,3) ;
98 R_loc.resize(nb,3) ;
99 r_int.resize(nb) ;
100 r_out.resize(nb) ;
101 delta_r.resize(nb) ;
102 delta_teta.resize(nb);
103
104 for (int i=0; i<nb; i++)
105 {
106
107 Cerr << " Lire plan de coupe en x " << finl;
108
109 is >> R_loc(i,0) ;
110
111 is >> R_loc(i,1) ;
112
113 is >> R_loc(i,2) ;
114
115 is >> C_trans(i,0) ;
116
117 is >> C_trans(i,1) ;
118
119 is >> C_trans(i,2) ;
120
121 is >> r_int(i);
122
123 is >> r_out(i);
124
125 is >> delta_r(i);
126
127 is >> delta_teta(i) ;
128 }
129
130 c_flux = 1 ;
131
132 is >> motbidon ;
133 if (motbidon != accfermee)
134 {
135 Cerr << " on attent } et pas " << motbidon << finl ;
136 }
137 break ;
138
139 }
140 case 1 :
141 {
142 Cerr << " Lire Richardson " << finl;
143 const Domaine_dis_base& zdis=mon_equation->inconnue().domaine_dis_base();
144 const Domaine_VEF& domaine_VEF=ref_cast(Domaine_VEF, zdis);
145 // const Probleme_base& pb = mon_equation->probleme();
146 const int nb_faces = domaine_VEF.nb_faces() ;
147
148 ch_ri.associer_domaine_dis_base(zdis);
149 ch_ri.nommer("Richardson");
150 ch_ri.fixer_nb_comp(1);
151 ch_ri.fixer_nb_valeurs_nodales(nb_faces);
152 ch_ri.fixer_unite("-");
153 c_Richardson = 1 ;
154 champs_compris_.ajoute_champ(ch_ri);
155 break ;
156
157 }
158 case 2 :
159 {
160 Cerr << " Lire Pression_porosite " << finl;
161 const Domaine_dis_base& zdis=mon_equation->inconnue().domaine_dis_base();
162 const Domaine_VEF& domaine_VEF=ref_cast(Domaine_VEF, zdis);
163 // const Probleme_base& pb = mon_equation->probleme();
164 const int nb_elem = domaine_VEF.nb_elem() ;
165
166 ch_p.associer_domaine_dis_base(zdis);
167 ch_p.nommer("Pression_porosite");
168 ch_p.fixer_nb_comp(1);
169 ch_p.fixer_nb_valeurs_nodales(nb_elem);
170 ch_p.fixer_unite("Pa.m3/kg");
171 c_pression = 1 ;
172 champs_compris_.ajoute_champ(ch_p);
173 break ;
174 }
175 default :
176 {
177 if (motlu == accfermee)
178 {
179 break ;
180 }
181 else
182 {
183 Cerr << "Erreur dans la lecture de Traitement_particulier_Brech_VEF" << finl;
184 Cerr << "Les mots cles possibles sont : calcul_flux ou Richardson " << finl;
185 Cerr << "Vous avez lu :" << motlu << finl;
186 exit();
187 break;
188 }
189 }
190 }
191 }
192 is >> motlu;
193 if (motlu != accfermee)
194 {
195 Cerr << "Erreur dans la lecture de Traitement_particulier_NS_Brech_VEF 1 ";
196 Cerr << "On attendait une } et pas " << motlu << finl;
197 exit();
198 }
199 // is >> motlu;
200 // if (motlu != accfermee)
201 // {
202 // Cerr << "Erreur dans la lecture de Traitement_particulier_NS_Brech_VEF 2 ";
203 // Cerr << "On attendait une } et pas " << motlu << finl;
204 // exit();
205 // }
206 // is >> motlu;
207 // if (motlu != accfermee)
208 // {
209 // Cerr << "Erreur dans la lecture de Traitement_particulier_NS_Brech_VEF 3 ";
210 // Cerr << "On attendait une } et pas " << motlu << finl;
211 // exit();
212 // }
213 }
214
215 return is;
216}
217
218
219void Traitement_particulier_NS_Brech_VEF::post_traitement_particulier()
220{
221
222 if (c_flux == 1 ) post_traitement_particulier_calcul_flux() ;
223 if (c_Richardson == 1 ) post_traitement_particulier_Richardson() ;
224 if (c_pression == 1 ) post_traitement_particulier_calcul_pression() ;
225}
226
227void Traitement_particulier_NS_Brech_VEF::post_traitement_particulier_calcul_flux()
228{
229
230 for (int ii =0; ii<nb; ii++)
231 {
232 int taille = 1 ;
233 DoubleTab coord_trace (taille,3) ;
234 DoubleTab Surf_trace(taille,3) ;
235 DoubleTab valeurs_(taille,3) ;
236
237 // Modifs VB pour prise en compte rho et calcul du flux enthalpique + Tmoy
238 OBS_PTR(Champ_base) rch1 ;
239 OBS_PTR(Champ_Inc_base) l_inco ;
240 const Probleme_base& pb = mon_equation->probleme();
241 /*
242 pb.a_pour_Champ_Inc("Temperature", rch1 ) ;
243 l_inco = ref_cast(Champ_Inc_base, rch1.valeur()) ;
244 const Champ_Inc_base& temp = l_inco.valeur();
245 */
246 rch1 = pb.get_champ("Temperature");
247 l_inco = ref_cast(Champ_Inc_base, rch1.valeur()) ;
248 const Champ_Inc_base& temp = l_inco.valeur();
249
250 DoubleTab temper_(taille,1);
251 DoubleTab rho_(taille,1) ;
252 DoubleTab cp_(taille,1) ;
253 // Fin Modifs VB
254 //
255 double rad, teta ;
256 double dr = ( r_out(ii) - r_int(ii) ) /delta_r(ii) ;
257 double dteta = 2.*3.14159/delta_teta(ii) ;
258
259 taille = (int)(delta_r(ii)) * (int)(delta_teta(ii)) ;
260 coord_trace.resize(taille,3) ;
261 Surf_trace.resize(taille,3) ;
262 valeurs_.resize(taille,3);
263
264 // Modifs VB pour prise en compte rho et calcul du flux enthalpique + Tmoy
265 rho_.resize(taille,1) ;
266 temper_.resize(taille,1);
267 cp_.resize(taille,1);
268 // Fin Modifs VB
269
270 double a,b,c ;
271
272 // Modifs VB pour calculer le flux sur une face inclinee
273 // anglex est l'angle que fait la normale a la face par rapport a l'axe des x (angle aigu)
274 double anglex = acos(C_trans(ii,0)/sqrt( C_trans(ii,0) * C_trans(ii,0)
275 + C_trans(ii,1) * C_trans(ii,1)
276 + C_trans(ii,2) * C_trans(ii,2) )) ;
277 // Fin Modifs VB
278
279 int i =-1 ;
280
281 for ( rad = (r_int(ii)+(dr/2.)) ; rad < r_out(ii) ; rad += dr )
282 for (teta = dteta/2. ; teta < 2.*3.14159 ; teta += dteta )
283 {
284 i++ ;
285 // Modifs VB pour calculer le flux sur une face inclinee
286 // a = 0. ;
287 // b = rad * cos(teta) ;
288 a = - rad * cos(teta) * sin(anglex) ;
289 b = rad * cos(teta) * cos(anglex) ;
290 // Fin Modifs VB
291 c = rad * sin(teta) ;
292 coord_trace(i,0) = a + R_loc(ii,0) ;
293 coord_trace(i,1) = b + R_loc(ii,1) ;
294 coord_trace(i,2) = c + R_loc(ii,2) ;
295 double ra = (dr/2.+rad) * (dr/2.+rad) ;
296 double ri = (rad - (dr/2.))* (rad - (dr/2.)) ;
297 Surf_trace(i,0) = 3.14159 * (ra -ri) / delta_teta(ii) ;
298 }
299
300 double flux_pos = 0. ;
301 double flux_neg = 0. ;
302 double flux ;
303 // Modifs VB pour prise en compte rho et calcul du flux enthalpique + Tmoy
304 double tempmoy = 0. ;
305 double massvol = 0. ;
306 double fluxE_pos = 0. ;
307 double fluxE_neg = 0. ;
308 double fluxE ;
309 // fin modifs VB
310
311 const Domaine_dis_base& zdis=mon_equation->inconnue().domaine_dis_base();
312 const Domaine& domaine=zdis.domaine();
313 IntVect les_polys(coord_trace.dimension(0));
314 domaine.chercher_elements(coord_trace, les_polys);
315 mon_equation->inconnue().valeur_aux_elems(coord_trace, les_polys, valeurs_);
316
317 // Modifs VB pour prise en compte rho et calcul du flux enthalpique + Tmoy
318 mon_equation->fluide().masse_volumique().valeur_aux_elems(coord_trace, les_polys, rho_);
319 mon_equation->fluide().capacite_calorifique().valeur_aux_elems(coord_trace, les_polys, cp_);
320 temp.valeur_aux_elems(coord_trace, les_polys, temper_);
321 // Fin Modifs VB
322
323 // On annulle la partie virtuelle
324 for (i=0 ; i<taille; i++)
325 if (les_polys(i)>=domaine.nb_elem())
326 {
327 for (int k=0; k<dimension; k++)
328 valeurs_(i,k)=0;
329 // Modifs VB pour prise en compte rho et calcul du flux enthalpique + Tmoy
330 rho_(i)=0.;
331 temper_(i)=0.;
332 cp_(i)=0.;
333 }
334 // Fin Modifs VB
335
336 for (i=0 ; i<taille; i++ )
337 {
338 // Modifs VB pour prise en compte rho et calcul du flux enthalpique + Tmoy
339 // flux = valeurs_(i,0)*Surf_trace(i,0) ;
340 massvol += rho_(i) ;
341 tempmoy += temper_(i) ;
342 flux = ( valeurs_(i,0)*cos(anglex) + valeurs_(i,1)*sin(anglex) )
343 * Surf_trace(i,0)*rho_(i) ;
344 fluxE = ( valeurs_(i,0)*cos(anglex) + valeurs_(i,1)*sin(anglex) )
345 * Surf_trace(i,0)*rho_(i)*cp_(i)*temper_(i) ;
346 flux=Process::mp_sum(flux);
347 fluxE=Process::mp_sum(fluxE);
348 if ( flux > 0. ) flux_pos += flux ;
349 if ( flux < 0. ) flux_neg += flux ;
350 if ( fluxE > 0. ) fluxE_pos += fluxE ;
351 if ( fluxE < 0. ) fluxE_neg += fluxE ;
352 }
353
354 tempmoy=Process::mp_sum(tempmoy);
355 massvol=Process::mp_sum(massvol);
356 tempmoy /= taille ;
357 massvol /= taille ;
358
359 if (je_suis_maitre())
360 {
361 Cout << finl ;
362 Cout << " ----------------------- " << finl ;
363 Cout << " Traitement_particulier flux en x (Plan Oyz) a: "<< R_loc(ii,0) << " "
364 << R_loc(ii,1) << " "
365 << R_loc(ii,2) << finl ;
366 Cout << " angle de la normale du plan par rapport a x = " << anglex/3.1415*180. << " (deg) " << finl ;
367 Cout << " flux_pos = " << flux_pos << " flux_neg = " << flux_neg << " (kg/s) " << finl ;
368 Cout << " fluxEnthalpique_pos = " << fluxE_pos << " fluxEnthalpique_neg = " << fluxE_neg << " (Watts) " << finl ;
369 Cout << " masse volumique moyenne sur le plan de coupe : " << massvol << finl ;
370 Cout << " temperature moyenne sur le plan de coupe : " << tempmoy << finl ;
371 }
372 // Fin Modifs VB
373 }
374}
375
376void Traitement_particulier_NS_Brech_VEF::post_traitement_particulier_Richardson()
377{
378 const Domaine_dis_base& zdis=mon_equation->domaine_dis();
379 const Domaine_VEF& domaine_VEF=ref_cast(Domaine_VEF, zdis);
380 const Domaine_Cl_VEF& domaine_Cl_VEF = ref_cast(Domaine_Cl_VEF,mon_equation->domaine_Cl_dis() );
381
382 OBS_PTR(Champ_base) rch1 ;
383 OBS_PTR(Champ_Inc_base) l_inco ;
384
385 const Probleme_base& pb = mon_equation->probleme();
386 /*
387 pb.a_pour_Champ_Inc("Temperature", rch1 ) ;
388
389 l_inco = ref_cast(Champ_Inc_base, rch1.valeur()) ;
390 */
391 rch1 = pb.get_champ("Temperature");
392 l_inco = ref_cast(Champ_Inc_base, rch1.valeur()) ;
393
394 const Champ_Inc_base& temp = l_inco.valeur();
395 const DoubleTab& temperature = temp.valeurs();
396
397 const DoubleTab& vitesse = mon_equation->inconnue().valeurs();
398 const DoubleVect& gravite = mon_equation->fluide().gravite().valeurs();
399 const DoubleVect& beta = mon_equation->fluide().beta_t().valeurs();
400
401 const int nb_faces = domaine_VEF.nb_faces();
402
403 DoubleVect P(nb_faces) ;
404 DoubleVect G(nb_faces) ;
405
406 DoubleVect beta_i(nb_faces) ;
407 beta_i = beta(0) ;
408
409 calculer_terme_production_K( domaine_VEF, domaine_Cl_VEF,
410 P, vitesse);
411
412 calculer_terme_destruction_K( domaine_VEF, domaine_Cl_VEF,
413 G, temperature,
414 beta_i, gravite) ;
415
416 DoubleVect& richard_loc = ch_ri.valeurs();
417 richard_loc = 0. ;
418
419 for ( int i=0; i<nb_faces; i++ )
420 {
421 if (std::fabs (P[i]) > 1.e-9 )richard_loc [i] = -G[i]/P[i] ;
422 else richard_loc [i] = 0. ;
423 }
424
425 ch_ri.changer_temps(mon_equation->inconnue().temps());
426
427}
428
429void Traitement_particulier_NS_Brech_VEF::post_traitement_particulier_calcul_pression()
430{
431 const Domaine_VEF& zvef=ref_cast(Domaine_VEF, mon_equation->domaine_dis());
432 const DoubleVect& porosite_face = mon_equation->milieu().porosite_face();
433 int i,comp;
434 int nb_face = zvef.nb_faces();
435 Operateur_Div divergence = mon_equation->operateur_divergence();
436 Operateur_Grad gradient = mon_equation->operateur_gradient();
437 SolveurSys solveur_pression_ = mon_equation->solveur_pression();
438
439 DoubleTab& pression=mon_equation->pression().valeurs();
440 DoubleTab& vitesse=mon_equation->vitesse().valeurs();
441 //DoubleTrav gradP(la_vitesse.valeurs());
442 //DoubleTrav gradP(vitesse);
443 //DoubleTrav inc_pre(pression);
444 DoubleTab gradP(vitesse);
445 DoubleTab inc_pre(pression);
446
447 // DoubleVect& inc_pre = ch_p.valeurs();
448 inc_pre = 0. ;
449 //DoubleTrav secmem(pression);
450 DoubleTab secmem(pression);
451
452 gradient.calculer(mon_equation->pression().valeurs(),gradP);
453 //gradient.calculer(la_pression,gradP);
454 // Cerr << "gradP " << gradP << finl;
455 //on veut BM-1Bt(psi*Pression)
456 mon_equation->solv_masse().appliquer(gradP);
457
458 DoubleTab grad_temp(vitesse);
459 for(i=0; i<nb_face; i++)
460 {
461 for (comp=0; comp<dimension; comp++)
462 grad_temp(i,comp) = gradP(i,comp)/porosite_face(i);
463 }
464
465 divergence.calculer(grad_temp, secmem);
466 secmem *= -1; // car div =-B
467 solveur_pression_.resoudre_systeme(mon_equation->matrice_pression().valeur(),secmem, inc_pre);
468
469 DoubleVect& la_pression_porosite = ch_p.valeurs();
470 la_pression_porosite = inc_pre ;
471 Cerr << "la_pression " << mon_equation->pression().valeurs() << finl;
472 Cerr << "ch_p.valeurs() " << ch_p.valeurs() << finl;
473}
474
475
476////////////////////////////////////////////////////////////////////////////
477//
478// Implementation des fonctions de la classe Calcul_Production_K_VEF
479//
480/////////////////////////////////////////////////////////////////////////////
481
482
483void Traitement_particulier_NS_Brech_VEF::
484calculer_terme_production_K(const Domaine_VEF& domaine_VEF,const Domaine_Cl_VEF& zcl_VEF,
485 DoubleVect& P,
486 const DoubleTab& vit) const
487{
488 // P est discretise comme K et Eps i.e au centre des elements
489 //
490 // P(elem) = -(2/3)*k(i)*div_U(i) + nu_t(i) * F(u,v,w)
491 //
492 // 2 2 2
493 // avec F(u,v,w) = 2[(du/dx) + (dv/dy) + (dw/dz) ] +
494 //
495 // 2 2 2
496 // (du/dy+dv/dx) + (du/dz+dw/dx) + (dw/dy+dv/dz)
497 //
498 // Rqs: On se place dans le cadre incompressible donc on neglige
499 // le terme (2/3)*k(i)*div_U(i)
500
501 P= 0;
502
503 // Calcul de F(u,v,w):
504
505 // const Domaine& domaine = domaine_VEF.domaine();
506 int nb_elem = domaine_VEF.nb_elem();
507 // const IntTab& elem_faces = domaine_VEF.elem_faces();
508 const IntTab& face_voisins = domaine_VEF.face_voisins();
509 // const int nb_faces = domaine_VEF.nb_faces();
510 const DoubleVect& volumes = domaine_VEF.volumes();
511 // int premiere_face_entier = domaine_VEF.premiere_face_int();
512 // const IntTab& les_Polys = domaine.les_elems();
513
514 int fac=0, poly1, poly2;
515 int nb_faces_ = domaine_VEF.nb_faces();
516
517 // (du/dx du/dy dv/dx dv/dy ...) pour un poly
518
519 ///////////////////////////////////////////////////////////////////////////////////////////////
520 // <
521 // calcul des gradients; < [ Ujp*np/vol(j) ]
522 // j
523 ////////////////////////////////////////////////////////////////////////////////////////////////
524 int n_bord;
525
526 DoubleTab gradient_elem(nb_elem,dimension,dimension);
527 Champ_P1NC::calcul_gradient(vit,gradient_elem,zcl_VEF);
528
529 // On a les gradient_elem par elements
530 ////////////////////////////////////////////////////////////////////////////////////
531
532 double du_dx;
533 double du_dy;
534 double du_dz;
535 double dv_dx;
536 double dv_dy;
537 double dv_dz;
538 double dw_dx;
539 double dw_dy;
540 double dw_dz;
541
542 /////////////////////////////////////////////////////////////////////////////////////////////////////////
543 // Calcul des du/dx dv/dy et des derivees croisees sur les faces de chaque elements dans le cas 2D
544 /////////////////////////////////////////////////////////////////////////////////////////////////////////
545
546 // Boucle sur les bords pour traiter les faces de bord
547 // en distinguant le cas periodicite
548
549 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
550 {
551 const Cond_lim& la_cl = zcl_VEF.les_conditions_limites(n_bord);
552 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
553 int ndeb = le_bord.num_premiere_face();
554 int nfin = ndeb + le_bord.nb_faces();
555
556 if (sub_type(Periodique,la_cl.valeur()))
557 {
558 for (fac=ndeb; fac<nfin; fac++)
559 {
560 poly1 = face_voisins(fac,0);
561 poly2 = face_voisins(fac,1);
562 double a=volumes(poly1)/(volumes(poly1)+volumes(poly2));
563 double b=volumes(poly2)/(volumes(poly1)+volumes(poly2));
564
565 du_dx =a*gradient_elem(poly1,0,0) + b*gradient_elem(poly2,0,0);
566 du_dy =a*gradient_elem(poly1,0,1) + b*gradient_elem(poly2,0,1);
567 dv_dx =a*gradient_elem(poly1,1,0) + b*gradient_elem(poly2,1,0);
568 dv_dy =a*gradient_elem(poly1,1,1) + b*gradient_elem(poly2,1,1);
569
570 // Determination du terme de production
571
572 P(fac) = ( 2*(du_dx *du_dx + dv_dy *dv_dy) + ((du_dy+dv_dx)*(du_dy+dv_dx)));
573
574 if (dimension == 3)
575 {
576 du_dz=a*gradient_elem(poly1,0,2) + b*gradient_elem(poly2,0,2);
577 dv_dz=a*gradient_elem(poly1,1,2) + b*gradient_elem(poly2,1,2);
578 dw_dx=a*gradient_elem(poly1,2,0) + b*gradient_elem(poly2,2,0);
579 dw_dy=a*gradient_elem(poly1,2,1) + b*gradient_elem(poly2,2,1);
580 dw_dz=a*gradient_elem(poly1,2,2) + b*gradient_elem(poly2,2,2);
581
582 // Determination du terme de production
583
584 P(fac) = (2*( du_dx*du_dx + dv_dy*dv_dy + dw_dz*dw_dz )
585 + ( (du_dy+dv_dx)*(du_dy+dv_dx)
586 + (du_dz+dw_dx)*(du_dz+dw_dx)
587 + (dw_dy+dv_dz)*(dw_dy+dv_dz)));
588 }
589 }
590 }
591 else
592 {
593 for (fac=ndeb; fac<nfin; fac++)
594 {
595 poly1 = face_voisins(fac,0);
596
597 du_dx=gradient_elem(poly1,0,0);
598 du_dy=gradient_elem(poly1,0,1);
599 dv_dx=gradient_elem(poly1,1,0);
600 dv_dy=gradient_elem(poly1,1,1);
601 // Determination du terme de production
602 P(fac) = ( 2*(du_dx*du_dx + dv_dy*dv_dy) + ((du_dy+dv_dx)*(du_dy+dv_dx)));
603
604
605 if (dimension == 3)
606 {
607 du_dz=gradient_elem(poly1,0,2);
608 dv_dz=gradient_elem(poly1,1,2);
609 dw_dx=gradient_elem(poly1,2,0);
610 dw_dy=gradient_elem(poly1,2,1);
611 dw_dz=gradient_elem(poly1,2,2);
612
613 // Determination du terme de production
614
615 P(fac) = (2*( du_dx*du_dx + dv_dy*dv_dy + dw_dz*dw_dz )
616 + ( (du_dy+dv_dx)*(du_dy+dv_dx)
617 + (du_dz+dw_dx)*(du_dz+dw_dx)
618 + (dw_dy+dv_dz)*(dw_dy+dv_dz)));
619 }
620 }
621 }
622 }
623
624
625 // Traitement des faces internes
626
627 for (; fac<nb_faces_; fac++)
628 {
629 poly1 = face_voisins(fac,0);
630 poly2 = face_voisins(fac,1);
631 double a=volumes(poly1)/(volumes(poly1)+volumes(poly2));
632 double b=volumes(poly2)/(volumes(poly1)+volumes(poly2));
633
634 du_dx=a*gradient_elem(poly1,0,0) + b*gradient_elem(poly2,0,0);
635 du_dy=a*gradient_elem(poly1,0,1) + b*gradient_elem(poly2,0,1);
636 dv_dx=a*gradient_elem(poly1,1,0) + b*gradient_elem(poly2,1,0);
637 dv_dy=a*gradient_elem(poly1,1,1) + b*gradient_elem(poly2,1,1);
638
639 // Determination du terme de production
640
641 P(fac) = ( 2*(du_dx*du_dx + dv_dy*dv_dy) + ((du_dy+dv_dx)*(du_dy+dv_dx)));
642
643 if (dimension == 3)
644 {
645 du_dz=a*gradient_elem(poly1,0,2) + b*gradient_elem(poly2,0,2);
646 dv_dz=a*gradient_elem(poly1,1,2) + b*gradient_elem(poly2,1,2);
647 dw_dx=a*gradient_elem(poly1,2,0) + b*gradient_elem(poly2,2,0);
648 dw_dy=a*gradient_elem(poly1,2,1) + b*gradient_elem(poly2,2,1);
649 dw_dz=a*gradient_elem(poly1,2,2) + b*gradient_elem(poly2,2,2);
650
651 // Determination du terme de production
652
653 P(fac) = (2*( du_dx*du_dx + dv_dy*dv_dy + dw_dz*dw_dz )
654 + ( (du_dy+dv_dx)*(du_dy+dv_dx)
655 + (du_dz+dw_dx)*(du_dz+dw_dx)
656 + (dw_dy+dv_dz)*(dw_dy+dv_dz)));
657 }
658 }
659}
660
661
662void Traitement_particulier_NS_Brech_VEF::
663calculer_terme_destruction_K(const Domaine_VEF& domaine_VEF,
664 const Domaine_Cl_VEF& zcl_VEF,
665 DoubleVect& G,const DoubleTab& temper,
666 const DoubleVect& beta,const DoubleVect& gravite) const
667{
668 //
669 // G est discretise comme K et Eps i.e au centre des faces
670 //
671 // --> ----->
672 // G(face) = beta alpha_t(face) G . gradT(face)
673 //
674 int nb_elem = domaine_VEF.nb_elem();
675 int nb_faces= domaine_VEF.nb_faces();
676 DoubleTrav u_teta(nb_faces,dimension);
677 DoubleTrav gradient_elem(nb_elem,dimension);
678 u_teta=0;
679 gradient_elem=0;
680
681 // const Domaine& le_dom=domaine_VEF.domaine();
682 // int nb_faces_elem = le_dom.nb_faces_elem();
683 // const IntTab& elem_faces = domaine_VEF.elem_faces();
684 const IntTab& face_voisins = domaine_VEF.face_voisins();
685 const DoubleTab& face_normales = domaine_VEF.face_normales();
686 int premiere_face_entier = domaine_VEF.premiere_face_int();
687 const DoubleVect& volumes = domaine_VEF.volumes();
688 int nb_faces_ = domaine_VEF.nb_faces();
689 int elem1,elem2,fac;
690
691 //DoubleVect coef(Objet_U::dimension);
692 // const IntTab& les_elem_faces = domaine_VEF.elem_faces();
693
694 // Calcul du gradient de temperature :
695
696 // On traite les bords
697 int n_bord;
698 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
699 {
700 const Cond_lim& la_cl = zcl_VEF.les_conditions_limites(n_bord);
701 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
702 int ndeb = le_bord.num_premiere_face();
703 int nfin = ndeb + le_bord.nb_faces();
704
705 if (sub_type(Periodique,la_cl.valeur()))
706 {
707 for (fac=ndeb; fac<nfin; fac++)
708 {
709 elem1=face_voisins(fac,0);
710 elem2=face_voisins(fac,1);
711 for (int i=0; i<dimension; i++)
712 {
713 gradient_elem(elem1,i) += 0.5*face_normales(fac,i)*temper(fac);
714 gradient_elem(elem2,i) -= 0.5*face_normales(fac,i)*temper(fac);
715 }
716 }
717 }
718 else
719 {
720 for (int fac2=ndeb; fac2<nfin; fac2++)
721 {
722 elem1=face_voisins(fac2,0);
723 for (int i=0; i<dimension; i++)
724 gradient_elem(elem1,i) += face_normales(fac2,i)*temper(fac2);
725 }
726 }
727 }
728
729 // On traite les faces internes
730 for (fac = premiere_face_entier ; fac<nb_faces_; fac++)
731 {
732 elem1=face_voisins(fac,0);
733 elem2=face_voisins(fac,1);
734 for (int i=0; i<dimension; i++)
735 {
736 gradient_elem(elem1,i) += face_normales(fac,i)*temper(fac);
737 gradient_elem(elem2,i) -= face_normales(fac,i)*temper(fac);
738 }
739 }
740
741 for (int elem=0; elem<nb_elem; elem++)
742 for (int i=0; i<dimension; i++)
743 gradient_elem(elem,i) /= volumes(elem);
744
745 // Calcul de u_teta :
746
747 // On traite les bords
748
749 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
750 {
751 const Cond_lim& la_cl = zcl_VEF.les_conditions_limites(n_bord);
752 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
753 int ndeb = le_bord.num_premiere_face();
754 int nfin = ndeb + le_bord.nb_faces();
755
756 if (sub_type(Periodique,la_cl.valeur()))
757 {
758 for (fac=ndeb; fac<nfin; fac++)
759 {
760 elem1 = face_voisins(fac,0);
761 elem2 = face_voisins(fac,1);
762 double a=volumes(elem1)/(volumes(elem1)+volumes(elem2));
763 double b=volumes(elem2)/(volumes(elem1)+volumes(elem2));
764 for (int i=0; i<dimension; i++)
765 u_teta(fac,i)=a*beta(elem1)*gradient_elem(elem1,i)
766 +b*beta(elem2)*gradient_elem(elem2,i);
767 }
768 }
769 else
770 {
771 for (fac=ndeb; fac< nfin; fac++)
772 {
773 elem1 = face_voisins(fac,0);
774 for (int i=0; i<dimension; i++)
775 u_teta(fac,i)=beta(elem1)*gradient_elem(elem1,i);
776 }
777 }
778 }
779
780 for (; fac<nb_faces_; fac++)
781 {
782 elem1 = face_voisins(fac,0);
783 elem2 = face_voisins(fac,1);
784 double a=volumes(elem1)/(volumes(elem1)+volumes(elem2));
785 double b=volumes(elem2)/(volumes(elem1)+volumes(elem2));
786 for (int i=0; i<dimension; i++)
787 u_teta(fac,i)=a*beta(elem1)*gradient_elem(elem1,i)+b*beta(elem2)*gradient_elem(elem2,i);
788 }
789 // -------> ------>
790 // Calcul de gravite . u_teta
791 //
792 G = 0;
793 for (fac=0; fac< nb_faces_; fac++)
794 {
795 if (dimension == 2)
796 G[fac] = gravite(0)*u_teta(fac,0) + gravite(1)*u_teta(fac,1);
797 else if (dimension == 3)
798 G[fac] = gravite(0)*u_teta(fac,0) + gravite(1)*u_teta(fac,1) + gravite(2)*u_teta(fac,2);
799 }
800}
801
Champ_Inc_base
Classe Champ_Inc_base.
Definition Champ_Inc_base.h:53
Champ_Inc_base::valeurs
DoubleTab & valeurs() override
Renvoie le tableau des valeurs du champ au temps courant.
Definition Champ_Inc_base.h:77
Champ_P1NC::calcul_gradient
static DoubleTab & calcul_gradient(const DoubleTab &, DoubleTab &, const Domaine_Cl_VEF &)
Definition Champ_P1NC.cpp:915
Champ_base
classe Champ_base Cette classe est la base de la hierarchie des champs.
Definition Champ_base.h:43
Champ_base::valeur_aux_elems
virtual DoubleTab & valeur_aux_elems(const DoubleTab &positions, const IntVect &les_polys, DoubleTab &valeurs) const
provoque une erreur ! doit etre surchargee par les classes derivees
Definition Champ_base.cpp:245
Domaine_Cl_VEF
Definition Domaine_Cl_VEF.h:40
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::face_normales
virtual double face_normales(int face, int comp) const
Definition Domaine_VF.h:47
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::nb_elem
int nb_elem() const
Definition Domaine_dis_base.h:49
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
Entree
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
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
Motcle
Une chaine de caractere (Nom) en majuscules.
Definition Motcle.h:26
Motcles
Un tableau d'objets de la classe Motcle.
Definition Motcle.h:63
Motcles::search
int search(const Motcle &t) const
Definition Motcle.cpp:321
Objet_U::Entree
friend class Entree
Definition Objet_U.h:76
Objet_U::dimension
static int dimension
Definition Objet_U.h:99
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
Probleme_base
classe Probleme_base C'est un Probleme_U qui n'est pas un couplage.
Definition Probleme_base.h:55
Probleme_base::get_champ
const Champ_base & get_champ(const Motcle &nom) const override
Definition Probleme_base.cpp:841
Process::mp_sum
static double mp_sum(double)
Calcule la somme de x sur tous les processeurs du groupe courant.
Definition Process.cpp:146
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
SolveurSys::resoudre_systeme
int resoudre_systeme(const Matrice_Base &matrice, const DoubleVect &secmem, DoubleVect &solution)
Definition SolveurSys.cpp:48
Sortie
Classe de base des flux de sortie.
Definition Sortie.h:52
Traitement_particulier_NS_Brech_VEF
classe Traitement_particulier_Brech_VEF Cette classe permet de faire les traitements particuliers
Definition Traitement_particulier_NS_Brech_VEF.h:35
Traitement_particulier_NS_Brech_VEF::lire
Entree & lire(Entree &) override
Definition Traitement_particulier_NS_Brech_VEF.cpp:62
Traitement_particulier_NS_Brech_VEF::Traitement_particulier_NS_Brech_VEF
Traitement_particulier_NS_Brech_VEF()
Definition Traitement_particulier_NS_Brech_VEF.cpp:26
Traitement_particulier_NS_Brech_VEF::ch_p
Champ_Fonc_P1NC ch_p
Definition Traitement_particulier_NS_Brech_VEF.h:46
Traitement_particulier_NS_Brech_VEF::post_traitement_particulier
void post_traitement_particulier() override
Definition Traitement_particulier_NS_Brech_VEF.cpp:219
Traitement_particulier_NS_Brech_VEF::ch_ri
Champ_Fonc_P1NC ch_ri
Definition Traitement_particulier_NS_Brech_VEF.h:45
Traitement_particulier_NS_VEF
classe Traitement_particulier_VEF Cette classe permet de ne faire aucyun traitement particulier
Definition Traitement_particulier_NS_VEF.h:31
Traitement_particulier_NS_base::OBS_PTR
OBS_PTR(Navier_Stokes_std) mon_equation
Traitement_particulier_NS_base::champs_compris_
Champs_compris champs_compris_
Definition Traitement_particulier_NS_base.h:68