1 : C YCVFPU SOURCE CHAT 05/01/13 04:16:39 5004 2 : SUBROUTINE YCVFPU(NEL,K0,NPF,IES,IAXI,NPTI,IPADI,NPTS,IPADS, 3 : & LEF,XCOOR, 4 : & VOLF, 5 : & UN,TK,TE, 6 : & F, 7 : & ANU,IKC,UET,YP) 8 : IMPLICIT INTEGER(I-N) 9 : IMPLICIT REAL*8 (A-H,O-Z) 10 : 11 : C************************************************************************ 12 : C 13 : C SYNTAXE : 14 : C 15 : C FPU (NU,UET,YP <,VPAROI> ) 16 : C 17 : C 1------2 18 : C (R1,AL1) LEF FLUIDE NOEUDS 1 2 19 : C 20 : C 21 : C ANU VISCOSITE CINEMATIQUE 22 : C UET U* 23 : C YP DISTANCE A LA PAROI 24 : C VPAROI VITESSE DE LA PAROI (PAR DEFAUT 0.D0) 25 : C 26 : C CAS TRIDIMENSIONNEL 27 : C 4 ________ 3 28 : C / FLUIDE / 29 : C 1 /________/2 30 : C 31 : C 32 : C************************************************************************ 33 : -INC CCOPTIO 34 : -INC CCREEL 35 : DIMENSION UN(NPTI,IES) 36 : DIMENSION F(NPTS,IES),TK(NPTS),TE(NPTS) 37 : 38 : DIMENSION UET(*) 39 : DIMENSION AF(2,2),CF(2,2) 40 : DIMENSION ANU(*),VOLF(*) 41 : DIMENSION XCOOR(*),XYZ(3,4) 42 : DIMENSION LEF(NPF,NEL),IPADI(*),IPADS(*) 43 : DIMENSION LE(4) 44 : DIMENSION UIX(4),UIY(4),UIZ(4),UFI(4) 45 : 46 : DIMENSION UPIU(101) 47 : 48 : C -------------------CAROLI----I------------------------------------ 49 : DATA UPIU( 1) /0.000000D+00/ 50 : DATA UPIU( 2) /0.999954D+00/ 51 : DATA UPIU( 3) /0.199858D+01/ 52 : DATA UPIU( 4) /0.298974D+01/ 53 : DATA UPIU( 5) /0.395957D+01/ 54 : DATA UPIU( 6) /0.488802D+01/ 55 : DATA UPIU( 7) /0.575502D+01/ 56 : DATA UPIU( 8) /0.654720D+01/ 57 : DATA UPIU( 9) /0.726016D+01/ 58 : DATA UPIU( 10) /0.789663D+01/ 59 : DATA UPIU( 11) /0.846323D+01/ 60 : DATA UPIU( 12) /0.896802D+01/ 61 : DATA UPIU( 13) /0.941906D+01/ 62 : DATA UPIU( 14) /0.982372D+01/ 63 : DATA UPIU( 15) /0.101884D+02/ 64 : DATA UPIU( 16) /0.105188D+02/ 65 : DATA UPIU( 17) /0.108193D+02/ 66 : DATA UPIU( 18) /0.110941D+02/ 67 : DATA UPIU( 19) /0.113464D+02/ 68 : DATA UPIU( 20) /0.115790D+02/ 69 : DATA UPIU( 21) /0.117944D+02/ 70 : DATA UPIU( 22) /0.119944D+02/ 71 : DATA UPIU( 23) /0.121809D+02/ 72 : DATA UPIU( 24) /0.123553D+02/ 73 : DATA UPIU( 25) /0.125189D+02/ 74 : DATA UPIU( 26) /0.126727D+02/ 75 : DATA UPIU( 27) /0.128178D+02/ 76 : DATA UPIU( 28) /0.129549D+02/ 77 : DATA UPIU( 29) /0.130849D+02/ 78 : DATA UPIU( 30) /0.132082D+02/ 79 : DATA UPIU( 31) /0.133256D+02/ 80 : DATA UPIU( 32) /0.134375D+02/ 81 : DATA UPIU( 33) /0.135443D+02/ 82 : DATA UPIU( 34) /0.136465D+02/ 83 : DATA UPIU( 35) /0.137444D+02/ 84 : DATA UPIU( 36) /0.138383D+02/ 85 : DATA UPIU( 37) /0.139285D+02/ 86 : DATA UPIU( 38) /0.140153D+02/ 87 : DATA UPIU( 39) /0.140988D+02/ 88 : DATA UPIU( 40) /0.141794D+02/ 89 : DATA UPIU( 41) /0.142573D+02/ 90 : DATA UPIU( 42) /0.143325D+02/ 91 : DATA UPIU( 43) /0.144052D+02/ 92 : DATA UPIU( 44) /0.144757D+02/ 93 : DATA UPIU( 45) /0.145440D+02/ 94 : DATA UPIU( 46) /0.146102D+02/ 95 : DATA UPIU( 47) /0.146746D+02/ 96 : DATA UPIU( 48) /0.147371D+02/ 97 : DATA UPIU( 49) /0.147979D+02/ 98 : DATA UPIU( 50) /0.148570D+02/ 99 : DATA UPIU( 51) /0.149147D+02/ 100 : DATA UPIU( 52) /0.149708D+02/ 101 : DATA UPIU( 53) /0.150256D+02/ 102 : DATA UPIU( 54) /0.150790D+02/ 103 : DATA UPIU( 55) /0.151311D+02/ 104 : DATA UPIU( 56) /0.151821D+02/ 105 : DATA UPIU( 57) /0.152319D+02/ 106 : DATA UPIU( 58) /0.152806D+02/ 107 : DATA UPIU( 59) /0.153282D+02/ 108 : DATA UPIU( 60) /0.153749D+02/ 109 : DATA UPIU( 61) /0.154206D+02/ 110 : DATA UPIU( 62) /0.154653D+02/ 111 : DATA UPIU( 63) /0.155092D+02/ 112 : DATA UPIU( 64) /0.155522D+02/ 113 : DATA UPIU( 65) /0.155944D+02/ 114 : DATA UPIU( 66) /0.156358D+02/ 115 : DATA UPIU( 67) /0.156765D+02/ 116 : DATA UPIU( 68) /0.157164D+02/ 117 : DATA UPIU( 69) /0.157556D+02/ 118 : DATA UPIU( 70) /0.157942D+02/ 119 : DATA UPIU( 71) /0.158321D+02/ 120 : DATA UPIU( 72) /0.158694D+02/ 121 : DATA UPIU( 73) /0.159060D+02/ 122 : DATA UPIU( 74) /0.159421D+02/ 123 : DATA UPIU( 75) /0.159776D+02/ 124 : DATA UPIU( 76) /0.160126D+02/ 125 : DATA UPIU( 77) /0.160470D+02/ 126 : DATA UPIU( 78) /0.160809D+02/ 127 : DATA UPIU( 79) /0.161143D+02/ 128 : DATA UPIU( 80) /0.161472D+02/ 129 : DATA UPIU( 81) /0.161797D+02/ 130 : DATA UPIU( 82) /0.162117D+02/ 131 : DATA UPIU( 83) /0.162433D+02/ 132 : DATA UPIU( 84) /0.162744D+02/ 133 : DATA UPIU( 85) /0.163051D+02/ 134 : DATA UPIU( 86) /0.163354D+02/ 135 : DATA UPIU( 87) /0.163653D+02/ 136 : DATA UPIU( 88) /0.163949D+02/ 137 : DATA UPIU( 89) /0.164240D+02/ 138 : DATA UPIU( 90) /0.164528D+02/ 139 : DATA UPIU( 91) /0.164813D+02/ 140 : DATA UPIU( 92) /0.165094D+02/ 141 : DATA UPIU( 93) /0.165371D+02/ 142 : DATA UPIU( 94) /0.165646D+02/ 143 : DATA UPIU( 95) /0.165917D+02/ 144 : DATA UPIU( 96) /0.166185D+02/ 145 : DATA UPIU( 97) /0.166450D+02/ 146 : DATA UPIU( 98) /0.166712D+02/ 147 : DATA UPIU( 99) /0.166971D+02/ 148 : DATA UPIU(100) /0.167228D+02/ 149 : DATA UPIU(101) /0.167481D+02/ 150 : C -----------------------CAROLI---------F----------------------- 151 : 152 : C CONSTANTES PHYSIQUES 153 : C 154 : DATA CNU,AKER /0.09D0,0.41D0/ 155 : 156 : C 157 : IF(IKC.NE.1)GO TO 9990 158 : R1=1.D0 159 : C 160 : C 161 : SQCNU=SQRT(CNU) 162 : KC=1 163 : C --------------------CAROLI-----------I----------------------- 164 : ANUL=ANU(KC) 165 : C --------------------CAROLI-----------F----------------------- 166 : 167 : IF(IES.EQ.3)GO TO 300 168 : 169 : C ********* 170 : C * 2D * 171 : C ********* 172 : 173 : CF(1,1)=R1/3.D0 174 : CF(2,1)=CF(1,1)/2.D0 175 : CF(1,2)=CF(2,1) 176 : CF(2,2)=CF(1,1) 177 : 178 : C IF(IPT.EQ.100)WRITE(IOIMP,*)' YCVFPU : ANUL,R1 ',ANUL,R1 179 : 180 : C write(6,*)' K0,NEL,IKC=',K0,NEL,IKC,' npf=',npf 181 : DO 50 K=1,NEL 182 : NK=K0+K 183 : KC=1+(1-IKC)*(NK-1) 183 : CALL RSETD(AF,CF,4) 184 : C 185 : C 186 : IP=LEF(1,K) 187 : XYZ(1,1)=XCOOR((IP-1)*(IES+1) +1) 188 : XYZ(2,1)=XCOOR((IP-1)*(IES+1) +2) 189 : IP=LEF(2,K) 190 : XYZ(1,2)=XCOOR((IP-1)*(IES+1) +1) 191 : XYZ(2,2)=XCOOR((IP-1)*(IES+1) +2) 192 : 193 : C WRITE(IOIMP,*)' LE.... ',(LE(MM),MM=1,4) 194 : 195 : DO 5 I=1,NPF 196 : NFU=IPADI(LEF(I,K)) 197 : UIX(I)= UN(NFU,1) 198 : UIY(I)= UN(NFU,2) 199 : 5 CONTINUE 200 : 201 : C 202 : C CALCUL DES LONGUEURS DE L'ELEMENT : AL1 203 : C 204 : IF (IAXI.EQ.0) THEN 205 : R1=1.D0 206 : ELSEIF (IAXI.EQ.1) THEN 207 : R1=(XYZ(2,1)+XYZ(2,2))*XPI 208 : ELSEIF (IAXI.EQ.2) THEN 209 : R1=(XYZ(1,1)+XYZ(1,2))*XPI 210 : ENDIF 211 : 212 : AL1X=XYZ(1,2)-XYZ(1,1) 213 : AL1Y=XYZ(2,2)-XYZ(2,1) 214 : AL1=SQRT(AL1X*AL1X+AL1Y*AL1Y) 215 : 216 : C 217 : C CALCUL DES COMPOSANTES TANGENTIELLES A LA PAROI DES VITESSES 218 : C ET DE LA MOYENNE DE LEURS DIFFERENCES FLUIDE-PAROI : UTM 219 : C 220 : 221 : UT1=UIX(1)*UIX(1) + UIY(1)*UIY(1) 222 : UT2=UIX(2)*UIX(2) + UIY(2)*UIY(2) 223 : UT1=SQRT(UT1) 224 : UT2=SQRT(UT2) 225 : UTM=(UT1+UT2)/2.D0 + 1.D-5 226 : 227 : C 228 : C CALCUL DE U ETOILE A PARTIR UTM 229 : C 230 : C --------------------CAROLI---I-------------------------- 231 : DO 57 I=1,4 232 : YPLUS=YP*UET(K)/ANUL 233 : YPLUS=ABS(YPLUS)+1.D-5 234 : IF(YPLUS.GT.100.0D0) THEN 235 : UPLUS=5.5D0+(LOG(YPLUS))/0.41D0 236 : ELSE 237 : IY=INT(YPLUS) 238 : RY=YPLUS-IY 239 : UPLUS=UPIU(IY+1)+RY*(UPIU(IY+2)-UPIU(IY+1)) 240 : ENDIF 241 : C ---- RELAX SU UET ------------------------------------ 242 : UET(NK)=0.2D0*UET(NK)+0.8D0*(ABS(UTM/(UPLUS+1.D-5))) 243 : 57 CONTINUE 244 : C --------------------CAROLI----F------------------ 245 : 246 : C************************************************************************ 247 : C CALCUL Q D M 248 : C************************************************************************ 249 : 250 : C CALCUL DU COEFFICIENT AK A PARTIR DE UTM 251 : 252 : AUTM=ABS(UTM) 253 : AK=0.D0 254 : IF(AUTM.GT.1.D-10)AK=UET(NK)*UET(NK)/AUTM 255 : 256 : C CALCUL DE LA MATRICE VITESSE 257 : DO 82 I=1,2 258 : DO 72 J=1,2 259 : AF(J,I)=AK*AL1*AF(J,I)*R1 260 : 72 CONTINUE 261 : 82 CONTINUE 262 : C 263 : C SATURATION VITESSE 264 : C 265 : DO 54 J=1,2 266 : AF(1,J)=AF(1,J)*UT1 267 : AF(2,J)=AF(2,J)*UT2 268 : 54 CONTINUE 269 : 270 : DO 65 I=1,NPF 271 : NF=IPADS(LEF(I,K)) 272 : FF1=(AF(1,I)+AF(2,I))*AL1X/AL1 273 : FF2=(AF(1,I)+AF(2,I))*AL1Y/AL1 274 : F(NF,1)=F(NF,1)-SIGN(FF1,UIX(I)) 275 : F(NF,2)=F(NF,2)-SIGN(FF2,UIY(I)) 276 : 65 CONTINUE 277 : 278 : C------------------------------------------------------------------------ 279 : C CALCUL Q D M FIN 280 : C------------------------------------------------------------------------ 281 : 282 : 283 : C 284 : C CALCUL DE KP ET DE EPSILONP 285 : C 286 : YPLUS2=-(YPLUS+0.01D0)*(YPLUS+0.01D0)*0.0017D0 287 : DO 66 I=1,NPF 288 : NF=IPADS(LEF(I,K)) 289 : TK(NF)=UET(NK)*UET(NK)/SQCNU 290 : TE(NF)=UET(NK)*UET(NK)*UET(NK)/ 291 : & (AKER*YP*(1.D0-EXP(YPLUS2))) 292 : 293 : 66 CONTINUE 294 : 295 : 50 CONTINUE 296 : RETURN 297 : 298 : C ********* 299 : C * 3D * 300 : C ********* 301 : 300 CONTINUE 302 : 303 : DO 350 K=1,NEL 304 : NK=K0+K 305 : KC=1+(1-IKC)*(NK-1) 306 : C 307 : UTM=0.D0 308 : 309 : DO 35 I=1,NPF 310 : NFA=IPADI(LEF(I,K)) 311 : UIX(I)=UN(NFA,1)+XPETIT 312 : UIY(I)=UN(NFA,2)+XPETIT 313 : UIZ(I)=UN(NFA,3)+XPETIT 314 : UFI(I)=UIX(I)*UIX(I)+UIY(I)*UIY(I)+UIZ(I)*UIZ(I) 315 : UFI(I)=SQRT(UFI(I)) 316 : UTM=UTM+UFI(I) 317 : 35 CONTINUE 318 : 319 : UTM=UTM/NPF + 1.D-5 320 : 321 : AIRF=VOLF(NK) 322 : C 323 : C CALCUL DE U ETOILE A PARTIR UTM 324 : C 325 : C --------------------CAROLI---I-------------------------- 326 : DO 357 I=1,4 327 : YPLUS=YP*UET(NK)/ANUL 328 : YPLUS=ABS(YPLUS)+1.D-5 329 : IF(YPLUS.GT.100.0D0) THEN 330 : UPLUS=5.5D0+(LOG(YPLUS))/0.41D0 331 : ELSE 332 : IY=INT(YPLUS) 333 : RY=YPLUS-IY 334 : UPLUS=UPIU(IY+1)+RY*(UPIU(IY+2)-UPIU(IY+1)) 335 : ENDIF 336 : C ---- RELAX SU UET ------------------------------------ 337 : UET(NK)=0.2D0*UET(NK)+0.8D0*(ABS(UTM/(UPLUS+1.D-5))) 338 : 357 CONTINUE 339 : C --------------------CAROLI----F------------------ 340 : 341 : C************************************************************************ 342 : C CALCUL Q D M 343 : C************************************************************************ 344 : 345 : C CALCUL DU COEFFICIENT AK A PARTIR DE UTM 346 : 347 : AUTM=ABS(UTM) 348 : AK=0.D0 349 : IF(AUTM.GT.1.D-10)AK=UET(NK)*UET(NK)/AUTM 350 : C 351 : DO 365 I=1,NPF 352 : NFA=IPADS(LEF(I,K)) 353 : F(NFA,1)=F(NFA,1)-AK/NPF*UIX(I)/UFI(I)*AIRF 354 : F(NFA,2)=F(NFA,2)-AK/NPF*UIY(I)/UFI(I)*AIRF 355 : F(NFA,3)=F(NFA,3)-AK/NPF*UIZ(I)/UFI(I)*AIRF 356 : 365 CONTINUE 357 : C------------------------------------------------------------------------ 358 : C CALCUL Q D M FIN 359 : C------------------------------------------------------------------------ 360 : C 361 : C CALCUL DE KP ET DE EPSILONP 362 : C 363 : YPLUS2=-(YPLUS+0.01D0)*(YPLUS+0.01D0)*0.0017D0 364 : DO 366 I=1,NPF 365 : NFA=IPADS(LEF(I,K)) 366 : TK(NFA)=UET(NK)*UET(NK)/SQCNU 367 : TE(NFA)=UET(NK)*UET(NK)*UET(NK)/ 368 : & (AKER*YP*(1.D0-EXP(YPLUS2))) 369 : 366 CONTINUE 370 : 371 : 350 CONTINUE 372 : 373 : RETURN 374 : 375 : 9990 CONTINUE 376 : WRITE(IOIMP,*)' IKC ',IKC 377 : WRITE(IOIMP,*) 378 : $ ' COEF D''UN MAUVAIS TYPE : ON VEUT (SCAL OU VECT DOMA)' 378 : CALL ARRET(0) 379 : RETURN 380 : 1002 FORMAT(10(1X,1PE11.4)) 381 : 1001 FORMAT(20(1X,I5)) 382 : END 383 : 384 : 385 : 386 : 387 : 388 : 389 :
© Cast3M 2003 - All rights reserved.
Disclaimer