Sources Esope | Cast3M

xcvfpu
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 &lt;,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&lt;- UN(NPTS)
C UETS(I) tableau local ; U* aux sommets de l'élément     &lt;- 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