Accueil | Cast3M

* fichier :  smithhutton.dgibi
************************************************************************
************************************************************************
*
*-------------------------------------------------------------------------------
* Convection/Diffusion : CAS SMITH ET HUTTON
* REFERENCE : NUMERICAL HEAT TRANSFER, VOL.5, p.439, 1982
*-------------------------------------------------------------------------------
* Les faces latérales et supérieure d'une boite rectangulaire sont
* imperméables. Le fluide rentre et sort de la boite par la face
* inférieure. Il rentre par la moitié gauche de la face inférieure
* et sort par la moitié droite de cette même face.
*
* Le champ de vitesse est connu et donné par :
*    v(x,y) = 2y(1-x2) i - 2x(1-y2) j
* la taille de la boite étant [-1,1]x[0,1]
*
* L'écoulement transporte un scalaire passif qui est injecté en entrée
* suivant le profil :
*    c(x,0) = 1 + tanh(10(2x+1))
* avec x variant de -1 à 0 (entrée du domaine).
* Sur les frontières imperméables du domaine, la concentration est
* constante et égale à 1-tanh(10).
*
* On cherche la solution stationnaire du transport par diffusion et
* convection du champ scalaire passif. On compare la solution en sortie
* avec la solution de référence. Plusieurs solutions sont calculées
* suivant le Peclet (rapport entre la convection et la diffusion).
*-------------------------------------------------------------------------------
*
'OPTI' 'DIME' 2 'ISOV' 'SULI' ;
*
* Pe : Peclet (convection sur diffusion)
* Les cas traités dans la référence sont 10, 100, 500, 1000 et 1000000
Pe  = 10 ;
*
* KSUPG   : Option de décentrement (CENTREE/SUPG/SUPGDC)
* GRAPH   : Booleen pour l'affichage des tracés à l'issue du calcul
* POST1   : Booleen pour affichage Résidu à chaque pas
* POST2   : Booleen pour affichage C(x,t) à chaque pas
* COMPLET : Booleen modifiant la finesse du maillage et la précision des calculs
* NX      : Nombre de maille suivant x
* NY      : Nombre de maille suivant y
* nbiter  : Nombre de pas de temps
* tolera  : Tolérance pour la comparaison avec la solution de référence
KSUPG   = 'SUPG' ;
GRAPH   = FAUX ;
POST1   = FAUX ;
POST2   = FAUX ;
COMPLET = faux ;
'OPTI' 'ELEM' 'QUA4' ;
'SI' ( COMPLET ) ;
     NY     = 20     ;
     NX     = 2 * NY ;
     nbiter = 2000   ;
     tolera = 0.01   ;
'SINON' ;
     NY     = 10     ;
     NX     = 2 * NY ;
     nbiter = 100    ;
     tolera = 0.2    ;
'FINSI' ;
*
*==================================================================
* Calcul du résidu en température et arrêt suivant un critère
*==================================================================
* E/  : RVX      : TABLE     : TABLE des données créées par EQEX
*                  ARG1      : Fréquence d'impression
*                  ARG2      : Critère d'arrêt
*  /S : MAT1     : MATRIK    : Objet vide
*  /S : CHP1     : CHPO      : Objet vide
*------------------------------------------------------------------
*
'DEBPROC' residu rvx*table ;
   RV   = rvx . 'EQEX' ;
   FREQ = RVX . 'ARG1' ;
   EPS0 = RVX . 'ARG2' ;
   NITER = RV . 'NITER' ;
   DD = RV . 'PASDETPS' . 'NUPASDT' ;
   NN = DD '/' FREQ ;
   L0 = 'EGA' (DD '-' (FREQ*NN)) 0 ;
   'SI' L0 ;
      RANG0 = RV . 'PASDETPS' . 'NUPASDT' ;
      TIME0 = RV . 'PASDETPS' . 'TPS'     ;
      CN0   = RV . 'INCO' . 'CN'  ;
      CNM0  = RV . 'INCO' . 'CN2' ;
      ERR0  = ('MAXIMUM' ('ABS' (CN0 '-' CNM0))) '+' 1.D-20 ;
      ERR10 = ('LOG' ERR0 ) '/' ('LOG' 10.) ;
      'MESSAGE' 'Résidu en Température au pas'
           RANG0 '(' TIME0 ') :' ERR0 ':' ERR10 ;
      RV . 'INCO' . 'IT'  = RV . 'INCO' . 'IT' 'ET' ('PROG' RANG0) ;
      RV . 'INCO' . 'TI'  = RV . 'INCO' . 'TI' 'ET' ('PROG' TIME0) ;
      RV . 'INCO' . 'ER'  = RV . 'INCO' . 'ER' 'ET' ('PROG' ERR10) ; 
      EV1 = 'EVOL' 'MANUEL' (RV . 'INCO' . 'IT') (RV . 'INCO' . 'ER') ;
      Y1 = ('LOG' EPS0) '/' ('LOG' 10) ;
      'SI' POST1 ;
           X1 = 0. ; X2 = RV . 'ITMA' ;
          'DESSIN' EV1 'YBOR' Y1 0. 'NCLK'
                  'TITR' 'Evolution du résidu' ;
      'FINSI' ;
      'SI' POST2 ;
           L1 = (PROG 0. PAS 100. 2000.)* 1.D-3 ;
           trace L1 cn0 DOMTOT CNT1 'TITR' 'Concentration' 'NCLK' ;
      'FINSI' ;
      'SI' ((ERR10 < Y1) 'ET' (DD > 10)) ;
         RV . 'TFINAL' = RV . 'PASDETPS' . 'TPS' ;
      'FINSI' ;
   'FINSI' ;
   RV . 'INCO' . 'CN2' = 'COPIER' RV . 'INCO' . 'CN' ;
   mat1 chp1 = 'KOPS' 'MATRIK' ;
