* fichier : darcy5.dgibi
*
******************** CAS TEST : darcy5.dgibi **********************
*
GRAPH = FAUX ;
*GRAPH = VRAI ;
'SAUT' 'PAGE' ;
*
*-------------------------------------------------------------------
* TEST DARCY5
*
* Validation de la résolution des équations de DARCY formulées en
* (VITESSE - PRESSION) par comparaison avec
* la résolution des équations de DARCY formulées en (VITESSE - CHARGE)
*
*
* (Résolution par une méthode d'éléments finis mixtes hybrides)
*
* Le maillage est composés de triangles et de quadrangles et contient
* des lignes de maillages courbes
*
* la massif est supposé isotrope, indéformable et saturé
*
* Pour permettre la comparaison, la densité de l'eau reste constante
*
*-------------------------------------------------------------------
*
OPTI ECHO 0 ;
*
'SAUT' 'PAGE' ;
*
*- Options générales de calcul
*
TITRE
'EFMY DARCY FORMULATION (VITESSE - PRESSION) 2D : darcy5.dgibi' ;
OPTI DIME 2 ELEM QUA4 ;
OPTI ISOV SURF ;
*
*- Données physiques
*
*- Densité de l'eau
RHO0 = 1000. ;
*
*- Pression athmosphérique
P0 = 1.013E5 ;
*
*- Composante du vecteur gravité
Gx0 = 0. ;
Gy0 = -9.81 ;
*
*- Perméabilité intrinsèque
LA0 = 1.e-9 ;
*
*======================================================================
*
*- Construction du maillage
*
*======================================================================
*
*- Discrétisation spatiale
*
ENX = 15 ;
ENY0 = 5 ;
ENY1 = 5 ;
ENY2 = 5 ;
*
*- Création des points
*
A0 = 0.d0 -2.d3 ; B0 = 5.d3 -2.d3 ;
*
A1 = 0.d0 -1.d3 ; B1 = 5.d3 -1.d3 ; C1 = 2.d3 -1.25d3 ;
*
A2 = 0.d0 0.0d3 ; B2 = 5.d3 0.0d3 ; C2 = 2.d3 0.25d3 ;
*
A3 = 0.d0 1.5d3 ; B3 = 5.d3 1.d3 ; C3 = 2.5d3 1.5d3 ;
*
*- Création des droites
*
AB0 = DROI ENX A0 B0 ;
AB1 = A1 cer3 C1 ENX B1 ;
AB2 = A2 cer3 C2 ENX B2 ;
AB3 = A3 cer3 C3 ENX B3 ;
*
*- Creation des SURFACES
*
NIVEAU0 = AB0 regle ENY0 AB1 ;
OPTI ELEM TRI3 ;
NIVEAU1 = AB1 regle ENY1 AB2 ;
OPTI ELEM QUA4 ;
NIVEAU2 = AB2 regle ENY2 AB3 ;
MASSIF0 = NIVEAU0 ET NIVEAU1 ET NIVEAU2 ;
HAUT = INVE AB3 ;
MASSIF0 = ORIENT MASSIF0 ;
QMASSIF = CHANGE MASSIF0 QUAF ;
QFHAUT = CHANGE HAUT QUAF ;
ELIM 0.01 (QMASSIF ET QFHAUT ) ;
*
*- Domaines
*
*HYTOT = DOMA MASSIF0 0.01 ;
*CHYB1 = DOMA HYTOT 'SURFACE' ;
*CHYB2 = DOMA HYTOT 'NORMALE' ;
*MCHYB = DOMA HYTOT 'ORIENTAT' ;
*HYHAUT = DOMA HAUT INCL HYTOT 0.01 ;
*
*- Modélisation
*
*MODHYB = MODL HYTOT DARCY ISOTROPE HYQ4 HYT3 ;
MODHYB = MODL QMASSIF DARCY ISOTROPE ;
MODHAU = MODL QFHAUT DARCY ISOTROPE ;
CHYB1 = DOMA MODHYB 'SURFACE' ;
CHYB2 = DOMA MODHYB 'NORMALE' ;
CEHAUT= DOMA MODHAU 'CENTRE' ;
CETOT= DOMA MODHYB 'CENTRE' ;
*
*- perméabilité intrinsèque
*
LMME = MANU CHPO CETOT 1 'K' la0 'NATURE' 'DIFFU' ;
*
*======================================================================
*
*- Résolution en en (VITESSE ; PRESSION)
*
*======================================================================
*
LMMA = KCHA MODHYB 'CHAM' LMME ;
MATI_P = MATERIAU MODHYB 'K' LMMA ;
*
**** Masse HYBride
CND1A_P = MHYBR MODHYB MATI_P ;
*
**** Masse FORce
M_P = MHYBR MODHYB 'MASSE' ;
*
**** MAtrice en Trace de pression TP
HND1A_P = MATP MODHYB CND1A_P ;
*
*- Conditions aux limites
*
BHAUT_P = BLOQUE CEHAUT 'TH' ;
*
*- TP imposee
*
PIMP = MANU CHPO CEHAUT 1 'TH' P0 'NATURE' 'DISCRET' ;
EHAUT_P = DEPI BHAUT_P PIMP ;
*
*- contribution des forces de volume au second membre
*
RHO = MANU CHPO CETOT 1 'SCAL' RHO0 'NATURE' 'DISCRET' ;
RCH = MANU CHPO CETOT 2 'FX' gx0 'FY' gy0 'NATURE' 'DISCRET' ;
RCH = RHO * RCH ;
GRAV2TP = SQTP MODHYB HND1A_P M_P RCH ;
*
*- Assemblage matrice et second membre
*
CCC1_P = HND1A_P ET BHAUT_P ;
* second membre du système matriciel en trace de charge
FFF1_P = GRAV2TP et EHAUT_P ;
*
*- Resolution en trace de pression
*
CHTER1_P = RESOUDRE CCC1_P FFF1_P ;
CHTER1_P = exco CHTER1_P 'TH' 'TH';
*
*- Calcul des pressions élémentaires
*
PCEN1 = HYBP MODHYB CND1A_P CHTER1_P M_P RCH ;
*
*- Calcul des débits
*
QFACE1_P = HDEBI MODHYB CND1A_P PCEN1 CHTER1_P M_P RCH;
VCENT1_P = HVIT MODHYB QFACE1_P ;
*VV = VECT VCENT1_P 0.5e9 'VX' 'VY' VERT ;
*titre 'VITESSES AUX CENTRES DES ELEMENTS : ECHL=1e9:2' ;
*TRAC VV CETOT ;
*
*=======================================================================
*
*- Résolution en (VITESSE ; CHARGE)
*
*=======================================================================
*
*
*- conductivité hydraulique
*
LMME = LMME * (ABS Gy0) * RHO0 ;
LMMA = KCHA MODHYB 'CHAM' LMME ;
MATI_H = MATERIAU MODHYB 'K' LMMA ;
*
**** Masse HYBride
CND1A_H = MHYBR MODHYB MATI_H ;
*
**** MAtrice en Trace de charge TH
HND1A_H = MATP MODHYB CND1A_H ;
*
*- Conditions aux limites
*
BHAUT_H = BLOQUE CEHAUT 'TH' ;
*
*- TH imposee
*
PIMP = NOMC 'TH' (CEHAUT COOR 2) ;
EHAUT_H = DEPI BHAUT_H PIMP ;
*
*- Assemblage matrice et second membre
*
CCC1_H = HND1A_H ET BHAUT_H ;
* second membre du système matriciel en trace de charge
FFF1_H = EHAUT_H ;
*
*- Resolution en trace de charge
*
CHTER1_H = RESOUDRE CCC1_H FFF1_H ;
CHTER1_H = exco CHTER1_H 'TH' 'TH';
*
*- Calcul de la charge
*
TCEN1 = HYBP MODHYB CND1A_H CHTER1_H ;
*
*- Calcul de la vitesse
*
QFACE1_H = HDEBI MODHYB CND1A_H TCEN1 CHTER1_H ;
VCENT1_H = HVIT MODHYB QFACE1_H ;
*VV = VECT VCENT1_H 0.5e9 'VX' 'VY' VERT ;
*titre 'CALCUL EN CHARGE : VITESSES AUX CENTRES DES ELTS : ECHL=1e9:2';
*TRAC VV CETOT ;
*
*=======================================================================
*
*- Comparaison des deux résolutions
*
*=======================================================================
*
*
*- comparaison des traces de charges
*
TP2 = (NOMC CHTER1_P SCAL) - P0 / (rho0 * gy0) ;
TP2 = ((DOMA MODHYB 'FACE') COOR 2) - TP2 ;
ER2 = (NOMC CHTER1_H SCAL)/TP2 - 1. ;
maxtp1 = ER2 ABS MAXI ;
*
*
*- comparaison des charges
*
PCEN2 = (NOMC PCEN1 SCAL) - P0 / (rho0 * gy0) ;
PCEN2 = CETOT COOR 2 - PCEN2 ;
ER2 = (NOMC TCEN1 SCAL)/PCEN2 - 1. ;
maxp1 = ER2 ABS MAXI ;
'SI' GRAPH ;
TITRE 'VARIATIONS RELATIVES DES CHARGES DH/H' ;
TRAC (KCHA MODHYB 'CHAM' er2) MODHYB (CONTOUR
(DOMA MODHYB 'MAILLAGE')) ;
'FINSI' ;
*
*- comparaison des vitesse
*
QFACE1_P = NOMC SCAL QFACE1_P ;
VFACE1_P = QFACE1_P * CHYB2 / CHYB1 ;
QFACE1_H = NOMC SCAL QFACE1_H ;
VFACE1_H = QFACE1_H * CHYB2 / CHYB1 ;
maxvf1 = (kops QFACE1_P '-' QFACE1_H) / (maxi QFACE1_H) ABS MAXI ;
MOT1 = 'MOTS' 'VX' 'VY' ;
VDVD = 'PSCA' VCENT1_H VCENT1_H MOT1 MOT1 ;
VD1 = VCENT1_P - VCENT1_H ;
VC1 = 'PSCA' VD1 VD1 MOT1 MOT1 ;
SDC1 = 'ABS' ( VC1 / VDVD ) ;
SDC1 = SDC1 '**' 0.5D0 ;
maxvc1 = MAXI SDC1 ;
'SI' GRAPH ;
TITRE ' VARIATIONS RELATIVES DES VITESSES |DV|/|V|' ;
TRAC ('KCHA' MODHYB 'CHAM' SDC1) MODHYB (CONTOUR
(DOMA MODHYB 'MAILLAGE')) ;
'FINSI' ;
*
'SAUT' 'PAGE' ;
'SAUT' 2 'LIGNE' ;
'MESS' ' VARIATIONS RELATIVES ' ;
'SAUT' 1 'LIGNE' ;
'MESS' ' TH H Vf Vc' ;
'SAUT' 1 'LIGNE' ;
'MESS' ' ' maxtp1 ' ' maxp1 ' ' maxvf1 ' ' maxvc1 ;
'SAUT' 2 'LIGNE' ;
*
EPS0 = 1.E-7 ;
LOG1 = MAXTP1 > EPS0 ;
LOG2 = MAXP1 > EPS0 ;
LOG3 = MAXVF1 > EPS0 ;
LOG4 = MAXVC1 > EPS0 ;
L0 = LOG1 'OU' LOG2 'OU' LOG3 'OU' LOG4 ;
*
'SI' ( L0 ) ;
'ERRE' 5 ;
'SINO' ;
'ERRE' 0 ;
'FINS' ;
FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer