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

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