* fichier :  rupt10.dgibi
************************************************************************
* Section : Mecanique Endommagement
************************************************************************

****************************************************
*                                                  *
*    VALIDATION DE LA METHODE DES DFEPLACEMENTS    *
*    DANS LE CAS D'UNE PLAQUE EN FLEXION PURE.     *
*            SOLUTION DE REFERENCE :               *
*     Compendium of STRESS INTENSITY FACTORS       *
*      by Rooke & Cartwright. Section Plates       *
*   and Shells "Central crack in a finite width    *
*         plate : uniform bending moment"          *
*                                                  *
*                                                  *
*  GEOMETRIE :                                     *
*                                                  *
*  longueur de la fissure : 2a = 2 x 17.5 MM       *
*  largeur de la plaque   : 2b = 2 x 70 MM         *
*  Hauteur de la plaque   : 2h = 2 x 280 MM        *
*  Epaisseur de la plaque : e  = 0.7    MM         *
*                                                  *
*  MATERIAU :                                      *
*                                                  *
*  Module d'Young         : E = 200000 MPa         *
*  Coefficient de poisson : NU = 0.3               *
*                                                  *
*  CHARGEMENT :                                    *
*                                                  *
*  Flexion pure           : M = 114.3 N.MM         *
*                      soit SIGf=10MPa             *
*                                                  *
*  SOLUTION ANALYTIQUE :                           *
*                                                  *
*  Fact d'int. de contr : K1 = 30.349 MPa(MM)**0.5 *
*                                                  *
****************************************************


opti echo 1 dime 2 elem TRI3 ;

*****************************************
*  SAISSIE DES PARAMETRES DE MAILLAGE   *
*****************************************
*
* Epaisseur de l'eprouvette

epa1 = .7;

* ---------------- A/W :

rap1 = 0.25 ;

* ---------------- H/W :

rap2 = 4. ;

* ---------------- W :

w = 70. ;

* largeur de bande deposee (pour un metal biphase)

ldep = 10. ;

*eprouvette SENB ou CCP ? (val = 0  ou  val = 1) ;

val = 1 ;

*taille de maille en pointe de fissure

tm = .2 ;

*nombre de decouge sur 45 degres ;

n = 4 ;

*introduction d'un pave de quadrangle en pointe de fissure ?
*(non : rep = 0     oui : rep = 1)

rep = 1 ;

***************************************
******* parametres intermediaires *****
***************************************

lcritw = (w * (1. - rap1)) ;
lcrita = w * rap1 ;
*demi-cote du carre
c = mini (prog (w / 5.) ldep lcrita lcritw) ;

*rayon du disque
r = c / 2. ;

*densite en peripherie de disque
*te = r / 10. ;
te = 1.5 * ((pi * r) / (4. * n)) ;

*densite au premier arc
si (ega rep 1) ;
 ta = tm * 2. ;
sinon ;
 ta = tm ;
finsi ;


*******************************************
*  MAILLAGE                               *
*******************************************

*----------------------------------------------------*
*** partie circulaire autour de la pointe de fissure *
*----------------------------------------------------*

p0 = 0. 0. ;
p1 = (rap1 * w) 0. ;
dens ta ;
si (rep ega 0) ;
 p2 = p1 plus (tm 0.) ;
sinon ;
 p2 = p1 plus ((tm * n * 2.) 0.) ;
 si ((tm * n * 2.) > r) ;
  nc =  r / 2. ;
  opti echo 0 ;
  mess 'Le produit (tm * n) doit etre inferieur a' nc ;
  opti echo 1 ;
  opti donn 5 ;
 finsi ;
finsi ;
dens te ;
p10 = p1 plus (r 0.) ;

ligam0 = p2 d p10 ;

su21 = ligam0 rota n 45. p1 ;
l21 = su21 cote 4 ;
c21 =  su21 cote 2 ;
su22 = (su21 cote 3) rota n 45. p1 ;
l22 = su22 cote 2 ;
c22 =  su22 cote 4 ;
su23 = (su22 cote 3) rota n 45. p1 ;
l23 = su23 cote 4 ;
c23 =  su23 cote 2 ;
su24 = (su23 cote 3) rota n 45. p1 ;
l24 = su24 cote 2 ;
c24 =  su24 cote 4 ;

rac = l21 et l22 et l23 et l24 ;
su2 = su21 et su22 et su23 et su24 ;

si (ega rep 0) ;
 su0 = cout p1 rac ;
 su0l = chan su0 ligne ;
 p01 = (p1 plus p2) / 2. ;
 lig0 =  su0l elem appuye strictement (p1 et (su0 poin proc p01)
   et p2)    ;
 ligam0 =  ligam0 et lig0 ;
 l1 = ligam0 ;
sinon ;
 dens tm ;
 qad = p1 plus ((tm * n) 0.) ;
 qag = p1 moin ((tm * n) 0.) ;
 basd = p1 d n qad ;
 basg = p1 d n qag ;
 bas = ordo (basg et basd) ;
 pave = bas trans n (0. (tm * n)) ;
 copav1 = pave cote 2 ;
 copav2 = pave cote 3 ;
 copav3 = pave cote 4 ;
 copav = copav1 et copav2 et copav3 ;
 dec = (enti (n / 1.2)) * (-1) ;
 joncd = qad d  p2 ;
 joncg = qag d  (su24 poin proc qag) ;
 sujon = dall copav joncd rac joncg plan ;
 su0 = sujon et pave ;
 lig0 = basd et joncd ;
 ligam0 =  ligam0 et lig0 ;
 l1 = ligam0 ;
finsi ;


*---------------------------------------------*
*** carre exterieur autour du fond de fissure *
*---------------------------------------------*

dens (2. * te) ;
pdb = p1 plus (c 0.) ;
pdh = p1 plus (c c) ;
pmh = p1 plus (0. c) ;
pgh = p1 plus (((-1) * c) c) ;
pgb = p1 plus (((-1) * c) 0.) ;

card = pdb d n pdh ;
carhd = pdh d n pmh ;
carhg = pmh d n pgh ;
carg = pgh d n pgb ;

decc = ((-1) * (enti ((c - r) / (1. * te)))) - 1 ;

diad = (c21 poin initial) d decc pdb ;
diadh = (c22 poin final) d decc pdh ;
diam = (c23 poin initial) d decc pmh ;
diagh = (c24 poin final) d decc pgh ;
diag = (c24 poin initial) d decc pgb ;

su3 = dall c21 diadh card diad plan;
su4 = (dall c22 diam carhd diadh plan) et
   (dall c23 diagh carhg diam plan) ;
su5 = dall c24 diag carg diagh plan ;

carre = su0 et su2 et su3 et su4 et su5 ;
ligam1 = ligam0 et diad ;


*-------------------------*
*** partie metal depose   *
*-------------------------*


clcritw =  (c < (lcritw + 0.01)) et (c > (lcritw - 0.01)) ;
clcrita =  (c < (lcrita + 0.01)) et (c > (lcrita - 0.01)) ;
cldep =    (c < (ldep + 0.01)) et (c >  (ldep - 0.01)) ;

si ((non clcritw ) et (non clcrita)) ;
* mess 'cas general' ;

 vtl = 0. (ldep - c) ;
 vtg = (((-1) * (lcrita - c)) 0.) ;
 vtd =  (lcritw - c) 0. ;
 nfl = enti ((ldep - c) / (mesu (carhd elem 1) long)) ;
 nfg = enti ((lcrita - c) / (mesu (carg elem 1) long)) ;
 nfd = enti ((lcritw - c) / (mesu (card elem 1) long)) ;
 si (nfg ega 0) ;
  nfg = 1 ;
 su7 = carg tran nfg vtg ;
 sinon ;
  densi = (mesu (carg elem 1) long) ;
  densf = densi * (1 + (nfg / 20.)) ;
  su7 = carg tran  'DINI' densi 'DFIN' densf vtg ;
 finsi ;
 si (nfd ega 0) ;
  nfd = 1 ;
  su6 = card tran nfd vtd ;
 sinon ;
  densi = (mesu (card elem 1) long) ;
  densf = densi * (1 + (nfd / 20.)) ;
  su6 = card tran  'DINI' densi 'DFIN' densf vtd ;
 finsi ;
 l8 = (su7 cote 4) et (su6 cote 2) et carhd et carhg ;
 l8 = ordo l8 ;
 si (cldep) ;
*  mess 'cas critique ldep' ;
  mdep = carre et su6 et su7 ;
  l8b = inve l8 ;
  ligam = ligam1 et (su6 cote 4) ;
  lsym1 = (su7 cote 3) ;
 sinon ;
  si (nfl ega 0) ;
   nfl = 1 ;
   su8 = l8 tran nfl vtl ;
  sinon ;
   densi = (mesu (carhd elem 1) long) ;
   densf = densi * (1 + (nfl / 10.)) ;
   su8 = l8 tran  'DINI' densi 'DFIN' densf vtl ;
  finsi ;
  mdep = carre et su6 et su7 et su8 ;
  l8b = su8 cote 3 ;
  ligam = ligam1 et (su6 cote 4) ;
  lsym1 =  (su7 cote 3) et (su8 cote 4) ;
 finsi ;

finsi ;

si (clcritw) ;
* mess 'cas critique W' ;
 vtl = 0. (ldep - c) ;
 vtt = (((-1) * (lcrita - c)) 0.) ;
 nfl = enti ((ldep - c) / (mesu (carhd elem 1) long)) ;
 nft = enti ((lcrita - c) / (mesu (carg elem 1) long)) ;
 si (nft ega 0) ;
  nft = 1 ;
  su7 = carg tran nft vtt ;
 sinon ;
  densi = (mesu (carg elem 1) long) ;
  densf = densi * (1 + (nft / 20.)) ;
  su7 = carg tran  'DINI' densi 'DFIN' densf vtt ;
 finsi ;
 l8 = (su7 cote 4) et carhd et carhg ;
 l8 = ordo l8 ;
 si (cldep) ;
*  mess 'cas critique ldep et lcritw' ;
  mdep = carre et su7 ;
  l8b = inve l8 ;
  ligam = ligam1 ;
  lsym1 = (su7 cote 3) ;
 sinon ;
  si (nfl ega 0) ;
   nfl = 1 ;
   su8 = l8 tran nfl vtl ;
  sinon ;
   densi = (mesu (carhd elem 1) long) ;
   densf = densi * (1 + (nfl / 10.)) ;
   su8 = l8 tran  'DINI' densi 'DFIN' densf vtl ;
  finsi ;
  mdep = carre et su7 et su8 ;
  l8b = su8 cote 3 ;
  ligam = ligam1 ;
  lsym1 =  (su7 cote 3) et (su8 cote 4) ;
 finsi ;
finsi ;

si (clcrita) ;
* mess 'cas critique A' ;
 vtl = 0. (ldep - c) ;
 vtt = (lcritw - c) 0. ;
 nfl = enti ((ldep - c) / (mesu (carhd elem 1) long)) ;
 nft = enti ((lcritw - c) / (mesu (card elem 1) long)) ;
 si (nft ega 0) ;
  nft = 1 ;
  su6 = card tran nft vtt ;
 sinon ;
  densi = (mesu (card elem 1) long) ;
  densf = densi * (1 + (nft / 20.)) ;
  su6 = card tran  'DINI' densi 'DFIN' densf vtt ;
 finsi ;
 l8 = (su6 cote 2) et carhd et carhg ;
 l8 = ordo l8 ;
 si (cldep) ;
*  mess 'cas critique ldep et lcrita' ;
  mdep = carre et su6 ;
  l8b = l8 ;
  ligam = ligam1  et (su6 cote 4) ; ;
  lsym1 = carg ;
 sinon ;
  si (nfl ega 0) ;
   nfl = 1 ;
   su8 = l8 tran nfl vtl ;
  sinon ;
   densi = (mesu (carhd elem 1) long) ;
   densf = densi * (1 + (nfl / 10.)) ;
   su8 = l8 tran  'DINI' densi 'DFIN' densf vtl ;
  finsi ;
  mdep = carre et su6 et su8 ;
  l8b = su8 cote 3 ;
  l8b = inve l8b ;
  ligam = ligam1 et (su6 cote 4) ;
  lsym1 = (su8 cote 2) et carg ;
 finsi ;
finsi ;

*-------------------------*
*** partie metal de base  *
*-------------------------*

l8b = ordo l8b ;
ppi = l8b poin initial ;pf = l8b poin final ;
lu = (mesu l8b long) / (nbel l8b) * 1.2 ;
fron1 = (ppi plus (0. lu)) d (enti ((nbel l8b) / 1.5))
   (pf plus (0. lu)) ;
lub = ((mesu fron1 long) / (nbel fron1)) * 1.2 ;
pib = fron1 poin initial ;pfb = fron1 poin final ;
fron2 =  (pib plus (0. lub)) d (enti ((nbel fron1) / 1.5))
   (pfb plus (0. lub)) ;

join1 = cout l8b fron1 ;
join2 = cout fron1 fron2 ;
join = join1 et join2 ;
joinl = chan join ligne ;
joinp = joinl poin droit p0 (p0 plus (0. 10.)) (tm / 10.) ;
lsym2 = joinl elem appuye strictement joinp ;

lrest = ((w * rap2) / 2.) - (ldep + lu + lub) ;
vt = 0. lrest ;
le = mesu (fron2 elem 1) long ;
corp = fron2 tran (enti (lrest / le)) vt ;
lsym3 = corp cote 2 ;

pinif = ligam poin initial;

*---------------------------*
*particularite pour un CCP  *
*---------------------------*

si (val ega 1) ;
 mbas = corp et join ;
*plan verticale de symetrie
 lsym = lsym1 et lsym2 et lsym3 ;
 lsym = lsym coul vert ;
*face de chargement
 schar = corp cote 3 ;
 schar = schar coul roug ;
*ligament
 ligam = ligam coul bleu ;
*maillage complet
 ccp = mbas et mdep ;
 titr 'eprouvette CCP' ;
* trac (ccp et ligam et lsym et schar) ;
finsi ;

*----------------------------*
*particularite pour un SENB  *
*----------------------------*

si (val ega 0) ;
 lresid = corp cote 3 ;
 resid = lresid tran 1 (0. le) ;
 mbas = corp et join et resid ;
*appui
 pap = p0 plus (0. ((rap2 / 2.) * w)) ;
 pap = mbas poin proc pap ;
 pap = manu pap poi1 vert ;
*point d'application de la force
 pchar = p0 plus (w 0.) ;
 pchar = mdep poin proc pchar ;
 pchar = manu pchar poi1 roug ;
*ligament
 ligam = ligam coul bleu ;
*maillage complet
 senb = mbas et mdep ;
 titr 'eprouvette SENB' ;
 trac (senb et ligam et pap et pchar) ;
finsi ;
*
*trac (cont ccp);
*
*  FIN DE LA DEFINITION DE LA GEOMETRIE
*
*
opti dime 3 mode trid;
*
*  Definition des modeles et du materiau
*
MO1 = MODELE CCP MECANIQUE ELASTIQUE COQ3;
MA0 = MATER MO1 YOUNG 2.E11 NU 0.3 EPAI epa1;
RIG0 = RIGI MO1 MA0;
*
*  Definition des conditions aux limites
*
BL1 = BLOQ UY RX ligam;
BL2 = BLOQ UX RY lsym;
BL3 = BLOQ UZ schar;
RIG1 = RIG0 ET BL1 ET BL2 ET BL3;
*
*  Chargement : Moment de flexion
*
pinich = schar poin initial;
pfinch = schar poin final;
schar=schar chan poi1;
PMILCH = schar diff (pinich et pfinch);
MOEXT = 10.*(w * (epa1** 2) / 3. / 2.) / ((NBNO schar) - 1);
MOMIL = MOEXT * (NBNO PMILCH);
FO1 = (MOME MX MOEXT (PINICH ET PFINCH)) ET (MOME MX MOMIL PMILCH);
*
*  Resolution
*
DEP1 = RESO RIG1 FO1;
*
*  Calcul du facteur d'intensite de contrainte
*
pfisin = lsym poin init;
ct1 = cont ccp;
psup = (ct1 diff (ligam et lsym)) poin droi pfisin p1;
ct2 = chan 'SEG2' ct1;
lifis1 = ct2 elem appu stri psup;
*
SUPTAB = TABLE;
SUPTAB.'MODMIXT' = FAUX;
SUPTAB.'LEVRE_1' = lifis1;
SUPTAB.'FRTFISS' = p1;
SUPTAB.'MEMBRANE' = FAUX;
SUPTAB.'FLEXION' = VRAI;
SUPTAB.'EPAI' = EPA1;
*
K1TH = 30.3486;
*
SIF SUPTAB MA0 DEP1;
K1MOY = SUPTAB.K1.TOTAL;
MESS '   K1 analytique      = 'K1TH;
MESS '   K1 calcule par SIF = 'K1MOY;
err1 = 100.*(K1TH - K1MOY)/K1TH;err1=abs err1;
mess '   Erreur : 'err1' en %';
*
SI (err1 < 2.) ;
  MESS 'PROCEDURE <SIF> ERR  0';
  ERRE 0;
SINON;
  MESS 'PROCEDURE <SIF> ERR  5';
  ERRE 5;
FINSI;
fin;