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
  • WT(I,LG) = FN(I,LG)
  • WS(I,LG) = FN(I,LG)
  • 140 CONTINUE
  • RETURN
  • ENDIF
  •  
  • IF(IDCEN.EQ.2)THEN
  • IF(NBNN.GT.27)CALL ARRET(0)
  • IF(NC.EQ.1)THEN
  • CALL RSETD(TT,TN,NBNN)
  • ELSE
  • DO 141 I=1,NBNN
  • U=0.D0
  • DO 142 N=1,NC
  • U=U+TN(I,N)*TN(I,N)
  • 142 CONTINUE
  • TT(I)=SQRT(U)
  • 141 CONTINUE
  • ENDIF
  • ENDIF
  • C
  • C- Calcul pour chaque point de Gauss de chaque élément de :
  • C- /L UML : Norme de la vitesse aux points de Gauss
  • C- 1/BM: module d'un temps caractéristique associé à la convection
  • C- XMB : caractéristique géométrique 1/2 He
  • C
  • XPETI=1.D-30
  •  
  • IF(KAL.EQV..FALSE.)THEN
  • CALL CALJTR(GR,XYZ,NES,IDIM,NBNN,NPG,AJT)
  • ENDIF
  •  
  • DO 144 LG=1,NPG
  • ANUK=V2(LG)/V1(LG)
  • C
  • UL=0.D0
  • DO 309 N=1,IDIM
  • UL=UL+VELCHE(LG+(N-1)*NPG)*VELCHE(LG+(N-1)*NPG)
  • 309 CONTINUE
  • UML=SQRT(UL)+XPETI
  • c.......................................................................
  • IF(KAL.EQV..FALSE.)THEN
  • BMI=0.D0
  • DO 310 N=1,IDIM
  • UHAT=0.D0
  • DO 311 M=1,IDIM
  • UHAT=UHAT+AJT(M,N,LG)*VELCHE(LG+(M-1)*NPG)
  • 311 CONTINUE
  • BMI=BMI+UHAT*UHAT
  • 310 CONTINUE
  • BM=SQRT(BMI) + XPETI
  • XMB=UML/BM
  • ELSE
  • c.......................................................................
  • XMB=0.D0
  • BM =0.D0
  • DO 320 N=1,IDIM
  • XH=0.D0
  • XB=0.D0
  • DO 322 M=1,IDIM
  • DO 321 K=1,NBNN
  • ci XB=XB+VELCHE(LG+(M-1)*NPG) *GR(M,K,LG)*XYZ(N,K)
  • ci XH=XH+VELCHE(LG+(M-1)*NPG)/UML*GR(M,K,LG)*XYZ(N,K)
  • XB=XB+VELCHE(LG+(M-1)*NPG) *HR(M,K,LG)*XREF(N,K)
  • ci XH=XH+VELCHE(LG+(M-1)*NPG)/UML*HR(M,K,LG)*XREF(N,K)
  • 321 CONTINUE
  • 322 CONTINUE
  • ci XMB=XMB+XH*XH
  • BM=BM+XB*XB
  • 320 CONTINUE
  • ci XMB=SQRT(XMB) + XPETI
  • ci BM=UML/XMB
  • BM=SQRT(BM) + XPETI
  • XMB = UML/BM
  • ENDIF
  • c.......................................................................
  • C
  • CALL INITD(UPIL,3,0.D0)
  • C
  • C- Calcul en chaque élément, pour chaque point de Gauss de
  • C- GRAD : vecteur unitaire porté par le gradient du champ scalaire
  • C- UP : projection de la vitesse sur la direction donnée par GRAD
  • C- UPIL : vecteur UP*GRAD aux points de Gauss
  • C UIL(KP,M,L) -> VELCHE(LG+(M-1)*NPG,K)
  • C- pour l'option SUPGDC
  • C
  • IF (IDCEN.EQ.2) THEN
  • DO 170 N=1,IDIM
  • GRAD(N)=0.D0
  • DO 170 I=1,NBNN
  • GRAD(N) = GRAD(N) + TT(I)*HR(N,I,LG)
  • 170 CONTINUE
  •  
  • AX=0.D0
  • DO 2301 M=1,IDIM
  • AX = AX + GRAD(M)*GRAD(M)
  • 2301 CONTINUE
  • AX = SQRT(AX) + XPETI
  •  
  • UPL=0.D0
  • DO 2302 N=1,IDIM
  • GRAD(N) = GRAD(N) / AX
  • UPL = UPL + GRAD(N) * VELCHE(LG+(N-1)*NPG)
  • 2302 CONTINUE
  •  
  • DO 2303 N=1,IDIM
  • UPIL(N) = GRAD(N) * UPL
  • 2303 CONTINUE
  •  
  • C
  • BPI=0.D0
  • DO 410 N=1,IDIM
  • UPHAT=0.D0
  • DO 411 M=1,IDIM
  • UPHAT=UPHAT+AJT(M,N,LG)*UPIL(M)
  • 411 CONTINUE
  •  
  • BPI=BPI+UPHAT*UPHAT
  • 410 CONTINUE
  • BP=SQRT(BPI) + XPETI
  • XPB=UPL/BP
  • ENDIF
  • C
  • C-----------------------------
  • C- DECENTREMENT suivant IDCEN
  • C-----------------------------
  • C On calcule dans chaque cas TO1 et TO2 ainsi que le tenseur
  • C associé à la viscosité numérique afin d'évaluer le pas de
  • C temps de stabilité pour les schémas explicites.
  • C
  • C---------
  • C SUPGDC : Base théorique dans : A New FE formulation for computational
  • C fluid dynamics, II Beyond SUPG, HUGHES et al., in Comp.Meth.Appl.Mech.
  • C Eng., vol 54, pp 341-355 (1986).
  • C---------
  • IF (IDCEN.EQ.2) THEN
  • C
  • C- Approximation "doublement asymptotique" basée sur la vitesse moyenne
  • C- HMK : Distance basé sur la vitesse moyenne
  • C- ALFA : Peclet de maille basé sur la vitesse moyenne
  • C
  • ALFA = UML * XMB / (ANUK+XPETI) / 3.D0
  • AKSI = MIN(1.D0,ALFA)
  • CCT = AKSI / BM * CMD
  • C
  • C- Approximation "doublement asymptotique" basée sur la projection de
  • C- la vitesse sur le gradient du champ scalaire
  • C- HMK : Distance basée sur la vitesse projetée
  • C- ALFA : Peclet de maille basé sur la vitesse projetée
  • C
  • ALFA = UPL * XPB / (ANUK+XPETI) / 3.D0
  • AKSI = MIN(1.D0,ALFA)
  • CCP = AKSI / BP
  • C
  • CPT = CCP - CCT
  • CC2 = MAX(0.D0,CPT) * CMD
  • C
  • TO1 = CCT
  • TO2 = CC2
  • SI1 = 1.D0
  • SI2 = 1.D0
  • IF(IKOMP.EQ.1)THEN
  • SI1 = -1.D0
  • SI2 = -1.D0
  • ENDIF
  • C-------
  • C SUPG :
  • C-------
  • ELSEIF (IDCEN.EQ.3.OR.IDCEN.EQ.19) THEN
  • C
  • C- Approximation "doublement asymptotique" basée sur la vitesse moyenne
  • C- HMK : Distance basé sur la vitesse moyenne
  • C- ALFA : Peclet de maille basé sur la vitesse moyenne
  • C
  • ALFA = UML * XMB / (ANUK+XPETI) / 3.D0
  • AKSI = MIN(1.D0,ALFA)
  • CCT = AKSI / BM * CMD
  • DYY=(Aire/XMB)**(IDIM-1)
  • Alg2=XMB/DYY
  • C cct2=cct*(alg2**1.5)
  • C cct2=cct*1.5
  • c cct2=cct
  •  
  • c if(ke.eq.1726)then
  • c if(ke.eq.7)then
  • c DX=XYZ(1,3)-XYZ(1,2)
  • c DY=XYZ(2,2)-XYZ(2,1)
  • c air=DX*DY
  • c Alg=DX/(aire**0.5)
  • c cct2=cct*(alg2**0.5)
  • c write(6,*)'---------------------------------------'
  • c write(6,*)' DXxDY=',air,' Aire=',Aire
  • c write(6,*)'DX =',DX,' DY =',DY,' Alg=',Alg,' Alg2=',alg2
  • c write(6,*)'U1x=',VELCHE(1),' U1y=',VELCHE(1+NPG)
  • cc write(6,1002)(XYZ(1,kk),kk=1,nbnn)
  • cc write(6,*)' coor y '
  • cc write(6,1002)(XYZ(2,kk),kk=1,nbnn)
  • c write(6,*)'UML= XMB= BM= CCT= CCT*U u*dx*cmd'
  • c write(6,1002)uml,xmb,bm,cct,(cct*uml*uml),(uml*xmb*cmd),
  • c &(cct2*uml*uml)
  • c write(6,*)'---------------------------------------'
  • c write(6,*)'---------------------------------------'
  • c endif
  • c cct=cct2
  • C
  • TO1 = CCT
  • TO2 = 0.D0
  • SI1 = 1.D0
  • SI2 = 1.D0
  • IF(IKOMP.EQ.1)THEN
  • SI1 = -1.D0
  • SI2 = -1.D0
  • ENDIF
  • C-------------------
  • C Tenseur Visqueux :
  • C-------------------
  • ELSEIF (IDCEN.EQ.4) THEN
  • DT19 = CMD * 0.5D0
  • TO1 = DT19
  • TO2 = 0.D0
  • SI1 = 1.D0
  • SI2 = 1.D0
  • IF(IKOMP.EQ.1)THEN
  • SI1 = -1.D0
  • SI2 = -1.D0
  • ENDIF
  • C-----------------------------
  • C Crank Nicholson généralisé :
  • C-----------------------------
  • ELSEIF (IDCEN.EQ.5) THEN
  • DT19 = CMD /3.D0
  • TO1 = DT19
  • TO2 = 0.D0
  • SI1 = -1.D0
  • SI2 = 1.D0
  • IF(IKOMP.EQ.1)THEN
  • SI1 = 1.D0
  • SI2 = -1.D0
  • ENDIF
  • ENDIF
  • C
  • C
  • C---------------------------------------------------
  • C Fonction test pour la formulation Petrov-Galerkin
  • C---------------------------------------------------
  • C Ce qui est diffusion numérique en explicite se transforme en
  • C ajoutant de la viscosité numérique (Tenseurs visqueux et CNG).
  • C WS : Fonction test pour la partie explicite
  • C
  • cw=0.
  • IF(IDIM.EQ.2)THEN
  • U1=VELCHE(LG)
  • SU1=SIGN(U1,UML)*cw
  • U2=VELCHE(LG+NPG)
  • SU2=SIGN(U2,UML)*cw
  • TU1=TO1*(U1+SU1)+TO2*UPIL(1)
  • TU2=TO1*(U2+SU2)+TO2*UPIL(2)
  • DO 2050 I=1,NBNN
  • W=TU1*HR(1,I,LG)+TU2*HR(2,I,LG)
  • WT(I,LG) = FN(I,LG) + SI1*W
  • WS(I,LG) = FN(I,LG) + SI2*W
  • 2050 CONTINUE
  • c write(6,*)'TO1 TO2 = ',TO1,TO2
  • C
  • ELSEIF(IDIM.EQ.3)THEN
  • TU1=TO1*VELCHE(LG)+TO2*UPIL(1)
  • TU2=TO1*VELCHE(LG+NPG)+TO2*UPIL(2)
  • TU3=TO1*VELCHE(LG+2*NPG)+TO2*UPIL(3)
  • DO 2051 I=1,NBNN
  • W=TU1*HR(1,I,LG)+TU2*HR(2,I,LG)+TU3*HR(3,I,LG)
  • WT(I,LG) = FN(I,LG) + SI1*W
  • WS(I,LG) = FN(I,LG) + SI2*W
  • 2051 CONTINUE
  • ENDIF
  • C
  • C- Si on est en conservatif, on rétablit les valeurs de la vitesse
  • C de transport qui ont été modifiées car elles sont utilisées
  • C dans d'autres subroutines ultérieures.
  • C
  • c IF (IKOMP .EQ. 1) THEN
  • c DO 235 N=1,IDIM
  • c UIL(KP,N,L) = XPETI
  • c DO 215 I=1,NBNN
  • c NF = IPADU(LE(I,K))
  • c UIL(KP,N,L) = UIL(KP,N,L) + UN(NF,N)*FN(I,L)
  • c215 CONTINUE
  • c235 CONTINUE
  • c ENDIF
  • c
  • 144 CONTINUE
  •  
  • C*************************************************************************
  • c write(6,*)' FIN CALDEC '
  • RETURN
  • 1001 FORMAT(20(1X,I5))
  • 1002 FORMAT(10(1X,1PE11.4))
  • END
  •  
  •  
  • © Cast3M 2003 - Tous droits réservés.
    Mentions légales