* fichier :  rupt7.dgibi
*opti echo 0;
*******************************************************
*                                                     *
*           Test Rupt7.dgibi: Jeux de données         *
*           ---------------------------------         *
*                                                     *
*******************************************************
*****   CALCUL DU FACTEUR D'INTENSITE DE          *****
*****   CONTRAINTES D'UNE PLAQUE AVEC FISSURE     *****
*****   RECTILIGNE INCLINEE                       *****
*****   CHARGEE EN TRACTION UNIFORME              *****
*****   PAR SIF (interpolation des deplacements)  *****
*****   ET PAR G_THETA (methode integrale)        *****
*******************************************************
*****   SOLUTION DE REFERENCE DONNE PAR :         *****
*****   [TADA, STRESS ANALYSIS HANDBOOK,1973]     *****
*****  +[ISIDA,IJF(7),301-316,1971] facteur de forme **
*******************************************************
*****   demi-LARGEUR DE LA PLAQUE : b=125 mm ;    *****
*****   LONGUEUR   : 2*h=250 mm                   *****
*****   LONGUEUR DE FISSURE : 50 mm   ;           *****
*****   ANGLE / OX : 45°                          *****
*****   CONTRAINTE UNIFORME : 10 N/m2             *****
*******************************************************
*
opti dime 2 elem qua8 mode plan defo echo 0;
*graph=faux;
graph=vrai;
opti trac psc eptr 5;
*
*--------------DEFINITION DU MAILLAGE-----------------;
*
*geometrie
b = 0.125;
h = 0.125;
alpha1 = 45.;
*alpha1 = 30.;
*alpha1 = 20.;
a = 0.050;
ax = 0.050 * (sin alpha1);
ay = 0.050 * (cos alpha1);
*points
p0=0. 0.;
p11=h 0.;
p12=p11 plus (0. 0.);
p2=(2.*h) 0.;
p3=0. b;
pf=p11 plus (ax ay);
p4=(p11 plus (ax 0.)) plus (0. b);p5=(2.*h) b;
*zone pointe de fissure
i=0;
pd = (prog 1. 2.5 5. 9. 13. 19. 25.) *1.E-3 ;
cim1=pf;cim2=pf;cim3=pf;cim41=pf;cim42=pf;
plfim1=pf;plfim2=pf;pim2=pf;
repeter bbbb 7;
   i= i + 1;
   v=extr pd i;
   pi1=pf moin (v 0.);pi2=pf plus (0. v);
   pi3=pf plus (v 0.);pi4=pf moin (0. v);
*   ci4=cerc  2 pi4 pf pi1;
*   plfi1=ci4 poin 2;
   plfi1 = pf moin ((v * (sin alpha1)) (v * (cos alpha1)));
   plfi2=plfi1 plus (0. 0.);
   n=2;si ((i ega 6 ) ou (i ega 7));n=4;finsi;
   si(alpha1 < 40.); m = n ;
   sinon;            m = n / 2;
   fins;
   ci1=cerc  n pi1 pf pi2;
   ci2=cerc  n pi2 pf pi3;
   ci3=cerc  n pi3 pf pi4;
   ci41=cerc  m pi4 pf plfi1;
   ci42=cerc  m plfi2 pf pi1;
   sui1=cout cim1 ci1;
   sui2=cout cim2 ci2;
   sui3=cout cim3 ci3;
   sui41=cout cim41 ci41;
   sui42=cout cim42 ci42;
   elim 1.E-6 (sui42 et sui1);
   elim 1.E-6 (sui1 et sui2 et sui3 et sui41);
   lfi1=(cont sui41) ELEM comp plfim1 plfi1;
   lfi2=(cont sui42) ELEM comp plfim2 plfi2;
   lhi=pim2 d 1 pi2;
   cim1=ci1;cim2=ci2;cim3=ci3;cim41=ci41;cim42=ci42;
   plfim1=plfi1;plfim2=plfi2;pim2=pi2;
   sut1=sui42 et sui1;
   sut2=sui2 et sui3 et sui41;
   si (i ega 1);
     sutt1=sut1;sutt2=sut2;lf1=lfi1;lf2=lfi2;lh=lhi;
   sinon;
     sutt1=sutt1 et sut1;
     sutt2=sutt2 et sut2;
     lf1=lf1 et lfi1;lf2=lf2 et lfi2;lh=lh et lhi;
   finsi;
*   trac (sutt1 et (coul sutt2 vert));
fin bbbb;
*zone exterieure gauche
lg=p0 d 5 p3;
lhc= p4 d 3 pi2;
lig1=lg et (p3 d 8 p4) et lhc;
lig1=lig1 et (inve ci1) et ci42;
lig1=lig1 et (plfi2 d 2 p12);
nbh = 6;
lig1=lig1 et (p12 d nbh p0);
su1=surf lig1 plan;
obj1=su1 et sutt1;
elim 1.E-6 obj1;
*zone exterieure droite
ld=p2 d 5 p5;
lig2=ld et (p5 d 8 p4) et lhc;
lig2=lig2 et ci2 et ci3 et ci41;
lig2=lig2 et (plfi1 d 2 p11);
lig2=lig2 et (p11 d nbh p2);
su2=surf lig2 plan;
obj2=su2 et sutt2;
elim 1.E-6 obj2;
*recup des levres de fissures
perobj1 = inve ( cont sutt1);
lifis1 = (perobj1 ELEM comp pf plfi2) coul rose;
elfis1 = lifis1 elem 'APPU' 'LARG' pf;
elfis1 = lifis1 elem 'APPU' 'LARG' elfis1;
elfis1 = ordo elfis1;
ptfis1 = poin elfis1 2;
vfis1 = pf moin ptfis1;
vfis1 = vfis1 / (norm vfis1);
vfis1ok = pf moin plfi2;
vfis1ok=vfis1ok / (norm vfis1ok);
psc1 = abs (psca vfis1 vfis1ok);
si (psc1 < 0.99);
    lifis1 = (perobj1 ELEM comp plfi2 pf) coul rose;
fins;
*recup des levres de fissures
perobj2 = inve ( cont sutt2);
lifis2 = (perobj2 ELEM comp pf plfi1) coul vert;
elfis2 = lifis2 elem 'APPU' 'LARG' pf;
elfis2 = lifis2 elem 'APPU' 'LARG' elfis2;
elfis2 = ordo elfis2;
ptfis2 = poin elfis2 2;
vfis2 = pf moin ptfis2;
vfis2 = vfis2 / (norm vfis2);
vfis2ok = pf moin plfi2;
vfis2ok=vfis2ok / (norm vfis2ok);
psc2 = abs (psca vfis2 vfis2ok);
si (psc2 < 0.99);
    lifis2 = (perobj2 ELEM comp plfi1 pf) coul rose;
fins;
*symetrie
obj1tmp obj2tmp lgtmp ldtmp
  = syme obj1 obj2 lg ld 'DROIT' p0 p2;
obj2s obj1s lds lgs
  = syme obj1tmp obj2tmp lgtmp ldtmp 'DROIT' p11 (p11 moin (0. b));
ELIM 1.E-6 obj1s obj1;
ELIM 1.E-6 obj2s obj2;
ELIM 1.E-6 lgs obj1s;
ELIM 1.E-6 lds obj2s;
ELIM 1.E-6 lgs lg;
ELIM 1.E-6 lds ld;
*assemblage
obj = (obj1 et obj2 et obj1s et obj2s);
ld= (ld et lds) coul roug;
lg= (lg et lgs) coul bleu;

*-----------DEFINITION DU MATERIAU-------------------;
*
* objaf = modl obj MECANIQUE ELASTIQUE ISOTROPE;
phi0 = PSIP obj lifis1 ;
phi0 = CHAN 'CHAM' phi0 obj;
obj1tot = phi0 ELEM 'SUPERIEUR' 1.E-4;
obj2tot = (obj DIFF obj1tot) coul BLEU;
trac (obj1tot et obj2tot);
mess (nbel obj) ((nbel obj1tot) + (nbel obj2tot));
mod1tot = mode obj1tot MECANIQUE ELASTIQUE ISOTROPE;
mod2tot = mode obj2tot MECANIQUE ELASTIQUE ISOTROPE;
objaf =mod1tot et mod2tot;
*
mat = mate objaf youn 20000 nu 0.3;
*
*-----------MATRICE DE RIGIDITE----------------------;
*
rig=rigi mat objaf;
*
*-----------DEFINITION DU CHARGEMENT-----------------;
*
cdl=(bloq ux uy p0) et (bloq uy p2);
sxx1 = 10.;
f1=pres mass objaf lg (-1.*sxx1);
f2=pres mass objaf ld (-1.*sxx1);
*
*-----------RESOLUTION ET CONTRAINTES----------------;
*
u=reso (rig et cdl) (f1 et f2);
sig=sigm mat objaf u;

si(graph);
  vf12 = vect (f1 et f2) 'FORC';
  trac vf12 obj ;
  trac sig objaf (defo u obj 20) ;
  trac (obj et lifis1 et lifis2) (defo u obj 20);
fins;

*opti trac x;
*----CALCUL DU FACTEUR D'INTENSITE DE CONTRAINTES----;
*----APPEL DE LA PROCEDURE SIF  ---------------------;
SIFTAB = TABL;
SIFTAB.'MODMIXTE' = VRAI;
SIFTAB.'LIFIS1' = lifis2;
SIFTAB.'LIFIS2' = lifis1;
SIFTAB.'FRTFISS' = pf;
SIF SIFTAB MAT U;
K1MOY = SIFTAB.K1;
K2MOY = SIFTAB.K2;

*----APPEL DE LA PROCEDURE G_THETA  ---------------------;
GTAB = TABL;
GTAB . 'OBJECTIF' = MOT 'DECOUPLAGE';
GTAB . 'FRONT_FISSURE'  = pf;
GTAB . 'LEVRE_SUPERIEURE' = lifis1;
GTAB . 'LEVRE_INFERIEURE' = lifis2;
GTAB . 'COUCHE' = 4;
GTAB . 'SOLUTION_RESO'  = u;
GTAB . 'CARACTERISTIQUES' = mat;
GTAB . 'MODELE' = objaf;
GTAB . 'CHARGEMENTS_MECANIQUES' = (f1 et f2) ;
G_THETA GTAB;
K1G = GTAB . 'RESULTATS' . 'I';
K2G = GTAB . 'RESULTATS' . 'II';

*-----DEPOUILLEMENT  ---------------------------------;
*pour a/b=0.4 et b/h=1, on a:
Fac1 = 1.216;
*solution infinie:
K0inf = sxx1 * ((a*pi)**0.5 ) * (cos alpha1);
*solution de reference
K1ref = K0inf * (cos alpha1) * Fac1;
K2ref = K0inf * (sin alpha1) * Fac1;

*ecarts relatif (%)
EC1 = (abs (1. - (abs (K1MOY / K1REF)))) * 100;
EC2 = (abs (1. - (abs (K2MOY / K2REF)))) * 100;
EC3 = (abs (1. - (abs (K1G / K1REF)))) * 100;
EC4 = (abs (1. - (abs (K2G / K2REF)))) * 100;

*messages
mess '------------------------------------------------';
mess (chai '                   K1           K2   ');
mess (chai 'reference     ' K1ref '  ' K2ref );
mess (chai '   SIF        ' K1MOY '  ' K2MOY );
mess (chai '   G_THETA    ' K1G   '  ' K2G   );
mess '------------------------------------------------';
mess (chai 'erreur relative (%)  K1        K2   ');
mess (chai '   SIF        ' EC1 '  ' EC2 );
mess (chai '   G_THETA    ' EC3 '  ' EC4 );
mess '------------------------------------------------';

*opti donn 5;
GTAB . 'MODELES_COMPOSITES' = TABL;
GTAB . 'MODELES_COMPOSITES' . 1 = mod1tot;
GTAB . 'MODELES_COMPOSITES' . 2 = mod2tot;
G_THETA GTAB;
K1Gc = GTAB . 'RESULTATS' . 'I';
K2Gc = GTAB . 'RESULTATS' . 'II';
mess (chai '   G_THETA composite    ' K1Gc   '  ' K2Gc   );

*OPTI DONN 5 trac X;
*-----------------------------------------------------;
*------------------MESSAGE D'ERREUR-------------------*
SI ((EC1 < 3.) ET (EC2 < 9.) et (EC3 < 1.) et (EC4 < 6.));
  MESS 'ERR   0';
  ERRE  0;
SINON;
  MESS 'ERR   5';
  ERRE  5;
FINSI;
*-----------------------------------------------------*
fin ;