* fichier : fluxtotalVF.dgibi * ********************* CAS TEST : transport1.dgibi ******************** * GRAPH = 'N' ; 'OPTI' 'ECHO' 1 ; CRIT1 = 1.D-6 ; 'SAUT' 'PAGE' ; * *--------------------------------------------------------------------- * TEST TRANSPORT1 * CALCUL DARCY ISOTROPE TRANSPORT. * Transport d'un front. * * dT * -- + div (uT - Kgrad(T)) = 0 * dt * * Ce test permet de vérifier le bon fonctionnement des opérateurs * utilisés afin de résoudre l'équation de transport par diffusion * convection d'un champ scalaire passif par la méthode d'éléments * finis mixtes hybrides implanté dans CASTEM2000. * * Le maillage ne comporte qu'un élément dans la direction transverse * à l'écoulement puisque le phénomène étudié est 1D. Dans le sens de * l'écoulement, le pas de discrétisation est constant. * * Les conditions initiales correspondent à l'entrée du front dans le * domaine. Comme le champ de vitesse choisi est positif, la trace de * concentration est égale à 1 à gauche du domaine. * * Les conditions aux limites imposent la trace de concentration à * gauche du domaine égale à un via une condition de flux total * (soit U C = U pour avoit C = 1) ; sur les autres frontières ce sont * les conditions aux limites naturelles. * * Le schéma en temps utilisé est le schéma d'Euler implicite * (teta-méthode avec teta=1.0). * * Données physiques : * VK : Coefficient de diffusion ................. 1.D-2 * VX1 : Vitesse suivant x......................... 1. * VY1 : Vitesse suivant y......................... 0. * T0 : Concentration initiale.................... 0. * T1 : Concentration en entrée................... 1. * * Données temporelles : * TMIN : Temps initial............................. 0. * TMAX : Temps final............................... 30. * DELTAT: Pas de temps.............................. 1. * T0 : Concentration initiale.................... 0. * T1 : Concentration en entrée................... 1. * * Parametres de maillage : * L : longueur du domaine ....................... 100. * H : hauteur du domaine ........................ 1. * INUMX : Nombre d'éléments en x .................... 100 * INUMY : Nombre d'éléments en y .................... 1 car 1D * *--------------------------------------------------------------------- * *------------------ * Options generales *------------------ * 'OPTI' 'DIME' 2 'ELEM' 'QUA4' ; 'OPTI' 'ISOV' 'SURF' ; * * *========= * MAILLAGE *========= * * *- Création des points supports du contour du domaine, et des droites *- passant par les centres et les faces pour le post-traitement. * L = 100.D0 ; LS2 = L / 2.D0 ; H = 10.D0 ; HS2 = H / 2.D0 ; X0 = 0.D0 ; X1 = X0 + L ; Y0 = 0.D0 ; Y1 = Y0 + H ; INUMX = 100 ; INUMY = 1 ; INUM1 = INUMX - 1 ; Y01 = Y0 + Y1 * 0.5D0 ; DX = X1 - X0 / INUMX ; DX1 = DX / 2.D0 ; XG = X0 + DX1 ; XD = X1 - DX1 ; * A1 = X0 Y0 ; A3 = X1 Y0 ; D1 = X0 Y1 ; D3 = X1 Y1 ; B1 = X0 Y01 ; B3 = X1 Y01 ; C1 = XG Y01 ; C3 = XD Y01 ; P6 = LS2 Y01 ; * *- Création des DROITES frontieres * DRBAS = A3 'DROI' INUMX A1 ; DRGAU = A1 'DROI' INUMY D1 ; DRHAU = D1 'DROI' INUMX D3 ; DRDRO = D3 'DROI' INUMY A3 ; PELIM = DX1 / (5. * INUMX) ; * *- Creation maillage GEOMETRIQUE * PTOT1 = 'DALL' DRBAS DRGAU DRHAU DRDRO ; PTOT2 = 'ORIE' PTOT1 ; * *- Creation maillage HYBRIDE y compris sous-objets (cond. limites) * DRMID = B1 'DROI' INUMX B3 ; DRMIC = C1 'DROI' INUM1 C3 ; EXT1 = 'MANU' 'POI1' B1 ; * QFTOT= 'CHAN' PTOT2 QUAF ; QFGAU= 'CHAN' DRGAU QUAF ; QFDRO= 'CHAN' DRDRO QUAF ; ELIM PELIM (QFTOT ET QFGAU ET QFDRO ET DRMID ET DRMIC ET EXT1) ; * *================ * INITIALISATIONS *================ * * ---------------- * = MODELISATION = * ---------------- MODHYB = MODE QFTOT 'DARCY' 'ANISOTROPE' ; MODDRO = MODE QFDRO 'DARCY' 'ANISOTROPE' ; MODGAU = MODE QFGAU 'DARCY' 'ANISOTROPE' ; CHYB1 = 'DOMA' MODHYB 'SURFACE' ; CHYB2 = 'DOMA' MODHYB 'NORMALE' ; * * --------------------- * = Donnees physiques = * --------------------- * T0 = 0.D0 ; CTT1 = 1.D0 ; T1 = (H ) / INUMY; VX1 = 1.D0 ; VY1 = 0.D0 ; VK = 1.D-2 ; MATI2 = 'MATE' MODHYB DIRECTION (1. 0.) 'K11' VK 'K21' 0.D0 'K22' VK; MATI2 = KCHA MODHYB MATI2 'CHPO'; * SPEED = 'MANU' 'CHPO' ('DOMA' MODHYB 'FACE') 2 'VX' VX1 'VY' VY1 'NATURE' 'DISCRET' ; * * ----------------------- * = Donnees transitoire = * ----------------------- * TETA : Parametre de le theta-méthode * TMAX : Temps final * TSUP : Temps pour conditions aux limites * DELTAT : Pas de temps * TETA = 0.00D0 ; TMIN = 0.D0 ; TMAX2 = 15.00D0 ; TMAX = 30.00D0 ; TSUP = 1.2D0 * TMAX ; DELTAT = 0.5D0 ; * LICALC = 'PROG' TMIN 'PAS' (DELTAT/2.D0) TMAX2 ; LICALC = LICALC ET ('PROG' (TMAX2 + DELTAT) 'PAS' DELTAT TMAX) ; LISAUV = 'PROG' TMAX ; * * ------------------------ * = Conditions initiales = * ------------------------ TFHY0 = 'MANU' 'CHPO' ('DOMA' MODHYB 'FACE') 1 'TH' T0 'NATURE' 'DISCRET' ; TFHY1 = 'MANU' 'CHPO' EXT1 1 'TH' T1 'NATURE' 'DISCRET' ; TFHYB = TFHY0 'ET' TFHY1 ; TCHYB = 'MANU' 'CHPO' ('DOMA' MODHYB 'CENTRE') 1 'H' T0 'NATURE' 'DISCRET' ; * * -------------- * = T imposée = * -------------- BBGAU = 'BLOQ' ('DOMA' MODGAU 'CENTRE') 'TH' ; EEGAU = 'DEPI' BBGAU T1 ; TT1 = 'KCHT' MODGAU SCAL CENTRE (-1.* T1); TTT1 = 'CHAR' TT1 ('EVOL' 'MANU' ('PROG' 0. TSUP) ('PROG' 1. 1.)) ; CHAGAU = 'CHAR' EEGAU ('EVOL' 'MANU' ('PROG' 0. TSUP) ('PROG' 1. 1.)) ; * * --------------------------- * = Table DARCY_TRANSITOIRE = * --------------------------- * * * *-- Table de transport : Transp = 'TABLE'; Transp . 'MODELE' = MODHYB ; Transp.'TEMPS' = 'TABLE'; Transp.'CONCENTRATION' = 'TABLE'; Transp.'FLUXDIFF' = 'TABLE'; Transp.'FLUXCONV' = 'TABLE'; Transp.'CARACTERISTIQUES' = MATI2; *Transp.'POROSITE' = MaPor; Transp.'CONVECTION' = SPEED ; Transp.'VITELEM' = kcht MODHYB VECT CENTRE (1.D0 0.); Transp.'VITELEM' = nomc (mots UX UY) (mots VX VY) Transp . 'VITELEM'; * Conditions initiales : Transp.'TEMPS'. 0 = TMIN ; Transp.'CONCENTRATION'. 0 = TCHYB; Transp.'FLUXDIFF'. 0 = 0.D0 * TFHY0; Transp.'FLUXCONV'. 0 = 0.D0 * TFHY0; * Conditions aux limites : *Transp . 'TRACE_IMPOSE' = TTT1 ; Transp . 'FLUXTOT_IMP' = TTT1 ; * Paramètres numériques : Transp.'THETA_DIFF' = 1.D0; Transp.'THETA_CONVECTION' = 1.0D0; Transp.'LUMP' = FAUX; Transp.'TYPDISCRETISATION' = 'VF'; Transp.'DECENTR' = FAUX; TABRES = table METHINV; TABRES . 'TYPINV' = 1; TABRES . 'PRECOND' = 3; Transp . 'METHINV' = TABRES; Transp.'TEMPS_CALCULES' = LiCalc; Transp.'TEMPS_SAUVES' = LiSauv; Transp . INTCONC = TABLE; Transp . INTCONC . 0 = 0.D0 * TCHYB; * ========== *Transp . 'DECROISSANCE' = 0.01D0; * | CALCUL | * ========== *======================= * Resolution transitoire *======================= * TRANSGEN TRANSP ; * * *================= * POST-TRAITEMENT *================= TRANS2 = TABLE TRANSP; * *-------------------------------------------------------------------- * Dans chaque cas on trace * La trace de concentration le long de DRMID * La concentration le long de DRMIC *-------------------------------------------------------------------- * Tests de NON-REGRESSION : * Principe du maximum * Position du centre de gravité du champ de concentration *-------------------------------------------------------------------- * * Critères numériques CFL = VX1 * DELTAT / DX ; PEK = VX1 * DX / (2. * VK) ; FOU = 2. * VK * DELTAT / (DX * DX) ; 'SAUT' 1 'LIGNE' ; 'MESS' ' Critères numériques ' ; 'MESS' ' CFL ' CFL ; 'MESS' ' PECLET ' PEK ; 'MESS' ' FOURIER ' FOU ; 'MESS' ' ' ; * XC YC = 'COOR' (DOMA MODHYB 'CENTRE') ; ISOR1 = INDEX ( TRANS2 . 'TEMPS') ; NTSOR = DIME ISOR1 ; NTSO1 = NTSOR - 1 ; IOK = FAUX ; IRESU = 1 ; * *----------------------- REPETER VISURESU NTSO1 ; *----------------------- * IRESU = IRESU + 1 ; INDI1 = ISOR1.IRESU ; TTRA = TRANS2 . 'TEMPS' . INDI1 ; * Principe du maximum EPR1 = TRANS2 . 'CONCENTRATION' . INDI1 ; PMA = 'MAXI' EPR1 ; PMI = 'MINI' EPR1 ; VTEST = PMA + PMI / 1.D0 ; 'SAUT' 1 'LIGNE' ; 'MESS' ' Principe du maximum ' ; 'MESS' ' Max et min théorique ' T1 T0 ; 'MESS' ' Max et min par maille ' PMA PMI ; 'MESS' ' Précision demandée ' CRIT1 ; 'MESS' ' ' ; * Centre de gravité M0 = 'KCHA' MODHYB 'CHAM' EPR1 ; M1 = 'INTG' MODHYB M0 ; * XV1 = XC * EPR1 ; MXV1 = 'KCHA' MODHYB 'CHAM' XV1 ; XV2 = 'INTG' MODHYB MXV1 ; XV3 = XV2 / M1 ; YV1 = YC * EPR1 ; MYV1 = 'KCHA' MODHYB 'CHAM' YV1 ; YV2 = 'INTG' MODHYB MYV1 ; YV3 = YV2 / M1 ; XTEST = TTRA / 2.D0 ; * 'SAUT' 1 'LIGNE' ; 'MESS' ' Centre de gravité du nuage ' ; 'MESS' ' THEORIQUE ' XTEST HS2 ; 'MESS' ' OBTENU ' XV3 YV3 ; 'MESS' ' TAILLE D UNE MAILLE ' DX ; 'MESS' ' ' ; 'SAUT' 1 'LIGNE' ; * Tests de non regression 'SI' ( VTEST 'NEG' CTT1 CRIT1 ) ; IOK = VRAI ; 'FINS' ; 'SI' ( XV3 'NEG' XTEST DX ) ; IOK = VRAI ; 'FINS' ; 'SI' ( YV3 'NEG' HS2 CRIT1 ) ; IOK = VRAI ; 'FINS' ; * 'SI' ('NEG' GRAPH 'N' ) ; LTI2 = 'CHAINE' 'Front 1D-h temps ' TTRA ; 'TITR' LTI2 ; AV2 = 'EVOL' 'ROUG' 'CHPO' EPR1 'H' DRMIC ; 'DESS' AV2 'MIMA' ; 'FINS' ; * *------------- FIN VISURESU ; *------------- * * 'SI' ( IOK ) ; mess 'erreur' ((VTEST - T1) / T1); mess 'erreur' ((XV3 - XTEST) / XTEST); mess 'erreur' ((YV3 - HS2) / HS2); 'ERRE' 5 ; 'SINO' ; 'ERRE' 0 ; 'FINSI' ; * 'FIN' ;