Télécharger darcy3EFMH.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : darcy3EFMH.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. *
  5. ************************** CAS TEST : darcy3.dgibi ******************
  6. *
  7. GRAPH = 'N' ;
  8. 'SAUT' 'PAGE' ;
  9. NITER = 40 ;
  10. *
  11. *-------------------------------------------------------------------
  12. * TEST DARCY3
  13. * CALCUL DARCY ORTHOTROPE 3D
  14. *
  15. * TEST TRANGEOL en EFMH
  16. * staionnaire comme limite d'un transitoire
  17. *
  18. * On effectue trois calculs sur un cube, maillé par des cubes
  19. * réguliers.
  20. * Les conditions aux limites varient suivant le cas considéré :
  21. * On impose le flux ou la charge sur les cotés du domaine.
  22. *
  23. * La solution analytique en charge est un polynome de degré un,
  24. * la vitesse est constante et la conductivité hydraulique orthotrope.
  25. * __ __
  26. * | |
  27. * | 1 0 0 |
  28. * K = | 0 3/4 0 |
  29. * | 0 0 1/2 |
  30. * |__ __|
  31. *
  32. * H(x,y,z) = -45 x -80 y -60z + 200.
  33. * V(x,y,z) = ( 45 ; 60 ; 30 )
  34. *
  35. * On s'attend à une précision de l'ordre de la précision machine.
  36. *
  37. *-------------------------------------------------------------------
  38. *
  39. 'SAUT' 'PAGE' ;
  40. *
  41. *- Options générales de calcul.
  42. *
  43. 'TITR' 'EFMH DARCY ORTHOTROPE 3D Lineaire : darcy3.dgibi' ;
  44. 'OPTI' 'DIME' 3 'ELEM' 'TET4' ;
  45. 'OPTI' 'ECHO' 1 ;
  46. *
  47. * ------------
  48. * = MAILLAGE =
  49. * ------------
  50. OEIL = 5.D0 6.D0 7.D0 ;
  51. VECX = 1.D0 0.D0 0.D0 ;
  52. VECY = 0.D0 1.D0 0.D0 ;
  53. VECZ = 0.D0 0.D0 1.D0 ;
  54. *
  55. ENX = 6 ;
  56. ENY = 6 ;
  57. ENZ = 6 ;
  58. DX = 1.D0 / ENX ;
  59. DY = 1.D0 / ENY ;
  60. DZ = 1.D0 / ENZ ;
  61. *
  62. *- Création des points
  63. *
  64. A0 = 0.D0 0.D0 0.D0 ;
  65. B0 = 1.D0 0.D0 0.D0 ;
  66. C0 = 1.D0 1.D0 0.D0 ;
  67. D0 = 0.D0 1.D0 0.D0 ;
  68. E0 = 0.D0 0.D0 1.D0 ;
  69. F0 = 1.D0 0.D0 1.D0 ;
  70. G0 = 1.D0 1.D0 1.D0 ;
  71. H0 = 0.D0 1.D0 1.D0 ;
  72. *
  73. *- Création des droites
  74. *
  75. AB = 'DROI' ENX A0 B0 ;
  76. AD = 'DROI' ENY A0 D0 ;
  77. AE = 'DROI' ENZ A0 E0 ;
  78. BC = 'DROI' ENY B0 C0 ;
  79. BF = 'DROI' ENZ B0 F0 ;
  80. CD = 'DROI' ENX C0 D0 ;
  81. CG = 'DROI' ENZ C0 G0 ;
  82. DH = 'DROI' ENZ D0 H0 ;
  83. EF = 'DROI' ENX E0 F0 ;
  84. EH = 'DROI' ENY E0 H0 ;
  85. FG = 'DROI' ENY F0 G0 ;
  86. GH = 'DROI' ENX G0 H0 ;
  87. *
  88. FE = 'INVE' EF ;
  89. EA = 'INVE' AE ;
  90. DC = 'INVE' CD ;
  91. HD = 'INVE' DH ;
  92. DA = 'INVE' AD ;
  93. HE = 'INVE' EH ;
  94. GF = 'INVE' FG ;
  95. FB = 'INVE' BF ;
  96. *
  97. *- Creation des faces du cube
  98. *
  99. SDRO = 'DALL' AB BF FE EA 'PLAN' ;
  100. SGAU = 'DALL' DC CG GH HD 'PLAN' ;
  101. SBAS = 'DALL' AB BC CD DA 'PLAN' ;
  102. SHAU = 'DALL' EF FG GH HE 'PLAN' ;
  103. SDEV = 'DALL' BC CG GF FB 'PLAN' ;
  104. SDER = 'DALL' AD DH HE EA 'PLAN' ;
  105. *
  106. *- Création maillage géométrique
  107. *
  108. ENXM = ENX + ENY + ENZ ;
  109. ELI0 = 1.D0 / ENXM / 10.D0 ;
  110. 'SI' ('EGA' ('VALEUR' 'ELEM') 'CUB8') ;
  111. CUBE1 = 'PAVE' SDER SDEV SBAS SHAU SGAU SDRO ;
  112. 'SINON' ;
  113. CUBE1 = 'VOLU' (SDER 'ET' SDEV 'ET' SBAS 'ET' SHAU
  114. 'ET' SGAU 'ET' SDRO) ;
  115. 'FINSI' ;
  116. QFTOT = CHANGE CUBE1 QUAF ;
  117. QFGAU = CHANGE SGAU QUAF ;
  118. QFDRO = CHANGE SDRO QUAF ;
  119. QFHAU = CHANGE SHAU QUAF ;
  120. QFBAS = CHANGE SBAS QUAF ;
  121. QFDEV = CHANGE SDEV QUAF ;
  122. QFDER = CHANGE SDER QUAF ;
  123. ELIM ELI0 (QFTOT ET QFGAU ET QFDRO ET QFHAU ET QFBAS ET QFDEV ET
  124. QFDER ) ;
  125. *
  126. *- Création maillage HYBRIDE et sous-objets (conditions aux limites)
  127. *
  128. MODHYB = MODE QFTOT 'DARCY' 'ANISOTROPE' ;
  129. MODGAU = MODE QFGAU 'DARCY' 'ANISOTROPE' ;
  130. MODDRO = MODE QFDRO 'DARCY' 'ANISOTROPE' ;
  131. MODHAU = MODE QFHAU 'DARCY' 'ANISOTROPE' ;
  132. MODBAS = MODE QFBAS 'DARCY' 'ANISOTROPE' ;
  133. MODDEV = MODE QFDEV 'DARCY' 'ANISOTROPE' ;
  134. MODDER = MODE QFDER 'DARCY' 'ANISOTROPE' ;
  135. C11 = 'DOMA' MODHYB 'VOLUME' ;
  136. CHYB1 = 'DOMA' MODHYB 'SURFACE' ;
  137. CHYB2 = 'DOMA' MODHYB 'NORMALE' ;
  138. CEGAU = 'DOMA' MODGAU 'CENTRE' ;
  139. CEDRO = 'DOMA' MODDRO 'CENTRE' ;
  140. CEHAU = 'DOMA' MODHAU 'CENTRE' ;
  141. CEBAS = 'DOMA' MODBAS 'CENTRE' ;
  142. CEDEV = 'DOMA' MODDEV 'CENTRE' ;
  143. CEDER = 'DOMA' MODDER 'CENTRE' ;
  144. *
  145. *- Solution analytique
  146. *
  147. XX YY ZZ = 'COOR' (DOMA MODHYB 'FACE' ) ;
  148. XXC YYC ZZC = 'COOR' (DOMA MODHYB 'CENTRE') ;
  149. *
  150. VKX = 1.D0 ;
  151. VKY = 0.75D0 ;
  152. VKZ = 0.5D0 ;
  153. VVKX = MANU 'CHPO' (DOMA MODHYB 'CENTRE') 'K11' 1.D0 ;
  154. VVKX = CHANGER ATTRIBUT VVKX 'NATU' DISCRET ;
  155. VVKY = MANU 'CHPO' (DOMA MODHYB 'CENTRE') 'K22' 0.75D0 ;
  156. VVKY = CHANGER ATTRIBUT VVKY 'NATU' DISCRET ;
  157. VVKZ = MANU 'CHPO' (DOMA MODHYB 'CENTRE') 'K33' 0.5D0 ;
  158. VVKZ = CHANGER ATTRIBUT VVKZ 'NATU' DISCRET ;
  159. AA = -45.D0 ;
  160. BB = -80.D0 ;
  161. CC = -60.D0 ;
  162. DD = 200.D0 ;
  163. AAA = -1.D0 * VKX * AA ;
  164. BBB = -1.D0 * VKY * BB ;
  165. CCC = -1.D0 * VKZ * CC ;
  166. *
  167. PANAF = (AA * XX) + (BB * YY) + (CC * ZZ) + DD ;
  168. PANAC = (AA * XXC) + (BB * YYC) + (CC * ZZC) + DD ;
  169. VANAC = 'MANU' 'CHPO' (DOMA MODHYB 'CENTRE') 3 'VX' AAA
  170. 'VY' BBB 'VZ' CCC ;
  171. VANAF = 'MANU' 'CHPO' (DOMA MODHYB 'FACE') 3 'VX' AAA
  172. 'VY' BBB 'VZ' CCC ;
  173. *
  174. * --------------
  175. * = RESOLUTION =
  176. * --------------
  177. *
  178. MATI3 = (NOMC 'K11' VVKX) et
  179. (NOMC 'K22' VVKY) et
  180. (NOMC 'K33' VVKZ) et
  181. (NOMC 'K21' (0.D0 * VVKX)) et
  182. (NOMC 'K31' (0.D0 * VVKX)) et
  183. (NOMC 'K32' (0.D0 * VVKX)) ;
  184. * ;
  185. *- Conditions aux limites
  186. *
  187. BBGAU = 'BLOQ' CEGAU 'TH' ;
  188. BBDRO = 'BLOQ' CEDRO 'TH' ;
  189. BBHAU = 'BLOQ' CEHAU 'TH' ;
  190. BBBAS = 'BLOQ' CEBAS 'TH' ;
  191. BBDEV = 'BLOQ' CEDEV 'TH' ;
  192. BBDER = 'BLOQ' CEDER 'TH' ;
  193. *
  194. *- TH imposée
  195. *
  196. TTIMP = 'REDU' PANAF CEGAU ;
  197. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH' ;
  198. EEGAU = 'DEPI' BBGAU TTIM2 ;
  199. TTIMP = 'REDU' PANAF CEDRO ;
  200. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH' ;
  201. EEDRO = 'DEPI' BBDRO TTIM2 ;
  202. TTIMP = 'REDU' PANAF CEBAS ;
  203. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH' ;
  204. EEBAS = 'DEPI' BBBAS TTIM2 ;
  205. TTIMP = 'REDU' PANAF CEHAU ;
  206. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH' ;
  207. EEHAU = 'DEPI' BBHAU TTIM2 ;
  208. TTIMP = 'REDU' PANAF CEDEV ;
  209. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH' ;
  210. EEDEV = 'DEPI' BBDEV TTIM2 ;
  211. TTIMP = 'REDU' PANAF CEDER ;
  212. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH' ;
  213. EEDER = 'DEPI' BBDER TTIM2 ;
  214. *
  215. *- Flux imposé
  216. *
  217. FLDRO = -60.D0 * (nomc 'FLUX' (doma moddro VOLUME));
  218. FLGAU = 60.D0 * (nomc 'FLUX' (doma modgau VOLUME));
  219. FLHAU = 30.D0 * (nomc 'FLUX' (doma modhau VOLUME));
  220. FLBAS = -30.D0 * (nomc 'FLUX' (doma modbas VOLUME));
  221. FLDEV = 45.D0 * (nomc 'FLUX' (doma moddev VOLUME));
  222. FLDER = -45.D0 * (nomc 'FLUX' (doma modder VOLUME));
  223. *
  224. *
  225. *- Assemblage et résolution en TH
  226. *
  227.  
  228.  
  229. h_lim = 'RESOUD' (BBGAU 'ET' BBDEV)
  230. (EEGAU 'ET' EEDEV) ;
  231. h_lim = 'EXCO' 'TH' h_lim;
  232. CHCLIM = TABLE;
  233. CHCLIM . 'NEUMANN' = ('NOMC' 'I35' ( FLHAU 'ET' FLBAS
  234. 'ET' FLDRO 'ET' FLDER));
  235. CHCLIM . 'DIRICHLET' = 'NOMC' 'I35' h_lim;
  236.  
  237.  
  238. GEOL1 = TABLE;
  239. GEOL1 . 'CONCENTRATION' = 'NOMC' 'I35' (0.D0 * PANAC) ;
  240. GEOL1 . 'LUMP' = FAUX ;
  241. GEOL1 . 'TYPDISCRETISATION' = 'EFMH' ;
  242. GEOL1 . 'THETA_DIFFUSION' = 1.0D0 ;
  243. GEOL1 . 'THETA_CONVECTION' = 1.0D0 ;
  244. GEOL1 . 'DECENTREMENT' = FAUX ;
  245. GEOL1 . 'DELTAT' = 10.D0 / NITER ;
  246. GEOL1 . 'DIFFUSIVITE' = mati3 ;
  247. GEOL1 . 'SOLVEUR' = 3 ;
  248. GEOL1 . 'PRECONDITIONNEUR' = 3 ;
  249. GEOL1 . 'POROSITE' = 'MANU' 'CHPO'
  250. (doma MODHYB CENTRE) 'CK' 1.D0 ;
  251. GEOL1 . 'CLIMITES' = CHCLIM ;
  252. GEOL1 . 'RECALCUL' = VRAI ;
  253. *
  254. GEOL1 GEOL2 = TRANGEOL Modhyb GEOL1;
  255.  
  256.  
  257.  
  258. * boucle transitoire
  259.  
  260.  
  261.  
  262. REPE blocc NITER ;
  263. mess 'boucle' &blocc;
  264. GEOL1 . 'RECALCUL' = FAUX ;
  265. GEOL1 GEOL2 = TRANGEOL Modhyb GEOL1 GEOL2;
  266. FIN blocc;
  267.  
  268.  
  269.  
  270.  
  271. CHTER3 = 'NOMC' 'TH' GEOL2 . 'TRACE_CONC';
  272. PCEN3 = 'NOMC' 'H' GEOL1 . 'CONCENTRATION' ;
  273. QFACE3 = 'NOMC' 'FLUX' GEOL1 . 'FLUXDIFF' ;
  274.  
  275. *
  276. *- Calcul de V
  277. *
  278. VCENT3 = 'HVIT' MODHYB QFACE3 ;
  279. QFACE3 = 'EXCO' QFACE3 'FLUX' 'SCAL' ;
  280. VFACE3 = QFACE3 * CHYB2 / CHYB1 ;
  281. *
  282. * -----------------
  283. * = Calcul ERREUR =
  284. * -----------------
  285. * ERReur relative en Trace de charge TH aux faces des éléments
  286. * ERReur relative en charge H au centre des éléments
  287. * Erreur relative sur la vitesse au centre des éléments
  288. *
  289. ERRTP3 = 'EXCO' CHTER3 'TH' 'SCAL' ;
  290. ERRTP3 = ERRTP3 - PANAF / PANAF ;
  291. ERRTP3 = 'ABS' ERRTP3 ;
  292. ERRP3 = 'EXCO' PCEN3 'H' 'SCAL' ;
  293. ERRP3 = ERRP3 - PANAC / PANAC ;
  294. ERRP3 = 'ABS' ERRP3 ;
  295. *
  296. MOT1 = 'MOTS' 'VX' 'VY' 'VZ' ;
  297. VDVD = 'PSCA' VANAC VANAC MOT1 MOT1 ;
  298. *
  299. VD3 = VANAC - VCENT3 ;
  300. VC3 = 'PSCA' VD3 VD3 MOT1 MOT1 ;
  301. SDC3 = 'ABS' ( VC3 / VDVD ) ;
  302. SDC3 = SDC3 '**' 0.5 ;
  303. *
  304. * -------------------
  305. * = Tracé resultats =
  306. * -------------------
  307. 'SI' ('NEG' GRAPH 'N') ;
  308. *
  309. *- Transformation des quantités aux centres en MCHAML constant.
  310. *
  311. ERRP3 = 'KCHA' MODHYB 'CHAM' ERRP3 ;
  312. SDS3 = 'KCHA' MODHYB 'CHAM' SDC3 ;
  313. *
  314. * L'erreur relative sur la charge au centre
  315. * L'erreur relative sur la Vitesse au centre
  316. *
  317. 'TITR' 'darcy3/3 : Erreur relative sur la charge' ;
  318. 'TRAC' MODHYB ERRP3 ;
  319. 'TITR' 'darcy3/3 : Erreur relative sur la vitesse' ;
  320. 'TRAC' MODHYB SDS3 ;
  321. *
  322. 'FINSI' ;
  323. *
  324. * -------------------
  325. * = Gestion ERREURS =
  326. * -------------------
  327. MAXTP3 = 'MAXI' ERRTP3 ;
  328. MAXP3 = 'MAXI' ERRP3 ;
  329. MAXV3 = 'MAXI' SDC3 ;
  330. *
  331. 'SAUT' 'PAGE' ;
  332. 'SAUT' 2 'LIGNE' ;
  333. 'MESS' ' ERREURS RELATIVES ' ;
  334. 'SAUT' 1 'LIGNE' ;
  335. 'MESS' ' cas test TH H V' ;
  336. 'SAUT' 1 'LIGNE' ;
  337. 'MESS' ' numero3 ' maxtp3 ' ' maxp3 ' ' maxv3 ;
  338. 'SAUT' 2 'LIGNE' ;
  339. *
  340. EPS0 = 1.D-7 ;
  341. LOG3 = MAXTP3 > EPS0 ;
  342. LOG6 = MAXP3 > EPS0 ;
  343. LOG9 = MAXV3 > EPS0 ;
  344. LTP0 = LOG3 ;
  345. LP0 = LOG6 ;
  346. LV0 = LOG9 ;
  347. L0 = LTP0 'OU' LP0 'OU' LV0 ;
  348. 'SI' ( L0 ) ;
  349. 'ERRE' 5 ;
  350. 'SINO' ;
  351. 'ERRE' 0 ;
  352. 'FINSI' ;
  353. *
  354. 'FIN' ;
  355.  
  356.  
  357.  
  358.  
  359.  
  360.  
  361.  
  362.  
  363.  
  364.  
  365.  
  366.  
  367.  
  368.  
  369.  

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