* fichier : darcy3_prisme_EFMH.dgibi ************************************************************************ * Section : Fluides Darcy ************************************************************************ * ************************** 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' 'PRI6' ; '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 = 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' ( 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 / 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-12 ; 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' ;