Télécharger la plaquette Cast3M


fr en

Télécharger caldec.eso

Retour à la liste

Numérotation des lignes : 

caldec

  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