1 : C FLOWM SOURCE KK2000 13/11/08 21:15:36 7860 2 : SUBROUTINE FLOWM(FC, FT, GAMA, ALFA, P, K0, A, AH, BH, CH, C, 3 : & SIGE0, KH, F0, FM, ETA) 4 : IMPLICIT INTEGER(I-N) 5 : IMPLICIT REAL*8(A-H,O-Z) 6 : 7 : -INC CCREEL 8 : 9 : C Subroutine de calcul de M = DF/DS 10 : 11 : C FC -> rc 12 : C Ft -> rt 13 : C Alfa -> AlfaT supposé 14 : 15 : REAL*8 SIGE0(6), KH, F0,PI3S6 16 : REAL*8 TRS, S(3,3), J2, RO, FC, FTB, FT 17 : REAL*8 KSIB, FM(6), RCB, REB, B0, B1, C0,C1, C2, C3 18 : REAL*8 C4, C5, J3, INTER, INTER2, TETA, DZERO, D1 19 : REAL*8 D2, Q, K0, KA, KSIH, KB, ROB 20 : REAL*8 DJ2DS(6), DROBDS(6), DKSIBD(6), DRCBDK 21 : REAL*8 DRCBDS(6), DB1DB0, DC0DB0, DC1DB0, DC2DB0 22 : REAL*8 DC3DB0, DC4DB0, DC5DB0, DD0DB0, DD1DB0, DD2DB0 23 : REAL*8 DQD0, DQDB0, DQD1, DQD2, DREBDK, DB0DKS 24 : REAL*8 DTEDJ2, DTEDJ3, DINTE(3,3), DJ3DS(6) 25 : REAL*8 DTETDS(6), D0DTET, D1DTET, D2DTET, DQDTET 26 : REAL*8 DQDS(6), DKBDKS, DFDS(6), DB0DS(6), M 27 : REAL*8 GAMA, ALFA, P, A, AH, BH, CH, C, INTER9 28 : REAL*8 DFDKH, DKBDK, ETA, DKDKH, INTER3, INTE10 29 : INTEGER I, J, K, MM 30 : 31 : 32 : C Calcul de la trace de la prédiction élastique 33 : TRS = 0.D0 34 : DO I = 1,3 35 : TRS = TRS + SIGE0(I) 36 : END DO 37 : 38 : 39 : C Calcul du déviateur de la prédiction 40 : S(1,1) = SIGE0(1) - 1.D0/3.D0 * TRS 41 : S(2,2) = SIGE0(2) - 1.D0/3.D0 * TRS 42 : S(3,3) = SIGE0(3) - 1.D0/3.D0 * TRS 43 : S(1,2) = SIGE0(4) 44 : S(1,3) = SIGE0(5) 45 : S(2,3) = SIGE0(6) 46 : S(2,1) = S(1,2) 47 : S(3,1) = S(1,3) 48 : S(3,2) = S(2,3) 49 : 50 : 51 : C Calcul de J2 52 : J2 = 0.D0 53 : DO I = 1,3 54 : DO J = 1,3 55 : J2 = J2 + 1.D0/2.D0 * S(I,J) * S(J,I) 56 : END DO 57 : END DO 58 : 59 : 60 : 61 : C*** condition de nullité de J2 : 62 : 63 : INTER2 = MAX(ABS(SIGE0(1)), ABS(SIGE0(2))) 64 : INTER2 = MAX(ABS(INTER2), ABS(SIGE0(3))) 65 : INTER2 = J2 / INTER2**2.D0 66 : 67 : IF (INTER2 .LT. 1.d-6) THEN 68 : J2 = 0.D0 69 : END IF 70 : 71 : 72 : RO = SQRT(2.D0)*SQRT(J2) 73 : ROB = RO / FC 74 : 75 : 76 : 77 : C**** Calcul de RCB 78 : 79 : FTB = FT / FC 80 : M = 3.D0 * (1.D0-(FTB)**(2.D0/GAMA))/(FTB + 81 : & 2.D0* (FTB)**(1.D0/GAMA)) 82 : KSIB = 1.D0/(SQRT(3.D0))*TRS/FC 83 : C = 1.d0 84 : 85 : INTER3 = M*M-12.D0*SQRT(3.D0)*M*KSIB+36.D0 86 : 87 : IF (INTER3 .LT. 0.D0) THEN 88 : write(6,*) 'INTER3', INTER3 89 : write(6,*) 'racine carre negative' 90 : END IF 91 : 92 : 93 : INTER3 = -M+SQRT(M*M-12.D0*SQRT(3.D0)*M*KSIB+36.D0) 94 : 95 : IF (INTER3 .LT. 0.D0 .AND. GAMA.LT.1.D0) THEN 96 : C write(6,*) 'Puissance GAMA (inf. 1) de nombre négatif' 97 : END IF 98 : 99 : C IF (INTER3 .LT. 0.) THEN 100 : C RCB = (1/6.)**GAMA*SQRT(2/3.)*(-M)**GAMA 101 : C WRITE(*,*) 'Expression alternative de RCB' 102 : C ELSE 103 : RCB =(1.D0/6.D0)**GAMA*SQRT(2.D0/3.D0)*(-M+ 104 : & SQRT(M*M-12.D0*SQRT(3.D0)*M*KSIB +36.D0))**GAMA 105 : 106 : C END IF 107 : 108 : 109 : 110 : C**** Calcul de Q 111 : 112 : INTER3 = M*M-3.D0*SQRT(3.D0)*M*KSIB+9.D0 113 : 114 : IF (INTER3 .LT. 0.D0) THEN 115 : write(6,*) 'racine carre negative' 116 : END IF 117 : 118 : INTER3 = -M+SQRT(M*M-3.D0*SQRT(3.D0)*M*KSIB+9.D0) 119 : 120 : IF (INTER3 .LT. 0.D0 .AND. GAMA.LT.1.D0) THEN 121 : C write(6,*) 'Puissance GAMA (inf. 1) de nombre négatif 2' 122 : END IF 123 : 124 : REB = (1.D0/3.D0)**GAMA*SQRT(2.D0/3.D0)*(-M+SQRT(M*M- 125 : & 3.D0*SQRT(3.D0)*M*KSIB +9.D0))**GAMA 126 : 127 : B0 = REB / RCB 128 : B1 = SQRT(3.D0)*(1.D0-ALFA)*B0/(1.D0+B0)+ 129 : & 2.D0*ALFA*B0/(SQRT(3.D0)) 130 : 131 : INTER9 = SQRT(3.D0)*B0/(1.D0+B0) 132 : 133 : INTE10 = 2.D0*B0/SQRT(3.D0) 134 : 135 : IF ((B0.GT.1.D0).OR.(B0.LE.0.5) 136 : & .OR.(B1.LE.INTER9).OR.(B1.GE.INTE10)) THEN 137 : write(6,*) 'ATENTION VALEUR B0 B1' 138 : write(6,*) B0, B1 139 : END IF 140 : 141 : C0 = (2.D0-SQRT(3.D0)*B1)*(2.D0*B0-SQRT(3.D0)*B1)/ 142 : & (B1*(1.D0+B0)- SQRT(3.D0)*B0)**2.D0 143 : C1 = 3.D0-C0*(1.D0+B0)**2.D0 144 : C2 = 1.D0+3.D0*C0*(1.D0-B0)**2.D0 145 : C3 = 2.D0*C0*SQRT(3.D0)*(1.D0-B0*B0) 146 : C4 = (1.D0+B0)*(1.D0-B0*C0) 147 : C5 = (1.D0-B0)*(1.D0-3.D0*B0*C0) 148 : C WRITE(6,*) 'B0,B1' 149 : C WRITE(6,*) B0,B1 150 : C WRITE(6,*) 'C0,C1,C2,C3,C4,C5' 151 : C WRITE(6,*) C0,C1,C2,C3,C4,C5 152 : 153 : J3 =0.D0 154 : DO I = 1,3 155 : DO J = 1,3 156 : DO K = 1,3 157 : J3 = J3 + 1.D0/3.D0 * S(I,J)*S(J,K)*S(K,I) 158 : ENDDO 159 : ENDDO 160 : ENDDO 161 : 162 : INTER = -3.D0*SQRT(3.D0)/2.D0 * J3 / (J2)**(1.5) 163 : 164 : 165 : C*** Correction sur le sin < 1 166 : C*** correction triaxial si J2 =0..., Teta =0 167 : 168 : IF (INTER .GT. 1.D0) THEN 169 : INTER = 1.D0 170 : END IF 171 : IF (INTER .LT. -1.D0) THEN 172 : INTER = -1.D0 173 : END IF 174 : 175 : INTER2 = MAX(ABS(SIGE0(1)), ABS(SIGE0(2))) 176 : INTER2 = MAX(ABS(INTER2), ABS(SIGE0(3))) 177 : INTER2 = J2 / INTER2**2.D0 178 : 179 : C** cas de l'essai triaxial confine : J2 environ 0 donc 0 180 : 181 : IF (INTER2 .LT. 1d-6) THEN 182 : TETA =0.D0 183 : ELSE 184 : TETA = 1.D0/3.D0 * DASIN(INTER) 185 : END IF 186 : 187 : 188 : C**** fin correction 189 : 190 : DZERO = C1*(COS(TETA))**2.D0 - C2*(SIN(TETA))**2.D0 191 : & + C3 * SIN(TETA) * COS(TETA) 192 : D1 = 2.D0*(C4*SQRT(3.D0)*COS(TETA) - C5 * SIN(TETA)) 193 : D2 = B0 * (4.D0-3.D0*B0*C0) 194 : 195 : IF (B0.LE.0.5) THEN 196 : C WRITE(6,*) 'Q a valeur limite.' 197 : Q = 0.5 198 : ELSE 199 : 200 : INTER = D1*D1-4.D0*DZERO*D2 201 : 202 : 203 : IF (INTER .LT. 0.D0) THEN 204 : write(6,*) 'RACINE CARREE Q negative', INTER 205 : END IF 206 : 207 : Q = (D1 - (SQRT(D1*D1 - 4.D0*DZERO*D2)))/(2.D0*DZERO) 208 : 209 : END IF 210 : 211 : C**** Calcul de KB 212 : 213 : KA = K0 + (1.D0-K0) * SQRT(KH*(2.D0-KH)) 214 : KSIH = A / (1.D0-KA) 215 : 216 : C write(6,*) 'KSIH', KSIH 217 : C write(6,*) 'KH', KH 218 : 219 : KB = KA**(2.D0*P)*(1.D0-KSIB**2.D0/KSIH**2.D0) 220 : 221 : DKDKH = (1.D0-K0)*(1.D0-KH)/SQRT(KH*(2.D0-KH)) 222 : 223 : C write(6,*) 'DKDKH', DKDKH 224 : 225 : DKBDK = (2.D0*P)*KA**(2.D0*P-1.D0)*(1.D0-(1.D0- 226 : & KA)**2.D0*KSIB**2.D0/A**2.D0) 227 : & + KA**(2.D0*P)*(2.D0*(1.D0-KA)*KSIB**2.D0/A**2.D0) 228 : 229 : C write(6,*) 'DKBDK', DKBDK 230 : 231 : DFDKH = -RCB**2.D0*Q**2.D0*DKBDK*DKDKH 232 : 233 : ETA=DFDKH 234 : 235 : IF (KH.GE.1.D0) THEN 236 : ETA = 0.D0 237 : END IF 238 : 239 : 240 : 241 : C write(6,*) 'ETA', ETA 242 : 243 : C** fonction seuil 244 : 245 : F0=ROB * ROB - KB * RCB*RCB * Q*Q 246 : 247 : 248 : C**** Deuxième partie : calcul du flux m = DFDS 249 : 250 : C**** calcul de dF/dSigma 251 : 252 : C****** calcul de DROBDS 253 : 254 : DJ2DS(1) = S(1,1) 255 : DJ2DS(2) = S(2,2) 256 : DJ2DS(3) = S(3,3) 257 : 258 : C*** * 2 ? 259 : 260 : DJ2DS(4) = 2.D0*S(1,2) 261 : DJ2DS(5) = 2.D0*S(1,3) 262 : DJ2DS(6) = 2.D0*S(2,3) 263 : 264 : C** fin ? 265 : 266 : C*** DROB**2/DS 267 : 268 : DO I = 1,6 269 : DROBDS(I) = 2.D0/(FC*FC)*DJ2DS(I) 270 : END DO 271 : 272 : 273 : C****** calcul de DRBCDS 274 : 275 : DO I = 1,6 276 : DKSIBD(I) = 0.D0 277 : END DO 278 : DO I = 1,3 279 : DKSIBD(I) = 1.D0/(SQRT(3.D0)*FC) 280 : ENDDO 281 : 282 : DRCBDK = -GAMA*((-M+SQRT(M*M-12.D0*SQRT(3.D0)*M*KSIB+ 283 : & 36.D0*C)) /6.D0)**(GAMA-1.D0)*(SQRT(2.D0)*M)/ 284 : & (SQRT(M*M- 12.D0*SQRT(3.D0) *M*KSIB+36.D0*C)) 285 : 286 : DO I = 1,6 287 : DRCBDS(I) = DRCBDK * DKSIBD(I) 288 : END DO 289 : 290 : C******* calcul de DQDS 291 : 292 : DB1DB0 = SQRT(3.D0)*(1.D0-ALFA)/(1.D0+B0)**2.D0 + 293 : & 2.D0*ALFA/(SQRT(3.D0)) 294 : 295 : DC0DB0 = 1.D0/(B1*(1.D0+B0)-SQRT(3.D0)*B0)**4.D0* 296 : & ((-SQRT(3.D0)*DB1DB0* 297 : & (2.D0*B0 -SQRT(3.D0)*B1)+(2.D0-SQRT(3.D0)*B1)* 298 : & (2.D0- SQRT(3.D0)*DB1DB0)) *(B1*(1.D0+B0)- 299 : & SQRT(3.D0)*B0)**2.D0 - 300 : & (2.D0-SQRT(3.D0)*B1) 301 : & * (2.D0*B0 - SQRT(3.D0)*B1) 302 : & * (2.D0*(B1*(1.D0+B0)-SQRT(3.D0)*B0) 303 : & *(DB1DB0*(1.D0+B0)+B1-SQRT(3.D0)))) 304 : 305 : DC1DB0 = -DC0DB0*(1.D0+B0)**2.D0-2.D0*C0*(1.D0+B0) 306 : DC2DB0 = 3.D0*DC0DB0*(1.D0-B0)**2.D0+3.D0*C0*2.D0* 307 : & (1.D0-B0)*(-1.D0) 308 : DC3DB0 = 2.D0*DC0DB0*SQRT(3.D0)*(1.D0-B0*B0)+ 309 : & 2.D0*C0*SQRT(3.D0)*(-2.D0*B0) 310 : DC4DB0 = (1.D0-B0*C0)+(1.D0+B0)*(-1.D0)*(C0+B0*DC0DB0) 311 : DC5DB0 = -(1.D0-3.D0*B0*C0)+(1.D0-B0)*(-3.D0)*(C0+B0*DC0DB0) 312 : 313 : DD0DB0 = DC1DB0 * COS(TETA)*COS(TETA)-DC2DB0*SIN(TETA)* 314 : & SIN(TETA) + DC3DB0*SIN(TETA)*COS(TETA) 315 : DD1DB0 = 2.D0*(DC4DB0*SQRT(3.D0)*COS(TETA)-DC5DB0*SIN(TETA)) 316 : DD2DB0 = (4.D0-3.D0*B0*C0)+B0*(-3.D0*(C0+B0*DC0DB0)) 317 : 318 : DQD0 = (4.D0*DZERO*D2/(2.D0*SQRT(D1*D1-4.D0*DZERO*D2))- 319 : & (D1-SQRT(D1*D1-4.D0*DZERO*D2))) 320 : & /(2.D0*DZERO*DZERO) 321 : DQD1 = 1.D0/(2.D0*DZERO)*(1.D0-(D1)/(SQRT(D1*D1-4.D0*DZERO*D2))) 322 : DQD2 = 1.D0/(SQRT(D1*D1-4.D0*DZERO*D2)) 323 : 324 : DQDB0 = DQD0 * DD0DB0 + DQD1*DD1DB0 + DQD2 * DD2DB0 325 : 326 : DREBDK = -GAMA*((-M+SQRT(M*M-3.D0*SQRT(3.D0)*M*KSIB+ 327 : & 9.D0*C))/3.D0) 328 : & **(GAMA-1.D0) * (M/SQRT(2.D0))/SQRT(M*M- 329 : & 3.D0*SQRT(3.D0)*M*KSIB+9.D0*C) 330 : 331 : DB0DKS = 1.D0/(RCB*RCB)*(DREBDK*RCB-REB*DRCBDK) 332 : 333 : DO I =1,6 334 : DB0DS(I) = DB0DKS*DKSIBD(I) 335 : END DO 336 : 337 : C****** eviter J2 = 0 (triaxial, définition des dérivées) 338 : PI3S6= XPI/6.D0 339 : C write(6,*) 'absolu', ABS(ABS(TETA)-PI3S6) 340 : IF (J2 .GT. 0.D0) THEN 341 : IF (ABS(ABS(TETA)-PI3S6) .GT. 1.d-6) THEN 342 : DTEDJ2 = 3.D0*SQRT(3.D0)*J3/(4.D0*COS(3.D0*TETA))/ 343 : & (J2)**2.5 344 : DTEDJ3 = - SQRT(3.D0)/(2.D0*COS(3.D0*TETA)*J2**1.5) 345 : ELSE 346 : DTEDJ2 =0.D0 347 : DTEDJ3 = 0.D0 348 : ENDIF 349 : ENDIF 350 : DO I = 1,3 351 : DO J = 1,3 352 : DINTE(I,J) = 0.D0 353 : END DO 354 : END DO 355 : 356 : DO I = 1, 3 357 : DO J = 1,I 358 : DO MM = 1,3 359 : DINTE(I,J) = DINTE(I,J) + S(I,MM)*S(MM,J) 360 : END DO 361 : IF (I .EQ. J) THEN 362 : DINTE(I,J) = DINTE(I,J) - 2.D0/3.D0 * J2 363 : END IF 364 : END DO 365 : END DO 366 : 367 : DJ3DS(1) = DINTE(1,1) 368 : DJ3DS(2) = DINTE(2,2) 369 : DJ3DS(3) = DINTE(3,3) 370 : 371 : C x2 ? 372 : C 373 : DJ3DS(4) = 2.D0*DINTE(2,1) 374 : DJ3DS(5) = 2.D0*DINTE(3,1) 375 : C DJ3DS(6) = DINTE(3,2) 376 : DJ3DS(6) = 2.D0*DINTE(3,2) 377 : 378 : DO I = 1,6 379 : DTETDS(I) = DTEDJ2*DJ2DS(I)+DTEDJ3*DJ3DS(I) 380 : END DO 381 : 382 : D0DTET = -2.D0*C1*COS(TETA)*SIN(TETA)-2.D0*C2*SIN(TETA)* 383 : & COS(TETA)+C3*((COS(TETA))**2.D0-(SIN(TETA))**2.D0) 384 : D1DTET = 2.D0*(-C4*SQRT(3.D0)*SIN(TETA)-C5*COS(TETA)) 385 : D2DTET = 0.D0 386 : 387 : DQDTET = DQD0 * D0DTET + DQD1*D1DTET + DQD2 * D2DTET 388 : 389 : DO I=1,6 390 : DQDS(I) = DQDB0*DB0DS(I) + DQDTET * DTETDS(I) 391 : END DO 392 : 393 : C**** Calcul de DKBDKISB 394 : 395 : DKBDKS = (KA**(2.D0*P))*(-2.D0)*KSIB*(1.D0-KA)**2.D0 396 : & / (A**2.D0) 397 : 398 : 399 : IF (KH.GE.1.D0) THEN 400 : DO I = 1,6 401 : DFDS(I) = DROBDS(I) - 2.D0*DRCBDS(I)*RCB*KB*Q*Q 402 : & - 2.D0*DQDS(I)*Q*KB*RCB*RCB 403 : FM(I) = DFDS(I) 404 : END DO 405 : ELSE 406 : DO I = 1,6 407 : 408 : DFDS(I) = DROBDS(I) - 2.D0*DRCBDS(I)*RCB*KB*Q*Q 409 : & - 2.D0*DQDS(I)*Q*KB*RCB*RCB - DKBDKS 410 : & * DKSIBD(I)*RCB*RCB*Q*Q 411 : FM(I) = DFDS(I) 412 : END DO 413 : END IF 414 : C write(6,*) 'DROBDS', DROBDS(4), DROBDS(5), DROBDS(6) 415 : C write(6,*) 'DQDS' , DQDS(4) , DQDS(5) , DQDS(6) 416 : C write(6,*) 'DTETDS', DTETDS(4), DTETDS(5), DTETDS(6) 417 : 418 : C DO I =1,6 419 : C (6,*), 'DFDS', I, DFDS(I) 420 : C END DO 421 : 422 : END 423 : 424 : 425 : 426 : 427 : 428 : 429 : 430 : 431 : 432 : 433 :
© Cast3M 2003 - All rights reserved.
Disclaimer