C FLOWM SOURCE CB215821 16/04/21 21:16:53 8920 SUBROUTINE FLOWM(FC, FT, GAMA, ALFA, P, K0, A, AH, BH, CH, C, & SIGE0, KH, F0, FM, ETA) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC CCREEL C Subroutine de calcul de M = DF/DS C FC -> rc C Ft -> rt C Alfa -> AlfaT supposé REAL*8 SIGE0(6), KH, F0,PI3S6 REAL*8 TRS, S(3,3), J2, RO, FC, FTB, FT REAL*8 KSIB, FM(6), RCB, REB, B0, B1, C0,C1, C2, C3 REAL*8 C4, C5, J3, INTER, INTER2, TETA, DZERO, D1 REAL*8 D2, Q, K0, KA, KSIH, KB, ROB REAL*8 DJ2DS(6), DROBDS(6), DKSIBD(6), DRCBDK REAL*8 DRCBDS(6), DB1DB0, DC0DB0, DC1DB0, DC2DB0 REAL*8 DC3DB0, DC4DB0, DC5DB0, DD0DB0, DD1DB0, DD2DB0 REAL*8 DQD0, DQDB0, DQD1, DQD2, DREBDK, DB0DKS REAL*8 DTEDJ2, DTEDJ3, DINTE(3,3), DJ3DS(6) REAL*8 DTETDS(6), D0DTET, D1DTET, D2DTET, DQDTET REAL*8 DQDS(6), DKBDKS, DFDS(6), DB0DS(6), M REAL*8 GAMA, ALFA, P, A, AH, BH, CH, C, INTER9 REAL*8 DFDKH, DKBDK, ETA, DKDKH, INTER3, INTE10 INTEGER I, J, K, MM C Calcul de la trace de la prédiction élastique TRS = 0.D0 DO I = 1,3 TRS = TRS + SIGE0(I) END DO C Calcul du déviateur de la prédiction S(1,1) = SIGE0(1) - 1.D0/3.D0 * TRS S(2,2) = SIGE0(2) - 1.D0/3.D0 * TRS S(3,3) = SIGE0(3) - 1.D0/3.D0 * TRS S(1,2) = SIGE0(4) S(1,3) = SIGE0(5) S(2,3) = SIGE0(6) S(2,1) = S(1,2) S(3,1) = S(1,3) S(3,2) = S(2,3) C Calcul de J2 J2 = 0.D0 DO I = 1,3 DO J = 1,3 J2 = J2 + 1.D0/2.D0 * S(I,J) * S(J,I) END DO END DO C*** condition de nullité de J2 : INTER2 = MAX(ABS(SIGE0(1)), ABS(SIGE0(2))) INTER2 = MAX(ABS(INTER2), ABS(SIGE0(3))) INTER2 = J2 / INTER2**2.D0 IF (INTER2 .LT. 1.d-6) THEN J2 = 0.D0 END IF RO = SQRT(2.D0)*SQRT(J2) ROB = RO / FC C**** Calcul de RCB FTB = FT / FC M = 3.D0 * (1.D0-(FTB)**(2.D0/GAMA))/(FTB + & 2.D0* (FTB)**(1.D0/GAMA)) KSIB = 1.D0/(SQRT(3.D0))*TRS/FC C = 1.d0 INTER3 = M*M-12.D0*SQRT(3.D0)*M*KSIB+36.D0 IF (INTER3 .LT. 0.D0) THEN write(6,*) 'INTER3', INTER3 write(6,*) 'racine carre negative' END IF INTER3 = -M+SQRT(M*M-12.D0*SQRT(3.D0)*M*KSIB+36.D0) IF (INTER3 .LT. 0.D0 .AND. GAMA.LT.1.D0) THEN C write(6,*) 'Puissance GAMA (inf. 1) de nombre négatif' END IF C IF (INTER3 .LT. 0.) THEN C RCB = (1/6.)**GAMA*SQRT(2/3.)*(-M)**GAMA C WRITE(*,*) 'Expression alternative de RCB' C ELSE RCB =(1.D0/6.D0)**GAMA*SQRT(2.D0/3.D0)*(-M+ & SQRT(M*M-12.D0*SQRT(3.D0)*M*KSIB +36.D0))**GAMA C END IF C**** Calcul de Q INTER3 = M*M-3.D0*SQRT(3.D0)*M*KSIB+9.D0 IF (INTER3 .LT. 0.D0) THEN write(6,*) 'racine carre negative' END IF INTER3 = -M+SQRT(M*M-3.D0*SQRT(3.D0)*M*KSIB+9.D0) IF (INTER3 .LT. 0.D0 .AND. GAMA.LT.1.D0) THEN C write(6,*) 'Puissance GAMA (inf. 1) de nombre négatif 2' END IF REB = (1.D0/3.D0)**GAMA*SQRT(2.D0/3.D0)*(-M+SQRT(M*M- & 3.D0*SQRT(3.D0)*M*KSIB +9.D0))**GAMA B0 = REB / RCB B1 = SQRT(3.D0)*(1.D0-ALFA)*B0/(1.D0+B0)+ & 2.D0*ALFA*B0/(SQRT(3.D0)) INTER9 = SQRT(3.D0)*B0/(1.D0+B0) INTE10 = 2.D0*B0/SQRT(3.D0) IF ((B0.GT.1.D0).OR.(B0.LE.0.5) & .OR.(B1.LE.INTER9).OR.(B1.GE.INTE10)) THEN write(6,*) 'ATENTION VALEUR B0 B1' write(6,*) B0, B1 END IF C0 = (2.D0-SQRT(3.D0)*B1)*(2.D0*B0-SQRT(3.D0)*B1)/ & (B1*(1.D0+B0)- SQRT(3.D0)*B0)**2.D0 C1 = 3.D0-C0*(1.D0+B0)**2.D0 C2 = 1.D0+3.D0*C0*(1.D0-B0)**2.D0 C3 = 2.D0*C0*SQRT(3.D0)*(1.D0-B0*B0) C4 = (1.D0+B0)*(1.D0-B0*C0) C5 = (1.D0-B0)*(1.D0-3.D0*B0*C0) C WRITE(6,*) 'B0,B1' C WRITE(6,*) B0,B1 C WRITE(6,*) 'C0,C1,C2,C3,C4,C5' C WRITE(6,*) C0,C1,C2,C3,C4,C5 J3 =0.D0 DO I = 1,3 DO J = 1,3 DO K = 1,3 J3 = J3 + 1.D0/3.D0 * S(I,J)*S(J,K)*S(K,I) ENDDO ENDDO ENDDO INTER = -3.D0*SQRT(3.D0)/2.D0 * J3 / (J2)**(1.5) C*** Correction sur le sin < 1 C*** correction triaxial si J2 =0..., Teta =0 IF (INTER .GT. 1.D0) THEN INTER = 1.D0 END IF IF (INTER .LT. -1.D0) THEN INTER = -1.D0 END IF INTER2 = MAX(ABS(SIGE0(1)), ABS(SIGE0(2))) INTER2 = MAX(ABS(INTER2), ABS(SIGE0(3))) INTER2 = J2 / INTER2**2.D0 C** cas de l'essai triaxial confine : J2 environ 0 donc 0 IF (INTER2 .LT. 1d-6) THEN TETA =0.D0 ELSE TETA = 1.D0/3.D0 * ASIN(INTER) END IF C**** fin correction DZERO = C1*(COS(TETA))**2.D0 - C2*(SIN(TETA))**2.D0 & + C3 * SIN(TETA) * COS(TETA) D1 = 2.D0*(C4*SQRT(3.D0)*COS(TETA) - C5 * SIN(TETA)) D2 = B0 * (4.D0-3.D0*B0*C0) IF (B0.LE.0.5) THEN C WRITE(6,*) 'Q a valeur limite.' Q = 0.5 ELSE INTER = D1*D1-4.D0*DZERO*D2 IF (INTER .LT. 0.D0) THEN write(6,*) 'RACINE CARREE Q negative', INTER END IF Q = (D1 - (SQRT(D1*D1 - 4.D0*DZERO*D2)))/(2.D0*DZERO) END IF C**** Calcul de KB KA = K0 + (1.D0-K0) * SQRT(KH*(2.D0-KH)) KSIH = A / (1.D0-KA) C write(6,*) 'KSIH', KSIH C write(6,*) 'KH', KH KB = KA**(2.D0*P)*(1.D0-KSIB**2.D0/KSIH**2.D0) DKDKH = (1.D0-K0)*(1.D0-KH)/SQRT(KH*(2.D0-KH)) C write(6,*) 'DKDKH', DKDKH DKBDK = (2.D0*P)*KA**(2.D0*P-1.D0)*(1.D0-(1.D0- & KA)**2.D0*KSIB**2.D0/A**2.D0) & + KA**(2.D0*P)*(2.D0*(1.D0-KA)*KSIB**2.D0/A**2.D0) C write(6,*) 'DKBDK', DKBDK DFDKH = -RCB**2.D0*Q**2.D0*DKBDK*DKDKH ETA=DFDKH IF (KH.GE.1.D0) THEN ETA = 0.D0 END IF C write(6,*) 'ETA', ETA C** fonction seuil F0=ROB * ROB - KB * RCB*RCB * Q*Q C**** Deuxième partie : calcul du flux m = DFDS C**** calcul de dF/dSigma C****** calcul de DROBDS DJ2DS(1) = S(1,1) DJ2DS(2) = S(2,2) DJ2DS(3) = S(3,3) C*** * 2 ? DJ2DS(4) = 2.D0*S(1,2) DJ2DS(5) = 2.D0*S(1,3) DJ2DS(6) = 2.D0*S(2,3) C** fin ? C*** DROB**2/DS DO I = 1,6 DROBDS(I) = 2.D0/(FC*FC)*DJ2DS(I) END DO C****** calcul de DRBCDS DO I = 1,6 DKSIBD(I) = 0.D0 END DO DO I = 1,3 DKSIBD(I) = 1.D0/(SQRT(3.D0)*FC) ENDDO DRCBDK = -GAMA*((-M+SQRT(M*M-12.D0*SQRT(3.D0)*M*KSIB+ & 36.D0*C)) /6.D0)**(GAMA-1.D0)*(SQRT(2.D0)*M)/ & (SQRT(M*M- 12.D0*SQRT(3.D0) *M*KSIB+36.D0*C)) DO I = 1,6 DRCBDS(I) = DRCBDK * DKSIBD(I) END DO C******* calcul de DQDS DB1DB0 = SQRT(3.D0)*(1.D0-ALFA)/(1.D0+B0)**2.D0 + & 2.D0*ALFA/(SQRT(3.D0)) DC0DB0 = 1.D0/(B1*(1.D0+B0)-SQRT(3.D0)*B0)**4.D0* & ((-SQRT(3.D0)*DB1DB0* & (2.D0*B0 -SQRT(3.D0)*B1)+(2.D0-SQRT(3.D0)*B1)* & (2.D0- SQRT(3.D0)*DB1DB0)) *(B1*(1.D0+B0)- & SQRT(3.D0)*B0)**2.D0 - & (2.D0-SQRT(3.D0)*B1) & * (2.D0*B0 - SQRT(3.D0)*B1) & * (2.D0*(B1*(1.D0+B0)-SQRT(3.D0)*B0) & *(DB1DB0*(1.D0+B0)+B1-SQRT(3.D0)))) DC1DB0 = -DC0DB0*(1.D0+B0)**2.D0-2.D0*C0*(1.D0+B0) DC2DB0 = 3.D0*DC0DB0*(1.D0-B0)**2.D0+3.D0*C0*2.D0* & (1.D0-B0)*(-1.D0) DC3DB0 = 2.D0*DC0DB0*SQRT(3.D0)*(1.D0-B0*B0)+ & 2.D0*C0*SQRT(3.D0)*(-2.D0*B0) DC4DB0 = (1.D0-B0*C0)+(1.D0+B0)*(-1.D0)*(C0+B0*DC0DB0) DC5DB0 = -(1.D0-3.D0*B0*C0)+(1.D0-B0)*(-3.D0)*(C0+B0*DC0DB0) DD0DB0 = DC1DB0 * COS(TETA)*COS(TETA)-DC2DB0*SIN(TETA)* & SIN(TETA) + DC3DB0*SIN(TETA)*COS(TETA) DD1DB0 = 2.D0*(DC4DB0*SQRT(3.D0)*COS(TETA)-DC5DB0*SIN(TETA)) DD2DB0 = (4.D0-3.D0*B0*C0)+B0*(-3.D0*(C0+B0*DC0DB0)) DQD0 = (4.D0*DZERO*D2/(2.D0*SQRT(D1*D1-4.D0*DZERO*D2))- & (D1-SQRT(D1*D1-4.D0*DZERO*D2))) & /(2.D0*DZERO*DZERO) DQD1 = 1.D0/(2.D0*DZERO)*(1.D0-(D1)/(SQRT(D1*D1-4.D0*DZERO*D2))) DQD2 = 1.D0/(SQRT(D1*D1-4.D0*DZERO*D2)) DQDB0 = DQD0 * DD0DB0 + DQD1*DD1DB0 + DQD2 * DD2DB0 DREBDK = -GAMA*((-M+SQRT(M*M-3.D0*SQRT(3.D0)*M*KSIB+ & 9.D0*C))/3.D0) & **(GAMA-1.D0) * (M/SQRT(2.D0))/SQRT(M*M- & 3.D0*SQRT(3.D0)*M*KSIB+9.D0*C) DB0DKS = 1.D0/(RCB*RCB)*(DREBDK*RCB-REB*DRCBDK) DO I =1,6 DB0DS(I) = DB0DKS*DKSIBD(I) END DO C****** eviter J2 = 0 (triaxial, définition des dérivées) PI3S6= XPI/6.D0 C write(6,*) 'absolu', ABS(ABS(TETA)-PI3S6) IF (J2 .GT. 0.D0) THEN IF (ABS(ABS(TETA)-PI3S6) .GT. 1.d-6) THEN DTEDJ2 = 3.D0*SQRT(3.D0)*J3/(4.D0*COS(3.D0*TETA))/ & (J2)**2.5 DTEDJ3 = - SQRT(3.D0)/(2.D0*COS(3.D0*TETA)*J2**1.5) ELSE DTEDJ2 =0.D0 DTEDJ3 = 0.D0 ENDIF ENDIF DO I = 1,3 DO J = 1,3 DINTE(I,J) = 0.D0 END DO END DO DO I = 1, 3 DO J = 1,I DO MM = 1,3 DINTE(I,J) = DINTE(I,J) + S(I,MM)*S(MM,J) END DO IF (I .EQ. J) THEN DINTE(I,J) = DINTE(I,J) - 2.D0/3.D0 * J2 END IF END DO END DO DJ3DS(1) = DINTE(1,1) DJ3DS(2) = DINTE(2,2) DJ3DS(3) = DINTE(3,3) C x2 ? C DJ3DS(4) = 2.D0*DINTE(2,1) DJ3DS(5) = 2.D0*DINTE(3,1) C DJ3DS(6) = DINTE(3,2) DJ3DS(6) = 2.D0*DINTE(3,2) DO I = 1,6 DTETDS(I) = DTEDJ2*DJ2DS(I)+DTEDJ3*DJ3DS(I) END DO D0DTET = -2.D0*C1*COS(TETA)*SIN(TETA)-2.D0*C2*SIN(TETA)* & COS(TETA)+C3*((COS(TETA))**2.D0-(SIN(TETA))**2.D0) D1DTET = 2.D0*(-C4*SQRT(3.D0)*SIN(TETA)-C5*COS(TETA)) D2DTET = 0.D0 DQDTET = DQD0 * D0DTET + DQD1*D1DTET + DQD2 * D2DTET DO I=1,6 DQDS(I) = DQDB0*DB0DS(I) + DQDTET * DTETDS(I) END DO C**** Calcul de DKBDKISB DKBDKS = (KA**(2.D0*P))*(-2.D0)*KSIB*(1.D0-KA)**2.D0 & / (A**2.D0) IF (KH.GE.1.D0) THEN DO I = 1,6 DFDS(I) = DROBDS(I) - 2.D0*DRCBDS(I)*RCB*KB*Q*Q & - 2.D0*DQDS(I)*Q*KB*RCB*RCB FM(I) = DFDS(I) END DO ELSE DO I = 1,6 DFDS(I) = DROBDS(I) - 2.D0*DRCBDS(I)*RCB*KB*Q*Q & - 2.D0*DQDS(I)*Q*KB*RCB*RCB - DKBDKS & * DKSIBD(I)*RCB*RCB*Q*Q FM(I) = DFDS(I) END DO END IF C write(6,*) 'DROBDS', DROBDS(4), DROBDS(5), DROBDS(6) C write(6,*) 'DQDS' , DQDS(4) , DQDS(5) , DQDS(6) C write(6,*) 'DTETDS', DTETDS(4), DTETDS(5), DTETDS(6) C DO I =1,6 C (6,*), 'DFDS', I, DFDS(I) C END DO END