169 const Domaine_dis_base& zdisbase=mon_equation->inconnue().domaine_dis_base();
171 int nb_faces = domaine_VDF.
nb_faces();
172 const double tps = mon_equation->schema_temps().temps_courant();
175 int nb_elems = domaine_VDF.
nb_elem();
177 int Ny = nb_elems/(nb_faces_y-nb_elems);
186 tau.resize(nb_faces);
187 U_RANS.resize(nb_faces);
188 umoy_spat.resize(nb_faces);
190 umoy_spat.resize(nb_faces);
191 umoy_spat.resize(nb_faces);
193 if (nom_pb_rans ==
"non_couple")
196 EFichier fic_vit(
"vitesse_RANS.dat");
198 for(
int num_face=0 ; num_face<nb_faces ; num_face++)
202 fic_vit >> U_RANS(face) ;
203 fic_verif << face <<
" " << U_RANS(face) << finl;
206 for(
int num_face=0 ; num_face<nb_faces ; num_face++)
208 tau(num_face) = alpha_tau;
222 utemp_gliss.resize(nb_faces);
223 utemp.resize(nb_faces);
227 utemp_gliss.resize(nb_faces);
228 utemp.resize(nb_faces);
229 utemp_sum.resize(nb_faces);
235 umoy.resize(nb_faces);
240 EFichier vit_umoy (
"utemp_sum.dat");
242 SFichier vit_reprise (
"LES.reprise");
243 SFichier vit_reprise2 (
"LES2.reprise");
247 Cerr <<
"trash = " << trash << finl;
249 Cerr <<
"trash = " << trash << finl;
251 for(
int num_face = 0 ; num_face<nb_faces ; num_face++)
253 vit_umoy >> utemp_sum(num_face);
254 vit_umoy2 >> umoy(num_face);
256 vit_reprise << num_face <<
" " << utemp_sum(num_face) << finl;
257 vit_reprise2 << num_face <<
" " << umoy(num_face) << finl;
264 Cerr <<
"Debut initialisation des tableaux de moyenne spatiale..." << finl;
270 const Domaine_dis_base& zdisbase=mon_equation->inconnue().domaine_dis_base();
273 const DoubleTab& xv = domaine_VDF.
xv();
275 int nb_faces = domaine_VDF.
nb_faces();
280 Ny = nb_elems/(nb_faces_y-nb_elems);
285 const IntVect& orientation = domaine_VDF.
orientation();
287 int num_face,j,ori,indic,indicv,indicw,trouve;
313 corresp_u.resize(nb_faces);
314 corresp_v.resize(nb_faces);
315 corresp_w.resize(nb_faces);
333 for (num_face=0; num_face<nb_faces; num_face++)
335 ori = orientation(num_face);
342 for (j=0; j<indic+1; j++)
346 corresp_u[num_face] = j;
355 corresp_u[num_face] = indic;
366 for (j=0; j<indicv+1; j++)
370 corresp_v[num_face] = j;
379 corresp_v[num_face] = indicv;
390 for (j=0; j<indicw+1; j++)
394 corresp_w[num_face] = j;
404 corresp_w[num_face] = indicw;
412 Cerr <<
"Cas de figure impossible!!" << finl;
422 Cerr <<
"Initialisation des tableaux de moyenne spatiale : OK" << finl;
428 const Domaine_VDF& domaine_VDF = le_dom_VDF.valeur();
429 int nb_elems = domaine_VDF.
nb_elem();
430 int nb_faces = domaine_VDF.
nb_faces();
431 const IntTab& face_voisins = domaine_VDF.
face_voisins();
432 const IntTab& elem_faces = domaine_VDF.
elem_faces();
433 DoubleVect norme_vit_elem(nb_elems);
434 norme_vit_elem.
resize(nb_faces);
435 double vit1, vit2, vit3;
436 DoubleTab norme(nb_faces);
438 for(
int num_elem = 0 ; num_elem<nb_elems ; num_elem++)
441 int face1 = elem_faces(num_elem,0);
442 int face2 = elem_faces(num_elem,3);
443 int face3 = elem_faces(num_elem,1);
444 int face4 = elem_faces(num_elem,4);
445 int face5 = elem_faces(num_elem,2);
446 int face6 = elem_faces(num_elem,5);
454 vit1 = 0.5*(U_RANS(face1)+U_RANS(face2));
455 vit2 = 0.5*(U_RANS(face3)+U_RANS(face4));
456 vit3 = 0.5*(U_RANS(face5)+U_RANS(face6));
458 norme_vit_elem(num_elem) = sqrt(vit1*vit1 + vit2*vit2 + vit3*vit3);
463 for(
int num_face = 0 ; num_face<nb_faces ; num_face++)
465 int elem0 = face_voisins(num_face,0);
466 int elem1 = face_voisins(num_face,1);
469 norme(num_face) = norme_vit_elem(elem1);
471 norme(num_face) = norme_vit_elem(elem0);
473 norme(num_face) = 0.5*(norme_vit_elem(elem0)+norme_vit_elem(elem1));
480 DoubleVect& champ_spat_y, DoubleVect& champ_spat_z,
481 DoubleVect& champ_spat_x_2, DoubleVect& champ_spat_y_2,
482 DoubleVect& champ_spat_z_2)
484 const Domaine_VDF& domaine_VDF = le_dom_VDF.valeur();
485 const IntVect& orientation = domaine_VDF.
orientation();
486 int nb_faces = domaine_VDF.
nb_faces();
494 for (
int num_face=0; num_face<nb_faces; num_face++)
496 int ori = orientation[num_face];
500 j=corresp_u[num_face];
501 champ_spat_x[j] += champ(num_face);
502 champ_spat_x_2[j] += champ(num_face)*champ(num_face);
508 j=corresp_v[num_face];
509 champ_spat_y[j] += champ(num_face);
510 champ_spat_y_2[j] += champ(num_face)*champ(num_face);
516 j=corresp_w[num_face];
517 champ_spat_z[j] += champ(num_face);
518 champ_spat_z_2[j] += champ(num_face)*champ(num_face);
525 champ_spat_x[j] = champ_spat_x[j]/compt_x[j];
526 champ_spat_x_2[j] = champ_spat_x_2[j]/compt_x[j];
531 champ_spat_y[j] = champ_spat_y[j]/compt_y[j];
532 champ_spat_y_2[j] = champ_spat_y_2[j]/compt_y[j];
537 champ_spat_z[j] = champ_spat_z[j]/compt_z[j];
538 champ_spat_z_2[j] = champ_spat_z_2[j]/compt_z[j];
613 const Domaine_VDF& domaine_VDF = le_dom_VDF.valeur();
615 const DoubleTab& vitesse = mon_equation->inconnue().valeurs();
616 const double dt = mon_equation->schema_temps().pas_de_temps();
617 const double tps = mon_equation->schema_temps().temps_courant();
618 const double dt_min = mon_equation->schema_temps().pas_temps_min();
620 const IntVect& orientation = domaine_VDF.
orientation();
621 int nb_elems = domaine_VDF.
nb_elem();
622 int nb_faces = domaine_VDF.
nb_faces();
635 SFichier fic_fface(
"force_face.dat", ios::app);
636 const DoubleTab& xv = domaine_VDF.
xv();
643 if(nom_pb_rans !=
"non_couple")
653 U_RANS = pb_rans->equation(0).inconnue().valeurs();
666 int num_face,Nx, Ny, Nz, j;
671 Nx = nb_elems/(nb_faces_x-nb_elems);
672 Ny = nb_elems/(nb_faces_y-nb_elems);
673 Nz = nb_elems/(nb_faces_z-nb_elems);
675 for (num_face=0; num_face<nb_faces; num_face++)
677 int ori = orientation[num_face];
678 vit = vitesse[num_face];
681 j=corresp_u[num_face];
682 umoy_x[j] += vit/((Nx+1)*Nz);
686 j=corresp_v[num_face];
687 umoy_y[j] += vit/(Nx*Nz);
691 j=corresp_w[num_face];
692 umoy_z[j] += vit/(Nx*(Nz+1));
698 for (i=0 ; i<Ny ; i++)
700 fic_moyx << i <<
" " << umoy_x(i) << finl;
702 for (i=0 ; i<Ny+1 ; i++)
704 fic_moyy << i <<
" " << umoy_y(i) << finl;
706 for (i=0 ; i<Ny ; i++)
708 fic_moyz << i <<
" " << umoy_z(i) << finl;
713 for(num_face = 0 ; num_face<nb_faces ; num_face++)
715 int ori = orientation(num_face);
719 j = corresp_u(num_face);
720 umoy(num_face) = umoy_x(j);
724 j = corresp_v(num_face);
725 umoy(num_face) = umoy_y(j);
729 j = corresp_w(num_face);
730 umoy(num_face) = umoy_z(j);
733 Cerr <<
"Erreur orientation dans la redistribution de la moy. spat. de SouCanRANS_LESVDFFa" << finl;
830 for (
int num_face=0; num_face<nb_faces; num_face++)
832 utemp_sum(num_face) = vitesse(num_face)*dt;
833 umoy(num_face) = vitesse(num_face);
837 else if (tps <= t_av)
839 for (
int num_face=0; num_face<nb_faces; num_face++)
841 utemp_sum(num_face) += vitesse(num_face)*dt;
842 umoy(num_face) = utemp_sum(num_face)/(tps-dt_min);
847 for (
int num_face=0; num_face<nb_faces; num_face++)
849 umoy(num_face) = (dt/t_av)*vitesse(num_face)
850 +(1-(dt/t_av))*umoy(num_face);
854 fic_f << tps <<
" " << dt <<
" " << tps-dt_min <<
" "<< utemp_sum(123);
855 fic_f <<
" " << umoy(123) <<
" " << vitesse(123) <<
" " << U_RANS(123) << finl;
899 for (
int num_face=0; num_face<nb_faces; num_face++)
901 utemp_sum(num_face) = vitesse(num_face)*dt;
902 umoy(num_face) = vitesse(num_face);
908 for (
int num_face=0; num_face<nb_faces; num_face++)
910 utemp_sum(num_face) += vitesse(num_face)*dt;
911 umoy(num_face) = utemp_sum(num_face)/(tps-dt_min);
914 fic_f << tps <<
" " << dt <<
" " << tps-dt_min <<
" "<< utemp_sum(123);
915 fic_f <<
" " << umoy(123) <<
" " << vitesse(123) <<
" " << U_RANS(123) << finl;
945 for (
int num_face=0; num_face<nb_faces; num_face++)
953 utemp_sum(num_face) = vitesse(num_face)*dt;
954 umoy(num_face) = vitesse(num_face);
960 for (
int num_face=0; num_face<nb_faces; num_face++)
962 utemp_sum(num_face) += vitesse(num_face)*dt;
963 umoy(num_face) = utemp_sum(num_face)/(tps-dt_min+dt);
970 Cerr <<
"moyenne = " << moyenne << finl;
971 Cerr <<
"Probleme pour le choix du type de moyenne" << finl;
979 DoubleVect force(nb_faces);
981 DoubleVect force_x(nb_faces);
982 DoubleVect force_y(nb_faces);
983 DoubleVect force_z(nb_faces);
985 DoubleVect force_x_2(nb_faces);
986 DoubleVect force_y_2(nb_faces);
987 DoubleVect force_z_2(nb_faces);
989 DoubleVect force_x_temp(nb_faces);
990 DoubleVect force_y_temp(nb_faces);
991 DoubleVect force_z_temp(nb_faces);
993 DoubleVect force_x_temp_2(nb_faces);
994 DoubleVect force_y_temp_2(nb_faces);
995 DoubleVect force_z_temp_2(nb_faces);
1000 for(
int num_face = 0 ; num_face<nb_faces ; num_face++)
1002 vol = volumes_entrelaces(num_face)*porosite_surf(num_face);
1004 force(num_face) = ((U_RANS(num_face)-umoy(num_face))/(dt))
1005 *(U_RANS(num_face)/norm(num_face))*vol;
1009 moy_spat(force, force_x, force_y, force_z, force_x_2, force_y_2, force_z_2);
1019 double i, j, epsilon = 1.e-8;
1020 double dt_post_inst=2.;
1021 double tps_stats=10.;
1023 modf((tps+dt)/dt_post_inst + epsilon, &i);
1024 modf(tps/dt_post_inst + epsilon, &j);
1030 force_x_temp_2, force_y_temp_2, force_z_temp_2, tps-tps_stats);
1038 static int cpt_sonde=0;
1041 for(
int face = 0 ; face<nb_faces ; face++)
1043 if((xv(face,1) > 1.1)&&(xv(face,1) < 0.9)&&(orientation(face)==0))
1046 fic_fface <<
"# " << num_face <<
" "<< orientation(num_face)
1047 <<
" " << xv(num_face,0) <<
" " << xv(num_face,1) <<
" " << xv(num_face,2) << finl;
1048 fic_f <<
"# " << num_face <<
" "<< orientation(num_face)
1049 <<
" " << xv(num_face,0) <<
" " << xv(num_face,1) <<
" " << xv(num_face,2) << finl;
1053 fic_fface << tps <<
" " << vitesse(num_face) <<
" " << force(num_face) << finl;
classe Domaine_Cl_dis_base Les objets Domaine_Cl_dis_base representent les conditions aux limites
int orientation(int) const override
inline DoubleVect& Domaine_VDF::porosite_face() {
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
int face_voisins(int num_face, int i) const
renvoie l'element voisin de numface dans la direction i.
classe Domaine_dis_base Cette classe est la base de la hierarchie des domaines discretisees.