* fichier : darcy1.dgibi
***************** CAS TEST : darcy1.dgibi *************************
*
GRAPH = 'N' ;
'SAUT' 'PAGE' ;
*
*-------------------------------------------------------------------
* TEST DARCY1
*
* CALCUL DARCY ISOTROPE
* 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 pour la résolution des équations de DARCY par une méthode
* d'éléments finis mixtes hybrides dans CASTEM2000.
*
* On effectue trois calculs sur un domaine carré, maillé par des
* quadrangles réguliers. Le phénomène étudié est monodimensionnel.
* Les conditions aux limites varient suivant le cas considéré :
* on impose le flux ou la charge sur les cotés du domaine.
*
* La solution analytique en charge étant un polynome de degré un
* en x, la vitesse est constante et horizontale.
* Si on note L la longueur du domaine, HD (resp. HG) la charge à
* droite du domaine (resp. à gauche) et dH=HD-HG on a :
* h(x) = (x-L) * dH/L + HD
* Et par suite
* V(x) = (-K*dH/L ; 0.)
*
* On s'attend à une précision de l'ordre de la précision machine.
*
*-------------------------------------------------------------------
*
'SAUT' 'PAGE' ;
*
*- Options générales de calcul.
*
'TITR' 'EFMH DARCY ISOTROPE 1D : darcy1.dgibi' ;
* Dimension 2, éléments à générer quadrangulaires ou triangulaires :
'OPTI' 'DIME' 2 'ELEM' 'QUA4' ;
* Pas d'affichage à l'écran des lignes de commandes lues :
'OPTI' 'ECHO' 1 ;
* ------------
* = MAILLAGE =
* ------------
*
*- Création des points supports des DROITES
*
* Dimensions :
L = 1.D0 ;
HA = 1.D0 ;
* Coordonnées :
XG = 0.D0 ;
XD = XG + L ;
YG = 0.D0 ;
YD = YG + HA ;
* Points
A1 = XG YG ;
A3 = XD YG ;
D1 = XG YD ;
D3 = XD YD ;
*
*- Création des DROITES frontières
*
* Discrétisation dans les deux directions :
INUMX = 9 ;
INUMY = 4 ;
* droites maillées :
DRBAS = A3 'DROI' INUMX A1 ;
DRGAU = A1 'DROI' INUMY D1 ;
DRHAU = D1 'DROI' INUMX D3 ;
DRDRO = D3 'DROI' INUMY A3 ;
*
*- Création maillage GEOMETRIQUE
*
* on va faire des parallèlépipèdes réguliers grâce à DALLER
* Les points définissant les quadrangles sont appellés 'SOMMET'
CARRE = 'DALLER' DRBAS DRGAU DRHAU DRDRO 'PLAN' ;
*
*- Création maillage HYBRIDE y compris sous-objets (cond. limites)
* ce sont les maillages avec les points sommets + les points faces +
* les points centres.
* on crée celui du maillage total
DOMHYB = 'CHANGER' CARRE 'QUAF' ;
* et ceux des droites sur lesquelles on voudra mettre des conditions
* aux limites (en fait les quatre côtés, ici) :
DOMGAU = 'CHANGER' DRGAU 'QUAF' ;
DOMDRO = 'CHANGER' DRDRO 'QUAF' ;
DOMHAU = 'CHANGER' DRHAU 'QUAF' ;
DOMBAS = 'CHANGER' DRBAS 'QUAF' ;
* On a créé à chaque fois les points faces et centres.
* Comme CARRE, DRGAU, DRDRO, DRHAU et DRBAS avaient des segments en
* commun, les opérations précédentes ont créé beaucoup de points en
* double. Il faut les éliminer pour qu'on parle bien de la même chose.
* critère d'élimination (inférieur à la plus petite dimension d'arête
ELI0 = L / INUMX / 10.D0 ;
* élimination proprement dite :
'ELIMINATION' ( DOMHYB ET DOMGAU ET DOMDRO ET DOMHAU ET DOMBAS) ELI0 ;
* ----------------
* = MODELISATION =
* ----------------
* Définition du modèle physique attaché au domaine d'étude :
* Il servira d'une part pour la résolution,
* d'autre part pour extraire les points faces, centres, et autres
* informations utiles
MODHYB = 'MODELE' DOMHYB 'DARCY' 'ISOTROPE' ;
* On voudra aussi disposer des points faces des frontières,
* c'est-à-dire des points centres des lignes constituant ces
* frontières. Pour cela, on attache à ces lignes une modélisation,
* inutile en soi, mais qui dit implicitement (comme c'est un modèle
* 'DARCY'), que l'on travaille en éléments finis mixte hybrides,
* et donc qu'on aura besoin, à un moment ou à un autre, des points
* centres, faces, et autres :
MODGAU = 'MODELE' DOMGAU 'DARCY' 'ISOTROPE' ;
MODDRO = 'MODELE' DOMDRO 'DARCY' 'ISOTROPE' ;
MODHAU = 'MODELE' DOMHAU 'DARCY' 'ISOTROPE' ;
MODBAS = 'MODELE' DOMBAS 'DARCY' 'ISOTROPE' ;
* On peut maintenant extraire via ces modèles, les points centre des
* maillages 'QUAF' des frontières, qui sont en fait les points faces
* en limite du domaine :
CEDRO = 'DOMA' MODDRO 'CENTRE' ;
CEHAU = 'DOMA' MODHAU 'CENTRE' ;
CEBAS = 'DOMA' MODBAS 'CENTRE' ;
CEGAU = 'DOMA' MODGAU 'CENTRE' ;
* Les autres petites informations, sur le domaine entier :
* Maillage des points centres :
HYCEN = 'DOMA' MODHYB 'CENTRE' ;
* Maillage des points faces :
HYFAC = 'DOMA' MODHYB 'FACE' ;
* Longueur des arêtes (on est en 2D) :
HYSUR = 'DOMA' MODHYB 'SURFACE' ;
* Surface des mailles :
HYVOL = 'DOMA' MODHYB 'VOLUME' ;
* Vecteurs normaux aux faces, champ-points à autant de composantes que
* la dimension de l'espace : 'UX' , 'UY' (, 'UZ').
HYNOR = 'DOMA' MODHYB 'NORMALE' ;
*
* ------------------------------------
* = Tenseur de perméabilité isotrope =
* ------------------------------------
* La grandeur physique correspondant au laplacien. C'est ici la
* perméabilité hydraulique.
VK = 0.75D0 ;
* On la stoque dans un objet de type matériau
MATI3 = 'MATERIAU' MODHYB 'K' VK ;
* -----------------------
* = Solution analytique =
* -----------------------
* Elle va servir à deux choses : définir les conditions aux limites,
* et vérifier la précision de la résolution.
*
* Trace de charge aux faces
* La charge étant linéaire, Th moyen est au centre de la face
*
XX YY = 'COOR' (DOMA MODHYB 'FACE') ;
HG = 100.D0 ;
HD = 20.D0 ;
DH = HD - HG ;
PANAF = (XX - L) * DH / L + HD ;
*
* Charge aux centres des éléments
* La charge étant linéaire, la moyenne est au centre de l'élément
*
XXC YYC = 'COOR' (DOMA MODHYB 'CENTRE' ) ;
PANAC = (XXC - L) * DH / L + HD ;
*
* Composantes scalaires de la vitesse
*
VXANA = ( -1.D0 ) * VK * DH / L ;
VYANA = 0.D0 ;
*
* Vitesse aux faces et aux centres
*
VANAF = 'MANU' 'CHPO' (DOMA MODHYB 'FACE') 2 'VX' VXANA 'VY' VYANA ;
VANAC = 'MANU' 'CHPO' (DOMA MODHYB 'CENTRE') 2 'VX' VXANA 'VY' VYANA ;
*
* -------------------------
* = matrice Masse HYBride =
* -------------------------
* C'est l'inverse du tenseur de perméabilité. Il ne faut donc jamais
* avoir mis une valeur de perméabilité nulle.
CND1A = 'MHYB' MODHYB MATI3 ;
* -------------------------
* = MAtrice globale en TH =
* -------------------------
* La matrice correspondant à div (K 'GRAD').
* C'est un objet de type RIGIDITE.
HND1A = 'MATP' MODHYB CND1A ;
* ----------------
* = TH imposée =
* ----------------
* Les conditions aux limites de Dirichlet portent sur la charge.
* On est en EFMH, donc ce sont les valeurs aux faces qu'on utilise.
* On les appelle trace de charge. C'est un champ-point appuyé aux
* face, de nom de composante 'TH'.
* Contribution au membre de gauche pour chaque sous-partie de la frontière
BBGAU = 'BLOQ' CEGAU 'TH' ;
BBDRO = 'BLOQ' CEDRO 'TH' ;
BBHAU = 'BLOQ' CEHAU 'TH' ;
BBBAS = 'BLOQ' CEBAS 'TH' ;
* Contribution au membre de droite pour chaque sous-partie de la frontière
TTIMP = 'REDU' PANAF CEGAU ;
TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH' ;
EEGAU = 'DEPI' BBGAU TTIM2 ;
*
TTIMP = 'REDU' PANAF CEDRO ;
TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH' ;
EEDRO = 'DEPI' BBDRO TTIM2 ;
*
TTIMP = 'REDU' PANAF CEBAS ;
TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH' ;
EEBAS = 'DEPI' BBBAS TTIM2 ;
*
TTIMP = 'REDU' PANAF CEHAU ;
TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH' ;
EEHAU = 'DEPI' BBHAU TTIM2 ;
* ----------------
* = Flux imposé =
* ----------------
* Conditions aux limites de Neuman.
* En discrétisation EFMH comme ici, ce sont des débits intégrés aux
* faces. C'est un champ-point de composante 'FLUX'
* Il faut tenir compte de l'orientation des normales aux faces,
* toujours sortantes aux frontières du domaine. Une valeur de débit
* positive correspond donc à un flux sortant.
* flux rentrant à gauche - (- K 'GRAD' h) :
FG = VK * DH / L ;
* champ de débit imposé à gauche, constant sur les faces de la
* frontière de gauche (il faut multiplier le flux par la longueur des
* arêtes pour obtenir le débit intégré) :
FGAU = 'MANU' 'CHPO' CEGAU 1 'FLUX' FG ;
LGAU = 'REDU' HYSUR CEGAU ;
FLGAU = FGAU * LGAU ;
* on en fait un champ de nature discret (diffus conviendrait tout aussi bien) :
* pour que sa nature ne reste pas 'INDETERMINE'. Cela permettra de le
* concaténer avec d'autres champs par la suite.
FLGAU = 'CHAN' FLGAU 'ATTRIBUT' 'NATURE' 'DISCRET' ;
*
* idem pour le flux sortant à droite :
FD = ( -1.D0 ) * FG ;
FDRO = 'MANU' 'CHPO' CEDRO 1 'FLUX' FD ;
LDRO = 'REDU' HYSUR CEDRO ;
FLDRO = FDRO * LDRO ;
FLDRO = 'CHAN' FLDRO 'ATTRIBUT' 'NATURE' 'DISCRET' ;
* ----------------
* = Assemblage =
* ----------------
* Création des membres de gauche et de droite de l'équation "A TH = B"
* pour les trois jeux de conditions aux limites choisis.
* CCC# correspond à A, FFF# correspond à B.
* Les CL de Neuman n'ont pas de contribution au terme de gauche.
* CL de Dirichlet partout :
CCC1 = HND1A 'ET' BBHAU 'ET' BBBAS 'ET' BBDRO 'ET' BBGAU ;
FFF1 = EEHAU 'ET' EEBAS 'ET' EEDRO 'ET' EEGAU ;
* CL de Dirichlet à gauche et à droite,
* et flux nul imposé en haut et en bas (c'est ce qui se passe par
* défaut, donc il suffit de ne rien mettre) :
CCC2 = HND1A 'ET' BBDRO 'ET' BBGAU ;
FFF2 = EEDRO 'ET' EEGAU ;
* CL de Dirichlet en hau et en bas,
* et flux imposés à gauche et à droite :
CCC3 = HND1A 'ET' BBHAU 'ET' BBBAS ;
FFF3 = EEHAU 'ET' EEBAS 'ET' FLDRO 'ET' FLGAU ;
* -----------------
* = Résolution TH =
* -----------------
* On résoud les trois systèmes d'équations (CL comprises)
* On obtient en sortie la trace de charge, valeur moyenne de la charge
* sur chaque arête.
CHTER1 = 'RESO' CCC1 FFF1 ;
CHTER2 = 'RESO' CCC2 FFF2 ;
CHTER3 = 'RESO' CCC3 FFF3 ;
* ----------------
* = Résolution H =
* ----------------
* On en déduit la charge moyenne par élément, sous la forme d'un
* champ-point de composante 'H' appuyé aux centres des éléments.
PCEN1 = 'HYBP' MODHYB CND1A CHTER1 ;
PCEN2 = 'HYBP' MODHYB CND1A CHTER2 ;
PCEN3 = 'HYBP' MODHYB CND1A CHTER3 ;
* ----------------
* = Résolution V =
* ----------------
* On en déduit aussi pour chaque jeu de conditions aux limites,
* les débits intégrés sur chaque arête, sous la forme d'un champ-point
* de composante 'FLUX' appuyé aux centre des faces :
QFACE1 = 'HDEB' MODHYB CND1A PCEN1 CHTER1 ;
* La vitesse moyenne dans chaque élément, sous la forme d'un champ-point
* de composantes 'VX', 'VY' (,'VZ') appuyé aux centres des éléments :
VCENT1 = 'HVIT' MODHYB QFACE1 ;
* La vitesse moyenne à travers chaque face, sous la forme d'un champ-point
* de composante 'SCAL' appuyé aux centres des faces :
VFACE1 = ('EXCO' QFACE1 'FLUX' 'SCAL') * HYNOR / HYSUR;
* de même pour le 2e jeux de CLs
QFACE2 = 'HDEB' MODHYB CND1A PCEN2 CHTER2 ;
VCENT2 = 'HVIT' MODHYB QFACE2 ;
VFACE2 = ('EXCO' QFACE2 'FLUX' 'SCAL') * HYNOR / HYSUR;
* de même pour le 3e jeux de CLs
QFACE3 = 'HDEB' MODHYB CND1A PCEN3 CHTER3 ;
VCENT3 = 'HVIT' MODHYB QFACE3 ;
VFACE3 = ('EXCO' QFACE3 'FLUX' 'SCAL') * HYNOR / HYSUR;
* -----------------
* = Calcul ERREUR =
* -----------------
* Comparaison calcul analytique, calcul numérique :
* Erreur relative en trace de charge (moyenne aux faces)
* Erreur relative en charge (moyenne aux centres des éléments)
*
* Jeu de CL n°1
ERRTP1 = 'EXCO' CHTER1 'TH' 'SCAL' ;
ERRTP1 = ERRTP1 - PANAF / PANAF ;
ERRTP1 = 'ABS' ERRTP1 ;
ERRP1 = 'EXCO' PCEN1 'H' 'SCAL' ;
ERRP1 = ERRP1 - PANAC / PANAC ;
ERRP1 = 'ABS' ERRP1 ;
* Jeu de CL n°2
ERRTP2 = 'EXCO' CHTER2 'TH' 'SCAL' ;
ERRTP2 = ERRTP2 - PANAF / PANAF ;
ERRTP2 = 'ABS' ERRTP2 ;
ERRP2 = 'EXCO' PCEN2 'H' 'SCAL' ;
ERRP2 = ERRP2 - PANAC / PANAC ;
ERRP2 = 'ABS' ERRP2 ;
* Jeu de CL n°3
ERRTP3 = 'EXCO' CHTER3 'TH' 'SCAL' ;
ERRTP3 = ERRTP3 - PANAF / PANAF ;
ERRTP3 = 'ABS' ERRTP3 ;
ERRP3 = 'EXCO' PCEN3 'H' 'SCAL' ;
ERRP3 = ERRP3 - PANAC / PANAC ;
ERRP3 = 'ABS' ERRP3 ;
*
* Erreur relative en vitesse aux centres des éléments
MOT1 = 'MOTS' 'VX' 'VY' ;
VDVD = 'PSCA' VANAC VANAC MOT1 MOT1 ;
* Jeu de CL n°1
VD1 = VANAC - VCENT1 ;
VC1 = 'PSCA' VD1 VD1 MOT1 MOT1 ;
SDC1 = 'ABS' ( VC1 / VDVD ) ;
SDC1 = SDC1 '**' 0.5D0 ;
* Jeu de CL n°2
VD2 = VANAC - VCENT2 ;
VC2 = 'PSCA' VD2 VD2 MOT1 MOT1 ;
SDC2 = 'ABS' ( VC2 / VDVD ) ;
SDC2 = SDC2 '**' 0.5D0 ;
* Jeu de CL n°3
VD3 = VANAC - VCENT3 ;
VC3 = 'PSCA' VD3 VD3 MOT1 MOT1 ;
SDC3 = 'ABS' ( VC3 / VDVD ) ;
SDC3 = SDC3 '**' 0.5D0 ;
*
* -------------------
* = Tracé resultats =
* -------------------
'SI' ('NEG' GRAPH 'N') ;
*
*- Transformation des quantités aux centres en MCHAML
*
ERRP1 = 'KCHA' MODHYB 'CHAM' ERRP1 ;
ERRP2 = 'KCHA' MODHYB 'CHAM' ERRP2 ;
ERRP3 = 'KCHA' MODHYB 'CHAM' ERRP3 ;
SDS1 = 'KCHA' MODHYB 'CHAM' SDC1 ;
SDS2 = 'KCHA' MODHYB 'CHAM' SDC2 ;
SDS3 = 'KCHA' MODHYB 'CHAM' SDC3 ;
*
* Dans chaque cas on trace
* L'erreur relative sur la charge moyenne aux centres des éléments
* L'erreur relative sur la Vitesse aux centres des éléments
*
'TITR' 'darcy1/1 : Erreur relative sur la charge' ;
'TRAC' MODHYB ERRP1 ;
'TITR' 'darcy1/1 : Erreur relative sur la vitesse' ;
'TRAC' MODHYB SDS1 ;
*
'TITR' 'darcy1/2 : Erreur relative sur la charge' ;
'TRAC' MODHYB ERRP2 ;
'TITR' 'darcy1/2 : Erreur relative sur la vitesse' ;
'TRAC' MODHYB SDS2 ;
*
'TITR' 'darcy1/3 : Erreur relative sur la charge' ;
'TRAC' MODHYB ERRP3 ;
'TITR' 'darcy1/3 : Erreur relative sur la vitesse' ;
'TRAC' MODHYB SDS3 ;
*
'FINSI' ;
*
* -------------------
* = Gestion ERREURS =
* -------------------
MAXTP1 = 'MAXI' ERRTP1 ;
MAXTP2 = 'MAXI' ERRTP2 ;
MAXTP3 = 'MAXI' ERRTP3 ;
MAXP1 = 'MAXI' ERRP1 ;
MAXP2 = 'MAXI' ERRP2 ;
MAXP3 = 'MAXI' ERRP3 ;
MAXV1 = 'MAXI' SDC1 ;
MAXV2 = 'MAXI' SDC2 ;
MAXV3 = 'MAXI' SDC3 ;
*
'SAUT' 'PAGE' ;
'SAUT' 2 'LIGNE' ;
'MESS' ' ERREURS RELATIVES ' ;
'SAUT' 1 'LIGNE' ;
'MESS' ' cas test TH H V' ;
'SAUT' 1 'LIGNE' ;
'MESS' ' numero1 ' maxtp1 ' ' maxp1 ' ' maxv1 ;
'SAUT' 1 'LIGNE' ;
'MESS' ' numero2 ' maxtp2 ' ' maxp2 ' ' maxv2 ;
'SAUT' 1 'LIGNE' ;
'MESS' ' numero3 ' maxtp3 ' ' maxp3 ' ' maxv3 ;
'SAUT' 2 'LIGNE' ;
*
EPS0 = 1.E-13 ;
LOG1 = MAXTP1 > EPS0 ; LOG2 = MAXTP2 > EPS0 ; LOG3 = MAXTP3 > EPS0 ;
LOG4 = MAXP1 > EPS0 ; LOG5 = MAXP2 > EPS0 ; LOG6 = MAXP3 > EPS0 ;
LOG7 = MAXV1 > EPS0 ; LOG8 = MAXV2 > EPS0 ; LOG9 = MAXV3 > EPS0 ;
LTP0 = LOG1 'OU' LOG2 'OU' LOG3 ;
LP0 = LOG4 'OU' LOG5 'OU' LOG6 ;
LV0 = LOG7 'OU' LOG8 'OU' LOG9 ;
L0 = LTP0 'OU' LP0 'OU' LV0 ;
* Compte-rendu de fin et sortie
'SI' L0 ;
'ERREUR' 5 ;
'SINON' ;
'ERREUR' 0 ;
'FINSI' ;
*
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer