Accueil | Cast3M

* fichier :  rupt29.dgibi
************************************************************************
************************************************************************
opti trac psc eptr 5 ;
 
GRAPH = 'O' ;
 
*******************************************************
*           Test rupt29.dgibi: Jeux de données        *
*           ---------------------------------         *
*******************************************************
 
*Cas test de validation pour le calcul de J sous plusieurs chargement 
*avec les procedures g_theta.procedur et g_calcul.procedur
*
*- chargement avec pression sur levres avec plusieurs types d elements
*  sur le maillage de la fissure
 
*Calcul en dimension 3 avec des elements CU20 et PRI15 sur un maillage 
*complet non symetrique
 
 
* DONNEES GEOMETRIQUES
*-------------------- 
 opti dime 3 elem cu20 echo 1;
 zz = 1e-6 ;
 
 aa = 2. ; Rsa = 5. ;
*
 R = Rsa * aa ;
 dda = aa / 40. ;
*
* ----- MAILLAGE
*
 pfond = aa 0. 0. ;
 np = 4 ; nt = 6 ;
 p1 = pfond plus ((np*dda) 0. 0.) ;
 p2 = p1 plus (0. 0. (np*dda)) ;
 p4 = pfond moin ((np*dda) 0. 0.) ;
 p3 = p4 plus (0. 0. (np*dda)) ;
 ss1 = (p1 d np pfond d np p4) tran np (0. 0. (np*dda)) ;
 p5 = pfond plus ((aa / 2.) 0. 0.) ;
 p6 = p5 plus (0. 0. (aa / 2.)) ;
 p8 = pfond moin ((aa / 2.) 0. 0.) ;
 p7 = p8 plus (0. 0. (aa / 2.)) ;
 ll1 = p1 d np p2 d (2*np) p3 d np p4 ;
 ll2 = p5 d np p6 d (2*np) p7 d np p8 ;
 ss2 = ll1 regl ll2 dini dda dfin (aa/8.) ;
 
 dx = R - (1.5 * aa) ;
 ss3 = (p5 d np p6) tran (dx 0. 0.) dini (aa/4.) dfin aa ;
 ss4 = ss1 et ss2 et ss3 ; elim zz ss4 ;
 vv1 = ss4 volu nt rota 90. (0. 0. 0.) (0. 0. 1.) ;
 
 pc = 0. 0. 0. ;
 p9 = 0. (aa / 2.) 0. ;
 nq = 4 ;
 ll3 = (p8 c nt pc p9) et (p9 d nq pc d nq p8) ;
 ss5 = surf ll3 plane ;
 vv2 = ss5 volu tran np (0. 0. (aa / 2.)) ;
 elim zz (vv1 et vv2) ;
 
 ss6 = (vv1 et vv2) poin plan p6 p7 (0. 1. (aa / 2.)) zz ;
 ss7 = (enve (vv1 et vv2)) elem appuye strictement ss6 ;
 dz = R - (0.5 * aa) ;
 vv3 = ss7 volu tran (0. 0. dz) dini (aa/4.) dfin aa ;
 
 vtot1 = vv1 et vv2 et vv3 ;
 
*---------- Surface de la fissure
 ss0 = vtot1 poin plan pc p8 p9 zz ;
 ss1 = (enve vtot1) elem appuye strictement ss0 ;
 chx = ss1 coor 1 ; chy = ss1 coor 2 ;
 chr = ((chx ** 2) + (chy ** 2)) ** 0.5 ;
 chrr = chan cham chr ss1 ;
 sulev1 = chrr elem infe (aa - zz) ;
 sulig = diff ss1 sulev1 ;
 ffis = (cont sulev1) elem appuye strictement sulig ;
 
*---------- Fin du modèle
 
 vtot2 = vtot1 syme plan pc p8 p9 ;
 sulev2 = sulev1 syme plan pc p8 p9 ;
 sulev2 = INVE sulev2 ;
 elim zz (vtot2 et sulev2 et sulig) ;
 
 vtot = vtot1 et vtot2 ;
 
*---------- Lignes et surfaces nommées
 sulev = sulev1 ET sulev2 ;
 
 ss0 = vtot poin plan pc p8 p7 zz ;
 symy = (enve vtot) elem appuye strictement ss0 ;
 
 ss0 = vtot poin plan pc p9 (0. 0. 1.) zz ;
 symx = (enve vtot) elem appuye strictement ss0 ;
 
 ss0 = vtot point plan (0. 0. R) (1. 0. R) (0. 1. R) zz ;
 ssup = (enve vtot) elem appuye strictement ss0 ;
 Rp = -1. * R ;
 ss0 = vtot point plan (0. 0. Rp) (1. 0. Rp) (0. 1. Rp) zz ;
 sinf = (enve vtot) elem appuye strictement ss0 ;
 
si (ega GRAPH 'O') ;
  trac vtot qual cach ;
fins ;
 
******************************************************************
 
