1 : C FPTAU SOURCE PV 08/09/11 21:15:54 6150 2 : SUBROUTINE FPTAU(IDIM,NP,NBEL,NUM,IPADT,XCOOR, 3 : & U,NPTI,NC,IPADS,NPTS,ANUT,AMUEF, 4 : & YP,AN,RO,IKR,AMU,IKM,UET,YPS,ITYPEL,TKE,NBINC) 5 : IMPLICIT INTEGER(I-N) 6 : IMPLICIT REAL*8 (A-H,O-Z) 7 : C....................................................................... 8 : C Ut(NP,NBEL) module vitesse tangente = |U - Un| 9 : C Un(NP,NBEL) module vitesse normale = Un 10 : C YP(NP) distance la paroi 11 : C UET(*) uet aux noeuds frontieres 12 : C AN(IDIM,NP) normale 13 : C ANN module de la normale 14 : C NUM(NP) connectivite element appuyé 15 : C NMS(NS) connectivite relative (frontiere) 16 : C NMY(NY) connectivite complementaire 17 : C RO(*) densite 18 : C AMUEF(NP,NBEL) Viscosite effective aux noeuds/elt 19 : C 20 : C 21 : C1/ ce qu'il faudrait faire c'est attribuer une normale pour chaque point 22 : C de l'element appuye -> c'est fait 23 : C....................................................................... 24 : 25 : -INC CCREEL 26 : DIMENSION LKKE(1000) 27 : real*8 ann(3),Xa(9),Ya(9),Za(9) 28 : DIMENSION XCOOR(*),XYZ(3,27) 29 : DIMENSION NUM(NP,NBEL),NMS(9),NMY(27) 30 : DIMENSION U(NPTI,NC),IPADT(*),VS(27,3),IPADS(*) 31 : DIMENSION ANUT(NPTI),AMUEF(NP,NBEL) 32 : DIMENSION YP(NP,NBEL),AN(IDIM,NP,NBEL) 33 : DIMENSION UT(27),UN(27) 34 : DIMENSION RO(*),AMU(*) 35 : DIMENSION UET(NPTS),YPS(NPTS),TKE(NPTS,*) 36 : DIMENSION NUMFA(27) 37 : DATA EPS/1.D-4/,NBIT/100/ 38 : c----------------------------------------------------------------------- 39 : c------- Caracterisation des elements ---------------------------------- 40 : c----------------------------------------------------------------------- 41 : DIMENSION NMF9(3,4),N1F9(3,4),N2F9(3,4) 42 : DATA NMF9/1,2,3,3,4,5,5,6,7,7,8,1/ 43 : DATA N1F9/8,9,4,2,9,6,4,9,8,6,9,2/,N2F9/7,6,5,1,8,7,3,2,1,5,4,3/ 44 : c----------------------------------------------------------------------- 45 : DIMENSION NMF7(3,3),N1F7(3,3),N2F7(1,3) 46 : DATA NMF7/1,2,3,3,4,5,5,6,1/ 47 : DATA N1F7/6,7,4,2,7,6,4,7,2/,N2F7/5,1,3/ 48 : c----------------------------------------------------------------------- 49 : DIMENSION NMF8(4,6),N1F8(4,6),N2F8(4,6) 50 : DATA NMF8/1,2,3,4,2,3,7,6,5,6,7,8,1,4,8,5,1,2,6,5,4,3,7,8/ 51 : c? DATA N1F8/6,7,4,2,7,6,4,7,2/,N2F8/5,1,3/ 52 : c----------------------------------------------------------------------- 53 : C CONSTANTES PHYSIQUES 54 : DATA CNU,AKER /0.09D0,0.41D0/ 55 : c----------------------------------------------------------------------- 56 : c----------------------------------------------------------------------- 57 : 58 : SQCNU=SQRT(CNU) 59 : 60 : c write(6,*)' DEBUT FPTAU : NBINC=',NBINC 61 : c write(6,*)'NES=',nes,' NPG=',npg,' idim=',idim,' np=',np 62 : c write(6,*)'NPTS=',NPTS 63 : 64 : DEUPI=1.D0 65 : cc IF(IAXI.NE.0)DEUPI=2.D0*XPI 66 : IERC=0 67 : ERRMAX=0.D0 68 : 69 : KT=0 70 : DO 2 K=1,NBEL 71 : c write(6,*)' DEBUT BCL sur K ',K 72 : NS=0 73 : NY=0 73 : CALL INITI(NMS,NP,0) 73 : c CALL INITI(MMS,NP,0) 73 : CALL INITI(NMY,NP,0) 74 : c write(6,*)' ELEMENT K=',K 75 : 76 : DO 20 I=1,NP 77 : J=NUM(I,K) 78 : JI=IPADT(NUM(I,K)) 79 : JM=(1-IKM)*(JI-1)+1 80 : c JS=IPADS(NUM(I,K)) 81 : c write(6,*)'I=',I,' num(i=',NUM(I,K),' Ipads()=', 82 : c & IPADS(J) 83 : IF(IPADS(J).NE.0)THEN 84 : NS=NS+1 85 : NMS(NS)=I 86 : c MMS(NS)=I 87 : c write(6,*)'I=',I,' num(i=',NUM(I,K),' Ipads()=', 88 : c & IPADS(J),' NS=',NS,' ',NMS(NS) 89 : ELSE 90 : NY=NY+1 91 : NMY(NY)=I 92 : ENDIF 93 : DO 10 N=1,IDIM 94 : XYZ(N,I) = XCOOR((J-1)*(IDIM+1)+N) 95 : VS(I,N)=U(JI,N) 96 : 10 CONTINUE 97 : c UETK(I,K)=VET(JS) 98 : c UETK(I,K)=0.D0 99 : c ROS(I) = RO(JR) 100 : c AMUS(I)= AMU(JM) 101 : 20 CONTINUE 102 : c write(6,*)'K=',K,' MMS=',(mms(ijk),ijk=1,nS) 103 : 104 : c write(6,*)'K=',K,' XYZ=',(XYZ(1,I),I=1,np) 105 : c write(6,*)'K=',K,' XYZ=',(XYZ(2,I),I=1,np) 106 : c write(6,*)'K=',K,' NP=',NP,' NS=',NS,' NY=',NY 107 : c write(6,*)'K=',K,' NUM ',(num(i,k),i=1,np) 108 : c write(6,*)'K=',K,' NMSO',(nms(i),i=1,nS) 109 : c write(6,*)'K=',K,' NMS ',(num(nms(i),k),i=1,nS) 110 : c write(6,*)'K=',K,' NMY ',(num(NMY(I),k),I=1,NY) 111 : 112 : IF(NY.EQ.1)GO TO 2 113 : IF(NS.EQ.1)GO TO 2 114 : IF(NS.LE.2.AND.IDIM.EQ.3)GO TO 2 115 : KT=KT+1 116 : c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. 117 : c._._ On Calcule la normale au SEG2 (NMS) et YP 118 : c._._ (La normale est la meme pour les deux pts) 119 : IF(IDIM.EQ.2.AND.NS.EQ.2)THEN 120 : 121 : Xa(1)=XYZ(1,NMS(1)) 122 : Ya(1)=XYZ(2,NMS(1)) 123 : Xa(2)=XYZ(1,NMS(2)) 124 : Ya(2)=XYZ(2,NMS(2)) 125 : c a=(Ya(1)-Ya(2))/(Xa(1)-Xa(2)) 126 : c a=1. 127 : aa=Ya(1)-Ya(2) 128 : bb=Xa(1)-Xa(2) 129 : AN(1,NMS(1),K)=aa 130 : AN(2,NMS(1),K)=(-1.)*bb 131 : ann(1)=((aa*aa)+(bb*bb))**0.5D0 132 : AN(1,NMS(1),K)=AN(1,NMS(1),K)/ann(1) 133 : AN(2,NMS(1),K)=AN(2,NMS(1),K)/ann(1) 134 : DO 241 I=1,NY 135 : Xq=XYZ(1,NMY(I)) 136 : Yq=XYZ(2,NMY(I)) 137 : YP(NMY(I),K)=abs((aa*Xq)-(bb*Yq)+(bb*Ya(1))-(aa*Xa(1)))/ann(1) 138 : 241 CONTINUE 139 : c la normale est la meme pour tous les pts (SEG2) 140 : DO 251 I=2,NS 141 : AN(1,NMS(I),K)=AN(1,NMS(1),K) 142 : AN(2,NMS(I),K)=AN(2,NMS(1),K) 143 : 251 CONTINUE 144 : DO 252 I=1,NY 145 : AN(1,NMY(I),K)=AN(1,NMS(1),K) 146 : AN(2,NMY(I),K)=AN(2,NMS(1),K) 147 : 252 CONTINUE 148 : c._._ FIN On Calcule la normale au SEG2 (NMS) et YP 149 : c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. 150 : 151 : c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. 152 : c._._ On Calcule la normale au SEG3 (NMS) et YP 153 : c._._ (La normale n'est pas la meme pour tous les pts) 154 : ELSEIF(IDIM.EQ.2.AND.NS.EQ.3)THEN 155 : 156 : C._._ Mais quelle est donc la face? 156 : call initi(numfa,np,0) 157 : do 2300 I=1,NS 158 : numfa(nms(I))=1 159 : 2300 continue 160 : nfa=numfa(2)+numfa(4)*2+numfa(6)*3+numfa(8)*4 161 : c write(6,*)'NFA=',NFA 162 : 163 : Xa(1)=XYZ(1,NMF9(1,NFA)) 164 : Ya(1)=XYZ(2,NMF9(1,NFA)) 165 : Xa(2)=XYZ(1,NMF9(2,NFA)) 166 : Ya(2)=XYZ(2,NMF9(2,NFA)) 167 : Xa(3)=XYZ(1,NMF9(3,NFA)) 168 : Ya(3)=XYZ(2,NMF9(3,NFA)) 169 : 170 : a=(Ya(1)-Ya(2))/(Xa(1)-Xa(2)) 171 : AN(1,NUM(NMF9(1,NFA),K),K)=a 172 : AN(2,NUM(NMF9(1,NFA),K),K)=-1.D0 173 : ann(1)=(a*a+1.D0)**0.5D0 174 : AN(1,NUM(NMF9(1,NFA),K),K)=AN(1,NUM(NMF9(1,NFA),K),K)/ann(1) 175 : AN(2,NUM(NMF9(1,NFA),K),K)=AN(2,NUM(NMF9(1,NFA),K),K)/ann(1) 176 : 177 : b=(Ya(1)-Ya(3))/(Xa(1)-Xa(3)) 178 : AN(1,NUM(NMF9(2,NFA),K),K)=b 179 : AN(2,NUM(NMF9(2,NFA),K),K)=-1.D0 180 : ann(2)=(b*b+1.D0)**0.5D0 181 : AN(1,NUM(NMF9(2,NFA),K),K)=AN(1,NUM(NMF9(2,NFA),K),K)/ann(2) 182 : AN(2,NUM(NMF9(2,NFA),K),K)=AN(2,NUM(NMF9(2,NFA),K),K)/ann(2) 183 : 184 : c=(Ya(3)-Ya(2))/(Xa(3)-Xa(2)) 185 : AN(1,NUM(NMF9(3,NFA),K),K)=c 186 : AN(2,NUM(NMF9(3,NFA),K),K)=-1.D0 187 : ann(3)=(c*c+1.D0)**0.5D0 188 : AN(1,NUM(NMF9(3,NFA),K),K)=AN(1,NUM(NMF9(3,NFA),K),K)/ann(3) 189 : AN(2,NUM(NMF9(3,NFA),K),K)=AN(2,NUM(NMF9(3,NFA),K),K)/ann(3) 190 : 191 : c write(6,*)'ann=',(ann(kk),kk=1,ns) 192 : DO 231 I=1,NS 193 : Xq=XYZ(1,N1F9(I,NFA)) 194 : Yq=XYZ(2,N1F9(I,NFA)) 195 : YP(N1F9(I,NFA),K)=abs(a*Xq-Yq+Ya(I)-a*Xa(I))/ann(I) 196 : Xq=XYZ(1,N2F9(I,NFA)) 197 : Yq=XYZ(2,N2F9(I,NFA)) 198 : YP(N2F9(I,NFA),K)=abs(a*Xq-Yq+Ya(I)-a*Xa(I))/ann(I) 199 : 231 CONTINUE 200 : 201 : c write(6,*)'NMY=',(NMY(I),I=1,NY) 202 : c write(6,*)'YP=',(YP(NMY(I),K),I=1,NY) 203 : c write(6,*)'YP=',(YP(NMS(I),K),I=1,NS) 204 : c la normale n'est pas la meme pour tous les pts (SEG3) 205 : DO 232 I=1,NS 206 : AN(1,N1F9(I,NFA),K)=AN(1,NMS(I),K) 207 : AN(2,N1F9(I,NFA),K)=AN(2,NMS(I),K) 208 : AN(1,N2F9(I,NFA),K)=AN(1,NMS(I),K) 209 : AN(2,N2F9(I,NFA),K)=AN(2,NMS(I),K) 210 : 232 CONTINUE 211 : c._._ On Calcule la normale au SEG3 (NMS) et YP 212 : c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. 213 : 214 : c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. 215 : c._._ On Calcule les normales aux CUB8,PRI6,TET4 et PYR5 216 : c._._ (La normale est la meme pour tous les pts) 217 : ELSEIF( ITYPEL.EQ.14 218 : & .OR.ITYPEL.EQ.16 219 : & .OR.ITYPEL.EQ.23 220 : & .OR.ITYPEL.EQ.25 )THEN 221 : 222 : Xa(1)=XYZ(1,NMS(1)) 223 : Ya(1)=XYZ(2,NMS(1)) 224 : Za(1)=XYZ(3,NMS(1)) 225 : 226 : Xa(2)=XYZ(1,NMS(2)) 227 : Ya(2)=XYZ(2,NMS(2)) 228 : Za(2)=XYZ(3,NMS(2)) 229 : 230 : Xa(3)=XYZ(1,NMS(3)) 231 : Ya(3)=XYZ(2,NMS(3)) 232 : Za(3)=XYZ(3,NMS(3)) 233 : 234 : AN(1,NMS(1),K)= 235 : & (Ya(2)-Ya(1))*(Za(3)-Za(1))-(Za(2)-Za(1))*(Ya(3)-Ya(1)) 236 : AN(2,NMS(1),K)= 237 : & (Za(2)-Za(1))*(Xa(3)-Xa(1))-(Xa(2)-Xa(1))*(Za(3)-Za(1)) 238 : AN(3,NMS(1),K)= 239 : & (Xa(2)-Xa(1))*(Za(3)-Za(1))-(Ya(2)-Ya(1))*(Xa(3)-Xa(1)) 240 : ann(1)=(AN(1,NMS(1),K)*AN(1,NMS(1),K) 241 : & +AN(2,NMS(1),K)*AN(2,NMS(1),K) 242 : & +AN(3,NMS(1),K)*AN(3,NMS(1),K))**0.5D0 243 : AN(1,NMS(1),K)=AN(1,NMS(1),K)/ann(1) 244 : AN(2,NMS(1),K)=AN(2,NMS(1),K)/ann(1) 245 : AN(3,NMS(1),K)=AN(3,NMS(1),K)/ann(1) 246 : 247 : DO 431 I=2,NS 248 : AN(1,NMS(I),K)=AN(1,NMS(1),K) 249 : AN(2,NMS(I),K)=AN(2,NMS(1),K) 250 : AN(3,NMS(I),K)=AN(3,NMS(1),K) 251 : 431 CONTINUE 252 : 253 : c write(6,*)'ann=',ann(1) 254 : DO 432 I=1,NY 255 : Xq=XYZ(1,NMY(I)) 256 : Yq=XYZ(2,NMY(I)) 257 : Zq=XYZ(3,NMY(I)) 258 : YP(NMY(I),K)=ABS( 259 : & (Xq-Xa(1))*AN(1,NMS(1),K)+(Yq-Ya(1))*AN(2,NMS(1),K) 260 : & +(Zq-Za(1))*AN(3,NMS(1),K)) 261 : c write(6,*)'(Xq-Xa(1))=',(Xq-Xa(1)),' (Yq-Ya(1))=',(Yq-Ya(1)), 262 : c &' (Zq-Za(1))=',(Zq-Za(1)) 263 : 432 CONTINUE 264 : 265 : c write(6,*)'AN ',AN(1,NMS(1),K),AN(2,NMS(1),K),AN(3,NMS(1),K) 266 : 267 : c write(6,*)'NMY=',(NMY(I),I=1,NY) 268 : c write(6,*)'YP=',(YP(NMY(I),K),I=1,NY) 269 : c write(6,*)'YP=',(YP(NMS(I),K),I=1,NS) 270 : c la normale n'est pas la meme pour tous les pts (SEG3) 271 : DO 433 I=1,NY 272 : AN(1,NMY(I),K)=AN(1,NMS(1),K) 273 : AN(2,NMY(I),K)=AN(2,NMS(1),K) 274 : AN(3,NMY(I),K)=AN(3,NMS(1),K) 275 : 433 CONTINUE 276 : c._._ FIN Calcul des normales aux CUB8,PRI6,TET4 et PYR5 277 : c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. 278 : 279 : c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. 280 : c._._ On Calcule les normales aux CU27,PR21,TE15 et PY19 281 : c._._ (La normale est la meme pour tous les pts) 282 : ELSEIF( ITYPEL.EQ.33 283 : & .OR.ITYPEL.EQ.34 284 : & .OR.ITYPEL.EQ.35 285 : & .OR.ITYPEL.EQ.36 )THEN 286 : 287 : Xa(1)=XYZ(1,NMS(1)) 288 : Ya(1)=XYZ(2,NMS(1)) 289 : Za(1)=XYZ(3,NMS(1)) 290 : 291 : Xa(2)=XYZ(1,NMS(2)) 292 : Ya(2)=XYZ(2,NMS(2)) 293 : Za(2)=XYZ(3,NMS(2)) 294 : 295 : Xa(3)=XYZ(1,NMS(3)) 296 : Ya(3)=XYZ(2,NMS(3)) 297 : Za(3)=XYZ(3,NMS(3)) 298 : 299 : Xa(4)=XYZ(1,NMS(NS)) 300 : Ya(4)=XYZ(2,NMS(NS)) 301 : Za(4)=XYZ(3,NMS(NS)) 302 : 303 : ANx1=(Ya(2)-Ya(1))*(Za(3)-Za(1))-(Za(2)-Za(1))*(Ya(3)-Ya(1)) 304 : ANy1=(Za(2)-Za(1))*(Xa(3)-Xa(1))-(Xa(2)-Xa(1))*(Za(3)-Za(1)) 305 : ANz1=(Xa(2)-Xa(1))*(Za(3)-Za(1))-(Ya(2)-Ya(1))*(Xa(3)-Xa(1)) 306 : ann(1)=(ANx1*ANx1+ANy1*ANy1+ANz1*ANz1)**0.5D0 307 : 308 : ANx2=(Ya(2)-Ya(1))*(Za(4)-Za(1))-(Za(2)-Za(1))*(Ya(4)-Ya(1)) 309 : ANy2=(Za(2)-Za(1))*(Xa(4)-Xa(1))-(Xa(2)-Xa(1))*(Za(4)-Za(1)) 310 : ANz2=(Xa(2)-Xa(1))*(Za(4)-Za(1))-(Ya(2)-Ya(1))*(Xa(4)-Xa(1)) 311 : ann(2)=(ANx2*ANx2+ANy2*ANy2+ANz2*ANz2)**0.5D0 312 : 313 : IF(ann(1).GT.ann(2))THEN 314 : AN(1,NMS(1),K)=ANx1/ann(1) 315 : AN(2,NMS(1),K)=ANy1/ann(1) 316 : AN(3,NMS(1),K)=ANz1/ann(1) 317 : ELSE 318 : AN(1,NMS(1),K)=ANx2/ann(2) 319 : AN(2,NMS(1),K)=ANy2/ann(2) 320 : AN(3,NMS(1),K)=ANz2/ann(2) 321 : ann(1)=ann(2) 322 : ENDIF 323 : 324 : DO 531 I=2,NS 325 : AN(1,NMS(I),K)=AN(1,NMS(1),K) 326 : AN(2,NMS(I),K)=AN(2,NMS(1),K) 327 : AN(3,NMS(I),K)=AN(3,NMS(1),K) 328 : 531 CONTINUE 329 : 330 : c write(6,*)'ann=',ann(1) 331 : DO 532 I=1,NY 332 : Xq=XYZ(1,NMY(I)) 333 : Yq=XYZ(2,NMY(I)) 334 : Zq=XYZ(3,NMY(I)) 335 : YP(NMY(I),K)=ABS( 336 : & (Xq-Xa(1))*AN(1,NMS(1),K)+(Yq-Ya(1))*AN(2,NMS(1),K) 337 : & +(Zq-Za(1))*AN(3,NMS(1),K)) 338 : c write(6,*)'(Xq-Xa(1))=',(Xq-Xa(1)),' (Yq-Ya(1))=',(Yq-Ya(1)), 339 : c &' (Zq-Za(1))=',(Zq-Za(1)) 340 : 532 CONTINUE 341 : 342 : c write(6,*)'AN ',AN(1,NMS(1),K),AN(2,NMS(1),K),AN(3,NMS(1),K) 343 : 344 : c write(6,*)'NMS=',(NMS(I),I=1,NS) 345 : c write(6,*)'NMY=',(NMY(I),I=1,NY) 346 : c write(6,*)'YP=',(YP(NMS(I),K),I=1,NS) 347 : c write(6,*)'YP=',(YP(NMY(I),K),I=1,NY) 348 : c la normale n'est pas la meme pour tous les pts (SEG3) 349 : DO 533 I=1,NY 350 : AN(1,NMY(I),K)=AN(1,NMS(1),K) 351 : AN(2,NMY(I),K)=AN(2,NMS(1),K) 352 : AN(3,NMY(I),K)=AN(3,NMS(1),K) 353 : 533 CONTINUE 354 : c._._ FIN Calcul des normales aux CU27,PR21,TE15 et PY19 355 : c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. 356 : ELSE 357 : 358 : write(6,*)'Type d elements non prevus dans FPU' 359 : write(6,*)'K=',K,' IDIM=',IDIM,' NS=',NS 360 : 361 : RETURN 362 : ENDIF 363 : 364 : c write(6,*)' FIN CALCUL NORMALES' 365 : c._._ On Calcule la vitesse normale Un et la vitesse tangentielle Ut (NMY) 366 : DO 244 I=1,NP 367 : U0=0.D0 368 : U1=0.D0 368 : call initd(un,np,0.D0) 369 : DO 243 N=1,IDIM 370 : U0=U0+VS(I,N)*AN(N,I,K) 371 : U1=U1+(VS(I,N) - UN(I)*AN(N,I,K))**2.D0 372 : 243 CONTINUE 373 : UN(I)=U0 374 : UT(I)=U1**0.5D0 + 1.D-10 375 : 244 CONTINUE 376 : c._._ On Calcule la vitesse normale Un et la vitesse tangentielle Ut (NMY) 377 : 378 : c write(6,*)'AN=',AN(1,NMS(1),K),AN(2,NMS(1),K) 379 : c write(6,*)'YP=',(YP(NMY(I),K),I=1,NY) 380 : c write(6,*)'UN=',(UN(NMY(I)),I=1,NY) 381 : c write(6,*)'UT=',(UT(NMY(I)),I=1,NY) 382 : c write(6,*)'VS= ' 383 : c write(6,1002)(VS(I,1),I=1,NP) 384 : c write(6,1002)(VS(I,2),I=1,NP) 385 : 385 : cc CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG, 386 : cc & NES,IDIM,NP,NPG,IAXI,AIRE,AJ,SGN) 387 : c._._ On Calcule la matrice masse diagonale maillage appuye 387 : cc CALL INITD(DG,NPTA,0.D0) 388 : cc DO 221 I=1,NP 389 : cc AI=0.D0 390 : cc DO 220 LG=1,NPG 391 : cc PDR=PGSQ(LG)*DEUPI*RPG(LG) 392 : cc AI=AI+FN(I,LG)*PDR 393 : c220 CONTINUE 394 : cc JD=IPADA(NUM(I,K)) 395 : cc DG(JD)=DG(JD)+AI 396 : c221 CONTINUE 397 : c._._ On Calcule la matrice masse diagonale maillage appuye 398 : 399 : UETM=0.D0 400 : YPM =0.D0 401 : DO 225 II=1,NY 402 : I=NMY(II) 403 : IS=NMS(II) 404 : YP(I,K)=YP(I,K)+1.D-20 405 : c write(6,*)'YP ',YP(I,K) 406 : 407 : JI=IPADT(NUM(I,K)) 408 : c JS=IPADS(NUM(IS,K)) 409 : JR=(1-IKR)*(JI-1)+1 410 : JM=(1-IKM)*(JI-1)+1 411 : c YPLUS=RO(JR)*YP(I,K)*UETK(I,K)/AMU(JM) 412 : YPLUS=RO(JR)*YP(I,K)*UT(I)/AMU(JM) 413 : YPLUS=YPLUS+1.D-10 414 : 415 : c UET0=UETK(I,K) 416 : UET0=1.D-20 417 : c write(6,*)'M=0 YPLUS=',YPLUS,UT(I),AMU(JM),UET0 418 : w=0.5D0 419 : 420 : DO 30 M=1,NBIT 421 : IF(YPLUS.LE.5.D0)THEN 422 : UPLUS=YPLUS 423 : UETN=UT(I)/UPLUS 424 : YPLUS=w*YPLUS + (1.D0-w)*RO(JR)*YP(I,K)*UETN/AMU(JM) 425 : ELSE 426 : c Reichardt 427 : UPLUS=1.D0/AKER*Log(1.D0+AKER*yplus)+ 428 : & 7.8D0*(1.D0-exp(-yplus/11.D0)-(yplus/11.D0)*exp(-yplus/3.D0)) 429 : UETN=UT(I)/UPLUS 430 : YPLUS=RO(JR)*YP(I,K)*UETN/AMU(JM) 431 : ENDIF 432 : resid=ABS(UETN-UET0)/UETN 433 : UET0=UETN 434 : IF(resid.LE.EPS)GO TO 31 435 : 30 CONTINUE 436 : write(6,*)'Operateur FPU, NON CONVERGENCE en 40 iterations: eps=', 437 : & eps,' residu ->',resid 438 : 31 CONTINUE 439 : 440 : c write(6,*)' Convergence en ',m,' iterations' 441 : c UETK(I,K)=UETN 442 : c VET(JS)=UETN 443 : UETM = (UETM*FLOAT(II-1)+UETN)/FLOAT(II) 444 : YPM = (YPM *FLOAT(II-1)+YP(I,K))/FLOAT(II) 445 : AMUEF(I,K)=AMU(JM) 446 : IF(YPLUS.GT.5.D0)THEN 447 : AMUEF(I,K)=(RO(JR)*UETN*UETN*YP(I,K)/UT(I))+AMU(JM) 448 : ANUT(IPADT(NUM(I,K)))=AMUEF(I,K) 449 : ENDIF 450 : 451 : 225 CONTINUE 452 : c write(6,*)' APRES 225 NY=',NY 453 : 454 : DO 226 I=1,NS 455 : JI=IPADT(NUM(I,K)) 456 : JM=(1-IKM)*(JI-1)+1 457 : AMUEF(I,K)=AMU(JM) 458 : ANUT(IPADT(NUM(NMS(I),K)))=AMUEF(I,K) 459 : UET(IPADS(NUM(NMS(I),K)))=UETM 460 : YPS(IPADS(NUM(NMS(I),K)))=YPM 461 : 226 CONTINUE 462 : 463 : IF(NBINC.EQ.3)THEN 464 : DO 227 I=1,NS 465 : JI=IPADT(NUM(I,K)) 466 : JM=(1-IKM)*(JI-1)+1 467 : TKE(IPADS(NUM(NMS(I),K)),1)=UETM*UETM/SQCNU 468 : TKE(IPADS(NUM(NMS(I),K)),2)=UETM*UETM*UETM/(AKER*YPM) 469 : 227 CONTINUE 470 : ENDIF 471 : 472 : c write(6,*)'NS=',NS 473 : c write(6,1001)(LKKE(I),I=1,NS) 474 : c write(6,1001)(NUM(I,K),I=1,NP) 475 : c write(6,*)'UET K=',K,(UETK(I,K),I=1,NP) 476 : c write(6,*)'AMUEF K=',K,(AMUEF(I,K),I=1,NP) 477 : 478 : c._._ On deduit aussi Tau aux sommets 479 : c._._ On deduit Nut aux points de Gauss 480 : 481 : C................................................................. 482 : 483 : C...... FIN Boucle sur les elements Indice K 484 : c write(6,*)'. FIN Boucle sur les elements Indice K ',K 485 : 2 CONTINUE 486 : 487 : c write(6,*)'UET=' 488 : c write(6,1001)(UET(i),i=1,npts) 489 : 490 : 491 : c write(6,*)' Nb d element traite KT=',KT 492 : c write(6,*)'RETOUR FPTAU' 493 : RETURN 494 : 1002 FORMAT(10(1X,1PE11.4)) 495 : 1001 FORMAT(20(1X,I5)) 496 : END 497 : 498 : 499 :
© Cast3M 2003 - All rights reserved.
Disclaimer