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