Télécharger xfem04.dgibi

Retour à la liste

Numérotation des lignes :

  1. optio norm auto;
  2. * fichier : xfem04.dgibi
  3. ************************************************************************
  4. ************************************************************************
  5. ************************************************************************
  6. *
  7. * Mots-clés : Mécanique de la rupture, XFEM, contact, frottement, LATIN
  8. *
  9. * Etude de reponse d'une plaque en compression avec fissure inclinee
  10. * via une formulation mixte a trois champs. Le contact frottant et la
  11. * resolution par la LATIN sont pris en compte par des procedures perso.
  12. * test des operateurs : RELA ACCRO FAIBLE, G_THETA
  13. * creation : 2013-02-13, Benoit Trolle (these SNCF-Lamcos)
  14. * integration dans Cast3m : 2013-12-17, Benoit Prabel (#)
  15. *
  16. ************************************************************************
  17.  
  18. OPTION DIME 2 MODE PLAN DEFO ELEM QUA4;
  19.  
  20.  
  21. * Masquer lignes excutées
  22. * 'OPTION' echo 0;
  23.  
  24. ************************************************************************
  25. * PROCEDURES INTERNES *
  26. ************************************************************************
  27.  
  28. *-----------------------------------------------------------------------
  29. * PLOTITER : Affichage des champs locaux au cours des itérations de la LATIN
  30. *----------------------------------------------------------------------
  31. 'DEBPROC' PLOTITER Tplus*CHPOINT Tmoins*CHPOINT Wplus*CHPOINT
  32. Wmoins*CHPOINT meshint*MAILLAGE;
  33. fac2 = 20;
  34. fac1 = 1E-12;
  35. VTplus = vect Tplus forc rouge;
  36. VTmoins = vect Tmoins forc vert;
  37. def1 = defo wplus meshint VTplus rouge;
  38. def2 = defo wmoins meshint VTmoins vert;
  39. def3 = def1 'ET' def2;
  40. 'TRACER' def3 'NCLK';
  41. 'FINPROC';
  42.  
  43.  
  44. *----------------------------------------------------------------------
  45. * LATIN : contient les étapes du solveur non linéaires LATIN
  46. *
  47. * ENTREE :
  48. * - crack : maillage de l'interface
  49. * - RHS0 : Chargement mécanique
  50. * - K1tot : matrice totale sans Kuwt (avec C.L.)
  51. * - Kuwt : matrice issu de RELA ACCRO FAIBLE
  52. * - Ktot1 : Kuwt et K1tot
  53. * - mod1tot : modele de tout le massif
  54. * - WPglo, WMglo, WPglo, WMglo : w+,w-,t+,t- solution du pas précédent
  55. * - f_crack : Coefficient de frottement entre les lèvres de la fissure
  56. * - ItMax : Nombre d'itération maximum
  57. * - eps1 : Précision demandée
  58. *
  59. * SORTIE :
  60. * - (U, W, T)
  61. * - erreur à chaque itération
  62. *----------------------------------------------------------------------
  63. 'DEBPROC' LATIN Tdisplay*TABLE crack*MAILLAGE K1tot*RIGIDITE
  64. mod1tot*MMODEL RHS0*CHPOINT IdPasT*ENTIER Ktot1*RIGIDITE Kuwt*RIGIDITE
  65. cht1*CHPOINT chn*CHPOINT ItMax*ENTIER eps1*FLOTTANT k*FLOTTANT
  66. f_crack*FLOTTANT mai1tot*MAILLAGE cht2/CHPOINT TPglo/CHPOINT
  67. TMglo/CHPOINT WPglo/CHPOINT WMglo/CHPOINT;
  68.  
  69.  
  70. * initialisation avec second menbre incomplet uniquement pour le premier pas de de temps
  71. * pour les suivants on initialise avec les champs locaux du pas précédents
  72.  
  73. *'SI' ((IdPasT 'EGA' 1) 'OU' ouv1);
  74. 'SI' (IdPasT 'EGA' 1);
  75.  
  76. *-----------------------------< Résolution initiale: étape globale
  77. ma1 = 'EXTRAIRE' mod1 'MAIL';
  78. cl2 = ('BLOQUE' 'AX'ma1) 'ET' ('BLOQUE' 'B1X'ma1) 'ET'
  79. ('BLOQUE' 'C1X'ma1) 'ET' ('BLOQUE' 'D1X'ma1) 'ET'
  80. ('BLOQUE' 'E1X'ma1) 'ET' ('BLOQUE' 'AY'ma1) 'ET'
  81. ('BLOQUE' 'B1Y'ma1) 'ET' ('BLOQUE' 'C1Y'ma1) 'ET'
  82. ('BLOQUE' 'D1Y'ma1) 'ET' ('BLOQUE' 'E1Y'ma1) 'ET'
  83. ('BLOQUE' 'AZ'ma1) 'ET' ('BLOQUE' 'B1Z'ma1) 'ET'
  84. ('BLOQUE' 'C1Z'ma1) 'ET' ('BLOQUE' 'D1Z'ma1) 'ET'
  85. ('BLOQUE' 'E1Z'ma1) ;
  86. ktot2 = ktot1 'ET' cl2;
  87. * u1 = RESO Ktot2 RHS0;
  88.  
  89. u1 = RESO Ktot1 RHS0;
  90.  
  91. *-----------------------------< extrait w et t fissure de u1
  92. * (dans leur formalisme global = moyenne + saut)
  93. wsauv fsauv = EXTCRACK u1 crack dim1;
  94.  
  95. *-----------------------------< recuperation champ + et -
  96. Wmoins Wplus Tmoins Tplus = SPLIT crack fsauv wsauv;
  97. 'SI' (Tdisplay.iteration);
  98. 'TITRE' 'apres etape GLOBALE INITIALE';
  99. PLOTITER Tplus Tmoins Wplus Wmoins crack;
  100. 'FINSI';
  101. *-----------------------------< Projection dans la base locale xxx_b = xxx_before (local step)
  102. 'SI' (('VALEUR' 'DIME') 'EGA' 3);
  103. TlocP_b TlocM_b WlocP_b WlocM_b = GLO2LOC cht1 chn Tplus Tmoins Wplus
  104. Wmoins cht2;
  105. 'SINON' ;
  106. TlocP_b TlocM_b WlocP_b WlocM_b = GLO2LOC cht1 chn Tplus Tmoins Wplus
  107. Wmoins ;
  108. 'FINSI';
  109.  
  110.  
  111. *-----------------------------< Etape locale pour la première itération
  112. TlocP TlocM WlocP WlocM = LOCSTEPI k f_crack crack TlocP_b TlocM_b
  113. WlocP_b WlocM_b;
  114.  
  115. *-----------------------------< Projection des w+- et t+- (i +1/2) dans le repère global
  116. 'SI' (('VALEUR' 'DIME') 'EGA' 3);
  117. WMglo WPglo TMglo TPglo = LOC2GLO cht1 chn TlocP TlocM WlocP WlocM
  118. crack cht2;
  119. 'SINON';
  120. WMglo WPglo TMglo TPglo = LOC2GLO cht1 chn TlocP TlocM WlocP WlocM
  121. crack;
  122. 'FINSI';
  123.  
  124. 'FINSI';
  125. *---> initialisation itérations
  126. it1 = 1 ;
  127.  
  128. * utlise solutions locales au pas de temps précédent pour initialiser
  129. 'SI' (IdPasT > 1);
  130.  
  131. * 'SI' (idpasT < 6); itmax = 10; 'FINSI';
  132.  
  133. *-----------------------------< Projection dans la base locale xxx_b = xxx_before (local step)
  134. 'SI' (('VALEUR' 'DIME') 'EGA' 3);
  135. TlocP_b TlocM_b WlocP_b WlocM_b = GLO2LOC cht1 chn TPglo TMglo
  136. WPglo WMglo cht2;
  137. 'SINON' ;
  138. TlocP_b TlocM_b WlocP_b WlocM_b = GLO2LOC cht1 chn TPglo TMglo
  139. WPglo WMglo ;
  140. 'FINSI';
  141.  
  142. TlocP = TlocP_b;
  143. TlocM = TlocM_b;
  144. WlocP = WlocP_b;
  145. WlocM = WlocM_b;
  146. 'FINSI';
  147. * sauvegarde champ locaux après étape locale
  148. ttlocp = table; ttlocm = table; twlocp = table; twlocm = table;
  149. ttlocp.it1 = tlocp_b; ttlocm.it1 = tlocm_b;
  150. twlocp.it1 = wlocp_b; twlocm.it1 = wlocm_b;
  151. 'SI' (Tdisplay.iteration);
  152. 'TITRE' 'apres etape LOCALE';
  153. PLOTITER TPglo TMglo WPglo WMglo crack;
  154. 'FINSI';
  155. *-----------------------------< Recombine les champs
  156. Wglo Tglo = LIPSRECO WMglo WPglo TMglo TPglo;
  157.  
  158. *-----------------------------< Mise à jour du second membre
  159. rhs1 = UPDATRHS Kuwt Tglo Wglo crack it1 rhs0;
  160.  
  161. *-----------------------------< Etape globale
  162. u1 = 'RESOUD' ktot1 rhs1;
  163.  
  164. *-----------------------------< extrait w et t fissure de u1
  165. * (dans leur formalisme global = moyenne + saut)
  166. wfiss fcon = EXTCRACK u1 crack dim1;
  167.  
  168. *-----------------------------< recuperation champ + et -
  169. Wmoins Wplus Tmoins Tplus = SPLIT Crack fcon wfiss;
  170. 'SI' (Tdisplay.iteration);
  171. 'TITRE' 'apres etape GLOBALE';
  172. PLOTITER Tplus Tmoins Wplus Wmoins crack;
  173. 'FINSI';
  174.  
  175.  
  176. *---< Projection dans la base locale
  177. 'SI' (('VALEUR' 'DIME') 'EGA' 3);
  178. TlocP_n TlocM_n WlocP_n WlocM_n = GLO2LOC cht1 chn Tplus Tmoins
  179. Wplus Wmoins cht2;
  180. 'SINON' ;
  181. TlocP_n TlocM_n WlocP_n WlocM_n = GLO2LOC cht1 chn Tplus Tmoins
  182. Wplus Wmoins ;
  183. 'FINSI';
  184.  
  185.  
  186. 'REPETER' blo7 (Itmax);
  187. it1 = it1 '+' 1;
  188.  
  189. TlocP_b = TlocP_n;
  190. TlocM_b = TlocM_n;
  191. WlocP_b = WlocP_n;
  192. WlocM_b = WlocM_n;
  193.  
  194. *-----------------------------< Etape locale (_a = after ; _b = before)
  195. TlocP_a TlocM_a WlocP_a WlocM_a ouv1 = LOCSTEP1
  196. k f_crack crack TlocP_b TlocM_b WlocP_b WlocM_b it1
  197. (twlocp.(it1-1)) (twlocm.(it1-1));
  198.  
  199. 'SI' (ouv1 'ET' (it1 > 15));
  200. 'QUITTER' blo7;
  201. 'FINSI';
  202.  
  203.  
  204. *-----------------------------< Projection des des w+- et t+- (i +1/2) dans le repère global
  205. 'SI' (('VALEUR' 'DIME') 'EGA' 3);
  206. WMglo WPglo TMglo TPglo = LOC2GLO cht1 chn TlocP_a TlocM_a WlocP_a
  207. WlocM_a crack cht2;
  208. 'SINON';
  209. WMglo WPglo TMglo TPglo = LOC2GLO cht1 chn TlocP_a TlocM_a WlocP_a
  210. WlocM_a crack;
  211. 'FINSI';
  212.  
  213. 'SI' (Tdisplay.iteration);
  214. 'TITRE' 'apres etape LOCALE';
  215. PLOTITER TPglo TMglo WPglo WMglo crack;
  216. 'FINSI';
  217. *-----------------------------< Recombine les champs
  218. Wglo Tglo = LIPSRECO WMglo WPglo TMglo TPglo;
  219.  
  220. *-----------------------------< Mise à jour du second membre
  221. rhs1 = UPDATRHS Kuwt Tglo Wglo crack it1 rhs0;
  222.  
  223. *-----------------------------< Etape globale
  224. u1 = 'RESOUD' ktot1 rhs1;
  225.  
  226. *-----------------------------< extrait w et t fissure de u1
  227. * (dans leur formalisme global = moyenne + saut)
  228. wfiss fcon = EXTCRACK u1 crack dim1;
  229.  
  230. *-----------------------------< recuperation champ + et -
  231. Wmoins Wplus Tmoins Tplus = SPLIT crack fcon wfiss;
  232. 'SI' (Tdisplay.iteration);
  233. 'TITRE' 'apres etape GLOABLE';
  234. PLOTITER Tplus Tmoins Wplus Wmoins crack;
  235. 'FINSI';
  236.  
  237. *-----------------------------< Projection dans la base locale
  238. 'SI' (('VALEUR' 'DIME') 'EGA' 3);
  239. TlocP_n TlocM_n WlocP_n WlocM_n = GLO2LOC cht1 chn Tplus Tmoins
  240. Wplus Wmoins cht2;
  241. 'SINON' ;
  242. TlocP_n TlocM_n WlocP_n WlocM_n = GLO2LOC cht1 chn Tplus Tmoins
  243. Wplus Wmoins ;
  244. 'FINSI';
  245.  
  246. * sauvegarge de l'histoire au cours des itérations
  247. * ttlocp.it1 = tlocp_n;
  248. * ttlocm.it1 = tlocm_n;
  249. * twlocp.it1 = wlocp_n;
  250. * twlocm.it1 = wlocm_n;
  251.  
  252. ttlocp.it1 = tlocp_a;
  253. ttlocm.it1 = tlocm_a;
  254. twlocp.it1 = wlocp_a;
  255. twlocm.it1 = wlocm_a;
  256.  
  257. 'SI' (it1 'EGA' 2);
  258. *-----------------------------< Calcul de l'indicateur d'erreur INITIAL
  259. TotalErr denN denT etaT etaN = ERRORINI TlocP_n TlocM_n WlocP_n
  260. WlocM_n TlocP TlocM WlocP WlocM k;
  261. lerror = 'PROG' TotalErr;letaT = 'PROG'etaT;letaN = 'PROG' etaN;
  262. 'SINON';
  263. *-----------------------------< Calcul de l'indicateur d'erreur (dénominateur fixé)
  264. TotalErr etaT etaN = GETERROR TlocP_n TlocM_n WlocP_n WlocM_n
  265. TlocP_a TlocM_a WlocP_a WlocM_a denN denT k;
  266. lerror = lerror'ET'TotalErr;letaT=letaT'ET'etaT;
  267. letaN=letaN'ET'etaN;
  268. 'FINSI';
  269. * 'MESSAGE' '';
  270. * 'MESSAGE' 'denN =' denN;
  271. * 'MESSAGE' 'denT =' denT;
  272. * 'MESSAGE' '';
  273.  
  274. SI (TotalErr < eps1) ;
  275. MESS 'La precision demandee a ete atteinte';
  276. QUITTER blo7 ;
  277. FINSI ;
  278.  
  279. SI (TotalErr > 1000000) ;
  280. 'MESSAGE' ' ';
  281. 'MESSAGE' ' ';
  282. 'MESSAGE' ' ';
  283. MESS '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!';
  284. MESS '!!! !!!';
  285. MESS '!!! LE CALCUL DIVERGE !! -> ARRET DU CALCUL !!!';
  286. MESS '!!! !!!';
  287. MESS '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!';
  288. 'MESSAGE' ' ';
  289. 'MESSAGE' ' ';
  290. 'MESSAGE' ' ';
  291. QUITTER blo7 ;
  292. FINSI ;
  293.  
  294. MESS ' ';
  295. MESS ' ';
  296. MESS '-----------> iteration :' it1 ;
  297. MESS 'Erreur :' totalerr;
  298.  
  299. 'FIN' blo7;
  300. * fin de la boucle sur les itération du solveur non-linéaire
  301.  
  302. 'MESSAGE' 'fin du cas fermeture';
  303. * 'FINSI';
  304. * fin du cas fermeture
  305.  
  306. * cas fissure complètement ouverte -> résoluion directe -> w = u sur la fissure et t = 0
  307. 'SI' ouv1;
  308. 'MESSAGE' ' ' ;
  309. 'MESSAGE' 'Fisssure COMPLETEMENT ouverte -> résolution DIRECTE';
  310. 'MESSAGE' ' ' ;
  311. u2 = 'RESOUD' k1tot RHS0;
  312. WMglo WPglo = XFEM 'FISS' crack u2 mod1tot;
  313. TPglo = 0.*Tplus; TMglo = 0.*Tplus;lerror = 'PROG' 0.;
  314. 'SINON';
  315. u2 = 'REDU' u1 mai1tot;
  316. 'FINSI';
  317.  
  318.  
  319. 'FINPROC' u2 TPglo TMglo WPglo WMglo lerror;
  320.  
  321.  
  322. *----------------------------------------------------------------------
  323. *
  324. * SAVFIELD : sauvegarde les champs dans les tables
  325. *
  326. *----------------------------------------------------------------------
  327.  
  328. 'DEBPROC' SAVFIELD u1*CHPOINT tplus*CHPOINT tmoins*CHPOINT wplus*CHPOINT
  329. wmoins*CHPOINT it1*ENTIER;
  330.  
  331. 'SI' (it1 'EGA' 0);
  332. ttplus = 'TABLE';
  333. ttmoins = 'TABLE';
  334. twplus = 'TABLE';
  335. twmoins = 'TABLE';
  336. tu = 'TABLE';
  337. 'FINSI';
  338. ttplus.it1 = tplus;
  339. ttmoins.it1 = tmoins;
  340. twplus.it1 = wplus;
  341. twmoins.it1 = wmoins;
  342. tu.it1 = u1;
  343. 'FINPROC' ttplus ttmoins twplus twmoins tu;
  344.  
  345.  
  346. *----------------------------------------------------------------------
  347. *
  348. * EXTCRACK : extrait les inconnus relatives à la fissure de u1
  349. *
  350. *----------------------------------------------------------------------
  351.  
  352. 'DEBPROC' EXTCRACK u1*CHPOINT meshint*MAILLAGE dim1*ENTIER;
  353.  
  354. 'SI' (dim1 'EGA' 2);
  355. lmot1 = 'MOTS' 'UX ' 'UY ' 'AX ' 'AY ';
  356. lmot2 = 'MOTS' 'FX ' 'FY ' 'FAX ' 'FAY ';
  357. 'SINON';
  358. lmot1 = 'MOTS' 'UX ' 'UY ' 'UZ ' 'AX ' 'AY ' 'AZ ' ;
  359. lmot2 = 'MOTS' 'FX ' 'FY ' 'FZ ' 'FAX ' 'FAY ' 'FAZ ' ;
  360. 'FINSI';
  361.  
  362. wsauv = 'ENLEVER' ('REDU' meshint u1) lmot2;
  363. fsauv = 'ENLEVER' ('REDU' meshint u1) lmot1;
  364.  
  365. 'FINPROC' wsauv fsauv;
  366.  
  367.  
  368. *----------------------------------------------------------------------
  369. *
  370. * SPLIT : Différencie les champs locaux sur les lèvres de la fissure
  371. * à partir des inconnus XFEM issus de l'étape globale
  372. *
  373. *----------------------------------------------------------------------
  374.  
  375. 'DEBPROC' SPLIT MeshInt*MAILLAGE fint1*CHPOINT uint1*CHPOINT;
  376.  
  377. 'SI' (('VALEUR' 'DIME') 'EGA' 2);
  378. lmot1 = 'MOTS' 'UX ' 'UY ';
  379. lmot2 = 'MOTS' 'AX ' 'AY ';
  380. lmot3 = 'MOTS' 'FX ' 'FY ';
  381. lmot4 = 'MOTS' 'FAX ' 'FAY ';
  382. 'SINON';
  383. lmot1 = 'MOTS' 'UX ' 'UY ' 'UZ ';
  384. lmot2 = 'MOTS' 'AX ' 'AY ' 'AZ ';
  385. lmot3 = 'MOTS' 'FX ' 'FY ' 'FZ ';
  386. lmot4 = 'MOTS' 'FAX ' 'FAY ' 'FAZ ';
  387. 'FINSI';
  388.  
  389. *-----------------------------< Déplacement
  390. umoy = 'EXCO' lmot1 uint1;
  391. uH = 'EXCO' lmot2 uint1;
  392. uH = 'NOMC' lmot2 lmot1 uH;
  393. usup = umoy '+' uH;
  394. uinf = umoy '-' uH;
  395. *-----------------------------< Effort
  396. tmoy = 'EXCO' lmot3 fint1;
  397. tH = 'EXCO' lmot4 fint1;
  398. tH = 'NOMC' lmot4 lmot3 tH;
  399. tsup = tmoy '+' tH;
  400. tinf = tmoy '-' tH;
  401.  
  402.  
  403. 'FINPROC' uinf usup tinf tsup;
  404.  
  405.  
  406. *----------------------------------------------------------------------
  407. *
  408. * LOCABA3D : Définit une lèvre plus et moins pour tous les éléments
  409. * d'interface à partir de la normale orienté dans le sens
  410. * de parcours des éléments d'interface
  411. *
  412. *----------------------------------------------------------------------
  413.  
  414. 'DEBPROC' LOCABA3D MeshInt*MAILLAGE mod_int*MMODEL;
  415.  
  416. chn = vsur mod_int 'NORM';
  417. mo1 = 'EXTRAIRE' chn 'COMP';
  418. mo2 = 'MOTS' 'NX ' 'NY ' 'NZ ';
  419. chn = 'CHANGER' 'CHPO' mod_int chn 'MOYE';
  420. chn = 'NOMC' chn mo1 mo2;
  421. cht1 = 0. * chn;
  422. * creation des tangentes arbitrairement orienté à partir de la normale
  423. xt1 = 'NOMC' 'T1X ' ('EXCO' 'NY ' chn);
  424. yt1 = 'NOMC' 'T1Y ' ('EXCO' 'NX ' chn);
  425. zt1 = 'NOMC' 'T1Z ' ('EXCO' 'NZ ' cht1);
  426. 'SI' ((('MAXIMUM' xt1) 'EGA' 0) 'ET' (('MAXIMUM' yt1) 'EGA' 0));
  427. xt1 = 'MANUEL' chpo ('EXTRAIRE' xt1 mail) 'T1X ' 1.;
  428. 'FINSI';
  429. cht1 = xt1 '+' yt1 '+' zt1;
  430. * on norme cht1
  431. m1 = 'EXTRAIRE' cht1 'COMP';
  432. cht1 = cht1 '/' (('PSCAL' cht1 cht1 m1 m1)**0.5);
  433.  
  434. * creation de la dernière tangente comme résultat du produit vectoriel
  435. xn1 = 'NOMC' 'NX ' ('EXCO' 'NX ' chn);
  436. yn1 = 'NOMC' 'NY ' ('EXCO' 'NY ' chn);
  437. zn1 = 'NOMC' 'NZ ' ('EXCO' 'NZ ' chn);
  438. xt1 = 'NOMC' 'T1X ' ('EXCO' 'T1X ' cht1);
  439. yt1 = 'NOMC' 'T1Y ' ('EXCO' 'T1Y ' cht1);
  440. zt1 = 'NOMC' 'T1Z ' ('EXCO' 'T1Z ' cht1);
  441. mxt1 = 'MOTS' 'T1X '; myt1 = 'MOTS' 'T1Y '; mzt1 = 'MOTS' 'T1Z ';
  442. mxt2 = 'MOTS' 'T2X '; myt2 = 'MOTS' 'T2Y '; mzt2 = 'MOTS' 'T2Z ';
  443. mxn1 = 'MOTS' 'NX '; myn1 = 'MOTS' 'NY '; mzn1 = 'MOTS' 'NZ ';
  444.  
  445. * pvec ne fonctionne pas avec les CHPO contrairement à ce qui est dit dans l'aide
  446. * on le fait à la main
  447. xt2 = (yn1 '*' zt1 myn1 mzt1 mxt2) '-' (zn1 '*' yt1 mzn1 myt1 mxt2);
  448. yt2 = (zn1 * xt1 mzn1 mxt1 myt2) '-' (xn1 * zt1 mxn1 mzt1 myt2);
  449. zt2 = (xn1 * yt1 mxn1 myt1 mzt2) '-' (yn1 * xt1 myn1 mxt1 mzt2);
  450. cht2 = xt2 '+' yt2 '+' zt2;
  451.  
  452.  
  453. 'FINPROC' chn cht1 cht2;
  454.  
  455.  
  456. *----------------------------------------------------------------------
  457. *
  458. * LOCABASE : Définit une lèvre plus et moins pour tous les éléments
  459. * d'interface à partir de la normale orienté dans le sens
  460. * de parcours des éléments d'interface
  461. * utiliser chn = vsur modele_interface 'NORM'; pour la création des normales
  462. *----------------------------------------------------------------------
  463.  
  464. 'DEBPROC' LOCABASE MeshInt*MAILLAGE;
  465.  
  466.  
  467. Nbe1em = 'NBEL' MeshInt;
  468.  
  469. * 'SI' (nbe1em '>EG' 6);
  470. 'SI' (nbe1em '>EG' 6E6);
  471. cht chn chb = @FRENET MeshInt;
  472.  
  473. nbn1 = 'NBNO' MeshInt;
  474. po1 = MeshInt 'POIN' 1;
  475. * Mise à zéro des normales pour le premier et dernier noeud
  476. chp1 = 'MANUEL' chpo po1 1 'SCAL' 0. natu discret;
  477. po2 = MeshInt 'POIN' nbn1;
  478. chp2 = 'MANUEL' chpo po2 1 'SCAL' 0. natu discret;
  479. chp3 = chp1 'ET' chp2;
  480.  
  481. idnode = 2;
  482. 'REPETER' blo1 (nbn1 '-' 2);
  483. po1 = MeshInt 'POIN' idnode;
  484. chp1 = 'MANUEL' chpo po1 1 'SCAL' 1. natu discret;
  485. chp3 = chp3 'ET' chp1;
  486. idnode = idnode '+' 1;
  487. 'FIN' blo1;
  488.  
  489. lmot1 = 'MOTS' 'SCAL';
  490. lmot2 = 'MOTS' 'UX ' 'UX ';
  491.  
  492. cht = chp3 * cht;
  493. chn = chp3 * chn;
  494.  
  495. * Cas du premier et dernier noeud de l'interface
  496.  
  497. * Premier noeud
  498. po1 = MeshInt 'POIN' 1;
  499. x1 = 'COORDONNEE' 1 po1;
  500. y1 = 'COORDONNEE' 2 po1;
  501.  
  502. po2 = MeshInt 'POIN' 2;
  503. x2 = 'COORDONNEE' 1 po2;
  504. y2 = 'COORDONNEE' 2 po2;
  505.  
  506. pmilieu = 'POIN' (x1 '+' ((x2'-' x1) '/' 2))
  507. (y1 '+' ((y2 '-' y1) '/' 2));
  508.  
  509. Lseg = (((x1 '-' x2) '**' 2) '+' ((y1 '-' y2) '**' 2 )) '**' (0.5);
  510. a = (x2 '-' x1) '/' Lseg;
  511. b = (y2 '-' y1) '/' Lseg;
  512. b2 = -1. * b;
  513.  
  514.  
  515. chp1 = 'MANUEL' chpo po1 2 'TX ' a 'TY ' b ;
  516. chp2 = 'MANUEL' chpo po1 2 'NX ' b2 'NY ' a;
  517.  
  518. cht = cht '+' chp1;
  519. chn = chn '+' chp2;
  520.  
  521.  
  522. * Dernier noeud
  523. po1 = MeshInt 'POIN' (nbn1 '-' 1);
  524. x1 = 'COORDONNEE' 1 po1;
  525. y1 = 'COORDONNEE' 2 po1;
  526.  
  527. po2 = MeshInt 'POIN' nbn1;
  528. x2 = 'COORDONNEE' 1 po2;
  529. y2 = 'COORDONNEE' 2 po2;
  530.  
  531. pmilieu = 'POIN' (x1 '+' ((x2'-' x1) '/' 2))
  532. (y1 '+' ((y2 '-' y1) '/' 2));
  533.  
  534. Lseg = (((x1 '-' x2) '**' 2) '+' ((y1 '-' y2) '**' 2 )) '**' (0.5);
  535. a = (x2 '-' x1) '/' Lseg;
  536. b = (y2 '-' y1) '/' Lseg;
  537. b2 = -1. * b;
  538.  
  539. chp1 = 'MANUEL' chpo po2 2 'TX ' a 'TY ' b ;
  540. chp2 = 'MANUEL' chpo po2 2 'NX ' b2 'NY ' a;
  541.  
  542. cht = cht '+' chp1;
  543. chn = chn '+' chp2;
  544.  
  545. co1 = 'EXTRAIRE' cht 'COMP';
  546. co2 = 'MOTS' 'T1X ' 'T1Y ';
  547.  
  548. cht = 'NOMC' cht co1 co2;
  549.  
  550. * cas avec moins de 5 éléments (@frenet ne fonctionne avec moins de 5 éléments)
  551. 'SINON';
  552.  
  553. cht = 'MANUEL' chpo MeshInt 'T1X ' 0. 'T1Y ' 0.;
  554. chn = 'MANUEL' chpo MeshInt 'NX ' 0. 'NY ' 0.;
  555. Nbe1em = 'NBEL' MeshInt;
  556. Nbnode = 'NBNO' MeshInt;
  557. idnode = 1;
  558.  
  559. 'REPETER' blo1 (Nbnode);
  560.  
  561. 'SI' ((idnode 'EGA' 1) 'OU' (idnode 'EGA' Nbnode));
  562. 'SI' (idnode 'EGA' 1);
  563. ele1 = 'ELEM' MeshInt 1;
  564. 'SINON';
  565. ele1 = 'ELEM' MeshInt nbe1em;
  566. 'FINSI';
  567. po1 = (ele1 point 1);
  568. x1 = 'COORDONNEE' 1 po1;
  569. y1 = 'COORDONNEE' 2 po1;
  570.  
  571. po2 = (ele1 point 2);
  572. x2 = 'COORDONNEE' 1 po2;
  573. y2 = 'COORDONNEE' 2 po2;
  574.  
  575. Lseg = (((x1 '-' x2) '**' 2) '+' ((y1 '-' y2) '**' 2 )) '**' (0.5);
  576. a = (x2 '-' x1) '/' Lseg;
  577. b = (y2 '-' y1) '/' Lseg;
  578. b2 = -1. * b;
  579.  
  580. 'SI' (idnode 'EGA' 1);
  581. chp1 = 'MANUEL' chpo po1 2 'T1X ' a 'T1Y ' b ;
  582. chp2 = 'MANUEL' chpo po1 2 'NX ' b2 'NY ' a;
  583. 'SINON';
  584. chp1 = 'MANUEL' chpo po2 2 'T1X ' a 'T1Y ' b ;
  585. chp2 = 'MANUEL' chpo po2 2 'NX ' b2 'NY ' a;
  586. 'FINSI';
  587.  
  588. 'SINON';
  589. po1 = (MeshInt point (idnode '-' 1));
  590. x1 = 'COORDONNEE' 1 po1;
  591. y1 = 'COORDONNEE' 2 po1;
  592. po2 = (MeshInt point idnode);
  593. x2 = 'COORDONNEE' 1 po2;
  594. y2 = 'COORDONNEE' 2 po2;
  595. po3 = (MeshInt point (idnode '+' 1));
  596. x3 = 'COORDONNEE' 1 po3;
  597. y3 = 'COORDONNEE' 2 po3;
  598.  
  599. Lseg1 = (((x1 '-' x2) '**' 2) '+' ((y1 '-' y2) '**' 2 )) '**' (0.5);
  600. Lseg2 = (((x3 '-' x2) '**' 2) '+' ((y3 '-' y2) '**' 2 )) '**' (0.5);
  601.  
  602. a1 = (x2 '-' x1) '/' Lseg1;
  603. b1 = (y2 '-' y1) '/' Lseg1;
  604. a2 = (x3 '-' x2) '/' Lseg2;
  605. b2 = (y3 '-' y2) '/' Lseg2;
  606. Lseg3 = (((a1 '+' a2)**2) '+' ((b1 '+' b2)**2)) '**' (0.5);
  607. chp1 = 'MANUEL' chpo po2 2 'T1X ' ((a1 + a2) '/' Lseg3) 'T1Y '
  608. ((b1 '+' b2) '/' Lseg3) ;
  609. chp2 = 'MANUEL' chpo po2 2 'NX '(-1 * (b2 '+' b1) '/' Lseg3)
  610. 'NY ' ((a1 + a2) '/' Lseg3);
  611. 'FINSI' ;
  612.  
  613. cht = cht '+' chp1;
  614. chn = chn '+' chp2;
  615. idnode = idnode '+' 1;
  616.  
  617. 'FIN' blo1;
  618.  
  619. 'FINSI';
  620.  
  621. 'FINPROC' chn cht;
  622.  
  623.  
  624. *----------------------------------------------------------------------
  625. *
  626. * GLO2LOC : Projection des champs de l'interface du repère global
  627. * vers le repère local
  628. *
  629. *----------------------------------------------------------------------
  630.  
  631. 'DEBPROC' GLO2LOC cht1*CHPOINT chn*CHPOINT Tplus*CHPOINT Tmoins*CHPOINT
  632. Wplus*CHPOINT Wmoins*CHPOINT cht2/CHPOINT;
  633.  
  634. fg1 = 'EXTRAIRE' tplus 'COMP';
  635. wg1 = 'EXTRAIRE' wplus 'COMP';
  636. t1 = 'EXTRAIRE' cht1 'COMP';
  637. n1 = 'EXTRAIRE' chn 'COMP' ;
  638. 'SI' (dim1 'EGA' 3);
  639. t2 = 'EXTRAIRE' cht2 'COMP';
  640. 'FINSI' ;
  641.  
  642. a = 'NOMC' 'FT1 ' (psca cht1 Tplus t1 fg1);
  643. b = 'NOMC' 'FN ' (psca chn Tplus n1 fg1);
  644. TlocP = a '+' b;
  645.  
  646. a = 'NOMC' 'FT1 ' (psca cht1 Tmoins t1 fg1);
  647. b = 'NOMC' 'FN ' (psca chn Tmoins n1 fg1);
  648. TlocM = a '+' b;
  649.  
  650. a = 'NOMC' 'UT1 ' (psca cht1 Wplus t1 wg1);
  651. b = 'NOMC' 'UN ' (psca chn Wplus n1 wg1);
  652. WlocP = a '+' b;
  653.  
  654. a = 'NOMC' 'UT1 ' (psca cht1 Wmoins t1 wg1);
  655. b = 'NOMC' 'UN ' (psca chn Wmoins n1 wg1);
  656. WlocM = a '+' b;
  657.  
  658. 'SI' (('DIME' fg1) 'EGA' 3);
  659. c = 'NOMC' 'FT2 ' (psca cht2 Tplus t2 fg1);
  660. TlocP = TlocP '+' c;
  661. c = 'NOMC' 'FT2 ' (psca cht2 Tmoins t2 fg1);
  662. TlocM = TlocM '+' c;
  663. c = 'NOMC' 'UT2 ' (psca cht2 Wplus t2 wg1);
  664. WlocP = WlocP '+' c;
  665. c = 'NOMC' 'UT2 ' (psca cht2 Wmoins t2 wg1);
  666. WlocM = WlocM '+' c;
  667. 'FINSI';
  668.  
  669. 'FINPROC' TlocP TlocM WlocP WlocM;
  670.  
  671.  
  672. *----------------------------------------------------------------------
  673. *
  674. * LOCSTEP1 : Etape locale pour le premier par de temps
  675. *
  676. *----------------------------------------------------------------------
  677.  
  678. 'DEBPROC' LOCSTEP1 Klatin*FLOTTANT f_crack*FLOTTANT MeshInt*MAILLAGE
  679. TlocP*CHPOINT TlocM*CHPOINT WlocP*CHPOINT WlocM*CHPOINT it1*ENTIER
  680. WlocPOld*CHPOINT WlocMOld*CHPOINT;
  681.  
  682. meshloc = 'CHANGER' POI1 meshint;
  683. wdiffold= wlocmold '-' wlocpold;
  684.  
  685. wdiff= wlocm '-' wlocp;
  686. tdiff= tlocm '-' tlocp;
  687. k_tan = klatin;
  688. k_nor = klatin;
  689.  
  690. TlocPnor = 'EXCO' 'FN ' TlocP;
  691. TlocPta1 = 'EXCO' 'FT1 ' TlocP;
  692. TlocMnor = 'EXCO' 'FN ' TlocM;
  693. TlocMta1 = 'EXCO' 'FT1 ' TlocM;
  694.  
  695. WlocPnor = 'EXCO' 'UN ' WlocP;
  696. WlocPta1 = 'EXCO' 'UT1 ' WlocP;
  697. WlocMnor = 'EXCO' 'UN ' WlocM;
  698. WlocMta1 = 'EXCO' 'UT1 ' WlocM;
  699.  
  700. w1old = 'EXCO' 'UT1 ' wdiffold;
  701. wm1old = 'EXCO' 'UT1 ' wlocmold;
  702. wp1old = 'EXCO' 'UT1 ' wlocpold;
  703.  
  704. * initialisation booleen pour détecter fissure complètement ouverte
  705. ouv1 = vrai;
  706.  
  707. 'SI' (('DIME' ('EXTRAIRE' tlocp 'COMP')) 'EGA' 3);
  708. TlocPta2 = 'EXCO' 'FT2 ' TlocP; TlocMta2 = 'EXCO' 'FT2 ' TlocM;
  709. WlocPta2 = 'EXCO' 'UT2 ' WlocP; WlocMta2 = 'EXCO' 'UT2 ' WlocM;
  710. TMta2 = 0. * TlocMta2;TPta2 = 0. * TlocPta2;
  711. WMta2 = 0. * WlocMta2;WPta2 = 0. * WlocPta2;
  712. w2old = 'EXCO' 'UT2 ' wdiffold;
  713. wm2old = 'EXCO' 'UT2 ' wlocmold;
  714. wp2old = 'EXCO' 'UT2 ' wlocpold;
  715. B3D = vrai;
  716. 'SINON';
  717. B3D = faux;
  718. 'FINSI';
  719. *
  720. * Cacul des indicateur de contact et de frottement initial
  721. ConIni = 0.5 * (('EXCO' 'UN 'wdiff)'-'(('EXCO' 'FN' tdiff)'/'k_nor));
  722.  
  723. * btrolle le 5 novembre pour coller avec ceux qui est écrit dans NewlocalStepOnInterface dans ELFE
  724. GliIni1 =(0.5*((K_tan*('EXCO' 'UT1 'wdiff))'-'('EXCO' 'FT1' tdiff)))'-'
  725. (0.5 * K_tan * w1old);
  726.  
  727. 'SI' B3D;
  728. GliIni2 =(0.5*((K_tan*('EXCO' 'UT2 'wdiff))'-'('EXCO' 'FT2' tdiff)))'-'
  729. (0.5 * K_tan * w2old);
  730. GliIni = ((gliIni1**2)'+'(gliIni2**2))'**'(0.5);
  731. 'SINON';
  732. GliIni = GliIni1;
  733. 'FINSI' ;
  734.  
  735. * seuil pour le glissement
  736. GSeuil = f_crack * ('ABS' (K_nor * ConIni));
  737.  
  738. * Initialisation des champs correspondant à l'itération i + 1/2
  739. TMnor = 0. * TlocMnor; TMta1 = 0. * TlocMta1;
  740. TPnor = 0. * TlocMnor; TPta1 = 0. * TlocPta1;
  741. WMnor = 0. * WlocMnor; WMta1 = 0. * WlocMta1;
  742. WPnor = 0. * WlocMnor; WPta1 = 0. * WlocPta1;
  743.  
  744. * initialisation des compteurs
  745. Ncon = 0;
  746. Nouv = 0;
  747.  
  748. *
  749. * Loi de comportement local i --> i + 1/2 avec Contact with friction (coulomb)
  750. *---
  751. Nbn1 = 'NBNO' Meshloc;
  752. idnode = 1;
  753. * boucle sur les noeuds de l'interface
  754. 'REPETER' blo1 (Nbn1);
  755.  
  756. ptemp = Meshloc 'POIN' idnode;
  757. ValCI = 'EXTRAIRE' ConIni 'SCAL' ptemp;
  758.  
  759. * Valeur à l'itération i
  760. wpn = 'EXTRAIRE' WlocPnor 'SCAL' ptemp;
  761. wmn = 'EXTRAIRE' WlocMnor 'SCAL' ptemp;
  762. wpt1 = 'EXTRAIRE' WlocPta1 'SCAL' ptemp;
  763. wmt1 = 'EXTRAIRE' WlocMta1 'SCAL' ptemp;
  764.  
  765. tpn = 'EXTRAIRE' TlocPnor 'SCAL' ptemp;
  766. tmn = 'EXTRAIRE' TlocMnor 'SCAL' ptemp;
  767. tpt1 = 'EXTRAIRE' TlocPta1 'SCAL' ptemp;
  768. tmt1 = 'EXTRAIRE' TlocMta1 'SCAL' ptemp;
  769.  
  770. wpo1 = 'EXTRAIRE' wp1old 'SCAL' ptemp;
  771. wmo1 = 'EXTRAIRE' wm1old 'SCAL' ptemp;
  772.  
  773. 'SI' B3D;
  774. tpt2 = 'EXTRAIRE' TlocPta2 'SCAL' ptemp;
  775. tmt2 = 'EXTRAIRE' TlocMta2 'SCAL' ptemp;
  776. wpt2 = 'EXTRAIRE' WlocPta2 'SCAL' ptemp;
  777. wmt2 = 'EXTRAIRE' WlocMta2 'SCAL' ptemp;
  778. wpo2 = 'EXTRAIRE' wp2old 'SCAL' ptemp;
  779. wmo2 = 'EXTRAIRE' wm2old 'SCAL' ptemp;
  780. 'FINSI';
  781.  
  782. * Si ouverture : >EG 0. pour le cas ou on bloque les sauts passe par là ;-)
  783. 'SI' (ValCI '>EG' 0.);
  784. Nouv = Nouv '+' 1;
  785. * 'MESSAGE' '---> Ouvert';
  786. * calcul des nouveaux déplacemtents (t= 0)
  787. chp1 = 'MANUEL' chpo ptemp 'SCAL' (wpn '-' (tpn '/' K_nor));
  788. chp2 = 'MANUEL' chpo ptemp 'SCAL' (wmn '-' (tmn '/' K_nor));
  789. chp3 = 'MANUEL' chpo ptemp 'SCAL' (wpt1 '-' (tpt1 '/' K_nor));
  790. chp4 = 'MANUEL' chpo ptemp 'SCAL' (wmt1 '-' (tmt1 '/' K_nor));
  791. WPnor = WPnor '+' chp1; WMnor = WMnor '+' chp2;
  792. WPta1 = WPta1 '+' chp3; WMta1 = WMta1 '+' chp4;
  793.  
  794. 'SI' B3D;
  795. chp5 = 'MANUEL' chpo ptemp 'SCAL' (wpt2 '-' (tpt2 '/' K_nor));
  796. chp6 = 'MANUEL' chpo ptemp 'SCAL' (wmt2 '-' (tmt2 '/' K_nor));
  797. WPta2 = WPta2 '+' chp5; WMta2 = WMta2 '+' chp6;
  798. 'FINSI';
  799.  
  800. * Si contact
  801. 'SINON';
  802. Ncon = Ncon '+' 1;
  803.  
  804. * Pas de test d'ouverture pour le premier et le dernier noeud.. il faudrait le faire juste pour le ou les fronts pour faire les choses correctements...
  805. 'SI' ((idnode > 1) 'ET' (idnode < Nbn1)); ouv1 = faux; 'FINSI';
  806. * problème normal
  807. * w
  808. chp1 = 'MANUEL' chpo ptemp 'SCAL' (0.5 * ((wpn '+' wmn ) '-'
  809. ((tpn '+' tmn) '/' K_nor)));
  810. WPnor = WPnor '+' chp1; WMnor = WMnor '+' chp1;
  811.  
  812. * t
  813. chp2 = 'MANUEL' chpo ptemp 'SCAL' (K_nor * ValCI);
  814. TPnor = TPnor '+' chp2; TMnor = TMnor '-' chp2;
  815.  
  816. * extrait les valeurs pour tester adhérence ou glissement
  817. ValGI = 'EXTRAIRE' GliIni 'SCAL' ptemp;
  818. G = 'EXTRAIRE' GSeuil 'SCAL' ptemp;
  819. ValGI1 = 'EXTRAIRE' GliIni1 'SCAL' ptemp;
  820. si B3D;
  821. ValGI2 = 'EXTRAIRE' GliIni2 'SCAL' ptemp;
  822. 'FINSI';
  823.  
  824. * Adhérence
  825. 'SI' (('ABS' ValGI) '<' G);
  826. * 'MESSAGE' '---> ADHERENCE';
  827. * 'MESSAGE' 'ValGI1 = ' ValGI1 'G =' G 'ValCI =' ValCI;
  828. * t
  829. chp31 = 'MANUEL' chpo ptemp 'SCAL' valGI1;
  830. TPta1 = TPta1 '+' chp31; TMta1 = TMta1 '-' chp31;
  831.  
  832. * w
  833. chp41='MANU' chpo ptemp'SCAL'((wpt1 '-' wpo1) '+'((valGI1'-'
  834. tpt1) '/' K_tan));
  835. WPta1 = WPta1 '+' chp41'+' ('MANU' chpo ptemp'SCAL' wpo1);
  836. WMta1 = WMta1 '+' chp41'+' ('MANU' chpo ptemp'SCAL' wmo1);
  837.  
  838. 'SI' B3D;
  839. * t
  840. chp32 = 'MANUEL' chpo ptemp 'SCAL' valGI2;
  841. TPta2 = TPta2 '+' chp32; TMta2 = TMta2 '-' chp32;
  842.  
  843. * w
  844. chp42='MANU' chpo ptemp 'SCAL' ((wpt2 '-' wpo2)'+'((valGI2'-'
  845. tpt2) '/' K_tan));
  846. WPta2 = WPta2 '+' chp42'+' ('MANU' chpo ptemp'SCAL' wpo2);
  847. WMta2 = WMta2 '+' chp42'+' ('MANU' chpo ptemp'SCAL' wmo2);
  848. 'FINSI';
  849.  
  850. * Glissement
  851. 'SINON';
  852. * 'LISTE' ptemp;
  853. * 'MESSAGE' '---> Glissement';
  854. * 'MESSAGE' 'ValGI = ' ValGI 'G =' G 'ValCI =' ValCI;
  855.  
  856. * t
  857. 'SI' (valGI 'EGA' 0.);
  858. * chp31 = 'MANUEL' chpo ptemp 'SCAL' 0.;
  859. 'SINON' ;
  860. chp31 = 'MANUEL' chpo ptemp 'SCAL'(G*valGI1'/'('ABS' valGI));
  861. * chp31 = 'MANUEL' chpo ptemp 'SCAL'(G*valGI1'/'(2*('ABS' valGI)));
  862. 'FINSI';
  863. TPta1 = TPta1 '+' chp31; TMta1 = TMta1 '-' chp31;
  864.  
  865. * w
  866. * newlocalstep oninterface
  867. chp41 = 'MANUEL' chpo ptemp 'SCAL'
  868. (wpt1 '-' wpo1 '+'( ( (G*valGI1'/'('ABS' valGI)) '-' tpt1) '/' K_tan));
  869. chp51 = 'MANUEL' chpo ptemp 'SCAL'
  870. (wmt1 '-' wmo1 '+'(((-1.*G*valGI1'/'('ABS' valGI))'-'tmt1) '/' K_tan));
  871. WPta1 = WPta1 '+' chp41 '+' ('MANU' chpo ptemp'SCAL' wpo1);
  872. WMta1 = WMta1 '+' chp51 '+' ('MANU' chpo ptemp'SCAL' wmo1);
  873.  
  874. 'SI' B3D;
  875. 'SI' (valGI 'EGA' 0.);
  876. chp32 = 'MANUEL' chpo ptemp 'SCAL' 0.;
  877. 'SINON' ;
  878. chp32 = 'MANUEL' chpo ptemp 'SCAL'(G*valGI2'/'('ABS' valGI));
  879. 'FINSI';
  880. chp42 = 'MANUEL' chpo ptemp 'SCAL'
  881. ((wpt2 '- 'wpo2) '+' ( ( (G*valGI2'/'('ABS' valGI))'-'tpt2)'/'K_tan));
  882. chp52 = 'MANUEL' chpo ptemp 'SCAL'
  883. ((wmt2 '-' wmo2)'+'(((-1.*G*valGI2'/'('ABS' valGI))'-'tmt2)'/'K_tan));
  884. TPta2 = TPta2 '+' chp32; TMta2 = TMta2 '-' chp32;
  885. WPta2 = WPta2 '+' chp42'+' ('MANU' chpo ptemp'SCAL' wpo2);
  886. WMta2 = WMta2 '+' chp52'+' ('MANU' chpo ptemp'SCAL' wmo2);
  887. 'FINSI';
  888.  
  889. 'FINSI';
  890.  
  891. 'FINSI';
  892.  
  893. idnode = idnode '+' 1;
  894.  
  895. 'FIN' blo1;
  896.  
  897. tpnor = 'NOMC' 'FN ' tpnor natu discret;
  898. tpta1 = 'NOMC' 'FT1 ' tpta1 natu discret;
  899. tmnor = 'NOMC' 'FN ' tmnor natu discret;
  900. tmta1 = 'NOMC' 'FT1 ' tmta1 natu discret;
  901. wpnor = 'NOMC' 'UN ' wpnor natu diffus;
  902. wpta1 = 'NOMC' 'UT1 ' wpta1 natu diffus;
  903. wmnor = 'NOMC' 'UN ' wmnor natu diffus;
  904. wmta1 = 'NOMC' 'UT1 ' wmta1 natu diffus;
  905.  
  906. TlocP_a = TPta1 'ET' TPnor;
  907. TlocM_a = TMta1 'ET' TMnor;
  908. WlocP_a = WPta1 'ET' WPnor;
  909. WlocM_a = WMta1 'ET' WMnor;
  910.  
  911. 'SI' B3D;
  912. tpta2 = 'NOMC' 'FT2 ' tpta2 natu discret;
  913. tmta2 = 'NOMC' 'FT2 ' tmta2 natu discret;
  914. wpta2 = 'NOMC' 'UT2 ' wpta2 natu diffus;
  915. wmta2 = 'NOMC' 'UT2 ' wmta2 natu diffus;
  916. TlocP_a = TlocP_a 'ET' Tpta2;
  917. TlocM_a = TlocM_a 'ET' Tmta2;
  918. WlocP_a = WlocP_a 'ET' WPta2;
  919. WlocM_a = WlocM_a 'ET' WPta2;
  920. 'FINSI' ;
  921.  
  922.  
  923. 'MESSAGE' ' ';
  924. 'MESSAGE' 'Nombre de noeud en contact = ' Ncon;
  925. 'MESSAGE' 'Nombre de noeud ouvert = ' Nouv;
  926. 'MESSAGE' ' ';
  927.  
  928. 'FINPROC' TlocP_a TlocM_a WlocP_a WlocM_a ouv1;
  929.  
  930.  
  931. *----------------------------------------------------------------------
  932. *
  933. * LOC2GLO : Projetion des champs à l'interface du repère local vers
  934. * le repère global
  935. *
  936. *----------------------------------------------------------------------
  937.  
  938. 'DEBPROC' LOC2GLO cht1*CHPOINT chn*CHPOINT TlocP*CHPOINT TlocM*CHPOINT
  939. WlocP*CHPOINT WlocM*CHPOINT MeshInt*MAILLAGE cht2/CHPOINT;
  940.  
  941. fl1 = 'EXTRAIRE' TlocM 'COMP' ;
  942. wl1 = 'EXTRAIRE' WlocM 'COMP' ;
  943. t1 = 'EXTRAIRE' cht1 'COMP' ;
  944. n1 = 'EXTRAIRE' chn 'COMP' ;
  945. vecx = 'MANUEL' chpo MeshInt 'UX ' 1. 'UY ' 0.;
  946. vecy = 'MANUEL' chpo MeshInt 'UX ' 0. 'UY ' 1.;
  947.  
  948. 'SI' (('DIME' fl1) 'EGA' 3);
  949. vecx = 'MANUEL' chpo MeshInt 'UX ' 1. 'UY ' 0. 'UZ ' 0. ;
  950. vecy = 'MANUEL' chpo MeshInt 'UX ' 0. 'UY ' 1. 'UZ ' 0. ;
  951. vecz = 'MANUEL' chpo MeshInt 'UX ' 0. 'UY ' 0. 'UZ ' 1. ;
  952. t2 = 'EXTRAIRE' cht2 'COMP';
  953. B3D = vrai;
  954. 'SINON';
  955. B3D = faux;
  956. 'FINSI';
  957.  
  958. mx1 = 'EXTRAIRE' vecx 'COMP';
  959.  
  960. a = 'NOMC' 'UT1 ' ('PSCA' vecx cht1 mx1 t1);
  961. b = 'NOMC' 'UN ' ('PSCA' vecx chn mx1 n1);
  962. chx = a '+' b;
  963. 'SI' B3D;
  964. c = 'NOMC' 'UT2 '('PSCA' vecx cht2 mx1 t2);
  965. chx = chx '+' c;
  966. 'FINSI';
  967.  
  968. a = 'NOMC' 'UT1 ' ('PSCA' vecy cht1 mx1 t1);
  969. b = 'NOMC' 'UN ' ('PSCA' vecy chn mx1 n1);
  970. chy = a '+' b;
  971. 'SI' B3D;
  972. c = 'NOMC' 'UT2 '('PSCA' vecy cht2 mx1 t2);
  973. chy = chy '+' c;
  974. a = 'NOMC' 'UT1 ' ('PSCA' vecz cht1 mx1 t1);
  975. b = 'NOMC' 'UN ' ('PSCA' vecz chn mx1 n1);
  976. c = 'NOMC' 'UT2 '('PSCA' vecz cht2 mx1 t2);
  977. chz = a '+' b '+' c;
  978. 'FINSI';
  979.  
  980. a = 'NOMC' 'FY ' (psca chy TlocP wl1 fl1);
  981. b = 'NOMC' 'FX ' (psca chx TlocP wl1 fl1);
  982. TPglo = b '+' a;
  983.  
  984. a = 'NOMC' 'FY ' (psca chy TlocM wl1 fl1);
  985. b = 'NOMC' 'FX ' (psca chx TlocM wl1 fl1);
  986. TMglo = b '+' a;
  987.  
  988. a = 'NOMC' 'UY ' (psca chy WlocP wl1 wl1);
  989. b = 'NOMC' 'UX ' (psca chx WlocP wl1 wl1);
  990. WPglo = b '+' a;
  991.  
  992. a = 'NOMC' 'UY ' (psca chy WlocM wl1 wl1);
  993. b = 'NOMC' 'UX ' (psca chx WlocM wl1 wl1);
  994. WMglo = b '+'a;
  995.  
  996. 'SI' B3D;
  997. c = 'NOMC' 'FZ ' (psca chz TlocP wl1 fl1);
  998. TPglo = TPglo '+' c;
  999. c = 'NOMC' 'FZ ' (psca chz TlocM wl1 fl1);
  1000. TMglo = TMglo '+' c;
  1001. c = 'NOMC' 'UZ ' (psca chz WlocP wl1 wl1);
  1002. WPglo = WPglo '+' c;
  1003. c = 'NOMC' 'UZ ' (psca chz WlocM wl1 wl1);
  1004. WMglo = WMglo '+' c;
  1005. 'FINSI';
  1006.  
  1007. 'FINPROC' WMglo WPglo TMglo TPglo;
  1008.  
  1009.  
  1010. *----------------------------------------------------------------------
  1011. *
  1012. * LIPSRECO : Recombine les champs locaux des 2 lèvres de la fissure
  1013. * issus de l'étape locale
  1014. * sur les maillages de l'étape globale en différenciant
  1015. * et en calculant les inconnus XFEM correspondantes
  1016. *
  1017. *----------------------------------------------------------------------
  1018.  
  1019. 'DEBPROC' LIPSRECO Wmoins*CHPOINT Wplus*CHPOINT Tmoins*CHPOINT
  1020. Tplus*CHPOINT;
  1021.  
  1022. Tmoy = (Tmoins '+' Tplus) '/' 2;
  1023. Wmoy = (Wmoins '+' Wplus) '/' 2;
  1024. tmo1 = 'EXTRAIRE' tmoy 'COMP';
  1025. wmo1 = 'EXTRAIRE' wmoy 'COMP';
  1026. 'SI' (('DIME' wmo1) 'EGA' 2);
  1027. tmo2 = 'MOTS' 'FAX ' 'FAY ';
  1028. wmo2 = 'MOTS' 'AX ' 'AY ';
  1029. 'SINON';
  1030. tmo2 = 'MOTS' 'FAX ' 'FAY ' 'FAZ ';
  1031. wmo2 = 'MOTS' 'AX ' 'AY ' 'AZ ';
  1032. 'FINSI';
  1033.  
  1034. Tsaut = 'NOMC' tmo1 tmo2 ((Tplus '-' Tmoy));
  1035. Wsaut = 'NOMC' wmo1 wmo2 ((Wplus '-' Wmoy));
  1036. Tglo = Tmoy '+' Tsaut;
  1037. Wglo = Wmoy '+' Wsaut;
  1038.  
  1039. 'FINPROC' Wglo Tglo;
  1040.  
  1041.  
  1042. *----------------------------------------------------------------------
  1043. *
  1044. * UPDATRHS : Mise à jour du second membre du sytème linéaire
  1045. *
  1046. *----------------------------------------------------------------------
  1047.  
  1048. 'DEBPROC' UPDATRHS Kuwt*RIGIDITE Tglo*CHPOINT Wglo*CHPOINT
  1049. meshint*MAILLAGE it1*ENTIER f1*CHPOINT;
  1050.  
  1051. lmot1 = 'EXTRAIRE' wglo 'COMP';
  1052. lmot2 = 'EXTRAIRE' Tglo 'COMP';
  1053. f2 = 'ENLEVER' ('REDU' (kuwt * wglo) MeshInt) lmot1;
  1054. f3 = 'ENLEVER' ('REDU' (kuwt * tglo) MeshInt) lmot1;
  1055. f4 = enle ('REDU' (kuwt * tglo) MeshInt) lmot2;
  1056. f5 = f3 '+' f2 '+' f4;
  1057. RHS1 = f1 '+' f5;
  1058.  
  1059. 'FINPROC' rhs1;
  1060.  
  1061.  
  1062. *----------------------------------------------------------------------
  1063. *
  1064. * ERRORINI : Calcul de l'indicateur d'erreur INITIALE
  1065. * basé sur les quantités w et t sur l'interface
  1066. * dans les bases locales
  1067. * entre la solition i + 1/2 et i + 1
  1068. * fixe les dénominateurs pour les appels suivant -> GETERROR
  1069. *----------------------------------------------------------------------
  1070.  
  1071. 'DEBPROC' ERRORINI TlocP_n*CHPOINT TlocM_n*CHPOINT WlocP_n*CHPOINT
  1072. WlocM_n*CHPOINT TlocP*CHPOINT TlocM*CHPOINT WlocP*CHPOINT WlocM*CHPOINT
  1073. Klatin*FLOTTANT;
  1074.  
  1075.  
  1076. *-> i + 1/2 : TlocP TlocM WlocP WlocM
  1077. TPnor = 'EXCO' 'FN ' TlocP;
  1078. TPta1 = 'EXCO' 'FT1 ' TlocP;
  1079. TMnor = 'EXCO' 'FN ' TlocM;
  1080. TMta1 = 'EXCO' 'FT1 ' TlocM;
  1081.  
  1082. WPnor = 'EXCO' 'UN ' WlocP;
  1083. WPta1 = 'EXCO' 'UT1 ' WlocP;
  1084. WMnor = 'EXCO' 'UN ' WlocM;
  1085. WMta1 = 'EXCO' 'UT1 ' WlocM;
  1086.  
  1087. *-> i + 1 : TlocP_n TlocM_n WlocP_n WlocM_n ( _n = new )
  1088. TPnor_n = 'EXCO' 'FN ' TlocP_n;
  1089. TPta1_n = 'EXCO' 'FT1 ' TlocP_n;
  1090. TMnor_n = 'EXCO' 'FN ' TlocM_n;
  1091. TMta1_n = 'EXCO' 'FT1 ' TlocM_n;
  1092.  
  1093. WPnor_n = 'EXCO' 'UN ' WlocP_n;
  1094. WPta1_n = 'EXCO' 'UT1 ' WlocP_n;
  1095. WMnor_n = 'EXCO' 'UN ' WlocM_n;
  1096. WMta1_n = 'EXCO' 'UT1 ' WlocM_n;
  1097.  
  1098.  
  1099. * ---> NUMERATEUR
  1100.  
  1101. * problème normal
  1102. numN = 'MAXIMUM' (
  1103. ((((TPnor_n '-' TPnor)**2)'+' ((TMnor_n '-' TMnor)**2)) '/' Klatin) '+'
  1104. ((((WPnor_n '-' WPnor)**2)'+' ((WMnor_n '-' WMnor)**2)) '*' Klatin));
  1105.  
  1106. * problème tangentiel
  1107. numT= 'MAXIMUM' (
  1108. ((((TPta1_n '-' TPta1)**2)'+' ((TMta1_n '-' TMta1)**2)) '/' Klatin) '+'
  1109. ((((WPta1_n '-' WPta1)**2)'+' ((WMta1_n '-' WMta1)**2)) '*' Klatin));
  1110.  
  1111.  
  1112. * ---> DENOMINATEUR
  1113.  
  1114. * problème normal
  1115. SNor_n = (((TPnor_n * TPnor_n) '+' (TMnor_n * TMnor_n) ) '/' Klatin) '+'
  1116. (Klatin * ((WPnor_n * WPnor_n) '+' (WMnor_n * WMnor_n)));
  1117.  
  1118. SNor = (((TPnor * TPnor) '+' (TMnor * TMnor) ) '/' Klatin) '+'
  1119. (Klatin * ((WPnor * WPnor) '+' (WMnor * WMnor)));
  1120.  
  1121. denN = ('MAXIMUM' SNor_n ) '+' ('MAXIMUM' SNor );
  1122.  
  1123.  
  1124. * problème tangentiel
  1125. STa1_n = (((TPta1_n * TPta1_n) '+' (TMta1_n * TMta1_n) ) '/' Klatin) '+'
  1126. (Klatin * ((WPta1_n * WPta1_n) '+' (WMta1_n * WMta1_n)));
  1127.  
  1128. STa1 = (((TPta1 * TPta1) '+' (TMta1 * TMta1) ) '/' Klatin) '+'
  1129. (Klatin * ((WPta1 * WPta1) '+' (WMta1 * WMta1)));
  1130.  
  1131. denT = ('MAXIMUM' STa1_n ) '+' ('MAXIMUM' STa1 );
  1132.  
  1133. * ---> ERREUR
  1134.  
  1135. *-> erreur normale
  1136. etaN = numN '/' denN;
  1137.  
  1138.  
  1139. *-> erreur tangentielle
  1140. etaT = numT '/' denT;
  1141.  
  1142. *-> erreur totale : ComputeLatinErrorNew dans ELFE-3D
  1143. totalErr = ((etaN) '+' (etaT )) '**' 0.5;
  1144.  
  1145. *-> erreur totale : publi récente (= manuscrit Epierres)
  1146. *ler1 = 'PROG' (etaT**0.5) (etaN**0.5);
  1147. *totalerr = 'MAXIMUM' ler1;
  1148. * 'SI' (etaT > etaN);
  1149. * totalErr = etaT;
  1150. * 'SINON';
  1151. * totalErr = etaN
  1152. * 'FINSI';
  1153.  
  1154. 'FINPROC' totalErr (denN) (denT) (etaT**0.5) (etaN**0.5);
  1155.  
  1156.  
  1157. *----------------------------------------------------------------------
  1158. *
  1159. * GETERROR : Calcul de l'indicateur d'erreur
  1160. * basé sur les quantités w et t sur l'interface
  1161. * dans les bases locales
  1162. * entre la solition i + 1/2 et i + 1
  1163. * Les DENOMINATEURS DOIVENT ETRE FOURNIT EN ENTREE !
  1164. *----------------------------------------------------------------------
  1165.  
  1166. 'DEBPROC' GETERROR TlocP_n*CHPOINT TlocM_n*CHPOINT WlocP_n*CHPOINT
  1167. WlocM_n*CHPOINT TlocP*CHPOINT TlocM*CHPOINT WlocP*CHPOINT WlocM*CHPOINT
  1168. den_N*FLOTTANT den_T*FLOTTANT Klatin*FLOTTANT ;
  1169.  
  1170. * exclusion des points anciennement et actuellement extrémités
  1171. * de la fissure pour le calcul de l'erreur
  1172.  
  1173.  
  1174. *-> i + 1/2 : TlocP TlocM WlocP WlocM
  1175. TPnor = 'EXCO' 'FN ' TlocP;
  1176. TPta1 = 'EXCO' 'FT1 ' TlocP;
  1177. TMnor = 'EXCO' 'FN ' TlocM;
  1178. TMta1 = 'EXCO' 'FT1 ' TlocM;
  1179.  
  1180. WPnor = 'EXCO' 'UN ' WlocP;
  1181. WPta1 = 'EXCO' 'UT1 ' WlocP;
  1182. WMnor = 'EXCO' 'UN ' WlocM;
  1183. WMta1 = 'EXCO' 'UT1 ' WlocM;
  1184.  
  1185. *-> i + 1 : TlocP_n TlocM_n WlocP_n WlocM_n ( _n = new )
  1186. TPnor_n = 'EXCO' 'FN ' TlocP_n;
  1187. TPta1_n = 'EXCO' 'FT1 ' TlocP_n;
  1188. TMnor_n = 'EXCO' 'FN ' TlocM_n;
  1189. TMta1_n = 'EXCO' 'FT1 ' TlocM_n;
  1190.  
  1191. WPnor_n = 'EXCO' 'UN ' WlocP_n;
  1192. WPta1_n = 'EXCO' 'UT1 ' WlocP_n;
  1193. WMnor_n = 'EXCO' 'UN ' WlocM_n;
  1194. WMta1_n = 'EXCO' 'UT1 ' WlocM_n;
  1195.  
  1196.  
  1197. * ---> NUMERATEUR
  1198.  
  1199. * problème normal
  1200. numN = 'MAXIMUM' (
  1201. ((((TPnor_n '-' TPnor)**2)'+' ((TMnor_n '-' TMnor)**2)) '/' Klatin) '+'
  1202. ((((WPnor_n '-' WPnor)**2)'+' ((WMnor_n '-' WMnor)**2)) '*' Klatin));
  1203.  
  1204.  
  1205. * problème tangentiel
  1206. numT= 'MAXIMUM' (
  1207. ((((TPta1_n '-' TPta1)**2)'+' ((TMta1_n '-' TMta1)**2)) '/' Klatin) '+'
  1208. ((((WPta1_n '-' WPta1)**2)'+' ((WMta1_n '-' WMta1)**2)) '*' Klatin));
  1209.  
  1210.  
  1211. * ---> ERREUR
  1212.  
  1213. *-> erreur normale
  1214. etaN = numN '/' den_N;
  1215.  
  1216.  
  1217. *-> erreur tangentielle
  1218. etaT = numT '/' den_T;
  1219.  
  1220. *-> erreur totale : ComputeLatinErrorNew dans ELFE-3D
  1221. totalErr = ((etaN) '+' (etaT )) '**' 0.5;
  1222.  
  1223. *-> erreur totale : publi récente (= manuscrit Epierres)
  1224. *ler1 = 'PROG' (etaT**0.5) (etaN**0.5);
  1225. *totalerr = ('MAXIMUM' ler1)'**'(0.5);
  1226. *'SI' (etaT > etaN);
  1227. * totalErr = etaT;
  1228. * 'SINON';
  1229. * totalErr = etaN
  1230. * 'FINSI';
  1231.  
  1232. 'FINPROC' totalerr (etaT**0.5) (etaN**0.5);
  1233.  
  1234.  
  1235. *----------------------------------------------------------------------
  1236. *
  1237. * LOCSTEPI : Etape locale pour le premier par de temps
  1238. *
  1239. *----------------------------------------------------------------------
  1240.  
  1241. 'DEBPROC' LOCSTEPI Klatin*FLOTTANT f_crack*FLOTTANT MeshInt*MAILLAGE
  1242. TlocP*CHPOINT TlocM*CHPOINT WlocP*CHPOINT WlocM*CHPOINT;
  1243.  
  1244. meshloc = 'CHANGER' POI1 meshint;
  1245. wdiff= wlocm '-' wlocp;
  1246. tdiff= tlocm '-' tlocp;
  1247. k_tan = klatin;
  1248. k_nor = klatin;
  1249.  
  1250. TlocP = 1.*TlocP;
  1251. TlocM = 1.*TlocM;
  1252. WlocP = 1.*WlocP;
  1253. WlocM = 1.*WlocM;
  1254.  
  1255. TlocPnor = 'EXCO' 'FN ' TlocP;
  1256. TlocPta1 = 'EXCO' 'FT1 ' TlocP;
  1257. TlocMnor = 'EXCO' 'FN ' TlocM;
  1258. TlocMta1 = 'EXCO' 'FT1 ' TlocM;
  1259.  
  1260. WlocPnor = 'EXCO' 'UN ' WlocP;
  1261. WlocPta1 = 'EXCO' 'UT1 ' WlocP;
  1262. WlocMnor = 'EXCO' 'UN ' WlocM;
  1263. WlocMta1 = 'EXCO' 'UT1 ' WlocM;
  1264.  
  1265. * initialisation booleen pour détecter fissure complètement ouverte
  1266. ouv1 = vrai;
  1267.  
  1268. 'SI' (('DIME' ('EXTRAIRE' tlocp 'COMP')) 'EGA' 3);
  1269. TlocPta2 = 'EXCO' 'FT2 ' TlocP; TlocMta2 = 'EXCO' 'FT2 ' TlocM;
  1270. WlocPta2 = 'EXCO' 'UT2 ' WlocP; WlocMta2 = 'EXCO' 'UT2 ' WlocM;
  1271. TMta2 = 0. * TlocMta2;TPta2 = 0. * TlocPta2;
  1272. WMta2 = 0. * WlocMta2;WPta2 = 0. * WlocPta2;
  1273. B3D = vrai;
  1274. 'SINON';
  1275. B3D = faux;
  1276. 'FINSI';
  1277. *
  1278. * Cacul des indicateur de contact et de frottement initial
  1279. ConIni = 0.5 * (('EXCO' 'UN 'wdiff)'-'(('EXCO' 'FN' tdiff)'/'k_nor));
  1280. GliIni1 = 0.5*((K_tan*('EXCO' 'UT1 'wdiff))'-'('EXCO' 'FT1' tdiff));
  1281.  
  1282. 'SI' B3D;
  1283. GliIni2 = 0.5*((K_tan*('EXCO' 'UT2 'wdiff))'-'('EXCO' 'FT2' tdiff));
  1284. GliIni = ((gliIni1**2)'+'(gliIni2**2))'**'(0.5);
  1285. 'SINON';
  1286. GliIni = GliIni1;
  1287. 'FINSI' ;
  1288.  
  1289. * seuil pour le glissement
  1290. GSeuil = f_crack * ('ABS' (K_nor * ConIni));
  1291.  
  1292. * Initialisation des champs correspondant à l'itération i + 1/2
  1293. TMnor = 0. * TlocMnor; TMta1 = 0. * TlocMta1;
  1294. TPnor = 0. * TlocMnor; TPta1 = 0. * TlocPta1;
  1295. WMnor = 0. * WlocMnor; WMta1 = 0. * WlocMta1;
  1296. WPnor = 0. * WlocMnor; WPta1 = 0. * WlocPta1;
  1297.  
  1298. * initialisation des compteurs
  1299. Ncon = 0;
  1300. Nouv = 0;
  1301. *
  1302. * Loi de comportement local i --> i + 1/2 avec Contact with friction (coulomb)
  1303. *---
  1304. Nbn1 = 'NBNO' Meshloc;
  1305. idnode = 0;
  1306.  
  1307. 'REPETER' blo1 (Nbn1);
  1308. idnode = idnode '+' 1;
  1309. ptemp = Meshloc 'POIN' idnode;
  1310. ValCI = 'EXTRAIRE' ConIni 'SCAL' ptemp;
  1311.  
  1312. * Valeur à l'itération i
  1313. wpn = 'EXTRAIRE' WlocPnor 'SCAL' ptemp;
  1314. wmn = 'EXTRAIRE' WlocMnor 'SCAL' ptemp;
  1315. wpt1 = 'EXTRAIRE' WlocPta1 'SCAL' ptemp;
  1316. wmt1 = 'EXTRAIRE' WlocMta1 'SCAL' ptemp;
  1317.  
  1318. tpn = 'EXTRAIRE' TlocPnor 'SCAL' ptemp;
  1319. tmn = 'EXTRAIRE' TlocMnor 'SCAL' ptemp;
  1320. tpt1 = 'EXTRAIRE' TlocPta1 'SCAL' ptemp;
  1321. tmt1 = 'EXTRAIRE' TlocMta1 'SCAL' ptemp;
  1322.  
  1323. 'SI' B3D;
  1324. tpt2 = 'EXTRAIRE' TlocPta2 'SCAL' ptemp;
  1325. tmt2 = 'EXTRAIRE' TlocMta2 'SCAL' ptemp;
  1326. wpt2 = 'EXTRAIRE' WlocPta2 'SCAL' ptemp;
  1327. wmt2 = 'EXTRAIRE' WlocMta2 'SCAL' ptemp;
  1328. 'FINSI';
  1329.  
  1330. * Si ouverture : >EG 0. pour le cas ou on bloque les sauts passe par là ;-)
  1331. 'SI' (ValCI '>EG' 0.);
  1332. * 'MESSAGE' 'OUVERTURE pour le noeud' idnode;
  1333. Nouv = Nouv '+' 1;
  1334. * calcul des noueaux déplacemtents (t= 0)
  1335. chp1 = 'MANUEL' chpo ptemp 'SCAL' (wpn '-' (tpn '/' K_nor));
  1336. chp2 = 'MANUEL' chpo ptemp 'SCAL' (wmn '-' (tmn '/' K_nor));
  1337. chp3 = 'MANUEL' chpo ptemp 'SCAL' (wpt1 '-' (tpt1 '/' K_nor));
  1338. chp4 = 'MANUEL' chpo ptemp 'SCAL' (wmt1 '-' (tmt1 '/' K_nor));
  1339. WPnor = WPnor '+' chp1; WMnor = WMnor '+' chp2;
  1340. WPta1 = WPta1 '+' chp3; WMta1 = WMta1 '+' chp4;
  1341.  
  1342. 'SI' B3D;
  1343. chp5 = 'MANUEL' chpo ptemp 'SCAL' (wpt2 '-' (tpt2 '/' K_nor));
  1344. chp6 = 'MANUEL' chpo ptemp 'SCAL' (wmt2 '-' (tmt2 '/' K_nor));
  1345. WPta2 = WPta2 '+' chp5; WMta2 = WMta2 '+' chp6;
  1346. 'FINSI';
  1347.  
  1348. * Si contact
  1349. 'SINON';
  1350. Ncon = Ncon '+' 1;
  1351. * 'MESSAGE' 'CONTACT pour le noeud' idnode;
  1352. 'SI' ((idnode > 1) 'ET' (idnode < Nbn1)); ouv1 = faux; 'FINSI';
  1353.  
  1354. * problème normal
  1355. * w
  1356. chp1 = 'MANUEL' chpo ptemp 'SCAL' (0.5 * ((wpn '+' wmn ) '-'
  1357. ((tpn '+' tmn) '/' K_nor)));
  1358. WPnor = WPnor '+' chp1; WMnor = WMnor '+' chp1;
  1359. * t
  1360. chp2 = 'MANUEL' chpo ptemp 'SCAL' (K_nor * ValCI);
  1361. TPnor = TPnor '+' chp2; TMnor = TMnor '-' chp2;
  1362.  
  1363. * extraction seuil adhérence
  1364. ValGI = 'EXTRAIRE' GliIni 'SCAL' ptemp;
  1365. G = 'EXTRAIRE' GSeuil 'SCAL' ptemp;
  1366. ValGI1 = 'EXTRAIRE' GliIni1 'SCAL' ptemp;
  1367. si B3D;
  1368. ValGI2 = 'EXTRAIRE' GliIni2 'SCAL' ptemp;
  1369. 'FINSI';
  1370.  
  1371. * Adhérence
  1372. 'SI' (('ABS' ValGI) '<' G);
  1373.  
  1374. * t
  1375. chp31 = 'MANUEL' chpo ptemp 'SCAL' valGI1;
  1376. TPta1 = TPta1 '+' chp31; TMta1 = TMta1 '-' chp31;
  1377.  
  1378. * w
  1379. chp41='MANU' chpo ptemp'SCAL'(wpt1'+'((valGI1'-'tpt1)'/'K_tan));
  1380. WPta1 = WPta1 '+' chp41; WMta1 = WMta1 '+' chp41;
  1381.  
  1382. 'SI' B3D;
  1383. chp32 = 'MANUEL' chpo ptemp 'SCAL' valGI2;
  1384. chp42='MANU' chpo ptemp'SCAL'(wpt2'+'((valGI2'-'tpt2)'/'K_tan));
  1385. TPta2 = TPta2 '+' chp32; TMta2 = TMta2 '-' chp32;
  1386. WPta2 = WPta2 '+' chp42; WMta2 = WMta2 '+' chp42;
  1387. 'FINSI';
  1388.  
  1389. * Glissement
  1390. 'SINON';
  1391. * 'LISTE' ptemp;
  1392. * 'MESSAGE' '---> Glissement';
  1393.  
  1394. * t
  1395. 'SI' (valGI 'EGA' 0.);
  1396. * chp31 = 'MANUEL' chpo ptemp 'SCAL' 0.;
  1397. 'SINON' ;
  1398. chp31 = 'MANUEL' chpo ptemp 'SCAL'(G*valGI1'/' ('ABS' valGI));
  1399. 'FINSI';
  1400. TPta1 = TPta1 '+' chp31; TMta1 = TMta1 '-' chp31;
  1401.  
  1402.  
  1403. * w
  1404. chp41 = 'MANUEL' chpo ptemp 'SCAL'
  1405. (wpt1 '+'( ( (G*valGI1'/'('ABS' valGI)) '-' tpt1) '/' K_tan));
  1406. chp51 = 'MANUEL' chpo ptemp 'SCAL'
  1407. (wmt1 '+'(((-1.*G*valGI1'/'('ABS' valGI))'-'tmt1) '/' K_tan));
  1408. WPta1 = WPta1 '+' chp41; WMta1 = WMta1 '+' chp51;
  1409.  
  1410. 'SI' B3D;
  1411. 'SI' (valGI 'EGA' 0.);
  1412. chp32 = 'MANUEL' chpo ptemp 'SCAL' 0.;
  1413. 'SINON' ;
  1414. chp32 = 'MANUEL' chpo ptemp 'SCAL'(G*valGI2'/'('ABS' valGI));
  1415. 'FINSI';
  1416. chp32 = 'MANUEL' chpo ptemp 'SCAL'(G*valGI2'/'('ABS' valGI));
  1417. chp42 = 'MANUEL' chpo ptemp 'SCAL'
  1418. (wpt2 '+'( ( (G*valGI2'/'('ABS' valGI)) '-' tpt2) '/' K_tan));
  1419. chp52 = 'MANUEL' chpo ptemp 'SCAL'
  1420. (wmt2 '+'(((-1.*G*valGI2'/'('ABS' valGI))'-'tmt2) '/' K_tan));
  1421. TPta2 = TPta2 '+' chp32; TMta2 = TMta2 '-' chp32;
  1422. WPta2 = WPta2 '+' chp42; WMta2 = WMta2 '+' chp52;
  1423. 'FINSI';
  1424.  
  1425. 'FINSI';
  1426.  
  1427. 'FINSI';
  1428.  
  1429. 'FIN' blo1;
  1430.  
  1431. tpnor = 'NOMC' 'FN ' tpnor natu discret;
  1432. tpta1 = 'NOMC' 'FT1 ' tpta1 natu discret;
  1433. tmnor = 'NOMC' 'FN ' tmnor natu discret;
  1434. tmta1 = 'NOMC' 'FT1 ' tmta1 natu discret;
  1435. wpnor = 'NOMC' 'UN ' wpnor natu diffus;
  1436. wpta1 = 'NOMC' 'UT1 ' wpta1 natu diffus;
  1437. wmnor = 'NOMC' 'UN ' wmnor natu diffus;
  1438. wmta1 = 'NOMC' 'UT1 ' wmta1 natu diffus;
  1439.  
  1440. TlocP_a = TPta1 'ET' TPnor;
  1441. TlocM_a = TMta1 'ET' TMnor;
  1442. WlocP_a = WPta1 'ET' WPnor;
  1443. WlocM_a = WMta1 'ET' WMnor;
  1444.  
  1445. 'SI' B3D;
  1446. tpta2 = 'NOMC' 'FT2 ' tpta2 natu discret;
  1447. tmta2 = 'NOMC' 'FT2 ' tmta2 natu discret;
  1448. wpta2 = 'NOMC' 'UT2 ' wpta2 natu diffus;
  1449. wmta2 = 'NOMC' 'UT2 ' wmta2 natu diffus;
  1450. TlocP_a = TlocP_a 'ET' Tpta2;
  1451. TlocM_a = TlocM_a 'ET' Tmta2;
  1452. WlocP_a = WlocP_a 'ET' WPta2;
  1453. WlocM_a = WlocM_a 'ET' WPta2;
  1454. 'FINSI' ;
  1455.  
  1456. 'MESSAGE' ' ';
  1457. 'MESSAGE' 'Nombre de noeud en contact = ' Ncon;
  1458. 'MESSAGE' 'Nombre de noeud ouvert = ' Nouv;
  1459. 'MESSAGE' ' ';
  1460.  
  1461. 'FINPROC' TlocP_a TlocM_a WlocP_a WlocM_a;
  1462.  
  1463.  
  1464. ************************************************************************
  1465. * AFFICHAGE *
  1466. ************************************************************************
  1467.  
  1468. OPTI TRAC PSC EPTR 6 POTR HELVETICA_16;
  1469. GRAPH = vrai;
  1470. COMPLET = faux;
  1471. * opti debug 1;
  1472.  
  1473. * Gestion de l'affichage au cours d'un calcul*
  1474. Tdisplay = tabl ;
  1475. * Affichage des levels set en début de cycle
  1476. Tdisplay . levelset = GRAPH;
  1477. * Affichage du repère locale
  1478. Tdisplay . rep_local = GRAPH;
  1479. * Affichage des déplacement, déformés, contraintes, au cours de itérations de la LATIN
  1480. * (pour chaque pas de temps !)
  1481. Tdisplay . Iteration = faux;
  1482. * Affichage des pressions de contact sur déforméee locale
  1483. Tdisplay . PasTemps = GRAPH;
  1484. * Affichage de la convergence de la LATIN
  1485. * (pour chaque pas de temps !)
  1486. Tdisplay . convergence = GRAPH;
  1487. * Affichage des FICS
  1488. Tdisplay . fics = GRAPH;
  1489. * Affichage du champ de déplacement global
  1490. Tdisplay . dep_global = vrai;
  1491.  
  1492.  
  1493.  
  1494. ************************************************************************
  1495. * PARAMETRES DU CALCUL *
  1496. ************************************************************************
  1497.  
  1498. dim1 = 2;
  1499.  
  1500. *----------------- < Précision pour solveur non-linéaire
  1501. si COMPLET;
  1502. * eps1 = 1E-6;
  1503. eps1 = 1E-8;
  1504. sino;
  1505. eps1 = 1E-3;
  1506. fins;
  1507.  
  1508. *----------------- < Nombre d'itéraion max
  1509. Itmax = 1000;
  1510.  
  1511. *----------------- < Materiau
  1512. Ey1 = 2.06E11;
  1513. nu1 = 0.3;
  1514. rho1 = 7800.0;
  1515.  
  1516. *----------------- < Frottement des levres de la Fissure
  1517. f_crack = 0.1;
  1518.  
  1519. *----------------- < Chargement
  1520. f_ext = -1.0E9;
  1521.  
  1522. *----------------- < facteur d'amplification pour les tracés
  1523. * effort
  1524. * fac1 = 0;
  1525. fac1 = 0.05 /(('ABS' f_ext));
  1526.  
  1527. * déformée
  1528. fac2 = 0.01*('ABS' (Ey1 '/' f_ext));
  1529.  
  1530. *----------------- < Paramètre géométrique Hammouda
  1531. L = 0.4;
  1532. W = 0.2;
  1533. a = 0.02;
  1534. beta = 45.;
  1535. *----------------- < Maillage structure
  1536. N1 = 10;
  1537. * taille élément
  1538. den4 = a '/' N1;
  1539. *----------------- < Taille élément interface
  1540. den3 = den4 '/' 1;
  1541.  
  1542.  
  1543. ************************************************************************
  1544. * MAILLAGE *
  1545. ************************************************************************
  1546.  
  1547. *----------------- < Maillage structure 1
  1548.  
  1549. *----< surface fine de propagation
  1550. p2 = (W '+' (2 * a)) (W '+' (2*a));
  1551. p3 = (W '+' (2 * a)) (W '-' (2*a));
  1552. p4 = (W '-' (2 * a)) (W '+' (2*a));
  1553. p5 = (W '-' (2 * a)) (W '-' (2*a));
  1554.  
  1555. l5 = 'DROIT' p2 p3 'DINI' den4 'DFIN' den4;
  1556. l6 = 'DROIT' p3 p5 'DINI' den4 'DFIN' den4;
  1557. l7 = 'DROIT' p5 p4 'DINI' den4 'DFIN' den4;
  1558. l8 = 'DROIT' p4 p2 'DINI' den4 'DFIN' den4;
  1559. lfi1 = l5 'ET' l6 'ET' l7 'ET' l8;
  1560. s2 = 'DALLER' l5 l6 l7 l8 'PLAN';
  1561.  
  1562. *----< contour grossier
  1563. p0 = 0. 0.;
  1564. p1 = (2*W) 0.;
  1565. p9 = (2*W) L;
  1566. p8 = 0. L;
  1567. l1 = 'DROIT' p0 N1 p1 ;
  1568. l2 = 'DROIT' p1 (N1) p9 ;
  1569. l3 = 'DROIT' p9 N1 p8 ;
  1570. l4 = 'DROIT' p8 (N1) p0 ;
  1571. lco1 = l1 'ET' l2 'ET' l3 'ET' l4;
  1572. lto1 = lco1 'ET' lfi1;
  1573. s1 = 'SURFACE' lto1 'PLANE';
  1574. linbas1 = l1;
  1575. linhaut = l3;
  1576. ling = l2;
  1577. lind = l4;
  1578. s3 = s1 'ET' s2;
  1579.  
  1580. *----------------- < Maillage interface (fissure) (si 1 seule front,
  1581. * définir tip1 comme le front)
  1582. tip2 = (W '+' (a * ('COS' beta))) (W '+' (a * ('SIN' beta)));
  1583. tip1 = (W '-' (a * ('COS' beta))) (W '-' (a * ('SIN' beta)));
  1584. crack = 'DROIT' tip2 tip1 'DINI' den3 'DFIN' den3 'COULEUR' viol;
  1585. stot1 = s3 'ET' crack;
  1586. si (GRAPH); TRAC stot1 'TITR' 'Maillage'; fins;
  1587.  
  1588. * Creation des Levels Sets
  1589. *---------------
  1590. psi0 phi0 = PSIPHI s2 crack 'DEUX' tip1 tip2;
  1591. si( Tdisplay.levelset);
  1592. TRAC phi0 s2 'TITR' 'Level Set \f';
  1593. TRAC psi0 s2 'TITR' 'Level Set \y';
  1594. finsi;
  1595.  
  1596.  
  1597. ************************************************************************
  1598. * MODELE, ENRICHISSEMENT, MATÉRIAU, MATRICE, RELATIONS *
  1599. ************************************************************************
  1600.  
  1601. *----------------- < Modele
  1602. mod1a = MODE s1 MECANIQUE ELASTIQUE ISOTROPE;
  1603. mod1b = MODE s2 MECANIQUE ELASTIQUE ISOTROPE xq4r;
  1604. * Il faut d'abord donner le modèle X-FEM pour qu'il soit défini en
  1605. * premier dans le modèle totale
  1606. mod1 = mod1b 'ET' mod1a;
  1607.  
  1608. *----------------- < Enrichissement
  1609. Che1X = tabl;
  1610. * Che1X . 0 = TRIE mod1 psi0 phi0 saut;
  1611. Che1X . 0 = TRIE mod1 psi0 phi0;
  1612.  
  1613. * constructionsion des blocages des ddl X-fem non actifs dans
  1614. * les éléments de transition.
  1615. Rel1X = tabl;
  1616. Rel1X . 0 = RELA mod1;
  1617.  
  1618. *----------------- < Matériau
  1619. mat1 = MATE mod1 'YOUN' Ey1 'NU' nu1 'RHO' rho1;
  1620.  
  1621. *----------------- < Latin
  1622. Klatin = 1*(Ey1/('MESURE' crack)) ;
  1623.  
  1624. *----------------- < Stabilisation
  1625. xi = 1*(-1) '/' (Klatin);
  1626. *'OPTION' donn 5;
  1627. *-----------------------------< Creation des differentes matrices blocs
  1628. * du système
  1629. kuwt = RELA 'ACCRO' 'FAIBL' 'PENA' Klatin 'STAB' xi 'NOMU' mod1 Crack;
  1630. * multiplication par facteur 2 à cause de l'utilisation d'enrichissement
  1631. * saut sur l'interface
  1632. kuwt = 2 * kuwt;
  1633.  
  1634. *-----------------------------< Rigidité formulation XFEM 1 champ
  1635. k1 = 'RIGIDITE' mod1 mat1;
  1636.  
  1637. *-----------------------------< Condition aux limites
  1638. cl1 = (BLOQ 'DEPL' (linbas1 'POIN' 'INITIAL')) 'ET'
  1639. (BLOQ 'DEPL' (linbas1 'POIN' 'FINAL'));
  1640. cl1 = cl1 'ET' ('BLOQUE' 'DEPL' 'DIRECTION' (0. 1.) linbas1);
  1641.  
  1642. *-----------------------------< Matrice globale 1 champ
  1643. K1tot = k1 'ET' cl1 'ET' (Rel1X.0);
  1644.  
  1645. *-----------------------------< Matrice globale 3 champs
  1646. Ktot1 = K1tot 'ET' Kuwt;
  1647.  
  1648. *-----------------------------< Effort imposé
  1649. RHS0 = PRES 'MASS' mod1 (-1.*f_ext) linhaut;
  1650. si (GRAPH); TRAC (vect RHS0 FORC roug) stot1 'TITR' 'chargement'; fins;
  1651.  
  1652. *-----------------------------< Construction repère local
  1653. mod_int = 'MODELISER' crack 'MECANIQUE' 'ZCO2';
  1654. chn cht1 = LOCABASE crack ;
  1655. chn = -1.*chn;
  1656. cht1 = -1.*cht1;
  1657.  
  1658. *-----------------------------< Affichage repère local
  1659. vn = vect chn 'NX ' 'NY ' jaune;
  1660. vt1 = vect cht1 'T1X ' 'T1Y ' blanc;
  1661. vloc = vn 'ET' vt1;
  1662. 'SI' Tdisplay.rep_local;
  1663. 'TRACER' vloc (crack coul bleu);
  1664. 'FINSI';
  1665.  
  1666.  
  1667. ************************************************************************
  1668. * CALCUL *
  1669. ************************************************************************
  1670.  
  1671. *-----------------------------< Calcul avec un unique pas de temps,
  1672. * seule la solution du problème de rupture nous interresse ici
  1673. IdPasT = 1;
  1674.  
  1675.  
  1676. *-----------------------------< Resolution
  1677. 'SI' (IdPasT 'EGA' 1);
  1678. u0 Tplus Tmoins Wplus Wmoins err = LATIN Tdisplay crack
  1679. K1tot mod1 RHS0 IdPasT ktot1 kuwt cht1 chn ItMax eps1 klatin
  1680. f_crack s3;
  1681. 'SINON' ;
  1682. u0 Tplus Tmoins Wplus Wmoins err = LATIN Tdisplay crack
  1683. K1tot mod1 RHS0 IdPasT ktot1 kuwt cht1 chn ItMax eps1 klatin
  1684. f_crack v1 cht2 Tplus Tmoins Wplus Wmoins;
  1685. 'FINSI';
  1686.  
  1687.  
  1688. ************************************************************************
  1689. * POST-TRAITEMENTS *
  1690. ************************************************************************
  1691.  
  1692. *-----------------------------< Affichage de la pressions de contact
  1693. * sur la déformée
  1694. 'SI'(Tdisplay.PasTemps);
  1695. b = xfem reco u0 mod1;
  1696. def0 = defo s3 b fac2;
  1697. VTplus = vect Tplus forc fac1 rouge;
  1698. VTmoins = vect Tmoins forc fac1 vert;
  1699. def1 = defo wplus crack fac2 VTplus rouge;
  1700. def2 = defo wmoins crack fac2 VTmoins vert;
  1701. def3 = def1 'ET' def2 'ET' def0;
  1702. 'TRACER' def3 TITR 'deformee';
  1703. 'TRACER' (def1 et def2) TITR 'deformee';
  1704. 'FINSI';
  1705. *-----------------------------< Convergence
  1706. 'SI' (Tdisplay.convergence);
  1707. * Réglage de la fenêtre de tracé
  1708. min1 = ('MINIMUM' err);
  1709. max1 = ('MAXI' err);
  1710. abs2 = 'MAXI' ('PROG' ('DIME' err));
  1711. abs1 = 'PROG' 1. PAS 1. abs2;
  1712.  
  1713. * Erreur
  1714. evo1 = 'EVOL' 'MANUEL' abs1 err;
  1715. DESS evo1
  1716. 'GRIL' 'POIN' 'GRIS'
  1717. 'LOGX' 'XBOR' 1 abs2 'TITX' 'Iterations' POSX CENT
  1718. 'LOGY' 'YBOR' min1 max1 'TITY' 'Erreur' POSY CENT
  1719. TITR 'Convergence de la LATIN';
  1720. 'FINSI';
  1721.  
  1722.  
  1723. *-----------------------------< Ajout des termes de contacts dans le
  1724. * calcul des FICs
  1725. wsaut = wplus '-' wmoins;
  1726. lmow = 'EXTRAIRE' wsaut 'COMP';
  1727. mow1 = 'MOTS' 'AX ' 'AY ';
  1728. wsaut = 'NOMC' lmow mow1 wsaut;
  1729.  
  1730. tpglo = EXCO (tplus) (mots 'FX' 'FY') 'NOID';
  1731. lmo1 = 'EXTRAIRE' tpglo 'COMP';
  1732. lmo2 = 'MOTS' 'SMX ' 'SMY ';
  1733. tpglo = 'NOMC' lmo1 lmo2 tpglo;
  1734.  
  1735. *-----------------------------< Calcul des FICs : front 1
  1736. 'MESSAGE' ' ';'MESSAGE' '----------------------------';
  1737. 'MESSAGE' 'Calcul des FICs : front 1';
  1738. 'MESSAGE' '----------------------------';'MESSAGE' ' ';
  1739. KTAB1 = TABL;
  1740. KTAB1 . 'OBJECTIF' = MOT 'DECOUPLAGE';
  1741. KTAB1 . 'PSI' = psi0;
  1742. KTAB1 . 'PHI' = phi0;
  1743. KTAB1 . 'FRONT_FISSURE' = tip1;
  1744. KTAB1 . 'MODELE' = mod1;
  1745. KTAB1 . 'CARACTERISTIQUES' = mat1;
  1746. KTAB1 . 'SOLUTION_RESO' = u0;
  1747. KTAB1 . 'CHARGEMENTS_MECANIQUES' = RHS0;
  1748. KTAB1 . 'COUCHE' = 2;
  1749. KTAB1 . 'MODELE_FISSURE' = mod_int;
  1750. KTAB1 . 'DEPLACEMENT_FISSURE' = wsaut;
  1751. KTAB1 . 'PRESSION_FISSURE' = TPglo ;
  1752.  
  1753. G_THETA KTAB1;
  1754.  
  1755. K1G_1 = KTAB1 . 'RESULTATS' . 'I';
  1756. K2G_1 = KTAB1 . 'RESULTATS' . 'II';
  1757. GK1K2_1 = Ey1 * ((K1G_1**2) + (K2G_1**2)) ;
  1758. * * pas de valeur négative pour KI
  1759. * 'SI' (K1G_1 < 0.);
  1760. * K1G_1 = 0.;
  1761. * 'FINSI' ;
  1762.  
  1763. *-----------------------------< Calcul des FICs : front 2
  1764. 'MESSAGE' ' ';'MESSAGE' '----------------------------';
  1765. 'MESSAGE' 'Calcul des FICs : front 2';
  1766. 'MESSAGE' '----------------------------';'MESSAGE' ' ';
  1767. KTAB2 = TABL;
  1768. KTAB2 . 'OBJECTIF' = MOT 'DECOUPLAGE';
  1769. KTAB2 . 'PSI' = psi0;
  1770. KTAB2 . 'PHI' = phi0;
  1771. KTAB2 . 'FRONT_FISSURE' = tip2;
  1772. KTAB2 . 'MODELE' = mod1;
  1773. KTAB2 . 'CARACTERISTIQUES' = mat1;
  1774. KTAB2 . 'SOLUTION_RESO' = u0;
  1775. KTAB2 . 'CHARGEMENTS_MECANIQUES' = RHS0;
  1776. KTAB2 . 'COUCHE' = 2;
  1777. KTAB2 . 'MODELE_FISSURE' = mod_int;
  1778. KTAB2 . 'DEPLACEMENT_FISSURE' = wsaut;
  1779. KTAB2 . 'PRESSION_FISSURE' = TPglo ;
  1780.  
  1781. G_THETA KTAB2;
  1782. K1G_2 = KTAB2 . 'RESULTATS' . 'I';
  1783. K2G_2 = KTAB2 . 'RESULTATS' . 'II';
  1784. GK1K2_2 = Ey1 * ((K1G_2**2) + (K2G_2**2)) ;
  1785. * * pas de valeur négative pour KI
  1786. * 'SI' (K1G_2 < 0.);
  1787. * K1G_2 = 0.;
  1788. * 'FINSI' ;
  1789.  
  1790. 'SI' (Tdisplay.fics);
  1791. * Affichage de la nouvelle fissure
  1792. q7v1 = VECT 2.E-3 (KTAB1 . CHAMP_THET1) 'VERT';
  1793. q7v2 = VECT 2.E-3 (KTAB2 . CHAMP_THET1) 'BLEU';
  1794. 'TRACER' (q7v1 et q7v2) (s2 et crack) TITR 'champs \q';
  1795. 'MESSAGE' ' ';'MESSAGE' ' ';'MESSAGE' ' ';
  1796. 'MESSAGE' '----------------------------';
  1797. 'MESSAGE' 'K1 front1 = ' K1G_1;
  1798. 'MESSAGE' 'K2 front1 = ' K2G_1;
  1799. 'MESSAGE' '----------------------------';
  1800. 'MESSAGE' ' ';'MESSAGE' ' ';'MESSAGE' ' ';
  1801. 'MESSAGE' '----------------------------';
  1802. 'MESSAGE' 'K1 front2 = ' K1G_2;
  1803. 'MESSAGE' 'K2 front2 = ' K2G_2;
  1804. 'MESSAGE' '----------------------------';
  1805. 'FINSI';
  1806.  
  1807.  
  1808. 'SI' (Tdisplay . dep_global);
  1809. *-----------------------------<AFfichage champ de déplacement
  1810. u2 = enle u0 'LX ';
  1811. u2 = 'REDU' u2 s3;
  1812. u2 = xfem reco u2 mod1;
  1813. ux2 = 'NOMC' 'SCAL' ('EXCO' 'UX ' u2);
  1814. uy2 = 'NOMC' 'SCAL'('EXCO' 'UY ' u2);
  1815. uxy2 = 1.*(((ux2 **2)'+'(uy2 **2)));
  1816. 'TRACER' uxy2 s3 10;
  1817. 'FINSI';
  1818.  
  1819. mess ' ';mess ' ';mess ' ';
  1820. 'MESSAGE' '----------------------';
  1821. 'MESSAGE' 'FIN DES POST-TRAITEMENTS';
  1822.  
  1823.  
  1824.  
  1825. ************************************************************************
  1826. *** TEST DE BON FONCTIONNEMENT ***
  1827. ************************************************************************
  1828.  
  1829. * valeur de reference obtenue en 2013 :
  1830. K2ref = -11.5E7;
  1831. xtol = 0.10;
  1832. err21 = abs ((K2G_1 - K2ref) / K2ref);
  1833. err22 = abs ((K2G_2 - K2ref) / K2ref);
  1834. mess 'ecart relatif sur K2:' err21 err22;
  1835. si ((err21 >eg xtol) ou (err22 >eg xtol));
  1836. ERRE 5;
  1837. fins;
  1838.  
  1839. * OPTI donn 5 trac X;
  1840. fin ;
  1841.  
  1842.  
  1843.  
  1844.  
  1845.  
  1846.  
  1847.  
  1848.  
  1849.  
  1850.  
  1851.  
  1852.  
  1853.  
  1854.  
  1855.  
  1856.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales