Télécharger mater2tg.eso

Retour à la liste

Numérotation des lignes :

  1. C MATER2TG SOURCE CB215821 16/04/21 21:17:44 8920
  2. SUBROUTINE mater2tg(cmat,ntot,lerror,idimpara3,
  3. * parahot3,npttot3,H66t,sigma6v,deps6v,lcp,
  4. * NOEL,NPT,Kinc)
  5.  
  6. c cmat is the name (characters) of the material
  7. c ntot is the integration point (global number for the whole structure): counter
  8. c idimpara3 is the number of parameters for parahot3
  9. c parahot3 is the array that contains the material parameters and state variables
  10. c npttot3 is the number of integration points (number for the considered material 3D)
  11. c and is thus the number of columns in parahot3
  12. c H66t is the tangent matrix
  13. c sigma6v is the nominal stress vector
  14. c deps6v is the incremental strains vector
  15.  
  16. c===========================================================================
  17.  
  18. c input
  19. c CMAT : name of the material
  20. c NTOT : # of the integration point ( global # for the entire structure )
  21. c can be NTOT1 or NTOT2 or NTOT3
  22.  
  23. c output
  24. c H66t : tangent stiffness matrix for 3D materials (THG add)
  25.  
  26. c input and output
  27. c PARAHOT3 : the first columns contain the material properties at elevated temperature,
  28. c the last columns contain strains, stresses, ....
  29.  
  30. c SIGMA6v : stress vector for 3D materials (THG add)
  31. c DEPS6V : incremental strain vector for 3D materials (THG add)
  32.  
  33. c=======================================================================
  34.  
  35.  
  36. c This routine is run for every iteration,
  37. c for every point of integration.
  38. c
  39. C At the integration point NTOT is the material described by CMAT.
  40. c This material has, at the point NTOT, the caracteristics PARAHOT
  41. c at elevated temperature.
  42.  
  43. c In multiaxial situation
  44. c We are going to calculate the stress vector sigm6v and the stiffness
  45. c matrix H66t as a function of the strain incr. vector deps6v,
  46. c calculated from the strain vector at the beginning of the time step.
  47.  
  48. c i o
  49. c a : df/dsigma
  50. c df/dsigmax, df/dsigmay, df/dtauxy, df/dsigmaz
  51. c * Ai : Hardening parameter
  52.  
  53. c THG add concrete model 3D
  54.  
  55. IMPLICIT REAL*8 (A-B,D-H,O-Z)
  56. implicit integer (I-K,M,N)
  57. implicit logical (L)
  58. implicit character*10 (C)
  59.  
  60. dimension parahot3(idimpara3,npttot3)
  61. dimension epstr6v(6), esigmac(6)
  62. dimension esigma6v(6),H66t(6,6),deps6v(6),sigma6v(6)
  63. dimension rt1(1,6),rc1(1,6),H66tsave(6,6),H66(6,6)
  64.  
  65. i1 = 1
  66. i2 = 2
  67. i3 = 3
  68. i4 = 4
  69. i5 = 5
  70. i6 = 6
  71. i7 = 7
  72. i8 = 8
  73. i9 = 9
  74. i10=10
  75. i11=11
  76. i12=12
  77. i13=13
  78. i14=14
  79. i15=15
  80. i16=16
  81. i17=17
  82. i18=18
  83. i19=19
  84. i20=20
  85. i21=21
  86. i22=22
  87. i23=23
  88. i24=24
  89. i25=25
  90. i26=26
  91. i27=27
  92. i28=28
  93. i29=29
  94. i30=30
  95. i31=31
  96. i32=32
  97. i33=33
  98. i34=34
  99. i35=35
  100. i36=36
  101. i37=37
  102. i38=38
  103. i39=39
  104.  
  105. r2 = 2.
  106. r3 = 3.
  107.  
  108. c TEMPERATURE
  109. T = parahot3(idimpara3-39,ntot)
  110. TMAX = parahot3(idimpara3-38,ntot)
  111.  
  112.  
  113. c Compressive strength
  114. fc20 = parahot3(i3,ntot)
  115. rloc = rKfcsi(t,tmax)
  116. parahot3(i3,ntot) = rloc * fc20
  117. fc = parahot3(i3,ntot)
  118.  
  119. c Tensile strength
  120. ft20 = parahot3(i4,ntot)
  121. rloc = rkftco(t,tmax)
  122. parahot3(i4,ntot) = rloc * ft20
  123. if (parahot3(i4,ntot).lt.(1.d5)) then
  124. c minimal value for ft = 0.10 MPa implemented for numerical stability
  125. c so if ft entered by the user < 0.10 MPa, the concrete is considered as
  126. c previously cracked with dt>0.96 and ft = 1.00 MPa
  127. parahot3(i4,ntot) = (1.d6)
  128. parahot3(idimpara3-41,ntot) =
  129. . max(0.96d0,parahot3(idimpara3-41,ntot))
  130. endif
  131. ft = parahot3(i4,ntot)
  132.  
  133. c Strain at peak stress
  134. c Strain at peak stress according to EC2: epsc1,ec2
  135. rloc = 1.0d0
  136. epsc1ec2 = EPSc1co(tmax,rloc)
  137. c Strain at peak stress according to ENV (minimal value): epsc1,env
  138. epsc1env = EPSc1co(tmax,-rloc)
  139. c Strain at peak stress for the instantaneous stress-strain relationship
  140. c epsc1,ETC = (2 * epsc1,min + epsc1,EC2)/3
  141. epsc1etc = (2.*epsc1env+epsc1ec2)/3.
  142. parahot3(i5,ntot) = epsc1etc
  143.  
  144. epsc1ec220 = EPSc1co(20.d0,rloc)
  145. epsc1env20 = EPSc1co(20.d0,-rloc)
  146. epsc1etc20 = (2.*epsc1env20+epsc1ec220)/3.
  147.  
  148. c Modulus at the origin
  149. parahot3(i1,ntot) = 2.0d0*parahot3(i3,ntot)/parahot3(i5,ntot)
  150.  
  151. c Compressive limit of elasticity
  152. parahot3(i6,ntot) = (0.3d0)*parahot3(i3,ntot)
  153. fc0 = parahot3(i6,ntot)
  154.  
  155. c Biaxial compression strength: considered as 1.20 x compressive strength
  156. parahot3(7,ntot) = (1.2d0)*parahot3(3,ntot)
  157.  
  158. c thermal strain
  159. c first, save thermal strain at the precedent step
  160. parahot3(idimpara3-31,ntot) = parahot3(idimpara3,ntot)
  161. c then, calculate new thermal strain
  162. parahot3(idimpara3,ntot) = epsthsi(t,tmax)
  163.  
  164. c Modulus used in tension
  165. parahot3(i12,ntot) = parahot3(i1,ntot)
  166.  
  167. c Material parameter dc: damage in compression at peak stress
  168. dco = parahot3(i9,ntot)
  169.  
  170. c Material parameter x,c: adimensional densitive crack energy in compression
  171. c This parameter is assumed to be constant with temperature
  172. rxc = parahot3(i10,ntot)
  173.  
  174. c Model parameter kappa,c1: accumulated compressive plastic strain at peak stress
  175. denom = (2.-2.*dco)*((parahot3(i8,ntot))-1.)
  176. parahot3(i13,ntot) = -parahot3(i5,ntot)*(1.-2.*dco)/(denom)
  177. rkappa1 = parahot3(i13,ntot)
  178.  
  179. c Model Parameter ac
  180. parahot3(i14,ntot) = ( -LOG(1-dco)/(rkappa1) )
  181.  
  182. c Model Parameter bc
  183. parahot3(i15,ntot) = 2*fc*rxc/((1-rxc)*((fc0*rkappa1)+
  184. . ((fc-fc0)*rkappa1*LOG(r2))))
  185.  
  186. c Parameter at
  187. if (ft20.ge.1.d6) then
  188. parahot3(i16,ntot) = (7.0d0*ft20)/(12.0d0*
  189. . parahot3(i11,ntot))
  190. else
  191. c for numerical reason at cannot be equal to 0 even if ft0 = 0
  192. parahot3(i16,ntot) = (7.0d0*1.0d6)/(12.0d0*
  193. . parahot3(i11,ntot))
  194. endif
  195.  
  196. c TRANSIENT CREEP STRAIN
  197. do iloc=1,6
  198. epstr6v(iloc) = 0.d0
  199. enddo
  200. c Phi function
  201. c First, phi(ti) is saved as phiprev
  202. phiprev = parahot3(17,ntot)
  203. c Then, phi(ti+1) is calculated and saved as phi
  204. c phi = 2*(epsc1,EC2 - epsc1,min) / (3*(fc/fck))
  205. phi = 2.*(epsc1ec2-epsc1env)/(3.*(fc/fc20))
  206. parahot3(17,ntot) = phi
  207. c Calculation of transient creep strain
  208. c If the temperature is bellow Tmax (T is decreasing or it is
  209. c not a first heating) the transient creep strain is unchanged
  210. if (T.lt.(tmax-1.0d-3)) then
  211. continue
  212. c Else, it is a first-time heating
  213. c so there is an increment of transient creep strain under compressive stress
  214. else
  215. c unless the concrete is in the descending branch of the material law
  216. c if (kappa,c < kappa,c1) (i.e. material in the ascending branch of the
  217. c compressive part of the material law)
  218. rkappac = parahot3(idimpara3-36,ntot)
  219. if (rkappac.le.rkappa1) then
  220. c epstr,i+1 = epstr,i + [phi(Ti+1) - phi(Ti)] * (1 - dc,i) * (sigma,eff-,i / fck)
  221. dc = parahot3(idimpara3-40,ntot)
  222. do iloc=i1,i6
  223. esigmac(iloc) = parahot3(29+iloc,ntot)
  224. epstr6v(iloc) = parahot3(17+iloc,ntot) + ((phi-phiprev)/
  225. . fc20)*(1.0d0 - dc) * esigmac(iloc)
  226. end do
  227. c epstr cannot decrease in absolute value
  228. do iloc=i1,i6
  229. if (epstr6v(iloc).lt.parahot3(17+iloc,ntot)) then
  230. parahot3(17+iloc,ntot) = epstr6v(iloc)
  231. endif
  232. end do
  233. endif
  234. endif
  235.  
  236. c CALCULATE THE INCREMENT IN MECHANICAL STRAIN DEPS
  237. c (the increment in thermal strain is subtracted)
  238. c Deps = Deps - (Epsth(n)-Epsth(n-1))
  239. deps6v(1) = deps6v(1) - parahot3(idimpara3,ntot)
  240. . + parahot3(idimpara3-31,ntot)
  241. deps6v(2) = deps6v(2) - parahot3(idimpara3,ntot)
  242. . + parahot3(idimpara3-31,ntot)
  243. if (.not.lcp) deps6v(3) = deps6v(3) - parahot3(idimpara3,ntot)
  244. . + parahot3(idimpara3-31,ntot)
  245.  
  246. lg = .false.
  247. lbroyden = .false.
  248. lfulldamage = .false.
  249.  
  250. n = 1
  251. 23 continue
  252. do istep = 1,n
  253. c To improve convergence, the step can be subdivided
  254.  
  255.  
  256. c ***********************************************************
  257. c *********************
  258. c ! 3D concrete model !
  259. c ! PLASTIC BOX !
  260. c *********************
  261.  
  262. call DPRAN3D(parahot3(i1,ntot),idimpara3,H66t,rt1,rc1,
  263. c *******
  264. . esigma6v,deps6v,lg,ntot,lerror,i1,i2,i3,i4,i5,
  265. . i6,n,lnoconv,istep,lcp,lbroyden,H66,lfulldamage)
  266. c Plastic Box of the 3D plastic-damage concrete model by T. Gernay
  267. c Yield functions of Drucker Prager in compression and Rankine in tension
  268. c Calculate the effective stresses and the effective elastoplastic tangent
  269. c modulus Dt (implicit process)
  270. c ***********************************************************
  271.  
  272. if (lnoconv) then
  273. n = n*2
  274. c the step is subdivided by factor 2
  275. c we come back to the state at the previous (converged) step: rkappac = rkappaic
  276. parahot3(idimpara3-35,ntot) = parahot3(idimpara3-37,ntot)
  277. parahot3(idimpara3-34,ntot) = parahot3(idimpara3-36,ntot)
  278. if (n.lt.32) go to 23
  279. c the step can be divided in maximum 16 "sub-steps"
  280. if (.not.lbroyden) then
  281. c if still no convergence when dividing the size of the step => try another algorithm
  282. lbroyden = .true.
  283. n = 1
  284. go to 23
  285. endif
  286. if (.not.lfulldamage) then
  287. c still no convergence => if completely damaged, "remove" the point (sigma = 0)
  288. dc = parahot3(idimpara3-40,ntot)
  289. dt = parahot3(idimpara3-41,ntot)
  290. if ((dt.gt.0.97).or.(dc.gt.0.97)) then
  291. lfulldamage = .true.
  292. n = 1
  293. go to 23
  294. endif
  295. endif
  296. c if we arrive here, all the strategies to improve convergence have been tested with
  297. c no success => lerror is true (send back an error from the MATERIAL subroutines)
  298. lerror = .true.
  299. return
  300. endif
  301. enddo
  302.  
  303. c if ((NOEL.eq.18).and.(NPT.eq.23)) then
  304. * print *,'SORTIE DPRAN3D contraintes = ',esigma6v
  305. c endif
  306.  
  307. c *************************************************************
  308. c *********************
  309. c ! 3D concrete model !
  310. c ! DAMAGE BOX !
  311. c *********************
  312. call DAMTG3D(parahot3(i1,ntot),idimpara3,H66t,rt1,rc1,
  313. c *******
  314. . esigma6v,lg,ntot,lerror,i1,i2,i3,i4,i5,i6,H66,
  315. . lfulldamage,lcp)
  316. c Damage Box of the 3D plastic-damage concrete model by T. Gernay
  317. c Calculate the nominal stresses and the nominal consistent algorithmic
  318. c tangent modulus Ct (explicit process)
  319. c *************************************************************
  320.  
  321. c if ((NOEL.eq.18).and.(NPT.eq.23)) then
  322. * print *,'SORTIE DAMTG3D contraintes = ', sigma6v
  323. c endif
  324.  
  325.  
  326. c The tangent matrix is symmetrized
  327. c because SAFIR uses a symmetric matrix for the resolution of the equilibrium
  328. c T
  329. c B = ( A + A ) / 2
  330. do iloc=2,6
  331. rloc = H66t(1,iloc) + H66t(iloc,1)
  332. H66t(1,iloc) = rloc/(2.d0)
  333. H66t(iloc,1) = H66t(1,iloc)
  334. enddo
  335. do iloc=3,6
  336. rloc = H66t(2,iloc) + H66t(iloc,2)
  337. H66t(2,iloc) = rloc/(2.d0)
  338. H66t(iloc,2) = H66t(2,iloc)
  339. enddo
  340. do iloc=4,6
  341. rloc = H66t(3,iloc) + H66t(iloc,3)
  342. H66t(3,iloc) = rloc/(2.d0)
  343. H66t(iloc,3) = H66t(3,iloc)
  344. enddo
  345. do iloc=5,6
  346. rloc = H66t(4,iloc) + H66t(iloc,4)
  347. H66t(4,iloc) = rloc/(2.d0)
  348. H66t(iloc,4) = H66t(4,iloc)
  349. enddo
  350. rloc = H66t(5,6) + H66t(6,5)
  351. H66t(5,6) = rloc/(2.d0)
  352. H66t(6,5) = H66t(5,6)
  353.  
  354. c Save the nominal stresses
  355. do iloc=i1,i6
  356. sigma6v(iloc) = PARAHOT3(idimpara3-13+iloc,ntot)
  357. end do
  358.  
  359. c if (lcp) then
  360. cc Plane stress model
  361. c parahot2(idimpara2-i3,ntot) = sigma6v(1)
  362. c parahot2(idimpara2-i2,ntot) = sigma6v(2)
  363. c parahot2(idimpara2-i1,ntot) = sigma6v(4)
  364. c PARAHOT2(idimpara2-i7,ntot)=parahot3(idimpara3-18,ntot)
  365. c PARAHOT2(idimpara2-i6,ntot)=parahot3(idimpara3-17,ntot)
  366. c PARAHOT2(idimpara2-i5,ntot)=parahot3(idimpara3-15,ntot)
  367. c endif
  368. END
  369.  
  370.  
  371.  
  372.  
  373.  

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