* fichier : cc2d3.dgibi
**
** --- 10 Novembre 1999 ---
**
** TEST CAVITE CARREE A PAROI DEFILANTE RE=400
**
** teste LAPL KONV (CENTREE) KBBT pression continue
** Algorithme de projection
** Elements P1-P1 Q1-Q1 2D stables si DT >= DT limite
** DT pas trop grand non plus sinon pas de convergence
** Avec l'algorithme VPI2 DT <= 0.3 alors que DT=1 marche avec VPI1
** TRI3 QUA4
** Algorithme VPI2
OPTION ISOV 'SULI' ;
OPTION TRACE X;
GRAPH = 'N' ;
COMPLET = FAUX;
*COMPLET = VRAI;
Rey = 400. ;
ds1=0.01 ; ds2=0.1 ;
Si COMPLET;
ITMAX = 50 ;
err1=1.e-2 ;
Sinon ;
ITMAX = 30 ;
err1=3.e-2 ;
Finsi ;
KPRESS = 'MSOMMET';
DISCR = 'LINE' ;
TYPK = 'QUA4' ;
KSUPG = 'CENTREE';
option dime 2 elem TYPK ;
p1= 0 0 ; p12=0.5 0. ; p2= 1 0 ;
ab=p1 d dini ds1 dfin ds2 p12 d dini ds2 dfin ds1 p2 ;
mt0 = ab trans dini ds1 dfin ds2 (0 0.5)
trans dini ds2 dfin ds1 (0 0.5) ;
mt1= ab trans dini ds1 dfin ds2 (0 0.5);
opti elem TRI3 ;
mt2= (cote 3 mt1) trans dini ds2 dfin ds1 (0 0.5) ;
mt2= inve mt2 ;
bc=cote 2 mt0; cd=cote 3 mt0; da=cote 4 mt0;
AA=bc moin (0.5 0.) ;
BB=ab plus (0. 0.5) ;
elim (mt0 et mt1 et mt2) 1.e-5 ;
mt = mt1 et mt2 ;
elim (AA et BB et mt) 1.e-5 ;
cnt= cont mt ;
mt= chan mt QUAF ;
$mt=mode mt 'NAVIER_STOKES' DISCR ;
MU =1. ;
RO= Rey ;
DT=0.3 ;
* La cavité est fermée il faut imposer la pression en 1 point !
prep1=doma $mt kpress;
bcp=elem prep1 POI1 (lect 1 ) ;
Mmt= doma $mt 'MAILLAGE' ;
doma $mt 'IMPR' ;
CD1= chan CD POI1 ;
CD = elem CD1 (lect 2 pas 1 ((nbel cd1) - 1) );
RV= eqex 'OMEGA' 1. 'NITER' 1 ITMA itmax 'FIDT' 1
'OPTI' 'EF' 'IMPL' KSUPG KPRESS 'DIV2'
ZONE $mt 'OPER' 'LAPN' (MU/RO) 'INCO' 'UN'
ZONE $mt 'OPER' 'KONV' 1. 'UNM' MU DT 'INCO' 'UN'
'OPTI' 'EF' 'BDF2'
ZONE $mt OPER DFDT 1. 'UNM' 'UNMM' DT 'UN' MU INCO 'UN'
;
RV= eqex RV CLIM
UN UIMP CD 1. UN VIMP CD 0. UN UIMP DA 0. UN VIMP DA 0.
UN UIMP AB 0. UN VIMP AB 0. UN UIMP BC 0. UN VIMP BC 0. ;
rv.'METHINV'.TYPINV=3 ;
rv.'METHINV'.IMPINV=0 ;
rv.'METHINV'.NITMAX=400;
rv.'METHINV'.PRECOND=3 ;
rv.'METHINV'.RESID =1.e-8 ;
rv. 'METHINV' . 'FCPRECT'=1 ;
rv. 'METHINV' . 'FCPRECI'=1 ;
RVP= EQEX
'OPTI' 'EF' KPRESS
ZONE $mt OPER KBBT (-1.) INCO 'UN' 'PRES'
CLIM PRES TIMP bcp 0.
;
rvp.'METHINV'.TYPINV=2 ;
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.'PROJ'= RVP ;
RV.'ITMA' = ITMAX ;
rv.inco= table inco ;
rv.inco.'UN' = kcht $mt vect sommet (1.e-5 1.e-5) ;
rv.inco.'UNM' = kcht $mt vect sommet (1.e-5 1.e-5) ;
rv.inco.'UNMM' = kcht $mt vect sommet (1.e-5 1.e-5) ;
rv.inco.'PRES'= kcht $mt scal kpress 0.;
rv.'TYPROJ'='VPI2';
exec rv ;
exec rv ;
AA = chan AA quaf ;
BB = chan BB quaf ;
elim (mt et aa et bb ) 1.e-3 ;
$AA=mode AA 'NAVIER_STOKES' DISCR ;
$BB=mode BB 'NAVIER_STOKES' DISCR ;
srti=doma $AA 'MAILLAGE' ;
srth=doma $BB 'MAILLAGE' ;
evolV = EVOL 'CHPO' (rv.'INCO'.'UN') UX (srti ) ;
evolH = EVOL 'CHPO' (rv.'INCO'.'UN') UY (srth ) ;
evx='EXTR' evolV 'ORDO' ;
evy='EXTR' evolV 'ABSC' ;
evolV= evol 'MANU' 'Vitesse' evx 'Hauteur' evy ;
rv.'EVOLV'=evolV ;
rv.'EVOLH'=evolH ;
si ('EGA' graph 'O' );
TAB1=TABLE;
TAB1.'TITRE'=TABLE ;
TAB1 . 1 ='MARQ REGU ' ;
TAB1.'TITRE' . 1 = mot 'Composante_UX ' ;
DESS evolV 'TITX' 'R (m)' 'TITY' 'V (m/s)' LEGE TAB1 ;
TAB1.'TITRE' . 1 = mot 'Composante_UY ' ;
DESS evolH 'TITX' 'R (m)' 'TITY' 'V (m/s)' LEGE TAB1 ;
c1=vect (rv.inco.'UN') 0.3 ux uy jaune ;
trace c1 mt ;
mtp= doma $mt 'MMAIL' ;
trace (rv.inco.pression) mtp (cont mtp) ;
finsi ;
lrr = prog
.35691E-38 1.83857E-02 3.87990E-02 6.18035E-02 8.78830E-02
.11787 .15314 .19540 .24501 .29606
.32867 .31443 .22963 .12329 9.19502E-03
8.18992E-02 .16068 .22291 .27060 .30865
.34925 .41007 .50336 .62601 .76143
.89040 1.0000;
lr = ABS( 'EXTRAIRE' (rv.'EVOLV') 'ABSC') ;
list lr ;
Dgsrti = doma $AA 'XXDIAGSI' ;
li= 'EXTR' (evol chpo Dgsrti srti) 'ORDO' ;
ER=( resu( ((lr - lrr)*(lr - lrr))*li) )** 0.5 ;
mess ' Ecart Q1-Q1 P1-P1 projection 2D ' er ;
Si ( er > err1 ) ; erreur 5 ; finsi ;
FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer