Télécharger ddsirank.eso

Retour à la liste

Numérotation des lignes :

ddsirank
  1. C DDSIRANK SOURCE CB215821 16/04/21 21:16:19 8920
  2. SUBROUTINE ddsigRank(sigma,d2gt,i1,i2,i3,i4,i5,i6,lcp,
  3. . parahot3,idimpara3)
  4.  
  5. c This subroutine calculates the second derivative of the Rankine
  6. c plastic function with respect to sigma
  7. c gt = ft = S1 = largest principal stress
  8.  
  9. IMPLICIT REAL*8 (A-B,D-H,O-Z)
  10. implicit integer (I-K,M,N)
  11. implicit logical (L)
  12. implicit character*10 (C)
  13.  
  14. **** dimension sigma(i6),dgt(i6),d2gt(i6,i6)
  15. dimension sigma(i6),dgt(6),d2gt(i6,i6)
  16. **** dimension rcossigmapr(i3,i3),P1(3,3),pi1(3)
  17. dimension rcossigmapr(3,3),P1(3,3),pi1(3)
  18. dimension sigmacp(3),vloc(3),vloc2(1,3),vloc1(1)
  19. dimension vloc3(3,3),d2gtcp(3,3),parahot3(idimpara3)
  20. dimension sigmaim(6),sigmail(6),rhi(6),dgtm(6),dgtl(6)
  21.  
  22.  
  23.  
  24. *****
  25. * MESSAGES D'ERREUR ( SUPPRESSION DES AUTOMATIC OBJECTS)
  26. IF(I3.GT.3) PRINT *, ' DDSIGRANK - ERREUR I2 = ', I2, ' > 2 '
  27. IF(I6.GT.6) PRINT *, ' DDSIGRANK - ERREUR I6 = ', I6, ' > 6 '
  28. *****
  29.  
  30.  
  31. if (lcp) then
  32. c In plane stress, the second derivative can be calculated analytically
  33. do jloc=1,6
  34. do iloc=1,6
  35. d2gt(iloc,jloc) = 0.d0
  36. enddo
  37. enddo
  38. P1(1,1) = 0.5d0
  39. P1(2,2) = 0.5d0
  40. P1(1,2) = -0.5d0
  41. P1(2,1) = -0.5d0
  42. do iloc=1,2
  43. P1(iloc,3) = 0.d0
  44. P1(3,iloc) = 0.d0
  45. enddo
  46. P1(3,3) = 2.d0
  47. pi1(1) = 1.d0
  48. pi1(2) = 1.d0
  49. pi1(3) = 0.d0
  50. sigmacp(1) = sigma(1)
  51. sigmacp(2) = sigma(2)
  52. sigmacp(3) = sigma(4)
  53. call mulAB(P1,sigmacp,vloc,3,3,3,1)
  54. call mulATB(sigmacp,vloc,vloc1,3,1,3,1)
  55. rloc = 0.5d0 * vloc1(1)
  56. psi1 = SQRT(rloc)
  57. call mulATB(sigmacp,P1,vloc2,3,1,3,3)
  58. call mulAB(vloc,vloc2,vloc3,3,1,1,3)
  59.  
  60. test = 0.5d0*(sigmacp(1)+sigmacp(2))
  61. psi1lim=max(1.E-7*test,1.d-3)
  62. if (ABS(psi1).le.psi1lim) then
  63. continue
  64. c d2gt = 0
  65. else
  66. do jloc=1,3
  67. do iloc=1,3
  68. d2gtcp(iloc,jloc)=(P1(iloc,jloc)/(2.0d0*psi1))
  69. . - (vloc3(iloc,jloc)/(4.0d0*psi1*psi1*psi1))
  70. enddo
  71. enddo
  72. do iloc=1,2
  73. d2gt(1,iloc) = d2gtcp(1,iloc)
  74. d2gt(2,iloc) = d2gtcp(2,iloc)
  75. d2gt(4,iloc) = d2gtcp(3,iloc)
  76. d2gt(iloc,4) = d2gtcp(iloc,3)
  77. end do
  78. d2gt(4,4) = d2gtcp(3,3)
  79. endif
  80. else
  81. c In full 3D, numerical differentiation
  82. do jloc=1,6
  83. do iloc=1,6
  84. d2gt(iloc,jloc) = 0.d0
  85. enddo
  86. enddo
  87. rh = 1.0d-5
  88. do iloc=1,6
  89. c ith component of sigma => ith column of d2gt
  90. c calculation of sigma + hi * ei and sigma - hi * ei
  91. do jloc=1,6
  92. rhi(jloc) = 0.d0
  93. enddo
  94. rhi(iloc) = 1.0d0 * rh
  95. do jloc=1,6
  96. sigmaim(jloc) = sigma(jloc) + rhi(jloc)
  97. sigmail(jloc) = sigma(jloc) - rhi(jloc)
  98. enddo
  99. c Calculation of q(x+hi*ei) = dft(x+hi*ei)/dsig
  100. call dsigRank(sigmaim,dgtm,3,6,parahot3,idimpara3,
  101. . rkappait,lcp)
  102. call dsigRank(sigmail,dgtl,3,6,parahot3,idimpara3,
  103. . rkappait,lcp)
  104. do jloc=1,6
  105. d2gt(jloc,iloc) = (dgtm(jloc) - dgtl(jloc))/(2.0d0*rh)
  106. enddo
  107. enddo
  108. c Calculation of f(x+hi*ei)
  109. endif
  110.  
  111. RETURN
  112. END
  113.  
  114.  
  115.  
  116.  
  117.  

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