Accueil | Cast3M

* fichier :  ccar5.dgibi
************************************************************************
************************************************************************
**
**  --- 6 JUIN 1999 ---
**
**   TEST CAVITE CARRE A PAROI DEFILANTE RE=400
**
**   Comparaison des algorithmes Implicite, projection
**   et semi explicite.
**   teste LAPL KONV NS (CENTREE) KBBT
**   Formulation MACRO CENTRE (Stabilisation) ELEMENT QUA9
**   La comparaison fine se fait avec COMPLET=VRAI
**   prevoir espace et temps cpu
 
OPTION ISOV 'SULI' ;
*OPTION TRACE PS;
GRAPH = 'N'                   ;
GRAPHI= 'N'                   ;
COMPLET = FAUX ;
*COMPLET = VRAI ;
Rey = 400.;
 
SI ( COMPLET ) ;
 
ds1=0.02; ds2=0.1 ;
ds1=0.01; ds2=0.08;
ITMAX  = 40 ;
ITMAI  = 30 ;
ITMAE  = 5000  ;
KSUPG  = 'CENTREE';
SINON        ;
 
err1=1.e-4 ;
 ds1=0.2; ds2=0.2 ;
ITMAX  = 30 ;
ITMAI  = 30 ;
ITMAE  = 2000;
ITMAX  = 10;
ITMAI  = 10;
ITMAE  = 100;
KSUPG  = 'SUPG';
 
FINSI ;
 
KPRESS = 'CENTRE' ;
DISCR  = 'LINE'  ;
TYPK   = 'QUA8'   ;
ALGO   = 'SPLT'   ;
 
  DEBPROC TEST KPRESS*MOT TYPK*MOT DISCR*MOT GRAPH*MOT ALGO*MOT
               ITMAX*ENTIER ITMAI*ENTIER ITMAE*ENTIER ;
 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 ;
 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 ;
 AA=bc moin (0.5 0.) ;
 BB=ab plus (0. 0.5) ;
 elim (AA et BB et mt) 1.e-3 ;
 
 mt= chan mt QUAF ;
 $mt=mode mt 'NAVIER_STOKES' DISCR ;
 
MU  =1. ;
RO= Rey  ;
EPSS=-1.e-8 ;
DT=1.    ;
 
* La cavité est fermée il faut imposer la pression en 1 point !
prep1=doma $mt kpress;
bcp=elem  prep1 POI1 (lect 15);
bcp = bcp coul rouge ;
Mmt= doma $mt 'MAILLAGE' ;
doma $mt 'IMPR' ;
Si (EGA GRAPH 'O') ; trace (Mmt et bcp); Finsi ;
*
 
Mab = chan quaf ab ;
Mcd = chan quaf cd ;
Mda = chan quaf da ;
Mbc = chan quaf bc ;
elim (Mmt et Mab et Mcd et Mda et Mbc) 1.e-5 ;
$ab = model Mab 'NAVIER_STOKES' DISCR ;
$cd = model Mcd 'NAVIER_STOKES' DISCR ;
$da = model Mda 'NAVIER_STOKES' DISCR ;
$bc = model Mbc 'NAVIER_STOKES' DISCR ;
ab = doma $ab maillage ;
cd = doma $cd maillage ;
da = doma $da maillage ;
bc = doma $bc maillage ;
 
 Si (EGA ALGO 'EXPL');
 
   RV= eqex  'OMEGA' 1. 'NITER' 1  ITMA itmax   'FIDT' 1
  'OPTI' 'EFM1' KSUPG
  ZONE  $mt    OPER  NS  (MU/RO)        INCO 'UN'
  'OPTI' 'EFM1' 'CENTREE'
  ZONE  $mt    OPER  DFDT 1. 'UN' 'DELTAT'  INCO 'UN' ;
 
 Sinon ;
 
   RV= eqex  'OMEGA' 1.    'NITER' 1  ITMA itmax   'FIDT' 1
  'OPTI' 'EF' 'IMPL' KSUPG
