Download coudep.dgibi

Back to the list

* fichier :  coudep.dgibi
* Ecoulement dans un Coude
* Test 3D pression continue methode de projection
* teste en 3D VNIMP FPU NS
* Ce test pose des problemes en QUAF ??
* Il faut un minimum de mailles pour converger a fort Reynolds
*
opti dime 2 elem QUA8 ;

 GRAPH=VRAI ;
 GRAPH=FAUX ;
 COMPLET = VRAI ;
 COMPLET = FAUX ;

 KPRES='MSOMMET' ;

 Si COMPLET ;
 R=0.5 ; nl = 4 ;n1= 5 ; lgd1=1.; lgd2=2.; n2=5 ;
 NU=1.e-5 ;
 NUT=5.e-3;
 YP=5.e-3 ;
 DT=1.    ;
 Sinon ;
 R=0.5 ; nl = 2 ;n1= 2 ; lgd1=1.; lgd2=2.; n2=2 ;
 NU=1.e-5 ;
 NUT=5.e-2;
 YP=5.e-2 ;
 DT=1.    ;
 Finsi ;

P0=0 0 ; p1= R 0 ;
p2 = p1 tour 90  p0 ;
p3 = p2 tour 90  p0 ;
p4 = p3 tour 90  p0 ;
mct= p1 c p0 n1 p2 c p0 n1 p3 c p0 n1 p4 c p0 n1 p1;
ct1= p1 c p0 n1 p2;
ct2= p2 c p0 n1 p3;
ct3= p3 c p0 n1 p4;
ct4= p4 c p0 n1 p1;

*ct2= p2 d  n1 p3;
*ct4= p4 d  n1 p1;

*mct= p1 c p0 n1 p2 d n1 p0 d n1 p1;
*mc= surf mct ;
mc= daller ct1 ct2 ct3 ct4 ;
*trace mc ;

opti dime 3 elem seg2 ;

q1= p0 plus (0 0 lgd1);
ax1 = q1 plus (1 0 0) ;
ax2 = q1 plus (1 1 0) ;
q2= q1 tour 90. ax1 ax2 ;
q3= q2 plus (lgd2 0 0);

lg= p0 d nl q1  c ax1 n2 q2 d nl q3 ;
*lg= p0 d nl q1  ;

opti elem CU20 ;

*trace lg ;
lg1= chan lg 'POI1' ;
nbp= (nbno lg1) - 2 ;
pi1= p0 ;
si1= mc ;
pi2= poin 2 lg1 ;
vh1=pi2 moins pi1 ;
si2= si1 plus vh1 ;
Entree = si1 ;
vc=si1 volu 2 si2;
*trace cache vc  ;

repeter bloc1 nbp ;
pi1=pi2;
si1=si2 ;
pi2= poin (&bloc1 + 2) lg1 ;
vh2=pi2 moins pi1 ;
vhr=vh2 ;

axe =vh1 pvec vhr ;
alfa= atg (norm(axe / (vh1 psca vhr) ));
Si ( alfa >EG 1.) ;
si2= si1 tour alfa pi1 (pi1 plus axe) ;
Sinon ;
si2= si1 ;
Finsi ;
si2=si2 plus vh2 ;
vh1=vh2 ;
vc = vc et (si1 volu 2 si2);
fin bloc1 ;
Sortie = si2 ;
oeil = 0. -1.e5 -100. ;
*trace cache vc oeil  ;
Sortie= Sortie coul rouge ;
Entree= Entree coul rouge ;

Mvc = chan vc QUAF ;
Mentree = chan entree QUAF ;
Msortie = chan sortie QUAF ;

ELIM (Mvc et Mentree et Msortie) 1.e-5 ;

DEBPROC TEST DISCR*MOT ;

$vc = mode Mvc 'NAVIER_STOKES' DISCR ;
$sortie = mode Msortie 'NAVIER_STOKES' DISCR ;
sortip=doma $sortie KPRES;
$entree = mode Mentree 'NAVIER_STOKES' DISCR ;
entree = doma $entree maillage ;
Menvlp = doma $vc 'ENVELOPPE' ;
Mparoi = diff Menvlp (Mentree et Msortie) ;

