Télécharger hdebit.eso

Retour à la liste

Numérotation des lignes :

  1. C HDEBIT SOURCE BP208322 17/04/18 21:15:07 9396
  2. SUBROUTINE HDEBIT
  3. C-----------------------------------------------------------------------
  4. C Calcul du débit à travers chaque face dans le cas d'une
  5. C formulation DARCY en éléments finis mixte hybride.
  6. C-----------------------------------------------------------------------
  7. C
  8. C Une face appartient en général à deux éléments. On ne définit qu'une
  9. C normale à la face. Le ddl aux faces est VN * S * signe(normale)
  10. C avec signe=1 si normale sortante et -1 sinon, S=Surface de la face et
  11. C VN=produit scalaire entre la vitesse à la face et la normale sortante.
  12. C
  13. C On obtient le débit en effectuant le produit du champoint de trace de
  14. C charge TH par une matrice - exprimée en fonction de la matrice de
  15. C rigidité de sous type DARCY et de la matrice DIV (matrice uniligne
  16. C composé de 1) - d'une part, par le MCHAML des orientations de normale
  17. C d'autre part.
  18. C
  19. C Comme nous utilisons des élémennts finis non conforme, il ne faut pas
  20. C assembler les quantitées calculées. Comme le flux est continu au signe
  21. C près, on ne stocke que le flux dans le sens de la normale calculée.
  22. C
  23. C -1 t (n) -1 (n)
  24. C Q = M1 * DIV * H - M1 * TH
  25. C
  26. C
  27. C avec M1 : Matrice massse hybride
  28. C DIV : Matrice uniligne representant la divergence
  29. C
  30. C
  31. C---------------------------
  32. C Phrase d'appel (GIBIANE) :
  33. C---------------------------
  34. C
  35. C CHP1 = HDEB MMODEL MRIGI1 CHP2 CHP3 (MRIGI2 CHPO4) (CHP5);
  36. C
  37. C------------------------
  38. C Opérandes et résultat :
  39. C------------------------
  40. C
  41. C CHP1 : CHAMPOINT contenant le débit porté par la normale à la face
  42. C MMODEL : Objet MMODEL spécifiant la formulation
  43. C MRIGI1 : RIGIDITE masses hybrides élémentaires de sous-type DARCY.
  44. C MRIGI2 : Matrices masses hybrides elementaires de sous-type MASSE.
  45. C CHP2 : CHPOINT centre contenant la charge H (ou concentration)
  46. C CHP3 : CHPOINT contenant la trace de charge TH
  47. C (ou les traces de concentration)
  48. C CHPO4 : CHPO centre de composante FX,FY(,FZ),
  49. C densite de force moyenne par élément
  50. C CHPO5 : CHPOIN face contenant le flux de la vitesse convective
  51. C
  52. C
  53. C-----------------------------------------------------------------------
  54. C
  55. C Langage : ESOPE + FORTRAN77
  56. C
  57. C Auteurs : 09/93 F.DABBENE - Cas permanent
  58. C 09/94 X.NOUVELLON - Extension au cas transitoire
  59. C 02/96 L.V.BENET - Prise en compte de forces de volume
  60. C 06/98 F.AURIOL - formulation prenant en compte les valeurs
  61. C des charges ou des concentrations au centre
  62. C des elements. Prise en compte de la vitesse
  63. C convective
  64. C
  65. C-----------------------------------------------------------------------
  66. IMPLICIT INTEGER(I-N)
  67. IMPLICIT REAL*8 (A-H,O-Z)
  68. *
  69. -INC CCOPTIO
  70. -INC SMCHAML
  71. -INC SMCHPOI
  72. -INC SMELEME
  73. -INC SMMODEL
  74. -INC SMRIGID
  75. -INC SMTABLE
  76. *
  77. SEGMENT IPMAHY
  78. INTEGER MAHYBR(NSOUS)
  79. ENDSEGMENT
  80. *
  81. LOGICAL LOGRE,LOGIN
  82. CHARACTER*4 NOMTOT(1)
  83. CHARACTER*4 NOMTO3(3)
  84. CHARACTER*4 NOMTH
  85. CHARACTER*6 CHAR6
  86. CHARACTER*8 TAPIND,TYPOBJ,CHARIN,CHARRE,LETYPE
  87. *
  88. * Initialisations
  89. *
  90. IVALIN = 0
  91. XVALIN = 0.D0
  92. LOGIN = .TRUE.
  93. IOBIN = 0
  94. TAPIND = 'MOT '
  95. *
  96. * Lecture du MMODEL et test de la formulation
  97. *
  98. CALL LIROBJ('MMODEL',IPMODE,1,IRET)
  99. IF (IERR.NE.0) RETURN
  100. MMODEL = IPMODE
  101. SEGACT MMODEL
  102. N1=KMODEL(/1)
  103. DO 7 I=1,N1
  104. IMODEL=KMODEL(I)
  105. SEGACT IMODEL
  106. IF(FORMOD(1).NE.'DARCY')THEN
  107. MOTERR(1:16) = 'DARCY '
  108. CALL ERREUR(719)
  109. RETURN
  110. ENDIF
  111. SEGDES IMODEL
  112. 7 CONTINUE
  113. SEGDES MMODEL
  114. *
  115. * Lecture de la TABLE domaine
  116. *
  117. IPTABL = 0
  118. CALL LEKMOD(MMODEL,IPTABL,IRET)
  119. IF (IERR.NE.0) RETURN
  120. CHARIN = 'MAILLAGE'
  121. TYPOBJ = 'MAILLAGE'
  122. CALL LEKTAB(IPTABL,CHARIN,IOBRE)
  123. IF (IERR.NE.0) RETURN
  124. IPGEOM = IOBRE
  125. CALL LEKTAB(IPTABL,'ELTFA',IOBRE)
  126. IF (IERR.NE.0) RETURN
  127. IELTFA = IOBRE
  128. CALL LEKTAB(IPTABL,'CENTRE',IOBRE)
  129. IF (IERR.NE.0) RETURN
  130. ICENTR = IOBRE
  131. CALL LEKTAB(IPTABL,'FACE',IOBRE)
  132. IF (IERR.NE.0) RETURN
  133. IPFACE = IOBRE
  134. *
  135. * récup. MCHAML orientation normale
  136. *
  137. CALL LEKTAB(IPTABL,'XXNORMAE',IORIE)
  138. IF (IERR.NE.0) RETURN
  139. *
  140. * récup. CHPO orientation normale
  141. *
  142. CALL LEKTAB(IPTABL,'XXNORMAF',INORM)
  143. IF (IERR.NE.0) RETURN
  144. *
  145. * récup. CHPO surface des faces
  146. *
  147. CALL LEKTAB(IPTABL,'XXSURFAC',ISURF)
  148. IF (IERR.NE.0) RETURN
  149. *
  150. * Lecture de RIGIDITE obligatoire
  151. *
  152. CALL LIROBJ('RIGIDITE',IPRIGI,1,IRET)
  153. IF (IERR.NE.0) RETURN
  154. RI1 = IPRIGI
  155. *
  156. * Lecture du CHAMPOINT des charges (ou concentrations au centre)
  157. *
  158. CALL LIROBJ('CHPOINT',ITPP,1,IRET)
  159. *
  160. ICHP1 = ITPP
  161. *
  162. * Lecture du CHAMPOINT de concentration aux facees
  163. *
  164. CALL LIROBJ('CHPOINT',ITPN,1,IRET)
  165. IF (IERR.NE.0) RETURN
  166. *
  167. * Lecture eventuelle du 1er CHAMPOINT optionnel
  168. *
  169. ICHP2 = 0
  170. CALL LIROBJ('CHPOINT',ICHP2,0,IRET2)
  171. IF (IERR.NE.0) RETURN
  172. *
  173. * Lecture eventuelle du 2em CHAMPOINT optionnel
  174. *
  175. ICHP3 = 0
  176. CALL LIROBJ('CHPOINT',ICHP3,0,IRET3)
  177. IF (IERR.NE.0) RETURN
  178. *
  179. * Lecture de la RIGIDITE complementaire
  180. *
  181. IPRIGI = 0
  182. CALL LIROBJ('RIGIDITE',IPRIGI,0,IRETR)
  183. IF (IERR.NE.0) RETURN
  184. RI2 = IPRIGI
  185. *
  186. * Test de la formulation contenu dans MMODEL
  187. *
  188. SEGACT MMODEL
  189. NSOUS = KMODEL(/1)
  190. SEGINI IPMAHY
  191. IDARCY = 0
  192. IPT1= IPGEOM
  193. IPT2= IPGEOM
  194. DO 10 ISOUS=1,NSOUS
  195. IF(NSOUS.GT.1) THEN
  196. SEGACT IPT2
  197. IPT1= IPT2.LISOUS(ISOUS)
  198. SEGDES IPT2
  199. ENDIF
  200. IMODEL = KMODEL(ISOUS)
  201. SEGACT IMODEL
  202. LETYPE = FORMOD(1)
  203. IF (LETYPE.EQ.'DARCY') THEN
  204. IDARCY = IDARCY + 1
  205. MAHYBR(ISOUS) = IPT1
  206. ENDIF
  207. SEGDES IMODEL
  208. 10 CONTINUE
  209. SEGDES MMODEL
  210. IF (IDARCY.EQ.0) THEN
  211. MOTERR = LETYPE
  212. CALL ERREUR(193)
  213. GOTO 100
  214. ENDIF
  215. *
  216. * Récuperation des pointeurs ELTFA pour les zones ou DARCY est defini
  217. *
  218. MELEME = IELTFA
  219. SEGACT MELEME
  220. LZONES = LISOUS(/1)
  221. IF (LZONES.EQ.0) LZONES = 1
  222. IPT4 = IPGEOM
  223. SEGACT IPT4
  224. DO 30 ISOUS=1,NSOUS
  225. IMAMEL = MAHYBR(ISOUS)
  226. IF (IMAMEL.NE.0) THEN
  227. DO 20 ISZ=1,LZONES
  228. IF (LZONES.EQ.1) THEN
  229. IPT2 = IPT4
  230. IPT3 = MELEME
  231. ELSE
  232. IPT2 = IPT4.LISOUS(ISZ)
  233. IPT3 = LISOUS(ISZ)
  234. ENDIF
  235. IF (IPT2.EQ.IMAMEL) THEN
  236. MAHYBR(ISOUS) = IPT3
  237. GOTO 30
  238. ENDIF
  239. 20 CONTINUE
  240. IF (IMAMEL.EQ.MAHYBR(ISOUS)) THEN
  241. INTERR(1) = ISOUS
  242. CALL ERREUR(698)
  243. SEGDES IPT4
  244. SEGDES MELEME
  245. GOTO 100
  246. ENDIF
  247. ENDIF
  248. 30 CONTINUE
  249. SEGDES IPT4
  250. SEGDES MELEME
  251. *
  252. * Test et tri du sous-type des matrices de rigiditée récupérées
  253. *
  254. IPFORC=0
  255. IPDARC=0
  256. SEGACT RI1
  257. LETYPE = RI1.MTYMAT
  258. IF(LETYPE.EQ.'DARCY')THEN
  259. IPDARC=RI1
  260. ELSE IF(LETYPE.EQ.'MASSE')THEN
  261. IPFORC=RI1
  262. ENDIF
  263. IF(RI2.NE.0)THEN
  264. SEGACT RI2
  265. LETYPE = RI2.MTYMAT
  266. IF(LETYPE.EQ.'MASSE')THEN
  267. IPFORC=RI2
  268. ELSE IF(LETYPE.EQ.'DARCY')THEN
  269. IPDARC=RI2
  270. ENDIF
  271. IF(IPFORC.EQ.0)THEN
  272. MOTERR(1:8) = 'RIGIDITE'
  273. MOTERR(9:16) = 'MASSE'
  274. CALL ERREUR(79)
  275. SEGDES RI1
  276. SEGDES RI2
  277. GOTO 100
  278. ENDIF
  279. ENDIF
  280. IF(IPDARC.EQ.0)THEN
  281. MOTERR(1:8) = 'RIGIDITE'
  282. MOTERR(9:16) = 'DARCY'
  283. CALL ERREUR(79)
  284. SEGDES RI1
  285. IF(IPFORC.NE.0)SEGDES RI2
  286. GOTO 100
  287. ENDIF
  288. *
  289. * Controle des pointeurs de MELEME de la rigidité
  290. *
  291. DO 40 ISOUS=1,NSOUS
  292. IMAMEL = MAHYBR(ISOUS)
  293. IF (IMAMEL.NE.0) THEN
  294. RI1=IPDARC
  295. IPTEST = RI1.IRIGEL(1,ISOUS)
  296. CALL CMPMAI(IMAMEL,IPTEST,LOGRE)
  297. IF(.NOT. LOGRE)THEN
  298. MOTERR(1:8) = 'DARCY'
  299. MOTERR(9:16) = ' ELTFA '
  300. INTERR(1) = ISOUS
  301. CALL ERREUR(698)
  302. SEGDES RI1
  303. GOTO 100
  304. ENDIF
  305. IF(IPFORC.NE.0)THEN
  306. RI2=IPFORC
  307. IPTEST = RI2.IRIGEL(1,ISOUS)
  308. CALL CMPMAI(IMAMEL,IPTEST,LOGRE)
  309. IF(.NOT. LOGRE)THEN
  310. MOTERR(1:8) = 'MASSE'
  311. MOTERR(9:16) = ' ELTFA '
  312. INTERR(1) = ISOUS
  313. CALL ERREUR(698)
  314. SEGDES RI1
  315. SEGDES RI2
  316. GOTO 100
  317. ENDIF
  318. ENDIF
  319. ENDIF
  320. 40 CONTINUE
  321. SEGDES RI1
  322. IF(IPFORC.NE.0)SEGDES RI2
  323. *
  324. * Test du CHAMPOINT de valeurs aux faces
  325. *
  326. IFORC=0
  327. MCHPOI=ITPN
  328. SEGACT MCHPOI
  329. NBSOUS=IPCHP(/1)
  330. SEGDES MCHPOI
  331. ICHP0= ITPN
  332. IITH=0
  333. IVCONV=0
  334. IF(NBSOUS.NE.1)THEN
  335. C c'est le cas en presence de multiplicateurs de Lagrange
  336. C on va extraire la composante TH
  337. NOMTH='TH '
  338. CALL EXCOPP(ITPN,NOMTH,NIFOUR,ITPTH,NOMTH,NIFOUR,0)
  339. IF(IERR.NE.0)GO TO 100
  340. ICHP0= ITPTH
  341. IITH=1
  342. ENDIF
  343. NOMTOT(1)=' '
  344. NBCOMP=-1
  345. CALL QUEPOI(ICHP0,IPFACE,INDIC,NBCOMP,NOMTOT)
  346. IF(IERR.NE.0)GO TO 100
  347. IF(INDIC.LT.0)THEN
  348. IRETOU=0
  349. IPOINT=0
  350. CHARIN=' '
  351. CALL SKNAME(ITPN,CHARIN,IRETOU,IPOINT)
  352. MOTERR(1:8)= CHARIN(1:8)
  353. IF(IRETOU.EQ.1)CALL ERREUR(-301)
  354. CALL ERREUR(609)
  355. GO TO 100
  356. ENDIF
  357. IF(IITH.EQ.0)THEN
  358. IF(NBSOUS.EQ.1)THEN
  359. MCHPOI=ICHP0
  360. SEGACT MCHPOI
  361. MSOUPO=IPCHP(1)
  362. SEGDES MCHPOI
  363. SEGACT MSOUPO
  364. NC=NOCOMP(/2)
  365. IF(NC.EQ.1)THEN
  366. IF(NOCOMP(1).EQ.'TH ')IITH=1
  367. ENDIF
  368. SEGDES MSOUPO
  369. ENDIF
  370. ENDIF
  371. *
  372. * Tri et test des CHAMPOINTs ARGUMENTS optionnels
  373. *
  374. IF(IPFORC.NE.0) THEN
  375. IF(IITH.NE.1)THEN
  376. IRETOU=0
  377. IPOINT=0
  378. CHARIN=' '
  379. CALL SKNAME(ITPN,CHARIN,IRETOU,IPOINT)
  380. MOTERR(1:4)= 'TH '
  381. MOTERR(5:12)= CHARIN(1:8)
  382. CALL ERREUR(77)
  383. ENDIF
  384. IF(IRET2.EQ.1)THEN
  385. IFORC=ICHP2
  386. NOMTOT(1)=' '
  387. NBCOMP=-1
  388. CALL QUEPOI(ICHP1,ICENTR,INDIC,NBCOMP,NOMTOT)
  389. IF(IERR.NE.0)GO TO 100
  390. IF(INDIC.LT.0)THEN
  391. CALL ERREUR(348)
  392. GO TO 100
  393. ENDIF
  394. ELSE
  395. CALL ERREUR(641)
  396. GO TO 100
  397. ENDIF
  398. NBCOMP = 0
  399. INDIC = 1
  400. NOMTO3(1) = ' '
  401. NOMTO3(2) = ' '
  402. NOMTO3(3) = ' '
  403. CALL QUEPOI(IFORC,ICENTR,INDIC,NBCOMP,NOMTO3)
  404. IF(NOMTO3(1).NE.'FX')THEN
  405. MOTERR(1:4)= NOMTO3(1)
  406. CALL ERREUR(197)
  407. ENDIF
  408. IF(IRET3.NE.0)THEN
  409. NBCOMP = -1
  410. INDIC = 0
  411. NOMTOT(1) = ' '
  412. CALL QUEPOI(ICHP3,IPFACE,INDIC,NBCOMP,NOMTOT)
  413. IF(IERR.NE.0)GO TO 100
  414. IF(INDIC.LT.0)THEN
  415. CALL ERREUR(348)
  416. GO TO 100
  417. ENDIF
  418. IVCONV=ICHP3
  419. ENDIF
  420. C Cas ou il n y a pas de rigidite RI2
  421. ELSE
  422. NOMTOT(1)=' '
  423. NBCOMP=-1
  424. CALL QUEPOI(ICHP1,ICENTR,INDIC,NBCOMP,NOMTOT)
  425. IF(IERR.NE.0)GO TO 100
  426. IF(INDIC.LT.0)THEN
  427. CALL ERREUR(348)
  428. GO TO 100
  429. ENDIF
  430. MCHPOI=ICHP0
  431. SEGACT MCHPOI
  432. MSOUPO=IPCHP(1)
  433. SEGDES MCHPOI
  434. SEGACT MSOUPO
  435. NC=NOCOMP(/2)
  436. MCHPO1= ICHP1
  437. SEGACT MCHPO1
  438. MSOUP1=MCHPO1.IPCHP(1)
  439. SEGDES MCHPO1
  440. SEGACT MSOUP1
  441. NC1=MSOUP1.NOCOMP(/2)
  442. IF(NC1.NE.NC)THEN
  443. MOTERR(1:40)=
  444. & 'pas le meme nombre de composantes '
  445. CALL ERREUR(-301)
  446. CALL ERREUR(21)
  447. GO TO 100
  448. ENDIF
  449. IF(IITH.EQ.1) THEN
  450. IF((MSOUP1.NOCOMP(1)).NE.'H ')THEN
  451. MOTERR(1:40)=
  452. & 'le nom de la composante doit etre H '
  453. CALL ERREUR(-301)
  454. CALL ERREUR(21)
  455. GO TO 100
  456. ENDIF
  457. IF(IRET2.NE.0)THEN
  458. NBCOMP = -1
  459. INDIC = 0
  460. NOMTOT(1) = ' '
  461. CALL QUEPOI(ICHP2,IPFACE,INDIC,NBCOMP,NOMTOT)
  462. IF(IERR.NE.0)GO TO 100
  463. IF(INDIC.LT.0)THEN
  464. CALL ERREUR(348)
  465. GO TO 100
  466. ENDIF
  467. IVCONV=ICHP2
  468. ENDIF
  469. ELSE
  470. NNC1=0
  471. DO 50 I=1,NC
  472. DO 35 J=1,NC1
  473. IF(NOCOMP(I).EQ.MSOUP1.NOCOMP(J))THEN
  474. NNC1=NNC1+1
  475. GO TO 32
  476. ENDIF
  477. 35 CONTINUE
  478. MOTERR(1:40)=
  479. & 'les noms de composantes sont differents '
  480. CALL ERREUR(-301)
  481. CALL ERREUR(21)
  482. GO TO 100
  483. 32 CONTINUE
  484. 50 CONTINUE
  485. IF(NNC1.NE.NC1)THEN
  486. MOTERR(1:40)=
  487. & 'un des CHPO a plus de comp. que l autre '
  488. CALL ERREUR(21)
  489. GO TO 100
  490. ENDIF
  491. IVCONV=ICHP2
  492. ENDIF
  493. SEGDES MSOUPO,MSOUP1
  494. ENDIF
  495. *
  496. * Construction du CHAMPOINT des debits aux faces
  497. *
  498. SEGDES IPMAHY
  499. CALL HDEBI2(IPMAHY,IPFACE,IPDARC,ICHP0,IORIE,
  500. S ICHP1,IVCONV,IITH,INORM,IORIE,ISURF,IPFORC,IFORC,IRET)
  501.  
  502.  
  503. IF (IRET.NE.0) CALL ECROBJ('CHPOINT',IRET)
  504. *
  505. * Ménage
  506. *
  507. 100 CONTINUE
  508. SEGSUP IPMAHY
  509. RETURN
  510. END
  511.  
  512.  
  513.  
  514.  
  515.  
  516.  
  517.  
  518.  
  519.  
  520.  
  521.  
  522.  
  523.  
  524.  
  525.  
  526.  
  527.  
  528.  
  529.  
  530.  
  531.  
  532.  

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