Download tube_scal_complet.dgibi
* 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 - All rights reserved.
Disclaimer