* fichier : darcy2.dgibi
*
******************** CAS TEST : darcy2.dgibi **********************
*
GRAPH = 'N' ;
'SAUT' 'PAGE' ;
*
*-------------------------------------------------------------------
* TEST DARCY2
* CALCUL DARCY ORTHOTROPE
* 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 par une méthode
* d'éléments finis mixtes hybrides dans CASTEM2000.
*
* On effectue un calcul sur un domaine carré, maillé 'aléatoirement'
* par des quadrangles et des triangles. Les conditions aux limites
* sont du type charge imposée.
*
* La solution analytique en charge étant un polynome de degré un,
* la vitesse est constante.
* H(x,y) = -45x -80y + 22.5
* V(x,y) = ( -45*KX ; -80*KY)
*
* On s'attend à peu d'erreur bien que le maillage soit 'aléatoire'.
*
*-------------------------------------------------------------------
*
'SAUT' 'PAGE' ;
*
*- Options générales de calcul.
*
'TITR' 'EFMY DARCY ORTHOTROPE : darcy2.dgibi' ;
'OPTI' 'DIME' 2 'ELEM' 'QUA4' ;
'OPTI' 'ECHO' 0 ;
*
* ------------
* = MAILLAGE =
* ------------
*- Création des points supports des DROITES
*
L = 0.1D0 ;
ML = -0.1D0 ;
A0 = ML ML ;
B0 = L ML ;
A1 = 0.D0 ML ;
B1 = L 0.D0 ;
C0 = L L ;
D0 = ML L ;
C1 = 0.05D0 L ;
D1 = ML -0.05D0 ;
*
*- Création et découpage des DROITES frontieres
*
INUMX = 8 ;
INUMY = 8 ;
INUMP = 2 ;
DRBA1 = A0 'DROI' INUMX A1 ;
DRBA2 = A1 'DROI' INUMP B0 ;
DRBAS = DRBA1 'ET' DRBA2 ;
DRDR1 = B0 'DROI' INUMP B1 ;
DRDR2 = B1 'DROI' INUMY C0 ;
DRDRO = DRDR1 'ET' DRDR2 ;
DRHA1 = C0 'DROI' INUMP C1 ;
DRHA2 = C1 'DROI' INUMX D0 ;
DRHAU = DRHA1 'ET' DRHA2 ;
DRGA1 = D0 'DROI' INUMP D1 ;
DRGA2 = D1 'DROI' INUMY A0 ;
DRGAU = DRGA1 'ET' DRGA2 ;
*
*- Creation maillage GEOMETRIQUE
*
BDHG = DRBAS 'ET' DRDRO 'ET' DRHAU 'ET' DRGAU ;
PTOT1 = 'SURF' BDHG ;
PTOT2 = 'ORIE' PTOT1 ;
QFTOT = CHANGE PTOT2 QUAF ;
*
*- Creation maillage HYBRIDE y compris sous-objets (cond. limites)
*
ELI0 = L / INUMX / 100.D0 ;
*HYTOT = 'DOMA' PTOT2 ELI0 ;
*CHYB1 = 'DOMA' HYTOT 'SURFACE' ;
*CHYB2 = 'DOMA' HYTOT 'NORMALE' ;
*MCHYB = 'DOMA' HYTOT 'ORIENTAT' ;
*HYGAU = 'DOMA' DRGAU 'INCL' HYTOT ELI0 ;
*HYDRO = 'DOMA' DRDRO 'INCL' HYTOT ELI0 ;
*HYHAU = 'DOMA' DRHAU 'INCL' HYTOT ELI0 ;
*HYBAS = 'DOMA' DRBAS 'INCL' HYTOT ELI0 ;
QFGAU = CHANGE DRGAU QUAF ;
QFDRO = CHANGE DRDRO QUAF ;
QFHAU = CHANGE DRHAU QUAF ;
QFBAS = CHANGE DRBAS QUAF ;
ELIM ELI0 ( QFTOT ET QFGAU ET QFDRO ET QFHAU ET QFBAS );
*
* ----------------
* = MODELISATION =
* ----------------
*
*MODHYB = 'MODL' HYTOT 'DARCY' 'ORTHOTROPE' 'HYQ4' 'HYT3' ;
MODHYB = 'MODL' QFTOT 'DARCY' 'ORTHOTROPE' ;
*option donn 5 ;
MODGAU = 'MODL' QFGAU 'DARCY' 'ORTHOTROPE' ;
MODDRO = 'MODL' QFDRO 'DARCY' 'ORTHOTROPE' ;
MODHAU = 'MODL' QFHAU 'DARCY' 'ORTHOTROPE' ;
MODBAS = 'MODL' QFBAS 'DARCY' 'ORTHOTROPE' ;
CEGAU= DOMA MODGAU 'CENTRE' ;
CEDRO= DOMA MODDRO 'CENTRE' ;
CEHAU= DOMA MODHAU 'CENTRE' ;
CEBAS= DOMA MODBAS 'CENTRE' ;
CETOT= DOMA MODHYB 'CENTRE' ;
CHYB1 = 'DOMA' MODHYB 'SURFACE' ;
CHYB2 = 'DOMA' MODHYB 'NORMALE' ;
*
* --------------------------------------
* = Tenseur de perméabilité orthotrope =
* --------------------------------------
*
VKX = 1.D0 ;
VKY = 0.75D0 ;
BZ = 1.D0 0.D0 ;
*
MATI3 = 'MATE' MODHYB 'DIRECTION' BZ 'PARALLELE' 'K1' VKX 'K2' VKY ;
*option donn 5 ;
*
* -----------------------
* = Solution analytique =
* -----------------------
*
* Charge aux faces
AA =-45.D0 ;
BB =-80.D0 ;
CC = 22.5D0 ;
AAA = ( -1.D0 ) * VKX * AA ;
BBB = ( -1.D0 ) * VKY * BB ;
*XX YY = 'COOR' HYTOT.'FACE' ;
XX YY = 'COOR' ( DOMA MODHYB 'FACE') ;
TANF1 = AA * XX ;
TANF2 = BB * YY ;
PANAF = TANF1 + TANF2 + CC ;
*
* 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') ;
TANC1 = AA * XXC ;
TANC2 = BB * YYC ;
PANAC = TANC1 + TANC2 + CC ;
*
* Vitesse aux faces et aux centres
*
VANAF = 'MANU' 'CHPO' ( DOMA MODHYB 'FACE') 2 'VX' AAA 'VY' BBB ;
VANAC = 'MANU' 'CHPO' (DOMA MODHYB 'CENTRE') 2 'VX' AAA 'VY' BBB ;
* -------------------------
* = matrice Masse HYBride =
* -------------------------
CND1A = 'MHYB' MODHYB MATI3 ;
*option donn 5 ;
* -------------------------
* = MAtrice globale en TH =
* -------------------------
HND1A = 'MATP' MODHYB CND1A ;
* ----------------
* = TH imposée =
* ----------------
*BBGAU = 'BLOQ' HYGAU.'CENTRE' 'TH' ;
BBGAU = 'BLOQ' CEGAU 'TH' ;
BBDRO = 'BLOQ' CEDRO 'TH' ;
BBHAU = 'BLOQ' CEHAU 'TH' ;
BBBAS = 'BLOQ' CEBAS 'TH' ;
*
*TTIMP = 'REDU' PANAF HYGAU.'CENTRE' ;
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 ;
* ----------------
* = Assemblage =
* ----------------
CCC1 = HND1A 'ET' BBHAU 'ET' BBBAS 'ET' BBDRO 'ET' BBGAU ;
FFF1 = EEHAU 'ET' EEBAS 'ET' EEDRO 'ET' EEGAU ;
* -----------------
* = Résolution Th =
* -----------------
CHTER1 = 'RESO' CCC1 FFF1 ;
*option donn 5 ;
* ----------------
* = Résolution H =
* ----------------
PCEN1 = 'HYBP' MODHYB CND1A CHTER1 ;
* ----------------
* = Résolution V =
* ----------------
QFACE1 = 'HDEB' MODHYB CND1A PCEN1 CHTER1 ;
VCENT1 = 'HVIT' MODHYB QFACE1 ;
QFACE1 = 'EXCO' QFACE1 'FLUX' 'SCAL' ;
VFACE1 = QFACE1 * CHYB2 / CHYB1 ;
* -----------------
* = Calcul ERREUR =
* -----------------
ERRTP1 = 'EXCO' CHTER1 'TH' 'SCAL' ;
ERRTP1 = ERRTP1 - PANAF / PANAF ;
ERRTP1 = 'ABS' ERRTP1 ;
*
ERRP1 = 'EXCO' PCEN1 'H' 'SCAL' ;
*option donn 5 ;
ERRP1 = ERRP1 - PANAC / PANAC ;
ERRP1 = 'ABS' ERRP1 ;
*
MOT1 = 'MOTS' 'VX' 'VY' ;
MOT2 = 'MOTS' 'UX' 'UY' ;
VAVN = 'PSCA' VANAF CHYB2 MOT1 MOT2 ;
VNN = VAVN * CHYB2 ;
DVDV = 'PSCA' VNN VNN MOT2 MOT2 ;
DVN1 = VNN - VFACE1 ;
DVN1N = 'PSCA' DVN1 DVN1 MOT2 MOT2 ;
SDF1 = 'ABS' ( DVN1N / DVDV ) ;
SDF1 = SDF1 '**' 0.5 ;
*
VDVD = 'PSCA' VANAC VANAC MOT1 MOT1 ;
VD1 = VANAC - VCENT1 ;
VC1 = 'PSCA' VD1 VD1 MOT1 MOT1 ;
SDC1 = 'ABS' ( VC1 / VDVD ) ;
SDC1 = SDC1 '**' 0.5 ;
* -------------------
* = TRACE RESULTATS =
* -------------------
* Dans chaque cas on trace :
* L'erreur relative sur la charge H
* L'erreur relative sur la vitesse V
*
'SI' ('NEG' GRAPH 'N') ;
* ERRP1 = 'KCHA' HYTOT 'CHAM' ERRP1 ;
ERRP1 = 'KCHA' MODHYB 'CHAM' ERRP1 ;
* SDS1 = 'KCHA' HYTOT 'CHAM' SDC1 ;
SDS1 = 'KCHA' MODHYB 'CHAM' SDC1 ;
'TITR' 'darcy2 : erreur relative charge' ;
'TRAC' MODHYB ERRP1 ;
'TITR' 'darcy2 : erreur relative vitesse' ;
'TRAC' MODHYB SDS1 ;
'FINSI' ;
* -------------------
* = Gestion ERREURS =
* -------------------
MAXTP1 = 'MAXI' ERRTP1 ;
MAXP1 = 'MAXI' ERRP1 ;
MAXVF1 = 'MAXI' SDF1 ;
MAXVC1 = 'MAXI' SDC1 ;
*
'SAUT' 'PAGE' ;
'SAUT' 2 'LIGNE' ;
'MESS' ' ERREURS RELATIVES ' ;
'SAUT' 1 'LIGNE' ;
'MESS' ' TH H Vf Vc' ;
'SAUT' 1 'LIGNE' ;
'MESS' ' ' maxtp1 ' ' maxp1 ' ' maxvf1 ' ' maxvc1 ;
'SAUT' 2 'LIGNE' ;
*
EPS0 = 1.E-12 ; EPS1 = 1.E-10 ;
LOG1 = MAXTP1 > EPS0 ;
LOG2 = MAXP1 > EPS0 ;
LOG3 = MAXVF1 > EPS1 ;
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