* ZONE  $mt    OPER  NS  (MU/RO)        INCO 'UN'
  ZONE  $mt    'OPER'  'LAPN'  MU  'INCO' 'UN'
  ZONE  $mt    'OPER'  'KONV' RO 'UN' MU DT 'INCO' 'UN'
  'OPTI' 'EFM1' 'CENTREE'
  ZONE  $mt    OPER  DFDT RO  'UN' DT INCO 'UN'
  ;
 
 Finsi ;
 
  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   ;
 
 'SI'  (EGA ALGO 'IMPL') ;
 
   RV.'ITMA' = ITMAI ;
   RV.'OMEGA' = 1.  ;
 
   RV= eqex RV
   'OPTI' 'EF' 'IMPL' KSUPG  KPRESS
   ZONE  $mt    OPER  KBBT 1.            INCO 'UN' 'PRES'
   CLIM PRES TIMP bcp 0. ;
 
 'SINON' ;
 
    Si (EGA ALGO 'EXPL') ; betastab =  1. ;
    Sinon ; betastab=  1.e2 ;
    Finsi ;
    Si (EGA TYPK 'TRI6') ; betastab = 1.e5 ; Finsi ;
 
    RVP= EQEX
    'OPTI' 'EFM1' KPRESS
    ZONE  $mt    OPER  KBBT  -1.   betastab   INCO 'UN' 'PRES'
    CLIM PRES TIMP bcp 0. ;
 
    rvp.'METHINV'.TYPINV=2 ;
    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 ;
 
    si (EGA ALGO 'SPLT') ;
    RV.'PROJ'= RVP ;
    RV.'ITMA' = ITMAX ;
    Finsi ;
 
    si (EGA ALGO 'EXPL') ;
    RV.'PROJ' = RVP ;
    RV.'TYPROJ'='PSCT';
    RV.'ITMA' = ITMAE ;
    RV.'ALFA' = 0.5   ;
    RV.'FIDT' = 50    ;
    Finsi ;
 
 'FINSI' ;
 
rv.inco= table inco ;
rv.inco.un =  kcht $mt  vect sommet (1.e-5 1.e-5) ;
rv.inco.'PRES'= kcht $mt scal kpress 0.;
 
 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' ;
 list evx ;
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 ;
 
Si (EGA ALGO 'IMPL') ;
pn= elno $mt (kcht $mt scal centre ( nomc scal (rv.inco.'PRES')));
sinon ;
pn= elno $mt (kcht $mt scal centre ( nomc scal (rv.inco.pression)));
finsi ;
trace pn mt ;
 
 mess ' MAX P ' (maxi pn) 'MIN P ' (mini pn) ;
finsi ;
FINPROC RV ;
 
 
 
 RV1= test 'CENTRE' TYPK 'MACRO' GRAPH 'SPLT' ITMAX ITMAI ITMAE ;
 
 evv= (rv1.'EVOLV');
 evh= (rv1.'EVOLH');
 
SI ( NON COMPLET ) ;
 
 lr = 'EXTRAIRE' (rv1.'EVOLV') 'ABSC' ;
 list lr ;
 
 lrr = prog
 -8.26684E-38 -2.74863E-02 -5.08159E-02 -7.49178E-02 -9.42076E-02
 -9.99882E-02 -8.10852E-02 -4.32674E-02  7.12863E-03  5.44891E-02
  0.10073      0.30942       1.0000;
 ER=SOMM( abs (lr - lrr) ) ;
 mess ' Ecart sur CENTRE QUA8 MACRO Projection ' er ;
 Si ( er > err1 ) ; erreur 5 ; finsi ;
FINSI ;
 
 
 RV2= test 'CENTRE' TYPK 'MACRO' GRAPH 'IMPL' ITMAX ITMAI ITMAE ;
 
 evv= evv et (rv2.'EVOLV');
 evh= evh et (rv2.'EVOLH');
 
SI ( NON COMPLET ) ;
 
 lr = 'EXTRAIRE' (rv2.'EVOLV') 'ABSC' ;
 list lr ;
 
 lrr=prog
  0.0000      -3.10947E-02 -5.38085E-02 -7.39196E-02 -8.54276E-02
 -8.30593E-02 -6.45006E-02 -3.23991E-02  3.13243E-03  3.90764E-02
  8.72317E-02  0.29477       1.0000;
 
 ER=SOMM( abs (lr - lrr) ) ;
 mess ' Ecart sur CENTRE QUA8 MACRO Implicite ' er ;
 Si ( er > err1 ) ; erreur 5 ; finsi ;
FINSI ;
 
 
 RV3= test 'CENTRE' TYPK 'MACRO' GRAPH 'EXPL' ITMAX ITMAI ITMAE ;
 list RV3.'TYPROJ' ;
 
 evv= evv et (rv3.'EVOLV');
 evh= evh et (rv3.'EVOLH');
 