*-- Modele mecanique et parametres materiau
EI20 = 197000. ; ALI20 = 15.54E-6 ;
mo2 = MODE vtot MECANIQUE ELASTIQUE ISOTROPE ;
ma2 = MATE mo2 'YOUN' EI20 'NU' 0.3 'ALPH' AlI20 ;
ri1 = RIGI mo2 ma2 ;
 
*--- Conditions aux limites
* Blocages
cd1 = BLOQ symx  'UX' ;
cd2 = BLOQ symy  'UY' ;
cd3 = BLOQ sinf  'UZ' ;
cdl = cd1 ET cd2 ET cd3 ;
*Traction uniaxiale (via un modele de pression)
pp = 100. ;
moph = MODE ssup 'CHARGEMENT' 'PRESSION' 'CONS' 'HAUT' ;
*Pression sur les levres (via un modele de pression)
mopl = MODE sulev 'CHARGEMENT' 'PRESSION' 'CONS' 'LEVRES' ;
 
* Construction des second membres
maph  = PRES moph 'PRES' pp ;
effN  = BSIG moph maph ;
mapl1 = PRES mopl 'PRES' pp ;
ff0   = BSIG mopl mapl1 ;
zer1  = ZERO mopl 'RIGIDITE' ;
chx chy chz = COOR zer1 ;
chr = (((chx ** 2) + (chy ** 2)) ** 0.5) / aa ;
*ch2 = pp * ((1. - (chr ** 2)) ** 0.5) ;
ch2 = pp * (chr ** 2) ;
mapl2 = PRES mopl ch2 ;
ff2  = BSIG mopl mapl2 ;
 
 
*--- 11 - RESOLUTION
depN dep0 dep2 = RESO (ri1 ET cdl) effN ff0 ff2 ;
 
 
*--- 12 - Calcul de J pour les trois conditions aux limites
*----------------------------------------------------------
EE = EI20 / 910. ; rpa = (3.1416 * aa / 1000.) ** 0.5 ;
* cas 1 : traction seule
SUPTAB = TABLE;
SUPTAB.'MODELE' = mo2 ET moph ;
SUPTAB.'CARACTERISTIQUES' = ma2 ;
SUPTAB.'PRESSION'         = maph ;
SUPTAB.'BLOCAGES_MECANIQUES' = cdl ;
SUPTAB.'OBJECTIF' = 'J' ;
SUPTAB.'LEVRE_SUPERIEURE' = sulev1 ;
SUPTAB.'LEVRE_INFERIEURE' = sulev2 ;
SUPTAB.'COUCHE' = 4 ;
SUPTAB.'FRONT_FISSURE' = ffis ;
SUPTAB.'SOLUTION_RESO' = depN ;
G_THETA SUPTAB ;
Jc0 = suptab.resultats.global ;
kc0 = (Jc0 * EE) ** 0.5 ;
ic0 = kc0 / rpa / pp ;
 
* cas 2 : pression uniforme sur les levres
SUPTAB.'MODELE' = mo2 ET mopl ;
SUPTAB.'CARACTERISTIQUES' = ma2 ;
SUPTAB.'PRESSION'         = mapl1 ;
SUPTAB.'SOLUTION_RESO' = dep0 ;
G_THETA SUPTAB ;
Jc1 = suptab.resultats.global ;
kc1 = (Jc1 * EE) ** 0.5 ;
ic1 = kc1 / rpa / pp ;
 
* cas 3 : pression parabolique sur les levres
SUPTAB.'MODELE' = mo2 ET mopl ;
SUPTAB.'CARACTERISTIQUES' = ma2 ;
SUPTAB.'PRESSION'         = mapl2 ;
SUPTAB.'SOLUTION_RESO' = dep2 ;
G_THETA SUPTAB ;
Jc2 = suptab.resultats.global ;
kc2 = (Jc2 * EE) ** 0.5 ;
ic2 = kc2 / rpa / pp ;
 
 
 
 
********************
* CALCUL DE L'ERREUR
********************
**---> Les 2 résultats doivent être identiques
**---> et converger vers la solution d'Irwin : 2/Pi
**---> Le résultat ic2 doit lui converger vers 4/(3.Pi)
err0 = ABS (100. * ((ic0 * pi / 2.) - 1.)) ;
err1 = ABS (100. * ((ic1 * pi / 2.) - 1.)) ;
err2 = ABS (100. * ((ic2 * 3. * pi / 4.) - 1.)) ;
 
 
****************************************
* AFFICHAGE DES RESULTATS ET DES ERREURS
****************************************
SAUT 5 'LIGNE' ;
 
mess 'Erreurs en % : ' err0 err1 err2 ;
mess ;
 
SI (err0 > 0.4) ;
  ERRE 'Erreur sur le calcul de Jc0' ;
FINS ;
SI (err1 > 0.4) ;
  ERRE 'Erreur sur le calcul de Jc1' ;
FINS ;
SI (err2 > 0.5) ;
  ERRE 'Erreur sur le calcul de Jc2' ;
FINS ;
 
 
FIN;