doma $vc 'IMPR' ;

* On extrait le champ des normales de la surface d'entree
* en prenant garde a l'orientation des normales
nentr= doma $entree 'NORMALEV';

Mcte  = doma $entree 'ENVELOPPE';
$cte= mode Mcte 'NAVIER_STOKES' DISCR ;
cte=doma $cte maillage ;

* On extrait les composantes nx ny etc
nx= exco 'UX' nentr;
ny= exco 'UY' nentr;
nz= exco 'UZ' nentr;
* on cree les componsantes cartesiennes du champ de vitesse
* qui seront imposees en condition limite type Dirichlet
U0=1. ;
usx=u0*nx ;
usy=u0*ny ;
usz=u0*nz ;

$envlp  = mode Menvlp  'NAVIER_STOKES' DISCR ;
 envlp=doma $envlp maillage ;
$paroi  = mode Mparoi  'NAVIER_STOKES' DISCR ;
paroi=doma $paroi maillage ;


rv= eqex 'ITMA' 5   'NITER' 1  'OMEGA' 1.
  OPTI EF 'SUPG' 'IMPL' KPRES
  'ZONE' $paroi 'OPER' 'FPU' 1. 'UN'  NU 'UET' YP 'INCO' 'UN'
  'ZONE' $vc 'OPER' 'NS' 1. 'UN' NUT 'INCO' 'UN'
  OPTI EF 'CENTREE'
  'ZONE' $vc 'OPER' 'DFDT' 1. 'UN' DT  'INCO' 'UN'
  'CLIM'
  'UN' UIMP entree  usx  'UN' VIMP entree usy 'UN' WIMP entree usz
* 'UN' UIMP paroi 0. 'UN' VIMP paroi 0. 'UN' WIMP paroi 0.
  ;

 rv.'METHINV'.TYPINV=3 ;
 rv.'METHINV'.IMPINV=0 ;
 rv.'METHINV'.NITMAX=100;
 rv.'METHINV'.PRECOND=3 ;
 rv.'METHINV'.RESID  =1.e-10 ;

rv.'INCO'= TABLE 'INCO' ;
rv.'INCO'.'UN'=kcht $vc vect sommet (0. 0. 0.) ;
rv.'INCO'.'PRES'=kcht $vc scal KPRES 0. ;

    RVP= EQEX 'OPTI' 'EF' KPRES
   'ZONE' $vc 'OPER' 'KBBT'  -1.  'INCO' 'UN' 'PRES'
   'ZONE' $paroi 'OPER' 'VNIMP' $vc 0.   'INCO' 'UN' 'PRES'
   'CLIM' 'PRES' TIMP sortip 0.
   ;

 rvp.'METHINV'.TYPINV=3 ;
 rvp.'METHINV'.IMPINV=0 ;
 rvp.'METHINV'.NITMAX=100;
 rvp.'METHINV'.PRECOND=3 ;
 rvp.'METHINV'.RESID  =1.e-10 ;


    RV.'PROJ'= RVP ;

 exec rv ;

 un=rv.inco.'UN' ;
 un1=redu un envlp ;
 ung=vect un1  0.1 ux uy uz jaune ;

Si GRAPH ;
 trace ung envlp oeil ;
Finsi;

 q= dbit un $envlp ;
 list q ;
 qe=dbit un $entree ;
 qs=dbit un $sortie ;
 qp=dbit un $paroi  ;
 mess 'qe=' qe 'qs=' qs 'qp=' qp 'qs+qp=' (qs+qp) ;


 yplus=yp*(rv.inco.'UET')/nu ;
 mess 'Mini yplus ' (mini yplus) ' Maxi yplus ' (maxi yplus) ;

 pn=exco rv.inco.'PRESSION' 'PRES' ;
 mtp= doma $vc 'MMAIL' ;

Si GRAPH ;
 trace pn mtp   ;
Finsi;

 mess 'Mini Pres ' (mini pn) ' Maxi Pres ' (maxi pn) ;

FINPROC ;

*test 'QUAF' ;

test 'MACRO';

test 'LINE' ;

FIN ;
 

 

 

 

© Cast3M 2003 - All rights reserved.
Disclaimer