SI ( NON COMPLET ) ;
 
 lr = 'EXTRAIRE' (rv3.'EVOLV') 'ABSC' ;
 list lr ;
 
 lrr=prog
  3.29048E-37 -1.67015E-02 -2.71125E-02 -3.64044E-02 -4.66237E-02
 -5.61821E-02 -5.96575E-02 -5.28560E-02 -3.42689E-02 -2.85057E-03
  5.39080E-02  0.27875       1.0000;
 
 ER=SOMM( abs (lr - lrr) ) ;
 mess ' Ecart sur CENTRE TRI6 MACRO Projection ' er ;
 Si ( er > err1 ) ; erreur 5 ; finsi ;
FINSI ;
*FIN ;
 
si ('EGA' graphi 'O' );
 
*  Solution de Ghia et al (1982) :
 
y = 'PROG' 0. 0.0547 0.0625 0.0703 0.1016 0.1719 0.2813 0.4531 0.5
           0.6172 0.7344 0.8516 0.9531 0.9609 0.9688 0.9766 1.;
 
'SI' (Rey 'EGA' 100.) ;
   u = 'PROG' 0. -0.03717 -0.04192 -0.04775 -0.06434 -0.10150 -0.15662
             -0.2109 -0.20581 -0.13641 0.00332 0.23151 0.68717 0.73722
              0.78871 0.84123 1.;
'FINSI' ;
 
'SI' (Rey 'EGA' 400.) ;
   u = 'PROG' 0. -0.08186 -0.09266 -0.10338 -0.14612 -0.24299 -0.32726
             -0.17119 -0.11477 0.02135 0.16256 0.29093 0.55892 0.61756
              0.68439 0.75837 1.;
'FINSI' ;
 
'SI' (Rey 'EGA' 1000.) ;
   u = 'PROG' 0. -0.18109 -0.20196 -0.22220 -0.29730 -0.38289 -0.27805
             -0.10648 -0.06080 0.05702 0.18719 0.33304 0.46604 0.51117
              0.57492 0.65928 1.;
'FINSI' ;
 
evolug = 'EVOL' 'MANU' u y;
 
TAB1=TABLE;
TAB1.'TITRE'=TABLE ;
TAB1 . 1 ='MARQ  CROI REGU '     ;
TAB1 . 2 ='MARQ  TRIA REGU '     ;
TAB1 . 3 ='MARQ  LOSA REGU '     ;
TAB1 . 4 ='TIRR MARQ PLUS REGU';
TAB1.'TITRE' . 1 = mot ' Proj '     ;
TAB1.'TITRE' . 2 = mot ' Impl '     ;
TAB1.'TITRE' . 3 = mot ' 1/2 Expl'  ;
TAB1.'TITRE' . 4 = MOT 'V Ghia et al (1982)';
 
evv = evv et evolug ;
 
DESS evv 'TITX' 'R (m)' 'TITY' 'V (m/s)' LEGE  TAB1
DATE MIMA TITR
 ' Cavite carree (Qua9 Macro stab comp ux) ; Comparaison Algorithmes ';
 
*  Solution de Ghia et al (1982) :
 
x = 'PROG' 0. 0.0625 0.0703 0.0781 0.0938 0.1563 0.2266 0.2344 0.5
           0.8047 0.8594 0.9063 0.9453 0.9531 0.9609 0.9688 1.;
 
'SI' (Rey 'EGA' 100.);
   v = 'PROG' 0. 0.09233 0.10091 0.10890 0.12317 0.16077 0.17507 0.17527
              0.05454 -0.24533 -0.22445 -0.16914 -0.10313 -0.08864
             -0.07391 -0.05906 0.;
'FINSI' ;
 
'SI' (Rey 'EGA' 400.);
   v = 'PROG' 0. 0.18360 0.19713 0.20920 0.22965 0.28124 0.30203 0.30174
              0.05186 -0.38598 -0.44993 -0.33827 -0.22847 -0.19254
             -0.15663 -0.12146 0.;
'FINSI' ;
'SI' (Rey 'EGA' 1000.);
   v = 'PROG' 0. 0.27485 0.29012 0.30353 0.32627 0.37095 0.33075 0.32235
              0.02526 -0.31966 -0.42665 -0.51550 -0.39188 -0.33714
             -0.27669 -0.21388 0.;
'FINSI' ;
 
evolvg = 'EVOL' 'MANU' x v;
 
evh = evh et evolvg ;
 
 
DESS evh 'TITX' 'R (m)' 'TITY' 'V (m/s)' LEGE  TAB1
DATE MIMA TITR
 ' Cavite carree (Qua9 Macro stab comp uy) ; Comparaison Algorithmes ';
finsi ;
 
*Si complet ;
*opti sauv 'ccar5.dgibi.sauv' ;
*sauv ;
*finsi ;
FIN ;