Télécharger dvisp.dgibi

Retour à la liste

* fichier :  dvisp.dgibi
************************************************************************
************************************************************************
* CAVITE CARREE ­ VAHL­DAVIS Methode de projection implicite *
* Malvina Renesson Aout 1999 * 
*************************************************************** 
GRAPH=VRAI ;
 GRAPH=FAUX ;
COMPLET=VRAI;
 COMPLET=FAUX ;
 
Si (NON Complet);
itmax= 80 ;
DISCR='LINE' ;
d1 = 0.04; d2 = 0.12; 
DT=.25;
Ra = 1.e6 ;
Sinon ;
*Ra = 1.e5 ;
*itmax= 600 ;
Ra = 1.e6 ;
itmax= 600 ;
DISCR='LINE' ;
d1 = 0.01; d2 = 0.07;
DT=.25;
 
Finsi ;
 
KSUPG ='CENTREE' ;
KPRESS='MSOMMET' ;
 
****************************** 
*ESTIMATION DE LA CONVERGENCE* 
****************************** 
DEBPROC CALCUL; 
ARGU RVX*'TABLE' MDZ/'MMODEL' MDZ/'TABLE';
RV = RVX.'EQEX'; 
DD = RV.PASDETPS.'NUPASDT'; 
NN = DD/5; 
LO = (DD-(5*NN)) EGA 0;
SI (LO); 
UN = RV.INCO.'UN'; 
UNM1 = RV.INCO.'UNM1'; 
unx = kcht $mt scal sommet (exco 'UX' un);
unm1x = kcht $mt scal sommet (exco 'UX' unm1);
uny = kcht $mt scal sommet (exco 'UY' un);
unm1y = kcht $mt scal sommet (exco 'UY' unm1);
ERRX = KOPS unx '-' unm1x;
ERRY = KOPS uny '-' unm1y;
ELIX = MAXI ERRX 'ABS'; 
ELIY = MAXI ERRY 'ABS'; 
ELIX = (LOG (ELIX + 1.0E-10))/(LOG 10.);
ELIY = (LOG (ELIY + 1.0E-10))/(LOG 10.);
MESSAGE 'ITER' RV.PASDETPS.'NUPASDT' 'ERREUR LINF' ELIX ELIY; 
RV.INCO.'UNM1' = KCHT $MT 'VECT' 'SOMMET' (RV.INCO.'UN'); 
IT = PROG RV.PASDETPS.'NUPASDT'; 
ER = PROG ELIY; 
RV.INCO.'IT' = (RV.INCO.'IT') ET IT; 
RV.INCO.'ER' = (RV.INCO.'ER') ET ER; 
FINSI; 
  as2 ama1  = 'KOPS' 'MATRIK' ;
 'FINPROC' as2 ama1 ;
 
************ 
* MAILLAGE * 
************ 
OPTI DIME 2; 
OPTI ELEM QUA8;
crit = 0.0000001; 
xxll = 1.;
p1 = 0. 0.; 
p15 = (xxll/2.) 0.; 
p2 = xxll 0.; 
p25 = xxll 0.5; 
p3 = xxll 1.; 
p35 = (xxll/2.) 1.; 
p4 = 0. 1.; 
p45 = 0. 0.5; 
 
bas1 = p1 d dini d1 dfin d2 p15; 
bas2 = p15 d dini d2 dfin d1 p2;
bas = (bas1 et bas2);elim crit bas; 
cdro = p2 d dini d1 dfin d2 p25 d dini d2 dfin d1 p3; 
haut = p3 d dini d1 dfin d2 p35 d dini d2 dfin d1 p4; 
cgau = p4 d dini d1 dfin d2 p45 d dini d2 dfin d1 p1; 
coupx = bas plus (0. 0.5) ;
coupy=cgau plus (0.5 0.) ;
cnt = bas et cdro et haut et cgau; 
mt = bas cdro haut cgau daller; 
bord = cont mt; 
*mt = surf bord; 
elim (mt et coupx) 1.e-5;
Si GRAPH ; trace mt ; Finsi ;
*option donn 5;
orienter mt; 
******** 
* DOMA * 
******** 
 
Mmt = chan QUAF mt ;
$mt = MODE Mmt 'NAVIER_STOKES' DISCR ;
mt = doma $mt maillage;
 
Mcgau=chan cgau QUAF ; 
$cgau=mode Mcgau 'NAVIER_STOKES' DISCR ;
 
Mcdro=chan cdro QUAF ; 
$cdro=mode Mcdro 'NAVIER_STOKES' DISCR ;
 
Mcnt=chan cnt QUAF ; 
$cnt=mode Mcnt 'NAVIER_STOKES' DISCR ;
 
