* fichier : cc3d1.dgibi
**
** --- 08 Novembre 1999 ---
**
** TEST CAVITE CUBIQUE
**
** teste KBBT NS en 3D + le Bi CG cc2d1.dgibi
** teste les elements a pression continue sur des tetrahedres
** et pyramides P1-P1 et Q1-Q1 algorithme de projection
** stabilises en prenant DT > DT limite
** Teste aussi la creation de volumes automatique
** Si on ne veut que des tetrahedres activer option elem tri3
GRAPH = 'N' ;
err1=5.e-2;
DISCR='LINE' ;
kpress='MSOMMET';
COMPLET = FAUX ;
SI(COMPLET) ;
ds1=0.05 ; ds2=0.08 ;
ITMAX=20 ;
SINON ;
ds1=0.1 ; ds2=0.1 ;
ITMAX=6;
FINSI ;
option dime 2 elem qua4 ;
*option dime 2 elem tri3 ;
opti isov suli ;
p1= 0 0 ; p12=0.5 0. ; p2= 1 0 ;
ab=p1 d dini ds1 dfin ds2 p12 d dini ds2 dfin ds1 p2 ;
ab12= p12 d dini ds1 dfin ds2 (0.5 0.5) d dini ds2
dfin ds1 (0.5 1.) ;
mt= ab trans dini ds1 dfin ds2 (0 0.5) trans dini ds2
dfin ds1 (0 0.5) ;
bc=cote 2 mt ;
cd=cote 3 mt ;
da=cote 4 mt ;
ct=ab et bc et cd et da ;
elim ct 1.e-3 ;
mt=ab bc cd da daller ;
mt=orie mt ;
opti dime 3 elem cub8;
prof=0.1;
mth=mt plus (0 0 prof) ;
ab= ab12 plus (0 0 prof) ;
oeil = 10 10 100 ;
nph=2 ;
cav=mt volu nph mth ;
*trace oeil cav;
f1=face 1 cav ;
f2=face 2 cav ;
f3=face 3 cav ;
psup=cd trans nph (0 0 prof) ;
psup=(CHAN TRI3 psup) ;
ff3= (CHAN TRI3 f3) ;
ELIM (PSUP ET FF3) 1.e-4 ;
env3= f1 et f2 et ff3 ;
cav = volu env3 ;
cav= chan cav quaf ;
ab=chan ab quaf ;
elim (cav et psup et ab et ff3 ) 1.e-5 ;
$mt=mode cav 'NAVIER_STOKES' DISCR ;
$AB=mode AB 'NAVIER_STOKES' DISCR ;
DOMA $mt IMPR ;
DT=2. ;
MU=1. ;
RO= 400. ;
prep1=doma $mt kpress;
bcp=elem prep1 POI1 (lect 1 ) ;
rv= eqex 'OMEGA' 1. 'NITER' 1 ITMA ITMAX
'OPTI' 'EF' 'IMPL' KPRESS 'SUPG'
ZONE $mt OPER NS 1. 'UN' (MU/RO) INCO 'UN'
'OPTI' 'EF' 'BDF2'
ZONE $mt OPER DFDT 1. 'UNM' 'UNMM' DT 'UN' MU INCO 'UN'
;
rv=eqex rv
CLIM
UN UIMP ff3 0. UN VIMP ff3 0. UN WIMP ff3 0.
UN WIMP (f1 et f2) 0.
UN UIMP PSUP 1.
;
rv.inco= table inco ;
rv.inco.tn = kcht $mt scal sommet 5. ;
rv.inco.un = kcht $mt vect sommet (1.e-5 1.e-5 1.e-5) ;
rv.inco.'UNM' = kcht $mt vect sommet (1.e-5 1.e-5 1.e-5) ;
rv.inco.'UNMM' = kcht $mt vect sommet (1.e-5 1.e-5 1.e-5) ;
rv.inco.pres = kcht $mt scal KPRESS 1.e-5 ;
rv.'METHINV'.TYPINV=3 ;
rv.'METHINV'.IMPINV=1 ;
rv.'METHINV'.NITMAX=1000;
rv.'METHINV'.PRECOND=3 ;
rv.'METHINV'.RESID =1.e-10 ;
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=1 ;
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 ;
exec rv ;
srti=doma $AB 'MAILLAGE' ;
evolV = EVOL 'CHPO' (rv.'INCO'.'UN') UX (srti ) ;
evx= extr evolV 'ORDO' ;
list evx ;
evy=extr evolV 'ABSC' ;
evolV= evol 'MANU' 'Vitesse' evx 'Hauteur' evy ;
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 ;
c1=vect (rv.inco.'UN') 0.3 ux uy uz jaune ;
*trace c1 mt ;
trace c1 cav ;
trace (rv.inco.'PRESSION') mt ;
finsi ;
lrr=prog
3.76000E-38 2.84554E-02 4.47889E-02 6.01943E-02 7.05196E-02
6.69563E-02 4.10891E-02 4.04425E-03 5.94186E-02 0.26219
1.0000 ;
evx = ABS evx ;
list evx ;
ER=resu( abs (evx - lrr) ) *0.1;
mess ' Ecart sur MSOMMET TET4 PYR5 ' er ;
Si ( er > 5.e-3) ; erreur 5 ; finsi ;
FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer