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