Télécharger la plaquette Cast3M


fr en

Télécharger darcy9.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  darcy1.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. *****************  CAS TEST : darcy1.dgibi  *************************
  5. *
  6. GRAPH = faux   ;
  7. 'SAUT' 'PAGE' ;
  8. *
  9. *-------------------------------------------------------------------
  10. * TEST DARCY9
  11. *
  12. * C'est une reprise de DARCY1
  13. * avec utilisation de conditions aux limites d'inégalité
  14. *
  15. *-  1 premier calcul comme DARCY1 pour vérifier que les modifs conservent
  16. *     le résultat
  17. *-  1 deuxième calcul avec CL d'inégalité et terme source pour l'exemple.
  18. *
  19. * Les modifications sont repérées par ------>
  20. *
  21. *  . multiplication par -1 de la matrice issue de MATP,
  22. *                          du terme source issu de SMTP
  23. *                          des débits au second membre,
  24. *  . on ne touche pas aux multiplicateurs de Lagrange (BLOQ et DEPI)
  25. *
  26. *-------------------------------------------------------------------
  27. *
  28. 'SAUT' 'PAGE' ;
  29. *
  30. *- Options générales de calcul.
  31. *
  32. 'TITR' 'EFMH DARCY CL inégalité : darcy9.dgibi' ;
  33. * Dimension 2, éléments à générer quadrangulaires ou triangulaires :
  34. 'OPTI' 'DIME' 2 'ELEM' 'TRI3'       ;
  35. * Pas d'affichage à l'écran des lignes de commandes lues :
  36. 'OPTI' 'ECHO' 0           ;
  37. *                                                     ------------
  38. *                                                     = MAILLAGE =
  39. *                                                     ------------
  40. *
  41. *- Création des points supports des DROITES
  42. *
  43. * Dimensions :
  44. L  = 1.D0    ;
  45. HA = 1.D0    ;
  46.  
  47. * Coordonnées :
  48. XG = 0.D0    ;
  49. XD = XG + L  ;
  50. YG = 0.D0    ;
  51. YD = YG + HA ;
  52.  
  53. * Points
  54. A1 = XG YG   ;
  55. A3 = XD YG   ;
  56. D1 = XG YD   ;
  57. D3 = XD YD   ;
  58. *
  59. *- Création des DROITES frontières
  60. *
  61. * Discrétisation dans les deux directions :
  62. INUMX = 9      ;
  63. INUMY = 4      ;
  64. * droites maillées :
  65. DRBAS = A3 'DROI' INUMX A1 ;
  66. DRGAU = A1 'DROI' INUMY D1 ;
  67. DRHAU = D1 'DROI' INUMX D3 ;
  68. DRDRO = D3 'DROI' INUMY A3 ;
  69. *
  70. *- Création maillage GEOMETRIQUE
  71. *
  72. * on va faire des parallèlépipèdes réguliers grâce à DALLER
  73. * Les points définissant les quadrangles sont appellés 'SOMMET'
  74. CARRE = 'DALLER' DRBAS DRGAU DRHAU DRDRO 'PLAN' ;
  75. *
  76. *- Création maillage HYBRIDE y compris sous-objets (cond. limites)
  77. * ce sont les maillages avec les points sommets + les points faces +
  78. * les points centres.
  79.  
  80. * on crée celui du maillage total
  81. DOMHYB = 'CHANGER' CARRE 'QUAF' ;
  82. * et ceux des droites sur lesquelles on voudra mettre des conditions
  83. * aux limites (en fait les quatre côtés, ici) :
  84. DOMGAU = 'CHANGER' DRGAU 'QUAF' ;
  85. DOMDRO = 'CHANGER' DRDRO 'QUAF' ;
  86. DOMHAU = 'CHANGER' DRHAU 'QUAF' ;
  87. DOMBAS = 'CHANGER' DRBAS 'QUAF' ;
  88.  
  89. * On a créé à chaque fois les points faces et centres.
  90. * Comme CARRE, DRGAU, DRDRO, DRHAU et DRBAS avaient des segments en
  91. * commun, les opérations précédentes ont créé beaucoup de points en
  92. * double. Il faut les éliminer pour qu'on parle bien de la même chose.
  93.  
  94. * critère d'élimination (inférieur à la plus petite dimension d'arête
  95. ELI0  = L / INUMX / 10.D0   ;
  96. * élimination proprement dite :
  97. 'ELIMINATION' ( DOMHYB ET DOMGAU ET DOMDRO ET DOMHAU ET DOMBAS) ELI0 ;
  98. *                                                 ----------------
  99. *                                                 = MODELISATION =
  100. *                                                 ----------------
  101. * Définition du modèle physique attaché au domaine d'étude :
  102. * Il servira d'une part pour la résolution,
  103. * d'autre part pour extraire les points faces, centres, et autres
  104. * informations utiles
  105. MODHYB = 'MODELE' DOMHYB 'DARCY' 'ISOTROPE'  ;
  106.  
  107. * On voudra aussi disposer des points faces des frontières,
  108. * c'est-à-dire des points centres des lignes constituant ces
  109. * frontières. Pour cela, on attache à ces lignes une modélisation,
  110. * inutile en soi, mais qui dit implicitement (comme c'est un modèle
  111. * 'DARCY'), que l'on travaille en éléments finis mixte hybrides,
  112. * et donc qu'on aura besoin, à un moment ou à un autre, des points
  113. * centres, faces, et autres :
  114. MODGAU = 'MODELE' DOMGAU 'DARCY' 'ISOTROPE'  ;
  115. MODDRO = 'MODELE' DOMDRO 'DARCY' 'ISOTROPE'  ;
  116. MODHAU = 'MODELE' DOMHAU 'DARCY' 'ISOTROPE'  ;
  117. MODBAS = 'MODELE' DOMBAS 'DARCY' 'ISOTROPE'  ;
  118.  
  119. * On peut maintenant extraire via ces modèles, les points centre des
  120. * maillages 'QUAF' des frontières, qui sont en fait les points faces
  121. * en limite du domaine :
  122. CEDRO  = 'DOMA' MODDRO 'CENTRE' ;
  123. CEHAU  = 'DOMA' MODHAU 'CENTRE' ;
  124. CEBAS  = 'DOMA' MODBAS 'CENTRE' ;
  125. CEGAU  = 'DOMA' MODGAU 'CENTRE' ;
  126.  
  127. * Les autres petites informations, sur le domaine entier :
  128. * Maillage des points centres :
  129. HYCEN = 'DOMA' MODHYB 'CENTRE' ;
  130. * Maillage des points faces :
  131. HYFAC = 'DOMA' MODHYB 'FACE' ;
  132. * Longueur des arêtes (on est en 2D) :
  133. HYSUR = 'DOMA' MODHYB 'SURFACE'  ;
  134. * Surface des mailles :
  135. HYVOL = 'DOMA' MODHYB 'VOLUME'  ;
  136. * Vecteurs normaux aux faces, champ-points à autant de composantes que
  137. * la dimension de l'espace : 'UX' , 'UY' (, 'UZ').
  138. HYNOR = 'DOMA' MODHYB 'NORMALE'  ;
  139. *
  140. *                             ------------------------------------
  141. *                             = Tenseur de perméabilité isotrope =
  142. *                             ------------------------------------
  143. * La grandeur physique correspondant au laplacien. C'est ici la
  144. * perméabilité hydraulique.
  145. VK    = 0.75D0   ;
  146. * On la stoque dans un objet de type matériau
  147. MATI3 = 'MATERIAU' MODHYB 'K' VK ;
  148.  
  149. * ==============
  150. * Premier Calcul
  151. * ==============
  152.  
  153. * Cas test 1D classique pour valider les modifications sur solution
  154. * analytique
  155.  
  156. *                                          = Solution analytique =
  157. *                                          -----------------------
  158. XX YY = 'COOR' (DOMA MODHYB 'FACE')    ;
  159. HG    = 100.D0     ;
  160. HD    =  20.D0     ;
  161. DH    = HD - HG    ;
  162. PANAF = (XX - L) * DH / L + HD ;
  163.  
  164. XXC YYC = 'COOR' (DOMA MODHYB 'CENTRE' )  ;
  165. PANAC   = (XXC - L) * DH / L + HD ;
  166.  
  167. VXANA = ( -1.D0 ) * VK * DH / L ;
  168. VYANA = 0.D0        ;
  169.  
  170. VANAF = 'MANU' 'CHPO' (DOMA MODHYB 'FACE')   2 'VX' VXANA 'VY' VYANA ;
  171. VANAC = 'MANU' 'CHPO' (DOMA MODHYB 'CENTRE') 2 'VX' VXANA 'VY' VYANA ;
  172. *
  173. *                                             = Modèle numérique =
  174. *                                             --------------------
  175. CND1A = 'MHYB' MODHYB  MATI3 ;
  176. *------> on multiplie par -1. pour qu'elle soit définie positive
  177. HND1A = -1.*('MATP' MODHYB  CND1A) ;
  178.  
  179. * CLs de Dirichlet
  180. * Contribution au membre de gauche pour chaque sous-partie de la frontière
  181. BBGAU = 'BLOQ' CEGAU 'TH' ;
  182. BBDRO = 'BLOQ' CEDRO 'TH' ;
  183. BBHAU = 'BLOQ' CEHAU 'TH' ;
  184. BBBAS = 'BLOQ' CEBAS 'TH' ;
  185.  
  186. * Contribution au membre de droite pour chaque sous-partie de la frontière
  187. TTIMP = 'REDU' PANAF CEGAU ;
  188. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  189. EEGAU = 'DEPI' BBGAU TTIM2          ;
  190. *
  191. TTIMP = 'REDU' PANAF CEDRO ;
  192. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  193. EEDRO = 'DEPI' BBDRO TTIM2          ;
  194. *
  195. TTIMP = 'REDU' PANAF CEBAS ;
  196. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  197. EEBAS = 'DEPI' BBBAS TTIM2          ;
  198. *
  199. TTIMP = 'REDU' PANAF CEHAU ;
  200. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  201. EEHAU = 'DEPI' BBHAU TTIM2          ;
  202.  
  203. * CL de Neuman
  204. FG    = VK * DH / L     ;
  205.  
  206. FGAU  = 'MANU' 'CHPO' CEGAU 1 'FLUX' FG   ;
  207. LGAU  = 'REDU' HYSUR  CEGAU   ;
  208. FLGAU = FGAU * LGAU           ;
  209. FLGAU = 'CHAN' FLGAU 'ATTRIBUT' 'NATURE' 'DISCRET' ;
  210.  
  211. FD    = ( -1.D0 ) * FG  ;
  212. FDRO  = 'MANU' 'CHPO' CEDRO 1 'FLUX' FD   ;
  213. LDRO  = 'REDU' HYSUR  CEDRO   ;
  214. FLDRO = FDRO * LDRO           ;
  215. FLDRO = 'CHAN' FLDRO 'ATTRIBUT' 'NATURE' 'DISCRET' ;
  216.  
  217. * Assemblage
  218. *------> les flux sont multipliés par -1. pour correspondre à HND1A
  219. CCC3  = HND1A 'ET' BBHAU 'ET' BBBAS  ;
  220. FFF3  =            EEHAU 'ET' EEBAS 'ET' (-1.* (FLDRO 'ET' FLGAU)) ;
  221.  
  222. CHTER3 = 'RESO' CCC3 FFF3 ;
  223.  
  224. PCEN3  = 'HYBP' MODHYB CND1A CHTER3 ;
  225. QFACE3 = 'HDEB' MODHYB CND1A PCEN3 CHTER3 ;
  226. VCENT3 = 'HVIT' MODHYB QFACE3 ;
  227. VFACE3 = ('EXCO' QFACE3 'FLUX' 'SCAL') * HYNOR / HYSUR;
  228.  
  229. * Calcul ERREUR
  230. ERRTP3 = 'EXCO' CHTER3 'TH' 'SCAL' ;
  231. ERRTP3 = ERRTP3 - PANAF / PANAF    ;
  232. ERRTP3 = 'ABS' ERRTP3  ;
  233. MAXTP3 = 'MAXI' ERRTP3 ;
  234.  
  235. ERRP3  = 'EXCO' PCEN3 'H' 'SCAL'   ;
  236. ERRP3  = ERRP3 - PANAC / PANAC     ;
  237. ERRP3  = 'ABS' ERRP3   ;
  238. MAXP3  = 'MAXI' ERRP3  ;
  239.  
  240. MOT1   = 'MOTS' 'VX' 'VY' ;
  241. VDVD   = 'PSCA' VANAC VANAC MOT1 MOT1 ;
  242. VD3    = VANAC - VCENT3           ;
  243. VC3    = 'PSCA' VD3 VD3 MOT1 MOT1 ;
  244. SDC3   = 'ABS' ( VC3 / VDVD )     ;
  245. SDC3   = SDC3 '**' 0.5D0          ;
  246. MAXV3  = 'MAXI' SDC3   ;
  247.  
  248. 'SI' (GRAPH) ;
  249.   'TITRE' 'Charge pour pb uniforme 1D' ;
  250.   'TRACER' ('KCHA' PCEN3 modhyb 'CHAM') modhyb ;
  251. 'FINSI' ;
  252.  
  253. * ===============
  254. * Deuxième Calcul
  255. * ===============
  256.  
  257. * Cas avec CL d'inégalité.
  258.  
  259. * Nouvelles CL
  260. * ------------
  261. * on avait Hgauche=HG=100, on met Hgauche >= HG
  262. BBGAU = 'BLOQ' CEGAU 'MINIMUM' 'TH' ;
  263. TTIMP = 'REDU' PANAF CEGAU ;
  264. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  265. EEGAU = 'DEPI' BBGAU TTIM2          ;
  266. *
  267. * on avait Hdroite=HD=20 , on met Hdroite <= HD
  268. BBDRO = 'BLOQ' CEDRO 'MAXIMUM' 'TH' ;
  269. TTIMP = 'REDU' PANAF CEDRO ;
  270. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  271. EEDRO = 'DEPI' BBDRO TTIM2          ;
  272.  
  273. * Résolution et post-traitement
  274. CCC4   = HND1A 'ET' BBHAU 'ET' BBBAS 'ET' BBGAU 'ET' BBDRO ;
  275. FFF4   =            EEHAU 'ET' EEBAS 'ET' EEGAU 'ET' EEDRO ;
  276. CHTER4 = 'RESO' CCC4 FFF4 ;
  277.  
  278. PCEN4  = 'HYBP' MODHYB CND1A CHTER4 ;
  279. QFACE4 = 'HDEB' MODHYB CND1A PCEN4 CHTER4 ;
  280. VCENT4 = 'HVIT' MODHYB QFACE4 ;
  281. VFACE4 = ('EXCO' QFACE4 'FLUX' 'SCAL') * HYNOR / HYSUR;
  282.  
  283. 'SI' (GRAPH) ;
  284.   'TITRE' 'Charge pour pb uniforme 1D avec CL unilatérales' ;
  285.   'TRACER' ('KCHA' PCEN4 modhyb 'CHAM') modhyb ;
  286. 'FINSI' ;
  287.  
  288.  
  289. * Calcul ERREUR
  290. ERRTP4 = 'EXCO' CHTER4 'TH' 'SCAL' ;
  291. ERRTP4 = ERRTP4 - PANAF / PANAF    ;
  292. ERRTP4 = 'ABS' ERRTP4  ;
  293. MAXTP4 = 'MAXI' ERRTP4 ;
  294.  
  295. ERRP4  = 'EXCO' PCEN4 'H' 'SCAL'   ;
  296. ERRP4  = ERRP4 - PANAC / PANAC     ;
  297. ERRP4  = 'ABS' ERRP4   ;
  298. MAXP4  = 'MAXI' ERRP4  ;
  299.  
  300. MOT1   = 'MOTS' 'VX' 'VY' ;
  301. VDVD   = 'PSCA' VANAC VANAC MOT1 MOT1 ;
  302. VD4    = VANAC - VCENT4           ;
  303. VC4    = 'PSCA' VD4 VD4 MOT1 MOT1 ;
  304. SDC4   = 'ABS' ( VC4 / VDVD )     ;
  305. SDC4   = SDC4 '**' 0.5D0          ;
  306. MAXV4  = 'MAXI' SDC4   ;
  307.  
  308. * Troisème Calcul
  309. * ===============
  310.  
  311. * Cas avec CL d'inégalité.
  312. * et un terme source pour illustrer son emploi.
  313.  
  314. * on l'impose au point centre proche du barycentre
  315. pt     = 'BARYCENTRE' carre ;
  316. pts    = 'POINT' hycen 'PROC' pt ;
  317. cs1    = 'MANU' 'CHPO' ('MANU' pts 'POI1') 1 'SOUR' 100. ;
  318. tersc1 = 'KCHT' modhyb 'SCAL' 'CENTRE' 'COMP' 'SOUR' 0. cs1 ;
  319. *------> les termes sources sont multipliés par -1. pour correspondre à HND1A
  320. SMTR   = -1.*('SMTP' modhyb CND1A tersc1) ;
  321.  
  322. * Résolution et post-traitement
  323. CCC5   = HND1A 'ET' BBHAU 'ET' BBBAS 'ET' BBGAU 'ET' BBDRO ;
  324. FFF5   =            EEHAU 'ET' EEBAS 'ET' EEGAU 'ET' EEDRO 'ET' SMTR ;
  325. CHTER5 = 'RESO' CCC5 FFF5 ;
  326.  
  327. PCEN5  = 'HYBP' MODHYB CND1A CHTER5 ;
  328. QFACE5 = 'HDEB' MODHYB CND1A PCEN5 CHTER5 ;
  329. VCENT5 = 'HVIT' MODHYB QFACE5 ;
  330. VFACE5 = ('EXCO' QFACE5 'FLUX' 'SCAL') * HYNOR / HYSUR;
  331.  
  332. 'SI' (GRAPH) ;
  333.   'TITRE' 'Charge pour pb uniforme 1D avec CL unilatérales' ;
  334.   'TRACER' ('KCHA' PCEN5 modhyb 'CHAM') modhyb ;
  335. 'FINSI' ;
  336.  
  337. * On vérifie que la charge a bien varié dans le bon sens (positif)
  338. vh4    = 'EXTR' pcen4 'H' pts ;
  339. vh5    = 'EXTR' pcen5 'H' pts ;
  340.  
  341. LOG3   = 'NON' (vh5 > vh4) ;
  342. *
  343. *                                                 -------------------
  344. *                                                 = Gestion ERREURS =
  345. *                                                 -------------------
  346. *
  347. 'SAUT' 'PAGE'                                                    ;
  348. 'SAUT' 2 'LIGNE'                                                 ;
  349. 'MESS' '                           ERREURS RELATIVES           ' ;
  350. 'SAUT' 1 'LIGNE'                                                 ;
  351. 'MESS' ' cas test        TH                  H                V' ;
  352. 'SAUT' 1 'LIGNE'                                                 ;
  353. 'MESS' ' numero3   ' maxtp3  '   ' maxp3 '   ' maxv3             ;
  354. 'SAUT' 1 'LIGNE'                                                 ;
  355. 'MESS' ' numero4   ' maxtp4  '   ' maxp4 '   ' maxv4             ;
  356. 'SAUT' 2 'LIGNE'                                                 ;
  357. *
  358. EPS0 = 1.E-13 ;
  359. LOG1 = MAXTP3 > EPS0 ; LOG2 = MAXTP4 > EPS0  ;
  360. LOG4 = MAXP3  > EPS0 ; LOG5 = MAXP4  > EPS0  ;
  361. LOG7 = MAXV3  > EPS0 ; LOG8 = MAXV4  > EPS0  ;
  362. LTP0 = LOG1 'OU' LOG2 ;
  363. LP0  = LOG4 'OU' LOG5 ;
  364. LV0  = LOG7 'OU' LOG8 ;
  365. L0   = LTP0 'OU' LP0  'OU' LV0 'OU' LOG3 ;
  366.  
  367. * Compte-rendu d'erreur
  368. 'SI' L0        ;
  369.    'ERREUR' 5  ;
  370. 'SINON'        ;
  371.    'ERREUR' 0  ;
  372. 'FINSI'        ;
  373.  
  374. 'FIN' ;
  375.  
  376.  
  377.  
  378.  
  379.  
  380.  
  381.  
  382.  
  383.  
  384.  
  385.  
  386.  
  387.  
  388.  
  389.  

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