$$$$ TRANSGEN * TRANSGEN PROCEDUR JC220346 12/09/12 21:15:09 7501 *TRANSGEN PROCEDUR GBM 04/01/03 'DEBP' TRANSGEN TRANSI*'TABLE'; * SUPPRIMER TOUS LES 'H', stocker composante de CHRG 'ET' * REMPLACER LES 'H' par cela. A faire. *--------------------------------------------------------------------- * Résolution de l'équation de darcy pour un problème d'écoulement ou * de transport par une méthode d'éléments * finis mixtes hybrides ou VF. Les inconnues du problème sont * - en EFMH, la concentration (H), la trace de * concentration (TH) et le débit diffusif (FLUX). * - en VF, la concentration (H) *--------------------------------------------------------------------- * *------------------------------ * Phrase d'appel (en GIBIANE) : *------------------------------ * * TRANSGEN TABLE ; * *---------------------------------- * Opérandes (à mettre dans TABLE) : *---------------------------------- * * ___________________________________________________________________ * | | * | Indice Contenu | * | | * ------------------------------------------------------------------- * | | * |------------------------------------------------ | * |Données physiques, géométriques et materielles : | * |------------------------------------------------ | * | | * |'MODELE' Objet modèle (MMODEL créé par MODE) | * | | * |'CARACTERISTIQUES' Données physiques et materielles : | * | diffusivité effective (CHAMELEM créé par MATE) | * | | * |'POROSITE' Valeur de la porosité (Type Champoint, Comp | * | 'SCAL', ou FLOTTANT) - Défaut 1. | * | | * |'DECROISSANCE' Valeur du terme de décroissance (Type FLOTTANT) | * | Tel que dC/dt = - Lambda * C - Défaut 0. | * | | * |'COEF_RETARD' Coefficient de retard linéaire dans le cas simple, | * | ou Pente à l'origine de la fonction F(C) dans le | * | cas d'isotherme non linéaire de Langmuir | * | ou Coefficient K de l'isotherme de Freundlich | * | (Type CHPO Centre 'SCAL', ou FLOTTANT) | * | | * |'LANGMUIR' Quantité maximale 'Fsat' adsorbée sur le solide | * | rapportée à l'unité de volume du fluide et exprimée| * | dans la meme unité que le soluté. | * | (Type CHPO Centre 'SCAL', ou FLOTTANT). | * | F = (R-1) C / [1 + ((R-1) C / Fsat)] | * | Si cet indice et le suivant sont absents, | * | l'équilibre d'adsorption est linéaire. Cet indice a| * | priorité sur l'indice FREUNDLICH. | * | | * |'FREUNDLICH' Exposant de la loi de Freundlich F = K (C ^ 1/n) | * | (Type FLOTTANT). | * | Dans ce cas (et si l'indice LANGMUIR n'existe pas),| * | l'indice 'COEF_RETARD' contient le coefficient | * | K ramené à une unité de volume de fluide. | * | - Non disponible pour l'instant - | * | | * |'LIMITE_SOLUBILITE' Limite de solubilité (Type chpoin, Comp 'H') | * | | * |'COEF_DISSOLUTION' Coef. de dissolution (Type CHPO Centre, Comp | * | 'SCAL'). Tel que dC/dt = Coef * (Csat - C) - Par | * | défaut, la dissolution est instantanée | * | | * |'CONVECTION' Débit intégré de la vitesse convective à travers | * | chaque face des éléments (Type CHPO Face, comp. | * | 'FLUX') | * | | * |---------------------- | * |Conditions initiales : | * |---------------------- | * | | * |'TEMPS' TABLE contenant à l'indice 0 la valeur du temps | * | initial (FLOTTANT) | * | | * | | * |'CONCENTRATION' TABLE contenant à l'indice 0 la concentration | * | (quantité d'élément par unité de volume d'eau) | * | (Type CHPO Centre, Comp 'H') | * | | * |'FLUX' TABLE contenant à l'indice 0 le flux diffusif | * | initial intégré sur chaque face (Type CHPO Face, | * | comp. 'FLUX') | * | | * |'PRECIPITE' TABLE contenant à l'indice 0 la quantité initiale | * | de précipité par unité de volume de milieu solide | * | (Type CHPO Centre, Comp 'H') | * | | * |'DISSOLUTION' TABLE contenant à l'indice 0 la quantité initiale | * | pour estimer la dissolution au premier pas de temps| * | (Type CHPO, Comp 'H'), voir plus loin. | * | | * |-------------------------------------- | * |Conditions aux limites / chargements : | * |-------------------------------------- | * | | * |'BLOCAGE' Contient les matrices de blocage (RIGIDITE) | * | | * |'TRACE_IMPOSE' Valeurs des traces imposées (charge ou concentra- | * | -tion) (CHARGEMENT 'TH' - Obligatoire si BLOCAGE) | * | | * |'FLUX_IMPOSE' Valeurs des flux imposés intégrés par face | * | (Type CHARGEMENT de CHPO Face, comp. 'FLUX'- | * | défaut 0.) | * | | * |'FLUXTOT_IMP' Valeurs des flux totaux imposés intégrés par face | * | (Type CHARGEMENT de CHPO Face, comp. 'FLUX'- | * | défaut 0.). Il s'agit du flux diffusif + convectif | * | | * |'MIXTES' Table : indice C contient les valeurs des flux | * | mixtes imposés intégrés par face | * | (Type CHARGEMENT de CHPO Face, comp. 'FLUX'- | * | défaut 0.) il est égal à A * flux diffusif + | * | B * Concentration | * | Indice A et B sont les réels A et B | * | | * |'SOURCE' Valeurs du terme source par maille et par unité de | * | temps (ex : puits, filiation) | * | Les valeurs à l'indice i sont les valeurs entre | * | les temps i-1 et i. | * | (CHARGEMENT de CHPO Centre, comp 'SOUR'- défaut 0.)| * | | * | | * |'DISSOLUTION_IMPOSEE' Valeurs des dissolutions imposées par unité| * | de temps et par maille. (Type CHARGEMENT de CHPO, | * | Comp 'H'). Les valeurs à l'indice i sont les | * | valeurs moyennes de dissolution par unité de temps | * | entre les temps i-1 et i. | * | Priorité de la dissolution imposée sur les | * | cinétiques. | * | | * |-------------------- | * |Données numériques : | * |-------------------- | * | | * | | * |'TEMPS_CALCULES' Valeur des temps calculés (LISTREEL) | * | Contient obligatoirement le temps final. | * | | * |'TEMPS_SAUVES' Valeur des temps sauvegardés (LISTREEL - défaut : | * | on sauve tous les pas de temps) | * | | * | | * | | * |'THETA_DIFF' Coefficient de relaxation compris entre 0. et 1. | * | (theta-méthode diffusion) ('FLOTTANT' - défaut 1.) | * | | * |'THETA_CONV' Idem pour la convection | * | ('FLOTTANT', Défaut = THETA_DIFF) | * |'THETA_DEC' Idem mais pour la décroissance | * | ('FLOTTANT' - défaut 1/2) | * | | * |'THETA_DISS' Idem mais pour la dissolution | * | ('FLOTTANT' - défaut 1.) | * | | * |'PENALISATION' Coefficient de pénalisation pour la prise en | * | compte de la limite de solubilité. La présence de | * | cet indice ou du suivant indique quel schéma a été | * | choisi. | * | (Type 'FLOTTANT') - Valeur conseillée 1.D7 | * | | * |'EPSI_LIM' Précision relative d'arrêt pour le shéma limite de | * | solubilité prédicteur-correcteur itératif | * | (Type FLOTTANT) - Valeur conseillée 5.D-3 | * | | * |'ITMAX_LIM' Nombre maxi d'itérations correspondant aux modules | * | de dissolution avant d'abandonner | * | (Type 'ENTIER') - Défaut 50 | * | | * |'EPSI_RET' Précision relative d'arrêt pour la résolution | * | itérative (Picard) de l'adsorption non linéaire | * | (Type FLOTTANT) - Défaut 1.D-4 | * | | * |'EPSI_COR' Petit saut de concentration pour calculer le coef. | * | de retard par la méthode de la corde lorsque le | * | retard est non-linéaire. | * | (Type FLOTTANT) - Défaut 1.D-4 | * | | * |'ITMAX_RET' Nombre maxi d'itérations correspondant au retard | * | non linéaire avant d'abandonner. | * | (Type 'ENTIER') - Défaut 20 | * |_________________________________________________________________| * * * *--------------------------------- * Résultats (stockés dans TABLE) : *--------------------------------- * * ___________________________________________________________________ * | | * | Indice Contenu | * | | * ------------------------------------------------------------------- * | | * | | * |'TEMPS' TABLE contenant les temps sauvegardés (FLOTTANT) | * | | * |'CONCENTRATION' TABLE contenant les concentrations | * | (Type CHPO Centre, Comp 'H') | * | | * |'FLUX' TABLE contenant les débits diffusifs intégrés | * | par face (Type CHPO Face, comp. 'FLUX') | * | | * |'PRECIPITE' TABLE contenant la quantité de précipité par maille| * | (Type CHPO Centre, Comp 'H') | * | | * |'DISSOLUTION' TABLE contenant la quantité de précipité dissoute | * | entre deux pas de temps par unité de volume et par | * | unité de temps. La valeur stockée à l'indice i, | * | est valable entre les temps i-1 et i | * | (Type CHPO, Comp 'H'). | * | ATTENTION, les valeurs de cette table résultat | * | n'ont aucun sens lorsque les temps sauvegardés ne | * | sont pas les memes que les temps calculés. Toute | * | tentative d'exploitation donnera alors des | * | résultats incohérents (erreurs de bilan) | * | | * |'RETARD' Si cet indice a été préalablement défini comme une | * | TABLE, alors il contient les valeurs du coefficient| * | de retard (Type 'CHPO' centre, Comp 'SCAL'). Sinon,| * | les valeurs du coefficient de retard ne sont pas | * | sauvegardées. | * |_________________________________________________________________| * * * ___________________________________________________________________ * | | * | Les tables résultats sont indicées par des entiers variant de 0 | * | à N . | * | A l'indice 0 on stocke les valeurs initiales, aux indices | * | suivants les champs correspondant au temps de sortie TEMPS.I . | * | Les champs servant en cas de reprise sont ceux correpondant au | * | dernier indice. | * |_________________________________________________________________| * * * TYPE DE RESOLUTION NOMINC = 'MOT' 'CONCENTRATION' ; NOMTETA = 'MOT' 'THETA_DIFF' ; NOMDDT = 'MOT' 'POROSITE' ; * *--------------------------------------------------------------------* * RECUPERATION DES DONNEES PHYSIQUES, GEOMETRIQUES ET MATERIELLES * *--------------------------------------------------------------------* * * MODELE 'SI' ( 'EXISTE' TRANSI 'MODELE' ) ; MODDARCY = TRANSI . 'MODELE' ; * METTRE LE PRECOND à 1 pour ne pas recalculer les préconditionnements * sur le maillage ('DOMA' MODDARCY TABLE) . 'PRECONDI' = 1 ; 'SINON' ; 'ERREUR' 'Il manque le modele.' ; 'QUITTER' TRANSGEN ; 'FINSI' ; * VOLUME VOLU1 = 'DOMA' MODDARCY 'VOLUME' ; * ORIENTATION MCHYB = 'DOMA' MODDARCY 'ORIENTAT' ; * MCENT = 'DOMA' MODDARCY 'CENTRE' ; * MFACE = 'DOMA' MODDARCY 'FACE' ; * MAILL = 'DOMA' MODDARCY 'MAILLAGE' ; * FILIATION 'SI' ('NON' ('EXISTE' TRANSI 'PERE')) ; TRANSI . 'PERE' = FAUX ; 'FINSI' ; * CARACTERISTIQUES 'SI' ( 'EXISTE' TRANSI 'CARACTERISTIQUES' ) ; 'SI' ('EGA' ('TYPE' TRANSI . 'CARACTERISTIQUES') 'CHPOINT') ; MAT1 = TRANSI . 'CARACTERISTIQUES' ; 'FINSI' ; * Si la diffusion est un chargement, le traitement à lieu dans le * corps du programme 'SINON' ; 'ERREUR' 'Indice CARACTERISTIQUES absent de la table de données.' ; 'QUITTER' TRANSGEN ; 'FINSI' ; * EMMAGASINEMENT OU POROSITE 'SI' ( 'EXISTE' TRANSI NOMDDT ) ; POROS = TRANSI . NOMDDT ; 'SINON' ; POROS = 1.D0 ; 'FINSI' ; 'SI' ('EGA' ('TYPE' POROS) 'FLOTTANT') ; POROS = 'MANU' 'CHPO' MCENT 'SCAL' POROS ; POROS = 'CHAN' 'ATTRIBUT' POROS 'NATURE' 'DISCRET' ; * VERRUE * POROS = 'MANU' 'CHPO' MCENT 'SCAL' 1.D0 NATURE DISCRET ; * 'FINSI' ; 'SI' ('NEG' ('TYPE' POROS) 'CHPOINT ') ; 'ERREUR' ('CHAINE' 'L indice ' NOMDDT ' doit etre de type ' 'FLOTTANT ou CHAMPOINT') ; 'SINON' ; POROS = 'NOMC' 'SCAL' POROS ; 'FINSI' ; *'SI' ((('MINI' POROS) < 0.D0) 'OU' (('MAXI' POROS) > 1.D0)) ; * 'ERREUR' 'La porosité doit etre comprise entre 0 et 1' ; *'FINSI' ; * RETARD 'SI' ( 'EXISTE' TRANSI 'COEF_RETARD' ) ; RETAR1 = TRANSI . 'COEF_RETARD' ; 'SI' ('EGA' ('TYPE' RETAR1) 'CHPOINT ') ; RETAR1 = NOMC 'SCAL' RETAR1 ; MINRET = ('MINI' RETAR1) ; 'SINON' ; MINRET = RETAR1 ; RETAR1 = 'MANU' 'CHPO' MCENT 'SCAL' MINRET ; RETAR1 = 'CHAN' 'ATTRIBUT' RETAR1 'NATURE' 'DISCRET' ; 'FINSI' ; 'SI' (MINRET < 1.D0) ; 'ERREUR' 'Le coefficient de retard doit etre supérieur à 1'; 'FINSI' ; 'SINON' ; * RETAR1 = kcht modhyb scal centre 1.D0 ; RETAR1 = 'MANU' 'CHPO' MCENT 'SCAL' 1.D0 ; RETAR1 = 'CHAN' 'ATTRIBUT' RETAR1 'NATURE' 'DISCRET' ; 'FINSI' ; RETAR1M1 = RETAR1 - 1. ; SAUVRET = 'EXISTE' TRANSI 'RETARD' ; * RETARD NON LINEAIRE RNONLINL = ('EXISTE' TRANSI 'COEF_RETARD') 'ET' ('EXISTE' TRANSI 'LANGMUIR') ; RNONLINF = ('EXISTE' TRANSI 'COEF_RETARD') 'ET' ('EXISTE' TRANSI 'FREUNDLICH') 'ET' ('NON' RNONLINL) ; RNONLIN = RNONLINL 'OU' RNONLINF ; * Cas Langmuir : 'SI' RNONLINL ; FSAT = TRANSI . 'LANGMUIR' ; 'SI' ('EGA' ('TYPE' FSAT) 'CHPOINT ') ; 'SI' (('MINI' FSAT) < 0.D0) ; 'ERREUR' 'La masse adsorbee a saturation doit etre positive' ; 'QUITTER' TRANSGEN ; 'FINSI' ; 'SINON' ; 'SI' (FSAT < 0.D0) ; 'ERREUR' 'La masse adsorbee a saturation doit etre positive' ; 'QUITTER' TRANSGEN ; 'FINSI' ; 'FINSI' ; RM1SURF = RETAR1M1 / FSAT ; 'FINSI' ; * Cas Freundlich : 'SI' RNONLINF ; UNSURN = TRANSI . 'FREUNDLICH' ; 'SI' ((UNSURN 1.D0)) ; 'ERREUR' 'L exposant de Freundlich doit etre compris entre 0' ' (exclu) et 1' ; 'QUITTER' TRANSGEN ; 'FINSI' ; 'FINSI' ; * DECROISSANCE DECRO = ( 'EXISTE' TRANSI 'DECROISSANCE' ) ; 'SI' DECRO ; LAMBD0 = TRANSI . 'DECROISSANCE' ; 'SI' (LAMBD0 < 0.D0) ; 'ERREUR' 'Le coefficient de décroissance doit etre positif.'; 'FINSI' ; 'SINON' ; LAMBD0 = 0. ; 'FINSI' ; * DISSOLUTION_IMPOSEE SOLUA = ('EXISTE' TRANSI 'DISSOLUTION_IMPOSEE' ) ; 'SI' SOLUA ; DISIMP = TRANSI . 'DISSOLUTION_IMPOSEE' ; 'FINSI' ; * LIMITE_SOLUBILITE * (Priorité de la dissolution imposée sur les autres processus) SOLUL = (('EXISTE' TRANSI 'LIMITE_SOLUBILITE' ) 'ET' ('NON' SOLUA)) ; 'SI' SOLUL ; 'SI' ('EGA' ('TYPE' TRANSI . 'LIMITE_SOLUBILITE') 'CHPOINT') ; LIMSOL = 'NOMC' 'H' TRANSI . 'LIMITE_SOLUBILITE' ; 'FINSI' ; 'SI' (('MINI' LIMSOL) < 0.D0) ; mess 'minimum de limsol' (MINI LIMSOL); 'ERREUR' 'La limite de solubilité doit etre positive.' ; 'FINSI' ; 'FINSI' ; * COEF_DISSOLUTION SOLUP = SOLUL 'ET' ('EXISTE' TRANSI 'COEF_DISSOLUTION' ) ; SOLUI = SOLUL 'ET' ('NON' SOLUP) ; 'SI' SOLUP ; codis = TRANSI . 'COEF_DISSOLUTION' ; 'SI' ('EGA' ('TYPE' CODIS) 'CHPOINT ') ; CODIS = 'NOMC' 'SCAL' TRANSI . 'COEF_DISSOLUTION' NATU DIFFUS ; 'SI' (('MINI' CODIS) < 0.D0) ; 'ERREUR' 'Le coefficient de dissolution doit etre positif.'; 'FINSI' ; 'SINON' ; CODIS = TRANSI . 'COEF_DISSOLUTION' ; CODIS = MANU CHPO MCENT SCAL CODIS ; 'SI' ((MINI CODIS) < 0.D0) ; 'ERREUR' 'Le coefficient de dissolution doit etre positif.'; 'FINSI' ; 'FINSI' ; 'FINSI' ; 'SI' (SOLUP) ; * Là où la précipitation est nulle (codis =0) on choisit * une limite de solubilité élevée à des fins d'optimisation * algorithmique 'SI' ('EGA' ('TYPE' TRANSI . 'LIMITE_SOLUBILITE') 'CHPOINT') ; dum = 'MASQUE' CODIS INFERIEUR 1.D-30 ; LIMSOL = ((1.D0 '-' dum) * LIMSOL) '+' ('NOMC' 'H' (1.D3*dum)) ; 'FINSI' ; 'FINSI' ; SOLUB = SOLUL 'OU' SOLUA ; * * ------------------------------------------------------------------ * * VERIFICATION DE LA BONNE STRUCTURE DES TABLES RESULTAT D'ORIGINE * * ------------------------------------------------------------------ * * * TEMPS 'SI' ( 'NON' ('EXISTE' TRANSI 'TEMPS' ) ) ; 'ERREUR' 'Indice TEMPS absent de la table de données.' ; 'QUITTER' TRANSGEN ; 'FINSI' ; * INCONNUE 'SI' ( 'NON' ('EXISTE' TRANSI NOMINC ) ) ; 'ERREUR' ('CHAINE' 'Indice ' NOMINC ' absent de la table de données.') ; 'QUITTER' TRANSGEN ; 'FINSI' ; * FLUXDIFF 'SI' ( 'NON' ('EXISTE' TRANSI 'FLUXDIFF' ) ) ; 'ERREUR' 'Indice FLUX diffusif absent de la table de données.' ; 'QUITTER' TRANSGEN ; 'FINSI' ; * FLUXCONV 'SI' ( 'NON' ('EXISTE' TRANSI 'FLUXCONV' ) ) ; 'ERREUR' 'Indice FLUX convectif absent de la table de données.' ; 'QUITTER' TRANSGEN ; 'FINSI' ; * DISSOLUTION, PRECIPITE 'SI' SOLUB ; * GBM mettre a jour dissolution . * 'SI' ( 'NON' ('EXISTE' TRANSI 'DISSOLUTION' ) ) ; * 'ERREUR' * 'Indice DISSOLUTION absent de la table de données.' ; * 'QUITTER' TRANSGEN ; * 'FINSI' ; 'SI' ( 'NON' ('EXISTE' TRANSI 'PRECIPITE' ) ) ; 'ERREUR' 'Indice PRECIPITE absent de la table de données.' ; 'QUITTER' TRANSGEN ; 'FINSI' ; 'FINSI' ; * TEST DES TAILLES DE TABLE IND1 = 'INDEX' ( TRANSI . 'TEMPS' ) ; IND3 = 'INDEX' ( TRANSI . NOMINC ) ; IND4 = 'INDEX' ( TRANSI . 'FLUXDIFF' ) ; LIN1 = 'DIME' IND1 ; LIN3 = 'DIME' IND3 ; LIN4 = 'DIME' IND4 ; 'SI' ( LIN1 'NEG' LIN3 ) ; 'ERREUR' ('CHAINE' 'Longueur des tables TEMPS et ' NOMINC ' différente.') ; 'QUITTER' TRANSGEN ; 'FINSI' ; 'SI' ( LIN1 'NEG' LIN4 ) ; 'ERREUR' 'Longueur des tables TEMPS et FLUXDIFF différente.' ; 'QUITTER' TRANSGEN ; 'FINSI' ; 'SI' SOLUB ; * GBM mettre a jour dissolution . * IND5 = 'INDEX' ( TRANSI . 'DISSOLUTION' ) ; IND6 = 'INDEX' ( TRANSI . 'PRECIPITE' ) ; * LIN5 = 'DIME' IND5 ; LIN6 = 'DIME' IND6 ; * 'SI' ( LIN1 'NEG' LIN5 ) ; * 'ERREUR' * 'Longueur des tables TEMPS et DISSOLUTION différente.' ; * 'QUITTER' TRANSGEN ; * 'FINSI' ; 'SI' ( LIN1 'NEG' LIN6 ) ; 'ERREUR' 'Longueur des tables TEMPS et PRECIPITE différente.' ; 'QUITTER' TRANSGEN ; 'FINSI' ; 'FINSI' ; * TEST DES INDICES DE TABLE IPO1 = 0 ; 'REPETER' BOU1 LIN1 ; IPO1 = IPO1 + 1 ; LAST1 = IND1 . IPO1 ; LAST3 = IND3 . IPO1 ; LAST4 = IND4 . IPO1 ; 'SI' SOLUB ; * GBM mettre a jour dissolution . * LAST5 = IND5 . IPO1 ; LAST6 = IND6 . IPO1 ; 'FINSI' ; 'SI' ( 'NEG' LAST1 LAST3 ) ; 'ERREUR' 'Indices des tables TEMPS et CONCENTRATION incohérents.' ; 'QUITTER' TRANSGEN ; 'FINSI' ; 'SI' ( 'NEG' LAST1 LAST4 ) ; 'ERREUR' 'Indices des tables TEMPS et FLUX incohérents.' ; 'QUITTER' TRANSGEN ; 'FINSI' ; 'SI' SOLUB ; * GBM mettre a jour dissolution . * 'SI' ( 'NEG' LAST1 LAST5 ) ; * 'ERREUR' * 'Indices des tables TEMPS et DISSOLUTION incohérents.'; * 'QUITTER' TRANSGEN ; * 'FINSI' ; 'SI' ( 'NEG' LAST1 LAST6 ) ; 'ERREUR' 'Indices des tables TEMPS et PRECIPITE incohérents.' ; 'QUITTER' TRANSGEN ; 'FINSI' ; 'FINSI' ; 'FIN' BOU1 ; * * ------------------------------------------------------------------ * * RECUPERATION DES CONDITIONS INITIALES OU DU DERNIER PAS SAUVE * * ------------------------------------------------------------------ * * TPSINI = TRANSI . 'TEMPS' . LAST1 ; *Momentanée pour la trace * On extrait les noms de composante des espèces NOMESP = 'EXTRAIRE' ('EXTRAIRE' TRANSI . NOMINC . LAST1 'COMP') 1 ; * nomespl est un nom d'espece locale à la procédure car * les EFMH veulent voir une composante 'H' NOMESPL = 'MOT' 'H' ; CHRG = 'NOMC' 'H' TRANSI . NOMINC . LAST1 ; *VERRUE FLU0 = 'NOMC' 'H' TRANSI . 'FLUXDIFF' . LAST1 ; FLUCO0 = 'NOMC' 'H' TRANSI . 'FLUXCONV' . LAST1 ; * *--------------------------------------------------------------------* * RECUPERATION DES THETA SCHEMAS DIFFUSION-CONVECTION * *--------------------------------------------------------------------* * THETA DIFFUSION 'SI' ( 'EXISTE' TRANSI 'THETA_DIFF' ) ; TETA = TRANSI . 'THETA_DIFF' ; 'SINON' ; TETA = 1.D0 ; 'FINSI' ; * THETA CONVECTION 'SI' ( 'EXISTE' TRANSI 'THETA_CONVECTION' ) ; TETAC = TRANSI . 'THETA_CONVECTION' ; 'SINON' ; TETAC = TETA ; 'FINSI' ; * *--------------------------------------------------------------------* * RECUPERATION DES CONDITIONS AUX LIMITES ET DES CHARGEMENTS * *--------------------------------------------------------------------* * * CL de type Dirichlet : BLOCAGE * TRACE_IMPOSE * On débranche les tests sur les longueurs de flux non nuls * en CL pour les cas ou ils ne sont pas intégrés en temps startflu = 1.D30; startflt = 1.D30; startmix = 1.D30; 'SI' ( 'EXISTE' TRANSI 'TRACE_IMPOSE') ; CHDIRI = TRANSI . 'TRACE_IMPOSE' ; 'FINSI' ; * CL de type Neumann : FLUX_IMPOSE 'SI' ( 'EXISTE' TRANSI 'FLUX_IMPOSE' ) ; FLUIMP = TRANSI . 'FLUX_IMPOSE' ; 'SI' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'VF') 'OU' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'EFMH') 'ET' ('EGA' (teta * tetac) 1.D0))) ; * On écrit les flux sous forme intégrale sauf en explicite et * kranck-Nickholson pour EFMH FLUIMP startflu = CHAMINT FLUIMP ; 'FINSI' ; 'FINSI' ; * CL de type flux total: FLUXTOT_IMP 'SI' ( 'EXISTE' TRANSI 'FLUXTOT_IMP' ) ; FLTOTIMP = TRANSI . 'FLUXTOT_IMP' ; 'SI' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'VF') 'OU' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'EFMH') 'ET' ('EGA' (teta * tetac) 1.D0))) ; * On écrit les flux sous forme intégrale sauf en explicite et * kranck-Nickholson pour EFMH FLTOTIMP startflt = CHAMINT FLTOTIMP ; 'FINSI' ; 'FINSI' ; * CL de type mixte : MIXTES 'SI' ( 'EXISTE' TRANSI 'MIXTES' ) ; FLUMIX = TRANSI . 'MIXTES' ; 'SI' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'VF') 'OU' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'EFMH') 'ET' ('EGA' (teta * tetac) 1.D0))) ; * On écrit les flux sous forme intégrale sauf en explicite et * kranck-Nickholson pour EFMH FLUMIX startmix = CHAMINT FLUMIX ; 'FINSI' ; 'FINSI' ; 'SI' ( 'EXISTE' TRANSI 'SOURCE' ) ; TERSOU = TRANSI . 'SOURCE' ; * on crée une evolution intégrée temporellement pour assurer la * conservation lors de liste de pas de temps de chargement differente * des temps de discrétisation TERSOU startsou = CHAMINT TERSOU ; 'FINSI' ; * *--------------------------------------------------------------------* * RECUPERATION DES données NUMERIQUES * *--------------------------------------------------------------------* * * PARAMETRES DE SOLUBILITE 'SI' (SOLUI 'OU' SOLUP) ; 'SI' ( 'EXISTE' TRANSI 'ITMAX_LIM') ; ITMAXI = TRANSI . 'ITMAX_LIM' ; 'SINON' ; ITMAXI = 50 ; 'FINSI' ; 'FINSI' ; 'SI' SOLUI ; 'SI' ( 'EXISTE' TRANSI 'PENALISATION') ; PENAL = TRANSI . 'PENALISATION' ; MET0 = 'PENALISATION' ; 'SINON' ; 'SI' ( 'EXISTE' TRANSI 'EPSI_LIM') ; EPS1 = TRANSI . 'EPSI_LIM' ; MET0 = 'PF DISSOLUTION' ; 'SINON' ; 'MESS' 'Il manque le paramètre ' ; 'ERREUR' 'du schéma numérique de limite de solubilité'; 'FINSI' ; 'FINSI' ; 'FINSI' ; * PARAMETRES DE RETARD NON LINEAIRE 'SI' RNONLIN ; 'SI' ( 'EXISTE' TRANSI 'EPSI_RET' ) ; EPSRNL = TRANSI . 'EPSI_RET' ; 'SINON' ; EPSRNL = 1.D-4 ; 'FINSI' ; 'SI' ( 'EXISTE' TRANSI 'ITMAX_RET' ) ; ITMAXRNL = TRANSI . 'ITMAX_RET' ; 'SINON' ; ITMAXRNL = 20 ; 'FINSI' ; * Petit saut de concentration minimal pour le calcul de dérivée 'SI' ( 'EXISTE' TRANSI 'EPSI_COR' ) ; EPSCORD = TRANSI 'EPSI_COR' ; 'SINON' ; 'SI' RNONLINL ; EPSCORD = ('MANU' 'CHPO' MCENT 1 'SCAL' 1.D-4) / RETAR1 * FSAT ; 'FINSI' ; 'SI' RNONLINF ; EPSCORD = 1.D-4 ; 'FINSI' ; 'FINSI' ; * Pour les isothermes linéaires la boucle de prise en compte de * retard non linéaire n'est parcourue qu'une fois : NPicard = ITMAXRNL ; 'SINON' ; NPicard = 1 ; 'FINSI' ; * THETA_DECROISSANCE 'SI' ( 'EXISTE' TRANSI 'THETA_DEC' ) ; BETA = TRANSI . 'THETA_DEC' ; 'SINON' ; BETA = 0.5 ; 'FINSI' ; * THETA_DISSOLUTION 'SI' ( 'EXISTE' TRANSI 'THETA_DIS' ) ; GAMMA = TRANSI . 'THETA_DIS' ; 'SINON' ; GAMMA = 1. ; 'FINSI' ; * TEMPS_CALCULES 'SI' ( 'EXISTE' TRANSI 'TEMPS_CALCULES' ) ; TPCAL = 'ORDO' (TRANSI . 'TEMPS_CALCULES') ; TRANSI . 'TEMPS_CALCULES' = TPCAL ; DCAL = 'DIME' TPCAL ; TPSFIN = 'EXTR' TPCAL DCAL ; TPS0 = 'EXTR' TPCAL 1 ; 'SI' ('NEG' DCAL 1) ; TPSANT = 'EXTR' TPCAL (DCAL - 1) ; TPS1 = 'EXTR' TPCAL 2 ; DT10 = TPS1 - TPS0 ; DT21 = TPSFIN - TPSANT ; PROPS = 'PROG' DT10 DT21 ; DTEPS = 'MINIMUM' PROPS ; 'SINON' ; DTEPS = TPSFIN - TPS0 ; 'FINSI' ; * EPS0 = DTEPS * 1.D-6 ; * ICAL est le premier indice des temps à calculer supérieur à * TpsIni + eps IOK1 = 0 ; 'REPETER' BOU2 DCAL ; ICAL = &BOU2 ; TEMS = 'EXTR' TPCAL ICAL ; 'SI' ( TPSINI '<' (TEMS - EPS0) ) ; IOK1 = 1 ; 'QUITTER' BOU2 ; 'FINSI' ; 'FIN' BOU2 ; 'SI' ( IOK1 'EGA' 0) ; 'MESS' 'Listreel des temps de calcul inférieur à tini.' ; 'MESS' 'On ne fait donc rien. ' ; 'QUITTER' TRANSGEN ; 'FINSI' ; DELTAT0 = TEMS - TPSINI ; 'SINON' ; 'ERREUR' 'Indice TEMPS_CALCULES absent de la table de données.' ; 'QUITTER' TRANSGEN ; 'FINSI' ; * TEMPS_SAUVES 'SI' ( 'EXISTE' TRANSI 'TEMPS_SAUVES' ) ; TPSOR = 'ORDO' (TRANSI . 'TEMPS_SAUVES') ; TRANSI . 'TEMPS_SAUVES' = TPSOR ; DSOR = 'DIME' TPSOR ; ISOR = 0 ; IOK1 = 0 ; 'REPETER' BOU3 DSOR ; ISOR = ISOR + 1 ; TEMS = 'EXTR' TPSOR ISOR ; 'SI' ( ( TPSINI '<' (TEMS - EPS0) ) 'ET' ( TEMS 'EG 1) ; 'MESS' ' ' NB0 ' itérations supp de solubilité' ; 'FINSI' ; * 'Fin module Dissolution d ordre 1' 'FINSI' ; *|-------------------------------------| *| Retard Non Linéaire | *|-------------------------------------| *- Calcul du nouveau coefficient de retard et de l'adsorbât 'SI' RNONLIN ; * 1. F(C) CC0 = 'NOMC' CHRG 'SCAL' ; 'SI' RNONLINL ; * Langmuir FF = RETAR1M1 * CC0 / (1.D0 + ( RM1SURF * CC0 )) ; 'SINON' ; * Freundlich * La bidouille utilisant les masques sert à obtenir le champ * des signes de CC0. a = 'MASQUE' CC0 'SUPERIEUR' 0. ; b = 'MASQUE' CC0 'INFERIEUR' 0. ; FF = RETAR1 * (('ABS' CC0) ** UNSURN) * (a - b) ; 'DETRUIT' a ; 'DETRUIT' b ; 'FINSI' ; FF0 = FF ; * 2. F(C + dC) CC1 = ('NOMC' CHA2 'SCAL') + EpsCord ; 'SI' RNONLINL ; * Langmuir FF1 = RETAR1M1 * CC1 / (1.D0 + ( RM1SURF * CC1 )) ; 'SINON' ; * Freundlich a = 'MASQUE' CC1 'SUPERIEUR' 0. ; b = 'MASQUE' CC1 'INFERIEUR' 0. ; FF1 = RETAR1 * (('ABS' CC1) ** UNSURN) * (a - b) ; 'DETRUIT' a ; 'DETRUIT' b ; 'FINSI' ; * 3. R = 1 + dF/dC RETAR2 = 1.D0 + ( (FF1 - FF0) / (CC1 - CC0) ) ; 'DETRUIT' CC1 ; 'DETRUIT' FF1 ; * *- Critère de sortie de la boucle de retard non linéaire : CRIT = ( 'ABS' (RETAR2 - RETARC) ) / RETARC ; CRIT = 'MAXIMUM' CRIT ; 'SI' (CRIT < EPSRNL) ; 'QUITTER' BOURNL ; 'SINON' ; * Le coef. de retard ne change que si le critère est mauvais. * Réévaluation des variables affectées devant DC/DT: RETARC = RETAR2 ; PORETSU = RETARC * POROS ; * coef devant dc/dt modifié TABMODI . 'POROSITE' = VRAI ; * Fonctions de volume réactualisées PORETSU = RETARC * POROS ; DECROI0 = 'EXCO' ((-1.) * LAMBD0 * PORETSU ) 'SCAL' 'SCAL'; DECROI0 = VOLU1 * DECROI0 ; 'FINSI' ; 'SINON' ; * NE sert a rien GBM * Retard linéaire F = (R - 1) C * zozo = 'NOMC' CHRG 'SCAL' ; * FF = RETAR1M1 * zozo ; * 'DETRUIT' zozo ; 'FINSI' ; * Si on atteint le nombre maximum d'itérations - en particulier si le * retard est linéaire - on quitte la boucle : 'SI' (IRNL 'EGA' ITMAXRNL) ; MESS ' sortie boucle de retard non linéaire sans convergence' ; 'QUITTER' BOURNL ; 'FINSI' ; *------------- 'FIN' BOURNL ; *------------- 'SI' (IRNL > 1 ) ; 'MESS' ' ' IRNL ' itérations de retard non linéaire.' ' Critère = ' crit ; 'FINSI' ; NBITRNL = NBITRNL + (IRNL - 1) ; *|-----------------------------------| *| Fin des modules de résolution | *|-----------------------------------| * *- Calcul de la dissolution par unité de volume et par unité de temps : * 'SI' SOLUA ; DIS2 = DIS1 / VOLU1 / DELTAT ; 'FINSI' ; * *- Archivage des resultats * * si le prochain temps à sauver tombe entre les deux * incréments de temps, on obtient les valeurs à sauver * par interpolation linéaire puis on les stocke. * *list ((PLACE) * 4. / 1000000.) ; 'SI' ( 'EXISTE' TRANSI 'TEMPS_SAUVES' ) ; * Sauvegarde de tous les temps intermédiaires, s'il y en a : 'REPETER' BOU9 ; TEMS = 'EXTR' TPSOR ISOR ; 'SI' ( ( PRECED '<' TEMS ) 'ET' ( TEMS '' DSOR ) ; ISOR = DSOR ; 'QUITTER' BOU9 ; 'FINSI' ; 'FIN' BOU9 ; 'SINON' ; * sinon, on sauvegarde tous les temps : LAST1 = LAST1 + 1 ; TRANSI . 'TEMPS' . LAST1 = TPS ; TRANSI . NOMINC . LAST1 = 'NOMC' NOMESP CHA2 ; TRANSI . 'FLUXDIFF' . LAST1 = 'NOMC' NOMESP FLU2 ; TRANSI . 'FLUXCONV' . LAST1 = 'NOMC' NOMESP FLUCO2 ; 'SI' SOLUB ; TRANSI . 'PRECIPITE' . LAST1 = 'NOMC' NOMESP PRE2 ; 'FINSI' ; 'SI' SAUVRET ; TRANSI . 'RETARD' . LAST1 = 'NOMC' NOMESP RETARC ; 'FINSI' ; 'FINSI' ; *On sauve l'intégrale de decroissance*concentration et * decroissance*précipité au cours du temps * stockée pour l'instant à tous les temps. Peut s'améliorer. * indice pour le stockage des intégrales de concentration LASTINT = LASTINT '+' 1; BM1 = 1.D0 '-' BETA ; * si l'espece fait l'objet d'une filiation on stocke les intégrales * du précipité et de la concentration. SI ((EXISTE TRANSI DECROISSANCE) 'ET' (TRANSI . 'PERE')) ; *TRANSI . 'INTFILI' . LASTINT = TRANSI . 'INTFILI' . (LASTINT '-' 1) * '+' (DELTAT * PORETSU * TRANSI . 'DECROISSANCE' * * ((BM1 * ('NOMC' NOMESP CHRG)) * '+' ((BETA) * ('NOMC' NOMESP CHA2)))) ; zozo = (DELTAT * TRANSI . 'DECROISSANCE') * PORETSU ; zaza = BM1 * CHRG ; zuzu = BETA * CHA2 ; zeze = 'KOPS' zaza '+' zuzu ; 'DETRUIT' zaza ; 'DETRUIT' zuzu ; zaza = 'KOPS' zozo * zeze ; 'DETRUIT' zozo ; 'DETRUIT' zeze ; zozo = 'NOMC' NOMESP zaza ; 'DETRUIT' zaza ; TRANSI . 'INTFILI' . LASTINT = TRANSI . 'INTFILI' . (LASTINT '-' 1) '+' zozo ; 'DETRUIT' zozo ; 'SI' SOLUB ; * TRANSI . 'INTFILI' . LASTINT = TRANSI . 'INTFILI' . LASTINT * '+' (DELTAT * (1.D0 '-' POROS) * TRANSI . 'DECROISSANCE' * * ((BM1 * ('NOMC' NOMESP PRCI)) * '+' ((BETA) * ('NOMC' NOMESP PRE2)))) ; zizi = (1.D0 '-' POROS) ; zozo = (DELTAT * TRANSI . 'DECROISSANCE') * zizi ; 'DETRUIT' zizi ; zaza = BM1 * PRCI ; zuzu = BETA * PRE2 ; zeze = 'KOPS' zaza '+' zuzu ; 'DETRUIT' zaza ; 'DETRUIT' zuzu ; zaza = 'KOPS' zozo * zeze ; 'DETRUIT' zozo ; 'DETRUIT' zeze ; zozo = 'NOMC' NOMESP zaza ; 'DETRUIT' zaza ; zaza = TRANSI . 'INTFILI' . LASTINT '+' zozo ; 'DETRUIT' TRANSI . 'INTFILI' . LASTINT ; 'DETRUIT' zozo ; TRANSI . 'INTFILI' . LASTINT = zaza ; 'OUBLIER' zaza ; 'FINSI' ; FINSI ; * Sauvegarde du dernier pas de temps si ça n'a pas déjà été fait. * Le dernier pas de temps ne peut pas faire l'objet d'interpolation * linéaire si le temps sauvegardé ne tombe pas pile dessus. 'SI' (IPAS 'EGA' DCAL) ; TPSDER = TRANSI. 'TEMPS' . LAST1 ; 'SI' ( TPSDER 'NEG' TPSFIN EPSDT ) ; LAST1 = LAST1 + 1 ; TRANSI . 'TEMPS' . LAST1 = TPS ; TRANSI . NOMINC . LAST1 = 'NOMC' NOMESP CHA2 ; TRANSI . 'FLUXDIFF' . LAST1 = 'NOMC' NOMESP FLU2 ; TRANSI . 'FLUXCONV' . LAST1 = 'NOMC' NOMESP FLUCO2 ; 'SI' SOLUB ; TRANSI . 'PRECIPITE' . LAST1 = 'NOMC' NOMESP PRE2 ; 'FINSI' ; 'SI' SAUVRET ; TRANSI . 'RETARD' . LAST1 = 'NOMC' NOMESP RETARC ; 'FINSI' ; 'FINSI' ; 'FINSI' ; * *- Initialisations pour le pas suivant * PRECED = TPS ; 'DETRUIT' CHRG ; CHRG = CHA2 ; 'DETRUIT' FLU0 ; FLU0 = FLU2 ; 'DETRUIT' FLUCO0 ; FLUCO0 = FLUCO2 ; 'OUBLIER' GEOL1 ; 'OUBLIER' GEOL2 ; GEOL1 = 'TABLE' GEOLPF1 ; GEOL2 = 'TABLE' GEOLPF2 ; 'DETRUIT' RETAR0 ; RETAR0 = RETARC ; DELOLD = DELTAT ; 'SI' SOLUB ; 'DETRUIT' PRCI ; PRCI = PRE2 ; 'FINSI' ; 'SI' ('EXISTE' TRANSI 'SOURCE') ; 'DETRUIT' TERSC2M1 ; TERSC2M1 = TERSC2 ; 'FINSI' ; 'SI' ('EXISTE' TRANSI 'FLUX_IMPOSE') ; 'DETRUIT' FLUIMPM1 ; FLUIMPM1 = FLUIMPO ; 'FINSI' ; 'SI' ('EXISTE' TRANSI 'FLUXTOT_IMP') ; 'DETRUIT' FLUTMPM1 ; FLUTMPM1 = FLUTMPO ; 'FINSI' ; 'SI' ('EXISTE' TRANSI 'MIXTES') ; 'DETRUIT' FLUMMPM1 ; FLUMMPM1 = FLUMMPO ; 'FINSI' ; *MESS 'TEMPS2' ; TEMPS; *'MENAGE' ; * destruction d'objets temporaires 'DETRUIT' TERSC3 ; *list ((PLACE) * 4. / 1000000.) ; *temps impr; * *============= 'FIN' BOUTPS ; COMM 'Boucle sur le temps' ; *============= * *- Indication du nombre de passages : * 'SI' (SOLUI 'ET' ('EGA' MET0 'PENALISATION')) ; 'MESS' 'Nombre total de matrices calculées : ' NBIT ; 'FINSI' ; 'SI' (SOLUI 'ET' ('EGA' MET0 'PF DISSOLUTION')) ; 'MESS' 'Nombre total d itérations supplémentaires de solubilité : ' NBIT ; 'FINSI' ; 'SI' SOLUP ; 'MESS' 'Nombre total d itérations supplémentaires de solubilité : ' NBIT; 'FINSI' ; 'SI' RNONLIN ; 'MESS' 'Nombre total d itérations supplémentaires de ' 'retard non linéaire : ' NBITRNL; 'FINSI' ; 'FINP' ; $$$$