* fichier : chimsour1d.dgibi
GRAPH = faux ;
******************* CAS TEST : chimsour1d.dgibi **********************
* --------------------------------------------------------------------
*
* TRANSPORT GEOCHIMIE AVEC SOURCE ( CAS 1D)
* la partie post-traitement avec graphique contient des exemples
* d'utilisation des procedures TRACHIS TRACHIT DESTRA
* avec des recherches d'identificateurs par NOESPCHI et NOCOMCHI
*
* --------------------------------------------------------------------
* emplacement de COMPOM
emp1 = 'MOT' '/u/castem/divers/COMPOM' ;
*emp1 = 'MOT' 'COMPOM' ;
*emp1 = 'MOT' '/export/home/castem2001/DGIBI/COMPOM' ;
'OPTION' 'ECHO' 0 ;
*
* Génération du maillage
* ----------------------
*
'OPTION' 'DIME' 2 'ELEM' 'QUA4' 'TRACER' 'PSC' ;
* Axe
ORIG = 0.0 0.0 ;
R10 = 1. 0.0;
XAXE = 'DROIT' 1 ORIG R10;
'ELIMINATION' XAXE 0.01;
* Points
Z90 = 0.0 90.0;
Z5 = 0.0 5.0;
Z10 = 0.0 10.0;
Z40 = 0.0 40.0 ;
Z390 = 0.0 390.0;
ZNAD = 0.0 -100.0;
ZM40 = 0.0 -40.0 ;
ZBOT = 0.0 -10.0;
ZB2 = 0.0 -5.0;
ZT2 = 0.0 5.0;
ZTOP = 0.0 10.0;
ZZEN = 0.0 400.0;
*
XAXE2 = 'PLUS' XAXE ZNAD;
*
ZONE1 = 'TRANSLATION' XAXE2 6 ('MOIN' ZM40 ZNAD) ;
XAXM40 = 'INVERSE' ('COTE' 3 ZONE1) ;
YN = 'INVERSE' ('COTE' 4 ZONE1) ;
LIGM40 = 'ELEM' ZONE1 'APPUYE' 'LARGEMENT' XAXM40 ;
ZONE2 = 'TRANSLATION' XAXM40 6 ('MOIN' ZBOT ZM40) ;
XAXBOT = 'INVERSE' ('COTE' 3 ZONE2) ;
YK = 'INVERSE' ('COTE' 4 ZONE2) ;
LIGBOT = 'ELEM' ZONE2 'APPUYE' 'LARGEMENT' XAXBOT ;
ZONE3 = 'TRANSLATION' XAXBOT 8 ('MOIN' ZTOP ZBOT) ;
XAXTOP = 'INVERSE' (COTE 3 ZONE3) ;
LIGTOP = 'ELEM' ZONE3 'APPUYE' 'LARGEMENT' XAXTOP ;
YB = 'INVERSE' (COTE 4 ZONE3) ;
ZONE4 = 'TRANSLATION' XAXTOP 6 ('MOIN' Z40 ZTOP) ;
XAXP40 = 'INVERSE' (COTE 3 ZONE4) ;
YS = 'INVERSE' (COTE 4 ZONE4) ;
ZONE5 = 'TRANSLATION' XAXP40 36 ('MOIN' ZZEN Z40) ;
LIGP40 = 'ELEM' ZONE5 'APPUYE' 'LARGEMENT' XAXP40 ;
XAXE1 = 'INVERSE' (COTE 3 ZONE5) ;
YZ = 'INVERSE' (COTE 4 ZONE5) ;
LIGAX1 = 'ELEM' ZONE5 'APPUYE' 'LARGEMENT' XAXE1 ;
*
YAXE = YN 'ET' YK 'ET' YB 'ET' YS 'ET' YZ;
X80 = 80.0 0.0;
*
MT = ZONE1 'ET' ZONE2 'ET' ZONE3 'ET' ZONE4 'ET' ZONE5 ;
*
'ELIMINATION' (MT 'ET' YAXE) 0.01;
'TASSER' MT;
SDOM = 'ELEM' MT 'CONTENANT' ORIG ;
SDOM = 'COULEUR' SDOM 'ROUGE' ;
AXELEM = 'ELEM' MT 'APPUYE' 'LARGEMENT' YAXE ;
*
* Calcul CHI1
* -----------
*
TABDON = 'TABLE' ;
TABDON.IDEN = 'LECT' 50 5 112 ;
TABDON.CLIM = 'TABLE' ;
TABDON.CLIM . TYP3 = 'LECT' 21442 ;
TABDON.CLIM . COMP3 = 'LECT' 112 ;
TABDON.CLIM . TYP6 = 'LECT' 21440 21441 12721 12722 ;
TB1 = 'CHI1' TABDON 'COMP' emp1 'LOGK' emp1 ;
*
*- Création des maillages hybrides
*
MFTOT = 'CHANER' MT 'QUAF' ;
MFSOUR = 'CHANER' SDOM 'QUAF' ;
MFBAS = 'CHANER' XAXE2 'QUAF' ;
MFAXE = 'CHANER' AXELEM 'QUAF' ;
MFP40 = 'CHANER' LIGP40 'QUAF' ;
MFTOP = 'CHANER' LIGTOP 'QUAF' ;
'ELIMINATION' 0.001
(MFTOT 'ET' MFSOUR 'ET' MFBAS 'ET' MFAXE 'ET' MFP40 'ET' MFTOP) ;
*
*- Modèle
*
MODHYB = 'MODELE' MFTOT 'DARCY' ANISOTROPE ;
MMSOUR = 'MODELE' MFSOUR 'DARCY' ANISOTROPE ;
MMBAS = 'MODELE' MFBAS 'DARCY' ANISOTROPE ;
MMAXE = 'MODELE' MFAXE 'DARCY' ANISOTROPE ;
MMP40 = 'MODELE' MFP40 'DARCY' ANISOTROPE ;
MMTOP = 'MODELE' MFTOP 'DARCY' ANISOTROPE ;
CHYB1 = 'DOMA' MODHYB 'SURFACE' ;
CHYB2 = 'DOMA' MODHYB 'NORMALE' ;
XVOLU = 'DOMA' MODHYB 'VOLUME' ;
CETOT = 'DOMA' MODHYB 'CENTRE' ;
MATOT = 'DOMA' MODHYB 'MAILLAGE ';
FATOT = 'DOMA' MODHYB 'FACE' ;
*
*- On génère la ligne SEGAXE pour le post-traitement
*
CEAXE = 'DOMA' MMAXE 'CENTRE' ;
NBCC = 'NBNO' CEAXE ;
PI1 = 'POINT' CEAXE 1 ;
PI0 = PI1 ;
I = 2 ;
'SI' (NBCC > 1) ;
PI2 = 'POINT' CEAXE I ;
SEGAXE = 'QUELCONQUE' 'SEG2' PI1 PI2 ;
SI (NBCC > 2) ;
NBCC2 = NBCC - 2 ;
'REPETER' BLOC6 nbcc2 ;
PI1 = PI2 ;
I = I + 1 ;
PI2 = 'POINT' CEAXE I ;
LILI = 'QUELCONQUE' 'SEG2' PI1 PI2 ;
SEGAXE = SEGAXE 'ET' LILI ;
'FIN' BLOC6 ;
'FINSI' ;
'FINSI' ;
*
* Données physiques
* -----------------
*
*- on entre la vitesse on en deduit le flux aux faces
*
V = 'MANU' CHPO FATOT 'NATURE' 'DISCRET' 2 'UX' 0. 'UY' 0.1;
VNCH = 'VECTEUR' V 1. 'UX' 'UY' 'ROUGE' ;
MOT1 = 'MOTS' 'UX' 'UY' ;
VAVN = 'PSCAL' V CHYB2 MOT1 MOT1 ;
VAVN = 'NOMC' 'SCAL' VAVN ;
QFACE = VAVN * CHYB1 ;
QFACE = 'NOMC' 'FLUX' QFACE ;
*
*- matériau pour le transport
*
DIR1 = 1. 0. ;
PORO = 0.5 ;
D11C = 1.D0 ;
D22C = 5.D0 ;
D11 = 'MANU' CHML MATOT SCAL D11C ;
D22 = 'MANU' CHML MATOT SCAL (D22C* poro) ;
D21 = 'MANU' CHML MATOT SCAL 1.D-15 ;
MAT1 = 'MATERIAU' MODHYB 'DIRECTION' DIR1 K11 D11 K21 D21 K22 D22 ;
*
*- pas de temps
*
*NBPAST nombre de pas de temps (200 pour un calcul significatif)
DELTAT= 1.25 ;
NBPAST= 10 ;
TETA = 1. ;
TFINAL = NBPAST * DELTAT ;
*
* Initialisation du système
* -------------------------
*
* CTOTC concentrations au centre
* CTOT concentrations aux faces
*
CTOT = 'MANU' 'CHPO' FATOT 3 X050 2.D-9 X005 1.D-9 X112 0.D0 ;
TT1 = 'EXTR' CTOT 'COMP' ;
NOCOMP = 'EXTR' TT1 1 ;
CTOTC = 'MANU' 'CHPO' CETOT 1 NOCOMP 0. ;
NBCOMP = 'DIME' TT1 ;
MASHYB = 'MHYB' MODHYB MAT1 ;
'REPETER' BOUC1 NBCOMP ;
NOCOMP = 'EXTR' TT1 &BOUC1 ;
TPR = 'EXCO' CTOT NOCOMP 'TH' ;
TPC = 'HYBP' MODHYB MASHYB TPR ;
TPC = 'NOMC' NOCOMP TPC ;
CTOTC = CTOTC + TPC ;
'FIN' BOUC1 ;
CLOGC = 'MANU' 'CHPO' CETOT 3 X050 -7. X005 -9. X112 -7. ;
FLOGC = 'MANU' 'CHPO' FATOT 3 X050 -7. X005 -9. X112 -7. ;
*
* Source
*
QELEM0 = 1.D-3;
QELEM = QELEM0 / DELTAT ;
SOURC0 = 'MANU' 'CHPO' CETOT 3 X050 0. X005 0. X112 0. ;
SOURC1 = 'MANU' 'CHPO' ('DOMA' MMSOUR 'CENTRE') 3
X050 (-1.* QELEM) X005 QELEM X112 0. ;
SOURC1 = SOURC1 + SOURC0 ;
SOURC = 'CHARGEMENT' SOURC1
('EVOL' 'MANU' ('PROG' 0. 'PAS' DELTAT (200. * DELTAT))
(PROG 1. 1. 199. * 0.) ) ;
*
* On initialise les concentrations des aqueux aux faces
* par un premier calcul CHI2
*
TBPAR2 = 'TABLE' ;
TBPARM = 'TABLE' ;
TBPARM.'SOUSTYPE' = 'DONNEES_CHIMIQUES' ;
TBPARM.TOT = CTOT ;
TBPARM.LOGC = FLOGC ;
TBPAR2.ITMAX = 25 ;
TBPAR2.EPS = 1.D-4 ;
TBPAR2.NFI = 2 ;
*
TB4 = 'CHI2' TB1 TBPAR2 TBPARM ;
*
TAQU = TB4.AQUE ;
*
* conditions aux limites
*
LIMBAS = 'REDU' TAQU ('DOMA' MMBAS 'CENTRE') ;
BBBAS = 'BLOQ' ('DOMA' MMBAS 'CENTRE') 'TH' ;
*
* Table de données de TRANSPORT
* -----------------------------
*
TABTRAN = 'TABLE' ;
TABTRAN.SOUSTYPE = 'MOT' 'GEOCHIMIE' ;
TABTRAN.'MODELE' = MODHYB ;
TABTRAN.'DIFFUSION' = MAT1 ;
TABTRAN.'POROSITE' = 'MANU' 'CHPO' CETOT 1 'CK' PORO ;
TABTRAN.'CONVECTION' = QFACE ;
TABTRAN.'CHIMI1' = TB1 ;
TABTRAN.'TAQU' = 'TABLE' ;
TABTRAN.'TAQU'. 0 = TAQU ;
TABTRAN.ITMAX = 25;
TABTRAN.EPS = 1.D-4 ;
TABTRAN.NFI = 1 ;
TABTRAN.LOGC = 'TABLE' ;
TABTRAN.LOGC. 0 = CLOGC ;
TABTRAN.TOT = 'TABLE' ;
TABTRAN.TOT. 0 = CTOTC ;
TABTRAN.PRECISION = 5.E-3 ;
TABTRAN.SORTIE = 'MOTS' 'FION' 'SOLU' ;
TABTRAN.'BLOCAGE' = BBBAS ;
TABTRAN.'TRACE_IMPOSE' = 'CHARGEMENT' LIMBAS
('EVOL' 'MANU' ('PROG' 0. 450.) ('PROG' 1. 1.)) ;
TABTRAN.'SOURCE' = SOURC ;
TABTRAN.'PAS_DE_TEMPS' = DELTAT ;
TABTRAN.'TEMPS_FINAL' = TFINAL ;
TABTRAN.'THETA' = 1.D0 ;
*
* Calcul couplé transport/géochimie
* ---------------------------------
*
CHITRNSP TABTRAN ;
*
* Controle des résultats
* ----------------------
*
* Calcul théorique 1D (pour le NA+)
*
VITCENT = 'HVIT' MODHYB QFACE ;
UY = 'EXCO' VY VITCENT SCAL ;
XX YY = 'COOR' CETOT ;
YY = 'NOMC' YY SCAL ;
PETITT = TFINAL ;
DENO1 = (4. * PI * PETITT * D22C) ** 0.5 ;
DENO2 = (4. * PETITT * D22C) ** 0.5 ;
CXYT1 = QELEM0 / PORO / DENO1 ;
CXYT21 = ((YY -((UY / PORO) * PETITT)) / DENO2) ** 2 ;
CXYT2 = CXYT1 * ('EXP' (-1. * CXYT21)) ;
SC1 = CXYT2 * XVOLU ;
SCXYT = 'RESULT' SC1 ;
tttt = 'CHAINE' 'Concentration en NA+ le long de l axe au temps '
petitt ;
'TITRE' TTTT ;
EV100 = 'EVOL' CHPOI CXYT2 'SCAL' SEGAXE ;
*
*- Calcul de l'intégrale de la concentration de NA+ et comparaison
* à la valeur théorique (on compare l'intégrale et la valeur maximale)
*
ECXN = 'EXCO' W002 TABTRAN.'SOLU'. NBPAST 'SCAL' ;
SECX1 = ECXN * XVOLU ;
SECXN = 'RESULT' SECX1 ;
SS1 = 'EXTR' SCXYT 'SCAL' ('POINT' ('EXTR' SCXYT 'MAIL') 1) ;
SS2 = 'EXTR' SECXN 'SCAL' ('POINT' ('EXTR' SECXN 'MAIL') 1) ;
DIFS = ('ABS' (SS1 - SS2) ) / SS1;
'SI' (DIFS < 1.D-2) ;
'ERREUR' 0 ;
'SINON' ;
'ERREUR' 5 ;
'FINSI' ;
MAECXN = 'MAXIMUM' ECXN ;
MACXYT2 = 'MAXIMUM' CXYT2 ;
DIFFNA = 'ABS' (ECXN - CXYT2) ;
MADIFF = 'MAXIMUM' DIFFNA ;
ERRDIF = MADIFF / MACXYT2 ;
'SI' (ERRDIF > 1.D-1) ;
'ERREUR' 5 ;
'FINSI' ;
*
* Post-traitement graphique
* -------------------------
*
'SI' GRAPH ;
EV51 = 'EVOL' 'VERT' CHPOI ECXN 'SCAL' SEGAXE ;
'DESSIN' (ev51 'ET' ev100) ;
MOTIA = 'CHAINE' 'Concentration le long de l axe'
' pour différents temps ';
LLTMP = 'LECT' 0 'PAS' 10 NBPAST ;
LISMD = 'MOTS' 'W002' ;
NE1 NE2 = NOESPCHI TB1 W002 ;
NO1 NO2 NO3 = NOCOMCHI TB1 'NUMCOMP' NE2 ;
TBDES = TRACHIS TABTRAN 'SOLU' LLTMP LISMD ('MOTS' NO1) SEGAXE ;
DESTRA TBDES 'MIMA' 'LOGO' 'DATE' 'TITRE' MOTIA ;
LISMD = 'MOTS' 'X050' ;
NO1 NO2 NO3 = NOCOMCHI TB1 'NOMINT' X050 ;
TBDES = TRACHIS TABTRAN 'AQUE' LLTMP LISMD ('MOTS' NO1) SEGAXE ;
DESTRA TBDES 'MIMA' 'LOGO' 'DATE' 'TITRE' MOTIA ;
LISMD = 'MOTS' 'X005' ;
NO1 NO2 NO3 = NOCOMCHI TB1 'NOMINT' X005 ;
LISMD = 'MOTS' 'X112' ;
NO1 NO2 NO3 = NOCOMCHI TB1 'NOMINT' X112 ;
LISMD = 'MOTS' 'W001' ;
NE1 NE2 = NOESPCHI TB1 W001 ;
NO1 NO2 NO3 = NOCOMCHI TB1 'NUMCOMP' NE2 ;
TBDES = TRACHIS TABTRAN 'SOLU' LLTMP LISMD ('MOTS' NO1) SEGAXE ;
DESTRA TBDES 'MIMA' 'LOGO' 'DATE' 'TITRE' MOTIA ;
LISMD = 'MOTS' 'W002' ;
NE1 NE2 = NOESPCHI TB1 W002 ;
NO1 NO2 NO3 = NOCOMCHI TB1 'NUMCOMP' NE2 ;
TBDES = TRACHIS TABTRAN 'SOLU' LLTMP LISMD ('MOTS' NO1) SEGAXE ;
DESTRA TBDES 'MIMA' 'LOGO' 'DATE' 'TITRE' MOTIA ;
TBN = 'TABLE' ;
NO1 NO2 NO3 = NOCOMCHI TB1 'NUMCOMP' 50 ;
TBN.1 = NO1;
NO1 NO2 NO3 = NOCOMCHI TB1 'NUMCOMP' 5 ;
TBN.2 = NO1;
NO1 NO2 NO3 = NOCOMCHI TB1 'NUMCOMP' 112 ;
TBN.3 = NO1;
MO1 NO1 = NOESPCHI TB1 50 ;
MO2 NO2 = NOESPCHI TB1 5 ;
MO3 NO3 = NOESPCHI TB1 112 ;
LISMD = 'MOTS' MO1 MO2 MO3 ;
CETOP = 'DOMA' MMTOP 'CENTRE' ;
XXTOP YYTOP= 'COOR' CETOP ;
YCTOP = 'EXTR' YYTOP 'SCAL' ('POINT' ('EXTR' YYTOP 'MAIL') 1) ;
TITOP = CHAIN 'Evolution en fonction du temps a '
'FORMAT' '(F6.3)' YCTOP 'metres' ;
TBDES = TRACHIT TABTRAN 'SOLU' LISMD TBN CETOP ;
DESTRA TBDES 'MIMA' 'LOGO' 'DATE' 'TITRE' TITOP ;
LISMD = 'MOTS' 'X050' 'X005' 'X112' ;
TBDES = TRACHIT TABTRAN 'AQUE' LISMD TBN CETOP ;
DESTRA TBDES 'MIMA' 'LOGO' 'DATE' 'TITRE' TITOP ;
'FINSI' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer