Télécharger darcy3_tetraedre_VF.dgibi
* fichier : darcy3_tetraedre_VF.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. * * * ------------ * = 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 * * * *- Creation des faces du cube * * *- Création maillage géométrique * ENXM = ENX + ENY + ENZ ; ELI0 = 1.D0 / ENXM / 10.D0 ; 'SINON' ; '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 ; QFDER ) ; * *- Création maillage HYBRIDE et sous-objets (conditions aux limites) * * *- Solution analytique * * VKX = 1.D0 ; VKY = 0.75D0 ; VKZ = 0.5D0 ; VVKX = CHANGER ATTRIBUT VVKX 'NATU' DISCRET ; VVKY = CHANGER ATTRIBUT VVKY 'NATU' DISCRET ; 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 ; 'VY' BBB 'VZ' CCC ; 'VY' BBB 'VZ' CCC ; * * -------------- * = RESOLUTION = * -------------- * * ; *- Conditions aux limites * * *- TH imposée * * *- Flux imposé * * * *- 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) ; CHCLIM = TABLE; GEOL1 = TABLE; GEOL1 . 'TYPDISCRETISATION' = 'VF' ; 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 ; * * deuxieme probleme h_lim = 'RESOUD' (BBHAU 'ET' BBBAS 'ET' BBGAU 'ET' BBDRO) (EEHAU 'ET' EEBAS 'ET' EEGAU 'ET' EEDRO) ; CHCLIM = TABLE; GEOL1 = TABLE; 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 ; * * troisieme probleme h_lim = 'RESOUD' (BBGAU 'ET' BBDEV) (EEGAU 'ET' EEDEV) ; CHCLIM = TABLE; 'ET' FLDRO 'ET' FLDER)); GEOL1 = TABLE; 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 ; * * *- Calcul de V * VFACE1 = QFACE1 * CHYB2 / CHYB1 ; * VFACE2 = QFACE2 * CHYB2 / CHYB1 ; * VFACE3 = QFACE3 * CHYB2 / CHYB1 ; * * ----------------- * = Calcul ERREUR = * ----------------- * ERReur relative en charge H au centre des éléments * Erreur relative sur la vitesse au centre des éléments * ERRP1 = (ERRP1 - PANAC) ** 2 ; ERRP1 = 'RESULT' ERRP1 ; ERRP1 = ERRP1 '/' ('RESULT' (PANAC ** 2)); ERRP1 = ('MAXIMUM' errp1)**0.5D0 ; * ERRP2 = (ERRP2 - PANAC) ** 2 ; ERRP2 = 'RESULT' ERRP2 ; ERRP2 = ERRP2 '/' ('RESULT' (PANAC ** 2)); ERRP2 = ('MAXIMUM' errp2)**0.5D0 ; * ERRP3 = (ERRP3 - PANAC) ** 2 ; ERRP3 = 'RESULT' ERRP3 ; ERRP3 = ERRP3 '/' ('RESULT' (PANAC ** 2)); ERRP3 = ('MAXIMUM' errp3)**0.5D0 ; * * * VD1 = VANAC - VCENT1 ; SDC1 = 'MAXIMUM' ( ('RESULT' VC1) / ('RESULT' VDVD )) ; SDC1 = SDC1 '**' 0.5 ; VD2 = VANAC - VCENT2 ; SDC2 = 'MAXIMUM' ( ('RESULT' VC2) / ('RESULT' VDVD) ) ; SDC2 = SDC2 '**' 0.5 ; VD3 = VANAC - VCENT3 ; 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. * * * Dans chaque cas on trace * L'erreur relative sur la charge au centre * L'erreur relative sur la Vitesse au centre * 'TRAC' MODHYB ERRP1 ; 'TRAC' MODHYB SDS1 ; * 'TRAC' MODHYB ERRP2 ; 'TRAC' MODHYB SDS2 ; * 'TRAC' MODHYB ERRP3 ; 'TRAC' MODHYB SDS3 ; * 'FINSI' ; * * ------------------- * = Gestion ERREURS = * ------------------- MAXP1 = ERRP1 ; MAXP2 = ERRP2 ; MAXP3 = ERRP3 ; MAXV1 = SDC1 ; MAXV2 = SDC2 ; MAXV3 = SDC3 ; * 'SAUT' 'PAGE' ; 'MESS' ' ERREURS RELATIVES ' ; 'MESS' ' cas test H V' ; 'MESS' ' numero1 ' maxp1 ' ' maxv1 ; 'MESS' ' numero2 ' maxp2 ' ' maxv2 ; 'MESS' ' numero3 ' maxp3 ' ' maxv3 ; * EPS0 = 1.D-8 ; LOG4 = MAXP1 > EPS0 ; LOG5 = MAXP2 > EPS0 ; LOG6 = MAXP3 > EPS0 ; LOG7 = MAXV1 > EPS0 ; LOG8 = MAXV2 > EPS0 ; LOG9 = MAXV3 > EPS0 ; LP0 = LOG4 'OU' LOG5 'OU' LOG6 ; LV0 = LOG7 'OU' LOG8 'OU' LOG9 ; L0 = LP0 'OU' LV0 ; 'SI' ( L0 ) ; 'SINO' ; 'FINSI' ; * 'FIN' ;
© Cast3M 2003 - Tous droits réservés.
Mentions légales