Télécharger la plaquette Cast3M


fr en

Télécharger darcy2.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  darcy2.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. *
  5. ********************  CAS TEST : darcy2.dgibi  **********************
  6. *
  7.  GRAPH = 'N'   ;
  8. 'SAUT' 'PAGE' ;
  9. *
  10. *-------------------------------------------------------------------
  11. * TEST DARCY2
  12. *                  CALCUL DARCY ORTHOTROPE
  13. *    Résolution par une méthode d'éléments finis mixtes hybrides.
  14. *
  15. * Ce test permet de vérifier le bon fonctionnement des opérateurs
  16. * utilisés afin de résoudre les équations de DARCY par une méthode
  17. * d'éléments finis mixtes hybrides dans CASTEM2000.
  18. *
  19. * On effectue un calcul sur un domaine carré, maillé 'aléatoirement'
  20. * par des quadrangles et des triangles. Les conditions aux limites
  21. * sont du type charge imposée.
  22. *
  23. * La solution analytique en charge étant un polynome de degré un,
  24. * la vitesse est constante.
  25. *       H(x,y) = -45x -80y + 22.5
  26. *       V(x,y) = ( -45*KX ; -80*KY)
  27. *
  28. * On s'attend à peu d'erreur bien que le maillage soit 'aléatoire'.
  29. *
  30. *-------------------------------------------------------------------
  31. *
  32. 'SAUT' 'PAGE' ;
  33. *
  34. *- Options générales de calcul.
  35. *
  36. 'TITR' 'EFMY DARCY ORTHOTROPE : darcy2.dgibi' ;
  37. 'OPTI' 'DIME' 2 'ELEM' 'QUA4'                 ;
  38.  'OPTI' 'ECHO' 0                               ;
  39. *
  40. *                                                     ------------
  41. *                                                     = MAILLAGE =
  42. *                                                     ------------
  43. *- Création des points supports des DROITES
  44. *
  45. L  =  0.1D0         ;
  46. ML = -0.1D0         ;
  47. A0 = ML     ML      ;
  48. B0 = L      ML      ;
  49. A1 = 0.D0   ML      ;
  50. B1 = L      0.D0    ;
  51. C0 = L      L       ;
  52. D0 = ML     L       ;
  53. C1 = 0.05D0 L       ;
  54. D1 = ML     -0.05D0 ;
  55. *
  56. *- Création et découpage des DROITES frontieres
  57. *
  58. INUMX = 8                  ;
  59. INUMY = 8                  ;
  60. INUMP = 2                  ;
  61. DRBA1 = A0 'DROI' INUMX A1 ;
  62. DRBA2 = A1 'DROI' INUMP B0 ;
  63. DRBAS = DRBA1 'ET' DRBA2   ;
  64. DRDR1 = B0 'DROI' INUMP B1 ;
  65. DRDR2 = B1 'DROI' INUMY C0 ;
  66. DRDRO = DRDR1 'ET' DRDR2   ;
  67. DRHA1 = C0 'DROI' INUMP C1 ;
  68. DRHA2 = C1 'DROI' INUMX D0 ;
  69. DRHAU = DRHA1 'ET' DRHA2   ;
  70. DRGA1 = D0 'DROI' INUMP D1 ;
  71. DRGA2 = D1 'DROI' INUMY A0 ;
  72. DRGAU = DRGA1 'ET' DRGA2   ;
  73. *
  74. *- Creation maillage GEOMETRIQUE
  75. *
  76. BDHG  = DRBAS 'ET' DRDRO 'ET' DRHAU 'ET' DRGAU ;
  77. PTOT1 = 'SURF' BDHG                            ;
  78. PTOT2 = 'ORIE' PTOT1                           ;
  79. QFTOT = CHANGE PTOT2 QUAF ;
  80. *
  81. *- Creation maillage HYBRIDE y compris sous-objets (cond. limites)
  82. *
  83. ELI0  = L / INUMX / 100.D0             ;
  84. *HYTOT = 'DOMA' PTOT2 ELI0              ;
  85. *CHYB1 = 'DOMA' HYTOT 'SURFACE'         ;
  86. *CHYB2 = 'DOMA' HYTOT 'NORMALE'         ;
  87. *MCHYB = 'DOMA' HYTOT 'ORIENTAT'        ;
  88. *HYGAU = 'DOMA' DRGAU 'INCL' HYTOT ELI0 ;
  89. *HYDRO = 'DOMA' DRDRO 'INCL' HYTOT ELI0 ;
  90. *HYHAU = 'DOMA' DRHAU 'INCL' HYTOT ELI0 ;
  91. *HYBAS = 'DOMA' DRBAS 'INCL' HYTOT ELI0 ;
  92.  QFGAU = CHANGE DRGAU QUAF ;
  93.  QFDRO = CHANGE DRDRO QUAF ;
  94.  QFHAU = CHANGE DRHAU QUAF ;
  95.  QFBAS = CHANGE DRBAS QUAF ;
  96.  ELIM ELI0 ( QFTOT ET QFGAU ET QFDRO ET QFHAU ET QFBAS );
  97. *
  98. *                                                 ----------------
  99. *                                                 = MODELISATION =
  100. *                                                 ----------------
  101. *
  102. *MODHYB = MODE HYTOT 'DARCY' 'ORTHOTROPE' 'HYQ4' 'HYT3' ;
  103. MODHYB = MODE QFTOT 'DARCY' 'ORTHOTROPE'  ;
  104. *option donn 5 ;
  105. MODGAU = MODE QFGAU 'DARCY' 'ORTHOTROPE'  ;
  106. MODDRO = MODE QFDRO 'DARCY' 'ORTHOTROPE'  ;
  107. MODHAU = MODE QFHAU 'DARCY' 'ORTHOTROPE'  ;
  108. MODBAS = MODE QFBAS 'DARCY' 'ORTHOTROPE'  ;
  109. CEGAU= DOMA MODGAU 'CENTRE' ;
  110. CEDRO= DOMA MODDRO 'CENTRE' ;
  111. CEHAU= DOMA MODHAU 'CENTRE' ;
  112. CEBAS= DOMA MODBAS 'CENTRE' ;
  113. CETOT= DOMA MODHYB 'CENTRE' ;
  114. CHYB1 = 'DOMA' MODHYB 'SURFACE'         ;
  115. CHYB2 = 'DOMA' MODHYB 'NORMALE'         ;
  116. *
  117. *                           --------------------------------------
  118. *                           = Tenseur de perméabilité orthotrope =
  119. *                           --------------------------------------
  120. *
  121. VKX = 1.D0      ;
  122. VKY = 0.75D0    ;
  123. BZ  = 1.D0 0.D0 ;
  124. *
  125. MATI3 = 'MATE' MODHYB 'DIRECTION' BZ 'PARALLELE' 'K1' VKX 'K2' VKY ;
  126. *option donn 5 ;
  127. *
  128. *                                          -----------------------
  129. *                                          = Solution analytique =
  130. *                                          -----------------------
  131. *
  132. *                                                 Charge aux faces
  133. AA    =-45.D0                ;
  134. BB    =-80.D0                ;
  135. CC    = 22.5D0               ;
  136. AAA   = ( -1.D0 ) * VKX * AA ;
  137. BBB   = ( -1.D0 ) * VKY * BB ;
  138. *XX YY = 'COOR' HYTOT.'FACE'  ;
  139. XX YY = 'COOR' ( DOMA MODHYB 'FACE')  ;
  140. TANF1 = AA * XX              ;
  141. TANF2 = BB * YY              ;
  142. PANAF = TANF1 + TANF2 + CC   ;
  143. *
  144. *                                 Charge aux centres des éléments
  145. * La charge étant linéaire, la moyenne est au centre de l'élément
  146. *
  147. XXC YYC = 'COOR' ( DOMA MODHYB 'CENTRE') ;
  148. TANC1   = AA * XXC              ;
  149. TANC2   = BB * YYC              ;
  150. PANAC   = TANC1 + TANC2 + CC    ;
  151. *
  152. *                                Vitesse aux faces et aux centres
  153. *
  154. VANAF = 'MANU' 'CHPO' ( DOMA MODHYB 'FACE')   2 'VX' AAA 'VY' BBB ;
  155. VANAC = 'MANU' 'CHPO' (DOMA MODHYB 'CENTRE') 2 'VX' AAA 'VY' BBB ;
  156. *                                        -------------------------
  157. *                                        = matrice Masse HYBride =
  158. *                                        -------------------------
  159. CND1A = 'MHYB' MODHYB  MATI3 ;
  160. *option donn 5 ;
  161. *                                        -------------------------
  162. *                                        = MAtrice globale en TH =
  163. *                                        -------------------------
  164. HND1A = 'MATP' MODHYB  CND1A ;
  165. *                                                 ----------------
  166. *                                                 =  TH imposée  =
  167. *                                                 ----------------
  168. *BBGAU = 'BLOQ' HYGAU.'CENTRE' 'TH' ;
  169. BBGAU = 'BLOQ' CEGAU 'TH' ;
  170. BBDRO = 'BLOQ' CEDRO 'TH' ;
  171. BBHAU = 'BLOQ' CEHAU 'TH' ;
  172. BBBAS = 'BLOQ' CEBAS 'TH' ;
  173. *
  174. *TTIMP = 'REDU' PANAF HYGAU.'CENTRE' ;
  175. TTIMP = 'REDU' PANAF CEGAU ;
  176. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  177. EEGAU = 'DEPI' BBGAU TTIM2          ;
  178. *
  179. TTIMP = 'REDU' PANAF CEDRO ;
  180. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  181. EEDRO = 'DEPI' BBDRO TTIM2          ;
  182. *
  183. TTIMP = 'REDU' PANAF CEBAS ;
  184. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  185. EEBAS = 'DEPI' BBBAS TTIM2          ;
  186. *
  187. TTIMP = 'REDU' PANAF CEHAU ;
  188. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  189. EEHAU = 'DEPI' BBHAU TTIM2          ;
  190. *                                                ----------------
  191. *                                                =  Assemblage  =
  192. *                                                ----------------
  193. CCC1 = HND1A 'ET' BBHAU 'ET' BBBAS 'ET' BBDRO 'ET' BBGAU ;
  194. FFF1 =            EEHAU 'ET' EEBAS 'ET' EEDRO 'ET' EEGAU ;
  195. *                                               -----------------
  196. *                                               = Résolution Th =
  197. *                                               -----------------
  198. CHTER1 = 'RESO' CCC1 FFF1 ;
  199. *option donn 5 ;
  200. *                                                    ----------------
  201. *                                                    = Résolution H =
  202. *                                                    ----------------
  203. PCEN1 = 'HYBP' MODHYB  CND1A CHTER1 ;
  204. *                                                ----------------
  205. *                                                = Résolution V =
  206. *                                                ----------------
  207. QFACE1 = 'HDEB' MODHYB  CND1A PCEN1 CHTER1  ;
  208. VCENT1 = 'HVIT' MODHYB  QFACE1        ;
  209. QFACE1 = 'EXCO' QFACE1 'FLUX' 'SCAL'             ;
  210. VFACE1 = QFACE1 * CHYB2 / CHYB1                  ;
  211. *                                                   -----------------
  212. *                                                   = Calcul ERREUR =
  213. *                                                   -----------------
  214. ERRTP1 = 'EXCO' CHTER1 'TH' 'SCAL' ;
  215.  
  216. ERRTP1 = ERRTP1 - PANAF / PANAF    ;
  217. ERRTP1 = 'ABS' ERRTP1              ;
  218. *
  219. ERRP1 = 'EXCO' PCEN1 'H' 'SCAL' ;
  220. *option donn 5 ;
  221. ERRP1 = ERRP1 - PANAC / PANAC   ;
  222. ERRP1 = 'ABS' ERRP1             ;
  223. *
  224. MOT1  = 'MOTS' 'VX' 'VY'             ;
  225. MOT2  = 'MOTS' 'UX' 'UY'             ;
  226. VAVN  = 'PSCA' VANAF CHYB2 MOT1 MOT2 ;
  227. VNN   = VAVN * CHYB2                 ;
  228. DVDV  = 'PSCA' VNN VNN MOT2 MOT2     ;
  229. DVN1  = VNN - VFACE1                 ;
  230. DVN1N = 'PSCA' DVN1 DVN1 MOT2 MOT2   ;
  231. SDF1  = 'ABS' ( DVN1N / DVDV )       ;
  232. SDF1  = SDF1 '**' 0.5                ;
  233. *
  234. VDVD  = 'PSCA' VANAC VANAC MOT1 MOT1 ;
  235. VD1   = VANAC - VCENT1               ;
  236. VC1   = 'PSCA' VD1 VD1 MOT1 MOT1     ;
  237. SDC1  = 'ABS' ( VC1 / VDVD )         ;
  238. SDC1  = SDC1 '**' 0.5                ;
  239. *                                                 -------------------
  240. *                                                 = TRACE RESULTATS =
  241. *                                                 -------------------
  242. *  Dans chaque cas on trace :
  243. *       L'erreur relative sur la charge H
  244. *       L'erreur relative sur la vitesse V
  245. *
  246. 'SI' ('NEG' GRAPH 'N') ;
  247. *    ERRP1 = 'KCHA' HYTOT 'CHAM' ERRP1         ;
  248.      ERRP1 = 'KCHA' MODHYB 'CHAM' ERRP1         ;
  249. *    SDS1  = 'KCHA' HYTOT 'CHAM' SDC1          ;
  250.      SDS1  = 'KCHA' MODHYB 'CHAM' SDC1          ;
  251.      'TITR' 'darcy2 : erreur relative charge'  ;
  252.      'TRAC' MODHYB ERRP1                       ;
  253.      'TITR' 'darcy2 : erreur relative vitesse' ;
  254.      'TRAC' MODHYB SDS1                        ;
  255. 'FINSI' ;
  256. *                                                 -------------------
  257. *                                                 = Gestion ERREURS =
  258. *                                                 -------------------
  259. MAXTP1 = 'MAXI' ERRTP1 ;
  260. MAXP1  = 'MAXI' ERRP1  ;
  261. MAXVF1 = 'MAXI' SDF1   ;
  262. MAXVC1 = 'MAXI' SDC1   ;
  263. *
  264. 'SAUT' 'PAGE'                                                        ;
  265. 'SAUT' 2 'LIGNE'                                                     ;
  266. 'MESS' '                        ERREURS RELATIVES                '   ;
  267. 'SAUT' 1 'LIGNE'                                                     ;
  268. 'MESS' '      TH                H                Vf              Vc' ;
  269. 'SAUT' 1 'LIGNE'                                                     ;
  270. 'MESS' ' ' maxtp1  ' ' maxp1 ' ' maxvf1 ' ' maxvc1                   ;
  271. 'SAUT' 2 'LIGNE'                                                     ;
  272. *
  273. EPS0 = 1.E-12 ; EPS1 = 1.E-10             ;
  274. LOG1 = MAXTP1 > EPS0                      ;
  275. LOG2 = MAXP1  > EPS0                      ;
  276. LOG3 = MAXVF1 > EPS1                      ;
  277. LOG4 = MAXVC1 > EPS0                      ;
  278. L0   = LOG1 'OU' LOG2 'OU' LOG3 'OU' LOG4 ;
  279. *
  280. 'SI' ( L0 ) ;
  281.    'ERRE' 5 ;
  282. 'SINO'      ;
  283.    'ERRE' 0 ;
  284. 'FINS'      ;
  285. *
  286. 'FIN' ;
  287.  
  288.  
  289.  
  290.  
  291.  
  292.  
  293.  
  294.  
  295.  
  296.  
  297.  

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