* 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 = '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' (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' ;