C ZCVFPU SOURCE CHAT 05/01/13 04:21:45 5004 SUBROUTINE ZCVFPU(NEL,K0,NPF,IES,NPTD,IAXI,IPADL, & LEF,XCOOR, & VOLF, & UN,TK,TE, & F, & DK,DE, & ANU,IKC,UET,YP) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C************************************************************************ C C SYNTAXE : C C FPU (NU,UET,YP <,VPAROI> ) C C 1------2 C (R1,AL1) LEF FLUIDE NOEUDS 1 2 C C C ANU VISCOSITE CINEMATIQUE C UET U* C YP DISTANCE A LA PAROI C VPAROI VITESSE DE LA PAROI (PAR DEFAUT 0.D0) C C CAS TRIDIMENSIONNEL C 4 ________ 3 C / FLUIDE / C 1 /________/2 C C C************************************************************************ -INC PPARAM -INC CCOPTIO -INC CCREEL DIMENSION UN(NPTD,IES),TK(*),TE(*) REAL*8 F(NPTD,IES) REAL*8 DK(*),DE(*) DIMENSION UET(*) DIMENSION AF(2,2),CF(2,2) DIMENSION ANU(*),VOLF(*) DIMENSION XCOOR(*),XYZ(3,4) DIMENSION LEF(NPF,NEL),IPADL(*) DIMENSION LE(4) DIMENSION UIX(4),UIY(4),UIZ(4),UFI(4) DIMENSION UPIU(101) C -------------------CAROLI----I------------------------------------ DATA UPIU( 1) /0.000000D+00/ DATA UPIU( 2) /0.999954D+00/ DATA UPIU( 3) /0.199858D+01/ DATA UPIU( 4) /0.298974D+01/ DATA UPIU( 5) /0.395957D+01/ DATA UPIU( 6) /0.488802D+01/ DATA UPIU( 7) /0.575502D+01/ DATA UPIU( 8) /0.654720D+01/ DATA UPIU( 9) /0.726016D+01/ DATA UPIU( 10) /0.789663D+01/ DATA UPIU( 11) /0.846323D+01/ DATA UPIU( 12) /0.896802D+01/ DATA UPIU( 13) /0.941906D+01/ DATA UPIU( 14) /0.982372D+01/ DATA UPIU( 15) /0.101884D+02/ DATA UPIU( 16) /0.105188D+02/ DATA UPIU( 17) /0.108193D+02/ DATA UPIU( 18) /0.110941D+02/ DATA UPIU( 19) /0.113464D+02/ DATA UPIU( 20) /0.115790D+02/ DATA UPIU( 21) /0.117944D+02/ DATA UPIU( 22) /0.119944D+02/ DATA UPIU( 23) /0.121809D+02/ DATA UPIU( 24) /0.123553D+02/ DATA UPIU( 25) /0.125189D+02/ DATA UPIU( 26) /0.126727D+02/ DATA UPIU( 27) /0.128178D+02/ DATA UPIU( 28) /0.129549D+02/ DATA UPIU( 29) /0.130849D+02/ DATA UPIU( 30) /0.132082D+02/ DATA UPIU( 31) /0.133256D+02/ DATA UPIU( 32) /0.134375D+02/ DATA UPIU( 33) /0.135443D+02/ DATA UPIU( 34) /0.136465D+02/ DATA UPIU( 35) /0.137444D+02/ DATA UPIU( 36) /0.138383D+02/ DATA UPIU( 37) /0.139285D+02/ DATA UPIU( 38) /0.140153D+02/ DATA UPIU( 39) /0.140988D+02/ DATA UPIU( 40) /0.141794D+02/ DATA UPIU( 41) /0.142573D+02/ DATA UPIU( 42) /0.143325D+02/ DATA UPIU( 43) /0.144052D+02/ DATA UPIU( 44) /0.144757D+02/ DATA UPIU( 45) /0.145440D+02/ DATA UPIU( 46) /0.146102D+02/ DATA UPIU( 47) /0.146746D+02/ DATA UPIU( 48) /0.147371D+02/ DATA UPIU( 49) /0.147979D+02/ DATA UPIU( 50) /0.148570D+02/ DATA UPIU( 51) /0.149147D+02/ DATA UPIU( 52) /0.149708D+02/ DATA UPIU( 53) /0.150256D+02/ DATA UPIU( 54) /0.150790D+02/ DATA UPIU( 55) /0.151311D+02/ DATA UPIU( 56) /0.151821D+02/ DATA UPIU( 57) /0.152319D+02/ DATA UPIU( 58) /0.152806D+02/ DATA UPIU( 59) /0.153282D+02/ DATA UPIU( 60) /0.153749D+02/ DATA UPIU( 61) /0.154206D+02/ DATA UPIU( 62) /0.154653D+02/ DATA UPIU( 63) /0.155092D+02/ DATA UPIU( 64) /0.155522D+02/ DATA UPIU( 65) /0.155944D+02/ DATA UPIU( 66) /0.156358D+02/ DATA UPIU( 67) /0.156765D+02/ DATA UPIU( 68) /0.157164D+02/ DATA UPIU( 69) /0.157556D+02/ DATA UPIU( 70) /0.157942D+02/ DATA UPIU( 71) /0.158321D+02/ DATA UPIU( 72) /0.158694D+02/ DATA UPIU( 73) /0.159060D+02/ DATA UPIU( 74) /0.159421D+02/ DATA UPIU( 75) /0.159776D+02/ DATA UPIU( 76) /0.160126D+02/ DATA UPIU( 77) /0.160470D+02/ DATA UPIU( 78) /0.160809D+02/ DATA UPIU( 79) /0.161143D+02/ DATA UPIU( 80) /0.161472D+02/ DATA UPIU( 81) /0.161797D+02/ DATA UPIU( 82) /0.162117D+02/ DATA UPIU( 83) /0.162433D+02/ DATA UPIU( 84) /0.162744D+02/ DATA UPIU( 85) /0.163051D+02/ DATA UPIU( 86) /0.163354D+02/ DATA UPIU( 87) /0.163653D+02/ DATA UPIU( 88) /0.163949D+02/ DATA UPIU( 89) /0.164240D+02/ DATA UPIU( 90) /0.164528D+02/ DATA UPIU( 91) /0.164813D+02/ DATA UPIU( 92) /0.165094D+02/ DATA UPIU( 93) /0.165371D+02/ DATA UPIU( 94) /0.165646D+02/ DATA UPIU( 95) /0.165917D+02/ DATA UPIU( 96) /0.166185D+02/ DATA UPIU( 97) /0.166450D+02/ DATA UPIU( 98) /0.166712D+02/ DATA UPIU( 99) /0.166971D+02/ DATA UPIU(100) /0.167228D+02/ DATA UPIU(101) /0.167481D+02/ C -----------------------CAROLI---------F----------------------- C CONSTANTES PHYSIQUES C DATA CNU,AKER /0.09D0,0.41D0/ C IF(IKC.NE.1)GO TO 9990 R1=1.D0 C C SQCNU=SQRT(CNU) KC=1 C --------------------CAROLI-----------I----------------------- ANUL=ANU(KC) C --------------------CAROLI-----------F----------------------- IF(IES.EQ.3)GO TO 300 C ********* C * 2D * C ********* CF(1,1)=R1/3.D0 CF(2,1)=CF(1,1)/2.D0 CF(1,2)=CF(2,1) CF(2,2)=CF(1,1) C IF(IPT.EQ.100)WRITE(IOIMP,*)' YCVFPU : ANUL,R1 ',ANUL,R1 DO 50 K=1,NEL NK=K0+K C WRITE(IOIMP,*)' XCVFPS NK=',NK KC=1+(1-IKC)*(NK-1) CALL RSETD(AF,CF,4) C C IP=LEF(1,K) XYZ(1,1)=XCOOR((IP-1)*(IES+1) +1) XYZ(2,1)=XCOOR((IP-1)*(IES+1) +2) IP=LEF(2,K) XYZ(1,2)=XCOOR((IP-1)*(IES+1) +1) XYZ(2,2)=XCOOR((IP-1)*(IES+1) +2) LE(1)=IPADL(LEF(1,K)) LE(2)=IPADL(LEF(2,K)) C WRITE(IOIMP,*)' LE.... ',(LE(MM),MM=1,4) DO 5 I=1,2 NF=LE(I) UIX(I)= UN(NF,1) UIY(I)= UN(NF,2) 5 CONTINUE C C CALCUL DES LONGUEURS DE L'ELEMENT : AL1 C IF (IAXI.EQ.0) THEN R1=1.D0 ELSEIF (IAXI.EQ.1) THEN R1=(XYZ(2,1)+XYZ(2,2))*XPI ELSEIF (IAXI.EQ.2) THEN R1=(XYZ(1,1)+XYZ(1,2))*XPI ENDIF AL1X=XYZ(1,2)-XYZ(1,1) AL1Y=XYZ(2,2)-XYZ(2,1) AL1=SQRT(AL1X*AL1X+AL1Y*AL1Y) C C CALCUL DES COMPOSANTES TANGENTIELLES A LA PAROI DES VITESSES C ET DE LA MOYENNE DE LEURS DIFFERENCES FLUIDE-PAROI : UTM C UT1=UIX(1)*UIX(1) + UIY(1)*UIY(1) UT2=UIX(2)*UIX(2) + UIY(2)*UIY(2) UT1=SQRT(UT1) UT2=SQRT(UT2) UTM=(UT1+UT2)/2.D0 C C CALCUL DE U ETOILE A PARTIR UTM C C --------------------CAROLI---I-------------------------- DO 57 I=1,2 YPLUS=YP*UET(K)/ANUL YPLUS=ABS(YPLUS)+1.D-5 IF(YPLUS.GT.100.0D0) THEN UPLUS=5.5D0+(LOG(YPLUS))/0.41D0 ELSE IY=INT(YPLUS) RY=YPLUS-IY UPLUS=UPIU(IY+1)+RY*(UPIU(IY+2)-UPIU(IY+1)) ENDIF C ---- RELAX SU UET ------------------------------------ UET(NK)=0.5D0*UET(NK)+0.5D0*(ABS(UTM/(UPLUS+1.D-5))) 57 CONTINUE C --------------------CAROLI----F------------------ C************************************************************************ C CALCUL Q D M C************************************************************************ C CALCUL DU COEFFICIENT AK A PARTIR DE UTM AUTM=ABS(UTM) AK=0.D0 IF(AUTM.GT.1.D-10)AK=UET(NK)*UET(NK)/AUTM C CALCUL DE LA MATRICE VITESSE DO 82 I=1,2 DO 72 J=1,2 AF(J,I)=AK*AL1*AF(J,I)*R1 72 CONTINUE 82 CONTINUE C C SATURATION VITESSE C C IF(IPAT.LE.2)WRITE(IOIMP,*)' UTM=',UTM,AL1 DO 54 J=1,2 AF(1,J)=AF(1,J)*UT1 AF(2,J)=AF(2,J)*UT2 54 CONTINUE DO 65 I=1,2 NF=LE(I) FF1=(AF(1,I)+AF(2,I))*AL1X/AL1 FF2=(AF(1,I)+AF(2,I))*AL1Y/AL1 C if(k.ge.10.AND.k.le.15)then C WRITE(IOIMP,*)' FF1,FF2=',ff1,ff2 C endif F(NF,1)=F(NF,1)+SIGN(FF1,UN(NF,1)) F(NF,2)=F(NF,2)+SIGN(FF2,UN(NF,2)) C------------------------------------------------------------------------ C CALCUL Q D M FIN C------------------------------------------------------------------------ 65 CONTINUE C C CALCUL DE KP ET DE EPSILONP C YPLUS2=-(YPLUS+0.01D0)*(YPLUS+0.01D0)*0.0017D0 DO 66 I=1,2 NF=LE(I) DK(NF)=XGRAND DE(NF)=XGRAND TK(NF)=UET(NK)*UET(NK)/SQCNU TE(NF)=UET(NK)*UET(NK)*UET(NK)/ & (AKER*YP*(1.D0-EXP(YPLUS2))) 66 CONTINUE 50 CONTINUE RETURN C ********* C * 3D * C ********* 300 CONTINUE DO 350 K=1,NEL NK=K0+K KC=1+(1-IKC)*(NK-1) C UTM=0.D0 DO 35 I=1,NPF NFA=IPADL(LEF(I,K)) UIX(I)=UN(NFA,1)+XPETIT UIY(I)=UN(NFA,2)+XPETIT UIZ(I)=UN(NFA,3)+XPETIT UFI(I)=UIX(I)*UIX(I)+UIY(I)*UIY(I)+UIZ(I)*UIZ(I) UFI(I)=SQRT(UFI(I)) UTM=UTM+UFI(I) 35 CONTINUE UTM=UTM/NPF AIRF=VOLF(NK) C C CALCUL DE U ETOILE A PARTIR UTM C C --------------------CAROLI---I-------------------------- DO 357 I=1,2 YPLUS=YP*UET(NK)/ANUL YPLUS=ABS(YPLUS)+1.D-5 IF(YPLUS.GT.100.0D0) THEN UPLUS=5.5D0+(LOG(YPLUS))/0.41D0 ELSE IY=INT(YPLUS) RY=YPLUS-IY UPLUS=UPIU(IY+1)+RY*(UPIU(IY+2)-UPIU(IY+1)) ENDIF C ---- RELAX SU UET ------------------------------------ UET(NK)=0.5D0*UET(NK)+0.5D0*(ABS(UTM/(UPLUS+1.D-5))) 357 CONTINUE C --------------------CAROLI----F------------------ C************************************************************************ C CALCUL Q D M C************************************************************************ C CALCUL DU COEFFICIENT AK A PARTIR DE UTM AUTM=ABS(UTM) AK=0.D0 IF(AUTM.GT.1.D-10)AK=UET(NK)*UET(NK)/AUTM C DO 365 I=1,NPF NFA=IPADL(LEF(I,K)) F(NFA,1)=F(NFA,1)+AK/NPF*UIX(I)/UFI(I)*AIRF F(NFA,2)=F(NFA,2)+AK/NPF*UIY(I)/UFI(I)*AIRF F(NFA,3)=F(NFA,3)+AK/NPF*UIZ(I)/UFI(I)*AIRF 365 CONTINUE C------------------------------------------------------------------------ C CALCUL Q D M FIN C------------------------------------------------------------------------ C C CALCUL DE KP ET DE EPSILONP C YPLUS2=-(YPLUS+0.01D0)*(YPLUS+0.01D0)*0.0017D0 DO 366 I=1,NPF NFA=IPADL(LEF(I,K)) DK(NFA)=XGRAND DE(NFA)=XGRAND TK(NFA)=UET(NK)*UET(NK)/SQCNU TE(NFA)=UET(NK)*UET(NK)*UET(NK)/ & (AKER*YP*(1.D0-EXP(YPLUS2))) 366 CONTINUE 350 CONTINUE RETURN 9990 CONTINUE WRITE(IOIMP,*)' IKC ',IKC WRITE(IOIMP,*) $ ' COEF D''UN MAUVAIS TYPE : ON VEUT (SCAL OU VECT DOMA)' CALL ARRET(0) RETURN 1002 FORMAT(10(1X,1PE11.4)) 1001 FORMAT(20(1X,I5)) END