Télécharger la plaquette Cast3M


fr en

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. * 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. * constructionsion des blocages des ddl X-fem non actifs dans 
  1613. * les éléments de transition.
  1614.  Rel1X = tabl;
  1615.  Rel1X . 0 = RELA mod1;
  1616.  
  1617. *----------------- < Matériau
  1618.  mat1 = MATE mod1 'YOUN' Ey1 'NU' nu1 'RHO' rho1;
  1619.  
  1620. *----------------- < Latin
  1621.  Klatin = 1*(Ey1/('MESURE' crack)) ;
  1622.  
  1623. *----------------- < Stabilisation
  1624.  xi = 1*(-1) '/' (Klatin);
  1625. *'OPTION' donn 5;
  1626. *-----------------------------< Creation des differentes matrices blocs 
  1627. *                               du système
  1628.  kuwt = RELA 'ACCRO' 'FAIBL' 'PENA' Klatin 'STAB' xi 'NOMU' mod1 Crack;
  1629. * multiplication par facteur 2 à cause de l'utilisation d'enrichissement 
  1630. * saut sur l'interface 
  1631.  kuwt = 2 * kuwt;
  1632.  
  1633. *-----------------------------< Rigidité formulation XFEM 1 champ
  1634.  k1 = 'RIGIDITE' mod1 mat1;
  1635.  
  1636. *-----------------------------< Condition aux limites 
  1637.  cl1 = (BLOQ  'DEPL' (linbas1 'POIN' 'INITIAL')) 'ET'
  1638.        (BLOQ  'DEPL' (linbas1 'POIN' 'FINAL'));
  1639.  cl1 = cl1 'ET' ('BLOQUE' 'DEPL' 'DIRECTION' (0. 1.) linbas1);
  1640.  
  1641. *-----------------------------< Matrice globale 1 champ
  1642.  K1tot  = k1  'ET' cl1 'ET' (Rel1X.0);
  1643.  
  1644. *-----------------------------< Matrice globale 3 champs
  1645.  Ktot1  = K1tot  'ET' Kuwt;
  1646.  
  1647. *-----------------------------< Effort imposé
  1648.  RHS0 = PRES 'MASS' mod1 (-1.*f_ext) linhaut;
  1649.  si (GRAPH); TRAC (vect RHS0 FORC roug) stot1 'TITR' 'chargement'; fins;
  1650.  
  1651. *-----------------------------< Construction repère local
  1652.  mod_int = 'MODELISER' crack 'MECANIQUE' 'ZCO2';
  1653.  chn cht1  = LOCABASE crack ;
  1654.  chn = -1.*chn;
  1655.  cht1 = -1.*cht1;
  1656.  
  1657. *-----------------------------< Affichage repère local
  1658.  vn = vect chn 'NX  ' 'NY  '  jaune; 
  1659.  vt1 = vect cht1 'T1X ' 'T1Y '  blanc; 
  1660.  vloc = vn 'ET' vt1;
  1661.  'SI' Tdisplay.rep_local;
  1662.    'TRACER' vloc (crack coul bleu);
  1663.  'FINSI';
  1664.  
  1665. LC1 = 'EXTR' mod_int 'FORC' ;
  1666. LC2 = 'EXTR' mod_int 'DEPL' ;
  1667. OPTI INCO LC1 LC2 ;
  1668. ************************************************************************
  1669. *                               CALCUL                                 *
  1670. ************************************************************************
  1671.  
  1672. *-----------------------------< Calcul avec un unique pas de temps, 
  1673. *     seule la solution du problème de rupture nous interresse ici
  1674.  IdPasT = 1;
  1675.  
  1676.  
  1677. *-----------------------------<  Resolution
  1678.  'SI' (IdPasT 'EGA' 1);
  1679.    u0 Tplus Tmoins Wplus Wmoins err = LATIN Tdisplay crack
  1680.    K1tot mod1 RHS0 IdPasT ktot1 kuwt cht1 chn ItMax eps1 klatin
  1681.    f_crack s3;
  1682.  'SINON' ;
  1683.    u0 Tplus Tmoins Wplus Wmoins err = LATIN Tdisplay crack
  1684.    K1tot mod1 RHS0 IdPasT ktot1 kuwt cht1 chn ItMax eps1 klatin
  1685.    f_crack v1 cht2 Tplus Tmoins Wplus Wmoins;
  1686.  'FINSI';
  1687.  
  1688.  
  1689. ************************************************************************
  1690. *                        POST-TRAITEMENTS                              *
  1691. ************************************************************************
  1692.  
  1693. *-----------------------------< Affichage de la pressions de contact 
  1694. *                               sur la déformée
  1695. 'SI'(Tdisplay.PasTemps);
  1696.    b = xfem reco u0 mod1;
  1697.    def0 =  defo s3 b fac2;
  1698.    VTplus = vect Tplus forc fac1 rouge;
  1699.    VTmoins = vect Tmoins forc fac1 vert;
  1700.    def1 = defo wplus crack fac2  VTplus rouge;                          
  1701.    def2 = defo wmoins crack fac2  VTmoins vert;
  1702.    def3 = def1 'ET' def2 'ET' def0;
  1703.    'TRACER' def3 TITR 'deformee';
  1704.    'TRACER' (def1 et def2) TITR 'deformee';
  1705.  'FINSI';
  1706. *-----------------------------< Convergence
  1707.  'SI' (Tdisplay.convergence);
  1708. *  Réglage de la fenêtre de tracé
  1709.    min1 = ('MINIMUM' err);  
  1710.    max1 = ('MAXI' err);
  1711.    abs2 = 'MAXI' ('PROG' ('DIME' err));
  1712.    abs1 = 'PROG' 1. PAS 1. abs2;
  1713.  
  1714. *  Erreur
  1715.    evo1 = 'EVOL' 'MANUEL' abs1 err;
  1716.    DESS  evo1 
  1717.    'GRIL' 'POIN' 'GRIS'
  1718.    'LOGX' 'XBOR' 1 abs2 'TITX' 'Iterations' POSX CENT
  1719.    'LOGY' 'YBOR' min1 max1 'TITY' 'Erreur' POSY CENT
  1720.    TITR 'Convergence de la LATIN';
  1721.  'FINSI';
  1722.  
  1723.  
  1724. *-----------------------------< Ajout des termes de contacts dans le 
  1725. *                               calcul des FICs
  1726.    wsaut = wplus '-' wmoins;
  1727.    lmow = 'EXTRAIRE' wsaut 'COMP';
  1728.    mow1 = 'MOTS' 'AX  ' 'AY  ';
  1729.    wsaut = 'NOMC' lmow mow1 wsaut;
  1730.  
  1731.    tpglo = EXCO (tplus) (mots 'FX' 'FY') 'NOID';
  1732.    lmo1 = 'EXTRAIRE' tpglo  'COMP';
  1733.    lmo2 = 'MOTS' 'SMX ' 'SMY  ';
  1734.    tpglo = 'NOMC' lmo1 lmo2 tpglo;
  1735.  
  1736. *-----------------------------< Calcul des FICs : front 1
  1737.    'MESSAGE' ' ';'MESSAGE' '----------------------------';
  1738.    'MESSAGE' 'Calcul des FICs : front 1';
  1739.    'MESSAGE' '----------------------------';'MESSAGE' ' ';
  1740.    KTAB1 = TABL;
  1741.    KTAB1 . 'OBJECTIF' = MOT 'DECOUPLAGE';
  1742.    KTAB1 . 'PSI' = psi0;
  1743.    KTAB1 . 'PHI' = phi0;
  1744.    KTAB1 . 'FRONT_FISSURE' = tip1;
  1745.    KTAB1 . 'MODELE' = mod1;
  1746.    KTAB1 . 'CARACTERISTIQUES' = mat1;
  1747.    KTAB1 . 'SOLUTION_RESO' = u0;
  1748.    KTAB1 . 'CHARGEMENTS_MECANIQUES' = RHS0;
  1749.    KTAB1 . 'COUCHE' = 2;
  1750.    KTAB1 . 'MODELE_FISSURE' = mod_int;
  1751.    KTAB1 . 'DEPLACEMENT_FISSURE' = wsaut;
  1752.    KTAB1 . 'PRESSION_FISSURE' = TPglo ;   
  1753.  
  1754.    G_THETA KTAB1;
  1755.  
  1756.    K1G_1 = KTAB1 . 'RESULTATS' . 'I';
  1757.    K2G_1 = KTAB1 . 'RESULTATS' . 'II';
  1758.    GK1K2_1 = Ey1 * ((K1G_1**2) + (K2G_1**2)) ;
  1759. * * pas de valeur négative pour KI 
  1760. *    'SI' (K1G_1 < 0.);
  1761. *       K1G_1 = 0.;
  1762. *    'FINSI' ;
  1763.  
  1764. *-----------------------------< Calcul des FICs : front 2
  1765.    'MESSAGE' ' ';'MESSAGE' '----------------------------';
  1766.    'MESSAGE' 'Calcul des FICs : front 2';
  1767.    'MESSAGE' '----------------------------';'MESSAGE' ' ';
  1768.    KTAB2 = TABL;
  1769.    KTAB2 . 'OBJECTIF' = MOT 'DECOUPLAGE';
  1770.    KTAB2 . 'PSI' = psi0;
  1771.    KTAB2 . 'PHI' = phi0;
  1772.    KTAB2 . 'FRONT_FISSURE' = tip2;
  1773.    KTAB2 . 'MODELE' = mod1;
  1774.    KTAB2 . 'CARACTERISTIQUES' = mat1;
  1775.    KTAB2 . 'SOLUTION_RESO' = u0;
  1776.    KTAB2 . 'CHARGEMENTS_MECANIQUES' = RHS0;
  1777.    KTAB2 . 'COUCHE' = 2;
  1778.    KTAB2 . 'MODELE_FISSURE' = mod_int;
  1779.    KTAB2 . 'DEPLACEMENT_FISSURE' = wsaut;
  1780.    KTAB2 . 'PRESSION_FISSURE' = TPglo ;   
  1781.  
  1782.    G_THETA KTAB2;
  1783.    K1G_2 = KTAB2 . 'RESULTATS' . 'I';
  1784.    K2G_2 = KTAB2 . 'RESULTATS' . 'II';
  1785.    GK1K2_2 = Ey1 * ((K1G_2**2) + (K2G_2**2)) ;
  1786. * * pas de valeur négative pour KI 
  1787. *    'SI' (K1G_2 < 0.);
  1788. *       K1G_2 = 0.;
  1789. *    'FINSI' ;
  1790.  
  1791.  'SI' (Tdisplay.fics);
  1792. * Affichage de la nouvelle fissure
  1793.    q7v1 = VECT 2.E-3 (KTAB1 . CHAMP_THET1) 'VERT';
  1794.    q7v2 = VECT 2.E-3 (KTAB2 . CHAMP_THET1) 'BLEU';
  1795.    'TRACER'  (q7v1 et q7v2) (s2 et crack) TITR 'champs \q';
  1796.    'MESSAGE' ' ';'MESSAGE' ' ';'MESSAGE' ' ';
  1797.    'MESSAGE' '----------------------------';
  1798.    'MESSAGE' 'K1 front1 = ' K1G_1;
  1799.    'MESSAGE' 'K2 front1 = ' K2G_1;
  1800.    'MESSAGE' '----------------------------';
  1801.    'MESSAGE' ' ';'MESSAGE' ' ';'MESSAGE' ' ';
  1802.    'MESSAGE' '----------------------------';
  1803.    'MESSAGE' 'K1 front2 = ' K1G_2;
  1804.    'MESSAGE' 'K2 front2 = ' K2G_2;
  1805.    'MESSAGE' '----------------------------';
  1806.  'FINSI';
  1807.  
  1808.  
  1809. 'SI' (Tdisplay . dep_global);
  1810. *-----------------------------<AFfichage champ de déplacement
  1811.  u2 = enle u0 'LX  ';
  1812.  u2 = 'REDU' u2 s3;
  1813.  u2 = xfem reco u2 mod1;
  1814.  ux2 =  'NOMC' 'SCAL' ('EXCO' 'UX  ' u2);
  1815.  uy2 =  'NOMC' 'SCAL'('EXCO' 'UY  ' u2);
  1816.  uxy2 = 1.*(((ux2 **2)'+'(uy2 **2)));
  1817.  'TRACER' uxy2 s3 10;
  1818. 'FINSI';
  1819.  
  1820. mess ' ';mess ' ';mess ' ';
  1821. 'MESSAGE' '----------------------';
  1822. 'MESSAGE' 'FIN DES POST-TRAITEMENTS';
  1823.  
  1824.  
  1825.  
  1826. ************************************************************************
  1827. ***                    TEST DE BON FONCTIONNEMENT                    ***
  1828. ************************************************************************
  1829.  
  1830. * valeur de reference obtenue en 2013 :
  1831. K2ref = -11.5E7;
  1832. xtol  = 0.10;
  1833. err21 = abs ((K2G_1 - K2ref) / K2ref);
  1834. err22 = abs ((K2G_2 - K2ref) / K2ref);
  1835. mess 'ecart relatif sur K2:' err21 err22;
  1836. si ((err21 >eg xtol) ou (err22 >eg xtol));
  1837.   ERRE 5;
  1838. fins;
  1839.  
  1840. * OPTI donn 5 trac X;
  1841. fin ;
  1842.  
  1843.  
  1844.  
  1845.  
  1846.  
  1847.  
  1848.  
  1849.  
  1850.  
  1851.  
  1852.  
  1853.  
  1854.  
  1855.  
  1856.  
  1857.  
  1858.  
  1859.  

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