API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Sep 8th, 2017
695
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Fortran 15.96 KB | None | 0 0
raw download clone embed print report
  1. *deck,usercreep       parallel  user                                    gal
  2.       SUBROUTINE usercreep (impflg, ldstep, isubst, matId , elemId,
  3.      &                      kDInPt, kLayer, kSecPt, nstatv, nprop,
  4.      &                      prop  , time  , dtime , temp  , dtemp ,
  5.      &                      toffst, Ustatev, creqv , pres  , seqv  ,
  6.      &                      delcr , dcrda)
  7. c*************************************************************************
  8. c
  9. c *** Implementation of API-579 Project Omega Creep Strain Formulas
  10. c     Written by Benjamin Turner, P.E.
  11. c     Revision: 7/18/2017
  12.  
  13. c     To-do: Add error checking for material parameter input
  14. c
  15. c
  16. c     input arguments
  17. c     ===============
  18. c      impflg   (in ,sc   ,i)             Explicit/implicit integration
  19. c                                         flag (currently not used)
  20. c      ldstep   (in ,sc   ,i)             Current load step
  21. c      isubst   (in ,sc   ,i)             Current sub step
  22. c      matId    (in ,sc   ,i)             number of material index
  23. c      elemId   (in ,sc   ,i)             Element number
  24. c      kDInPt   (in ,sc   ,i)             Material integration point
  25. c      kLayer   (in ,sc   ,i)             Layer number
  26. c      kSecPt   (in ,sc   ,i)             Section point
  27. c      nstatv   (in ,sc   ,i)             Number of state variables
  28. c      nprop    (in ,sc   ,i)             size of mat properties array
  29. c      prop     (dp ,ar(*),i)             mat properties array
  30. c      time                               Current time
  31. c      dtime                              Current time increment
  32. c      temp                               Current temperature
  33. c      dtemp                              Current temperature increment
  34. c      toffst   (dp, sc,   i)             temperature offset from absolute zero
  35. c      seqv     (dp ,sc  , i)             equivalent effective stress
  36. c      creqv    (dp ,sc  , i)             equivalent effective creep strain
  37. c      pres     (dp ,sc  , i)             hydrostatic pressure stress, -(Sxx+Syy+Szz)/3
  38. c
  39. c      Material Constants for Input -- from API-579
  40. c      prop(1)                            A0
  41. c      prop(2)                            A1
  42. c      prop(3)                            A2
  43. c      prop(4)                            A3
  44. c      prop(5)                            A4
  45. c      prop(6)                            B0
  46. c      prop(7)                            B1
  47. c      prop(8)                            B2
  48. c      prop(9)                            B3
  49. c      prop(10)                           B4
  50. c
  51. c      Constants -- from API-579
  52. c      prop(11)                           beta_omega: Prager factor = 0.33
  53. c      prop(12)                           d_cd: Creep ductility adjustment factor = 0
  54. c      prop(13)                           d_sr: Creep strain adjustment factor = 0
  55. c      Convergence criteria
  56. c      prop(14)                           ds_inc: Stres increment
  57. c      prop(15)                           de_inc: Creep increment
  58. c      Logging information
  59. c      prop(16)                           dolog: For dumping output: x>0 is on
  60. c      Minimum Temp/Stress
  61. c      prop(17)                           minT: The minimum temp to calc creep
  62. c      prop(18)                           minS: The minimum stress to calc creep
  63. c
  64. c     input output arguments              input desc     / output desc
  65. c     ======================              ==========       ===========
  66. c      Ustatev  (dp,ar(*), i/o)           user defined iinternal state variables at
  67. c                                         time 't' / 't+dt'.
  68. c                                         This array will be passed in containing the
  69. c                                         values of these variables at start of the
  70. c                                         time increment. They must be updated in this
  71. c                                         subroutine to their values at the end of
  72. c                                         time increment, if any of these internal
  73. c                                         state variables are associated with the
  74. c                                         creep behavior.
  75. c
  76. c      Ustatev(1)                         D_c    : damage factor: initialize to 0
  77. c      Ustatev(2)                         Ddot_c : damage factor rate: initialize to 0
  78. c      Ustatev(3)                         e_c    : creep strain: initialize to 0
  79. c      Ustatev(4)                         edot_c : creep strain: rate initialize to 0
  80. c
  81. c     output arguments
  82. c     ================
  83. c      delcr    (dp ,sc  , o)             incremental creep strain
  84. c      dcrda    (dp,ar(*), o)             output array
  85. c                                         dcrda(1) - derivative of incremental creep
  86. c                                                    strain to effective stress
  87. c                                         dcrda(2) - derivative of incremental creep
  88. c                                                    strain to creep strain
  89. c
  90. c
  91. c     local variables
  92. c     ===============
  93. c      A0,A1,A2,A3,A4 (dp, sc, l)           Material Constants
  94. c      B0,B1,B2,B3,B4 (dp, sc, l)           Material Constants
  95. c      beta           (dp, sc, l)           Prager constant
  96. c      d_cd,d_sr      (dp, sc, l)           Adjustment factors
  97. c      omega          (dp, sc, l)           Omega uniaxial damage parameter
  98. c      lomega         (dp, sc, l)           Log of omega uniaxial damage parameter
  99. c      omega_m        (dp, sc, l)           Omega multiaxial damage parameter
  100. c      omega_n        (dp, sc, l)           Minimum Omega uniaxial damage parameter
  101. c      S_l            (dp, sc, l)           Stress parameter
  102. c      S_l2, S_l3     (dp, sc, l)           Stress parameter (exponentiated)
  103. c      s_1, s_2, s_3  (dp, sc, l)           Principal stresses
  104. c      s_e            (dp, sc, l)           Effective stresses
  105. c      D_c            (dp, sc, l)           Cumulative damage
  106. c      Ddot_c         (dp, sc, l)           Cumulative damage rate
  107. c      e_c            (dp, sc, l)           Creep strain
  108. c      edot_c         (dp, sc, l)           Creep strain rate
  109. c      edot_oc        (dp, sc, l)           Creep strain parameter
  110. c      ledot_oc       (dp, sc, l)           Log of creep strain parameter
  111. c      n              (dp, sc, l)           Intermediate value
  112. c      delta          (dp, sc, l)           Intermediate stress value
  113. c      T              (dp, sc, l)           Temperature including offset
  114. c      de_inc         (dp, sc, l)           Strain increment for differentiation
  115. c      ds_inc         (dp, sc, l)           Stress increment for differentiation
  116. c      s_vec          (dp, sc, l)           Stress vector from element
  117. c      ncomp          (dp, sc, l)           Number of stress components
  118. c      delcr_e        (dp, sc, l)           Creep strain wrt strain
  119. c      delcr_s        (dp, sc, l)           Creep strain wrt stress
  120. c      dolog          (dp, sc, l)           Any integer > 0 turns on logging output
  121. c      minT           (dp, sc, l)           Minimum temp (F) to calc creep
  122. c      minS           (dp, sc, l)           Min triaxial stress (ksi) to calc creep
  123. c
  124. c*************************************************************************
  125. c
  126. c --- parameters
  127. c
  128. #include "impcom.inc"
  129. #include "ansysdef.inc"
  130.  
  131.       EXTERNAL         get_ElmInfo,
  132.      &                 get_ElmData,
  133.      &                 prinst,
  134.      &                 wrinqr
  135.  
  136.       DOUBLE PRECISION ZERO
  137.       PARAMETER        (ZERO = 0.0d0)
  138. c
  139. c --- argument list
  140. c
  141.  
  142.       INTEGER          ldstep, isubst, matId , elemId,
  143.      &                 kDInPt, kLayer, kSecPt, nstatv,
  144.      &                 impflg, nprop
  145.       DOUBLE PRECISION dtime , time  , temp  , dtemp , toffst,
  146.      &                 creqv , seqv  , pres, delcr
  147.       DOUBLE PRECISION prop(*), dcrda(*), Ustatev(nstatv)
  148.  
  149. c -- for debug output
  150.       INTEGER          wrinqr, iott, ntrg
  151.  
  152. c
  153. c --- local variables
  154. c
  155.       DOUBLE PRECISION A0    , A1    , A2     , A3   , A4   ,
  156.      &                 B0    , B1    , B2     , B3   , B4   ,
  157.      &                 beta  ,
  158.      &                 d_cd  , d_sr  ,
  159.      &                 omega , omega_m, omega_n,
  160.      &                 lomega,
  161.      &                 edot_c, edot_oc,
  162.      &                 ledot_oc,
  163.      &                 S_l   ,
  164.      &                 S_l2  , S_l3,
  165.      &                 s_1   , s_2   , s_3    ,
  166.      &                 s_e   ,
  167.      &                 D_c   , Ddot_c,
  168.      &                 e_c   ,
  169.      &                 n     , delta ,
  170.      &                 T     ,
  171.      &                 de_inc, ds_inc,
  172.      &                 delcr_e, delcr_s,
  173.      &                 minT, minS
  174.       INTEGER          ncomp , dolog
  175.  
  176. c --- for the stress vector
  177.       INTEGER          MSVAR
  178.       parameter (MSVAR = 11)
  179.       DOUBLE PRECISION s_vec(MSVAR)
  180. c
  181. c*************************************************************************
  182. c
  183. c *** skip when stress and creep strain are all zero
  184.       if (seqv .le. ZERO .and. creqv .le. ZERO) go to 990
  185. c *** API-579 Material Constants
  186.       A0      = prop(1)
  187.       A1      = prop(2)
  188.       A2      = prop(3)
  189.       A3      = prop(4)
  190.       A4      = prop(5)
  191.       B0      = prop(6)
  192.       B1      = prop(7)
  193.       B2      = prop(8)
  194.       B3      = prop(9)
  195.       B4      = prop(10)
  196. c      prop(11)                           beta_omega: Prager factor = 0.33
  197. c      prop(12)                           d_cd: Creep ductility adjustment factor = 0
  198. c      prop(13)                           d_sr: Creep strain adjustment factor = 0
  199.       beta    = prop(11)
  200.       d_cd = prop(12)
  201.       d_sr = prop(13)
  202.  
  203. c *** Numerical differentiation (may need to adjust)
  204. c     ds_inc is delta(stress) while de_inc is delta(strain)
  205. c      ds_inc   = 2.0d0**(-5.0d0)
  206. c      de_inc   = 2.0d0**(-15.0d0)
  207.       ds_inc   = prop(14)
  208.       de_inc   = prop(15)
  209. c *** Is logging turned on?
  210.       dolog    = prop(16)
  211.  
  212. c *** Min temp and min stress
  213.       minT    = prop(17)
  214.       minS    = prop(18)
  215.  
  216. c *** Absolute temperature
  217.       T        = temp + toffst
  218.       if(T .le. ZERO) go to 990
  219.  
  220. c *** Min temp
  221.       if(temp .le. minT) go to 2000
  222.  
  223. c *** Get number of stress components
  224.       call get_ElmInfo ('NCOMP',ncomp)
  225.  
  226. c *** Get stress vector
  227.       call get_ElmData ('SIG ',elemId, kDInPt, ncomp, s_vec)
  228.  
  229. c *** Get principal stresses -- Maybe need to put in storage var first??
  230.       call prinst(s_vec(1))
  231. c    primary function:  computes principal stresses and stress intensity
  232. c    secondary functions:  none
  233. c *** Notice - This file contains ANSYS Confidential information ***
  234.  
  235. c  input arguments:
  236. c     variable (typ,siz,intent)    description
  237. c     s        (dp,ar(11),inout)   - stress vector
  238. c                  s(1)=sx
  239. c                  s(2)=sy
  240. c                  s(3)=sz
  241. c                  s(4)=sigxy
  242. c                  s(5)=sigyz
  243. c                  s(6)=sigzx
  244.  
  245. c
  246. c  output arguments:
  247. c     variable (typ,siz,intent)    description
  248. c     s        (dp,ar(11),inout)   - stress vector
  249. c                  s(7)=sig1
  250. c                  s(8)=sig2
  251. c                  s(9)=sig3
  252. c                  s(10)=s.i.
  253. c                  s(11)=sige
  254.  
  255. c ********* note: all changes to this routine must be made in
  256. c                 post1  (paprst)
  257. c
  258. c *** Save principal stresses and equivalent
  259. c *** Make sure they're in KSI! KSI!
  260.  
  261.       s_1 = s_vec(7) / 1000
  262.       s_2 = s_vec(8) / 1000
  263.       s_3 = s_vec(9) / 1000
  264.       s_e = seqv / 1000
  265.  
  266. c *** Min triaxial stress
  267.       if(s_1+s_2+s_3 .le. minS) go to 2000
  268.  
  269. c *** Compute intermediate values
  270.  
  271.       delta = beta * (((s_1+s_2+s_3)/s_e) - 1)
  272.       S_l = LOG10 (s_e)
  273.       S_l2 = S_l**2
  274.       S_l3 = S_l**3
  275.  
  276. c *** Compute omega axial damage parameters
  277.  
  278.       lomega = ((B0+d_cd) + ((B1 + B2*S_l + B3*S_l2 + B4*S_l3)/T))
  279.       omega = 10**lomega
  280.       n = -1 * (A2 + 2*A3*S_l + 3*A4*S_l2)/T
  281.       omega_n = omega - n
  282.       if ( omega_n .le. 3 ) omega_n = 3
  283.       omega_m = omega_n ** (delta + 1)
  284.  
  285. c *** Compute intermediate creep strain rate parameters
  286.  
  287.       ledot_oc = -1* ((A0+d_sr) + ((A1 + A2*S_l + A3*S_l2 + A4*S_l3)/T))
  288.       edot_oc = 10**ledot_oc
  289.  
  290. c *** Compute creep strain rate
  291.       e_c = creqv
  292. c     if(e_c .le. TINY) e_c = sqrt(TINY)
  293.       edot_c = edot_oc*exp(omega_m*e_c)
  294.  
  295. c *** Compute incremental creep strain
  296.       delcr = edot_c * dtime
  297.  
  298. c *** Compute damage factors
  299.       D_c = Ustatev(1)
  300.       Ddot_c = omega_m * edot_oc
  301.  
  302. c *** Update damage factors
  303.       Ustatev(1) = D_c + Ddot_c*dtime
  304.       Ustatev(2) = Ddot_c
  305.       Ustatev(3) = creqv
  306.       Ustatev(4) = edot_c
  307.  
  308. c *** DONE DONE DONE DONE DONE WITH FIRST PASS
  309. c *** NOW CALCULATE NUMERICAL DERIVATIVES
  310. c *** DONE DONE DONE DONE DONE WITH FIRST PASS
  311. c *** Compute relative to effective strain
  312. c *** Compute creep strain rate (dstrain)
  313.       e_c = creqv+de_inc
  314. c      if(e_c .le. TINY) e_c = sqrt(TINY)
  315.       edot_c = edot_oc*exp(omega_m*e_c)
  316.  
  317. c *** Compute incremental creep strain (dstrain)
  318.       delcr_e = edot_c * dtime
  319.  
  320. c *** DONE DONE DONE DONE DONE WITH SECOND PASS
  321. c *** DONE DONE DONE DONE DONE WITH SECOND PASS
  322.  
  323. c *** Compute relative to effective stress
  324.       s_e = (seqv+ds_inc)/1000
  325.  
  326. c *** Compute intermediate values (dstress)
  327.  
  328.       delta = beta * (((s_1+s_2+s_3)/s_e) - 1)
  329.       S_l = LOG10 (s_e)
  330.       S_l2 = S_l**2
  331.       S_l3 = S_l**3
  332.  
  333. c *** Compute omega axial damage parameters (dstress)
  334.  
  335.       lomega = ((B0+d_cd) + ((B1 + B2*S_l + B3*S_l2 + B4*S_l3)/T))
  336.       omega = 10**lomega
  337.       n = -1 * (A2 + 2*A3*S_l + 3*A4*S_l2)/T
  338.       omega_n = omega - n
  339.       if ( omega_n .le. 3 ) omega_n = 3
  340.       omega_m = omega_n ** (delta + 1)
  341.  
  342. c *** Compute intermediate creep strain rate parameters (dstress)
  343.  
  344.       ledot_oc = -1* ((A0+d_sr) + ((A1 + A2*S_l + A3*S_l2 + A4*S_l3)/T))
  345.       edot_oc = 10**ledot_oc
  346.  
  347. c *** Compute creep strain rate (dstress)
  348.       e_c = creqv
  349. c      if(e_c .le. TINY) e_c = sqrt(TINY)
  350.       edot_c = edot_oc*exp(omega_m*e_c)
  351.  
  352. c *** Compute incremental creep strain (dstress)
  353.       delcr_s = edot_c * dtime
  354. c *** Derivative of incremental creep strain to effective stress
  355.       dcrda(1) = (delcr_s - delcr) / ds_inc
  356. c *** derivative of incremental creep strain to effective creep strain
  357.       dcrda(2) = (delcr_e - delcr) / de_inc
  358.  
  359. c *** save incremental creep strain as variable 2
  360.       Ustatev(4) = delcr
  361.  
  362. c *** trigger if we find the nodes
  363.       ntrg = 0
  364.       if (elemId .eq. 1337) then
  365.         ntrg = 1
  366.       elseif (elemId .eq. 1336) then
  367.         ntrg = 1
  368.       elseif (elemId .eq. 3797) then
  369.         ntrg = 1
  370.       else
  371.         ntrg = 0
  372.       end if
  373.  
  374. c *** debug output because of failing node Convergence
  375.  
  376.       if (dolog .gt. 0 .and. ntrg .gt. 0) then
  377.  
  378.         iott = wrinqr(WR_OUTPUT)
  379.         write(iott, 1000) 'ElemID: ', elemId, ' MatID: ', matId,
  380.      & ' Temp: ', T
  381. 1000  format(A,I,A,I,A,F)
  382.  
  383.         write(iott, 1010), 'A0: ', A0, ' A1: ', A1, ' A2: ', A2,
  384.      & ' A3: ', A3, ' A4: ', A4, ' B0: ', B0, ' B1: ', B1, ' B2: ',
  385.      & B2, ' B3: ', B3, ' B4: ', B4
  386. 1010  format(10(A, F))
  387.  
  388.         write(iott, 1020), 'Beta: ', beta, ' d_cd: ', d_cd, ' d_sr: ',
  389.      & d_sr, ' omega: ', omega, ' omega_m: ', omega_m, ' omega_n: ',
  390.      & omega_n, ' n: ', n, ' lomega: ', lomega
  391. 1020  format(8(A, F))
  392.  
  393.        write(iott, 1030), 'S_l: ', S_l, ' S_l2: ', S_l2, ' S_l3: ',
  394.      & S_l3, ' s_1: ', s_1, ' s_2: ', s_2, ' s_3: ', s_3, ' s_e: ',
  395.      & s_e, ' D_c: ', D_c, ' Ddot_c: ', Ddot_c, ' delta: ', delta
  396. 1030  format(10(A, F))
  397.  
  398.        write(iott, 1040), 'edot_c: ', edot_c, ' edot_oc: ', edot_oc,
  399.      & ' ledot_oc: ', ledot_oc, ' e_c: ', e_c, ' de_inc: ', de_inc,
  400.      & ' ds_inc: ', ds_inc, ' delcr_e: ', delcr_e, ' delcr_s: ',
  401.      & delcr_s, ' dcrda(1): ', dcrda(1), ' dcrda(2) ', dcrda(2)
  402. 1040  format(10(A, F))
  403.       end if
  404.  
  405.  990  continue
  406.       return
  407.  
  408. c *** Set all the creep and such to zero because we don't need to calculate it
  409. c *** Set incremental creep strain to zero
  410. 2000  delcr = 0
  411. c *** Set derivatives to zero
  412.       dcrda(1) = 0
  413.       dcrda(2) = 0
  414.       continue
  415.       return
  416.       end
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!