Télécharger caldec.eso

Retour à la liste

Numérotation des lignes :

  1. C CALDEC SOURCE MAGN 10/07/08 21:15:00 6709
  2. SUBROUTINE CALDEC(WT,WS,XYZ,GR,HR,FN,NES,IDIM,NBNN,NPG,AJT,
  3. &IDCEN,CMD,V1,V2,VELCHE,TN,NC,IKOMP,XREF,AIRE,KE)
  4. IMPLICIT INTEGER(I-N)
  5. IMPLICIT REAL*8 (A-H,O-Z)
  6. C***********************************************************************
  7. C
  8. C Ce Sp calcul les fonctions de forme Wt
  9. C----------------------------------------------------------------------
  10. C HISTORIQUE : 20/10/01 : Création
  11. C
  12. C HISTORIQUE :
  13. C
  14. C
  15. C---------------------------
  16. C Paramètres Entrée/Sortie :
  17. C---------------------------
  18. C
  19. C S/WT : Fonctions de forme Tilde (Formulation Petrov Galerkin)
  20. C S/WS : Fonctions de forme Tilde (Formulation Petrov Galerkin)
  21. C pour le cas CNG uniquement
  22. C E/ XYZ : Coordommées des noeud de l'élément
  23. C E/ GR : Gradient des fonctions de forme sur l'élément de référence
  24. C E/ HR : Gradient des fonctions de forme sur l'élément courant
  25. C E/ FN : Fonctions de forme
  26. C E/ NES : dimension espace de l'élément de référence
  27. C E/ IDIM : dimension espace calcul
  28. C E/ NBNN : nombre de noeuds de l'élément
  29. C E/ NPG : nombre de points de Gauss
  30. C E/ AJT : Jacobien Transposé
  31. C E/ IDCEN : Entier indiquant le type de décentrement souhaité
  32. C IDCEN 1-> CENTREE 2-> SUPGDC 3-> SUPG 4-> TVISQUEU 5-> CNG
  33. C E/ CMD : Coefficient multiplicateur du décentrement
  34. C Si IDCEN=4 ou =5 CMD=DT
  35. C E/ V1 V2 : Coefficient de l'équation : dans l'ordre Ro et Mu
  36. C E/ VELCHE : Champ de vitesse aux points de Gauss
  37. C E/ TN : Grandeur transportée
  38. C E/ NC : Nombre de composantes de cette grandeur
  39. C E/ IKOMP : =0 formulation non conservative =1 formulation conservative
  40. C----------------------------------------------------------------------
  41. C************************************************************************
  42. DIMENSION XYZ(IDIM,NBNN),AJT(IDIM,IDIM,NPG),XREF(NES,NBNN)
  43. DIMENSION WT(NBNN,NPG),WS(NBNN,NPG)
  44. DIMENSION FN(NBNN,NPG),GR(NES,NBNN,NPG),HR(IDIM,NBNN,NPG)
  45. DIMENSION V1(NPG),V2(NPG),VELCHE(8),TN(NBNN,NC),TT(27)
  46. LOGICAL*1 KAL
  47. DIMENSION UPIL(3),GRAD(3)
  48.  
  49. C*****************************************************************************
  50. CCALDEC
  51. c write(6,*)' DEBUT CALDEC ',IDCEN,CMD
  52. IF(IDCEN.EQ.19)THEN
  53. KAL=.TRUE.
  54. ELSE
  55. KAL=.FALSE.
  56. ENDIF
  57.  
  58. C----------
  59. C CENTREE :
  60. C----------
  61. C
  62. IF (IDCEN.EQ.0.OR.IDCEN.EQ.1) THEN
  63. DO 140 LG=1,NPG
  64. DO 140 I =1,NBNN
  65. WT(I,LG) = FN(I,LG)
  66. WS(I,LG) = FN(I,LG)
  67. 140 CONTINUE
  68. RETURN
  69. ENDIF
  70.  
  71. IF(IDCEN.EQ.2)THEN
  72. IF(NBNN.GT.27)CALL ARRET(0)
  73. IF(NC.EQ.1)THEN
  74. CALL RSETD(TT,TN,NBNN)
  75. ELSE
  76. DO 141 I=1,NBNN
  77. U=0.D0
  78. DO 142 N=1,NC
  79. U=U+TN(I,N)*TN(I,N)
  80. 142 CONTINUE
  81. TT(I)=SQRT(U)
  82. 141 CONTINUE
  83. ENDIF
  84. ENDIF
  85. C
  86. C- Calcul pour chaque point de Gauss de chaque élément de :
  87. C- /L UML : Norme de la vitesse aux points de Gauss
  88. C- 1/BM: module d'un temps caractéristique associé à la convection
  89. C- XMB : caractéristique géométrique 1/2 He
  90. C
  91. XPETI=1.D-30
  92.  
  93. IF(KAL.EQV..FALSE.)THEN
  94. CALL CALJTR(GR,XYZ,NES,IDIM,NBNN,NPG,AJT)
  95. ENDIF
  96.  
  97. DO 144 LG=1,NPG
  98. ANUK=V2(LG)/V1(LG)
  99. C
  100. UL=0.D0
  101. DO 309 N=1,IDIM
  102. UL=UL+VELCHE(LG+(N-1)*NPG)*VELCHE(LG+(N-1)*NPG)
  103. 309 CONTINUE
  104. UML=SQRT(UL)+XPETI
  105. c.......................................................................
  106. IF(KAL.EQV..FALSE.)THEN
  107. BMI=0.D0
  108. DO 310 N=1,IDIM
  109. UHAT=0.D0
  110. DO 311 M=1,IDIM
  111. UHAT=UHAT+AJT(M,N,LG)*VELCHE(LG+(M-1)*NPG)
  112. 311 CONTINUE
  113. BMI=BMI+UHAT*UHAT
  114. 310 CONTINUE
  115. BM=SQRT(BMI) + XPETI
  116. XMB=UML/BM
  117. ELSE
  118. c.......................................................................
  119. XMB=0.D0
  120. BM =0.D0
  121. DO 320 N=1,IDIM
  122. XH=0.D0
  123. XB=0.D0
  124. DO 322 M=1,IDIM
  125. DO 321 K=1,NBNN
  126. ci XB=XB+VELCHE(LG+(M-1)*NPG) *GR(M,K,LG)*XYZ(N,K)
  127. ci XH=XH+VELCHE(LG+(M-1)*NPG)/UML*GR(M,K,LG)*XYZ(N,K)
  128. XB=XB+VELCHE(LG+(M-1)*NPG) *HR(M,K,LG)*XREF(N,K)
  129. ci XH=XH+VELCHE(LG+(M-1)*NPG)/UML*HR(M,K,LG)*XREF(N,K)
  130. 321 CONTINUE
  131. 322 CONTINUE
  132. ci XMB=XMB+XH*XH
  133. BM=BM+XB*XB
  134. 320 CONTINUE
  135. ci XMB=SQRT(XMB) + XPETI
  136. ci BM=UML/XMB
  137. BM=SQRT(BM) + XPETI
  138. XMB = UML/BM
  139. ENDIF
  140. c.......................................................................
  141. C
  142. CALL INITD(UPIL,3,0.D0)
  143. C
  144. C- Calcul en chaque élément, pour chaque point de Gauss de
  145. C- GRAD : vecteur unitaire porté par le gradient du champ scalaire
  146. C- UP : projection de la vitesse sur la direction donnée par GRAD
  147. C- UPIL : vecteur UP*GRAD aux points de Gauss
  148. C UIL(KP,M,L) -> VELCHE(LG+(M-1)*NPG,K)
  149. C- pour l'option SUPGDC
  150. C
  151. IF (IDCEN.EQ.2) THEN
  152. DO 170 N=1,IDIM
  153. GRAD(N)=0.D0
  154. DO 170 I=1,NBNN
  155. GRAD(N) = GRAD(N) + TT(I)*HR(N,I,LG)
  156. 170 CONTINUE
  157.  
  158. AX=0.D0
  159. DO 2301 M=1,IDIM
  160. AX = AX + GRAD(M)*GRAD(M)
  161. 2301 CONTINUE
  162. AX = SQRT(AX) + XPETI
  163.  
  164. UPL=0.D0
  165. DO 2302 N=1,IDIM
  166. GRAD(N) = GRAD(N) / AX
  167. UPL = UPL + GRAD(N) * VELCHE(LG+(N-1)*NPG)
  168. 2302 CONTINUE
  169.  
  170. DO 2303 N=1,IDIM
  171. UPIL(N) = GRAD(N) * UPL
  172. 2303 CONTINUE
  173.  
  174. C
  175. BPI=0.D0
  176. DO 410 N=1,IDIM
  177. UPHAT=0.D0
  178. DO 411 M=1,IDIM
  179. UPHAT=UPHAT+AJT(M,N,LG)*UPIL(M)
  180. 411 CONTINUE
  181.  
  182. BPI=BPI+UPHAT*UPHAT
  183. 410 CONTINUE
  184. BP=SQRT(BPI) + XPETI
  185. XPB=UPL/BP
  186. ENDIF
  187. C
  188. C-----------------------------
  189. C- DECENTREMENT suivant IDCEN
  190. C-----------------------------
  191. C On calcule dans chaque cas TO1 et TO2 ainsi que le tenseur
  192. C associé à la viscosité numérique afin d'évaluer le pas de
  193. C temps de stabilité pour les schémas explicites.
  194. C
  195. C---------
  196. C SUPGDC : Base théorique dans : A New FE formulation for computational
  197. C fluid dynamics, II Beyond SUPG, HUGHES et al., in Comp.Meth.Appl.Mech.
  198. C Eng., vol 54, pp 341-355 (1986).
  199. C---------
  200. IF (IDCEN.EQ.2) THEN
  201. C
  202. C- Approximation "doublement asymptotique" basée sur la vitesse moyenne
  203. C- HMK : Distance basé sur la vitesse moyenne
  204. C- ALFA : Peclet de maille basé sur la vitesse moyenne
  205. C
  206. ALFA = UML * XMB / (ANUK+XPETI) / 3.D0
  207. AKSI = MIN(1.D0,ALFA)
  208. CCT = AKSI / BM * CMD
  209. C
  210. C- Approximation "doublement asymptotique" basée sur la projection de
  211. C- la vitesse sur le gradient du champ scalaire
  212. C- HMK : Distance basée sur la vitesse projetée
  213. C- ALFA : Peclet de maille basé sur la vitesse projetée
  214. C
  215. ALFA = UPL * XPB / (ANUK+XPETI) / 3.D0
  216. AKSI = MIN(1.D0,ALFA)
  217. CCP = AKSI / BP
  218. C
  219. CPT = CCP - CCT
  220. CC2 = MAX(0.D0,CPT) * CMD
  221. C
  222. TO1 = CCT
  223. TO2 = CC2
  224. SI1 = 1.D0
  225. SI2 = 1.D0
  226. IF(IKOMP.EQ.1)THEN
  227. SI1 = -1.D0
  228. SI2 = -1.D0
  229. ENDIF
  230. C-------
  231. C SUPG :
  232. C-------
  233. ELSEIF (IDCEN.EQ.3.OR.IDCEN.EQ.19) THEN
  234. C
  235. C- Approximation "doublement asymptotique" basée sur la vitesse moyenne
  236. C- HMK : Distance basé sur la vitesse moyenne
  237. C- ALFA : Peclet de maille basé sur la vitesse moyenne
  238. C
  239. ALFA = UML * XMB / (ANUK+XPETI) / 3.D0
  240. AKSI = MIN(1.D0,ALFA)
  241. CCT = AKSI / BM * CMD
  242. DYY=(Aire/XMB)**(IDIM-1)
  243. Alg2=XMB/DYY
  244. C cct2=cct*(alg2**1.5)
  245. C cct2=cct*1.5
  246. c cct2=cct
  247.  
  248. c if(ke.eq.1726)then
  249. c if(ke.eq.7)then
  250. c DX=XYZ(1,3)-XYZ(1,2)
  251. c DY=XYZ(2,2)-XYZ(2,1)
  252. c air=DX*DY
  253. c Alg=DX/(aire**0.5)
  254. c cct2=cct*(alg2**0.5)
  255. c write(6,*)'---------------------------------------'
  256. c write(6,*)' DXxDY=',air,' Aire=',Aire
  257. c write(6,*)'DX =',DX,' DY =',DY,' Alg=',Alg,' Alg2=',alg2
  258. c write(6,*)'U1x=',VELCHE(1),' U1y=',VELCHE(1+NPG)
  259. cc write(6,1002)(XYZ(1,kk),kk=1,nbnn)
  260. cc write(6,*)' coor y '
  261. cc write(6,1002)(XYZ(2,kk),kk=1,nbnn)
  262. c write(6,*)'UML= XMB= BM= CCT= CCT*U u*dx*cmd'
  263. c write(6,1002)uml,xmb,bm,cct,(cct*uml*uml),(uml*xmb*cmd),
  264. c &(cct2*uml*uml)
  265. c write(6,*)'---------------------------------------'
  266. c write(6,*)'---------------------------------------'
  267. c endif
  268. c cct=cct2
  269. C
  270. TO1 = CCT
  271. TO2 = 0.D0
  272. SI1 = 1.D0
  273. SI2 = 1.D0
  274. IF(IKOMP.EQ.1)THEN
  275. SI1 = -1.D0
  276. SI2 = -1.D0
  277. ENDIF
  278. C-------------------
  279. C Tenseur Visqueux :
  280. C-------------------
  281. ELSEIF (IDCEN.EQ.4) THEN
  282. DT19 = CMD * 0.5D0
  283. TO1 = DT19
  284. TO2 = 0.D0
  285. SI1 = 1.D0
  286. SI2 = 1.D0
  287. IF(IKOMP.EQ.1)THEN
  288. SI1 = -1.D0
  289. SI2 = -1.D0
  290. ENDIF
  291. C-----------------------------
  292. C Crank Nicholson généralisé :
  293. C-----------------------------
  294. ELSEIF (IDCEN.EQ.5) THEN
  295. DT19 = CMD /3.D0
  296. TO1 = DT19
  297. TO2 = 0.D0
  298. SI1 = -1.D0
  299. SI2 = 1.D0
  300. IF(IKOMP.EQ.1)THEN
  301. SI1 = 1.D0
  302. SI2 = -1.D0
  303. ENDIF
  304. ENDIF
  305. C
  306. C
  307. C---------------------------------------------------
  308. C Fonction test pour la formulation Petrov-Galerkin
  309. C---------------------------------------------------
  310. C Ce qui est diffusion numérique en explicite se transforme en
  311. C ajoutant de la viscosité numérique (Tenseurs visqueux et CNG).
  312. C WS : Fonction test pour la partie explicite
  313. C
  314. cw=0.
  315. IF(IDIM.EQ.2)THEN
  316. U1=VELCHE(LG)
  317. SU1=SIGN(U1,UML)*cw
  318. U2=VELCHE(LG+NPG)
  319. SU2=SIGN(U2,UML)*cw
  320. TU1=TO1*(U1+SU1)+TO2*UPIL(1)
  321. TU2=TO1*(U2+SU2)+TO2*UPIL(2)
  322. DO 2050 I=1,NBNN
  323. W=TU1*HR(1,I,LG)+TU2*HR(2,I,LG)
  324. WT(I,LG) = FN(I,LG) + SI1*W
  325. WS(I,LG) = FN(I,LG) + SI2*W
  326. 2050 CONTINUE
  327. c write(6,*)'TO1 TO2 = ',TO1,TO2
  328. C
  329. ELSEIF(IDIM.EQ.3)THEN
  330. TU1=TO1*VELCHE(LG)+TO2*UPIL(1)
  331. TU2=TO1*VELCHE(LG+NPG)+TO2*UPIL(2)
  332. TU3=TO1*VELCHE(LG+2*NPG)+TO2*UPIL(3)
  333. DO 2051 I=1,NBNN
  334. W=TU1*HR(1,I,LG)+TU2*HR(2,I,LG)+TU3*HR(3,I,LG)
  335. WT(I,LG) = FN(I,LG) + SI1*W
  336. WS(I,LG) = FN(I,LG) + SI2*W
  337. 2051 CONTINUE
  338. ENDIF
  339. C
  340. C- Si on est en conservatif, on rétablit les valeurs de la vitesse
  341. C de transport qui ont été modifiées car elles sont utilisées
  342. C dans d'autres subroutines ultérieures.
  343. C
  344. c IF (IKOMP .EQ. 1) THEN
  345. c DO 235 N=1,IDIM
  346. c UIL(KP,N,L) = XPETI
  347. c DO 215 I=1,NBNN
  348. c NF = IPADU(LE(I,K))
  349. c UIL(KP,N,L) = UIL(KP,N,L) + UN(NF,N)*FN(I,L)
  350. c215 CONTINUE
  351. c235 CONTINUE
  352. c ENDIF
  353. c
  354. 144 CONTINUE
  355.  
  356. C*************************************************************************
  357. c write(6,*)' FIN CALDEC '
  358. RETURN
  359. 1001 FORMAT(20(1X,I5))
  360. 1002 FORMAT(10(1X,1PE11.4))
  361. END
  362.  
  363.  

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