Télécharger la plaquette Cast3M


fr en

Télécharger ktest-calp.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  ktest-calp.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ***********************************************************
  5. *                                                         *
  6. * ktest-calp.dgibi                                        *
  7. *                                                         *
  8. * Test de fonctionnement de l'opérateur CALP sur un       *
  9. * simple exemple de plaque carrée en flexion pure.        *
  10. * Plusieurs types de coques sont testées en même temps.   *
  11. *                                                         *
  12. * Auteur : Michel Bulik                                   *
  13. *                                                         *
  14. * Date : Août 1995                                        *
  15. *                                                         *
  16. ***********************************************************
  17.  
  18. *** Options ...
  19.  
  20.     opti dime 3 mode trid elem seg2 echo 1 ;
  21.  
  22. *** Paramètres ...
  23.  
  24. * longueur et largeur de la plaque 
  25.     L = 1. ;
  26.  
  27. * nombre de divisions le long de ses côtés
  28.     nbdiv = 10 ;
  29.  
  30. * module de YOUNG
  31.     valyoun = 2.e+11 ;
  32.  
  33. * coeff. de Poisson nul pour qu'il n'y ait pas 
  34. * de courbure anticlastique
  35.     valnu   = 0. ;
  36.  
  37. * épaisseur de la plaque
  38.     ep = 1.e-2 ;
  39.  
  40. * calcul linéaire ou non 
  41.     calclin = vrai ;
  42.  
  43. * affichage graphique ou non
  44.     graph   = faux ;
  45.  
  46. * types d'éléments testés
  47.     typdel = mots COQ4 DKT DST COQ3 ;
  48.  
  49. *** Points ...
  50.  
  51.      dens (L / nbdiv) ;
  52.  
  53.      p1 = 0 0 0 ;
  54.      p2 = L 0 0 ;
  55.  
  56.      vectz  = 0 0 L ;
  57.      vectxz = L 0 L ;
  58.  
  59. *** Lignes ...
  60.  
  61.     li1 = p1 d nbdiv p2 ;
  62.  
  63. *** Surface ...
  64.  
  65.     opti elem qua4 ;
  66.  
  67.     su1 = li1 tran nbdiv vectz ;
  68.  
  69. *** Début de la boucle sur les types d'éléments ...
  70.  
  71.     nbfois = dime typdel ;
  72.     repeter surelem nbfois ;
  73.  
  74.     elemact = extr typdel &surelem ;
  75.  
  76.     si(ega &surelem 2) ;
  77.         mess 'SU1 passe au type TRI3' ;
  78.         su1 = chan tri3 su1 ;
  79.     finsi ;
  80.  
  81.     mess '--------------------------------------' ;
  82.     mess 'On analyse actuellement l élément ' elemact ;
  83.  
  84. *** Modèle ...
  85.  
  86.     mo1 = mode su1 mecanique elastique elemact ;
  87.  
  88. *** Matériau ...
  89.  
  90.     ma1 = mate mo1 YOUN valyoun NU valnu ;
  91.  
  92.     ca1 = cara mo1 EPAI ep ;
  93.  
  94.     ma1 = ma1 et ca1 ;
  95.  
  96. *** Rigidité + Blocages ...
  97.  
  98.     si(calclin) ;
  99.        ri1 = rigi mo1 ma1 ;
  100.     finsi ;
  101.  
  102.     cl1 = (bloq uy uz li1) et
  103.           (bloq rx    li1) et
  104.           (bloq ux    p1 ) ;
  105.  
  106. *** Chargement ...
  107.  
  108. * un champ de forces obtenus par l'intégration avec les
  109. * fonctions de forme - on va le modifier pour obtenir
  110. * le chargement réparti uniformément sur le bord supérieur
  111. * de la plaque
  112.  
  113.     MOP = 'MODE' SU1 'CHARGEMENT' 'PRESSION' elemact ;
  114.     MAP = 'PRES' MOP 'PRES' 10000. ;
  115.     fopres = 'BSIG' MOP MAP ;
  116.  
  117.     mafo1 = su1 poin droi vectz vectxz ;
  118.     fo1 = exco fy (redu fopres mafo1) mx ;
  119. * moment appliqué par unité de longueur 
  120.     momtot = (abs (maxi (resu fo1))) / L ;
  121.     mess 'Moment appliqué par (par unité de longueur) = '
  122.           momtot ;
  123.  
  124. *** L'affichage éventuel du chargement ...
  125.  
  126.     vecchar = faux ;
  127.     si(exis (extr fo1 comp) FX) ;
  128.        titr 'Chargement : forces nodales' ;
  129.        vecfo1 = vect fo1 FX FY FZ (1 / (maxi (abs fo1))) ;
  130.        vecchar =  vrai ;
  131.     finsi ;
  132.     si (exis (extr fo1 comp) MX) ;
  133.        titr 'Chargement : couples nodaux' ;
  134.        vecfo1 = vect fo1 MX MY MZ (1 / (maxi (abs fo1))) ;
  135.        vecchar =  vrai ;
  136.     finsi ;
  137.     si(graph et vecchar) ;
  138.        trac vecfo1 su1 ;
  139.     finsi ;
  140.  
  141. *** Résolution ...
  142.  
  143.     si(calclin) ;
  144.        de1 = reso (ri1 et cl1) fo1 ;
  145.     sinon ;
  146.        tabnl = table ;
  147.        tabnl . MODELE = mo1 ;
  148.        tabnl . CARACTERISTIQUES = ma1 ;
  149.        tabnl . BLOCAGES_MECANIQUES = cl1 ;
  150.        lr1 = prog 0. 1. ;
  151.        ev1 = evol manu 't' lr1 'f(t)' lr1 ;
  152.        tabnl . TEMPS_CALCULES = prog 0.0 pas 1.e-1 1.0 ;
  153.        tabnl . TEMPS_SAUVES = lr1 ;
  154.        tabnl . GRANDS_DEPLACEMENTS = vrai ;
  155.        tabn1 . 'HYPOTHESE_DEFORMATIONS'='LINEAIRE' ;
  156.  
  157.        tabnl . MTOL = 1.e+30 ;
  158.  
  159.        pasapas tabnl ;
  160.  
  161.        de1 = tabnl . DEPLACEMENTS . 1 ;
  162.     finsi ;
  163.  
  164. *** Affichage éventuel de la déformée et des DDL ...
  165.  
  166.     defo0 = defo su1 de1 0. ;
  167.     titr 'La deformee' ;
  168.     defo1 = defo su1 de1 1. ;
  169.     si(graph) ;
  170.        trac (defo0 et defo1) ;
  171.  
  172. *       lescomp  = extr (redu de1 su1) comp ;
  173.        lescomp = mots UY RX ;
  174.        nbcomp = dime lescomp ;
  175.        repe surcomp nbcomp ;
  176.           lacomp = extr lescomp &surcomp ;
  177.           titr 'Element' elemact ' : Degre de liberte' 
  178.                lacomp ;
  179.           trac (exco lacomp de1) su1 ;
  180.        fin surcomp ;
  181.  
  182.     finsi ;
  183.  
  184. *** Contraintes ...
  185.  
  186.     si(calclin) ;
  187.        si1 = sigm mo1 ma1 de1 'LINE' ;
  188.     sinon ;
  189.        si1 = sigm mo1 ma1 de1 'QUAD' ;
  190.     finsi ;
  191.     si1 = rten si1 mo1 (1 0 0) (0 0 1) ;
  192.  
  193. *** L'affichage éventuel des contraintes ...
  194.  
  195.     si(graph) ;
  196.  
  197.        lescomp  = extr si1 comp ;
  198.        nbcomp = dime lescomp ;
  199.        repe surcom2 nbcomp ;
  200.           lacomp = extr lescomp &surcom2 ;
  201.           titr 'Element' elemact ' : Contrainte generalisee' 
  202.                 lacomp ;
  203.           trac (exco lacomp si1) mo1 ;
  204.        fin surcom2 ;
  205.  
  206.     finsi ;
  207.  
  208. *** Vérif de la valeur de la contrainte ...
  209.  
  210.     valmax = maxi (exco m22 si1) ;
  211.     valmin = mini (exco m22 si1) ;
  212.     si( (valmax - valmin) > 1.e-6 ) ;
  213.         mess 'Attention ! Résultat non homogène (M22)' ;
  214. *        valm22 = (valmin + valmax) / 2 ;
  215.          valm22 = (intg (exco m22 si1) mo1) / (L * L) ;
  216.     sinon ;
  217.         valm22 = valmax ;
  218.     finsi ;
  219.     ecartrel = 100 * ((valm22 - momtot) / momtot) ;
  220.     mess 'Ecart relatif (M22) = ' ecartrel '%' ;
  221.     si(ecartrel > 1.e-5) ; erreur 5 ; finsi ;
  222.  
  223. *** Application de CALP ...
  224.  
  225.     si1sup = calp si1    ca1   mo1 SUPE ;
  226.     si1inf = calp si1    ca1   mo1 INFE ;
  227.  
  228.     valtheor = (6 * momtot) / (ep * ep) ;
  229.  
  230.     valsup = (intg (exco smyy si1sup) mo1) / (L * L) ;
  231.     valinf = (intg (exco smyy si1inf) mo1) / (L * L) ;
  232.  
  233.     ecsup = 100 * ((valsup - valtheor) / valtheor) ;
  234.     mess 'Ecart relatif (peau supérieure) = ' ecsup '%' ;
  235.     si(ecsup > 1.e-5) ; erreur 5 ; finsi ;
  236.  
  237.     ecinf = 100 * ((valinf + valtheor) / valtheor) ;
  238.     mess 'Ecart relatif (peau inférieure) = ' ecinf '%' ;
  239.     si(ecinf > 1.e-5) ; erreur 5 ; finsi ;
  240.  
  241.  
  242. *** Résultats analytiques pour le déplacement
  243. *** et la rotation ...
  244.  
  245. * inertie par unité de longueur ...
  246.     I = (ep * ep * ep) / 12 ;
  247.  
  248.     depanal = (momtot * L * L) / (2 * valyoun * I) ;
  249.     depcal  = (maxi (exco uy de1)) ;
  250.     ecdep = 100 * ((depcal - depanal) / depanal) ;
  251.     mess 'Écart relatif (déplacement maxi) = ' ecdep '%' ;
  252.     si(ecdep > 0.5) ; erreur 5 ; finsi ;
  253.  
  254.     rotanal = (momtot * L) / (valyoun * I) ;
  255.     rotcal  = abs (mini (exco rx de1)) ;
  256.     ecrot   = 100 * ((rotcal - rotanal) / rotanal) ;
  257.     mess 'Écart relatif (rotation maxi) = ' ecrot '%' ;
  258.     si(ecrot > 2.0) ; erreur 5 ; finsi ;
  259.  
  260. *** Déformations généralisées ...
  261.  
  262.     si (calclin) ;
  263.        ep1 = epsi mo1 de1 ma1 'LINE' ;
  264.     sinon ;
  265.        ep1 = epsi mo1 de1 ma1 'QUAD' ;
  266.     finsi ;
  267.     ep1 = rten ep1 mo1 (1 0 0) (0 0 1) ;
  268.  
  269. *** et sur les peaux ...
  270.  
  271.     ep1sup = calp ep1    ca1   mo1 SUPE ;
  272.     ep1inf = calp ep1    ca1   mo1 INFE ;
  273.  
  274. *** Comparaison avec les valeurs analytiques ...
  275.  
  276. * RTTT
  277.     valmax = maxi (exco rttt ep1) ;
  278.     valmin = mini (exco rttt ep1) ;
  279.     si( (valmax - valmin) > 1.e-6 ) ;
  280.         mess 'Attention ! Résultat non homogène (RTTT)' ;
  281.         valrttt = (intg (exco rttt ep1) mo1) / (L * L) ;
  282.     sinon ;
  283.         valrttt = valmax ;
  284.     finsi ;
  285.     valanal = momtot / (valyoun * I) ;
  286.     ecartrel = 100 * ((valrttt - valanal) / valanal) ;
  287.     mess 'Ecart relatif (RTTT) = ' ecartrel '%' ;
  288.     si(ecartrel > 1.e-5) ; erreur 5 ; finsi ;
  289.  
  290. * EPYY (peau supérieure) 
  291.     valanal = (ep * valrttt) / 2 ;
  292.  
  293.     valmax = maxi (exco epyy ep1sup) ;
  294.     valmin = mini (exco epyy ep1sup) ;
  295.     si( (valmax - valmin) > 1.e-6 ) ;
  296.         mess 'Attention ! Résultat non homogène (EPYY sup)' ;
  297.         valepyys = (intg (exco epyy ep1sup) mo1) / (L * L) ;
  298.     sinon ;
  299.         valepyys = valmax ;
  300.     finsi ;
  301.     ecartyys = 100 * ((valepyys - valanal) / valanal) ;
  302.     mess 'Ecart relatif (EPYY sup) = ' ecartyys '%' ;
  303.     si(ecartyys > 1.e-5) ; erreur 5 ; finsi ;
  304.  
  305.     valmax = maxi (exco epyy ep1inf) ;
  306.     valmin = mini (exco epyy ep1inf) ;
  307.     si( (valmax - valmin) > 1.e-6 ) ;
  308.         mess 'Attention ! Résultat non homogène (EPYY inf)' ;
  309.         valepyyi = (intg (exco epyy ep1inf) mo1) / (L * L) ;
  310.     sinon ;
  311.         valepyyi = valmax ;
  312.     finsi ;
  313.     ecartyyi = 100 * ((valepyyi + valanal) / valanal) ;
  314.     mess 'Ecart relatif (EPYY inf) = ' ecartyyi '%' ;
  315.     si(ecartyyi > 1.e-5) ; erreur 5 ; finsi ;
  316.  
  317. *** Fin de la boucle sur les types d'éléments ...
  318.  
  319.     fin surelem ;
  320.  
  321. *** Bye ...
  322.  
  323.     fin ;
  324.  
  325.  
  326.  
  327.  
  328.  
  329.  
  330.  
  331.  
  332.  
  333.  
  334.  
  335.  
  336.  

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