C CALDEC SOURCE MAGN 10/07/08 21:15:00 6709 SUBROUTINE CALDEC(WT,WS,XYZ,GR,HR,FN,NES,IDIM,NBNN,NPG,AJT, &IDCEN,CMD,V1,V2,VELCHE,TN,NC,IKOMP,XREF,AIRE,KE) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C*********************************************************************** C C Ce Sp calcul les fonctions de forme Wt C---------------------------------------------------------------------- C HISTORIQUE : 20/10/01 : Création C C HISTORIQUE : C C C--------------------------- C Paramètres Entrée/Sortie : C--------------------------- C C S/WT : Fonctions de forme Tilde (Formulation Petrov Galerkin) C S/WS : Fonctions de forme Tilde (Formulation Petrov Galerkin) C pour le cas CNG uniquement C E/ XYZ : Coordommées des noeud de l'élément C E/ GR : Gradient des fonctions de forme sur l'élément de référence C E/ HR : Gradient des fonctions de forme sur l'élément courant C E/ FN : Fonctions de forme C E/ NES : dimension espace de l'élément de référence C E/ IDIM : dimension espace calcul C E/ NBNN : nombre de noeuds de l'élément C E/ NPG : nombre de points de Gauss C E/ AJT : Jacobien Transposé C E/ IDCEN : Entier indiquant le type de décentrement souhaité C IDCEN 1-> CENTREE 2-> SUPGDC 3-> SUPG 4-> TVISQUEU 5-> CNG C E/ CMD : Coefficient multiplicateur du décentrement C Si IDCEN=4 ou =5 CMD=DT C E/ V1 V2 : Coefficient de l'équation : dans l'ordre Ro et Mu C E/ VELCHE : Champ de vitesse aux points de Gauss C E/ TN : Grandeur transportée C E/ NC : Nombre de composantes de cette grandeur C E/ IKOMP : =0 formulation non conservative =1 formulation conservative C---------------------------------------------------------------------- C************************************************************************ DIMENSION XYZ(IDIM,NBNN),AJT(IDIM,IDIM,NPG),XREF(NES,NBNN) DIMENSION WT(NBNN,NPG),WS(NBNN,NPG) DIMENSION FN(NBNN,NPG),GR(NES,NBNN,NPG),HR(IDIM,NBNN,NPG) DIMENSION V1(NPG),V2(NPG),VELCHE(8),TN(NBNN,NC),TT(27) LOGICAL*1 KAL DIMENSION UPIL(3),GRAD(3) C***************************************************************************** CCALDEC c write(6,*)' DEBUT CALDEC ',IDCEN,CMD IF(IDCEN.EQ.19)THEN KAL=.TRUE. ELSE KAL=.FALSE. ENDIF C---------- C CENTREE : C---------- C IF (IDCEN.EQ.0.OR.IDCEN.EQ.1) THEN DO 140 LG=1,NPG 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