Télécharger gernay.eso

Retour à la liste

Numérotation des lignes :

  1. C GERNAY SOURCE PV 17/12/08 21:17:25 9660
  2.  
  3. SUBROUTINE GERNAY(WRK52,WRK53,WRK54,WRK3,IB,IGAU,NBPGAU)
  4. *
  5. *______________________________________________________________________
  6. *
  7. * MODELE GERNAY
  8. *______________________________________________________________________
  9. *
  10. IMPLICIT INTEGER(I-N)
  11. IMPLICIT REAL*8(A-H,O-Z)
  12.  
  13.  
  14. -INC PPARAM
  15. -INC CCOPTIO
  16. -INC DECHE
  17. *
  18. SEGMENT WRK3
  19. REAL*8 WORK(LW),WORK2(LW2)
  20. ENDSEGMENT
  21. *
  22. PARAMETER(UN=1.D0,DEUX=2.D0,UNDEMI=.5D0)
  23. *
  24. REAL*8 YOUN, NU, BULK3, G2, LAM,
  25. & TETA0, DTETA, TETAF, ALPH, TREF, DYDT, DNDT,
  26. & YOUN0, NU0, BULK30, G20, LAM0,
  27. & YOUNF, NUF, BULK3F, G2F, LAMF,
  28. & DELDSE(6,6), ETHERM(6), DTHERM(6),
  29. & YREF, NREF, EPSE0(6), DEPSE(6),DHOOK(6,6),
  30. & DSIG2(3,3),DEPST2(3,3),TENS(3,3),DSIG(6),
  31. & V1(3),V2(3),V3(3),V(3,3),P(3,3)
  32.  
  33. REAL*8 PARAHOT3(83,1),sigma6v(6), deps6v(6),STRESS(6),DSTRAN(6)
  34. REAL*8 STRAN(6),STATEV(32),PROPS(10),DDSDDE(6,6)
  35. character*10 cmat
  36. logical lerror,lcp
  37. *
  38. * ON EXPLORE LE DECHE
  39. *
  40. NPT=1
  41. NOEL=1
  42. KERRE=0
  43. NMAT = XMAT(/1)
  44.  
  45. NVARI = VAR0(/1)
  46. NSTRS = SIG0(/1)
  47.  
  48. CALL ZERO(STRESS,6,1)
  49. CALL ZERO(STRAN,6,1)
  50. CALL ZERO(DSTRAN,6,1)
  51. CALL ZERO(STATEV,32,1)
  52. CALL ZERO(PARAHOT3,83,1)
  53. *
  54. * TRANSFER DE L'ETAT INITIAL, DU DEPST, ET DU MATERIAU
  55. *
  56. DO I=1,NSTRS
  57. STRESS(I) = SIG0(I)
  58. STRAN(I) = EPST0(I)
  59. DSTRAN(I)= DEPST(I)
  60. ENDDO
  61. *
  62. DO I=1,NVARI
  63. STATEV(I) = VAR0(I)
  64. ENDDO
  65. *
  66. PROPS( 1) = XMAT( 2)
  67. PROPS( 2) = XMAT( 5)
  68. PROPS( 3) = XMAT( 6)
  69. PROPS( 4) = XMAT( 7)
  70. PROPS( 5) = XMAT( 8)
  71. PROPS( 6) = XMAT( 9)
  72. PROPS( 7) = XMAT(10)
  73. PROPS( 8) = XMAT(11)
  74. PROPS( 9) = XMAT(12)
  75. PROPS(10) = XMAT( 1)
  76.  
  77. C
  78. C-----------------------------------------------------------------------
  79. C Constantes de materiau
  80. C ----------------------
  81. C PROPS(1) : NU Coefficient de Poisson
  82. c PROPS(2) : RC Resistance à la compression
  83. c PROPS(3) : RT Resistance à la traction
  84. c PROPS(4) : RB Resistance à la compression biaxiale
  85. c PROPS(5) : EPS1 Deformation au pic en compression
  86. c PROPS(6) : DILA Parametre de dilatance alpha_g
  87. c PROPS(7) : DAM1 Endommagement de compression au pic
  88. c PROPS(8) : XC Paramètre adimensionnel, rapport de l'energie de rupture
  89. c en compression dissipée avant le pic sur l'énergie totale
  90. c PROPS(9) : GT Densite d'energie de rupture en traction
  91. c PROPS(10) : YOUN Young pour le calcul de la rigidite elastique
  92. c
  93. C-----------------------------------------------------------------------
  94. C
  95.  
  96. * print *,' '
  97. * print *, ' POINT DINTEGRATION IGAU = ', IGAU
  98. * print *, ' stress', stress
  99. * print *, 'dstran', dstran
  100. * print *,' STATEV 1 à 4' , (statev(i),i=1,4)
  101. * print *,' STATEV 5 à 10 epspl ', (statev(i),i=5,10)
  102. * print *,' STATEV 11 à 16 epstr ', (statev(i),i=11,16)
  103. * print *,' STATEV 17 à 22 sig,eff,c ', (statev(i),i=17,22)
  104. * print *,' STATEV 24 à 29 sig,eff ', (statev(i),i=24,29)
  105. * print *,' STATEV 30 à 32 tmax,phi,epsth ',(statev(i),i=30,32)
  106.  
  107. cmat = 'SICO_ETC3D'
  108. idimpara3 = 83
  109. c ATTENTION si idimpara3 est modifié il faut modifier aussi la dimension
  110. c de parahot3 (plus haut)
  111. lerror = .false.
  112. **** npttot3=30
  113. npttot3=1
  114. *
  115. lcp = .false.
  116. *
  117. * CAS CONTRAINTES PLANES
  118. *
  119. if (IFOUR.EQ.-2) then
  120. lcp = .true.
  121. endif
  122.  
  123.  
  124. * TEMPERATURES EN DEBUT ET FIN DE PAS
  125.  
  126. TETA1 = ture0(1)
  127. TETA2 = turef(1)
  128.  
  129. C Recuperation des constantes de materiau
  130. C
  131. PARAHOT3(2,NPT) = PROPS(1)
  132. PARAHOT3(3,NPT) = PROPS(2)
  133. PARAHOT3(4,NPT) = PROPS(3)
  134. PARAHOT3(7,NPT) = PROPS(4)
  135. PARAHOT3(5,NPT) = PROPS(5)
  136. PARAHOT3(8,NPT) = PROPS(6)
  137. PARAHOT3(9,NPT) = PROPS(7)
  138. PARAHOT3(10,NPT)= PROPS(8)
  139. PARAHOT3(11,NPT)= PROPS(9)
  140. C
  141. c remplissage des variables internes
  142. c accumulated plastic strain kappa at the beginning of the step
  143. c in tension and in compression
  144. PARAHOT3(idimpara3-37,NPT) = STATEV(1)
  145. PARAHOT3(idimpara3-36,NPT) = STATEV(2)
  146.  
  147. c damage variable in tension and in compression
  148. PARAHOT3(idimpara3-41,NPT) = STATEV(3)
  149. PARAHOT3(idimpara3-40,NPT) = STATEV(4)
  150.  
  151. c plastic strains
  152. do i=1,6
  153. PARAHOT3(idimpara3-25+i,NPT) = STATEV(4+i)
  154. end do
  155.  
  156. c transient creep strain at the beginning of the step
  157. do i=1,6
  158. PARAHOT3(17+i,NPT) = STATEV(10+i)
  159. end do
  160.  
  161. c effective stress in compression (used for transient creep strain calculation)
  162. do i=1,6
  163. PARAHOT3(23+i,NPT) = STATEV(16+i)
  164. end do
  165. do i=1,6
  166. PARAHOT3(29+i,NPT) = PARAHOT3(23+i,NPT)
  167. end do
  168.  
  169. c effective stress after convergence (used for sigf3d forward calculation)
  170. do i=1,6
  171. PARAHOT3(idimpara3-7+i,NPT) = STATEV(23+i)
  172. end do
  173.  
  174. c remplissage des contraintes
  175. c remplissage du vecteur contraintes et du vecteur
  176. c increment de deformation
  177. if (.not.lcp) then
  178. do i=1,6
  179. PARAHOT3(idimpara3-13+i,NPT) = STRESS(i)
  180. sigma6v(i)= STRESS(i)
  181. deps6v(i)= DSTRAN(i)
  182. end do
  183. else
  184. do i=1,4
  185. PARAHOT3(idimpara3-13+i,NPT) = STRESS(i)
  186. sigma6v(i)= STRESS(i)
  187. deps6v(i)= DSTRAN(i)
  188. end do
  189. PARAHOT3(idimpara3-8,NPT)=0.d0
  190. PARAHOT3(idimpara3-7,NPT)=0.d0
  191. sigma6v(5) = 0.d0
  192. sigma6v(6) = 0.d0
  193. deps6v(5) = 0.d0
  194. deps6v(6) = 0.d0
  195. endif
  196.  
  197. c remplissage des deformations au debut du pas t
  198. c et en fin de pas
  199. if (.not.lcp) then
  200. do i=1,6
  201. PARAHOT3(idimpara3-31+i,NPT) = STRAN(i)
  202. PARAHOT3(idimpara3-19+i,NPT) = STRAN(i) + DSTRAN(i)
  203. end do
  204. else
  205. do i=1,4
  206. PARAHOT3(idimpara3-31+i,NPT) = STRAN(i)
  207. PARAHOT3(idimpara3-19+i,NPT) = STRAN(i) + DSTRAN(i)
  208. end do
  209. endif
  210.  
  211. c TEMPERATURES
  212.  
  213. c ********************************************************************************
  214. c probleme pour vecteur STATEV - STATEV(30,31,32) ne passe pas d'un pas à l'autre
  215. c => il faut passer le vecteur STATEV à 32 composantes au lieu de 29
  216. c ********************************************************************************
  217. c actualisation de la temperature au pas courant
  218. **** TEMP = TEMP + DTEMP
  219. **** PARAHOT3(idimpara3-39,NPT) = TEMP
  220. PARAHOT3(idimpara3-39,NPT) = TETA2
  221. c actualisation éventuelle de la temperature maximale atteinte
  222. **** if (TEMP.gt.STATEV(30)) STATEV(30) = TEMP
  223. if (TETA2.gt.STATEV(30)) STATEV(30) = TETA2
  224. PARAHOT3(idimpara3-38,NPT) = STATEV(30)
  225.  
  226. c phi function
  227. PARAHOT3(17,NPT) = STATEV(31)
  228.  
  229. c Thermal strain
  230. PARAHOT3(idimpara3,NPT) = STATEV(32)
  231.  
  232. c appel a mater2
  233. NTOT = 1
  234.  
  235. CALL MATER2TG (cmat,NTOT, lerror,idimpara3, parahot3, npttot3,
  236. & ddsdde, sigma6v, deps6v,lcp,NOEL,NPT,Kinc)
  237. * print *, ' sortie de mater2tg'
  238.  
  239. CALL MATER3TG (cmat,npt,idimpara3,parahot3,npttot3,lcp,
  240. . NOEL,NPT)
  241. * print *, ' sortie de mater3tg'
  242.  
  243. C sauvetage des resultats
  244.  
  245. c remplissage des contraintes
  246. if (.not.lcp) then
  247. do i=1,6
  248. STRESS(i) = sigma6v(i)
  249. end do
  250. else
  251. do i=1,4
  252. STRESS(i) = sigma6v(i)
  253. end do
  254. endif
  255. c if ((NOEL.eq.18).and.(NPT.eq.23)) then
  256. * print *,'GERNAY CONTRAINTES EN FIN DE PAS stress = ', stress
  257. c endif
  258.  
  259. c actualisation des variables internes
  260. c accumulated plastic strain kappa at the beginning of the step
  261. c in tension and in compression
  262. STATEV(1) = PARAHOT3(idimpara3-37,NPT)
  263. STATEV(2) = PARAHOT3(idimpara3-36,NPT)
  264.  
  265. c damage variable in tension and in compression
  266. STATEV(3) = PARAHOT3(idimpara3-41,NPT)
  267. STATEV(4) = PARAHOT3(idimpara3-40,NPT)
  268.  
  269. c plastic strains
  270. do i=1,6
  271. STATEV(4+i) = PARAHOT3(idimpara3-25+i,NPT)
  272. end do
  273.  
  274. c transient creep strain at the beginning of the step
  275. do i=1,6
  276. STATEV(10+i) = PARAHOT3(17+i,NPT)
  277. end do
  278.  
  279. c effective stress in compression (used for transient creep strain calculation)
  280. do i=1,6
  281. STATEV(16+i) = PARAHOT3(23+i,NPT)
  282. end do
  283.  
  284. c equivalent strain
  285. STATEV(23) = STATEV(1) + STATEV(2)
  286.  
  287. c effective stress after convergence (used for sigf3d forward calculation)
  288. do i=1,6
  289. STATEV(23+i) = PARAHOT3(idimpara3-7+i,NPT)
  290. end do
  291.  
  292. c maximum temperature
  293. STATEV(30) = PARAHOT3(idimpara3-38,NPT)
  294.  
  295. c phi function
  296. STATEV(31) = PARAHOT3(17,NPT)
  297.  
  298. c Thermal strain
  299. STATEV(32) = PARAHOT3(idimpara3,NPT)
  300.  
  301. c do iloc=1,NTENS
  302. do iloc=1,4
  303. test = ABS(stress(iloc))
  304. if (test.gt.10.d20) then
  305. print *,'>>>>>>>>>>>>>>>>>>>> problem stress = ',
  306. . stress
  307. print *,'in NOEL = , NPT = ', IB , IGAU
  308. KINC = -2
  309. endif
  310. enddo
  311. if (KINC.eq.-2) then
  312. print *,'in NOEL = , NPT = ', IB , IGAU
  313. endif
  314. *
  315. * TRANSFERT DE L'ETAT FINAL
  316. *
  317. DO I=1,NSTRS
  318. SIGF(I) = STRESS(I)
  319. ENDDO
  320. *
  321. DO I=1,NVARI
  322. VARF(I) = STATEV(I)
  323. ENDDO
  324. *
  325. RETURN
  326. END
  327.  
  328.  
  329.  
  330.  
  331.  
  332.  
  333.  
  334.  

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