Accueil | Cast3M

*
************************************************************************
************************************************************************
* ancien nom de fichier :  test_sif_3d.dgibi
 
 
************************************************************************
* Verification & Validation des procedures SIF et G_THETA
************************************************************************
* Geometrie 3D : 
* -------------
* fissure en forme de disque de rayon l2 (penny-shaped crack)
* dans un cylindre de rayon l1 et hauteur 2*l1 (milieu quasi-infini : l1>>l2)
* soumise a une contrainte sig, dans la direction orthogonale a la fissure (traction)
*
* Validation : 
* -----------
* Test sur le calcul du KI avec SIF et G_THETA options J et DECOUPLAGE
* Valeur theorique K_I = 2 * sig * (l2 / pi)**0.5
*
************************************************************************
 
************************************************************************
* OPTIONS ET PARAMETRES
************************************************************************
 
* Options generales
OPTI 'DIME' 3 'ELEM' 'CU20' ;
* GRAPH   = VRAI;  OPTI 'TRAC' 'PSC' 'EPTR' 8 'POTR' 'HELVETICA_16';
GRAPH   = FAUX;
COMPLET = FAUX;
 
* Parametres
l1  = 10.   ;
l2  = 1.    ;
sig = 42.E6 ;
 
 
************************************************************************
* MAILLAGE 
************************************************************************
 
* --densites (loin et pres de la fissure)
den1 = l1 / 2. ;
den2 = l2 / 6. ;
den3 = l2 / 12. ;
 
* --fissure levre superieure
OPTI 'DENS' den2 ;
p1  = 0. 0. 0. ;
OPTI 'DENS' den3 ;
p2  = l2 0. 0. ;
l12 = DROI p1 p2 'DINI' den2 'DFIN' den3;
s1  = l12 ROTA 90 p1 (p1 PLUS (0. 0. 1.)) ;
lff = s1 COTE 2 ;
ELIM s1 1.E-10 ;
s1s = REGE s1 ;
s1i = s1s PLUS (0. 0. 0.) ;
ELIM s1i lff 1.E-10 ;
 
* --section non fisuree
OPTI 'DENS' den3 ;
nbel3 = 4; 
p3  = (l2+(nbel3*den3)) 0. 0. ;
p4  = 0. (l2+(nbel3*den3)) 0. ;
l34 = CERC (NBEL lff) p3 p1 p4 ;
s2a  = REGL lff nbel3 l34;
OPTI 'DENS' den1 ;
p5  = l1 0. 0. ;
p6  = 0. l1 0. ;
l56 = CERC (NBEL lff) p5 p1 p6 ;
s2b  = REGL l34 l56 'DINI' den3 'DFIN' den1;
s2 = s2a et s2b;
 
* --volumes superieurs et inferieurs
v1   = VOLU 'TRAN' (s1s ET s2) 'DINI' den3 'DFIN' den2 (0. 0. (0.05 * l1 ));
v1   = VOLU 'TRAN' v1          'DINI' den3 'DFIN' den1 (0. 0. (0.95 * l1 ));
v2   = VOLU 'TRAN' (s1i ET s2) 'DINI' den3 'DFIN' den2 (0. 0. (-0.05 * l1));
v2   = VOLU 'TRAN' v2          'DINI' den3 'DFIN' den1 (0. 0. (-0.95 * l1));
vtot = v1 ET v2 ;
etot = ENVE vtot ;
atot = ARET etot ;
 
* trace et information sur le maillage
MESS 'nbel='(NBEL vtot) ' nbno=' (NBNO vtot);
SI GRAPH; 
  trac ((coul s1s ORAN) et s2) 'TITRE' 'Maillage en z=0';
  trac (vtot et (coul ORAN s1s))  'TITRE' 'Maillage (la fissure en orange)';
  trac vtot 'CACH' 'TITRE' 'Maillage';
FINSI;
 
 
************************************************************************
* MODELE
************************************************************************
 
* Modele et materiau
mo = MODE vtot 'MECANIQUE' ;
you = 200.E9;
nu = 0.3;
ma = MATE mo   'YOUN' you 'NU' nu;
 
* Blocages
xtot ytot ztot = COOR vtot ;
xtol = 1.E-8;
px0   = etot POIN 'PLAN' p1 p4 (0. 0. 1.) xtol ;
py0   = etot POIN 'PLAN' p1 p3 (0. 0. 1.) xtol ;
*px0   = xtot POIN 'MINI' ;
*py0   = ytot POIN 'MINI' ;
bl1   = BLOQ 'UX' px0 ;
bl2   = BLOQ 'UY' py0 ;
bl3   = BLOQ 'UZ' l56 ;
 
* Chargement
shaut = v1 FACE 2 ;
sbas  = v2 FACE 2 ;
f1    = PRES 'MASS' mo shaut (-1. * sig) ;
f2    = PRES 'MASS' mo sbas  (-1. * sig) ;
f     = f1 ET f2 ;
SI GRAPH; 
  vec_f = VECT f 'FORC' 'ROUG';
  trac vec_f vtot 'CACH' 'TITRE' 'Maillage';
FINSI;
 
************************************************************************
* RESOLUTION
************************************************************************
 
* Resolution lineaire
r0  = RIGI mo ma ;
rig = r0 ET bl1 ET bl2 ET bl3 ;
u   = RESO rig f ;
 
************************************************************************
* POST TRAITEMENT
************************************************************************
 
* post-traitement graphique
SI GRAPH; 
    svm1 = VMIS (sigm u mo ma) mo ma;
    def1 = defo u vtot svm1 mo;
    izo1 = prog 4.E7 PAS 0.5E7 1.E8;
    trac def1 izo1; 
    trac def1 'BOITE' (s1s et s2a) izo1; 
FINSI;
 
* Solution de reference
k1ref = 2. * sig * ((l2 / pi) ** 0.5) ;
 
* Solution par SIF
ts = TABL ;
ts . 'MAILLAGE' = vtot ;
ts . 'FRTFISS'  = lff  ;
ts . 'LEVRE_1'  = s1s  ;
ts . 'MODMIXTE' = VRAI ;
ts . 'LEVRE_2'  = s1i  ;
SIF ma u ts ;
k1s  = ts . 'K1' . p2 ;
errs = 100. * (k1s - k1ref) / k1ref ;
 
* Solution par G_THETA DECOUPLAGE
tk = TABL ;
tk . 'OBJECTIF' = MOT 'DECOUPLAGE' ;
tk . 'COUCHE' =  3 ;
tk . 'FRONT_FISSURE' = lff ;
tk . 'LEVRE_SUPERIEURE' = s1s ;
tk . 'LEVRE_INFERIEURE' = s1i ;
tk . 'SOLUTION_RESO' = u ;
tk . 'MODELE' = mo ;
tk . 'CARACTERISTIQUES' = ma ;
tk . 'CHARGEMENTS_MECANIQUES' = f ;
tk . 'BLOCAGES_MECANIQUES' = bl1 ET bl2 ET bl3 ;
 
G_THETA tk ;
 
k1  = tk . 'RESULTATS' . 'I' . 'GLOBAL' ;
errk1 = 100. * (k1 - k1ref) / k1ref ;
* trace de K1 le long du front 
evk1 = EXTR (EVOL ROSE 'CHPO' tk . 'CHPO_RESULTATS' 'K1' lff) 'COUR' 1;
* dess evk1 'YBOR' 0 80.E6 'YGRA' 10.E6;
 
* Solution par G_THETA J
tg = TABL ;
tg . 'OBJECTIF' = MOT 'J' ;
tg . 'COUCHE' =  3 ;
tg . 'FRONT_FISSURE' = lff ;
tg . 'LEVRE_SUPERIEURE' = s1s ;
tg . 'LEVRE_INFERIEURE' = s1i ;
tg . 'SOLUTION_RESO' = u ;
tg . 'MODELE' = mo ;
tg . 'CARACTERISTIQUES' = ma ;
tg . 'CHARGEMENTS_MECANIQUES' = f ;
tg . 'BLOCAGES_MECANIQUES' = bl1 ET bl2 ET bl3 ;
 
G_THETA tg ;
 
g1  = tg . 'RESULTATS' . 'GLOBAL' ;
Estar = you / (1. - (nu**2));
k1g = (Estar * g1)**0.5;
errk1g = 100. * (k1g - k1ref) / k1ref ;
* trace de K1g le long du front 
evk1g = EXTR (EVOL AZUR 'CHPO' tg . 'CHPO_RESULTATS' 'J' lff) 'COUR' 1;
evk1g = (Estar * evk1g)**0.5; 
SI GRAPH; 
* dess evk1g 'YBOR' 0 80.E6 'YGRA' 10.E6;
    tdess = TABL;
    tdess . 'TITRE' = tabl;
    tdess . 'TITRE' . 1 = mot 'G-\q : K';
    tdess . 'TITRE' . 2 = mot 'G-\q : J';
    dess (evk1 et evk1g) 'TITY' 'K_{I}'
    'YBOR' 0 80.E6 'YGRA' 10.E6 'LEGE' 'SO' tdess;
FINSI;
 
* verification de l'independace du contour pour G
SI COMPLET ;
    g1prog = prog;
    ng = 6;
    repe Bg ng;
        tg . 'COUCHE' =  &Bg ;
        G_THETA tg ;
        g1prog  = g1prog et tg . 'RESULTATS' . 'GLOBAL' ;
    fin Bg;
    k1prog = (Estar * g1prog)**0.5;
* evGcon = EVOL 'AZUR' 'MANU' 'LEGE' 'G-\t : J' 
*           'contour' (prog 1 pas 1 ng)  'G' g1prog;
* DESS evGcon XBOR 1 ng XGRA 1 ;
    evK1con = EVOL 'AZUR' 'MANU' 'LEGE' 'G-\t : J' 
            'contour' (prog 1 pas 1 ng)  'K_{I}' k1prog;
    SI GRAPH; 
    DESS evK1con XBOR 1 ng XGRA 1 ; 
    FINSI;
FINSI; 
 
************************************************************************
* TEST DE BON FONCTIONNEMENT
************************************************************************
 
* Affichage
OPTI 'ECHO' 0 ;
SAUT 3 'LIGNE' ;
MESS 'BILAN, CALCUL DE KI' ;
MESS ;
MESS 'Solution de reference                       :' ' ' k1ref  ;
MESS 'Calcul avec SIF                + erreur (%) :' ' ' k1s errs ;
MESS 'Calcul avec G_THETA decouplage + erreur (%) :' ' ' k1 errk1 ;
MESS 'Calcul avec G_THETA J          + erreur (%) :' ' ' k1g errk1g ;
 
* opti echo 1 donn 5 TRAC 'X';
 
* Test d'erreur
logerr = FAUX ;
SI ((ABS errs) > 4.E0) ;
  MESS '--> Erreur dans le calcul du KI avec SIF'  ;
  logerr = VRAI ;
FINSI ;
SI ((ABS errk1) > 5.E0) ;
  MESS '--> Erreur dans le calcul du KI avec G_THETA'  ;
  logerr = VRAI ;
FINSI ;
SI logerr ;
  ERRE 'ECHEC  DU CAS TEST !' ;
SINON ;
  MESS 'SUCCES DU CAS TEST !' ;
FINSI ;
 
 
 
FIN ;