* fichier : darcy3.dgibi
*
************************** CAS TEST : darcy3.dgibi ******************
*
GRAPH = 'N' ;
'SAUT' 'PAGE' ;
*
*-------------------------------------------------------------------
* TEST DARCY3
* CALCUL DARCY ORTHOTROPE 3D
* 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é afin de résoudre les équations de DARCY par une méthode
* d'éléments finis mixtes hybrides dans CASTEM2000.
*
* On effectue trois calculs sur un cube, maillé par des cubes
* réguliers.
* 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 est un polynome de degré un,
* la vitesse est constante et la conductivité hydraulique orthotrope.
* __ __
* | |
* | 1 0 0 |
* K = | 0 3/4 0 |
* | 0 0 1/2 |
* |__ __|
*
* H(x,y,z) = -45 x -80 y -60z + 200.
* V(x,y,z) = ( 45 ; 60 ; 30 )
*
* 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 ORTHOTROPE 3D Lineaire : darcy3.dgibi' ;
'OPTI' 'DIME' 3 'ELEM' 'CUB8' ;
'OPTI' 'ECHO' 0 ;
*
* ------------
* = MAILLAGE =
* ------------
OEIL = 5.D0 6.D0 7.D0 ;
VECX = 1.D0 0.D0 0.D0 ;
VECY = 0.D0 1.D0 0.D0 ;
VECZ = 0.D0 0.D0 1.D0 ;
*
ENX = 5 ;
ENY = 10 ;
ENZ = 5 ;
DX = 1.D0 / ENX ;
DY = 1.D0 / ENY ;
DZ = 1.D0 / ENZ ;
*
*- Création des points
*
A0 = 0.D0 0.D0 0.D0 ;
B0 = 1.D0 0.D0 0.D0 ;
C0 = 1.D0 1.D0 0.D0 ;
D0 = 0.D0 1.D0 0.D0 ;
E0 = 0.D0 0.D0 1.D0 ;
F0 = 1.D0 0.D0 1.D0 ;
G0 = 1.D0 1.D0 1.D0 ;
H0 = 0.D0 1.D0 1.D0 ;
*
*- Création des droites
*
AB = 'DROI' ENX A0 B0 ;
AD = 'DROI' ENY A0 D0 ;
AE = 'DROI' ENZ A0 E0 ;
BC = 'DROI' ENY B0 C0 ;
BF = 'DROI' ENZ B0 F0 ;
CD = 'DROI' ENX C0 D0 ;
CG = 'DROI' ENZ C0 G0 ;
DH = 'DROI' ENZ D0 H0 ;
EF = 'DROI' ENX E0 F0 ;
EH = 'DROI' ENY E0 H0 ;
FG = 'DROI' ENY F0 G0 ;
GH = 'DROI' ENX G0 H0 ;
*
FE = 'INVE' EF ;
EA = 'INVE' AE ;
DC = 'INVE' CD ;
HD = 'INVE' DH ;
DA = 'INVE' AD ;
HE = 'INVE' EH ;
GF = 'INVE' FG ;
FB = 'INVE' BF ;
*
*- Creation des faces du cube
*
SDRO = 'DALL' AB BF FE EA 'PLAN' ;
SGAU = 'DALL' DC CG GH HD 'PLAN' ;
SBAS = 'DALL' AB BC CD DA 'PLAN' ;
SHAU = 'DALL' EF FG GH HE 'PLAN' ;
SDEV = 'DALL' BC CG GF FB 'PLAN' ;
SDER = 'DALL' AD DH HE EA 'PLAN' ;
*
*- Création maillage géométrique
*
ENXM = ENX + ENY + ENZ ;
ELI0 = 1.D0 / ENXM / 10.D0 ;
CUBE1 = 'PAVE' SDER SDEV SBAS SHAU SGAU SDRO ;
QFTOT = CHANGE CUBE1 QUAF ;
QFGAU = CHANGE SGAU QUAF ;
QFDRO = CHANGE SDRO QUAF ;
QFHAU = CHANGE SHAU QUAF ;
QFBAS = CHANGE SBAS QUAF ;
QFDEV = CHANGE SDEV QUAF ;
QFDER = CHANGE SDER QUAF ;
ELIM ELI0 (QFTOT ET QFGAU ET QFDRO ET QFHAU ET QFBAS ET QFDEV ET
QFDER ) ;
*
*- Création maillage HYBRIDE et sous-objets (conditions aux limites)
*
*HYTOT = 'DOMA' CUBE1 ELI0 ;
*C11 = 'DOMA' HYTOT 'VOLUME' ;
*CHYB1 = 'DOMA' HYTOT 'SURFACE' ;
*CHYB2 = 'DOMA' HYTOT 'NORMALE' ;
*MCHYB = 'DOMA' HYTOT 'ORIENTAT' ;
*
*HYGAU = 'DOMA' SGAU 'INCL' HYTOT ELI0 ;
*HYDRO = 'DOMA' SDRO 'INCL' HYTOT ELI0 ;
*HYHAU = 'DOMA' SHAU 'INCL' HYTOT ELI0 ;
*HYBAS = 'DOMA' SBAS 'INCL' HYTOT ELI0 ;
*HYDEV = 'DOMA' SDEV 'INCL' HYTOT ELI0 ;
*HYDER = 'DOMA' SDER 'INCL' HYTOT ELI0 ;
MODHYB = 'MODL' QFTOT 'DARCY' 'ORTHOTROPE' ;
MODGAU = 'MODL' QFGAU 'DARCY' 'ORTHOTROPE' ;
MODDRO = 'MODL' QFDRO 'DARCY' 'ORTHOTROPE' ;
MODHAU = 'MODL' QFHAU 'DARCY' 'ORTHOTROPE' ;
MODBAS = 'MODL' QFBAS 'DARCY' 'ORTHOTROPE' ;
MODDEV = 'MODL' QFDEV 'DARCY' 'ORTHOTROPE' ;
MODDER = 'MODL' QFDER 'DARCY' 'ORTHOTROPE' ;
C11 = 'DOMA' MODHYB 'VOLUME' ;
CHYB1 = 'DOMA' MODHYB 'SURFACE' ;
CHYB2 = 'DOMA' MODHYB 'NORMALE' ;
CEGAU = 'DOMA' MODGAU 'CENTRE' ;
CEDRO = 'DOMA' MODDRO 'CENTRE' ;
CEHAU = 'DOMA' MODHAU 'CENTRE' ;
CEBAS = 'DOMA' MODBAS 'CENTRE' ;
CEDEV = 'DOMA' MODDEV 'CENTRE' ;
CEDER = 'DOMA' MODDER 'CENTRE' ;
*
*- Solution analytique
*
XX YY ZZ = 'COOR' (DOMA MODHYB 'FACE' ) ;
XXC YYC ZZC = 'COOR' (DOMA MODHYB 'CENTRE') ;
*
VKX = 1.D0 ;
VKY = 0.75D0 ;
VKZ = 0.5D0 ;
AA = -45.D0 ;
BB = -80.D0 ;
CC = -60.D0 ;
DD = 200.D0 ;
AAA = -1.D0 * VKX * AA ;
BBB = -1.D0 * VKY * BB ;
CCC = -1.D0 * VKZ * CC ;
*
PANAF = (AA * XX) + (BB * YY) + (CC * ZZ) + DD ;
PANAC = (AA * XXC) + (BB * YYC) + (CC * ZZC) + DD ;
VANAC = 'MANU' 'CHPO' (DOMA MODHYB 'CENTRE') 3 'VX' AAA
'VY' BBB 'VZ' CCC ;
VANAF = 'MANU' 'CHPO' (DOMA MODHYB 'FACE') 3 'VX' AAA
'VY' BBB 'VZ' CCC ;
*
* --------------
* = RESOLUTION =
* --------------
*
*MODHYB = 'MODL' HYTOT 'DARCY' 'ORTHOTROPE' 'HYC8' ;
MATI3 = 'MATE' MODHYB 'DIRECTION' VECX VECY 'PARALLELE'
'K1' VKX 'K2' VKY 'K3' VKZ ;
CND1A = 'MHYB' MODHYB MATI3 ;
HND1A = 'MATP' MODHYB CND1A ;
* ;
*- Conditions aux limites
*
BBGAU = 'BLOQ' CEGAU 'TH' ;
BBDRO = 'BLOQ' CEDRO 'TH' ;
BBHAU = 'BLOQ' CEHAU 'TH' ;
BBBAS = 'BLOQ' CEBAS 'TH' ;
BBDEV = 'BLOQ' CEDEV 'TH' ;
BBDER = 'BLOQ' CEDER 'TH' ;
*
*- TH imposée
*
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 ;
TTIMP = 'REDU' PANAF CEDEV ;
TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH' ;
EEDEV = 'DEPI' BBDEV TTIM2 ;
TTIMP = 'REDU' PANAF CEDER ;
TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH' ;
EEDER = 'DEPI' BBDER TTIM2 ;
*
*- Flux imposé
*
FLUD = -60.D0 * DX * DZ ;
FLUG = -1.D0 * FLUD ;
FLUH = 30.D0 * DX * DY ;
FLUB = -1.D0 * FLUH ;
FLUV = 45.D0 * DY * DZ ;
FLUR = -1.D0 * FLUV ;
*
FLGAU = 'MANU' 'CHPO' CEGAU 1 'FLUX' FLUG 'NATURE' 'DISCRET' ;
FLDRO = 'MANU' 'CHPO' CEDRO 1 'FLUX' FLUD 'NATURE' 'DISCRET' ;
FLBAS = 'MANU' 'CHPO' CEBAS 1 'FLUX' FLUB 'NATURE' 'DISCRET' ;
FLHAU = 'MANU' 'CHPO' CEHAU 1 'FLUX' FLUH 'NATURE' 'DISCRET' ;
FLDEV = 'MANU' 'CHPO' CEDEV 1 'FLUX' FLUV 'NATURE' 'DISCRET' ;
FLDER = 'MANU' 'CHPO' CEDER 1 'FLUX' FLUR 'NATURE' 'DISCRET' ;
*
*- Assemblage et résolution en TH
*
CCC1 = HND1A 'ET' BBHAU 'ET' BBBAS 'ET' BBDRO 'ET'
BBGAU 'ET' BBDEV 'ET' BBDER ;
FFF1 = EEHAU 'ET' EEBAS 'ET' EEDRO 'ET'
EEGAU 'ET' EEDEV 'ET' EEDER ;
*
CCC2 = HND1A 'ET' BBHAU 'ET' BBBAS 'ET' BBGAU 'ET' BBDRO ;
FFF2 = FLDEV 'ET' FLDER 'ET' EEHAU 'ET' EEBAS 'ET' EEGAU 'ET' EEDRO ;
*
CCC3 = HND1A 'ET' BBGAU 'ET' BBDEV ;
FFF3 = FLHAU 'ET' FLBAS 'ET' FLDRO 'ET' FLDER 'ET' EEGAU 'ET' EEDEV ;
*
CHTER1 = 'RESO' CCC1 FFF1 ;
CHTER2 = 'RESO' CCC2 FFF2 ;
CHTER3 = 'RESO' CCC3 FFF3 ;
*
*- Calcul de H
*
PCEN1 = 'HYBP' MODHYB CND1A CHTER1 ;
PCEN2 = 'HYBP' MODHYB CND1A CHTER2 ;
PCEN3 = 'HYBP' MODHYB CND1A CHTER3 ;
*
*- Calcul de V
*
QFACE1 = 'HDEB' MODHYB CND1A PCEN1 CHTER1 ;
VCENT1 = 'HVIT' MODHYB QFACE1 ;
QFACE1 = 'EXCO' QFACE1 'FLUX' 'SCAL' ;
VFACE1 = QFACE1 * CHYB2 / CHYB1 ;
*
QFACE2 = 'HDEB' MODHYB CND1A PCEN2 CHTER2 ;
VCENT2 = 'HVIT' MODHYB QFACE2 ;
QFACE2 = 'EXCO' QFACE2 'FLUX' 'SCAL' ;
VFACE2 = QFACE2 * CHYB2 / CHYB1 ;
*
QFACE3 = 'HDEB' MODHYB CND1A PCEN3 CHTER3 ;
VCENT3 = 'HVIT' MODHYB QFACE3 ;
QFACE3 = 'EXCO' QFACE3 'FLUX' 'SCAL' ;
VFACE3 = QFACE3 * CHYB2 / CHYB1 ;
*
* -----------------
* = Calcul ERREUR =
* -----------------
* ERReur relative en Trace de charge TH aux faces des éléments
* ERReur relative en charge H au centre des éléments
* Erreur relative sur la vitesse au centre des éléments
*
ERRTP1 = 'EXCO' CHTER1 'TH' 'SCAL' ;
ERRTP1 = ERRTP1 - PANAF / PANAF ;
ERRTP1 = 'ABS' ERRTP1 ;
ERRP1 = 'EXCO' PCEN1 'H' 'SCAL' ;
ERRP1 = ERRP1 - PANAC / PANAC ;
ERRP1 = 'ABS' ERRP1 ;
*
ERRTP2 = 'EXCO' CHTER2 'TH' 'SCAL' ;
ERRTP2 = ERRTP2 - PANAF / PANAF ;
ERRTP2 = 'ABS' ERRTP2 ;
ERRP2 = 'EXCO' PCEN2 'H' 'SCAL' ;
ERRP2 = ERRP2 - PANAC / PANAC ;
ERRP2 = 'ABS' ERRP2 ;
*
ERRTP3 = 'EXCO' CHTER3 'TH' 'SCAL' ;
ERRTP3 = ERRTP3 - PANAF / PANAF ;
ERRTP3 = 'ABS' ERRTP3 ;
ERRP3 = 'EXCO' PCEN3 'H' 'SCAL' ;
ERRP3 = ERRP3 - PANAC / PANAC ;
ERRP3 = 'ABS' ERRP3 ;
*
MOT1 = 'MOTS' 'VX' 'VY' 'VZ' ;
VDVD = 'PSCA' VANAC VANAC MOT1 MOT1 ;
*
VD1 = VANAC - VCENT1 ;
VC1 = 'PSCA' VD1 VD1 MOT1 MOT1 ;
SDC1 = 'ABS' ( VC1 / VDVD ) ;
SDC1 = SDC1 '**' 0.5 ;
VD2 = VANAC - VCENT2 ;
VC2 = 'PSCA' VD2 VD2 MOT1 MOT1 ;
SDC2 = 'ABS' ( VC2 / VDVD ) ;
SDC2 = SDC2 '**' 0.5 ;
VD3 = VANAC - VCENT3 ;
VC3 = 'PSCA' VD3 VD3 MOT1 MOT1 ;
SDC3 = 'ABS' ( VC3 / VDVD ) ;
SDC3 = SDC3 '**' 0.5 ;
*
* -------------------
* = Tracé resultats =
* -------------------
'SI' ('NEG' GRAPH 'N') ;
*
*- Transformation des quantités aux centres en MCHAML constant.
*
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 au centre
* L'erreur relative sur la Vitesse au centre
*
'TITR' 'darcy3/1 : Erreur relative sur la charge' ;
'TRAC' MODHYB ERRP1 ;
'TITR' 'darcy3/1 : Erreur relative sur la vitesse' ;
'TRAC' MODHYB SDS1 ;
*
'TITR' 'darcy3/2 : Erreur relative sur la charge' ;
'TRAC' MODHYB ERRP2 ;
'TITR' 'darcy3/2 : Erreur relative sur la vitesse' ;
'TRAC' MODHYB SDS2 ;
*
'TITR' 'darcy3/3 : Erreur relative sur la charge' ;
'TRAC' MODHYB ERRP3 ;
'TITR' 'darcy3/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.D-7 ;
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 ;
'SI' ( L0 ) ;
'ERRE' 5 ;
'SINO' ;
'ERRE' 0 ;
'FINSI' ;
*
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer