* fichier : dvispqt.dgibi
* CAVITE CARREE VAHLDAVIS 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;
d1 = 0.05; d2 = 0.2;
DT=.5;
Ra = 1.e5 ;
*Ra = 1.e6 ;
Sinon ;
Ra = 1.e5 ;
*Ra = 1.e6 ;
itmax= 600 ;
DISCR='LINE' ;
d1 = 0.01; d2 = 0.07;
DT=.5;
Finsi ;
KSUPG ='CENTREE' ;
KPRESS='MSOMMET' ;
******************************
*ESTIMATION DE LA CONVERGENCE*
******************************
DEBPROC CALCUL;
ARGU RVX*'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;
mt = surf cnt ;
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=5 ;
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.);
corsx = coor 1 (doma $MT SOMMET) ;
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' ;
**pas defini avec ce maillage: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.e5 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.e5 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.e5 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 - 3.7)/3.7 ; ernm=abs((MAXI Nusg) - 6.4)/6.4 ;
eru=abs((MAXI (exco un 'UY')) - 79.)/79.;
erp=abs((MAXI PSI) - 10.)/10.;
mess 'ern = ' ern ' ernm = ' ernm ' eru = ' eru ' erp = ' erp ;
Si (ern > 0.05 ); erreur 5 ; Finsi ;
Si (ernm > 0.1 ); erreur 5 ; Finsi ;
Si (eru > 0.1 ); erreur 5 ; Finsi ;
Si (erp > 0.05 ); erreur 5 ; Finsi ;
Finsi ;
FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer