Télécharger darcy5.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : darcy5.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. *
  5. ******************** CAS TEST : darcy5.dgibi **********************
  6. *
  7. GRAPH = FAUX ;
  8. *GRAPH = VRAI ;
  9. 'SAUT' 'PAGE' ;
  10. *
  11. *-------------------------------------------------------------------
  12. * TEST DARCY5
  13. *
  14. * Validation de la résolution des équations de DARCY formulées en
  15. * (VITESSE - PRESSION) par comparaison avec
  16. * la résolution des équations de DARCY formulées en (VITESSE - CHARGE)
  17. *
  18. *
  19. * (Résolution par une méthode d'éléments finis mixtes hybrides)
  20. *
  21. * Le maillage est composés de triangles et de quadrangles et contient
  22. * des lignes de maillages courbes
  23. *
  24. * la massif est supposé isotrope, indéformable et saturé
  25. *
  26. * Pour permettre la comparaison, la densité de l'eau reste constante
  27. *
  28. *-------------------------------------------------------------------
  29. *
  30. OPTI ECHO 0 ;
  31. *
  32. 'SAUT' 'PAGE' ;
  33. *
  34. *- Options générales de calcul
  35. *
  36. TITRE
  37. 'EFMY DARCY FORMULATION (VITESSE - PRESSION) 2D : darcy5.dgibi' ;
  38.  
  39. OPTI DIME 2 ELEM QUA4 ;
  40.  
  41. *
  42. *- Données physiques
  43. *
  44. *- Densité de l'eau
  45. RHO0 = 1000. ;
  46. *
  47. *- Pression athmosphérique
  48. P0 = 1.013E5 ;
  49. *
  50. *- Composante du vecteur gravité
  51. Gx0 = 0. ;
  52. Gy0 = -9.81 ;
  53. *
  54. *- Perméabilité intrinsèque
  55. LA0 = 1.e-9 ;
  56. *
  57. *======================================================================
  58. *
  59. *- Construction du maillage
  60. *
  61. *======================================================================
  62. *
  63. *- Discrétisation spatiale
  64. *
  65. ENX = 15 ;
  66. ENY0 = 5 ;
  67. ENY1 = 5 ;
  68. ENY2 = 5 ;
  69. *
  70. *- Création des points
  71. *
  72. A0 = 0.d0 -2.d3 ; B0 = 5.d3 -2.d3 ;
  73. *
  74. A1 = 0.d0 -1.d3 ; B1 = 5.d3 -1.d3 ; C1 = 2.d3 -1.25d3 ;
  75. *
  76. A2 = 0.d0 0.0d3 ; B2 = 5.d3 0.0d3 ; C2 = 2.d3 0.25d3 ;
  77. *
  78. A3 = 0.d0 1.5d3 ; B3 = 5.d3 1.d3 ; C3 = 2.5d3 1.5d3 ;
  79. *
  80. *- Création des droites
  81. *
  82. AB0 = DROI ENX A0 B0 ;
  83. AB1 = A1 cer3 C1 ENX B1 ;
  84. AB2 = A2 cer3 C2 ENX B2 ;
  85. AB3 = A3 cer3 C3 ENX B3 ;
  86. *
  87. *- Creation des SURFACES
  88. *
  89. NIVEAU0 = AB0 regle ENY0 AB1 ;
  90. OPTI ELEM TRI3 ;
  91. NIVEAU1 = AB1 regle ENY1 AB2 ;
  92. OPTI ELEM QUA4 ;
  93. NIVEAU2 = AB2 regle ENY2 AB3 ;
  94. MASSIF0 = NIVEAU0 ET NIVEAU1 ET NIVEAU2 ;
  95. HAUT = INVE AB3 ;
  96. MASSIF0 = ORIENT MASSIF0 ;
  97. QMASSIF = CHANGE MASSIF0 QUAF ;
  98. QFHAUT = CHANGE HAUT QUAF ;
  99. ELIM 0.01 (QMASSIF ET QFHAUT ) ;
  100. *
  101. *- Domaines
  102. *
  103. *HYTOT = DOMA MASSIF0 0.01 ;
  104. *CHYB1 = DOMA HYTOT 'SURFACE' ;
  105. *CHYB2 = DOMA HYTOT 'NORMALE' ;
  106. *MCHYB = DOMA HYTOT 'ORIENTAT' ;
  107. *HYHAUT = DOMA HAUT INCL HYTOT 0.01 ;
  108. *
  109. *- Modélisation
  110. *
  111. *MODHYB = MODE HYTOT DARCY ISOTROPE HYQ4 HYT3 ;
  112. MODHYB = MODE QMASSIF DARCY ISOTROPE ;
  113. MODHAU = MODE QFHAUT DARCY ISOTROPE ;
  114. CHYB1 = DOMA MODHYB 'SURFACE' ;
  115. CHYB2 = DOMA MODHYB 'NORMALE' ;
  116. CEHAUT= DOMA MODHAU 'CENTRE' ;
  117. CETOT= DOMA MODHYB 'CENTRE' ;
  118. *
  119. *- perméabilité intrinsèque
  120. *
  121. LMME = MANU CHPO CETOT 1 'K' la0 'NATURE' 'DIFFU' ;
  122. *
  123. *======================================================================
  124. *
  125. *- Résolution en en (VITESSE ; PRESSION)
  126. *
  127. *======================================================================
  128. *
  129. LMMA = KCHA MODHYB 'CHAM' LMME ;
  130. MATI_P = MATERIAU MODHYB 'K' LMMA ;
  131. *
  132. **** Masse HYBride
  133. CND1A_P = MHYBR MODHYB MATI_P ;
  134. *
  135. **** Masse FORce
  136. M_P = MHYBR MODHYB 'MASSE' ;
  137. *
  138. **** MAtrice en Trace de pression TP
  139. HND1A_P = MATP MODHYB CND1A_P ;
  140. *
  141. *- Conditions aux limites
  142. *
  143. BHAUT_P = BLOQUE CEHAUT 'TH' ;
  144. *
  145. *- TP imposee
  146. *
  147. PIMP = MANU CHPO CEHAUT 1 'TH' P0 'NATURE' 'DISCRET' ;
  148. EHAUT_P = DEPI BHAUT_P PIMP ;
  149. *
  150. *- contribution des forces de volume au second membre
  151. *
  152. RHO = MANU CHPO CETOT 1 'SCAL' RHO0 'NATURE' 'DISCRET' ;
  153. RCH = MANU CHPO CETOT 2 'FX' gx0 'FY' gy0 'NATURE' 'DISCRET' ;
  154. RCH = RHO * RCH ;
  155. GRAV2TP = SQTP MODHYB HND1A_P M_P RCH ;
  156. *
  157. *- Assemblage matrice et second membre
  158. *
  159. CCC1_P = HND1A_P ET BHAUT_P ;
  160. * second membre du système matriciel en trace de charge
  161. FFF1_P = GRAV2TP et EHAUT_P ;
  162. *
  163. *- Resolution en trace de pression
  164. *
  165. CHTER1_P = RESOUDRE CCC1_P FFF1_P ;
  166. CHTER1_P = exco CHTER1_P 'TH' 'TH';
  167. *
  168. *- Calcul des pressions élémentaires
  169. *
  170. PCEN1 = HYBP MODHYB CND1A_P CHTER1_P M_P RCH ;
  171. *
  172. *- Calcul des débits
  173. *
  174. QFACE1_P = HDEBI MODHYB CND1A_P PCEN1 CHTER1_P M_P RCH;
  175. VCENT1_P = HVIT MODHYB QFACE1_P ;
  176. *VV = VECT VCENT1_P 0.5e9 'VX' 'VY' VERT ;
  177. *titre 'VITESSES AUX CENTRES DES ELEMENTS : ECHL=1e9:2' ;
  178. *TRAC VV CETOT ;
  179. *
  180. *=======================================================================
  181. *
  182. *- Résolution en (VITESSE ; CHARGE)
  183. *
  184. *=======================================================================
  185. *
  186. *
  187. *- conductivité hydraulique
  188. *
  189. LMME = LMME * (ABS Gy0) * RHO0 ;
  190. LMMA = KCHA MODHYB 'CHAM' LMME ;
  191. MATI_H = MATERIAU MODHYB 'K' LMMA ;
  192. *
  193. **** Masse HYBride
  194. CND1A_H = MHYBR MODHYB MATI_H ;
  195. *
  196. **** MAtrice en Trace de charge TH
  197. HND1A_H = MATP MODHYB CND1A_H ;
  198. *
  199. *- Conditions aux limites
  200. *
  201. BHAUT_H = BLOQUE CEHAUT 'TH' ;
  202. *
  203. *- TH imposee
  204. *
  205. PIMP = NOMC 'TH' (CEHAUT COOR 2) ;
  206. EHAUT_H = DEPI BHAUT_H PIMP ;
  207. *
  208. *- Assemblage matrice et second membre
  209. *
  210. CCC1_H = HND1A_H ET BHAUT_H ;
  211. * second membre du système matriciel en trace de charge
  212. FFF1_H = EHAUT_H ;
  213. *
  214. *- Resolution en trace de charge
  215. *
  216. CHTER1_H = RESOUDRE CCC1_H FFF1_H ;
  217. CHTER1_H = exco CHTER1_H 'TH' 'TH';
  218. *
  219. *- Calcul de la charge
  220. *
  221. TCEN1 = HYBP MODHYB CND1A_H CHTER1_H ;
  222. *
  223. *- Calcul de la vitesse
  224. *
  225. QFACE1_H = HDEBI MODHYB CND1A_H TCEN1 CHTER1_H ;
  226. VCENT1_H = HVIT MODHYB QFACE1_H ;
  227. *VV = VECT VCENT1_H 0.5e9 'VX' 'VY' VERT ;
  228. *titre 'CALCUL EN CHARGE : VITESSES AUX CENTRES DES ELTS : ECHL=1e9:2';
  229. *TRAC VV CETOT ;
  230. *
  231. *=======================================================================
  232. *
  233. *- Comparaison des deux résolutions
  234. *
  235. *=======================================================================
  236. *
  237. *
  238. *- comparaison des traces de charges
  239. *
  240. TP2 = (NOMC CHTER1_P SCAL) - P0 / (rho0 * gy0) ;
  241. TP2 = ((DOMA MODHYB 'FACE') COOR 2) - TP2 ;
  242. ER2 = (NOMC CHTER1_H SCAL)/TP2 - 1. ;
  243. maxtp1 = ER2 ABS MAXI ;
  244. *
  245. *
  246. *- comparaison des charges
  247. *
  248. PCEN2 = (NOMC PCEN1 SCAL) - P0 / (rho0 * gy0) ;
  249. PCEN2 = CETOT COOR 2 - PCEN2 ;
  250. ER2 = (NOMC TCEN1 SCAL)/PCEN2 - 1. ;
  251. maxp1 = ER2 ABS MAXI ;
  252.  
  253. 'SI' GRAPH ;
  254. TITRE 'VARIATIONS RELATIVES DES CHARGES DH/H' ;
  255. TRAC (KCHA MODHYB 'CHAM' er2) MODHYB (CONTOUR
  256. (DOMA MODHYB 'MAILLAGE')) ;
  257. 'FINSI' ;
  258. *
  259. *- comparaison des vitesse
  260. *
  261. QFACE1_P = NOMC SCAL QFACE1_P ;
  262. VFACE1_P = QFACE1_P * CHYB2 / CHYB1 ;
  263. QFACE1_H = NOMC SCAL QFACE1_H ;
  264. VFACE1_H = QFACE1_H * CHYB2 / CHYB1 ;
  265. maxvf1 = (kops QFACE1_P '-' QFACE1_H) / (maxi QFACE1_H) ABS MAXI ;
  266.  
  267. MOT1 = 'MOTS' 'VX' 'VY' ;
  268. VDVD = 'PSCA' VCENT1_H VCENT1_H MOT1 MOT1 ;
  269. VD1 = VCENT1_P - VCENT1_H ;
  270. VC1 = 'PSCA' VD1 VD1 MOT1 MOT1 ;
  271. SDC1 = 'ABS' ( VC1 / VDVD ) ;
  272. SDC1 = SDC1 '**' 0.5D0 ;
  273. maxvc1 = MAXI SDC1 ;
  274.  
  275. 'SI' GRAPH ;
  276. TITRE ' VARIATIONS RELATIVES DES VITESSES |DV|/|V|' ;
  277. TRAC ('KCHA' MODHYB 'CHAM' SDC1) MODHYB (CONTOUR
  278. (DOMA MODHYB 'MAILLAGE')) ;
  279. 'FINSI' ;
  280.  
  281. *
  282. 'SAUT' 'PAGE' ;
  283. 'SAUT' 2 'LIGNE' ;
  • 'MESS' ' VARIATIONS RELATIVES ' ;
  • 'SAUT' 1 'LIGNE' ;
  • 'MESS' ' TH H Vf Vc' ;
  • 'SAUT' 1 'LIGNE' ;
  • 'MESS' ' ' maxtp1 ' ' maxp1 ' ' maxvf1 ' ' maxvc1 ;
  • 'SAUT' 2 'LIGNE' ;
  • *
  • EPS0 = 1.E-7 ;
  • LOG1 = MAXTP1 > EPS0 ;
  • LOG2 = MAXP1 > EPS0 ;
  • LOG3 = MAXVF1 > EPS0 ;
  • LOG4 = MAXVC1 > EPS0 ;
  • L0 = LOG1 'OU' LOG2 'OU' LOG3 'OU' LOG4 ;
  • *
  • 'SI' ( L0 ) ;
  • 'ERRE' 5 ;
  • 'SINO' ;
  • 'ERRE' 0 ;
  • 'FINS' ;
  •  
  • FIN ;
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  • © Cast3M 2003 - Tous droits réservés.
    Mentions légales