Télécharger la plaquette Cast3M


fr en

Télécharger flowm.eso

Retour à la liste

Numérotation des lignes : 

flowm

  1. C FLOWM     SOURCE    CB215821  16/04/21    21:16:53     8920
  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 * ASIN(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.  
  434.  
  435.  
  436.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales