1 : C XCVFPU SOURCE MAGN 06/04/11 21:17:42 5397 2 : SUBROUTINE XCVFPU(FN,GR,PG,XYZ,HR,PGSQ,RPG, 3 : & NES,IDIM,NP,NPG,IAXI,AJ,NBEL,NUM,XCOOR, 4 : & NPTI,IPADI,UN,NPTS,IPADS,UET,F,AMU,IKM,RO,IKR, 5 : & YP,DS,TK,TE,NBINC) 6 : IMPLICIT INTEGER(I-N) 7 : IMPLICIT REAL*8 (A-H,O-Z) 8 : 9 : C************************************************************************ 10 : C 11 : C SYNTAXE : 12 : C 13 : C FPU (RO,UN,MU,UET,YP <,VPAROI> ) VPAROI non fait 14 : C 15 : C 16 : C 17 : C RO cst ou au sommet ROS(NPTS) 18 : C ANU VISCOSITE CINEMATIQUE 19 : C UET U* 20 : C YP DISTANCE A LA PAROI 21 : C VPAROI VITESSE DE LA PAROI (PAR DEFAUT 0.D0) 22 : C 23 : C IDIM : dime espace NES : dime espace elt 24 : C NP : nb pt par elt, NPG : nb pt d'intégration 25 : C NPTS : nb pt de la frontière NPTI : nb pt total 26 : C NBEL : nb d'éléments 27 : C NUM(NP,NBEL) connectivités de l'élément (par type d'élément) 28 : C IPADS: Indirections maillage frontière 29 : C IPADI: Indirections maillage total 30 : C 31 : C 32 : C VS(I,3) tableau local ; vitesse aux sommets de l'élément<- UN(NPTS) 33 : C UETS(I) tableau local ; U* aux sommets de l'élément <- UET(NPTS) 34 : C UXL,UYL,UZL ux,uy,uz aux points de gauss 35 : C 36 : C DS(NP,NPTS) valeurs de U* intégrés aux noeuds : Somme_S Na U*a dS 37 : C TK(NP,NPTS) valeurs de Kp intégrés aux noeuds : 38 : C TE(NP,NPTS) valeurs de Ep intégrés aux noeuds : 39 : C sert à projeter les champs des points de gauss au sommet en divisant 40 : C ensuite par la matrice masse diagonale (Fait dans yfpu) 41 : C 42 : C 43 : C 44 : C 45 : C 46 : C************************************************************************ 47 : -INC CCREEL 48 : DIMENSION FN(NP,NPG),GR(IDIM,NP,NPG),PG(NPG),AJ(IDIM,IDIM,NPG) 49 : DIMENSION XYZ(IDIM,NP),HR(NES,NP,NPG),PGSQ(NPG),RPG(NPG) 50 : DIMENSION XCOOR(*),DS(NPTS) 51 : DIMENSION NUM(NP,NBEL) 52 : 53 : DIMENSION UN(NPTI,IDIM),RO(*),AMU(*) 54 : DIMENSION UET(NPTS),F(NPTS,IDIM),TK(NPTS),TE(NPTS) 55 : DIMENSION IPADI(*),IPADS(*) 56 : 57 : 58 : C CONSTANTES PHYSIQUES 59 : C 60 : DATA CNU,AKER /0.09D0,0.41D0/ 61 : DIMENSION VS(16,3),UETS(16),ROS(16),AMUS(16) 62 : DATA NITMAX/100/,W/0.7D0/,EPS/1.D-4/ 63 : 64 : SQCNU=SQRT(CNU) 65 : DEUPI=1.D0 66 : IF(IAXI.NE.0)DEUPI=2.D0*XPI 66 : c?? CALL INITD(DS,NPTS,0.D0) 67 : IERC=0 68 : ERRMAX=0.D0 69 : NBP=0 70 : NBIM=0 71 : DO 2 K=1,NBEL 72 : DO 20 I=1,NP 73 : J=NUM(I,K) 74 : JI=IPADI(NUM(I,K)) 75 : JR=(1-IKR)*(JI-1)+1 76 : JM=(1-IKM)*(JI-1)+1 77 : JS=IPADS(NUM(I,K)) 78 : 79 : DO 10 N=1,IDIM 80 : XYZ(N,I) = XCOOR((J-1)*(IDIM+1)+N) 81 : VS(I,N)=UN(JI,N) 82 : 10 CONTINUE 83 : UETS(I)=UET(JS) 84 : ROS(I) = RO(JR) 85 : AMUS(I)= AMU(JM) 86 : 87 : 20 CONTINUE 88 : 88 : CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG, 89 : & NES,IDIM,NP,NPG,IAXI,AIRE,AJ,SGN) 90 : 91 : DO 60 J=1,NP 92 : JS=IPADS(NUM(J,K)) 93 : FX=0.D0 94 : FY=0.D0 95 : FZ=0.D0 96 : FU=0.D0 97 : FK=0.D0 98 : FE=0.D0 99 : 100 : DO 50 L=1,NPG 101 : 102 : UXL=0.D0 103 : UYL=0.D0 104 : UZL=0.D0 105 : UETL=0.D0 106 : AMUL=0.D0 107 : ROL =0.D0 108 : 109 : DO 5 I=1,NP 110 : UXL=UXL+FN(I,L)*VS(I,1) 111 : UYL=UYL+FN(I,L)*VS(I,2) 112 : IF(IDIM.EQ.3)UZL=UZL+FN(I,L)*VS(I,3) 113 : UETL=UETL+FN(I,L)*UETS(I) 114 : AMUL=AMUL+FN(I,L)*AMUS(I) 115 : ROL =ROL +FN(I,L)*ROS(I) 116 : 5 CONTINUE 117 : UNL=UXL*AJ(1,IDIM,L)+UYL*AJ(2,IDIM,L) 118 : IF(IDIM.EQ.3)UNL=UNL+UZL*AJ(3,IDIM,L) 119 : UXNL=UNL*AJ(1,IDIM,L) 120 : UYNL=UNL*AJ(2,IDIM,L) 121 : UXTL=UXL-UXNL 122 : UYTL=UYL-UYNL 123 : * Y AURAIT IL UNE ERREUR UYX ET UYTL ET NON QUE DES X GBM 124 : UTL=(UXTL*UXTL+UYTL*UYTL) 125 : 126 : IF(IDIM.EQ.3)THEN 127 : UZNL=UNL*AJ(3,IDIM,L) 128 : UZTL=UZL-UZNL 129 : UTL=UTL+(UZTL*UZTL) 130 : ENDIF 131 : UTL=UTL**0.5D0+1.D-30 132 : TXL=UXTL/UTL 133 : TYL=UYTL/UTL 134 : TZL=UZTL/UTL 135 : 136 : UTM=ABS(UTL) 137 : ANUL = AMUL/ROL 138 : 139 : C 140 : C CALCUL DE U ETOILE A PARTIR UTM 141 : C 142 : C --------------------CAROLI---I-------------------------- 143 : UETLN=UETL 144 : DO 57 MI=1,NITMAX 145 : YPLUS=YP*UETLN/ANUL 146 : YPLUS=ABS(YPLUS) 147 : c Reichardt 148 : UPLUS=1.D0/AKER*Log(1.+AKER*yplus)+ 149 : & 7.8D0*(1.D0-exp(-yplus/11.D0)-(yplus/11.D0)*exp(-yplus/3.D0)) 150 : 151 : 152 : C ---- RELAX SU UET ------------------------------------ 153 : UETI=ABS(UTM/(UPLUS+1.D-5)) 154 : UETLNT=W*UETLN+(1.D0-W)*UETI 155 : ERR=ABS(UETLNT-UETLN)/MAX(UETLN,UETLNT) 156 : UETLN0=UETLN 157 : UETLN=UETLNT 158 : IF(ERR.LT.EPS)GO TO 58 159 : 160 : 57 CONTINUE 161 : 162 : IF(ERR.GT.ERRMAX)THEN 163 : C IF(L.EQ.1.AND.J.EQ.1.AND.K.LE.5)write(6,*)'UETLN',UETLN,UETLN0 164 : ERRMAX=ERR 165 : NBP=NBP+1 166 : IERC=1 167 : ENDIF 168 : 169 : 170 : 58 CONTINUE 171 : IF(NBIM.LT.MI)NBIM=MI 172 : C --------------------CAROLI----F------------------ 173 : 174 : C************************************************************************ 175 : C CALCUL Q D M 176 : C************************************************************************ 177 : 178 : C CALCUL DU COEFFICIENT AK A PARTIR DE UTM 179 : 180 : UETLN=UETLN+(ANUL/YP) 181 : AK=0.D0 182 : IF(UTM.GT.1.D-10)AK=UETLN*UETLN*FN(J,L)*PGSQ(L)*DEUPI*RPG(L) 183 : FX=FX-AK*TXL*ROL 184 : FY=FY-AK*TYL*ROL 185 : IF(IDIM.EQ.3)FZ=FZ-AK*TZL*ROL 186 : FU=FU+UETLN*FN(J,L)*PGSQ(L)*DEUPI*RPG(L) 187 : YPLUS2=-(YPLUS+0.01D0)*(YPLUS+0.01D0)*0.0017D0 188 : 189 : IF(YPLUS.LE.0.1)THEN 190 : Ret= 300. 191 : FKP=(0.057+0.05*(sqrt(Ret/1600.)))*(yplus*yplus) 192 : & *exp(-yplus/7.3) 193 : FKP=FKP+4.6*(1.-exp(-yplus/20.))/(4.*yplus/Ret+1.) 194 : & *(1.-exp(-((yplus/3.)**2.))) 195 : FKE = 1./AKER/((YPLUS**4.+15.**4.)**0.25) 196 : FK=FK+((UETLN*UETLN/SQCNU*FKP)*FN(J,L)*PGSQ(L)*DEUPI*RPG(L)) 197 : FE=FE+((ROL*UETLN*UETLN*UETLN*UETLN/ANUL*FKE) 198 : & *FN(J,L)*PGSQ(L)*DEUPI*RPG(L)) 199 : 200 : ELSE 201 : 202 : FK=FK+((UETLN*UETLN/SQCNU)*FN(J,L)*PGSQ(L)*DEUPI*RPG(L)) 203 : FE=FE+((UETLN*UETLN*UETLN/(AKER*YP*(1.D0-EXP(YPLUS2)))) 204 : C FE=FE+((UETLN*UETLN*UETLN/(AKER*YP )) 205 : & *FN(J,L)*PGSQ(L)*DEUPI*RPG(L)) 206 : 207 : ENDIF 208 : 50 CONTINUE 209 : 210 : F(JS,1)=F(JS,1)+FX 211 : F(JS,2)=F(JS,2)+FY 212 : IF(IDIM.EQ.3)F(JS,3)=F(JS,3)+FZ 213 : DS(JS)=DS(JS)+FU 214 : 215 : IF(NBINC.EQ.3)THEN 216 : C 217 : C CALCUL DE KP ET DE EPSILONP 218 : C 219 : TK(JS)=TK(JS)+FK 220 : TE(JS)=TE(JS)+FE 221 : ENDIF 222 : 223 : 60 CONTINUE 224 : 225 : 2 CONTINUE 226 : IF(IERC.NE.0) 227 : & write(6,*)' NON CONVERGENCE EN ',NITMAX,'ITERATIONS',ERRMAX,NBP 228 : 229 : write(6,*)' FPU : Nb max d iterations ',nbim 230 : 231 : RETURN 232 : 1002 FORMAT(10(1X,1PE11.4)) 233 : 1001 FORMAT(20(1X,I5)) 234 : END 235 : 236 : 237 : 238 : 239 : 240 : 241 :
© Cast3M 2003 - All rights reserved.
Disclaimer