Télécharger ylap13.eso

Retour à la liste

Numérotation des lignes :

ylap13
  1. C YLAP13 SOURCE OF166741 24/12/13 21:17:41 12097
  2. SUBROUTINE YLAP13(MU,KAPPA,CV,IROC,IVITC,IGRVC,
  3. & IGRTC,IVIMP,ITAUIM,IQIMP,
  4. & MELEMC,MELEMF,MELEFL,ISURF,INORM,IDIAM,
  5. & ICHFLU,DT)
  6. C
  7. C**** Traitement des conditions aux bords
  8. C
  9. C**** 'VIMP' : la vitesse imposé n'importe ou!
  10. C
  11. C 'QIMP' : flux de chaleur imposé n'importe ou
  12. C
  13. C 'TAUI' : tenseur de contraintes visqueux imposé n'importe ou
  14. C
  15.  
  16. IMPLICIT INTEGER(I-N)
  17.  
  18. -INC PPARAM
  19. -INC CCOPTIO
  20. -INC CCREEL
  21. -INC SMCHPOI
  22. -INC SMELEME
  23. -INC SMCOORD
  24. -INC SMLENTI
  25. -INC SMLMOTS
  26. C
  27. POINTEUR MPROC.MPOVAL, MPVITC.MPOVAL, MPGRVF.MPOVAL,
  28. & MPGRTF.MPOVAL,
  29. & MPVIMP.MPOVAL, MPTAUI.MPOVAL, MPQIMP.MPOVAL,
  30. & MPSURF.MPOVAL, MPNORM.MPOVAL, MPDIAM.MPOVAL,
  31. & MPFLUX.MPOVAL
  32. C
  33. POINTEUR MELEMC.MELEME,MELEMF.MELEME,MELEFL.MELEME
  34. POINTEUR MLCENT.MLENTI,MLEVIM.MLENTI,MLEQIM.MLENTI,MLETAI.MLENTI
  35. C
  36. INTEGER IROC,IVITC,IGRVC,IGRTC,IVIMP,ITAUIM,IQIMP
  37. & ,ISURF,INORM,IDIAM,ICHFLU
  38. & ,NFAC, NLCF, NGCF, NGCF1, NGCEG, NGCED
  39. & ,NLCEG,NLCED,NLFVI,NLFTI,NLFQI
  40. & , ICOORX, IGEOM
  41.  
  42. REAL*8 MU,KAPPA,CV,LAMBRO,DT, UNSDT
  43. & ,UXG,UYG
  44. & ,XG,YG,XFMXG,YFMYG,DRG
  45. & ,UXD,UYD
  46. & ,XD,YD,XFMXD,YFMYD,DRD,ALPHA,UMALPH
  47. & ,UXF,UYF,DUXXF,DUXYF,DUYXF,DUYYF,DTXF,DTYF
  48. & ,DSTDU,TAUXX,TAUXY,TAUYY,QX,QY,XF,YF
  49. & ,CNX, CNY, ORIENT, RO, DIAM, DIAM2, CELL, SURF
  50. & ,ZF,ZG,ZD,ZFMZD,UZF,UZG,UZD,TAUZZ,TAUXZ,TAUYZ
  51. & ,DUXZF,DUYZF,DUZXF,DUZYF,DUZZF,QZ,DTZF,CNZ,ZFMZG
  52. C
  53. CHARACTER*8 TYPE
  54. C
  55. C**** Initialisation de 1/DT
  56. C
  57. UNSDT = 0.0D0
  58. LAMBRO = KAPPA / CV
  59. C
  60. C**** KRIPAD pour la correspondance global/local de centre
  61. C
  62. CALL KRIPAD(MELEMC,MLCENT)
  63. C
  64. C EN KRIPAD
  65. C SEGACT MELEMC
  66. C SEGACT MLCENT
  67. C
  68. CALL LICHT(IROC,MPROC,TYPE,IGEOM)
  69. CALL LICHT(IVITC,MPVITC,TYPE,IGEOM)
  70. CALL LICHT(IGRVC,MPGRVF,TYPE,IGEOM)
  71. CALL LICHT(IGRTC,MPGRTF,TYPE,IGEOM)
  72. CALL LICHT(ISURF,MPSURF,TYPE,IGEOM)
  73. CALL LICHT(INORM,MPNORM,TYPE,IGEOM)
  74. CALL LICHT(IDIAM,MPDIAM,TYPE,IGEOM)
  75. CALL LICHT(ICHFLU,MPFLUX,TYPE,IGEOM)
  76. C
  77. C EN LICHT
  78. C SEGACT*MOD MPROC
  79. C SEGACT*MOD MPVITC
  80. C SEGACT*MOD MPGRVF
  81. C SEGACT*MOD MPGRTF
  82. C SEGACT*MOD MPSURF
  83. C SEGACT*MOD MPNORM
  84. C SEGACT*MOD MPDIAM
  85. C SEGACT*MOD MPFLUX
  86. C
  87. IF(IVIMP .GT. 0)THEN
  88. CALL LICHT(IVIMP,MPVIMP,TYPE,IGEOM)
  89. C SEGACT*MOD MPVIMP
  90. CALL KRIPAD(IGEOM,MLEVIM)
  91. C SEGACT IGEOM
  92. C SEGACT MLEVIM
  93. MELEME = IGEOM
  94. SEGDES MELEME
  95. ENDIF
  96. IF(ITAUIM .GT. 0)THEN
  97. CALL LICHT(ITAUIM,MPTAUI,TYPE,IGEOM)
  98. C SEGACT*MOD MPTAUI
  99. CALL KRIPAD(IGEOM,MLETAI)
  100. C SEGACT IGEOM
  101. C SEGACT MLETAI
  102. MELEME = IGEOM
  103. SEGDES MELEME
  104. ENDIF
  105. IF(IQIMP .GT. 0)THEN
  106. CALL LICHT(IQIMP,MPQIMP,TYPE,IGEOM)
  107. C SEGACT*MOD MPQIMP
  108. CALL KRIPAD(IGEOM,MLEQIM)
  109. C SEGACT IGEOM
  110. C SEGACT MLEQIM
  111. MELEME = IGEOM
  112. SEGDES MELEME
  113. ENDIF
  114. C
  115. SEGACT MELEFL
  116. SEGACT MELEMF
  117. NFAC = MELEMF.NUM(/2)
  118. C
  119. C**** Boucle sur les faces
  120. C
  121. DO NLCF = 1, NFAC, 1
  122. C
  123. C******* NLCF = numero local du centre de facel
  124. C NGCF = numero global du centre de facel
  125. C NLCF1 = numero local du centre de face
  126. C NGCEG = numero global du centre ELT "gauche"
  127. C NLCEG = numero local du centre ELT "gauche"
  128. C NGCED = numero global du centre ELT "droite"
  129. C NLCED = numero local du centre ELT "droite"
  130. C
  131. NGCF = MELEMF.NUM(1,NLCF)
  132. NGCF1 = MELEFL.NUM(2,NLCF)
  133. IF(NGCF .NE. NGCF1)THEN
  134. MOTERR(1:40)= 'FACEL et FACE = ? '
  135. WRITE(IOIMP,*) MOTERR(1:40)
  136. MOTERR(1:40)=
  137. & 'LAPN , subroutine ylap13.eso. '
  138. WRITE(IOIMP,*) MOTERR(1:40)
  139. CALL ERREUR(5)
  140. GOTO 9999
  141. ENDIF
  142. C
  143. NGCEG = MELEFL.NUM(1,NLCF)
  144. NGCED = MELEFL.NUM(3,NLCF)
  145. NLCEG = MLCENT.LECT(NGCEG)
  146. NLCED = MLCENT.LECT(NGCED)
  147. C
  148. C******* On controlle si sur NGCF on impose de CL
  149. C
  150. C NLFVI = numero local du centre de face sul le maillage des
  151. C "vitesses" "imposées"
  152. C
  153. C NLFTI = numero local du centre de face sul le maillage des
  154. C "tau" "imposés"
  155. C
  156. C NLFQI = numero local du centre de face sul le maillage des
  157. C "q" "imposés"
  158. C
  159. IF(IVIMP .GT. 0)THEN
  160. NLFVI = MLEVIM.LECT(NGCF)
  161. ELSE
  162. NLFVI = 0
  163. ENDIF
  164. C
  165. IF(ITAUIM .GT. 0)THEN
  166. NLFTI = MLETAI.LECT(NGCF)
  167. ELSE
  168. NLFTI = 0
  169. ENDIF
  170. C
  171. IF(IQIMP .GT. 0)THEN
  172. NLFQI = MLEQIM.LECT(NGCF)
  173. ELSE
  174. NLFQI = 0
  175. ENDIF
  176. C
  177. IF(NGCEG .NE. NGCED)THEN
  178. C
  179. C********** Parametres geometriques
  180. C
  181. ICOORX = ((IDIM + 1) * (NGCF - 1))+1
  182. XF = MCOORD.XCOOR(ICOORX)
  183. YF = MCOORD.XCOOR(ICOORX+1)
  184. ZF = MCOORD.XCOOR(ICOORX+2)
  185. C
  186. ICOORX = ((IDIM + 1) * (NGCEG - 1))+1
  187. XG = MCOORD.XCOOR(ICOORX)
  188. YG = MCOORD.XCOOR(ICOORX+1)
  189. ZG = MCOORD.XCOOR(ICOORX+2)
  190. XFMXG = XF - XG
  191. YFMYG = YF - YG
  192. ZFMZG = ZF - ZG
  193. DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG)+(ZFMZG*ZFMZG))
  194. C
  195. ICOORX = ((IDIM + 1) * (NGCED - 1))+1
  196. XD = MCOORD.XCOOR(ICOORX)
  197. YD = MCOORD.XCOOR(ICOORX+1)
  198. ZD = MCOORD.XCOOR(ICOORX+2)
  199. XFMXD = XF - XD
  200. YFMYD = YF - YD
  201. ZFMZD = ZF - ZD
  202. DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD)+(ZFMZD*ZFMZD))
  203. C
  204. C********** F=G -> DRG = 0 -> ALPHA = 0
  205. C
  206. ALPHA=DRG/(DRG+DRD)
  207. UMALPH=1.0D0 - ALPHA
  208. C
  209. C********** Les valeurs à l'interface
  210. C
  211. C DRG=0 -> F=G
  212. C
  213. C
  214. C********** Tenseur de contraintes visqueux
  215. C
  216. DUXXF = MPGRVF.VPOCHA(NLCF,1)
  217. DUXYF = MPGRVF.VPOCHA(NLCF,2)
  218. DUXZF = MPGRVF.VPOCHA(NLCF,3)
  219. DUYXF = MPGRVF.VPOCHA(NLCF,4)
  220. DUYYF = MPGRVF.VPOCHA(NLCF,5)
  221. DUYZF = MPGRVF.VPOCHA(NLCF,6)
  222. DUZXF = MPGRVF.VPOCHA(NLCF,7)
  223. DUZYF = MPGRVF.VPOCHA(NLCF,8)
  224. DUZZF = MPGRVF.VPOCHA(NLCF,9)
  225. C
  226. IF (NLFTI .GT. 0) THEN
  227. TAUXX = MPTAUI.VPOCHA(NLFTI,1)
  228. TAUYY = MPTAUI.VPOCHA(NLFTI,2)
  229. TAUZZ = MPTAUI.VPOCHA(NLFTI,3)
  230. TAUXY = MPTAUI.VPOCHA(NLFTI,4)
  231. TAUXZ = MPTAUI.VPOCHA(NLFTI,5)
  232. TAUYZ = MPTAUI.VPOCHA(NLFTI,6)
  233. ELSE
  234. DSTDU = 2.0D0/3.0D0 * (DUXXF + DUYYF + DUZZF)
  235. TAUXX = MU * (2.0D0 * DUXXF - DSTDU)
  236. TAUYY = MU * (2.0D0 * DUYYF - DSTDU)
  237. TAUZZ = MU * (2.0D0 * DUZZF - DSTDU)
  238. TAUXY = MU * (DUXYF + DUYXF)
  239. TAUXZ = MU * (DUXZF + DUZXF)
  240. TAUYZ = MU * (DUZYF + DUYZF)
  241. ENDIF
  242. C
  243. C********** Vitesse
  244. C
  245. IF( NLFVI .GT. 0) THEN
  246. UXF = MPVIMP.VPOCHA(NLFVI,1)
  247. UYF = MPVIMP.VPOCHA(NLFVI,2)
  248. UZF = MPVIMP.VPOCHA(NLFVI,3)
  249. ELSE
  250. C************* Gauche
  251. UXG = MPVITC.VPOCHA(NLCEG,1)
  252. UYG = MPVITC.VPOCHA(NLCEG,2)
  253. UZG = MPVITC.VPOCHA(NLCEG,3)
  254. C************* Droite
  255. UXD = MPVITC.VPOCHA(NLCED,1)
  256. UYD = MPVITC.VPOCHA(NLCED,2)
  257. UZD = MPVITC.VPOCHA(NLCED,3)
  258. C************* Face
  259. UXF = UMALPH * UXG + ALPHA * UXD
  260. UYF = UMALPH * UYG + ALPHA * UYD
  261. UZF = UMALPH * UZG + ALPHA * UZD
  262. C************* Correction de la vitesse lineaire exacte
  263. UXF = UXF +
  264. & (DUXXF * ((XFMXG * UMALPH)+ (XFMXD * ALPHA))) +
  265. & (DUXYF * ((YFMYG * UMALPH)+ (YFMYD * ALPHA))) +
  266. & (DUXZF * ((ZFMZG * UMALPH)+ (ZFMZD * ALPHA)))
  267. UYF = UYF +
  268. & (DUYXF * ((XFMXG * UMALPH)+ (XFMXD * ALPHA))) +
  269. & (DUYYF * ((YFMYG * UMALPH)+ (YFMYD * ALPHA))) +
  270. & (DUYZF * ((ZFMZG * UMALPH)+ (ZFMZD * ALPHA)))
  271. UZF = UZF +
  272. & (DUZXF * ((XFMXG * UMALPH)+ (XFMXD * ALPHA))) +
  273. & (DUZYF * ((YFMYG * UMALPH)+ (YFMYD * ALPHA))) +
  274. & (DUZZF * ((ZFMZG * UMALPH)+ (ZFMZD * ALPHA)))
  275. ENDIF
  276. C
  277. C********** Flux de chaleur
  278. C
  279. IF(NLFQI .GT. 0)THEN
  280. QX = MPQIMP.VPOCHA(NLFQI,1)
  281. QY = MPQIMP.VPOCHA(NLFQI,2)
  282. QZ = MPQIMP.VPOCHA(NLFQI,3)
  283. ELSE
  284. C************* Gauche
  285. DTXF = MPGRTF.VPOCHA(NLCF,1)
  286. DTYF = MPGRTF.VPOCHA(NLCF,2)
  287. DTZF = MPGRTF.VPOCHA(NLCF,3)
  288. C
  289. QX = -1.0D0 * KAPPA * DTXF
  290. QY = -1.0D0 * KAPPA * DTYF
  291. QZ = -1.0D0 * KAPPA * DTZF
  292. C
  293. ENDIF
  294. ELSE
  295. C
  296. C********** MURS
  297. C
  298. C Etat a gauche = Etat droite
  299. C
  300. ALPHA=0.0D0
  301. UMALPH=1.0D0
  302. C
  303. C********** Parametres geometriques
  304. C
  305. ICOORX = ((IDIM + 1) * (NGCF - 1))+1
  306. XF = MCOORD.XCOOR(ICOORX)
  307. YF = MCOORD.XCOOR(ICOORX+1)
  308. ZF = MCOORD.XCOOR(ICOORX+2)
  309. C
  310. ICOORX = ((IDIM + 1) * (NGCEG - 1))+1
  311. XG = MCOORD.XCOOR(ICOORX)
  312. YG = MCOORD.XCOOR(ICOORX+1)
  313. ZG = MCOORD.XCOOR(ICOORX+2)
  314. XFMXG = XF - XG
  315. YFMYG = YF - YG
  316. ZFMZG = ZF - ZG
  317. C
  318. C********** Tenseur de contraintes visqueux
  319. C
  320. DUXXF = MPGRVF.VPOCHA(NLCF,1)
  321. DUXYF = MPGRVF.VPOCHA(NLCF,2)
  322. DUXZF = MPGRVF.VPOCHA(NLCF,3)
  323. DUYXF = MPGRVF.VPOCHA(NLCF,4)
  324. DUYYF = MPGRVF.VPOCHA(NLCF,5)
  325. DUYZF = MPGRVF.VPOCHA(NLCF,6)
  326. DUZXF = MPGRVF.VPOCHA(NLCF,7)
  327. DUZYF = MPGRVF.VPOCHA(NLCF,8)
  328. DUZZF = MPGRVF.VPOCHA(NLCF,9)
  329. C
  330. IF (NLFTI .GT. 0) THEN
  331. TAUXX = MPTAUI.VPOCHA(NLFTI,1)
  332. TAUYY = MPTAUI.VPOCHA(NLFTI,2)
  333. TAUZZ = MPTAUI.VPOCHA(NLFTI,3)
  334. TAUXY = MPTAUI.VPOCHA(NLFTI,4)
  335. TAUXZ = MPTAUI.VPOCHA(NLFTI,5)
  336. TAUYZ = MPTAUI.VPOCHA(NLFTI,6)
  337. ELSE
  338. DSTDU = 2.0D0/3.0D0 * (DUXXF + DUYYF + DUZZF)
  339. TAUXX = MU * (2.0D0 * DUXXF - DSTDU)
  340. TAUYY = MU * (2.0D0 * DUYYF - DSTDU)
  341. TAUZZ = MU * (2.0D0 * DUZZF - DSTDU)
  342. TAUXY = MU * (DUXYF + DUYXF)
  343. TAUXZ = MU * (DUXZF + DUZXF)
  344. TAUYZ = MU * (DUZYF + DUYZF)
  345. ENDIF
  346. C
  347. C********** Vitesse
  348. C
  349. IF( NLFVI .GT. 0) THEN
  350. UXF = MPVIMP.VPOCHA(NLFVI,1)
  351. UYF = MPVIMP.VPOCHA(NLFVI,2)
  352. UZF = MPVIMP.VPOCHA(NLFVI,3)
  353. ELSE
  354. UXF = MPVITC.VPOCHA(NLCEG,1)
  355. UYF = MPVITC.VPOCHA(NLCEG,2)
  356. UZF = MPVITC.VPOCHA(NLCEG,3)
  357. C************* Correction de la vitesse lineaire exacte
  358. UXF = UXF +
  359. & (DUXXF * XFMXG ) +
  360. & (DUXYF * YFMYG ) +
  361. & (DUXZF * ZFMZG )
  362. UYF = UYF +
  363. & (DUYXF * XFMXG ) +
  364. & (DUYYF * YFMYG ) +
  365. & (DUYZF * ZFMZG )
  366. UZF = UZF +
  367. & (DUZXF * XFMXG ) +
  368. & (DUZYF * YFMYG ) +
  369. & (DUZZF * ZFMZG )
  370. ENDIF
  371. C
  372. C********** Flux de chaleur
  373. C
  374. IF(NLFQI .GT. 0)THEN
  375. QX = MPQIMP.VPOCHA(NLFQI,1)
  376. QY = MPQIMP.VPOCHA(NLFQI,2)
  377. QZ = MPQIMP.VPOCHA(NLFQI,3)
  378. ELSE
  379. C************* Gauche
  380. DTXF = MPGRTF.VPOCHA(NLCF,1)
  381. DTYF = MPGRTF.VPOCHA(NLCF,2)
  382. DTZF = MPGRTF.VPOCHA(NLCF,3)
  383. C
  384. QX = -1.0D0 * KAPPA * DTXF
  385. QY = -1.0D0 * KAPPA * DTYF
  386. QZ = -1.0D0 * KAPPA * DTZF
  387. C
  388. ENDIF
  389. ENDIF
  390. C
  391. C******* On calcule le sign du pruduit scalare
  392. C (Normales de Castem) * (vecteur "gauche" -> "centre")
  393. C
  394. CNX = MPNORM.VPOCHA(NLCF,1)
  395. CNY = MPNORM.VPOCHA(NLCF,2)
  396. CNZ = MPNORM.VPOCHA(NLCF,3)
  397. ORIENT = CNX * XFMXG + CNY * YFMYG + CNZ * ZFMZG
  398. ORIENT = SIGN(1.0D0,ORIENT)
  399. IF(ORIENT .NE. 1.0D0)THEN
  400. MOTERR(1:40)=
  401. & 'LAPN , subroutine ylap13.eso. '
  402. WRITE(IOIMP,*) MOTERR(1:40)
  403. MOTERR(1:40)=
  404. & 'Orientation normales. '
  405. WRITE(IOIMP,*) MOTERR(1:40)
  406. CALL ERREUR(5)
  407. GOTO 9999
  408. ENDIF
  409. C
  410. C******* Le flux aux interfaces
  411. C
  412. SURF = MPSURF.VPOCHA(NLCF,1)
  413. MPFLUX.VPOCHA(NLCF,1) = ((TAUXX * CNX) + (TAUXY * CNY) +
  414. & (TAUXZ * CNZ))
  415. & * SURF * (-1.0D0)
  416. MPFLUX.VPOCHA(NLCF,2) = ((TAUXY * CNX) + (TAUYY * CNY) +
  417. & (TAUYZ * CNZ))
  418. & * SURF * (-1.0D0)
  419. MPFLUX.VPOCHA(NLCF,3) = ((TAUXZ * CNX) + (TAUYZ * CNY) +
  420. & (TAUZZ * CNZ))
  421. & * SURF * (-1.0D0)
  422. MPFLUX.VPOCHA(NLCF,4) = (
  423. & ((TAUXX * UXF + TAUXY * UYF + TAUXZ * UZF - QX) * CNX) +
  424. & ((TAUXY * UXF + TAUYY * UYF + TAUYZ * UZF - QY) * CNY) +
  425. & ((TAUXZ * UXF + TAUYZ * UYF + TAUZZ * UZF - QZ) * CNZ))
  426. & * SURF * (-1.0D0)
  427. C
  428. C****** Le pas de temps
  429. C
  430. RO=UMALPH*MPROC.VPOCHA(NLCEG,1) +
  431. & ALPHA*MPROC.VPOCHA(NLCED,1)
  432. DIAM = UMALPH*MPDIAM.VPOCHA(NLCEG,1) +
  433. & ALPHA*MPDIAM.VPOCHA(NLCED,1)
  434. DIAM2=DIAM*DIAM
  435. CELL = 6.0D0*MU / (DIAM2*RO)
  436. CELL = MAX(CELL, (6.0D0*LAMBRO/(DIAM2*RO)))
  437. C
  438. IF(CELL .GT. UNSDT)THEN
  439. UNSDT = CELL
  440. ENDIF
  441. C
  442. ENDDO
  443. C
  444. C
  445. DT = 1.0D0 / (UNSDT + XPETIT)
  446. C
  447. SEGDES MELEFL
  448. SEGDES MELEMF
  449. SEGDES MELEMC
  450. SEGDES MPSURF
  451. SEGDES MPNORM
  452. SEGDES MPDIAM
  453. SEGSUP MLCENT
  454. C
  455. SEGDES MPROC
  456. SEGDES MPVITC
  457. SEGDES MPGRVF
  458. SEGDES MPGRTF
  459. SEGDES MPFLUX
  460. C
  461. IF(IVIMP .GT. 0) THEN
  462. SEGDES MPVIMP
  463. SEGSUP MLEVIM
  464. ENDIF
  465. IF(ITAUIM .GT. 0)THEN
  466. SEGDES MPTAUI
  467. SEGSUP MLETAI
  468. ENDIF
  469. IF(IQIMP .GT. 0)THEN
  470. SEGDES MPQIMP
  471. SEGDES MLEQIM
  472. ENDIF
  473. C
  474. 9999 CONTINUE
  475. RETURN
  476. END
  477.  
  478.  
  479.  
  480.  
  481.  
  482.  
  483.  
  484.  
  485.  

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