Télécharger la plaquette Cast3M


fr en

Télécharger couplage_thermique.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  couplage_thermique.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. *
  5. *------------ Cas test couplage_thermique.dgibi --------------------
  6. *
  7. *
  8. * Tests des opérateurs de Castem-fluide
  9. *
  10. *
  11. *- Options générales
  12. *
  13. OPTI DIME 2 ELEM QUA8 ECHO 0 ;
  14. GRAPH = 'N' ;
  15. EPS0  = 1.D-6 ;
  16. *
  17. *-------------------------------------------------------------------
  18. * Couplage thermique entre deux domaines solides l'échange entre les
  19. * deux domaines étant réalisés via un coefficient d'échange : on
  20. * compare à l'état stationnaire solutions calculées et analytique.
  21. *-------------------------------------------------------------------
  22. * La géométrie est un barreau constitué de deux parties de longueur
  23. * respective L1 et L2 et de diffusivité thermique D1 et D2. On note
  24. * h le coefficient d'échange thermique à la frontière commune et
  25. * T1(x) (resp. T2(x)) la température dans le premier barreau (resp.
  26. * dans le deuxième).
  27. *
  28. * Chargement : Flux nul sur les frontières haute et basse (1D) et
  29. * température imposées aux deux extrémités (notées T1 et T2).
  30. *
  31. * Solution stationnaire : Soit Tg (resp. Td) la température au niveau
  32. * de l'interface entre les deux barreaux vue du premier barreau (resp.
  33. * du deuxième). A l'état stationnaire, l'égalité des gradients donne :
  34. *    Tg = T1 + (T2-T1)(a2*h)/(a1*h + a2*h + a1*a2)
  35. *    Td = T2 + (T1-T2)(a1*h)/(a1*h + a2*h + a1*a2)
  36. * avec a1=D1/L1 et a2=D2/L2.
  37. * On a alors
  38. *    T1(x) = T1 + (Tg-T1) x/L1 pour x in [0,L1]
  39. *    T2(x) = T2 + (Td-T2) (L1+L2-x)/L2 pour x in [L1,L1+L2]
  40. *
  41. * Application numérique : D1=D2=h=100. ; L1=L2=0.5 ; T1=40. ; T2=20.
  42. *    T1(x) = 40 - 10x
  43. *    T2(x) = 30 - 10x
  44. * Les valeurs considérées permettent de se placer dans la meme config
  45. * que celle prise par Trio_U dans son dossier de test de validation
  46. * (rapport DTP/SMTH/LMTL/96-22 de Pierre LEDAC).
  47. *-------------------------------------------------------------------
  48. * Auteur : F.Dabbene                                           01/99
  49. *-------------------------------------------------------------------
  50. * 06/06  : Correction d'une erreur dans le post-traitement
  51. *          Tracé de la solution via un mchaml au lieu d'un chpo
  52. *          Formulations (DISCR) LINE, MACRO et QUAF possibles
  53. *          Ajout du cas permanent
  54. *-------------------------------------------------------------------
  55. *
  56. *--------------------------------------- Début de la procédure VFRON
  57. DEBP VFRON ;
  58. ARGU RVX*TABLE ;
  59. *
  60. * Restriction des champs à la frontière
  61. *
  62. RV   = RVX . 'EQEX' ;
  63. TINC = RV . 'INCO' ;
  64. LPT1 = 'DOMA' RVX.'DOMZ' 'SOMMET' ;
  65. TINC . 'T1F' = 'KCHT' RVX.'DOMZ' 'SCAL' 'SOMMET' TINC.'T1' ;
  66. TINC . 'T2F' = 'KCHT' RVX.'DOMZ' 'SCAL' 'SOMMET' TINC.'T2' ;
  67. s2 ma1 = 'KOPS' 'MATRIK' ;
  68. 'RESPROC' s2 ma1 ;
  69. 'FINPROC' ;
  70. *--------------------------------------- Fin de la procédure VFRON
  71. *
  72. *= MAILLAGE
  73. *
  74. XMIN = 0. ; L1 = 0.5 ; X1 = XMIN + L1 ; L2 = L1 ; X2 = X1 + L2 ;
  75. YMIN = 0. ; DY = 0.5 ; Y1 = YMIN + DY ; Y2 = Y1 + DY ;
  76. NX = 20 ; NY = 1 ;
  77. *
  78. *------------------------------------------ Points du premier domaine
  79. P1 = XMIN YMIN ;
  80. P2 = X1   YMIN ;
  81. P3 = X1   Y2   ;
  82. P4 = XMIN Y2   ;
  83. *------------------------------------------ Points du deuxième domaine
  84. PD1 = P2 ;
  85. PD2 = X2  YMIN ;
  86. PD3 = X2  Y2   ;
  87. PD4 = P3       ;
  88. *------------------------------------------ Lignes du premier domaine
  89. P1P2 = P1 'DROIT' NX P2 ;
  90. P2P3 = P2 'DROIT' NY P3 ;
  91. P3P4 = P3 'DROIT' NX P4 ;
  92. P4P1 = P4 'DROIT' NY P1 ;
  93. *------------------------------------------ Lignes du deuxième domaine
  94. PD1PD2 = PD1 'DROIT' NX PD2 ;
  95. PD2PD3 = PD2 'DROIT' NY PD3 ;
  96. PD3PD4 = PD3 'DROIT' NX PD4 ;
  97. PD4PD1 = 'INVERSE' P2P3 ;
  98. *------------------------------------------ Maillages
  99. *
  100. * Afin d'éviter toute confusion de maillages, ces objets sont à écraser 
  101. * après la création des modèles (à remplacer par le maillage qui
  102. * sous-tend chaque modèle)
  103. DOM1 = 'DALLER' P1P2 P2P3 P3P4 P4P1 ;
  104. DOM2 = 'DALLER' PD1PD2 PD2PD3 PD3PD4 PD4PD1 ;
  105. DGAU = P4P1 ;
  106. DDRO = PD2PD3 ;
  107. DBA1 = P1P2 ;
  108. DBA2 = PD1PD2 ;
  109. FRON = P2P3 ;
  110. *
  111. *= Création des MODELES
  112. *
  113. DDOM1 = 'CHANGER' DOM1 'QUAF' ;
  114. DDOM2 = 'CHANGER' DOM2 'QUAF' ;
  115. DDOMT = DDOM1 'ET' DDOM2 ;
  116. DDBA1 = 'CHANGER' DBA1 'QUAF' ;
  117. DDBA2 = 'CHANGER' DBA2 'QUAF' ;
  118. DDGAU = 'CHANGER' DGAU 'QUAF' ;
  119. DDDRO = 'CHANGER' DDRO 'QUAF' ;
  120. DFRON = 'CHANGER' FRON 'QUAF' ;
  121. *
  122. DISCR = 'LINE' ;
  123. $DOM1 = 'MODELISER' DDOM1 'NAVIER_STOKES' DISCR ;
  124. $DOM2 = 'MODELISER' DDOM2 'NAVIER_STOKES' DISCR ;
  125. $DOMT = 'MODELISER' DDOMT 'NAVIER_STOKES' DISCR ;
  126. $DBA1 = 'MODELISER' DDBA1 'NAVIER_STOKES' DISCR ;
  127. $DBA2 = 'MODELISER' DDBA2 'NAVIER_STOKES' DISCR ;
  128. $DGAU = 'MODELISER' DDGAU 'NAVIER_STOKES' DISCR ;
  129. $DDRO = 'MODELISER' DDDRO 'NAVIER_STOKES' DISCR ;
  130. $FRON = 'MODELISER' DFRON 'NAVIER_STOKES' DISCR ;
  131. *
  132. DOM1 = 'DOMA' $DOM1 'MAILLAGE' ;
  133. DOM2 = 'DOMA' $DOM2 'MAILLAGE' ;
  134. DOMT = 'DOMA' $DOMT 'MAILLAGE' ;
  135. DGAU = 'DOMA' $DGAU 'MAILLAGE' ;
  136. DDRO = 'DOMA' $DDRO 'MAILLAGE' ;
  137. DBA1 = 'DOMA' $DBA1 'MAILLAGE' ;
  138. DBA2 = 'DOMA' $DBA2 'MAILLAGE' ;
  139. FRON = 'DOMA' $FRON 'MAILLAGE' ;
  140. *
  141. *= Données numériques
  142. *
  143. D1 = 100. ; D2 = D1 ; H = 100. ;
  144. T1 = 40.  ; T2 = 20. ; T12 = T1 + T2 / 2. ;
  145. *
  146. *= Solution analytique
  147. *
  148. DEN1 = D1/L1*H + (D2/L2*H) + (D2*D1/L1/L2) ;
  149. TG   = T2 - T1 * D2 / L2 * H /DEN1 + T1 ;
  150. TD   = T1 - T2 * D1 / L1 * H /DEN1 + T2 ;
  151. D1X = 'EXTR' ('EVOL' 'CHPO' DBA1 ('COOR' 1 DBA1)) 'ORDO' 1 ;
  152. D2X = 'EXTR' ('EVOL' 'CHPO' DBA2 ('COOR' 1 DBA2)) 'ORDO' 1 ;
  153. X1X = D1X / L1 ;
  154. NBPT = 'NBEL' ('CHAN' 'POI1' DBA1) ;
  155. X2X = D2X - ('PROG' NBPT*L1) / L2 ;
  156. T1X = (TG-T1)*X1X + ('PROG' NBPT*T1) ;
  157. T2X = (T2-TD)*X2X + ('PROG' NBPT*TD) ;
  158. EV1 = 'EVOL' 'ROUG' 'MANU' D1X T1X ;
  159. EV2 = 'EVOL' 'ROUG' 'MANU' D2X T2X ;
  160. EV3 = EV1 'ET' EV2 ;
  161. *
  162. *
  163. *
  164. *                  PREMIERE MODELISATION DU PROBLEME
  165. *
  166. *-------------------------------------------------------------------
  167. * L'état stationnaire est obtenu en tant que régime asymptotique d'un
  168. * transitoire. La résolution est implicite MAIS l'échange surfacique
  169. * à l'interface est traité par 'ECHI' en semi-explicite. De ce fait,
  170. * seule la solution convergée a un sens, le transitoire est faux.
  171. *-------------------------------------------------------------------
  172. *
  173. *= Description du problème (table EQEX)
  174. *
  175. RV = 'EQEX' $DOMT 'ITMA' 20 'ALFA' 1.
  176.             'OPTI' 'EF' 'IMPL'
  177.             'ZONE' $DOM1 'OPER' 'LAPN' D1 'INCO' 'T1'
  178.             'ZONE' $DOM2 'OPER' 'LAPN' D2 'INCO' 'T2' ;
  179. RV = 'EQEX' RV
  180.             'OPTI' 'EF' 'IMPL'
  181.             'ZONE' $FRON 'OPER' 'VFRON'
  182.             'ZONE' $FRON 'OPER' 'ECHI' H 'T2F' 'INCO' 'T1'
  183.             'ZONE' $FRON 'OPER' 'ECHI' H 'T1F' 'INCO' 'T2' ;
  184. RV = 'EQEX' RV
  185.             'OPTI' 'EF' 'IMPL' 'CENTREE'
  186.             'ZONE' $DOM1 'OPER' 'DFDT' 1. 'T1' 10. 'INCO' 'T1'
  187.             'ZONE' $DOM2 'OPER' 'DFDT' 1. 'T2' 10. 'INCO' 'T2' ;
  188. RV = 'EQEX' RV
  189.             'CLIM' 'T1' 'TIMP' DGAU T1
  190.                    'T2' 'TIMP' DDRO T2 ;
  191. *
  192. *= Initialisation de la table INCO et résolution
  193. *
  194. RV . 'INCO' . 'T1' = 'KCHT' $DOM1 SCAL SOMMET T12 ;
  195. RV . 'INCO' . 'T2' = 'KCHT' $DOM2 SCAL SOMMET T12 ;
  196. EXEC RV ;
  197. *
  198. *= Post-traitement de la solution calculée
  199. *
  200. EVC1 = 'EVOL' 'CHPO' DBA1 RV.'INCO'.'T1' ;
  201. EVC2 = 'EVOL' 'CHPO' DBA2 RV.'INCO'.'T2' ;
  202. CT1  = 'EXTR' EVC1 'ORDO' 1 ;
  203. CT2  = 'EXTR' EVC2 'ORDO' 1 ;
  204. EVC3 = ('EVOL' 'VERT' 'MANU' D1X CT1)
  205.   'ET' ('EVOL' 'VERT' 'MANU' D2X CT2) ;
  206. *
  207. *= Tracés
  208. *
  209. 'SI' ('EGA' GRAPH 'O') ;
  210. CHAM1 = 'CHAN' 'CHAM' RV.'INCO'.'T1' DOM1 ;
  211. CHAM2 = 'CHAN' 'CHAM' RV.'INCO'.'T2' DOM2 ;
  212. 'TRAC' $DOMT (CHAM1 'ET' CHAM2) ;
  213. *
  214. TAB1 = 'TABLE' ;
  215. TAB1 . 'TITRE' = 'TABLE' ;
  216. TAB1 . 'TITRE' . 2 = 'Solution castem (1)' ;
  217. TAB1 . 'TITRE' . 1 = '---------------' ;
  218. TAB1 . 'TITRE' . 4 = 'Solution exacte' ;
  219. TAB1 . 'TITRE' . 3 = '---------------' ;
  220. TAB1 . 1 = 'MARQ CROI NOLI' ;
  221. TAB1 . 2 = 'MARQ CROI NOLI' ;
  222. 'DESS' (EVC3 'ET' EV3) 'TITR' 'Temperature en fonction de x'
  223.                        'TITX' 'x' 'TITY' 'T'
  224.                        'LEGE' TAB1 ;
  225. 'FINSI' ;
  226. *
  227. *= Controle erreur
  228. *
  229. SOM1 = 'ABS' (('MAXI' ('INTG' EV3)) - ('MAXI' ('INTG' EVC3))) ;
  230. TEST = SOM1 '<' EPS0 ;
  231. 'SI' TEST ;
  232.      'ERRE' 0 ;
  233. 'SINO' ;
  234.      'MESS'  'Analytic Integration -  Computed = Difference '  ;
  235.      'MESS'  ('MAXI' ('INTG' EV3))'-' ('MAXI' ('INTG' EVC3)) '=' SOM1 ;
  236.      'ERRE' 5 ;
  237. 'FINSI' ;
  238. *
  239. *
  240. *
  241. *                  DEUXIEME MODELISATION DU PROBLEME
  242. *
  243. *-------------------------------------------------------------------
  244. * L'état stationnaire est obtenu en tant que régime asymptotique d'un
  245. * transitoire. La résolution est implicite, y compris au niveau de
  246. * l'interface (l'échange surfacique est traité par 'MDIA' et non pas
  247. * par 'ECHI'). Le transitoire a donc aussi un sens physique, le flux
  248. * échangé étant continu pendant le transitoire.
  249. *-------------------------------------------------------------------
  250. *
  251. *= Description du problème (table EQEX)
  252. *
  253. RV = 'EQEX' $DOMT 'ITMA' 5 'ALFA' 1.
  254.             'OPTI' EF IMPL
  255.             'ZONE' $DOM1 'OPER' 'LAPN' D1 'INCO' 'T1'
  256.             'ZONE' $DOM2 'OPER' 'LAPN' D2 'INCO' 'T2' ;
  257. RV = 'EQEX' RV
  258.             'OPTI' EFM1 IMPL
  259.             'ZONE' $FRON 'OPER' 'MDIA'      H  'INCO' 'T1' 'T1'
  260.             'ZONE' $FRON 'OPER' 'MDIA' (-1.*H) 'INCO' 'T2' 'T1'
  261.             'ZONE' $FRON 'OPER' 'MDIA'      H  'INCO' 'T2' 'T2'
  262.             'ZONE' $FRON 'OPER' 'MDIA' (-1.*H) 'INCO' 'T1' 'T2' ;
  263. RV = 'EQEX' RV
  264.             'OPTI' EF IMPL CENTREE
  265.             'ZONE' $DOM1 'OPER' 'DFDT' 1. 'T1' 10. 'INCO' 'T1'
  266.             'ZONE' $DOM2 'OPER' 'DFDT' 1. 'T2' 10. 'INCO' 'T2' ;
  267. RV = 'EQEX' RV
  268.             'CLIM' 'T1' 'TIMP' DGAU T1
  269.                    'T2' 'TIMP' DDRO T2 ;
  270. *
  271. *= Initialisation de la table INCO et résolution
  272. *
  273. RV . 'INCO' . 'T1' = 'KCHT' $DOM1 SCAL SOMMET T12 ;
  274. RV . 'INCO' . 'T2' = 'KCHT' $DOM2 SCAL SOMMET T12 ;
  275. EXEC RV ;
  276. *
  277. *= Post-traitement de la solution calculée
  278. *
  279. EVC1 = 'EVOL' 'CHPO' DBA1 RV.'INCO'.'T1' ;
  280. EVC2 = 'EVOL' 'CHPO' DBA2 RV.'INCO'.'T2' ;
  281. CT1  = 'EXTR' EVC1 'ORDO' 1 ;
  282. CT2  = 'EXTR' EVC2 'ORDO' 1 ;
  283. EVC3 = ('EVOL' 'VERT' 'MANU' D1X CT1)
  284.   'ET' ('EVOL' 'VERT' 'MANU' D2X CT2) ;
  285. *
  286. *= Tracés
  287. *
  288. 'SI' ('EGA' GRAPH 'O') ;
  289. CHAM1 = 'CHAN' 'CHAM' RV.'INCO'.'T1' DOM1 ;
  290. CHAM2 = 'CHAN' 'CHAM' RV.'INCO'.'T2' DOM2 ;
  291. 'TRAC' $DOMT (CHAM1 'ET' CHAM2) ;
  292. *
  293. TAB1 = 'TABLE' ;
  294. TAB1 . 'TITRE' = 'TABLE' ;
  295. TAB1 . 'TITRE' . 2 = 'Solution castem (2)' ;
  296. TAB1 . 'TITRE' . 1 = '---------------' ;
  297. TAB1 . 'TITRE' . 4 = 'Solution exacte' ;
  298. TAB1 . 'TITRE' . 3 = '---------------' ;
  299. TAB1 . 1 = 'MARQ CROI NOLI' ;
  300. TAB1 . 2 = 'MARQ CROI NOLI' ;
  301. 'DESS' (EVC3 'ET' EV3) 'TITR' 'Temperature en fonction de x'
  302.                        'TITX' 'x' 'TITY' 'T'
  303.                        'LEGE' TAB1 ;
  304. 'FINSI' ;
  305. *
  306. *= Controle erreur
  307. *
  308. SOM1 = 'ABS' (('MAXI' ('INTG' EV3)) - ('MAXI' ('INTG' EVC3))) ;
  309. TEST = SOM1 '<' EPS0 ;
  310. 'SI' TEST ;
  311.      'ERRE' 0 ;
  312. 'SINO' ;
  313.      'MESS'  'Analytic Integration -  Computed = Difference '  ;
  314.      'MESS'  ('MAXI' ('INTG' EV3))'-' ('MAXI' ('INTG' EVC3)) '=' SOM1 ;
  315.      'ERRE' 5 ;
  316. 'FINSI' ;
  317. *
  318. *
  319. *
  320. *                  TROISIEME MODELISATION DU PROBLEME
  321. *
  322. *-------------------------------------------------------------------
  323. * On resout directement le problème stationnaire en full implicite
  324. *-------------------------------------------------------------------
  325. *
  326. *= Description du problème (table EQEX)
  327. *
  328. RV = 'EQEX' $DOMT 'ITMA' 1 'ALFA' 1.
  329.             'OPTI' EF IMPL
  330.             'ZONE' $DOM1 'OPER' 'LAPN' D1 'INCO' 'T1'
  331.             'ZONE' $DOM2 'OPER' 'LAPN' D2 'INCO' 'T2' ;
  332. RV = 'EQEX' RV
  333.             'OPTI' EFM1 IMPL
  334.             'ZONE' $FRON 'OPER' 'MDIA'      H  'INCO' 'T1' 'T1'
  335.             'ZONE' $FRON 'OPER' 'MDIA' (-1.*H) 'INCO' 'T2' 'T1'
  336.             'ZONE' $FRON 'OPER' 'MDIA'      H  'INCO' 'T2' 'T2'
  337.             'ZONE' $FRON 'OPER' 'MDIA' (-1.*H) 'INCO' 'T1' 'T2' ;
  338. RV = 'EQEX' RV
  339.             'CLIM' 'T1' 'TIMP' DGAU T1
  340.                    'T2' 'TIMP' DDRO T2 ;
  341. *
  342. *= Initialisation de la table INCO et résolution
  343. *
  344. RV . 'INCO' . 'T1' = 'KCHT' $DOM1 SCAL SOMMET T12 ;
  345. RV . 'INCO' . 'T2' = 'KCHT' $DOM2 SCAL SOMMET T12 ;
  346. EXEC RV ;
  347. *
  348. *= Post-traitement de la solution calculée
  349. *
  350. EVC1 = 'EVOL' 'CHPO' DBA1 RV.'INCO'.'T1' ;
  351. EVC2 = 'EVOL' 'CHPO' DBA2 RV.'INCO'.'T2' ;
  352. CT1  = 'EXTR' EVC1 'ORDO' 1 ;
  353. CT2  = 'EXTR' EVC2 'ORDO' 1 ;
  354. EVC3 = ('EVOL' 'VERT' 'MANU' D1X CT1)
  355.   'ET' ('EVOL' 'VERT' 'MANU' D2X CT2) ;
  356. *
  357. *= Tracés
  358. *
  359. 'SI' ('EGA' GRAPH 'O') ;
  360. CHAM1 = 'CHAN' 'CHAM' RV.'INCO'.'T1' DOM1 ;
  361. CHAM2 = 'CHAN' 'CHAM' RV.'INCO'.'T2' DOM2 ;
  362. 'TRAC' $DOMT (CHAM1 'ET' CHAM2) ;
  363. *
  364. TAB1 = 'TABLE' ;
  365. TAB1 . 'TITRE' = 'TABLE' ;
  366. TAB1 . 'TITRE' . 2 = 'Solution castem (3)' ;
  367. TAB1 . 'TITRE' . 1 = '---------------' ;
  368. TAB1 . 'TITRE' . 4 = 'Solution exacte' ;
  369. TAB1 . 'TITRE' . 3 = '---------------' ;
  370. TAB1 . 1 = 'MARQ CROI NOLI' ;
  371. TAB1 . 2 = 'MARQ CROI NOLI' ;
  372. 'DESS' (EVC3 'ET' EV3) 'TITR' 'Temperature en fonction de x'
  373.                        'TITX' 'x' 'TITY' 'T'
  374.                        'LEGE' TAB1 ;
  375. 'FINSI' ;
  376. *
  377. *= Controle erreur
  378. *
  379. SOM1 = 'ABS' (('MAXI' ('INTG' EV3)) - ('MAXI' ('INTG' EVC3))) ;
  380. TEST = SOM1 '<' EPS0 ;
  381. 'SI' TEST ;
  382.      'ERRE' 0 ;
  383. 'SINO' ;
  384.      'MESS'  'Analytic Integration -  Computed = Difference '  ;
  385.      'MESS'  ('MAXI' ('INTG' EV3))'-' ('MAXI' ('INTG' EVC3)) '=' SOM1 ;
  386.      'ERRE' 5 ;
  387. 'FINSI' ;
  388. 'FIN' ;
  389.  
  390.  
  391.  
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  

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