C XCVFPU SOURCE GF238795 18/02/05 21:16:22 9726 SUBROUTINE XCVFPU(FN,GR,PG,XYZ,HR,PGSQ,RPG, & NES,IDIM,NP,NPG,IAXI,AJ,NBEL,NUM,XCOOR, & NPTI,IPADI,UN,NPTS,IPADS,UET,F,AMU,IKM,RO,IKR, & YP,DS,TK,TE,NBINC) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C************************************************************************ C C SYNTAXE : C C FPU (RO,UN,MU,UET,YP <,VPAROI> ) VPAROI non fait C C C C RO cst ou au sommet ROS(NPTS) 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 IDIM : dime espace NES : dime espace elt C NP : nb pt par elt, NPG : nb pt d'intégration C NPTS : nb pt de la frontière NPTI : nb pt total C NBEL : nb d'éléments C NUM(NP,NBEL) connectivités de l'élément (par type d'élément) C IPADS: Indirections maillage frontière C IPADI: Indirections maillage total C C C VS(I,3) tableau local ; vitesse aux sommets de l'élément<- UN(NPTS) C UETS(I) tableau local ; U* aux sommets de l'élément <- UET(NPTS) C UXL,UYL,UZL ux,uy,uz aux points de gauss C C DS(NP,NPTS) valeurs de U* intégrés aux noeuds : Somme_S Na U*a dS C TK(NP,NPTS) valeurs de Kp intégrés aux noeuds : C TE(NP,NPTS) valeurs de Ep intégrés aux noeuds : C sert à projeter les champs des points de gauss au sommet en divisant C ensuite par la matrice masse diagonale (Fait dans yfpu) C C C C C C************************************************************************ -INC CCREEL DIMENSION FN(NP,NPG),GR(IDIM,NP,NPG),PG(NPG),AJ(IDIM,IDIM,NPG) DIMENSION XYZ(IDIM,NP),HR(NES,NP,NPG),PGSQ(NPG),RPG(NPG) DIMENSION XCOOR(*),DS(NPTS) DIMENSION NUM(NP,NBEL) DIMENSION UN(NPTI,IDIM),RO(*),AMU(*) DIMENSION UET(NPTS),F(NPTS,IDIM),TK(NPTS),TE(NPTS) DIMENSION IPADI(*),IPADS(*) C CONSTANTES PHYSIQUES C DATA CNU,AKER /0.09D0,0.41D0/ DIMENSION VS(16,3),UETS(16),ROS(16),AMUS(16) DATA NITMAX/100/,W/0.7D0/,EPS/1.D-4/ REAL*8 UZTL UZTL=0.D0 SQCNU=SQRT(CNU) DEUPI=1.D0 IF(IAXI.NE.0)DEUPI=2.D0*XPI c?? CALL INITD(DS,NPTS,0.D0) IERC=0 ERRMAX=0.D0 NBP=0 NBIM=0 DO 2 K=1,NBEL DO 20 I=1,NP J=NUM(I,K) JI=IPADI(NUM(I,K)) JR=(1-IKR)*(JI-1)+1 JM=(1-IKM)*(JI-1)+1 JS=IPADS(NUM(I,K)) DO 10 N=1,IDIM XYZ(N,I) = XCOOR((J-1)*(IDIM+1)+N) VS(I,N)=UN(JI,N) 10 CONTINUE UETS(I)=UET(JS) ROS(I) = RO(JR) AMUS(I)= AMU(JM) 20 CONTINUE CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG, & NES,IDIM,NP,NPG,IAXI,AIRE,AJ,SGN) DO 60 J=1,NP JS=IPADS(NUM(J,K)) FX=0.D0 FY=0.D0 FZ=0.D0 FU=0.D0 FK=0.D0 FE=0.D0 DO 50 L=1,NPG UXL=0.D0 UYL=0.D0 UZL=0.D0 UETL=0.D0 AMUL=0.D0 ROL =0.D0 DO 5 I=1,NP UXL=UXL+FN(I,L)*VS(I,1) UYL=UYL+FN(I,L)*VS(I,2) IF(IDIM.EQ.3)UZL=UZL+FN(I,L)*VS(I,3) UETL=UETL+FN(I,L)*UETS(I) AMUL=AMUL+FN(I,L)*AMUS(I) ROL =ROL +FN(I,L)*ROS(I) 5 CONTINUE UNL=UXL*AJ(1,IDIM,L)+UYL*AJ(2,IDIM,L) IF(IDIM.EQ.3)UNL=UNL+UZL*AJ(3,IDIM,L) UXNL=UNL*AJ(1,IDIM,L) UYNL=UNL*AJ(2,IDIM,L) UXTL=UXL-UXNL UYTL=UYL-UYNL * Y AURAIT IL UNE ERREUR UYX ET UYTL ET NON QUE DES X GBM UTL=(UXTL*UXTL+UYTL*UYTL) IF(IDIM.EQ.3)THEN UZNL=UNL*AJ(3,IDIM,L) UZTL=UZL-UZNL UTL=UTL+(UZTL*UZTL) ENDIF UTL=UTL**0.5D0+1.D-30 TXL=UXTL/UTL TYL=UYTL/UTL TZL=UZTL/UTL UTM=ABS(UTL) ANUL = AMUL/ROL C C CALCUL DE U ETOILE A PARTIR UTM C C --------------------CAROLI---I-------------------------- UETLN=UETL DO 57 MI=1,NITMAX YPLUS=YP*UETLN/ANUL YPLUS=ABS(YPLUS) c Reichardt UPLUS=1.D0/AKER*Log(1.+AKER*yplus)+ & 7.8D0*(1.D0-exp(-yplus/11.D0)-(yplus/11.D0)*exp(-yplus/3.D0)) C ---- RELAX SU UET ------------------------------------ UETI=ABS(UTM/(UPLUS+1.D-5)) UETLNT=W*UETLN+(1.D0-W)*UETI ERR=ABS(UETLNT-UETLN)/MAX(UETLN,UETLNT) UETLN0=UETLN UETLN=UETLNT IF(ERR.LT.EPS)GO TO 58 57 CONTINUE IF(ERR.GT.ERRMAX)THEN C IF(L.EQ.1.AND.J.EQ.1.AND.K.LE.5)write(6,*)'UETLN',UETLN,UETLN0 ERRMAX=ERR NBP=NBP+1 IERC=1 ENDIF 58 CONTINUE IF(NBIM.LT.MI)NBIM=MI C --------------------CAROLI----F------------------ C************************************************************************ C CALCUL Q D M C************************************************************************ C CALCUL DU COEFFICIENT AK A PARTIR DE UTM UETLN=UETLN+(ANUL/YP) AK=0.D0 IF(UTM.GT.1.D-10)AK=UETLN*UETLN*FN(J,L)*PGSQ(L)*DEUPI*RPG(L) FX=FX-AK*TXL*ROL FY=FY-AK*TYL*ROL IF(IDIM.EQ.3)FZ=FZ-AK*TZL*ROL FU=FU+UETLN*FN(J,L)*PGSQ(L)*DEUPI*RPG(L) YPLUS2=-(YPLUS+0.01D0)*(YPLUS+0.01D0)*0.0017D0 IF(YPLUS.LE.0.1)THEN Ret= 300. FKP=(0.057+0.05*(sqrt(Ret/1600.)))*(yplus*yplus) & *exp(-yplus/7.3) FKP=FKP+4.6*(1.-exp(-yplus/20.))/(4.*yplus/Ret+1.) & *(1.-exp(-((yplus/3.)**2.))) FKE = 1./AKER/((YPLUS**4.+15.**4.)**0.25) FK=FK+((UETLN*UETLN/SQCNU*FKP)*FN(J,L)*PGSQ(L)*DEUPI*RPG(L)) FE=FE+((ROL*UETLN*UETLN*UETLN*UETLN/ANUL*FKE) & *FN(J,L)*PGSQ(L)*DEUPI*RPG(L)) ELSE FK=FK+((UETLN*UETLN/SQCNU)*FN(J,L)*PGSQ(L)*DEUPI*RPG(L)) FE=FE+((UETLN*UETLN*UETLN/(AKER*YP*(1.D0-EXP(YPLUS2)))) C FE=FE+((UETLN*UETLN*UETLN/(AKER*YP )) & *FN(J,L)*PGSQ(L)*DEUPI*RPG(L)) ENDIF 50 CONTINUE F(JS,1)=F(JS,1)+FX F(JS,2)=F(JS,2)+FY IF(IDIM.EQ.3)F(JS,3)=F(JS,3)+FZ DS(JS)=DS(JS)+FU IF(NBINC.EQ.3)THEN C C CALCUL DE KP ET DE EPSILONP C TK(JS)=TK(JS)+FK TE(JS)=TE(JS)+FE ENDIF 60 CONTINUE 2 CONTINUE IF(IERC.NE.0) & write(6,*)' NON CONVERGENCE EN ',NITMAX,'ITERATIONS',ERRMAX,NBP write(6,*)' FPU : Nb max d iterations ',nbim RETURN 1002 FORMAT(10(1X,1PE11.4)) 1001 FORMAT(20(1X,I5)) END