'FINP' mat1 chp1 ;
*------------------------------------------------------------------
*
*
*- MAILLAGE 
*
*
* Points
A1 = -1.0 0.0 ;
A2 =  1.0 0.0 ;
A3 =  1.0 1.0 ;
A4 = -1.0 1.0 ;
A0 =  0.0 0.0 ;
*
* Lignes
LIN  = A1 'DROI' NY A0 ;
LOUT = A0 'DROI' NY A2 ;
FBAS = LIN 'ET' LOUT  ;
FDRO = A2 'DROI' NY A3 ;
FHAU = A3 'DROI' NX A4 ;
FGAU = A4 'DROI' NY A1 ;
LIMP = FDRO 'ET' FHAU 'ET' FGAU ;
*
* Maillage
DOMTOT = 'DALL' FBAS FDRO FHAU FGAU 'PLAN' ;
CNT1   = 'CONT' DOMTOT ;
*
* Modèles et sous-modèles
DOM2     = 'CHAN' 'QUAF' DOMTOT ;
LIN2     = 'CHAN' 'QUAF' LIN ;
LIMP2    = 'CHAN' 'QUAF' LIMP ;
$DOMTOT  = 'MODE' DOM2  'NAVIER_STOKES' 'LINE' ;
$LIN     = 'MODE' LIN2  'NAVIER_STOKES' 'LINE' ;
$LIMP    = 'MODE' LIMP2 'NAVIER_STOKES' 'LINE' ;
*
* Récupération des maillages et fusion des supports
DOMTOT   = 'DOMA' $DOMTOT 'MAILLAGE' ;
LIN      = 'DOMA' $LIN    'MAILLAGE' ;
LIMP     = 'DOMA' $LIMP   'MAILLAGE' ;
FDOMTOT  = 'DOMA' $DOMTOT 'FACE' ;
CLIN     = 'DOMA' $LIN  'CENTRE'  ;
CLIMP    = 'DOMA' $LIMP 'CENTRE'  ;
'ELIM' FDOMTOT (CLIN 'ET' CLIMP) 1.D-3 ; 
*
* Champ de vitesse
XX YY = 'COOR' DOMTOT ;
VXSH  = (2.*YY)*(1.0-(XX*XX)) ;
VYSH  = (-2.*XX)*(1.0-(YY*YY)) ;
VX0   = 'NOMC' 'UX' VXSH ;
VY0   = 'NOMC' 'UY' VYSH ;
VX    = 'KCHT' $DOMTOT 'SCAL' 'SOMMET' 'COMP' 'UX' VX0 ;
VY    = 'KCHT' $DOMTOT 'SCAL' 'SOMMET' 'COMP' 'UY' VY0 ;
CHVIT = 'KCHT' $DOMTOT 'VECT' 'SOMMET' 'COMP' 'UX' 'UY' (VX 'ET' VY) ;
*
* Diffusion (1/Pe)
DIF = 1.0 / ('FLOT' Pe) ;
*
* Profil de concentration à l'entrée
XBAS = 'COOR' 1 LIN ;
TOTO = 2.0*XBAS + 1.0 * 10. ;
SOLUTION = 'TANH' TOTO ;
SOLUTION = 1.0 + SOLUTION ;
CHP1 = 'KCHT' $LIN 'SCAL' 'SOMMET' 0. SOLUTION ;
CHP1 = 'NOMC' 'CN' CHP1 ;
*
* Conditions aux limites en concentration sur les frontières imperméables
C1 = (1.0 - (TANH 10.0)) ;
*
* Description du problème de transport
RV1 = 'EQEX' $DOMTOT 'ITMA' nbiter 'ALFA' 0.7
      'ZONE' $DOMTOT 'OPER' residu 1 1.D-7
      'OPTI' KSUPG 'EFM1' 'EXPL'
      'ZONE' $DOMTOT 'OPER' 'TSCAL' DIF 'VITESSE' 0. 'INCO' 'CN'
      'OPTI' 'CENTREE'
      'ZONE' $DOMTOT 'OPER' 'DFDT' 1. 'CN' 'DELTAT' 'INCO' 'CN' ;
*
* Description des conditions aux limites
RV1 = 'EQEX' RV1
      'CLIM' 'CN' 'TIMP' LIN  CHP1
      'CLIM' 'CN' 'TIMP' LIMP C1 ;
*
* Description des conditions initiales
RV1 . 'INCO' = TABLE 'INCO' ;
RV1 . 'INCO' . 'CN'      = 'KCHT' $DOMTOT 'SCAL' 'SOMMET' 0. ;
*
* Autres data
RV1 . 'INCO' . 'VITESSE' = CHVIT ;
RV1 . 'INCO' . 'CN2'     = 'KCHT' $DOMTOT 'SCAL' 'SOMMET' 0. ;
RV1 . 'INCO' . 'IT'      = 'PROG' ;
RV1 . 'INCO' . 'TI'      = 'PROG' ;
RV1 . 'INCO' . 'ER'      = 'PROG' ;
*
*
*- CALCUL
*
*
EXEC RV1 ;
*
*
*- ANALYSE DES RESULTATS
*
*
EVOL1 = 'EVOL' 'CHPO' (RV1 . 'INCO' . 'CN') 'SCAL' FBAS ;
EVOL2 = 'EVOL' 'CHPO' XX 'SCAL' FBAS ;
LIX   = 'EXTR' EVOL2 'ORDO' ;
LIU   = 'EXTR' EVOL1 'ORDO' ;
EVOL3 = 'EVOL' 'MANU' 'X' LIX 'C(X,0)' LIU ;
*
* Solutions de référence
LIXT = PROG 0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 ;
SI ( Pe 'EGA' 10 ) ;
LIUT = PROG 1.989 1.402 1.146 0.946 0.775 0.621 0.480 0.349
            0.227 0.111 0.000 ;
FINSI ;
SI ( Pe 'EGA' 100 ) ;
LIUT = PROG 2.000 1.940 1.836 1.627 1.288 0.869 0.480 0.209
            0.070 0.017 0.000 ;
FINSI ;
SI ( Pe 'EGA' 500 ) ;
LIUT = PROG 2.000 2.000 1.998 1.965 1.702 0.947 0.242 0.023
            0.001 0.000 0.000 ;
FINSI ;
SI ( Pe 'EGA' 1000 ) ;
LIUT = PROG 2.000 2.000 2.000 1.985 1.841 0.951 0.154 0.001
            0.000 0.000 0.000 ;
FINSI ;
SI ( Pe 'EGA' 1000000 ) ;
LIUT = PROG 2.000 2.000 2.000 1.999 1.964 1.000 0.036 0.001
            0.000 0.000 0.000 ;
FINSI ;
 
EVOL4 = EVOL 'MANU' 'X' LIXT 'Uref(X,0)' LIUT ;
 
LIUC = IPOL LIXT LIX LIU ;
NP = DIME LIXT ;
ERR0 = 0. ;
REPETER BLOC1 NP ;
    UCAL = EXTRAIRE LIUC &BLOC1 ;
    UREF = EXTRAIRE LIUT &BLOC1 ;
    ERR0 = ERR0 + ((UCAL-UREF)*(UCAL-UREF)) ;
FIN BLOC1 ;
ERR0 = ERR0/NP ;
ERR0 = ERR0 '**' 0.5 ;
MESSAGE 'ERREUR ' ERR0 ;
*
*
*- POST-TRAITEMENT
*
*
'SI' GRAPH ;
*
* Maillage
'TRAC' DOMTOT 'TITR' 'Maillage' ;
*
* Vitesse
UNCH = 'VECT' CHVIT 0.1 'UX' 'UY' 'ROUGE' ;
'TRAC' UNCH DOMTOT CNT1 'TITR' 'Vitesse transportante' ;
*
* Concentration
'TRAC' DOMTOT (RV1 . 'INCO' . 'CN') CNT1 'TITR' 'Concentration' ;
*
* Convergence vers la solution stationnaire
EV1 = 'EVOL' 'MANU' 'ITERATIONS' (RV1 . 'INCO' . 'IT')
                    'LOG|E|inf'  (RV1 . 'INCO' . 'ER') ;
'DESS' EV1 'XBOR' 0. ('FLOT' NBITER) 'YBOR' -20.0 0.0
           'TITR' 'Evolution du résidu' ;
*
* Comparaison avec la solution analytique
TAB1 = TABLE ;
TAB1 . TITRE = TABLE ;
TAB1. TITRE . 1 = 'MOT' 'CAST3M' ;
TAB1. TITRE . 2 = 'MOT' 'REFERENCE' ;
TAB1. 2 = 'MARQ LOSA NOLI' ;
'TITR' 'Comparaison CAST3M/Référence' ;
'DESS' (EVOL3 ET EVOL4) 'LEGE' TAB1
       'MIMA' 'XBOR' -1.0 1.0 'YBOR' -1.0 3.0 ;
*
FINSI ;
*
* Arret si problème
*
SI ( ERR0 > tolera ) ;
    'MESS' 'Erreur de ' err0 ' > ' tolera ;
    ERREUR 5 ;
FINSI ; 
'MESS' 'Erreur de ' err0 ' < ' tolera ;
FIN ;