Download darcy3_tetraedre_EFMH.dgibi
* fichier : darcy3_tetraedre_EFMH.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' 'TET4' ;
'OPTI' 'ECHO' 1 ;
*
* ------------
* = 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 = 10 ;
ENY = 10 ;
ENZ = 10 ;
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 ;
'SI' ('EGA' ('VALEUR' 'ELEM') 'CUB8') ;
CUBE1 = 'PAVE' SDER SDEV SBAS SHAU SGAU SDRO ;
'SINON' ;
CUBE1 = 'VOLU' (SDER 'ET' SDEV 'ET' SBAS 'ET' SHAU
'ET' SGAU 'ET' SDRO) ;
'FINSI' ;
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)
*
MODHYB = 'MODL' QFTOT 'DARCY' 'ANISOTROPE' ;
MODGAU = 'MODL' QFGAU 'DARCY' 'ANISOTROPE' ;
MODDRO = 'MODL' QFDRO 'DARCY' 'ANISOTROPE' ;
MODHAU = 'MODL' QFHAU 'DARCY' 'ANISOTROPE' ;
MODBAS = 'MODL' QFBAS 'DARCY' 'ANISOTROPE' ;
MODDEV = 'MODL' QFDEV 'DARCY' 'ANISOTROPE' ;
MODDER = 'MODL' QFDER 'DARCY' 'ANISOTROPE' ;
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 ;
VVKX = MANU 'CHPO' (DOMA MODHYB 'CENTRE') 'K11' 1.D0 ;
VVKX = CHANGER ATTRIBUT VVKX 'NATU' DISCRET ;
VVKY = MANU 'CHPO' (DOMA MODHYB 'CENTRE') 'K22' 0.75D0 ;
VVKY = CHANGER ATTRIBUT VVKY 'NATU' DISCRET ;
VVKZ = MANU 'CHPO' (DOMA MODHYB 'CENTRE') 'K33' 0.5D0 ;
VVKZ = CHANGER ATTRIBUT VVKZ 'NATU' DISCRET ;
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 =
* --------------
*
MATI3 = (NOMC 'K11' VVKX) et
(NOMC 'K22' VVKY) et
(NOMC 'K33' VVKZ) et
(NOMC 'K21' (0.D0 * VVKX)) et
(NOMC 'K31' (0.D0 * VVKX)) et
(NOMC 'K32' (0.D0 * VVKX)) ;
* ;
*- 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�
*
FLDRO = -60.D0 * (nomc 'FLUX' (doma moddro VOLUME));
FLGAU = 60.D0 * (nomc 'FLUX' (doma modgau VOLUME));
FLHAU = 30.D0 * (nomc 'FLUX' (doma modhau VOLUME));
FLBAS = -30.D0 * (nomc 'FLUX' (doma modbas VOLUME));
FLDEV = 45.D0 * (nomc 'FLUX' (doma moddev VOLUME));
FLDER = -45.D0 * (nomc 'FLUX' (doma modder VOLUME));
*
*
*- Assemblage et r�solution en TH
*
h_lim = 'RESOUD' ( BBHAU 'ET' BBBAS 'ET' BBDRO 'ET'
BBGAU 'ET' BBDEV 'ET' BBDER)
(EEHAU 'ET' EEBAS 'ET' EEDRO 'ET'
EEGAU 'ET' EEDEV 'ET' EEDER) ;
h_lim = 'EXCO' 'TH' h_lim;
CHCLIM = TABLE;
CHCLIM . 'DIRICHLET' = 'NOMC' 'I35' h_lim;
GEOL1 = TABLE;
GEOL1 . 'CONCENTRATION' = 'NOMC' 'I35' (0.D0 * PANAC) ;
GEOL1 . 'LUMP' = FAUX ;
GEOL1 . 'TYPDISCRETISATION' = 'EFMH' ;
GEOL1 . 'THETA_DIFFUSION' = 1.0D0 ;
GEOL1 . 'THETA_CONVECTION' = 1.0D0 ;
GEOL1 . 'DECENTREMENT' = FAUX ;
GEOL1 . 'DELTAT' = 1.D15 ;
GEOL1 . 'DIFFUSIVITE' = mati3 ;
GEOL1 . 'SOLVEUR' = 2 ;
GEOL1 . 'PRECONDITIONNEUR' = 3 ;
GEOL1 . 'POROSITE' = 0.D0 * PANAC ;
GEOL1 . 'CLIMITES' = CHCLIM ;
GEOL1 . 'RECALCUL' = VRAI ;
*
GEOL1 GEOL2 = TRANGEOL Modhyb GEOL1;
CHTER1 = 'NOMC' 'TH' GEOL2 . 'TRACE_CONC';
PCEN1 = 'NOMC' 'H' GEOL1 . 'CONCENTRATION' ;
QFACE1 = 'NOMC' 'FLUX' GEOL1 . 'FLUXDIFF' ;
* deuxieme probleme
h_lim = 'RESOUD' (BBHAU 'ET' BBBAS 'ET' BBGAU 'ET' BBDRO)
(EEHAU 'ET' EEBAS 'ET' EEGAU 'ET' EEDRO) ;
h_lim = 'EXCO' 'TH' h_lim;
CHCLIM = TABLE;
CHCLIM . 'NEUMANN' = ('NOMC' 'I35' (FLDEV 'ET' FLDER));
CHCLIM . 'DIRICHLET' = 'NOMC' 'I35' h_lim;
GEOL1 = TABLE;
GEOL1 . 'CONCENTRATION' = 'NOMC' 'I35' (0.D0 * PANAC) ;
GEOL1 . 'LUMP' = FAUX ;
GEOL1 . 'TYPDISCRETISATION' = 'EFMH' ;
GEOL1 . 'THETA_DIFFUSION' = 1.0D0 ;
GEOL1 . 'THETA_CONVECTION' = 1.0D0 ;
GEOL1 . 'DECENTREMENT' = FAUX ;
GEOL1 . 'DELTAT' = 1.D15 ;
GEOL1 . 'DIFFUSIVITE' = mati3 ;
GEOL1 . 'SOLVEUR' = 2 ;
GEOL1 . 'PRECONDITIONNEUR' = 3 ;
GEOL1 . 'POROSITE' = 0.D0 * PANAC ;
GEOL1 . 'CLIMITES' = CHCLIM ;
GEOL1 . 'RECALCUL' = VRAI ;
*
GEOL1 GEOL2 = TRANGEOL Modhyb GEOL1;
CHTER2 = 'NOMC' 'TH' GEOL2 . 'TRACE_CONC';
PCEN2 = 'NOMC' 'H' GEOL1 . 'CONCENTRATION' ;
QFACE2 = 'NOMC' 'FLUX' GEOL1 . 'FLUXDIFF' ;
* troisieme probleme
h_lim = 'RESOUD' (BBGAU 'ET' BBDEV)
(EEGAU 'ET' EEDEV) ;
h_lim = 'EXCO' 'TH' h_lim;
CHCLIM = TABLE;
CHCLIM . 'NEUMANN' = ('NOMC' 'I35' ( FLHAU 'ET' FLBAS
'ET' FLDRO 'ET' FLDER));
CHCLIM . 'DIRICHLET' = 'NOMC' 'I35' h_lim;
GEOL1 = TABLE;
GEOL1 . 'CONCENTRATION' = 'NOMC' 'I35' (0.D0 * PANAC) ;
GEOL1 . 'LUMP' = FAUX ;
GEOL1 . 'TYPDISCRETISATION' = 'EFMH' ;
GEOL1 . 'THETA_DIFFUSION' = 1.0D0 ;
GEOL1 . 'THETA_CONVECTION' = 1.0D0 ;
GEOL1 . 'DECENTREMENT' = FAUX ;
GEOL1 . 'DELTAT' = 1.D15 ;
GEOL1 . 'DIFFUSIVITE' = mati3 ;
GEOL1 . 'SOLVEUR' = 2 ;
GEOL1 . 'PRECONDITIONNEUR' = 3 ;
GEOL1 . 'POROSITE' = 0.D0 * PANAC ;
GEOL1 . 'CLIMITES' = CHCLIM ;
GEOL1 . 'RECALCUL' = VRAI ;
*
GEOL1 GEOL2 = TRANGEOL Modhyb GEOL1;
GEOL1 GEOL2 = TRANGEOL Modhyb GEOL1 GEOL2;
CHTER3 = 'NOMC' 'TH' GEOL2 . 'TRACE_CONC';
PCEN3 = 'NOMC' 'H' GEOL1 . 'CONCENTRATION' ;
QFACE3 = 'NOMC' 'FLUX' GEOL1 . 'FLUXDIFF' ;
*
*- Calcul de V
*
VCENT1 = 'HVIT' MODHYB QFACE1 ;
QFACE1 = 'EXCO' QFACE1 'FLUX' 'SCAL' ;
VFACE1 = QFACE1 * CHYB2 / CHYB1 ;
*
VCENT2 = 'HVIT' MODHYB QFACE2 ;
QFACE2 = 'EXCO' QFACE2 'FLUX' 'SCAL' ;
VFACE2 = QFACE2 * CHYB2 / CHYB1 ;
*
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) ** 2 ;
ERRTP1 = 'RESULT' ERRTP1 ;
ERRTP1 = ERRTP1 '/' ('RESULT' (PANAF ** 2));
ERRTP1 = ('MAXIMUM' errtp1)**0.5D0 ;
ERRP1 = 'EXCO' PCEN1 'H' 'SCAL' ;
ERRP1 = (ERRP1 - PANAC) ** 2 ;
ERRP1 = 'RESULT' ERRP1 ;
ERRP1 = ERRP1 '/' ('RESULT' (PANAC ** 2));
ERRP1 = ('MAXIMUM' errp1)**0.5D0 ;
*
ERRTP2 = 'EXCO' CHTER2 'TH' 'SCAL' ;
ERRTP2 = (ERRTP2 - PANAF) ** 2 ;
ERRTP2 = 'RESULT' ERRTP2 ;
ERRTP2 = ERRTP2 '/' ('RESULT' (PANAF ** 2));
ERRTP2 = ('MAXIMUM' errtp2)**0.5D0 ;
ERRP2 = 'EXCO' PCEN2 'H' 'SCAL' ;
ERRP2 = (ERRP2 - PANAC) ** 2 ;
ERRP2 = 'RESULT' ERRP2 ;
ERRP2 = ERRP2 '/' ('RESULT' (PANAC ** 2));
ERRP2 = ('MAXIMUM' errp2)**0.5D0 ;
*
*
ERRTP3 = 'EXCO' CHTER3 'TH' 'SCAL' ;
ERRTP3 = (ERRTP3 - PANAF) ** 2 ;
ERRTP3 = 'RESULT' ERRTP3 ;
ERRTP3 = ERRTP3 '/' ('RESULT' (PANAF ** 2));
ERRTP3 = ('MAXIMUM' errtp3)**0.5D0 ;
ERRP3 = 'EXCO' PCEN3 'H' 'SCAL' ;
ERRP3 = (ERRP3 - PANAC) ** 2 ;
ERRP3 = 'RESULT' ERRP3 ;
ERRP3 = ERRP3 '/' ('RESULT' (PANAC ** 2));
ERRP3 = ('MAXIMUM' errp3)**0.5D0 ;
*
*
MOT1 = 'MOTS' 'VX' 'VY' 'VZ' ;
VDVD = 'PSCA' VANAC VANAC MOT1 MOT1 ;
*
VD1 = VANAC - VCENT1 ;
VC1 = 'PSCA' VD1 VD1 MOT1 MOT1 ;
SDC1 = 'MAXIMUM' ( ('RESULT' VC1) / ('RESULT' VDVD )) ;
SDC1 = SDC1 '**' 0.5 ;
VD2 = VANAC - VCENT2 ;
VC2 = 'PSCA' VD2 VD2 MOT1 MOT1 ;
SDC2 = 'MAXIMUM' ( ('RESULT' VC2) / ('RESULT' VDVD) ) ;
SDC2 = SDC2 '**' 0.5 ;
VD3 = VANAC - VCENT3 ;
VC3 = 'PSCA' VD3 VD3 MOT1 MOT1 ;
SDC3 = 'MAXIMUM' ( ('RESULT' VC3) / ('RESULT' 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 = ERRTP1 ;
MAXTP2 = ERRTP2 ;
MAXTP3 = ERRTP3 ;
MAXP1 = ERRP1 ;
MAXP2 = ERRP2 ;
MAXP3 = ERRP3 ;
MAXV1 = SDC1 ;
MAXV2 = SDC2 ;
MAXV3 = 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-8 ;
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