Mcoupx = chan coupx 'QUAF' ;
$coupx=mode Mcoupx 'NAVIER_STOKES' DISCR ;
 
Mcoupy = chan coupy 'QUAF' ;
$coupy=mode Mcoupy 'NAVIER_STOKES' DISCR ;
 
Elim (Mmt et Mcgau et Mcdro et Mcnt et Mcoupx
et Mcoupy) 1.e-5 ;
cgau=doma $cgau maillage ;
cdro=doma $cdro maillage ;
cnt=doma $cnt maillage ;
coupx=doma $coupx maillage ;
coupy=doma $coupy maillage ;
 
*option donn 5; 
************** 
* PARAMETRES * 
************** 
Pr = 0.71; 
Gr = Ra/Pr; 
NU = 1/(Gr**0.5); 
ALF= NU/Pr; 
gb = kcht $mt vect centre (0. -1.);
uref = 0.2 ;
 
*********************** 
* CREATION DES TABLES *
*********************** 
RV = EQEX  'OMEGA' 1.    'NITER' 1  ITMA itmax   'FIDT' 1
  'OPTI' 'EF' 'IMPL' KSUPG
'ZONE' $MT 'OPER' CALCUL
'ZONE' $MT 'OPER' 'NS' 1. 'UN' NU 'GB' 'TN' 0.5 'INCO' 'UN'
'ZONE' $MT 'OPER' 'TSCAL' 1. 'UN' ALF 0. 'INCO' 'TN'
'OPTI'  'CENTREE'
 ZONE  $mt    OPER  DFDT 1. 'UN' DT  INCO 'UN'
 ZONE  $mt    OPER  DFDT 1. 'TN' DT  INCO 'TN'
; 
****
 
RV = EQEX RV 
'CLIM' 'UN' 'UIMP' cnt 0. 'UN' 'VIMP' cnt 0. 
'TN' 'TIMP' cgau 1.
'TN' 'TIMP' cdro 0.;
 
 rv.'METHINV'.TYPINV=3 ;
 rv.'METHINV'.IMPINV=0 ;
 rv.'METHINV'.NITMAX=400;
 rv.'METHINV'.PRECOND=3 ;
 rv.'METHINV'.RESID  =1.e-10;
 rv. 'METHINV' . 'FCPRECT'=1   ;
 rv. 'METHINV' . 'FCPRECI'=1   ;
 
 
****
* La cavité est fermée il faut imposer la pression en 1 point !
prep1=doma $mt kpress;
bcp=elem  prep1 POI1 (lect 1 ) ;
 
RVP = EQEX
'OPTI' 'EF' KPRESS
ZONE  $mt    OPER  KBBT  -1.  INCO 'UN' 'PRES'
CLIM PRES TIMP bcp 0. ;
 
    rvp.'METHINV'.TYPINV=1 ;
    rvp.'METHINV'.IMPINV=0 ;
    rvp.'METHINV'.NITMAX=300;
    rvp.'METHINV'.PRECOND=3 ;
    rvp.'METHINV'.RESID  =1.e-8 ;
    rvp.'METHINV' . 'FCPRECT'=100 ;
    rvp.'METHINV' . 'FCPRECI'=100 ;
 
 
 
RV.INCO = TABLE INCO;
RV.INCO.'UN' = kcht $MT VECT SOMMET (0. 0.); 
*RV.INCO.'TN' = kcht $MT SCAL SOMMET (1.-corsx);
RV.INCO.'TN' = kcht $MT SCAL SOMMET 0.5 ;
RV.INCO.'PRES' = kcht $MT SCAL KPRESS 0. ;
RV.INCO.'UNM1' = kcht $MT VECT SOMMET (1.E-5 1.E-5);
RV.INCO.'IT' = PROG 1; 
RV.INCO.'ER' = PROG 0.; 
RV.INCO.'GB' =  gb ;
RV.'PROJ' = RVP;
 
EXEC RV; 
 
************* 
* RESULTATS * 
************* 
*option sauv 'SAV';
*sauv Pr Ra Gr NU ALF gb uref RV;
*************
OPTI ISOV SULI; 
 
 
un=(RV.INCO.'UN') * (Pr*(Gr**0.5)) ;
unch = vect un (0.3/(Pr*(Gr**0.5))) UX UY JAUNE; 
 
u1y = exco un 'UX' ;
u2x = exco un 'UY' ;
evu=evol 'CHPO'  u2x  coupx ;
evv=evol 'CHPO'  u1y  coupy ;
 
EVOL4 = EVOL 'MANU' 'ITERATIONS' (RV.INCO.'IT') 'LOG---E---inf' 
(RV.INCO.'ER');
 
