Télécharger xfem04.dgibi

Retour à la liste

Numérotation des lignes :

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

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