* fichier : chimsour1d.dgibi ************************************************************************ * Section : Fluides Transport ************************************************************************ GRAPH = faux ; * repertoire des fichiers "divers" DIVERS = VENV 'CASTEM_DIVERS'; * ******************* 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 = 'CHAINE' 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' ;