Télécharger la plaquette Cast3M


fr en

Télécharger source2.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier : source2.dgibi
  2. * section : thermique
  3. *----------------------------------------------------------------------*
  4. *                           SOURCE2.DGIBI                              *
  5. *----------------------------------------------------------------------*
  6. *
  7. * Objet :
  8. * -------
  9. *
  10. * Verfication / validation d'un modele de source de chaleur.
  11. * Cas d'une SOURCE GAUSSIENNE SPHERIQUE.
  12. *
  13. * Description :
  14. * -------------
  15. * Comparaison des flux nodaux equivalents obtenus avec le modele a
  16. * ceux obtenus en construisant le champ "a la main", puis en l'integrant
  17. * avec l'operateur SOURCE.
  18. *
  19. * Type de calcul : Aucun
  20. * Mode de calcul : 2D PLAN, AXIS et 3D
  21. * Type d'element : TRI3 TRI6 QUA4 QUA8 TET4 TE10 PYR5 PY13 PRI6 PR15
  22. *                  CUB8 CU20
  23. * Objectifs      : Ecart relatif entre flux integres < 1.e-12
  24. *                  Ecart relatif entre la chaleur totale fournie (QTOT)
  25. *                  et la resultante sur le champ integre < 1.e-12
  26. *
  27. *----------------------------------------------------------------------*
  28. *
  29. * IG1 vrai : traces actives
  30. IG1      = faux ;
  31. *
  32. *------------------------ 2D ELEMENTS LINEAIRES -----------------------*
  33. *
  34. opti dime 2 elem qua4 ;
  35.  
  36. * Parametres du maillage :
  37. lo1      = 50.e-3 ;
  38. ep1      = 20.e-3 ;
  39.  
  40. * Parametres de la source :
  41. QT1      = 1.e3 ;
  42. OR1      = 30.e-3 ep1 ;
  43. RG1      = 5.e-3 ;
  44.  
  45. * Maillage :
  46. l1       = (0 0) droi 10 (0 ep1) ;
  47. s1       = l1    tran 50 (lo1 0) ;
  48. l2       = s1 cote 3 ;
  49. s2       = l2    tran 50 (lo1 ep1) ;
  50. s2       = s2 chan tri3 ;
  51. s0       = s1 et s2 ;
  52.  
  53. si IG1 ;
  54.   ltyp1    = s0 elem type ;
  55.   mot1     = mot 'ELEM =' ;
  56.   repe bm1 (dime ltyp1) ;
  57.     mot1     = chai mot1 ' ' (ltyp1 extr &bm1) ;
  58.   fin bm1 ;
  59.   mot1     = chai ' Maillage test en 2D (' mot1 ')' ;
  60.   trac s0 titr mot1 ;
  61. fins;
  62.  
  63. * MODE PLAN :
  64. * -----------
  65. *
  66. opti mode plan ;
  67. *
  68. moq1     = mode s0 thermique source gaussienne ;
  69. maq1     = mate moq1 'QTOT' QT1 'ORIG' OR1 'RGAU' RG1 ;
  70. sq1      = sour moq1 maq1 ;
  71.  
  72. qtmoq1   = sq1 resu maxi ;
  73.  
  74. mo1      = mode s0 thermique ;
  75. ma1      = mate mo1 k 25. ;
  76.  
  77. x1 y1    = s0 coor ;
  78. x1 y1    = (x1 - (or1 coor 1)) (y1 - (or1 coor 2)) ;
  79. x1 y1    = (chan cham x1 mo1 stresses) (chan cham y1 mo1 stresses) ;
  80. xqx1     = 4. / PI / RG1 / RG1 ;
  81. sqref1   = exp (-2. * ((x1 * x1) + (y1 * y1)) / RG1 / RG1) * QT1 * xqx1 ;
  82. sqref1   = sour mo1 sqref1 ;
  83.  
  84. si IG1 ;
  85.   trac sq1 S0 titr ' Flux nodaux du champ de source Gaussienne du MODELE ';
  86.   trac sqref1 s0 titr ' Flux nodaux du champ de source Gaussienne REFERENCE' ;
  87.   trac (sq1 - sqref1) s0 titr ' Difference des 2 champs' ;
  88. fins ;
  89.  
  90. err1     = (maxi abs (sq1 - sqref1)) / (maxi abs sqref1) ;
  91. list err1 ;
  92. err2     = (abs (qtmoq1 - qt1)) / (abs qt1) ;
  93. list err2 ;
  94. si ((err1 > 1.e-12) ou (err2 > 1.e-4)) ;
  95.   ERRE ' Probleme MODE SOURCE' ;
  96. fins ;
  97.  
  98. * MODE AXIS :
  99. * -----------
  100. *
  101. opti mode axis ;
  102. *
  103. moq1     = mode s0 thermique source gaussienne spherique ;
  104. maq1     = mate moq1 'QTOT' QT1 'ORIG' OR1 'RGAU' RG1 ;
  105. sq1      = sour moq1 maq1 ;
  106.  
  107. qtmoq1   = sq1 resu maxi ;
  108.  
  109. mo1      = mode s0 thermique ;
  110. ma1      = mate mo1 k 25. ;
  111.  
  112. x0 y0    = (or1 coor 1) (or1 coor 2) ;
  113. x1 y1    = s0 coor ;
  114. x1 y1    = (x1 - (or1 coor 1)) (y1 - (or1 coor 2)) ;
  115. x1 y1    = (chan cham x1 mo1 stresses) (chan cham y1 mo1 stresses) ;
  116. rc2      = 2. ** 0.5 ;
  117. rcpi1    = pi ** 0.5 ;
  118. xik1     = 0.25 * RG1 * RG1 * (exp (-2. * x0 * x0 / RG1 / RG1)) ;
  119. xik2     = 0.5 * RG1 * x0 * rcpi1 / rc2 ;
  120. xik3     = RG1 * x0 * (erf (x0 * rc2 / RG1)) * rcpi1 / (rc2 * rc2 * rc2) ;
  121. xqx1     = PI * RG1 * rcpi1 / rc2 ;
  122. xqx1     = xqx1 * (xik1 + xik2 + xik3) ;
  123. sqref1   = exp (-2. * ((x1 * x1) + (y1 * y1)) / RG1 / RG1) * QT1 / xqx1 ;
  124. sqref1   = sour mo1 sqref1 ;
  125.  
  126. si IG1 ;
  127.   trac sq1 S0 titr ' Flux nodaux du champ de source Gaussienne du MODELE ';
  128.   trac sqref1 s0 titr ' Flux nodaux du champ de source Gaussienne REFERENCE' ;
  129.   trac (sq1 - sqref1) s0 titr ' Difference des 2 champs' ;
  130. fins ;
  131.  
  132. err1     = (maxi abs (sq1 - sqref1)) / (maxi abs sqref1) ;
  133. list err1 ;
  134. err2     = (abs (qtmoq1 - qt1)) / (abs qt1) ;
  135. list err2 ;
  136. si ((err1 > 1.e-12) ou (err2 > 1.e-4)) ;
  137.   ERRE ' Probleme MODE SOURCE' ;
  138. fins ;
  139.  
  140. *---------------------- 2D ELEMENTS QUADRATIQUES ----------------------*
  141.  
  142. * Passage en EF quadratiques :
  143. s0       = s0 chan quad ;
  144.  
  145. si IG1 ;
  146.   ltyp1    = s0 elem type ;
  147.   mot1     = mot 'ELEM =' ;
  148.   repe bm1 (dime ltyp1) ;
  149.     mot1     = chai mot1 ' ' (ltyp1 extr &bm1) ;
  150.   fin bm1 ;
  151.   mot1     = chai ' Maillage test en 2D (' mot1 ')' ;
  152.   trac s0 titr mot1 ;
  153. fins;
  154.  
  155. * MODE PLAN :
  156. * -----------
  157. *
  158. opti mode plan ;
  159. *
  160. moq1     = mode s0 thermique source gaussienne spherique ;
  161. maq1     = mate moq1 'QTOT' QT1 'ORIG' OR1 'RGAU' RG1 ;
  162. sq1      = sour moq1 maq1 ;
  163.  
  164. qtmoq1   = sq1 resu maxi ;
  165.  
  166. mo1      = mode s0 thermique ;
  167. ma1      = mate mo1 k 25. ;
  168.  
  169. x1 y1    = s0 coor ;
  170. x1 y1    = (x1 - (or1 coor 1)) (y1 - (or1 coor 2)) ;
  171. x1 y1    = (chan cham x1 mo1 stresses) (chan cham y1 mo1 stresses) ;
  172. xqx1     = 4. / PI / RG1 / RG1 ;
  173. sqref1   = exp (-2. * ((x1 * x1) + (y1 * y1)) / RG1 / RG1) * QT1 * xqx1 ;
  174. sqref1   = sour mo1 sqref1 ;
  175.  
  176. si IG1 ;
  177.   trac sq1 S0 titr ' Flux nodaux du champ de source Gaussienne du MODELE ';
  178.   trac sqref1 s0 titr ' Flux nodaux du champ de source Gaussienne REFERENCE' ;
  179.   trac (sq1 - sqref1) s0 titr ' Difference des 2 champs' ;
  180. fins ;
  181.  
  182. err1     = (maxi abs (sq1 - sqref1)) / (maxi abs sqref1) ;
  183. list err1 ;
  184. err2     = (abs (qtmoq1 - qt1)) / (abs qt1) ;
  185. list err2 ;
  186. si ((err1 > 1.e-12) ou (err2 > 1.e-4)) ;
  187.   ERRE ' Probleme MODE SOURCE' ;
  188. fins ;
  189.  
  190. * MODE AXIS :
  191. * -----------
  192. *
  193. opti mode axis ;
  194. *
  195. moq1     = mode s0 thermique source gaussienne ;
  196. maq1     = mate moq1 'QTOT' QT1 'ORIG' OR1 'RGAU' RG1 ;
  197. sq1      = sour moq1 maq1 ;
  198.  
  199. qtmoq1   = sq1 resu maxi ;
  200.  
  201. mo1      = mode s0 thermique ;
  202. ma1      = mate mo1 k 25. ;
  203.  
  204. x0 y0    = (or1 coor 1) (or1 coor 2) ;
  205. x1 y1    = s0 coor ;
  206. x1 y1    = (x1 - (or1 coor 1)) (y1 - (or1 coor 2)) ;
  207. x1 y1    = (chan cham x1 mo1 stresses) (chan cham y1 mo1 stresses) ;
  208. rc2      = 2. ** 0.5 ;
  209. rcpi1    = pi ** 0.5 ;
  210. xik1     = 0.25 * RG1 * RG1 * (exp (-2. * x0 * x0 / RG1 / RG1)) ;
  211. xik2     = 0.5 * RG1 * x0 * rcpi1 / rc2 ;
  212. xik3     = RG1 * x0 * (erf (x0 * rc2 / RG1)) * rcpi1 / (rc2 * rc2 * rc2) ;
  213. xqx1     = PI * RG1 * rcpi1 / rc2 ;
  214. xqx1     = xqx1 * (xik1 + xik2 + xik3) ;
  215. sqref1   = exp (-2. * ((x1 * x1) + (y1 * y1)) / RG1 / RG1) * QT1 / xqx1 ;
  216. sqref1   = sour mo1 sqref1 ;
  217.  
  218. si IG1 ;
  219.   trac sq1 S0 titr ' Flux nodaux du champ de source Gaussienne du MODELE ';
  220.   trac sqref1 s0 titr ' Flux nodaux du champ de source Gaussienne REFERENCE' ;
  221.   trac (sq1 - sqref1) s0 titr ' Difference des 2 champs' ;
  222. fins ;
  223.  
  224. err1     = (maxi abs (sq1 - sqref1)) / (maxi abs sqref1) ;
  225. list err1 ;
  226. err2     = (abs (qtmoq1 - qt1)) / (abs qt1) ;
  227. list err2 ;
  228. si ((err1 > 1.e-12) ou (err2 > 1.e-4)) ;
  229.   ERRE ' Probleme MODE SOURCE' ;
  230. fins ;
  231.  
  232. *------------------------ 3D ELEMENTS LINEAIRES -----------------------*
  233. *
  234. opti dime 3 elem cub8 ;
  235. s0       = s0 chan line ;
  236.  
  237. * Maillage du volume :
  238. v0       = s0 volu tran 5 (0 0 ep1) ;
  239. x1       = v0 coor 1 ;
  240.  
  241. * Ajour de tetraedres :
  242. v3       = v0 elem appu stri (x1 poin supe (1.5 * lo1)) ;
  243. v0       = v0 diff v3 ;
  244. v0       = v0 et (v3 chan tet4) ;
  245.  
  246. * Ajout d'un element pyramide :
  247. el1      = (v0 elem cub8) elem 1 ;
  248. pe1      = el1 chan poi1 ;
  249. ps1      = pe1 poin 8 ;
  250. py1      = manu pyr5 (pe1 poin 1) (pe1 poin 2) (pe1 poin 3) (pe1 poin 4) ps1 ;
  251. py2      = manu pyr5 (pe1 poin 1) (pe1 poin 2) (pe1 poin 6) (pe1 poin 5) ps1 ;
  252. py3      = manu pyr5 (pe1 poin 2) (pe1 poin 3) (pe1 poin 7) (pe1 poin 6) ps1 ;
  253. py0      = py1 et py2 et py3 ;
  254. v0       = (v0 diff el1) et py0 ;
  255.  
  256. si IG1 ;
  257.   ltyp1    = v0 elem type ;
  258.   mot1     = mot 'ELEM =' ;
  259.   repe bm1 (dime ltyp1) ;
  260.     mot1     = chai mot1 ' ' (ltyp1 extr &bm1) ;
  261.   fin bm1 ;
  262.   mot1     = chai ' Maillage test en 3D (' mot1 ')' ;
  263.   trac v0 titr mot1 ;
  264. fins;
  265.  
  266. * Parametres de la source en 3D 
  267. OR1      = 30.e-3 (0.5 * ep1) ep1 ;
  268.  
  269. * Modeles :
  270. moq1     = mode V0 thermique source gaussienne ;
  271. maq1     = mate moq1 'QTOT' QT1 'ORIG' OR1 'RGAU' RG1 ;
  272. sq1      = sour moq1 maq1 ;
  273.  
  274. qtmoq1   = sq1 resu maxi ;
  275.  
  276. mo1      = mode V0 thermique ;
  277. ma1      = mate mo1 k 25. ;
  278.  
  279. x1 y1 z1 = V0 coor ;
  280. x1 y1 z1 = (x1 - (or1 coor 1)) (y1 - (or1 coor 2)) (z1 - (or1 coor 3)) ;
  281. x1 y1 z1 = (chan cham x1 mo1 stresses) (chan cham y1 mo1 stresses) (chan cham z1 mo1 stresses) ;
  282. xqx1     = (32. / (PI * PI * PI)) ** 0.5 / RG1 / RG1 / RG1 ;
  283. sqref1   = exp (-2. *((x1 * x1) + (y1 * y1) + (z1 * z1)) / RG1 / RG1) * QT1 * xqx1 ;
  284. sqref1   = sour mo1 sqref1 ;
  285.  
  286. si IG1 ;
  287.   trac sq1 V0 titr ' Flux nodaux du champ de source Gaussienne du MODELE ' coup or1 (0 1 0) (0 0 1) ;
  288.   trac sqref1 V0 titr ' Flux nodaux du champ de source Gaussienne REFERENCE' coup or1 (0 1 0) (0 0 1) ;
  289.   trac (sq1 - sqref1) V0 titr ' Difference des 2 champs' coup or1 (0 1 0) (0 0 1) ;
  290. fins ;
  291.  
  292. err1     = (maxi abs (sq1 - sqref1)) / (maxi abs sqref1) ;
  293. list err1 ;
  294. err2     = (abs (qtmoq1 - qt1)) / (abs qt1) ;
  295. list err2 ;
  296. si ((err1 > 1.e-12) ou (err2 > 2.e-4)) ;
  297.   ERRE ' Probleme MODE SOURCE' ;
  298. fins ;
  299.  
  300. *---------------------- 3D ELEMENTS QUADRATIQUES ----------------------*
  301. *
  302. v0       = v0 chan quad ;
  303.  
  304. si IG1 ;
  305.   ltyp1    = v0 elem type ;
  306.   mot1     = mot 'ELEM =' ;
  307.   repe bm1 (dime ltyp1) ;
  308.     mot1     = chai mot1 ' ' (ltyp1 extr &bm1) ;
  309.   fin bm1 ;
  310.   mot1     = chai ' Maillage test en 3D (' mot1 ')' ;
  311.   trac v0 titr mot1 ;
  312. fins;
  313.  
  314. * Modeles :
  315. moq1     = mode V0 thermique source gaussienne spherique ;
  316. maq1     = mate moq1 'QTOT' QT1 'ORIG' OR1 'RGAU' RG1 ;
  317. sq1      = sour moq1 maq1 ;
  318.  
  319. qtmoq1   = sq1 resu maxi ;
  320.  
  321. mo1      = mode V0 thermique ;
  322. ma1      = mate mo1 k 25. ;
  323.  
  324. x1 y1 z1 = V0 coor ;
  325. x1 y1 z1 = (x1 - (or1 coor 1)) (y1 - (or1 coor 2)) (z1 - (or1 coor 3)) ;
  326. x1 y1 z1 = (chan cham x1 mo1 stresses) (chan cham y1 mo1 stresses) (chan cham z1 mo1 stresses) ;
  327. xqx1     = (32. / (PI * PI * PI)) ** 0.5 / RG1 / RG1 / RG1 ;
  328. sqref1   = exp (-2. *((x1 * x1) + (y1 * y1) + (z1 * z1)) / RG1 / RG1) * QT1 * xqx1 ;
  329. sqref1   = sour mo1 sqref1 ;
  330.  
  331. si IG1 ;
  332.   trac sq1 V0 titr ' Flux nodaux du champ de source Gaussienne du MODELE ' ;
  333.   trac sqref1 V0 titr ' Flux nodaux du champ de source Gaussienne REFERENCE' ;
  334.   trac (sq1 - sqref1) V0 titr ' Difference des 2 champs' ;
  335. fins ;
  336.  
  337. err1     = (maxi abs (sq1 - sqref1)) / (maxi abs sqref1) ;
  338. list err1 ;
  339. err2     = (abs (qtmoq1 - qt1)) / (abs qt1) ;
  340. list err2 ;
  341. si ((err1 > 1.e-12) ou (err2 > 1.e-4)) ;
  342.   ERRE ' Probleme MODE SOURCE' ;
  343. fins ;
  344.  
  345. *-------------------------------- FIN ---------------------------------*
  346. fin ;
  347.  
  348.  
  349.  
  350.  
  351.  
  352.  
  353.  

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