Télécharger ksupg.eso

Retour à la liste

Numérotation des lignes :

  1. C KSUPG SOURCE CHAT 05/01/13 01:09:19 5004
  2. C KSUPG SOURCE MAGN 97/07/19 21:18:18 2765
  3. SUBROUTINE KSUPG(FN,GR,PG,XYZ,HR,PGSQ,RPG,
  4. & NES,NP,NPG,IAXI,XCOOR,
  5. & WT,WS,HK,PGSK,RPGK,AIRE,
  6. & UMJ,DUMJ,KDEB,KFIN,LRV,NPX,NPGX,
  7. & TN,IPADT,UN,IPADU,NPTD,NEL,ANUK,
  8. & IDCEN,LE,
  9. & AL,AH,AP,
  10. & DTM1,DT,DTT1,DTT2,DIAEL,NUEL)
  11. C--------------------------------------------------------------------
  12. C Calcul des fonctions tests associées à la formulation Eléments
  13. C Finis Petrov-Galerkin afin de stabiliser les termes de convection.
  14. C On en profite pour calculer des temps caractéristiques.
  15. C--------------------------------------------------------------------
  16. C Cette subroutine étant écrite en fortran pur, certains arguments
  17. C servent uniquement au dimensionnement des tableaux.
  18. C Afin de doper les calculs, on decoupe la boucle sur les éléments
  19. C par paquets afin de bénéficier à fond de la vectorisation.
  20. C--------------------------------------------------------------------
  21. C
  22. C--------------------------
  23. C Paramètre Entrée/Sortie :
  24. C--------------------------
  25. C
  26. C /S FN : Fonction de forme associé à la transformation géométrique
  27. C /S GR : Gradient de FN dans l'élément de référence
  28. C /S PG : Poids d'intégration associé à chaque point de Gauss
  29. C /S XYZ : Coordonnées des sommets de l'élément
  30. C /S HR : Gradient des fonctions de forme dans l'élément courant
  31. C /S PGSQ : Produit du poids d'intégration par detJ
  32. C /S RPG : Rayon associé à chaque point de Gauss (cas axisymétrique)
  33. C
  34. C E/ NES : Dimension en espace du problème traité (2 en 2D et 3 en 3D)
  35. C E/ NP : Nombre de points support de ddl pour l'élément considéré
  36. C E/ NPG : Nombre de points de Gauss
  37. C E/ IAXI : Précise l'axe de symétrie dans le cas axisymétrique
  38. C (en fait l'axe de symétrie est toujours y -> IAXI=2)
  39. C E/ XCOOR : Coordonnées de l'ensemble des points (voir SMCOORD.INC)
  40. C Petrov-Galerkin choisi pour stabiliser la convection
  41. C Petrov-Galerkin choisi pour stabiliser la convection
  42. C /S HK : Gradient de FN dans l'élément courant
  43. C /S PGSK : Produit du poids d'intégration par detJ
  44. C /S RPGK : Rayon associé à chaque point de Gauss (cas axisymétrique)
  45. C /S AIRE : Volume de chaque élément
  46. C /S UMJ : Vitesse servant à la convection en chaque point de Gauss
  47. C /S DUMJ : Gradient de Vitesse servant en formulation conservative
  48. C E/ KDEB : Indice du début de la fenetre
  49. C E/ KFIN : Indice de fin de la fenetre
  50. C E/ LRV : Dimension de la fenetre
  51. C E/ NPX : Nombre de points maximum par élément
  52. C E/ NPGX : Nombre de points de Gauss maximum
  53. C E/ TN : Inconnue transportée au pas de temps précédent
  54. C E/ IPADT : Correspondance numérotation locale et globale pour les
  55. C points SOMMET : J=LECT(I) : le point numéro I est rangé
  56. C en Jème position pour le CHAMP TN
  57. C E/ UN : Champ de vitesse transportant aux sommets des éléments
  58. C E/ IPADU : Correspondance numérotation locale et globale pour les
  59. C points SOMMET : J=LECT(I) : le point numéro I est rangé
  60. C en Jème position pour le CHAMP UN
  61. C E/ NPTD : Nombre de points sous-tendant le spg associé à UN
  62. C E/ NEL : Nombre de points sous-tendant le spg associé à XXDXDY
  63. C E/ ANUK : Champ contenant le coefficient de diffusion (diffusivité
  64. C pour l'énergie, rapport viscosité sur densité pour qdm)
  65. C E/ IDCEN : Contient le schéma associé au terme convectif
  66. C (1=CENTREE 2=SUPGDC 3=SUPG 4=TVISQUEUX 5=CNG)
  67. C E/ LE : Maillage de la zone considérée
  68. C E/ AL : Première composante de XXDXDY (longueur carac. en x)
  69. C E/ AH : Deuxième composante de XXDXDY (longueur carac. en y)
  70. C E/ AP : Troiième composante de XXDXDY (longueur carac. en z)
  71. C E/ DTM1 : Pas de temps imposé
  72. C E/S DT : Pas de temps
  73. C /S DTT1 : Temps caractéristique associé à la diffusion/convection
  74. C /S DTT2 : Temps caractéristique associé à la diffusion
  75. C /S DIAEL : Longueur carac. de l'élément le plus pénalisant pour dt
  76. C /S NUEL : Numéro de l'élément le plus pénalisant pour dt
  77. C
  78. C--------------------------------------------------------------------
  79. IMPLICIT INTEGER(I-N)
  80. IMPLICIT REAL*8 (A-H,O-Z)
  81. -INC CCOPTIO
  82. -INC CCREEL
  83. DIMENSION IPADT(*),IPADU(*),LE(NP,NEL)
  84.  
  85. DIMENSION XCOOR(*)
  86. DIMENSION UN(NPTD,IDIM),TN(*)
  87. DIMENSION ANUK(LRV)
  88.  
  89. DIMENSION FN(NP,NPG),GR(IDIM,NP,NPG),PG(NPG),XYZ(IDIM,NP)
  90. DIMENSION HR(IDIM,NP,NPG),PGSQ(NPG),RPG(NPG)
  91. DIMENSION WT(LRV,NPX,NPGX),WS(LRV,NPX,NPGX),HK(LRV,3,NPX,NPGX)
  92. DIMENSION PGSK(LRV,NPGX),RPGK(LRV,NPGX),AIRE(LRV)
  93. C
  94. PARAMETER (LRV1=64,NPGX1=64)
  95. DIMENSION AL(LRV),AH(LRV),AP(LRV)
  96. DIMENSION UM(LRV1),UMI(LRV1,3),BMI(LRV1),TOI1(LRV1),TOI2(LRV1)
  97. DIMENSION GRAD(LRV1,3,NPGX1),UMJ(LRV,3,NPG),DUMJ(LRV,3,NPG)
  98. DIMENSION UP(LRV1,NPGX1),UPI(LRV1,3,NPGX1)
  99. DIMENSION BM(LRV1,NPGX1),BP(LRV1,NPGX1)
  100. DIMENSION TO1(LRV1,NPGX1),TO2(LRV1,NPGX1)
  101. DIMENSION CXT(LRV1,NPGX1),CYT(LRV1,NPGX1),CZT(LRV1,NPGX1)
  102. DIMENSION CXY(LRV1,NPGX1),CXZ(LRV1,NPGX1),CYZ(LRV1,NPGX1)
  103. DIMENSION AL2(LRV1),AH2(LRV1),AP2(LRV1),XMB(LRV1)
  104. C
  105. CC0 = 1.D0
  106. IF (LRV.NE.LRV1) CALL ARRET(0)
  107. IF (NPG.GT.NPGX1) CALL ARRET(0)
  108. C
  109. C- Pour les éléments de la fenetre KDEB:KFIN initialisations des
  110. C- données associées à l'élément fini utilisé : fonctions de forme,
  111. C- gradients des fonctions de forme, aire de l'élément, poids des
  112. C- points de Gauss et rayon des points de Gauss pour le cas axi.
  113. C
  114. DO 100 K=KDEB,KFIN
  115. KP = K - KDEB + 1
  116. DO 20 I=1,NP
  117. J = LE(I,K)
  118. DO 10 N=1,IDIM
  119. XYZ(N,I) = XCOOR((J-1)*(IDIM+1)+N)
  120. 10 CONTINUE
  121. 20 CONTINUE
  122. CALL CALJBC(FN,GR,PG,XYZ,HR,PGSQ,RPG,
  123. & NES,IDIM,NP,NPG,IAXI,AIRE(KP))
  124. C
  125. DO 50 L=1,NPG
  126. DO 40 N=1,IDIM
  127. DO 30 I=1,NP
  128. HK(KP,N,I,L) = HR(N,I,L)
  129. 30 CONTINUE
  130. 40 CONTINUE
  131. PGSK(KP,L) = PGSQ(L)
  132. RPGK(KP,L) = RPG(L)
  133. 50 CONTINUE
  134. 100 CONTINUE
  135. C
  136. C- Initialisations des variables
  137. C
  138. DO 120 N=1,IDIM
  139. DO 110 K=KDEB,KFIN
  140. KP = K - KDEB + 1
  141. UMI(KP,N) = XPETIT
  142. 110 CONTINUE
  143. 120 CONTINUE
  144. C
  145. DO 150 L=1,NPG
  146. DO 140 N=1,IDIM
  147. DO 130 K=KDEB,KFIN
  148. KP = K - KDEB + 1
  149. UMJ(KP,N,L) = XPETIT
  150. DUMJ(KP,N,L) = XPETIT
  151. UPI(KP,N,L) = XPETIT
  152. GRAD(KP,N,L) = XPETIT
  153. 130 CONTINUE
  154. 140 CONTINUE
  155. 150 CONTINUE
  156. C
  157. C- Calcul en chaque élément, pour chaque point de Gauss
  158. C- UMJ : Vitesse
  159. C- GRAD : Gradient de l'inconnue scalaire au pas de temps précédent
  160. C- Le gradient n'est à calculer que pour l'option SUPGDC
  161. C
  162. IF (IDCEN.EQ.2) THEN
  163. DO 190 N=1,IDIM
  164. DO 180 L=1,NPG
  165. DO 170 I=1,NP
  166. DO 160 K=KDEB,KFIN
  167. KP = K - KDEB + 1
  168. NF = IPADU(LE(I,K))
  169. NT = IPADT(LE(I,K))
  170. UMJ(KP,N,L) = UMJ(KP,N,L) + UN(NF,N)*FN(I,L)
  171. DUMJ(KP,N,L) = DUMJ(KP,N,L) + UN(NF,N)*HK(KP,N,I,L)
  172. GRAD(KP,N,L) = GRAD(KP,N,L) + TN(NT) *HK(KP,N,I,L)
  173. 160 CONTINUE
  174. 170 CONTINUE
  175. 180 CONTINUE
  176. 190 CONTINUE
  177. ELSE
  178. DO 230 N=1,IDIM
  179. DO 220 L=1,NPG
  180. DO 210 I=1,NP
  181. DO 200 K=KDEB,KFIN
  182. KP = K - KDEB + 1
  183. NF = IPADU(LE(I,K))
  184. UMJ(KP,N,L) = UMJ(KP,N,L) + UN(NF,N)*FN(I,L)
  185. DUMJ(KP,N,L) = DUMJ(KP,N,L) + UN(NF,N)*HK(KP,N,I,L)
  186. 200 CONTINUE
  187. 210 CONTINUE
  188. 220 CONTINUE
  189. 230 CONTINUE
  190. ENDIF
  191. C
  192. C- Calcul pour chaque élément de
  193. C- UMI : Vitesse moyenne
  194. C- pour les options SUPGDC et SUPG
  195. C
  196. IF (IDCEN.EQ.2.OR.IDCEN.EQ.3) THEN
  197. DO 260 L=1,NPG
  198. DO 250 N=1,IDIM
  199. DO 240 K=KDEB,KFIN
  200. KP = K - KDEB + 1
  201. UMI(KP,N) = UMI(KP,N) + UMJ(KP,N,L)*PGSK(KP,L)/AIRE(KP)
  202. 240 CONTINUE
  203. 250 CONTINUE
  204. 260 CONTINUE
  205. ENDIF
  206. C
  207. C
  208. C====================================================================
  209. C Cas 2D
  210. C====================================================================
  211. C
  212. C
  213. IF (IDIM.EQ.2) THEN
  214. C
  215. C- Calcul pour chaque élément de
  216. C- UM : Norme de la vitesse moyenne
  217. C- BMI : Inverse d'un temps caractéristique associé à la convection
  218. C- pour l'option SUPGDC et SUPG
  219. C
  220. IF (IDCEN.EQ.2.OR.IDCEN.EQ.3) THEN
  221. DO 1000 K=KDEB,KFIN
  222. KP = K - KDEB + 1
  223. UM(KP) = UMI(KP,1)*UMI(KP,1) + UMI(KP,2)*UMI(KP,2)
  224. UM(KP) = SQRT(UM(KP))
  225. BMIX = UMI(KP,1) / AL(KP)
  226. BMIY = UMI(KP,2) / AH(KP)
  227. BMI(KP) = BMIX*BMIX + BMIY*BMIY
  228. BMI(KP) = SQRT(BMI(KP)) + XPETIT
  229. 1000 CONTINUE
  230. ENDIF
  231. C
  232. C- Calcul de facteur géométrique associé à chaque élément
  233. C
  234. DO 1010 K=KDEB,KFIN
  235. KP=K-KDEB+1
  236. AL2(KP) = 1.D0 / AL(KP) / AL(KP)
  237. AH2(KP) = 1.D0 / AH(KP) / AH(KP)
  238. XMB(KP) = (AL(KP)+AH(KP)) / 2.D0
  239. 1010 CONTINUE
  240. C
  241. C- Calcul en chaque élément, pour chaque point de Gauss de
  242. C- GRAD : vecteur unitaire porté par le gradient du champ scalaire
  243. C- UP : projection de la vitesse sur la direction donnée par GRAD
  244. C- UPI : vecteur UP*GRAD
  245. C- BM : inverse d'un temps caractéristique basé sur la vitesse
  246. C- BP : idem BP mais en considérant la projection de V sur GRAD
  247. C- pour l'option SUPGDC
  248. C
  249. IF (IDCEN.EQ.2) THEN
  250. DO 1030 L=1,NPG
  251. DO 1020 K=KDEB,KFIN
  252. KP = K - KDEB + 1
  253. AX = GRAD(KP,1,L)*GRAD(KP,1,L) + GRAD(KP,2,L)*GRAD(KP,2,L)
  254. AX = SQRT(AX) + XPETIT
  255. GRAD(KP,1,L) = GRAD(KP,1,L) / AX
  256. GRAD(KP,2,L) = GRAD(KP,2,L) / AX
  257. UP(KP,L) = GRAD(KP,1,L) * UMJ(KP,1,L)
  258. & + GRAD(KP,2,L) * UMJ(KP,2,L)
  259. UPI(KP,1,L) = GRAD(KP,1,L) * UP(KP,L)
  260. UPI(KP,2,L) = GRAD(KP,2,L) * UP(KP,L)
  261. 1020 CONTINUE
  262. 1030 CONTINUE
  263. DO 1050 L=1,NPG
  264. DO 1040 K=KDEB,KFIN
  265. KP = K - KDEB + 1
  266. BMX = UMJ(KP,1,L) / AL(KP)
  267. BMY = UMJ(KP,2,L) / AH(KP)
  268. BPX = UPI(KP,1,L) / AL(KP)
  269. BPY = UPI(KP,2,L) / AH(KP)
  270. BM(KP,L) = BMX*BMX + BMY*BMY
  271. BP(KP,L) = BPX*BPX + BPY*BPY
  272. BM(KP,L) = SQRT(BM(KP,L)) + XPETIT
  273. BP(KP,L) = SQRT(BP(KP,L)) + XPETIT
  274. 1040 CONTINUE
  275. 1050 CONTINUE
  276. ENDIF
  277. C
  278. C-----------------------------
  279. C- DECENTREMENT suivant IDCEN
  280. C-----------------------------
  281. C On calcule dans chaque cas TO1 et TO2 ainsi que le tenseur
  282. C associé à la viscosité numérique afin d'évaluer le pas de
  283. C temps de stabilité pour les schémas explicites.
  284. C
  285. C----------
  286. C CENTREE :
  287. C----------
  288. IF (IDCEN.EQ.1) THEN
  289. SI1 = 0.D0
  290. SI2 = 0.D0
  291. DO 1110 L=1,NPG
  292. DO 1100 K=KDEB,KFIN
  293. KP = K - KDEB + 1
  294. TO1(KP,L) = 0.D0
  295. TO2(KP,L) = 0.D0
  296. CXT(KP,L) = 0.D0
  297. CYT(KP,L) = 0.D0
  298. CXY(KP,L) = 0.D0
  299. 1100 CONTINUE
  300. 1110 CONTINUE
  301. C---------
  302. C SUPGDC : Base théorique dans : A New FE formulation for computational
  303. C fluid dynamics, II Beyond SUPG, HUGHES et al., in Comp.Meth.Appl.Mech.
  304. C Eng., vol 54, pp 341-355 (1986).
  305. C---------
  306. ELSEIF (IDCEN.EQ.2) THEN
  307. SI1 = 1.D0
  308. SI2 = 1.D0
  309. DO 1130 L=1,NPG
  310. DO 1120 K=KDEB,KFIN
  311. KP = K - KDEB + 1
  312. C
  313. C- Approximation "doublement asymptotique" basée sur la vitesse moyenne
  314. C- HMK : Distance basé sur la vitesse moyenne
  315. C- ALFA : Peclet de maille basé sur la vitesse moyenne
  316. C
  317. HMK = 2.D0 * UM(KP) / BM(KP,L)
  318. ALFA = UM(KP) * HMK / ANUK(KP) / 2.D0
  319. IF (ALFA.GT.3.D0) THEN
  320. AKSI = 1.D0
  321. ELSE
  322. AKSI = ALFA / 3.D0
  323. ENDIF
  324. CCT = AKSI / BM(KP,L)
  325. XMB(KP) = HMK
  326. C
  327. C- Approximation "doublement asymptotique" basée sur la projection de
  328. C- la vitesse sur le gradient du champ scalaire
  329. C- HMK : Distance basée sur la vitesse projetée
  330. C- ALFA : Peclet de maille basé sur la vitesse projetée
  331. C
  332. HPK = 2.D0 * UP(KP,L) / BP(KP,L)
  333. ALFA = UP(KP,L) * HPK / ANUK(KP) / 2.D0
  334. IF (ALFA.GT.3.D0) THEN
  335. AKSI = 1.D0
  336. ELSE
  337. AKSI = ALFA / 3.D0
  338. ENDIF
  339. CCP = AKSI / BP(KP,L)
  340. C
  341. CPT = CCP - CCT
  342. IF (CPT.GE.0.D0) then
  343. CC2 = CPT * CC0
  344. ELSE
  345. CC2 = 0.D0
  346. ENDIF
  347. C
  348. TO1(KP,L) = CCT
  349. TO2(KP,L) = CC2
  350. CXT(KP,L) = UMJ(KP,1,L)*UMJ(KP,1,L)*CCT
  351. & + UPI(KP,1,L)*UPI(KP,1,L)*CC2
  352. CYT(KP,L) = UMJ(KP,2,L)*UMJ(KP,2,L)*CCT
  353. & + UPI(KP,2,L)*UPI(KP,2,L)*CC2
  354. CXY(KP,L) = UMJ(KP,1,L)*UMJ(KP,2,L)*CCT
  355. & + UPI(KP,1,L)*UPI(KP,2,L)*CC2
  356.  
  357. 1120 CONTINUE
  358. 1130 CONTINUE
  359. C-------
  360. C SUPG :
  361. C-------
  362. ELSEIF (IDCEN.EQ.3) THEN
  363. SI1 = 1.D0
  364. SI2 = 1.D0
  365. DO 1150 L=1,NPG
  366. DO 1140 K=KDEB,KFIN
  367. KP = K - KDEB + 1
  368. C
  369. C- Approximation "doublement asymptotique" basée sur la vitesse moyenne
  370. C- HMK : Distance basé sur la vitesse moyenne
  371. C- ALFA : Peclet de maille basé sur la vitesse moyenne
  372. C
  373. HMK = 2.D0 * UM(KP) / BMI(KP)
  374. ALFA = UM(KP) * HMK / ANUK(KP) / 2.D0
  375. IF (ALFA.GT.3.D0) THEN
  376. AKSI = 1.D0
  377. ELSE
  378. AKSI = ALFA / 3.D0
  379. ENDIF
  380. CCT = AKSI / BMI(KP)
  381. XMB(KP) = HMK
  382. C
  383. TO1(KP,L) = CCT
  384. TO2(KP,L) = 0.D0
  385. CXT(KP,L) = UMJ(KP,1,L)*UMJ(KP,1,L)*CCT
  386. CYT(KP,L) = UMJ(KP,2,L)*UMJ(KP,2,L)*CCT
  387. CXY(KP,L) = UMJ(KP,1,L)*UMJ(KP,2,L)*CCT
  388. 1140 CONTINUE
  389. 1150 CONTINUE
  390. C-------------------
  391. C Tenseur Visqueux :
  392. C-------------------
  393. ELSEIF (IDCEN.EQ.4) THEN
  394. SI1 =-1.D0
  395. SI2 = 1.D0
  396. DT19 = DTM1 * 0.5D0
  397. DO 1170 L=1,NPG
  398. DO 1160 K=KDEB,KFIN
  399. KP = K - KDEB + 1
  400. TO1(KP,L) = DT19
  401. TO2(KP,L) = 0.D0
  402. CXT(KP,L) = UMJ(KP,1,L)*UMJ(KP,1,L)*DT19
  403. CYT(KP,L) = UMJ(KP,2,L)*UMJ(KP,2,L)*DT19
  404. CXY(KP,L) = UMJ(KP,1,L)*UMJ(KP,2,L)*DT19
  405. 1160 CONTINUE
  406. 1170 CONTINUE
  407. C-----------------------------
  408. C Crank Nicholson généralisé :
  409. C-----------------------------
  410. ELSEIF (IDCEN.EQ.5) THEN
  411. SI1 =-1.D0
  412. SI2 = 1.D0
  413. DT19 = DTM1/6.D0
  414. DO 1190 L=1,NPG
  415. DO 1180 K=KDEB,KFIN
  416. KP = K - KDEB + 1
  417. TO1(KP,L) = DT19
  418. TO2(KP,L) = 0.D0
  419. 1180 CONTINUE
  420. 1190 CONTINUE
  421. ENDIF
  422. C
  423. C---------------------------
  424. C Pas de temps de stabilité
  425. C---------------------------
  426. C
  427. C La viscosité utilisée dans l'évaluation des pas de temps de stabilité
  428. C explicite ajoute à la viscosité physique la viscosité numérique :
  429. C DT1 : Correspond à une CFL de 1 et un Peclet de maille de 1
  430. C CFL=1 -> dx=Vdt |
  431. C |-> dt=2D/V2
  432. C Pe=1 -> dx=2D/V |
  433. C DT2 : Correspond au pas de temps de stabilité lié à la diffusion
  434. C FOU=1 -> dt=0.5 dx2/D
  435. C
  436. IF (IDCEN.NE.5) THEN
  437. DO 2010 L=1,NPG
  438. DO 2000 K=KDEB,KFIN
  439. KP = K - KDEB + 1
  440. DT0 = DT
  441. DT1 = 2.D0 /
  442. & ( UMJ(KP,1,L)*UMJ(KP,1,L)/(CXT(KP,L)+ANUK(KP))
  443. & + UMJ(KP,2,L)*UMJ(KP,2,L)/(CYT(KP,L)+ANUK(KP)) )
  444. DT2 = 0.5D0 /
  445. & ( (CXT(KP,L)+ANUK(KP))*AL2(KP)
  446. & + (CYT(KP,L)+ANUK(KP))*AH2(KP) )
  447. IF (DT1.LT.DT) DT=DT1
  448. IF (DT2.LT.DT) DT=DT2
  449. IF (DT.NE.DT0) THEN
  450. DTT1 = DT1
  451. DTT2 = DT2
  452. DIAEL = XMB(KP)
  453. NUEL = K
  454. ENDIF
  455. 2000 CONTINUE
  456. 2010 CONTINUE
  457. ENDIF
  458. C
  459. C---------------------------------------------------
  460. C Fonction test pour la formulation Petrov-Galerkin
  461. C---------------------------------------------------
  462. C Ce qui est diffusion numérique en explicite se transforme en
  463. C ajoutant de la viscosité numérique (Tenseurs visqueux et CNG).
  464. C WS : Fonction test pour la partie explicite
  465. C
  466. DO 2050 L=1,NPG
  467. DO 2050 I=1,NP
  468. DO 2050 K=KDEB,KFIN
  469. KP = K - KDEB + 1
  470. WT(KP,I,L) = FN(I,L) + SI1 *
  471. & (TO1(KP,L)*(UMJ(KP,1,L)*HK(KP,1,I,L)+UMJ(KP,2,L)*HK(KP,2,I,L))
  472. & +TO2(KP,L)*(UPI(KP,1,L)*HK(KP,1,I,L)+UPI(KP,2,L)*HK(KP,2,I,L)))
  473. WS(KP,I,L) = FN(I,L) + SI2 *
  474. & (TO1(KP,L)*(UMJ(KP,1,L)*HK(KP,1,I,L)+UMJ(KP,2,L)*HK(KP,2,I,L))
  475. & +TO2(KP,L)*(UPI(KP,1,L)*HK(KP,1,I,L)+UPI(KP,2,L)*HK(KP,2,I,L)))
  476. 2050 CONTINUE
  477. C
  478. C
  479. C====================================================================
  480. C Cas 3D
  481. C====================================================================
  482. C
  483. C
  484. ELSEIF (IDIM.EQ.3) THEN
  485. IF (IDCEN.EQ.2.OR.IDCEN.EQ.3) THEN
  486. DO 5000 K=KDEB,KFIN
  487. KP = K - KDEB + 1
  488. UM(KP) = UMI(KP,1)*UMI(KP,1)
  489. & + UMI(KP,2)*UMI(KP,2)
  490. & + UMI(KP,3)*UMI(KP,3)
  491. UM(KP) = SQRT(UM(KP))
  492. BMIX = UMI(KP,1) / AL(KP)
  493. BMIY = UMI(KP,2) / AH(KP)
  494. BMIZ = UMI(KP,3) / AP(KP)
  495. BMI(KP) = BMIX*BMIX + BMIY*BMIY + BMIZ*BMIZ
  496. BMI(KP) = SQRT(BMI(KP)) + XPETIT
  497. 5000 CONTINUE
  498. ENDIF
  499. *
  500. DO 5010 K=KDEB,KFIN
  501. KP=K-KDEB+1
  502. AL2(KP) = 1.D0 / AL(KP) / AL(KP)
  503. AH2(KP) = 1.D0 / AH(KP) / AH(KP)
  504. AP2(KP) = 1.D0 / AP(KP) / AP(KP)
  505. XMB(KP) = (AL(KP)+AH(KP)+AP(KP)) / 3.D0
  506. 5010 CONTINUE
  507. *
  508. IF (IDCEN.EQ.2) THEN
  509. DO 5030 L=1,NPG
  510. DO 5020 K=KDEB,KFIN
  511. KP = K-KDEB+1
  512. AX = GRAD(KP,1,L)*GRAD(KP,1,L)
  513. & + GRAD(KP,2,L)*GRAD(KP,2,L)
  514. & + GRAD(KP,3,L)*GRAD(KP,3,L)
  515. AX = SQRT(AX) + XPETIT
  516. GRAD(KP,1,L) = GRAD(KP,1,L) / AX
  517. GRAD(KP,2,L) = GRAD(KP,2,L) / AX
  518. GRAD(KP,3,L) = GRAD(KP,3,L) / AX
  519. UP(KP,L) = GRAD(KP,1,L) * UMJ(KP,1,L)
  520. & + GRAD(KP,2,L) * UMJ(KP,2,L)
  521. & + GRAD(KP,3,L) * UMJ(KP,3,L)
  522. UPI(KP,1,L) = GRAD(KP,1,L) * UP(KP,L)
  523. UPI(KP,2,L) = GRAD(KP,2,L) * UP(KP,L)
  524. UPI(KP,3,L) = GRAD(KP,3,L) * UP(KP,L)
  525. 5020 CONTINUE
  526. 5030 CONTINUE
  527. DO 5050 L=1,NPG
  528. DO 5040 K=KDEB,KFIN
  529. KP = K-KDEB+1
  530. BMX = UMJ(KP,1,L) / AL(KP)
  531. BMY = UMJ(KP,2,L) / AH(KP)
  532. BMZ = UMJ(KP,3,L) / AP(KP)
  533. BPX = UPI(KP,1,L) / AL(KP)
  534. BPY = UPI(KP,2,L) / AH(KP)
  535. BPZ = UPI(KP,3,L) / AP(KP)
  536. BM(KP,L) = BMX*BMX + BMY*BMY + BMZ*BMZ
  537. BP(KP,L) = BPX*BPX + BPY*BPY + BPZ*BPZ
  538. BM(KP,L) = SQRT(BM(KP,L)) + XPETIT
  539. BP(KP,L) = SQRT(BP(KP,L)) + XPETIT
  540. 5040 CONTINUE
  541. 5050 CONTINUE
  542. ENDIF
  543. C
  544. C-----------------------------
  545. C- DECENTREMENT suivant IDCEN
  546. C-----------------------------
  547. C
  548. C----------
  549. C CENTREE :
  550. C----------
  551. IF (IDCEN.EQ.1) THEN
  552. SI1= 0.D0
  553. SI2= 0.D0
  554. DO 5110 L=1,NPG
  555. DO 5100 K=KDEB,KFIN
  556. KP = K - KDEB + 1
  557. TO1(KP,L) = 0.D0
  558. TO2(KP,L) = 0.D0
  559. CXT(KP,L) = 0.D0
  560. CYT(KP,L) = 0.D0
  561. CZT(KP,L) = 0.D0
  562. CXY(KP,L) = 0.D0
  563. CXZ(KP,L) = 0.D0
  564. CYZ(KP,L) = 0.D0
  565. 5100 CONTINUE
  566. 5110 CONTINUE
  567. C---------
  568. C SUPGDC :
  569. C---------
  570. ELSEIF (IDCEN.EQ.2) THEN
  571. SI1= 1.D0
  572. SI2= 1.D0
  573. DO 5130 L=1,NPG
  574. DO 5120 K=KDEB,KFIN
  575. KP = K - KDEB + 1
  576. HMK = 2.D0 * UM(KP) / BM(KP,L)
  577. ALFA = UM(KP)*HMK / ANUK(KP) / 2.D0
  578. IF (ALFA.GT.3.D0) THEN
  579. AKSI = 1.D0
  580. ELSE
  581. AKSI = ALFA / 3.D0
  582. ENDIF
  583. CCT = AKSI / BM(KP,L)
  584. XMB(KP) = HMK
  585. *
  586. HPK = 2.D0 * UP(KP,L) / BP(KP,L)
  587. ALFA = UP(KP,L)*HPK / ANUK(KP) / 2.D0
  588. IF (ALFA.GT.3.D0) THEN
  589. AKSI = 1.D0
  590. ELSE
  591. AKSI = ALFA / 3.D0
  592. ENDIF
  593. CCP = AKSI / BP(KP,L)
  594. *
  595. CPT=CCP-CCT
  596. IF (CPT.GE.0.D0) then
  597. CC2 = CPT * CC0
  598. ELSE
  599. CC2 = 0.D0
  600. ENDIF
  601. *
  602. TO1(KP,L) = CCT
  603. TO2(KP,L) = CC2
  604. CXT(KP,L) = UMJ(KP,1,L)*UMJ(KP,1,L)*CCT
  605. & + UPI(KP,1,L)*UPI(KP,1,L)*CC2
  606. CYT(KP,L) = UMJ(KP,2,L)*UMJ(KP,2,L)*CCT
  607. & + UPI(KP,2,L)*UPI(KP,2,L)*CC2
  608. CZT(KP,L) = UMJ(KP,3,L)*UMJ(KP,3,L)*CCT
  609. & + UPI(KP,3,L)*UPI(KP,3,L)*CC2
  610. CXY(KP,L) = UMJ(KP,1,L)*UMJ(KP,2,L)*CCT
  611. & + UPI(KP,1,L)*UPI(KP,2,L)*CC2
  612. CXZ(KP,L) = UMJ(KP,1,L)*UMJ(KP,3,L)*CCT
  613. & + UPI(KP,1,L)*UPI(KP,3,L)*CC2
  614. CYZ(KP,L) = UMJ(KP,2,L)*UMJ(KP,3,L)*CCT
  615. & + UPI(KP,2,L)*UPI(KP,3,L)*CC2
  616.  
  617. 5120 CONTINUE
  618. 5130 CONTINUE
  619. C-------
  620. C SUPG :
  621. C-------
  622. ELSEIF (IDCEN.EQ.3) THEN
  623. SI1= 1.D0
  624. SI2= 1.D0
  625. DO 5150 L=1,NPG
  626. DO 5140 K=KDEB,KFIN
  627. KP = K - KDEB + 1
  628. HMK = 2.D0 * UM(KP) / BM(KP,L)
  629. ALFA = UM(KP)*HMK / ANUK(KP) / 2.D0
  630. IF (ALFA.GT.3.D0) THEN
  631. AKSI = 1.D0
  632. ELSE
  633. AKSI = ALFA / 3.D0
  634. ENDIF
  635. CCT = AKSI / BM(KP,L)
  636. XMB(KP) = HMK
  637. TO1(KP,L) = CCT
  638. TO2(KP,L) = 0.D0
  639. CXT(KP,L) = UMJ(KP,1,L)*UMJ(KP,1,L)*CCT
  640. CYT(KP,L) = UMJ(KP,2,L)*UMJ(KP,2,L)*CCT
  641. CZT(KP,L) = UMJ(KP,3,L)*UMJ(KP,3,L)*CCT
  642. CXY(KP,L) = UMJ(KP,1,L)*UMJ(KP,2,L)*CCT
  643. CXZ(KP,L) = UMJ(KP,1,L)*UMJ(KP,3,L)*CCT
  644. CYZ(KP,L) = UMJ(KP,2,L)*UMJ(KP,3,L)*CCT
  645. 5140 CONTINUE
  646. 5150 CONTINUE
  647. C-------------------
  648. C Tenseur Visqueux :
  649. C-------------------
  650. ELSEIF (IDCEN.EQ.4) THEN
  651. SI1 =-1.D0
  652. SI2 = 1.D0
  653. DT19 = DTM1 * 0.5D0
  654. DO 5170 L=1,NPG
  655. DO 5160 K=KDEB,KFIN
  656. KP = K - KDEB + 1
  657. TO1(KP,L) = DT19
  658. TO2(KP,L) = 0.D0
  659. CXT(KP,L) = UMJ(KP,1,L)*UMJ(KP,1,L)*DT19
  660. CYT(KP,L) = UMJ(KP,2,L)*UMJ(KP,2,L)*DT19
  661. CZT(KP,L) = UMJ(KP,3,L)*UMJ(KP,3,L)*DT19
  662. CXY(KP,L) = UMJ(KP,1,L)*UMJ(KP,2,L)*DT19
  663. CXZ(KP,L) = UMJ(KP,1,L)*UMJ(KP,3,L)*DT19
  664. CYZ(KP,L) = UMJ(KP,2,L)*UMJ(KP,3,L)*DT19
  665. 5160 CONTINUE
  666. 5170 CONTINUE
  667. C-----------------------------
  668. C Crank Nicholson généralisé :
  669. C-----------------------------
  670. ELSEIF(IDCEN.EQ.5)THEN
  671. SI1 =-1.D0
  672. SI2 = 1.D0
  673. DT19 = DTM1/6.D0
  674. DO 5190 L=1,NPG
  675. DO 5180 K=KDEB,KFIN
  676. KP = K - KDEB + 1
  677. TO1(KP,L) = DT19
  678. TO2(KP,L) = 0.D0
  679. 5180 CONTINUE
  680. 5190 CONTINUE
  681. ENDIF
  682. C
  683. C---------------------------
  684. C Pas de temps de stabilité
  685. C---------------------------
  686. C
  687. IF (IDCEN.NE.5) THEN
  688. DO 5210 L=1,NPG
  689. DO 5200 K=KDEB,KFIN
  690. KP = K - KDEB + 1
  691. DT0 = DT
  692. DT1 = 2.D0 /
  693. & ( UMJ(KP,1,L)*UMJ(KP,1,L)/(CXT(KP,L)+ANUK(KP))
  694. & + UMJ(KP,2,L)*UMJ(KP,2,L)/(CYT(KP,L)+ANUK(KP))
  695. & + UMJ(KP,3,L)*UMJ(KP,3,L)/(CZT(KP,L)+ANUK(KP)) )
  696.  
  697. DT2 = 0.5D0 /
  698. & ( (CXT(KP,L)+ANUK(KP))*AL2(KP)
  699. & + (CYT(KP,L)+ANUK(KP))*AH2(KP)
  700. & + (CZT(KP,L)+ANUK(KP))*AP2(KP) )
  701.  
  702. IF (DT1.LT.DT) DT=DT1
  703. IF (DT2.LT.DT) DT=DT2
  704. IF (DT.NE.DT0) THEN
  705. DTT1 = DT1
  706. DTT2 = DT2
  707. DIAEL = XMB(KP)
  708. NUEL = K
  709. ENDIF
  710. 5200 CONTINUE
  711. 5210 CONTINUE
  712. END IF
  713. C
  714. C---------------------------------------------------
  715. C Fonction test pour la formulation Petrov-Galerkin
  716. C---------------------------------------------------
  717. C
  718. DO 5220 L=1,NPG
  719. DO 5220 I=1,NP
  720. DO 5220 K=KDEB,KFIN
  721. KP = K - KDEB + 1
  722. WT(KP,I,L)= FN(I,L) + SI1
  723. & * ( TO1(KP,L) * ( UMJ(KP,1,L)*HK(KP,1,I,L)
  724. & + UMJ(KP,2,L)*HK(KP,2,I,L)
  725. & + UMJ(KP,3,L)*HK(KP,3,I,L) )
  726. & + TO2(KP,L) * ( UPI(KP,1,L)*HK(KP,1,I,L)
  727. & + UPI(KP,2,L)*HK(KP,2,I,L)
  728. & + UPI(KP,3,L)*HK(KP,3,I,L) ))
  729. WS(KP,I,L)= FN(I,L) + SI2
  730. & * ( TO1(KP,L) * ( UMJ(KP,1,L)*HK(KP,1,I,L)
  731. & + UMJ(KP,2,L)*HK(KP,2,I,L)
  732. & + UMJ(KP,3,L)*HK(KP,3,I,L))
  733. & + TO2(KP,L) * ( UPI(KP,1,L)*HK(KP,1,I,L)
  734. & + UPI(KP,2,L)*HK(KP,2,I,L)
  735. & + UPI(KP,3,L)*HK(KP,3,I,L) ))
  736. 5220 CONTINUE
  737. ENDIF
  738. RETURN
  739. END
  740.  
  741.  
  742.  
  743.  

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