* fichier : darcy4.dgibi ************************************************************************ * Section : Fluides Darcy ************************************************************************ * ************************** CAS TEST : darcy4.dgibi ****************** * GRAPH = 'N' ; 'SAUT' 'PAGE' ; * *--------------------------------------------------------------------- * TEST DARCY4 * CALCUL DARCY ISOTROPE TRANSITOIRE * AVEC TERME SOURCE * Résolution par une méthode d'éléments finis mixtes hybrides. * * Ce test permet de vérifier le bon fonctionnement des opérateurs * utilisés afin de résoudre les équations de DARCY transitoire par * une méthode d'éléments finis mixtes hybrides dans CASTEM2000. * * On effectue un calcul sur un domaine 2D maillé par des * quadrangles réguliers. Le phénomène étudié est monodimensionnel. * Les conditions aux limites sont de type trace de charge imposée * à gauche du domaine, flux imposé à droite et flux nul en haut et * en bas. * * La solution analytique en charge est un polynome de degré 2 * en x et en t, les variables x et t étant séparées. * La vitesse est donc linéaire en x et indépendante du temps. * Le terme source dépend de x et de t. * * Solution analytique : * H(x,t) = a * t**2 + b * t + c * x**2 + d * x + e * V(x,t) = - K * ( 2. * c * x + d ) * S(x,t) = 2 * a * f * t + b * f - 2 * c * K * * avec K : Conductivité hydraulique ................. 1. * a : Coefficient .............................. -1. * b : Coefficient .............................. 20. * c : Coefficient .............................. -0.01 * d : Coefficient .............................. 2. * e : Coefficient .............................. 10. * f : Coefficient d'emmagasinnement ............ 2. * * Remarque : Ne pas prendre t>40 car h devient négatif et on peut * planter le calcul d'erreur relative. * * Parametres de maillage : * L : longueur du domaine ....................... 100. * H : hauteur du domaine ........................ 1. * INUMX : Nombre d'éléments en x .................... 25 * INUMY : Nombre d'éléments en y .................... 1 car 1D * *---------------------------------------------------------------- * *------------------ * Options generales *------------------ * 'OPTI' 'DIME' 2 'ELEM' 'QUA4' ; 'OPTI' 'ISOV' 'SURF' ; 'OPTI' 'ECHO' 0 ; * * *========= * MAILLAGE *========= * * * d1------------d3 * | | Post-traitement via les lignes b1b3 et c1c3 * b1 c1------c3 b3 b1b3 : Extrémités de la droite joignant les faces * | | c1c3 : Extrémités de la droite joignant les centres * a1------------a3 * * *- Création des points supports des DROITES * L = 100.D0 ; H = 1.D0 ; X0 = 0.D0 ; X1 = L ; Y0 = 0.D0 ; Y1 = Y0 + H ; INUMX = 25 ; INUMY = 1 ; INUM1 = INUMX - 1 ; DL = L / INUMX ; DH = H / INUMY ; Y01 = Y0 + Y1 * 0.5D0 ; DX1 = X1 - X0 / INUMX / 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 ; * *- 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 ; DRMIF = B1 'DROI' INUMX B3 ; DRMIC = C1 'DROI' INUM1 C3 ; ELI0 = DH / 10.D0 ; * *- Creation maillage GEOMETRIQUE * PTOT1 = 'DALL' DRBAS DRGAU DRHAU DRDRO ; PTOT2 = 'ORIE' PTOT1 ; QFTOT = CHANGE PTOT2 QUAF ; QFGAU = CHANGE DRGAU QUAF ; QFDRO = CHANGE DRDRO QUAF ; ELIM ELI0 ( QFTOT ET QFGAU ET QFDRO ET DRMIF ET DRMIC) ; * * * * * ---------------- * = MODELISATION = * ---------------- MODHYB = MODE QFTOT 'DARCY' 'ISOTROPE' ; MODGAU = MODE QFGAU 'DARCY' 'ISOTROPE' ; MODDRO = MODE QFDRO 'DARCY' 'ISOTROPE' ; CHYB1 = 'DOMA' MODHYB 'SURFACE' ; CHYB2 = 'DOMA' MODHYB 'NORMALE' ; CETOT = DOMA MODHYB 'CENTRE' ; CEDRO = DOMA MODDRO 'CENTRE' ; CEGAU = DOMA MODGAU 'CENTRE' ; * *===================================== * CONDITIONS INITIALES : CAS PERMANENT *===================================== * * ----------------------- * = Solution analytique = * ----------------------- A = -1.D0 ; B = 20.D0 ; C = -0.01D0 ; D = 2.D0 ; E = 10.D0 ; F = 2.D0 ; VK = 1.D0 ; * TRACE DE CHARGE XX YY = 'COOR' (DOMA MODHYB 'FACE' ) ; PANAF0 = C * XX * XX + ( D * XX ) + E ; PANAF0 = 'EXCO' PANAF0 'SCAL' 'TH' ; PANAF0 = 'CHAN' 'ATTRIBUT' PANAF0 'NATURE' 'DISCRET' ; * CHARGE XXC YYC = 'COOR' CETOT ; PANAC0 = C * XXC * XXC + ( D * XXC ) + E ; PANAC0 = 'EXCO' PANAC0 'SCAL' 'H' ; PANAC0 = 'CHAN' 'ATTRIBUT' PANAC0 'NATURE' 'DISCRET' ; * VITESSE VXANA = -1.D0 * VK * ( 2.D0 * C * XXC + D ) ; VXANA = 'EXCO' VXANA 'SCAL' 'VX' ; VXANA = 'CHAN' VXANA 'ATTRIBUT' 'NATURE' 'DISCRET' ; VYANA = 'MANU' 'CHPO' CETOT 1 'VY' 0.D0 'NATURE' 'DISCRET' ; VANAC = VXANA 'ET' VYANA ; MOT1 = 'MOTS' 'VX' 'VY' ; * ---------------- * = Tenseur K = * ---------------- MATI1 = 'MATE' MODHYB 'K' VK ; * ------------------------- * = matrice Masse HYBride = * ------------------------- CND1A = 'MHYB' MODHYB MATI1 ; * ------------------------- * = MAtrice globale en TH = * ------------------------- HND1A = 'MATP' MODHYB CND1A ; * ----------------------- * = Second membre en TH = * ----------------------- AIRE = DL * DH ; SR0 = -2.D0 * VK * C * AIRE ; SRCVA0 = 'MANU' 'CHPO' CETOT 1 'SOUR' SR0 'NATURE' 'DISCRET' ; SMSOUR = 'SMTP' MODHYB CND1A SRCVA0 ; * ----------------- * = TH imposée = * ----------------- BBGAU = 'BLOQ' CEGAU 'TH' ; BBDRO = 'BLOQ' CEDRO 'TH' ; EEGAU = 'DEPI' BBGAU PANAF0 ; EEDRO = 'DEPI' BBDRO PANAF0 ; * ----------------- * = Résolution TH = * * ----------------- CCC1 = HND1A 'ET' BBDRO 'ET' BBGAU ; FFF1 = SMSOUR 'ET' EEDRO 'ET' EEGAU ; CHTER1 = 'RESO' CCC1 FFF1 ; TH1 = 'EXCO' CHTER1 'TH' 'TH' ; * ----------------------------- * = Post-traitement Q, H et V = * ----------------------------- PCEN1 = 'HYBP' MODHYB CND1A TH1 SRCVA0 ; QFACE1 = 'HDEB' MODHYB CND1A PCEN1 TH1 ; VCENT1 = 'HVIT' MODHYB QFACE1 ; * * * *=================================== * DONNEES SPECIFIQUES AU TRANSITOIRE *=================================== * * * TINI = 0.D0 ; TFINAL = 20.D0 ; TSUP = 1.5D0 * TFINAL ; DELTAT = 1.0D0 ; LICALC = 'PROG' 0.D0 'PAS' DELTAT TFINAL ; LISAUV = 'PROG' 0.D0 'PAS' 5.D0 TFINAL ; COEMM = 'MANU' 'CHML' (DOMA MODHYB 'MAILLAGE') 'CK' F ; * ------------------------------------ * = Nouvelles conditions aux limites = * ------------------------------------ * * Définition d'un listreel contenant les valeurs des temps de calcul * Ce LISTREEL sert à définir les chargements variables en temps * TDEF = PROG 0.D0 PAS 0.001D0 TSUP ; * * Valeur imposée à gauche = a*t*t + b*t + e * BLCGE = BBGAU ; VALGA0 = 'MANU' 'CHPO' CEGAU 1 'TH' E ; VALGA0 = 'DEPI' BBGAU VALGA0 ; VALI0 = 'CHAR' VALGA0 ('EVOL' 'MANU' ('PROG' 0.D0 TSUP) ('PROG' 1.D0 1.D0)) ; VALGA1 = 'MANU' 'CHPO' CEGAU 1 'TH' B ; VALGA1 = 'DEPI' BBGAU VALGA1 ; VALI1 = 'CHAR' VALGA1 ('EVOL' 'MANU' ('PROG' 0.D0 TSUP) ('PROG' 0.D0 TSUP)) ; VALGA2 = 'MANU' 'CHPO' CEGAU 1 'TH' A ; VALGA2 = 'DEPI' BBGAU VALGA2 ; TDEF2 = TDEF * TDEF ; VALI2 = 'CHAR' VALGA2 ('EVOL' 'MANU' TDEF TDEF2) ; TPGAU = VALI0 'ET' VALI1 'ET' VALI2 ; * * Flux imposé à droite = -K * ( 2. * c * x + d ) * dh * FLVAL = 2.D0 * C * X1 + D * (-1.D0) * VK * DH ; FLDRO = 'MANU' 'CHPO' CEDRO 1 'FLUX' FLVAL ; FLXDRO = 'CHAR' FLDRO ('EVOL' 'MANU' ('PROG' 0.D0 TSUP) ('PROG' 1.D0 1.D0) ) ; * * Définition de la source moyenne par élément sur un pas de temps * s(t) = (2 * a * f * tmoy + b * f - 2 * c * K) * aire * avec aire = aire des éléments * tmoy = t(n+1) - deltat/2 = t(n) + deltat/2 * SRCE01 = B * F - ( 2.D0 * VK * C ) ; SRCE02 = (-1.D0) * F * A * DELTAT ; SR0 = SRCE01 + SRCE02 * AIRE ; SRCV10 = 'MANU' 'CHPO' CETOT 1 'SOUR' SR0 'NATURE' 'DISCRET' ; SRCVA0 = SRCV10 ; NTDEF = 'DIME' TDEF ; LIST1 = 'PROG' NTDEF * 1.D0 ; CHSRC0 = 'CHAR' SRCVA0 ('EVOL' 'MANU' TDEF LIST1) ; PRSR1 = 2.D0 * A * F * AIRE * TDEF ; SRCVA1 = 'MANU' 'CHPO' CETOT 1 'SOUR' 1.D0 'NATURE' 'DISCRET' ; CHSRC1 = 'CHAR' SRCVA1 ('EVOL' 'MANU' TDEF PRSR1 ) ; CHSRCE = CHSRC0 'ET' CHSRC1 ; * * --------------------------- * = Table DARCY_TRANSITOIRE = * --------------------------- TRANS2 = 'TABLE' ; TRANS2. 'TEMPS' = 'TABLE' ; TRANS2. 'TRACE_CHARGE' = 'TABLE' ; TRANS2. 'CHARGE' = 'TABLE' ; TRANS2. 'FLUX' = 'TABLE' ; * TRANS2.'SOUSTYPE' = 'DARCY_TRANSITOIRE' ; TRANS2.'MODELE' = MODHYB ; TRANS2.'CARACTERISTIQUES' = MATI1 ; TRANS2.'EMMAGASINEMENT' = COEMM ; TRANS2.'DOMAINE' = DOMA MODHYB 'TABLE' ; * TRANS2. 'TEMPS' . 0 = TINI ; TRANS2. 'TRACE_CHARGE' . 0 = TH1 ; TRANS2. 'CHARGE' . 0 = PCEN1 ; TRANS2. 'FLUX' . 0 = QFACE1 ; * TRANS2.'BLOCAGE' = BLCGE ; TRANS2.'TRACE_IMPOSE' = TPGAU ; TRANS2.'FLUX_IMPOSE' = FLXDRO ; TRANS2.'SOURCE' = CHSRCE ; * TRANS2.'THETA' = 0.5D0 ; TRANS2.'TEMPS_CALCULES' = LICALC ; TRANS2.'TEMPS_SAUVES' = LISAUV ; * *======================= * Resolution transitoire *======================= * DARCYTRA TRANS2 ; * * *================= * POST-TRAITEMENT *================= * -------------------------- * = Calcul ERREUR ET TRACE = * -------------------------- *-------------------------------------------------------------------- * Dans chaque cas on trace (si GRAPH<>N) *-------------------------------------------------------------------- * Les traces de charge le long de DRMIF (ligne b1b3) * L'erreur relative et l'erreur absolue en Trace de charge * La charge le long de DRMIC (ligne c1c3) * L'erreur relative et l'erreur absolue sur la charge * La Vitesse aux centres des elements (calculé et analytique) * L'erreur relative et l'erreur absolue sur la vitesse au centre *-------------------------------------------------------------------- * * VDVD = 'PSCA' VANAC VANAC MOT1 MOT1 ; VCANA = 'VECT' VANAC 0.01D0 VX VY 'ROUG' ; *-------------------------------------------------- * * VISURESU * Calcul des solutions analytiques et des erreurs * Visualisation sur DRMIF et DRMIC des courbes y=f(x) à t fixé * * ISOR1 : Liste des indices de la table TRANS2.'TEMPS' * NTSOR : Nombre de pas de temps sauvegardés * IOK : Booleen indiquant si le critère de précision * demandé est vérifiée (Gestion des erreurs) * IRESU : Compteur interne permettant de récuperer les * indices de table via l'index ISOR1. * ISOR1 = 'INDEX' ( TRANS2 . 'TEMPS') ; NTSOR = 'DIME' ISOR1 ; IOK = FAUX ; IRESU = 0 ; *-------------------------------------------------- 'REPETER' VISURESU NTSOR ; * IRESU = IRESU + 1 ; INDI1 = ISOR1 . IRESU ; TTRA = TRANS2 . 'TEMPS' . INDI1 ; * Erreurs en trace de charge PANAF1 = A*TTRA*TTRA + (B*TTRA) + (C*XX*XX) + (D*XX) + E ; PANAD1 = 'REDU' PANAF1 DRMIF ; CALCF1 = TRANS2 . 'TRACE_CHARGE' . INDI1 ; CALCF1 = 'EXCO' CALCF1 'TH' 'SCAL' ; CALCD1 = 'REDU' CALCF1 DRMIF ; ETPAB1 = CALCD1 - PANAD1 ; ETPAB1 = 'ABS' ETPAB1 ; ETPRE1 = ETPAB1 / PANAD1 ; ETPRE1 = 'ABS' ETPRE1 ; * Erreurs en charge K1 = C * AIRE * AIRE /12.D0 + E ; PANAC1 = A*TTRA*TTRA + (B*TTRA) + (C*XXC*XXC) + (D*XXC) + K1 ; CALC1 = TRANS2 . 'CHARGE' . INDI1 ; CALC1 = 'EXCO' CALC1 'H' 'SCAL' ; EPAB1 = CALC1 - PANAC1 ; EPAB1 = 'ABS' EPAB1 ; EPRE1 = EPAB1 / PANAC1 ; EPRE1 = 'ABS' EPRE1 ; * Erreurs en vitesse EFL1 = TRANS2 . 'FLUX' . INDI1 ; EVI1 = 'HVIT' MODHYB EFL1 ; VCAL1 = 'VECT' EVI1 0.01D0 'VX' 'VY' 'VERT' ; ERRV1 = EVI1 - VANAC ; EVAB1 = 'PSCA' ERRV1 ERRV1 MOT1 MOT1 ; EVRE1 = 'ABS' ( EVAB1 / VDVD ) ; EVAB1 = EVAB1 '**' 0.5D0 ; EVRE1 = EVRE1 '**' 0.5D0 ; * Evolutions des solutions analytiques AV1 = 'EVOL' 'ROUG' 'CHPO' PANAD1 'SCAL' DRMIF ; AV3 = 'EVOL' 'ROUG' 'CHPO' PANAC1 'SCAL' DRMIC ; * Tracés 'SI' ('NEG' GRAPH 'N' ) ; LTITRE = 'TEXTE' 'darcy4 Th temps ' TTRA ; 'TITR' LTITRE ; EV1 = 'EVOL' 'VERT' 'CHPO' CALCD1 'SCAL' DRMIF ; 'DESS' (EV1 ET AV1) ; LTITRE = 'TEXTE' 'darcy4 erreur relative Th temps ' TTRA ; 'TITR' LTITRE ; EV2 = 'EVOL' 'VERT' 'CHPO' ETPRE1 'SCAL' DRMIF ; 'DESS' EV2 ; LTITRE = 'TEXTE' 'darcy4 erreur absolue Th temps ' TTRA ; 'TITR' LTITRE ; EV3 = 'EVOL' 'VERT' 'CHPO' ETPAB1 'SCAL' DRMIF ; 'DESS' EV3 ; LTITRE = 'TEXTE' 'darcy4 h temps ' TTRA ; 'TITR' LTITRE ; EV4 = 'EVOL' 'VERT' 'CHPO' CALC1 'SCAL' DRMIC ; 'DESS' (EV4 ET AV3) ; LTITRE = 'TEXTE' 'darcy4 erreur relative h temps' TTRA ; 'TITR' LTITRE ; EV5 = 'EVOL' 'VERT' 'CHPO' EPRE1 'SCAL' DRMIC ; 'DESS' EV5 ; LTITRE = 'TEXTE' 'darcy4 erreur absolue h temps' TTRA ; 'TITR' LTITRE ; EV6 = 'EVOL' 'VERT' 'CHPO' EPAB1 'SCAL' DRMIC ; 'DESS' EV6 ; LTITRE = 'TEXTE' 'darcy4 erreur relat V temps ' TTRA ; 'TITR' LTITRE ; EV7 = 'EVOL' 'VERT' 'CHPO' EVRE1 'SCAL' DRMIC ; 'DESS' EV7 ; LTITRE = 'TEXTE' 'darcy4 erreur absolue V temps ' TTRA ; 'TITR' LTITRE ; EV8 = 'EVOL' 'VERT' 'CHPO' EVAB1 'SCAL' DRMIC ; 'DESS' EV8 ; 'FINSI' ; * *- Gestion ERREUR * MAXTP1 = 'MAXI' ETPRE1 ; MAXP1 = 'MAXI' EPRE1 ; MAXVC1 = 'MAXI' EVRE1 ; LOGRES = IRESU 'EGA' 1 ; 'SI' (LOGRES) ; 'SAUT' 'PAGE' ; 'SAUT' 2 'LIGNE' ; 'MESS' ' ERREURS RELATIVES '; 'SAUT' 1 'LIGNE' ; 'MESS' ' temps TH H Vc'; 'SAUT' 1 'LIGNE' ; 'FINSI' ; 'MESS' ttra ' ' maxtp1 ' ' maxp1 ' ' maxvc1 ; EPS0 = 1.E-04 ; LOG1 = MAXTP1 > EPS0 ; LOG2 = MAXP1 > EPS0 ; LOG3 = MAXVC1 > EPS0 ; L0 = LOG1 'OU' LOG2 'OU' LOG3 ; 'SI' ( L0 ) ; IOK = VRAI ; 'FINSI' ; 'FIN' VISURESU ; * ********************************* *- En cas de tracé uniquement ********************************* * 'SI' ('NEG' GRAPH 'N' ) ; *-------------------------------------------------- * * CONSTAB * (préparation pour VISPTS) * Transformation des champoints TH=f(t) en listreels * * TABEVOL : Table contenant TH=f(t) sous forme de LISTREEL * IRESU : Compteur interne permettant de récuperer les * indices de table via l'index ISOR1. * TABEVOL = 'TABLE' ; IRESU = 0 ; * *-------------------------------------------------- 'REPETER' CONSTAB NTSOR ; IRESU = IRESU + 1 ; INDI1 = ISOR1 . IRESU ; TTRA = TRANS2 . 'TEMPS' . INDI1 ; ETP1 = TRANS2 . 'TRACE_CHARGE' . INDI1 ; ETP1 = 'EXCO' ETP1 'TH' 'SCAL' ; EV1 = 'EVOL' 'VERT' 'CHPO' ETP1 'SCAL' DRMIF ; ORD1 = 'EXTR' EV1 'ORDO' ; TABEVOL . IRESU = ORD1 ; 'SI' ( IRESU 'EGA' 1) ; LITRAC = 'PROG' TTRA ; 'SINO' ; LITRAC = LITRAC 'ET' ('PROG' TTRA) ; 'FINSI' ; 'FIN' CONSTAB ; *-------------------------------------------------- * * VISPTS * Visualisation des courbes TH=f(t) pour un x fixé * * NTTRA : Nombre de valeurs en temps à considérer * LISTX : Listreel des abscisses * NXTOT : Nombre de points sur la ligne DRMIF * IX : Position dans LISTX du premier point à dépouiller * PASX : Incrément pour le pas suivant à dépouiller * * NTTRA = 'DIME' LITRAC ; LISTX = 'EXTR' EV1 'ABSC' ; NXTOT = 'NBNO' DRMIF ; IX = 5 ; PASX = 10 ; * *-------------------------------------------------- 'REPETER' VISTPS ; * Solution analytique VALX = 'EXTR' LISTX IX ; PPROG0 = 'PROG' NTTRA * (C*VALX*VALX + (D*VALX) + E ) ; PPROG1 = A * LITRAC * LITRAC + (B * LITRAC) ; PEVOL1 = PPROG0 + PPROG1 ; EVOL1 = 'EVOL' 'ROUG' 'MANU' 'Temps' LITRAC 'Th' PEVOL1 ; * Solution calculée ICRB = 0 ; 'REPETER' CONSTR NTSOR ; ICRB = ICRB + 1 ; VAL1 = 'EXTR' TABEVOL.ICRB IX ; 'SI' ( ICRB 'EGA' 1 ) ; PROG1 = 'PROG' VAL1 ; 'SINO' ; PROG1 = PROG1 'ET' ('PROG' VAL1 ) ; 'FINSI' ; 'FIN' CONSTR ; EVOL2 = 'EVOL' 'VERT' 'MANU' 'Temps' LITRAC 'Th' PROG1 ; * Tracés LTITRE = 'TEXTE' 'darcy4 Th en x=' VALX ; 'TITR' LTITRE ; 'DESS' ( EVOL1 ET EVOL2 ) ; LTITRE = 'TEXTE' 'erreur relative Th x=' VALX ; 'TITR' LTITRE ; ERRTP1 = PEVOL1 - PROG1 / PEVOL1 ; EVOL3 = 'EVOL' 'VERT' 'MANU' 'Temps' LITRAC 'Th' ERRTP1 ; 'DESS' EVOL3 ; * Test de sortie de boucle 'DETR' PROG1 ; IX = IX + PASX ; 'SI' (IX > NXTOT ) ; 'QUITTER' VISTPS ; 'FINSI' ; 'FIN' VISTPS ; * 'FINSI' ; * **************************** * FIN des tracés **************************** * 'SAUT' 4 'LIGNE' ; * 'SI' ( IOK ) ; 'ERRE' 5 ; 'SINO' ; 'ERRE' 0 ; 'FINSI' ; * 'FIN' ;