Accueil | Cast3M

* fichier :  darcy3EFMH.dgibi
************************************************************************
************************************************************************
*
************************** CAS TEST : darcy3.dgibi ******************
*
 GRAPH = 'N'   ;
'SAUT' 'PAGE' ;
NITER = 40 ;
*
*-------------------------------------------------------------------
* TEST DARCY3
*                  CALCUL DARCY ORTHOTROPE 3D
*
* TEST TRANGEOL en EFMH
* staionnaire comme limite d'un transitoire
*
* 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 = 6          ;
ENY = 6          ;
ENZ = 6          ;
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 = MODE QFTOT 'DARCY' 'ANISOTROPE'     ;
MODGAU = MODE QFGAU 'DARCY' 'ANISOTROPE'     ;
MODDRO = MODE QFDRO 'DARCY' 'ANISOTROPE'     ;
MODHAU = MODE QFHAU 'DARCY' 'ANISOTROPE'     ;
MODBAS = MODE QFBAS 'DARCY' 'ANISOTROPE'     ;
MODDEV = MODE QFDEV 'DARCY' 'ANISOTROPE'     ;
MODDER = MODE 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' (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'            = 10.D0 / NITER         ;
GEOL1 . 'DIFFUSIVITE'       = mati3                 ;
GEOL1 . 'SOLVEUR'           = 3                     ;
GEOL1 . 'PRECONDITIONNEUR'  = 3                     ;
GEOL1 . 'POROSITE'          = 'MANU' 'CHPO' 
                (doma MODHYB CENTRE) 'CK' 1.D0 ;
GEOL1 . 'CLIMITES'          = CHCLIM                ;
GEOL1 . 'RECALCUL'          = VRAI                  ;
*
GEOL1 GEOL2 = TRANGEOL Modhyb GEOL1;
 
 
 
* boucle transitoire
 
 
 
REPE blocc NITER   ;
   mess 'boucle' &blocc;
   GEOL1 . 'RECALCUL'          = FAUX                  ;
   GEOL1 GEOL2 = TRANGEOL Modhyb GEOL1 GEOL2;
FIN blocc;
 
 
 
 
CHTER3 = 'NOMC' 'TH' GEOL2 . 'TRACE_CONC';
PCEN3 = 'NOMC' 'H' GEOL1 . 'CONCENTRATION' ;
QFACE3 = 'NOMC' 'FLUX' GEOL1 . 'FLUXDIFF' ;
 
*
*- Calcul de V
*
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
*
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 ;
*
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.
*
ERRP3 = 'KCHA' MODHYB 'CHAM' ERRP3 ;
SDS3  = 'KCHA' MODHYB 'CHAM' SDC3  ;
*
*   L'erreur relative sur la charge au centre
*   L'erreur relative sur la Vitesse au centre
*
'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 =
*                                                 -------------------
MAXTP3 = 'MAXI' ERRTP3 ;
MAXP3  = 'MAXI' ERRP3  ;
MAXV3  = 'MAXI' SDC3   ;
*
'SAUT' 'PAGE'                                                       ;
'SAUT' 2 'LIGNE'                                                    ;
'MESS' '                             ERREURS RELATIVES            ' ;
'SAUT' 1 'LIGNE'                                                    ;
'MESS' ' cas test          TH                  H                 V' ;
'SAUT' 1 'LIGNE'                                                    ;
'MESS' ' numero3     ' maxtp3  '   ' maxp3 '   ' maxv3              ;
'SAUT' 2 'LIGNE'                                                    ;
*
EPS0 = 1.D-7 ;
LOG3 = MAXTP3 > EPS0 ;
LOG6 = MAXP3  > EPS0 ;
LOG9 = MAXV3  > EPS0 ;
LTP0 = LOG3 ;
LP0  = LOG6 ;
LV0  = LOG9 ;
L0   = LTP0 'OU' LP0  'OU' LV0  ;
'SI' ( L0 ) ;
   'ERRE' 5 ;
'SINO'      ;
   'ERRE' 0 ;
'FINSI'     ;
*
'FIN' ;