Télécharger tube_scal_complet.dgibi

Retour à la liste

* fichier :  tube_scal_complet.dgibi
************************************************************************
************************************************************************
******************************************************************
* Transport d'un scalaire dans un tube                           *
*                                                                *
* FORMULATION VF COMPRESSIBLE EXPLICITE/IMPLICITE                *
* SOLVEURS: UPWIND/CENTERED                                      *
*                                                                *
* A. BECCANTINI LTMF    NOVEMBRE 2001                            *
******************************************************************
 
 'OPTION'  'DIME' 2  'ELEM' 'QUA4'  'ISOV' 'SULI' 
   'ECHO' 1 'TRAC' 'X' ;
 
 GRAPH   = FAUX ;
*GRAPH   = VRAI ;
 
******************* 
******************* 
****   MESH   *****
*******************
******************* 
 
 RAF = 16 ;
 NY = 4 ;
 NX = 12 '*' RAF ;
 
 L   = 1.0D0  ;
 DX = L '/' NX '/' 2.0D0;
 H   = NY '*' DX  ;
 
 xD =  0.5D0 '*' L ;
 xG = -1.0D0 '*' xD  ;
 yH =  0.5D0 '*' H   ;
 yB = -1.0D0 '*' yH ;
 
 A1 = xG yB     ;
 A2 = 0.0D0 yB  ;
 A3 = xD yB     ;
 A4 = xD yH     ;
 A5 = 0.0D0 yH  ;
 A6 = xG yH     ;
 VECTG = XG 0.0D0 ;
 VECTD = XD 0.0D0 ;
 
 xBG = xG '-' DX;
 XBD = xD '+' DX;
 
 B1 = xBG yB;
 B2 = xBG yH;
 B3 = xBD yB;
 B4 = xBD yH;
 
 BAS1  = A1 'DROI' NX A2 ;
 BAS2  = A2 'DROI' NX A3 ;
 HAU2  = A4 'DROI' NX A5 ;
 HAU1  = A5 'DROI' NX A6 ;
 LATG  = B1 'DROI' NY B2 ;
 LAT1  = A1 'DROI' NY A6 ;
 LAT12 = A2 'DROI' NY A5 ;
 LAT2  = A3 'DROI' NY A4 ;
 LATD  = B3 'DROI' NY B4 ;
 
 
 DOM1  = LAT12 'TRANSLATION' NX VECTG ;
 DOM2  = LAT12 'TRANSLATION' NX VECTD ;
 
 VECTFG = (-1.0D0 '*' DX) 0.0D0 ;
 VECTFD = DX 0.0D0 ;
 
 FRONTG = LAT1 'TRANSLATION' 1 VECTFG ;
 FRONTG = FRONTG 'COUL' ROUG ;
 
 
*
*** Rotation
*
 
 ANGLE = 30.0D0;
 ORIG = 0.0D0 0.0D0;
 
 DOM1   = DOM1   'TOURNER' ANGLE ORIG;
 DOM2   = DOM2   'TOURNER' ANGLE ORIG;
 FRONTG = FRONTG 'TOURNER' ANGLE ORIG;
 
 DOMINT = DOM1 'ET' DOM2 ;
 'ELIMINATION' DOMINT 1D-6;
 
 DOMTOT = DOMINT 'ET' FRONTG ;
 'ELIMINATION' DOMTOT 1D-6;
 
*
**** EULER model to make the domaine table
*
 
 $DOMINT   = 'MODE' DOMINT 'EULER' ;
 $FRONTG   = 'MODE' FRONTG 'EULER' ;
 $DOMTOT   = 'MODE' DOMTOT  'EULER' ;
 
 TDOMTOT = 'DOMA' $DOMTOT 'VF'  ;
 TDOMINT = 'DOMA' $DOMINT 'VF'  ;
 TFRONTG = 'DOMA' $FRONTG 'VF'  ;
 
 QDOMINT = TDOMINT . 'QUAF' ;
 QFRONTG = TFRONTG . 'QUAF' ;
 QDOMTOT = TDOMTOT . 'QUAF' ;
 'ELIMINATION' QDOMTOT 1D-6 QDOMINT ; 
 'ELIMINATION' QDOMTOT 1D-6 QFRONTG ; 
 
*
******* Post-Treatment line
*
 
 'OPTION' 'ELEM' 'QUA4' ;
 XINIT = XG '+' (0.5D0 '*' DX) ;
 YINIT = YB '+' (0.5D0 '*' DX) ;
 XFIN  = XD '-' (0.5D0 '*' DX) ;
 YFIN = YINIT ;
 PINI = XINIT YINIT;
 PFIN = XFIN YFIN ;
 IAUX = (2 '*' NX) '-' 1 ;
 COURB = PINI 'DROIT' IAUX PFIN;
 COURB = COURB 'TOURNER' ANGLE ORIG;
 COURB = COURB 'COULEUR' 'VERT';
 'ELIMINATION'  0.00001 Courb ('DOMA' $DOMTOT 'CENTRE') ;
 
 'SI' GRAPH ;
    'TRACER' (DOMTOT 'ET' COURB) 'TITRE' ('CHAINE' 'Maillage ');
 'FINSI' ;
 
**********************************
**********************************
****   INITIAL CONDITIONS     ****
**********************************
**********************************
 
*
*** Left state
*
 
 rog   = 1.0  ;
 
*
*** Right state
*
 
 rod   =  1.250D-1 ;
 
* Speed
 
 unscal = 1.0 ;
 
 UX = unscal '*' ('COS' ANGLE) ;
 UY = unscal '*' ('SIN' ANGLE) ;
 
 RO_f  = 'MANUEL' 'CHPO' ('DOMA' $FRONTG 'CENTRE')  1 'SCAL' rog
         'NATU' 'DISCRET';
 
 RO_i  = 'MANUEL' 'CHPO' ('DOMA' $DOMINT 'CENTRE')  1 'SCAL' rod
         'NATU' 'DISCRET';
 
 RN  = RO_f et RO_i ;
 
 UN  = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'FACE') 2 'UX' UX 'UY' UY ;
 
*
*** I.C. graphics
*
 
 MOD1     =  'MODELISER' ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE' ;
 
 'SI' GRAPH ;
    CHM_RN   =  'KCHA' $DOMTOT 'CHAM' RN ;
    'TRACER' CHM_RN  MOD1 'TITR'  ('CHAINE'  'RO at t=' 0.0);
 'FINSI' ;
 
*******************************************
*******************************************
******** EXECUTION ************************
*******************************************
*******************************************
 
 RN0 = 'COPIER' RN ;
 
*
**** Computation parameters
*
 
 TPS = 0.0 ;
 SF = 2. ;
* 
* SF = safety factor
* In a regular mesh SF = 2 <-> CFL = 1
*
 TFINAL  = 0.5 ;
 METO = 'UPWIND' ;
 TEMPDIS = 'EXPLICIT' ;
 
*
* METO = 'CENTERED'
*        'UPWIND'
*
* TEMPDIS = 'EXPLICIT' -> Forward Euler scheme
*           'IMPLICIT' -> Backward Euler scheme
*           'CN'       -> Cranck-Nicholson
*           'SDIRK'    -> Singly Diagonal Implicit Runge Kutta
*
 
*
**** Temporal loop
*
 
 'MESSAGE' ;
 'MESSAGE' ('CHAINE' 'SF      = ' SF) ;
 'MESSAGE' ;
 
*
**** Jacobian and global matrix
*
 
 RF =  'PRET' 'CLAUDEIS' 'FACE' 1  $DOMTOT RN ;
 JACO = 'KONV' 'VF' 'CLAUDEIS' 'FACE' 'JACO' METO
                    $DOMTOT RF UN ;
 RESI DELTAT = 'KONV' 'VF' 'CLAUDEIS' 'FACE' 'RESI' METO
                    $DOMTOT RF UN ;
 
 DT_CON = SF '*' DELTAT ;
 
 'SI' ('EGA' TEMPDIS 'IMPLICIT') ;
    MJACO = 'KOPS' 'MULT' -1. JACO ;
    MATIDE = 'KOPS' 'MATIDE' ('MOTS' 'SCAL')
             ('DOMA' $DOMTOT 'CENTRE') 'MATRIK' ;
    MATASS = ('KOPS' 'MULT' (1. '/' DT_CON) MATIDE) 'ET'
              MJACO ;
 'FINSI' ;
 
 'SI' ('EGA' TEMPDIS 'CN');
    MJACO = 'KOPS' 'MULT' -0.5 JACO ;
    MATIDE = 'KOPS' 'MATIDE' ('MOTS' 'SCAL')
             ('DOMA' $DOMTOT 'CENTRE') 'MATRIK' ;
    MATASS = ('KOPS' 'MULT' (1. '/' DT_CON) MATIDE) 'ET'
              MJACO ;
 'FINSI' ;
 
 'SI' ('EGA' TEMPDIS 'SDIRK') ;
    A11 = (3. '+' (3 '**' 0.5)) '/' 6. ;
    A22 = A11 ;
    A21 = 1. '-' (2. '*' A11) ;
    B1 = 0.5 ;
    B2 = 0.5 ;
    MJACO = 'KOPS' 'MULT' (-1. '*' A11) JACO ;
    MATIDE = 'KOPS' 'MATIDE' ('MOTS' 'SCAL')
             ('DOMA' $DOMTOT 'CENTRE') 'MATRIK' ;
    MATASS = ('KOPS' 'MULT' (1. '/' DT_CON) MATIDE) 'ET'
              MJACO ;
 'FINSI' ;
 
 
 CLIM = 'MANUEL' 'CHPO' ('DOMA' $FRONTG 'CENTRE') 1 'SCAL' 0.0 ;
 
*********************************** 
**** Matrix inversion *************
***********************************
 
 MATTAB = 'TABLE' 'METHINV' ;
 MATTAB .  'TYPINV' = 5 ;
 MATTAB .  'IMPINV' = 0 ;
 
* Methode de numerotation de DDL
 MATTAB . 'TYRENU' = 'SLOA' ;
 MATTAB . 'PCMLAG' = 'APR2' ;
 MATTAB . 'GMRESTRT' = 50 ;
* ILU
 MATTAB . 'PRECOND' = 3 ;
 MATTAB . 'NITMAX' = 1000 ;
 MATTAB . 'RESID' = 1.D-10 ;
 DELTAR = RN '*' 0.0 ;
 
 MATTAB . 'MATASS' = MATASS ;
 MATTAB . 'MAPREC' = MATASS ;
 MATTAB . 'XINIT' = DELTAR ;
 
******************************
**** Temporal loop     *******
******************************
 
 'REPETER' BL1 -1 ;
 
    RF =  'PRET' 'CLAUDEIS' 'FACE' 1  $DOMTOT RN ;
 
    RESI DELTAT = 'KONV' 'VF' 'CLAUDEIS' 'FACE' 'RESI' METO
                  $DOMTOT RF UN ;
 
    DT_CON = SF '*' DELTAT ;
 
    DTTPS = 'MAXIMUM' ('PROG' ((TFINAL '-' TPS) '*' (1.0001))
          1.0E-15) ;
 
    DTMIN = 'MINIMUM' ('PROG' DT_CON DTTPS) ;
 
 
    'SI' ('EGA' TEMPDIS 'IMPLICIT') ;
*
******* Implicit
*
       'SI' (DTMIN < DT_CON) ;
           MATASS = ('KOPS' 'MULT' (1. '/' DTMIN) MATIDE) 'ET'
                 MJACO ;
           MATTAB . 'MATASS' = MATASS ;
           MATTAB . 'MAPREC' = MATASS ;
       'FINSI' ;
*
       DELTAR = 'KRES' MATASS
                'TYPI' (MATTAB)
                'CLIM' CLIM
                'SMBR' RESI
                'IMPR' 0 ;
    'FINSI' ;
*     
    'SI' ('EGA' TEMPDIS 'CN') ;
       'SI' (DTMIN < DT_CON) ;
           MATASS = ('KOPS' 'MULT' (1. '/' DTMIN) MATIDE)
                    'ET' MJACO ;
           MATTAB . 'MATASS' = MATASS ;
           MATTAB . 'MAPREC' = MATASS ;
       'FINSI' ;
*        
       DELTAR  = 'KRES' MATASS
                 'TYPI' (MATTAB)
                 'CLIM' CLIM
                 'SMBR' RESI
                 'IMPR' 0 ;
    'FINSI' ;
*
    'SI' ('EGA' TEMPDIS 'SDIRK') ;
       'SI' (DTMIN < DT_CON) ;
          MATASS = ('KOPS' 'MULT' (1. '/' DTMIN) MATIDE) 'ET'
              MJACO ;
           MATTAB . 'MATASS' = MATASS ;
           MATTAB . 'MAPREC' = MATASS ;
       'FINSI' ;
*        
       DELTAR1  = 'KRES' MATASS
                  'TYPI' (MATTAB)
                  'CLIM' CLIM
                  'SMBR' RESI
                  'IMPR' 0 ;
 
       RF =  'PRET' 'CLAUDEIS' 'FACE' 1  $DOMTOT (RN '+' (A21 '*'
             DELTAR1)) ;
 
       RESI DELTAT = 'KONV' 'VF' 'CLAUDEIS' 'FACE' 'RESI' METO
                    $DOMTOT RF UN ;
 
       DELTAR2  = 'KRES' MATASS
                  'TYPI' (MATTAB)
                  'CLIM' CLIM
                  'SMBR' RESI
                  'IMPR' 0 ;
 
       DELTAR = (B1 '*' DELTAR1) '+' (B2 '*' DELTAR2) ;
 
    'FINSI' ;
*
    'SI' ('EGA' TEMPDIS 'EXPLICIT') ;
       DELTAR = RESI * DTMIN ;
    'FINSI' ;
*
**** Increment of the variables (convection)
*
 
    RN = RN '+' DELTAR ;
 
    TPS  = TPS '+' DTMIN ;
 
    'SI' (((&BL1 '/' 20) '*' 20) 'EGA' &BL1) ;
       'MESSAGE' ('CHAINE' 'ITER = ' &BL1 '  TPS = ' TPS) ;
    'FINSI' ;
 
    'SI' (TPS > TFINAL) ;
      'QUITTER' BL1 ;
    'FINSI' ;
 
 'FIN' BL1 ;
 
*
*** Solutions 
*
 
 'SI' GRAPH ;
    CHM_RN   =  'KCHA' $DOMTOT 'CHAM' RN ;
    'TRACER' CHM_RN  MOD1 'TITR'  ('CHAINE'  'RO at t=' TFINAL);
 'FINSI' ;
 
*
*** Objects evolutions
*
 
 xx yy = 'COORDONNEE' Courb;
 ss =  'KOPS' ('KOPS' ('COS' ANGLE) '*' xx) '+'
             ('KOPS' ('SIN' ANGLE) '*' yy);
 
 evxx = 'EVOL' 'CHPO' ss Courb ;
 lxx = 'EXTRAIRE'  evxx 'ORDO' ;
 
 x0 = 'MINIMUM' lxx;
 x1 = 'MAXIMUM' lxx ;
 
* ro
 
 evro = 'EVOL'  'CHPO' RN Courb ;
 lro  = 'EXTRAIRE'  evro 'ORDO' ;
 evro = 'EVOL' 'MANU' 'x' lxx 'ro' lro;
 tro  = 'CHAINE'  '1D ' METO  ', ' TEMPDIS ' : RO  '
      '  tmps ' TFINAL  ' elem ' 'QUA4' ;
 
 LROAN = 'PROG' ;
 XPOS = (TPS '*' UNSCAL) '-' (L '/' 2.) ;
 
 'REPETER' BL1 ('DIME' lxx) ;
     Xcel = 'EXTRAIRE' lxx &BL1 ;
     'SI' (XCEL < XPOS) ;
        LROAN = LROAN 'ET' ('PROG' rOG) ;
     'SINON' ;
        LROAN = LROAN 'ET' ('PROG' rOD) ;
     'FINSI' ;
 'FIN' BL1 ;
  evroan = 'EVOL' 'MANU' 'x' lxx 'ro' lroan;
 
 'SI' GRAPH ;
    'DESSIN' (EVRO 'ET' EVROAN) 'TITRE' tro ;
 'FINSI' ;
 
 ERRO = 'MAXIMUM' ('EXTRAIRE' (EVRO '-' EVROAN) 'ORDO') 'ABS' ;
 
 'SI' (ERRO > 1.0D-6) ;
    'MESSAGE' 'Transport of passive scalars: problem' ;
    'ERREUR' 5 ;
 'FINSI' ;
 
 'FIN' ;
 
 
 
 
 
 
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales