Télécharger xlap13.eso

Retour à la liste

Numérotation des lignes :

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

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