Télécharger la plaquette Cast3M


fr en

Télécharger poudre2.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  poudre2.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. complet = faux;
  5. *
  6. *************************************************
  7. *  pour calcul complet mettre complet à : vrai;
  8. ***************************************************
  9. *
  10. ************************************************************
  11. *                                                          *
  12. *    Cas test pour la loi poudre_A                         *
  13. *                                                          *
  14. *    par Christophe DELLIS   (CEREM)                       *
  15. *                                                          *
  16. ************************************************************
  17. *  cas oedometrique
  18. *  un cylindre est densifié par mise en pression
  19. *  la densite initial est uniforme la densite finale
  20. *  est connue analytiquement
  21. *
  22. opti echo 0 ;
  23. *
  24. option dime 2 elem qua8 mode axis ;
  25. *
  26. *    parametres du maillage
  27. *
  28. rayon_1 = 10.0 ;
  29. *
  30. haute_1 = 15.0 ;
  31. *
  32. * definition des coordonnees des points
  33. *
  34. xa0 = 0.0 ;
  35. ya0 = 0.0 ;
  36. *
  37. xa1 = rayon_1 ;
  38. ya1 = ya0 ;
  39. *
  40. xa2 = xa1 ;
  41. ya2 = ya1 + haute_1 ;
  42. *
  43. xa3 = 0.0 ;
  44. ya3 = ya2 ;
  45. *
  46. * definition des densites
  47. *
  48. dena0 = 10.0 ;
  49. dena1 = 10.0 ;
  50. dena2 = 10.0 ;
  51. dena3 = 10.0 ;
  52. *
  53. * definition des points
  54. *
  55. densite dena0 ;
  56. a0 = xa0 ya0 ;
  57. *
  58. densite dena1 ;
  59. a1 = xa1 ya1 ;
  60. *
  61. densite dena2 ;
  62. a2 = xa2 ya2 ;
  63. *
  64. densite dena3 ;
  65. a3 = xa3 ya3 ;
  66. *
  67. * definition des lignes
  68. *
  69. si complet ;
  70.  n2 = 2;
  71. sinon;
  72.  n2= 2;
  73. finsi;
  74. la0a1 = a0 droi N2 a1 ;
  75. la1a2 = a1 droi N2 a2 ;
  76. la2a3 = a2 droi N2 a3 ;
  77. la3a0 = a3 droi N2 a0 ;
  78. *
  79. la3a2 = inve la2a3 ;
  80. la2a1 = inve la1a2 ;
  81. *
  82. * lignes pour les conditions limites
  83. *
  84. l_symz = la0a1 ;
  85. l_symx = la3a0 ;
  86. l_ext1 = la1a2 ;
  87. l_ext2 = la2a3 ;
  88. *
  89. * definition des surfaces
  90. *
  91. elmat1 = surf ( la0a1 et la1a2 et la2a3 et la3a0 ) plan ;
  92. *
  93. piece_1 = coul elmat1 bleu ;
  94. *
  95. tot = piece_1 ;
  96. elim (tot et a2) 0.001 ;
  97. *
  98. *
  99. ************************************************************
  100. *                                                          *
  101. *       calcul                                             *
  102. *                                                          *
  103. ************************************************************
  104. *
  105. T0TEMPE = 925.0 ;
  106. T1TEMPE = 925.0 ;
  107. T0TEMPS = 0.0 ;
  108. T1TEMPS = 1200.0 ;
  109. T2TEMPS = 2400.0 ;
  110. T3TEMPS = 3600.0 ;
  111. *
  112. *******************************************
  113. * Champ-point de temperature
  114. *******************************************
  115. *
  116. clth0001 = BLOQ 'T' (l_ext1 et l_ext2) ;
  117. p_bloth1 = DEPI clth0001 1. ;
  118. p_temps1 = prog T0TEMPS T1TEMPS T2TEMPS T3TEMPS ;
  119. p_theta1 = prog T0TEMPE T1TEMPE T1TEMPE T0TEMPE ;
  120. ev_ther1 = EVOL MANU TEMPS p_temps1 'T' p_theta1 ;
  121. ch_ther1 = CHAR TIMP p_bloth1 ev_ther1 ;
  122. *
  123. ************************************************************
  124. *                                                          *
  125. *       materiau                                           *
  126. *                                                          *
  127. ************************************************************
  128. P1MODEL = MODE piece_1 MECANIQUE ELASTIQUE
  129. VISCOPLASTIQUE POUDRE_A CONS DEUX;
  130. *
  131. V1YOUNG = 118.0E+03 ;
  132. V1NU = 0.33 ;
  133. V1ALPHA = 10.0E-10 ;
  134. V1A = 7.76e-15 ;
  135. V1N = 4.55 ;
  136. V1QSRT = 0.0 ;
  137. V1RHOR = 0.68 ;
  138. V1F0 = 0.1098 ;
  139. V1F1 = 1.0168 ;
  140. V1F2 = -0.2591 ;
  141. V1F3 = 0.0 ;
  142. V1F4 = 0.0 ;
  143. V1F5 = 0.66 ;
  144. V1C0 = 2.10851 ;
  145. V1C1 = 1.21767 ;
  146. V1C2 = -0.43081 ;
  147. V1C3 = 0.0 ;
  148. V1C4 = 0.0 ;
  149. V1C5 = 50.0 ;
  150. *
  151. P1MATER = MATER P1MODEL YOUNG V1YOUNG NU V1NU ALPH V1ALPHA TALP 20. TREF 0.
  152.  A V1A N V1N QSRT V1QSRT F0 V1F0 F1 V1F1 F2 V1F2 F3 V1F3
  153.  F4 V1F4 F5 V1F5 C0 V1C0 C1 V1C1 C2 V1C2 C3 V1C3 C4 V1C4
  154.  C5 V1C5 RHOR V1RHOR ;
  155. *
  156. P2MODEL = MODE piece_1 THERMIQUE ISOTROPE CONS UN ;
  157. *
  158. P2MATER = MATER P2MODEL 'C' 460 'K' 14.6 'RHO' 4.54 ;
  159. *
  160. ************************************************************
  161. *                                                          *
  162. *          chargement                                      *
  163. *                                                          *
  164. ************************************************************
  165. *
  166. p0bloca = bloque uz l_symz ;
  167. p1bloca = bloque ur (l_symx et l_ext1) ;
  168. *
  169. p_bloca = p0bloca et p1bloca ;
  170. *
  171. t0pres = 0.0 ;
  172. t1pres = 100.0 ;
  173. *
  174. li_tps = prog t0temps t1temps t2temps t3temps ;
  175. li_pres = prog t0pres t1pres t1pres t0pres ;
  176. ev_pres = evol manu temps li_tps pression li_pres ;
  177. *dess ev_pres titr 'Evolution du chargement en pression' nclk ;
  178. *
  179. p3force = pres mass p1model 1.0 l_ext2 ;
  180. titr 'Force appliquees' ;
  181. vec1 = vect p3force 'FR' 'FZ' 0.1;
  182. *trac vec1 piece_1 nclk ;
  183. *
  184. p3charg = char p3force ev_pres 'MECA';
  185. *
  186. p_charg = p3charg ;
  187. *
  188. ************************************************************
  189. *                                                          *
  190. *        calcul                                            *
  191. *                                                          *
  192. ************************************************************
  193. *
  194. p_model = p1model et p2model ;
  195. p_mater = p1mater et p2mater ;
  196. *
  197. * Initialisation de la table pour la procedure PASAPAS
  198. *
  199. si complet ;
  200. tfin = 2400.;
  201. sinon;
  202. tfin = 400;
  203. finsi;
  204. TAB1 = TABLE ;
  205. TAB1 . 'CHARGEMENT' = (p_charg et ch_ther1) ;
  206. TAB1 . 'MODELE' = p_model ;
  207. TAB1 . 'CARACTERISTIQUES' = p_mater ;
  208. TAB1 . 'BLOCAGES_MECANIQUES' = p_bloca ;
  209. TAB1 . 'BLOCAGES_THERMIQUES' = clth0001 ;
  210. TAB1 . 'TEMPS_CALCULES' = PROG 0.0
  211.         PAS 50.0  tfin;
  212. TAB1 . 'TEMPS_SAUVES' = PROG tfin;
  213. TAB1 . 'MAXITERATION' = 49 ;
  214. TAB1 . 'MAXISOUSPAS' = 500 ;
  215. TAB1 . 'MOVA' = 'MOT' 'VHOR' ;
  216. TAB1 . 'PRECISION' = 1.0E-4 ;
  217. TAB1 . 'PROCEDURE_THERMIQUE' = 'MOT' 'LINEAIRE' ;
  218. *
  219. PASAPAS TAB1 ;
  220. *
  221. ************************************************************
  222. *                                                          *
  223. *        Test de validation                                *
  224. *                                                          *
  225. ************************************************************
  226. *
  227. ind1 = index (TAB1 . 'VARIABLES_INTERNES') ;
  228. nb1 = dime ind1 ;
  229. *
  230. *   densite finale
  231. *
  232. vm1 = TAB1 . 'VARIABLES_INTERNES' .(ind1.nb1) ;
  233. vm2 = exco VHOR vm1 ;
  234. vm3 = redu vm2 piece_1 ;
  235. rhof1 = maxi vm3 ;
  236. *
  237. * Calcul theorique analytique de la densite finale :
  238. *
  239. si complet;
  240. rhof2 = 0.7269266 ;
  241. sinon;
  242. rhof2 = .68004;
  243. finsi;
  244. *
  245. * Calcul theorique analytique de la contrainte radiale finale :
  246. *
  247. si complet;
  248. sr1 = -77.91412;
  249. sinon;
  250. sr1 = -20.074 ;
  251. finsi;
  252. *
  253. * Lecture des deplacements du noeud au coin de l'echantillon :
  254. *
  255. ind2 = index (TAB1 . 'DEPLACEMENTS' ) ;
  256. nb2 = dime ind2 ;
  257. uz1 = (extr TAB1 .'DEPLACEMENTS'.(ind2.nb2) uz a2) ;
  258. *
  259. *
  260. * Lecture de la contrainte radiale du noeud au coin de l'echantillon :
  261. *
  262. ind1 = index (TAB1 . 'CONTRAINTES' ) ;
  263. nb1 = dime ind1 ;
  264. sr2 = maxi (exco TAB1 . 'CONTRAINTES'.(ind1.nb1) smrr ) ;
  265. list sr2 ;
  266. *
  267. * Calcul analytique des deplacements pour une densification isotrope :
  268. *
  269. si complet;
  270. uz2 = haute_1 * ((v1rhor/rhof1)**(1./1.) - 1.) ;
  271. sinon;
  272. uz2 = -4.00401E-03;
  273. finsi;
  274. list uz2 ;list uz1;
  275. *
  276. * Calcul des erreurs sur les deplacements et la densite :
  277. *
  278. erho1 = abs ((rhof1 - rhof2 )/rhof2);
  279. euz1 = abs ((uz1 - uz2) / uz2) ;
  280. esr1 = abs ((sr1 - sr2) / sr2) ;
  281. *
  282. * Ecriture des erreurs a l'ecran :
  283. *
  284. opti echo 0 ;
  285. sauter 3 lignes ;
  286. mess 'Calcul en matrice fermee : ' ;
  287. mess 'Erreur sur le deplacement vertical :  '  euz1 ;
  288. mess 'Erreur sur la contrainte radiale :    '  esr1 ;
  289. mess 'Erreur sur la densite finale:         '  erho1 ;
  290. sauter 3 lignes ;
  291.  
  292. *
  293.  
  294. si ( (euz1 + erho1 ) >eg 0.05 ) ;
  295.    mess '---------RESULTATS INCORRECTS-------------' ;
  296.    sauter 3 lignes ;
  297.    erreur 5;
  298. sinon ;
  299.    mess '---------RESULTATS CORRECTS-------------' ;
  300.    sauter 3 lignes ;
  301. finsi ;
  302. *
  303. fin ;
  304.  
  305.  
  306.  
  307.  
  308.  
  309.  
  310.  
  311.  
  312.  
  313.  
  314.  
  315.  

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