GRADT = KOPS RV.INCO.'TN' 'GRAD' $MT; 
DTDX = KCHT $MT 'SCAL' 'CENTRE' (EXCO 'UX' GRADT); 
DTDX = (-1)*(ELNO $MT DTDX); 
 
EVOLd = EVOL 'CHPO' DTDX SCAL cdro ;
EVOLg = EVOL 'CHPO' DTDX SCAL cgau ;
evol1=evold et evolg ;
LISTE EVOL1; 
TAB1 = TABLE; 
TAB1.1 = 'MARQ LOSA';
TAB1.2 = 'MARQ CROI'; 
 
Si GRAPH ; 
DESS EVOL1 'XBOR' 0. 1. 'GRIL' TAB1 'TITRE' ' Nusselt' ; 
trac RV.INCO.'TN' mt (cont mt) 14 'TITRE' ' Isothermes'  ;  
trace unch mt (cont mt) 'TITRE' ' Vitesses' ; 
dess evu  'TITRE' 'Coupe ox à y=1/2 (uy)' ;
dess evv  'TITRE' 'Coupe oy à x=1/2 (ux)' ;
dess EVOL4 'XBOR' 0. itmax 'YBOR' -10.0 0.0;
Finsi;
 
Nusd = (redu DTDX cdro);
Nusg = (redu DTDX cgau);
 
Numg = somt (Nusg*(doma $cgau 'XXDIAGSI')) ;
Numd = somt (Nusd*(doma $cdro 'XXDIAGSI')) ;
 
un2=KOPS un '*' (Pr*(Gr**0.5)); 
sw = KOPS un 'ROT' $mt ; 
rk = EQEX $mt 'OPTI' 'EF' 'IMPL' 
ZONE $mt OPER LAPN -1. INCO 'PSI' 
ZONE $mt OPER FIMP sw INCO 'PSI' 
'CLIM' 'PSI' 'TIMP' cnt 0.; 
rk.'INCO'=table 'INCO' ; 
rk.'INCO'.'PSI'=kcht $mt scal sommet 0. ; 
EXec rk ;
psi=rk.'INCO'.'PSI'; 
psi2=kops psi '*' (Pr*(Gr**0.5)); 
 
Si GRAPH ; 
trac psi mt CNT 14 'TITRE' 'Fonction de courant '; 
Finsi;
 
opti 'ECHO' 0 ;
MESS '----------- Ra = ' Ra '------------------------------------';
MESS '                     NUSSELT MOYEN   NUSSELT MAX   NUSSELT MIN' ;
MESS 'CAST3M : Gauche     ' NUMG '  ' (MAXI Nusg) '  ' (MINI Nusg)    ;
MESS 'CAST3M : Droite     ' NUMD '  ' (MAXI Nusd) '  ' (MINI Nusd)    ;
Si (EGA Ra 1.e5);
MESS 'Vahl Davis : Ra 1.e6    4.519              7.717         0.729';
Finsi ;
Si (EGA Ra 1.e6);
MESS 'Vahl Davis : Ra 1.e6    8.8               17.925         0.989';
Finsi ;
MESS '---------------------------------------------------------------';
MESS '                         Umax            Vmax   ' ;
MESS 'CAST3M      :          '
 (MAXI (exco un 'UX')) (MAXI (exco un 'UY'));
Si (EGA Ra 1.e5);
MESS 'Vahl Davis  : Ra 1.e6    34.73            68.59  ' ;
Finsi ;
Si (EGA Ra 1.e6);
MESS 'Vahl Davis  : Ra 1.e6    64.63           219.36  ' ;
Finsi ;
MESS '---------------------------------------------------------------';
MESS '                          Psi max         ' ;
MESS 'CAST3M      :           '  (MAXI PSI)       ;
Si (EGA Ra 1.e5);
MESS 'Vahl Davis  : Ra 1.e6      9.612          ' ;
Finsi;
Si (EGA Ra 1.e6);
MESS 'Vahl Davis  : Ra 1.e6     16.75           ' ;
Finsi;
MESS '---------------------------------------------------------------';
opti 'ECHO' 1 ; 
 
Si (NON COMPLET) ;
 
ern=abs( NUMD - 7.25)/7.25 ; ernm=abs((MAXI Nusg) - 11.5)/11.5 ;
eru=abs((MAXI (exco un 'UY')) - 277.)/277.;
erp=abs((MAXI PSI) - 25.3)/25.3;
mess 'ern = ' ern ' ernm = ' ernm ' eru = ' eru  ' erp = ' erp ;
Si (ern  > 0.005); erreur 5 ; Finsi ;
Si (ernm > 0.005); erreur 5 ; Finsi ;
Si (eru  > 0.005); erreur 5 ; Finsi ;
Si (erp  > 0.005); erreur 5 ; Finsi ;
 
Finsi ;
 
FIN ; 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales