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

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