Télécharger dpran3d.eso

Retour à la liste

Numérotation des lignes :

  1. C DPRAN3D SOURCE CB215821 16/04/21 21:16:31 8920
  2. SUBROUTINE dpran3d(parahot3,idimpara3,H66t,rt1,rc1,esigma,
  3. . deps,lg,ntot,lerror,i1,i2,i3,i4,i5,i6,n,lnoconv,
  4. . istep,lcp,lbroyden,H66,lfulldamage)
  5.  
  6. c parahot3 is the array that contains the material parameters and state variables
  7. c idimpara3 is the number of parameters for parahot3
  8. c H66t is the tangent matrix
  9. c rt1 and rc1 are the vectors used in DAMTG3D for the calculation of the tangent matrix
  10. c esigma is the effective stress vector
  11. c deps is the incremental strains vector
  12. c lg and lerror are logical
  13. c ntot is the integration point (global number for the whole structure): counter
  14. c n, istep, lnoconv, lbroyden and lfulldamage are used for the algorithmic strategy (step subdivider)
  15. c H66 is the elastic matrix
  16.  
  17. c ===========================================================================================
  18. c ! Plastic Box of the 3D plastic-damage concrete model by T. Gernay !
  19. c ! Yield functions of Drucker Prager in compression and Rankine in tension !
  20. c ! Calculate the effective stresses and the effective elastoplastic tangent !
  21. c ! modulus Dt (implicit process) !
  22. c ! !
  23. c ! input !
  24. c ! DEPS : Incremental strain vector (6 components) !
  25. c ! !
  26. c ! output !
  27. c ! H66t : tangent stiffness matrix (6x6) !
  28. c ! !
  29. c ! input and output !
  30. c ! PARAHOT3 : the first columns contain the material properties at elevated temperature !
  31. c ! the last columns contain strains, stresses, .... !
  32. c ! !
  33. c ! ESIGMA : Effective stress vector (6 components) !
  34. c ! !
  35. c ===========================================================================================
  36.  
  37. IMPLICIT REAL*8 (A-B,D-H,O-Z)
  38. implicit integer (I-K,M,N)
  39. implicit logical (L)
  40. implicit character*10 (C)
  41.  
  42. dimension deps(i6)
  43. c Incremental strain since the end of the previous time step (input)
  44. c i.e. since the last converged point
  45. dimension H66t(i6,i6)
  46. c Tangent matrix (output)
  47. dimension parahot3(idimpara3)
  48. c parahot3: Material prop. and various info. at elevated temp. (input and output)
  49. c parahot3(idimpara3-37,ntot): equivalent plastic strain in tension at the last converged step
  50. c parahot3(idimpara3-35,ntot): equivalent plastic strain in tension at the end of the current step
  51. c parahot3(idimpara3-36,ntot): equivalent plastic strain in compression at the last converged step
  52. c parahot3(idimpara3-34,ntot): equivalent plastic strain in compression at the end of the current step
  53. c parahot3(idimpara3-30:idimpara3-25,ntot) = strain at the end of the previous time step
  54. c parahot3(idimpara3-24:idimpara3-19,ntot) = plastic strains at the beginning of dpran3D
  55. c parahot3(idimpara3-18:idimpara3-13,ntot) = strain at the end of dpran3D
  56. c parahot3(idimpara3- 6:idimpara3- 1,ntot) = effective stress at the end of dpran3D
  57. c parahot3(idimpara3,ntot) = thermal strain
  58. dimension esigma(i6)
  59. c Effective stress (output)
  60. **** dimension EPSMEC(i6),EPSSIG(i6),esigmae(i6),esigmaf(i6)
  61. dimension EPSMEC(6),EPSSIG(6),esigmae(6),esigmaf(6)
  62. **** dimension H66(i6,i6),x(i2),d2gt(i6,i6),trav66b(i6,i6),trav6(i6)
  63. dimension H66(6,6),x(2),d2gt(6,6),trav66b(6,6),trav6(6)
  64. **** dimension trav66(i6,i6),Da(i6,i6),dft(i6)
  65. dimension trav66(6,6),Da(6,6),dft(6)
  66. **** dimension dfc(i6),dgc(i6),dgt(i6)
  67. dimension dfc(6),dgc(6),dgt(6)
  68. **** dimension trav1(1),DftD(i1,i6),DfcD(i1,i6),Ddgt(i6),Ddgc(i6)
  69. dimension trav1(1),DftD(1,6),DfcD(1,6),Ddgt(6),Ddgc(6)
  70. **** dimension trav16(i1,i6),rn1(i1,i6),rn2(i1,i6),rn3(i1,i6)
  71. dimension trav16(1,6),rn1(1,6),rn2(1,6),rn3(1,6)
  72. **** dimension rt1(i1,i6),rc1(i1,i6),vloc1(i1),vloc(i6)
  73. dimension rt1(1,6),rc1(1,6),vloc1(1),vloc(6)
  74. **** dimension Ddgtn1(i6,i6),Ddgtn2(i6,i6),Ddgcn3(i6,i6)
  75. dimension Ddgtn1(6,6),Ddgtn2(6,6),Ddgcn3(6,6)
  76. **** dimension Ddgcn4(i6,i6),rn4(i1,i6),dplas(6),d2gc(6,6),vloc2(6)
  77. dimension Ddgcn4(6,6),rn4(1,6),dplas(6),d2gc(6,6),vloc2(6)
  78.  
  79.  
  80. c Deviatoric projection operator P2 (stressdeviator = P2 * stress)
  81. c --------------------------------
  82. dimension P2(6,6)
  83. data P2 / 2.0d0, -1.0d0, -1.0d0, 0.0d0, 0.0d0, 0.0d0,
  84. . -1.0d0, 2.0d0, -1.0d0, 0.0d0, 0.0d0, 0.0d0,
  85. . -1.0d0, -1.0d0, 2.0d0, 0.0d0, 0.0d0, 0.0d0,
  86. . 0.0d0, 0.0d0, 0.0d0, 6.0d0, 0.0d0, 0.0d0,
  87. . 0.0d0, 0.0d0, 0.0d0, 0.0d0, 6.0d0, 0.0d0,
  88. . 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 6.0d0 /
  89.  
  90.  
  91. c Matrix of unity I, or U
  92. c -----------------------
  93. dimension U(6,6)
  94. data U / 1.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0,
  95. . 0.0d0, 1.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0,
  96. . 0.0d0, 0.0d0, 1.0d0, 0.0d0, 0.0d0, 0.0d0,
  97. . 0.0d0, 0.0d0, 0.0d0, 1.0d0, 0.0d0, 0.0d0,
  98. . 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1.0d0, 0.0d0,
  99. . 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1.0d0 /
  100.  
  101. c Vector pi2
  102. c -----------------------
  103. dimension pi2(6)
  104. data pi2 / 1.0d0,
  105. . 1.0d0,
  106. . 1.0d0,
  107. . 0.0d0,
  108. . 0.0d0,
  109. . 0.0d0 /
  110.  
  111.  
  112.  
  113. *****
  114. * MESSAGES D'ERREUR ( SUPPRESSION DES AUTOMATIC OBJECTS)
  115. IF(I1.GT.1) PRINT *, ' DPRAND3D - ERREUR I1 = ', I1, ' > 1 '
  116. IF(I6.GT.6) PRINT *, ' DPRAND3D - ERREUR I6 = ', I6, ' > 6 '
  117. *****
  118.  
  119. lerror = .false.
  120. lnoconv = .false.
  121.  
  122. P2(1,1) = 2.0d0
  123. P2(2,2) = 2.0d0
  124. P2(3,3) = 2.0d0
  125. i7 = 7
  126. i8 = 8
  127. i9 = 9
  128. i10=10
  129. i11=11
  130. i12=12
  131. i13=13
  132. i14=14
  133. i15=15
  134. i16=16
  135. i17=17
  136. i18=18
  137. i19=19
  138. i20=20
  139. i21=21
  140. i22=22
  141. i23=23
  142. i24=24
  143. i25=25
  144. i26=26
  145. i27=27
  146. i28=28
  147. i29=29
  148. i30=30
  149. i31=31
  150. i32=32
  151. i33=33
  152. i34=34
  153. i35=35
  154. i36=36
  155. i37=37
  156. i38=38
  157. i39=39
  158.  
  159. r0 = 0.
  160. r1 = 1.
  161. r2 = 2.
  162. r10 = 10.
  163.  
  164.  
  165. Eo = parahot3(i1)
  166. poison = parahot3(i2)
  167. fc = parahot3(i3)
  168. ft = parahot3(i4)
  169. c Strain at peak stress
  170. epsc1 = parahot3(i5)
  171. c Compressive limit of elasticity
  172. fc0 = parahot3(i6)
  173. c Biaxial compressive strength
  174. fb = parahot3(i7)
  175. c Dilatancy parameter
  176. ag = parahot3(i8)
  177. c Damage to peak stress
  178. dc = parahot3(i9)
  179. c Compressive dissipated energy parameter
  180. rxc = parahot3(i10)
  181. c Crack energy in tension
  182. gt = parahot3(i11)
  183. c Model parameters
  184. rkappa1 = parahot3(i13)
  185. ac = parahot3(i14)
  186. bc = parahot3(i15)
  187.  
  188. c small strain for starting the return mapping algorithm
  189. tol = 0.0001*epsc1
  190.  
  191. c Hardening parameter for tension at the beginning of the step
  192. rkappait= parahot3(idimpara3-37)
  193. c Hardening parameter for compression at the beginning of the step
  194. rkappaic= parahot3(idimpara3-36)
  195.  
  196. c ===================================================================
  197. c ! !
  198. c ! CALCULATION OF THE MATRIX OF ELASTICITY !
  199. c ! ======================================= !
  200. c ===================================================================
  201.  
  202. c Matrix of elasticity, (6x6, full 3D)
  203. c --------------------
  204. c [Eo*(1-nu)/((1+nu)*(1-2*nu)) Eo*nu/((1+nu)*(1-2*nu)) Eo*nu/((1+nu)*(1-2*nu)) 0 0 0 ]
  205. c [Eo*nu/((1+nu)*(1-2*nu)) Eo*(1-nu)/((1+nu)*(1-2*nu)) Eo*nu/((1+nu)*(1-2*nu)) 0 0 0 ]
  206. c H66 = [Eo*nu/((1+nu)*(1-2*nu)) Eo*nu/((1+nu)*(1-2*nu)) Eo*(1-nu)/((1+nu)*(1-2*nu)) 0 0 0 ]
  207. c [0 0 0 Eo/(1+nu) 0 0 ]
  208. c [0 0 0 0 Eo/(1+nu) 0 ]
  209. c [0 0 0 0 0 Eo/(1+nu) ]
  210.  
  211. c Matrix of elasticity, (3x3, plane stress)
  212. c --------------------
  213. c [ 1 nu 0 0 0 0]
  214. c H66 = Eo/(1-nu*nu) [nu 1 0 0 0 0]
  215. c [ 0 0 0 0 0 0]
  216. c [ 0 0 0 (1-nu)/2 0 0]
  217. c [ 0 0 0 0 0 0]
  218. c [ 0 0 0 0 0 0]
  219. c++
  220.  
  221. do jloc=i1,i6
  222. do iloc=i1,i6
  223. H66(iloc,jloc) = r0
  224. end do
  225. end do
  226.  
  227. if (.not.lcp) then
  228. c 3D model
  229. rloc = Eo/((r1+poison)*(r1-2.0d0*poison))
  230. H66(i1,i1) = rloc*(r1-poison)
  231. H66(i2,i2) = rloc*(r1-poison)
  232. H66(i3,i3) = rloc*(r1-poison)
  233. H66(i1,i2) = rloc*poison
  234. H66(i1,i3) = rloc*poison
  235. H66(i2,i3) = rloc*poison
  236. H66(i2,i1) = rloc*poison
  237. H66(i3,i1) = rloc*poison
  238. H66(i3,i2) = rloc*poison
  239. rloc2 = Eo/(2.0d0*(r1+poison))
  240. H66(i4,i4) =rloc2
  241. H66(i5,i5) =rloc2
  242. H66(i6,i6) =rloc2
  243. else
  244. c PLANE STRESS
  245. rloc = Eo/(r1-poison*poison)
  246. H66(i1,i1) = rloc
  247. H66(i2,i1) = poison * (rloc)
  248. H66(i1,i2) = poison * (rloc)
  249. H66(i2,i2) = rloc
  250. H66(4,4) = rloc * ((1.0d0-poison)/(2.0d0))
  251. deps(3) = 0.d0
  252. deps(5) = 0.d0
  253. deps(6) = 0.d0
  254. endif
  255.  
  256. c =============================================================
  257. c ! !
  258. c ! ELASTIC STRESS PREDICTOR !
  259. c ! ======================== !
  260. c =============================================================
  261.  
  262. c Total mechanical strain: EPSMEC[i,j] = epsmec[i-1,conv] + deps[i,j] - epspl[i-1,j]
  263. c epsplastic has to be subtracted because Sigma initial = 0.
  264. do iloc=i1,i6
  265. EPSMEC(iloc) = parahot3(idimpara3-31+iloc)
  266. . + ( (DEPS(iloc) * istep)/n )
  267. . - parahot3(idimpara3-25+iloc)
  268. c Calculation of the INSTANTANEOUS STRESS-DEPENDENT STRAIN
  269. EPSSIG(iloc) = EPSMEC(iloc) - parahot3(17+iloc)
  270. end do
  271. c The effective stresses are calculated as functions of the instantaneous
  272. c stress-dependent strain eps,sig because the transient creep strain is explicitly
  273. c calculated apart from these constitutive relationships. Still in parahot, the mechanical
  274. c strains are saved, not the instantaneous stress-dependent strains.
  275.  
  276.  
  277. if ( (.not.lnoconv).and.(n.gt.1) ) then
  278. c We are in a subdivision of the step (for algorithmic reasons)
  279. c The plastic strain considered in epsmec has to be adjusted
  280. rkappait = parahot3(idimpara3-35)
  281. rkappaic = parahot3(idimpara3-34)
  282. x1 = rkappait - parahot3(idimpara3-37)
  283. x2 = rkappaic - parahot3(idimpara3-36)
  284. call dplas3D(x1,x2,P2,pi2,ag,fc,fb,parahot3,idimpara3,
  285. . dplas,lcp)
  286. do iloc=i1,i6
  287. EPSSIG(iloc) = EPSMEC(iloc)-dplas(iloc)-parahot3(17+iloc)
  288. end do
  289. endif
  290.  
  291. c Elastic stress predictor: eff sigma[i,j] = Hel x EPSSIG[i,j]
  292. call mulab(H66,EPSSIG,ESIGMAE,i6,i6,i6,i1)
  293.  
  294. *
  295. * print *,' prediction contraintes = ',esigmae
  296. *
  297.  
  298. c ==============================================================
  299. c ! !
  300. c ! CONCRETE COMPLETELY DAMAGED !
  301. c ! =========================== !
  302. c ==============================================================
  303. dc = parahot3(idimpara3-40)
  304. dt = parahot3(idimpara3-41)
  305.  
  306. if ((dc.gt.0.98).or.(lfulldamage).or.(dt.gt.0.98)) then
  307. c Concrete completely damaged
  308. do iloc=1,6
  309. ESIGMA(iloc) = 0.d0
  310. end do
  311. do iloc=1,6
  312. parahot3(idimpara3-13+iloc) = esigma(iloc)
  313. end do
  314. do jloc=1,6
  315. do iloc=1,6
  316. H66t(iloc,jloc) = 0.d0
  317. end do
  318. end do
  319. c store the mechanical strain at the end of this iteration
  320. c mechanical strain = epsmec + epspl
  321. do iloc=i1,i6
  322. parahot3(idimpara3-19+iloc) = EPSMEC(iloc) +
  323. . parahot3(idimpara3-25+iloc)
  324. end do
  325. return
  326. endif
  327.  
  328. c ==========================================================================
  329. c ! !
  330. c ! CALCULATE THE PLASTICITY FUNCTIONS WITH THE ELASTIC STRESS PREDICTOR !
  331. c ! ==================================================================== !
  332. c ==========================================================================
  333.  
  334. c Calculation of the hardening functions
  335. tau_compr = fhardcc(rkappaic,fc,fc0,rkappa1,
  336. . ac,bc,lerror)
  337. tau_tension = fhardct(rkappait,ft)
  338. if (lerror) then
  339. write(2,*) '>Error returned from compressive hardening
  340. . to DPRAN3d'
  341. return
  342. endif
  343.  
  344. c Calculation of the yield functions for the elastic stress predictor
  345. alpha = (fb-fc)/(2.*fb-fc)
  346. f_compr = DPfunc(esigmae,alpha)
  347.  
  348. f_tension = Ranfunc(esigmae,i3,i6,parahot3,
  349. . idimpara3,lcp)
  350.  
  351. c++
  352. f_compr = f_compr - ( (1.0d0-alpha) * tau_compr )
  353. f_tension = f_tension - tau_tension
  354.  
  355. c =======================================================
  356. c ! !
  357. c ! TEST ELASTIC OR PLASTIC STEP !
  358. c ! ============================ !
  359. c =======================================================
  360.  
  361. c The plasticity functions must be satisfied with a precision of 10 Pa, i.e. 0.000010 N/mm²
  362. if ((f_compr.lt.-r10).and.(f_tension.lt.-r10)) then
  363.  
  364. C ***************
  365. c ! ELASTICITY !
  366. C ***************
  367.  
  368. do iloc=i1,i6
  369. ESIGMA(iloc) = ESIGMAE(iloc)
  370. end do
  371.  
  372. do jloc=i1,i6
  373. do iloc=i1,i6
  374. H66T(iloc,jloc) = H66(iloc,jloc)
  375. end do
  376. end do
  377. c store the stress at the end of this iteration
  378. do iloc=i1,i6
  379. parahot3(idimpara3-7+iloc) = ESIGMAE(iloc)
  380. end do
  381.  
  382. c The equivalent plastic strain in tension has not changed
  383. parahot3(idimpara3-35) = parahot3(idimpara3-37)
  384. c The equivalent plastic strain in compression has not changed
  385. parahot3(idimpara3-34) = parahot3(idimpara3-36)
  386.  
  387. c vectors for the calculation of H66t in subroutine damtg3d
  388. do iloc=i1,i6
  389. rt1(i1,iloc) = r0
  390. rc1(i1,iloc) = r0
  391. end do
  392.  
  393. else
  394.  
  395. C ***************
  396. c ! PLASTICITY !
  397. C ***************
  398.  
  399. * PRINT *, ' >>>>>>>> ON ENTRE EN PLASTICITE'
  400.  
  401. c RETURN MAPPING
  402. c --------------
  403. c Calculation of Dlambdac and Dlambdat, the increments of the eq. plastic strain which brings
  404. c back on the plasticity surfaces
  405. c We have to solve the following non linear system of equations
  406. c rkt * function_t(Dlambdac,Dlambdat) + (1-rkt) * Dlambdat = 0
  407. c rkc * function_c(Dlambdac,Dlambdat) + (1-rkc) * Dlambdac = 0
  408. c with rkc and rkt = 1 or 0, depending whether the related surface of plasticity is violated
  409.  
  410. if (f_compr.lt.-r10) then
  411. c only the tension surface is violated by the trial stress
  412. Dlambdac = r0
  413. rkc = r0
  414. rkt = r1
  415.  
  416. else if (f_tension.lt.-r10) then
  417. c only the compression surface is violated by the trial stress
  418. Dlambdat = r0
  419. rkc = r1
  420. rkt = r0
  421.  
  422. else
  423. c Both surfaces are violated by the trial stress
  424. rkc = r1
  425. rkt = r1
  426.  
  427. endif
  428.  
  429. c X=0 at the beginning of every SAFIR iteration, X being the vector (Dlambdat,Dlambdac)
  430. X(i1) = r0
  431. X(i2) = r0
  432.  
  433. c ===================================================================================
  434. c ! !
  435. c ! CALL TO SUBROUTINE ALGOPLAS TO SOLVE THE NONLINEAR SYSTEM IN DLAMBDAT, DLAMBDAC !
  436. c ! =============================================================================== !
  437. c ===================================================================================
  438.  
  439. if (.not.lcp) then
  440. c 3D MODEL
  441. if (.not.lbroyden) then
  442. call Algoplas1(rkt,rkc,esigmae,H66,P2,pi2,alpha,
  443. . rkappait,rkappaic,x,fc,fb,ft,Eo,poison,lg,ntot,
  444. . fc0,rkappa1,ac,bc,ag,lerror,i1,i2,i3,i4,i5,i6,
  445. . parahot3,idimpara3,deps,lnoconv,lcp,U)
  446. else
  447. call Algoplas2(rkt,rkc,esigmae,H66,P2,pi2,alpha,
  448. . rkappait,rkappaic,x,fc,fb,ft,Eo,poison,lg,ntot,
  449. . fc0,rkappa1,ac,bc,ag,lerror,i1,i2,i3,i4,i5,i6,
  450. . parahot3,idimpara3,deps,lnoconv,lcp,U)
  451. endif
  452. if (lnoconv) then
  453. return
  454. endif
  455.  
  456. if (lerror) then
  457. write(2,*) 'lerror is returned from ALGOPLAS to
  458. . DPRAN3D. Ntot =',ntot
  459. return
  460. endif
  461.  
  462. c The increment of plastic multipliers Dlambdat, Dlambdac have been found
  463.  
  464. c ==========================================================
  465. c ! !
  466. c ! UPDATE OF THE HARDENING PARAMETERS AND THE STRESSES !
  467. c ! =================================================== !
  468. c ==========================================================
  469.  
  470. Dlambdat = x(i1)
  471. Dlambdac = x(i2)
  472.  
  473. if (n.gt.1) then
  474. c we are currently in a step subdivision (for algorithmic reason)
  475. c we take the value saved at the end of the previous step-subdivision
  476. rkappait = parahot3(idimpara3-35)
  477. rkappaic = parahot3(idimpara3-34)
  478. else
  479. c we take the value at the previous (converged) time step
  480. rkappait= parahot3(idimpara3-37)
  481. rkappaic= parahot3(idimpara3-36)
  482. endif
  483.  
  484. c update of the hardening parameters
  485. rkappat = rkappait+Dlambdat
  486. rkappac = rkappaic+Dlambdac
  487.  
  488. c Update of the stresses: esigmaf = F(Dlambdat,Dlambdac (and esigmae))
  489. ag = parahot3(8)
  490. call sigf3D_implicit(esigmae,Dlambdat,Dlambdac,esigmaf,H66,
  491. c ***************
  492. . P2,pi2,ag,i6,fc,fb,parahot3,idimpara3,lerror,lcp,
  493. . esigmae,U)
  494.  
  495. if (lerror) then
  496. return
  497. endif
  498.  
  499. else
  500.  
  501. c PLANE STRESS MODEL
  502. if (.not.lbroyden) then
  503. call Algoplas1cp(rkt,rkc,esigmae,H66,P2,pi2,alpha,
  504. . rkappait,rkappaic,x,fc,fb,ft,Eo,poison,lg,ntot,
  505. . fc0,rkappa1,ac,bc,ag,lerror,i1,i2,i3,i4,i5,i6,
  506. . parahot3,idimpara3,deps,esigmaf,rkappat,rkappac,
  507. . lnoconv,lcp,U)
  508. else
  509. call Algoplas2cp(rkt,rkc,esigmae,H66,P2,pi2,alpha,
  510. . rkappait,rkappaic,x,fc,fb,ft,Eo,poison,lg,ntot,
  511. . fc0,rkappa1,ac,bc,ag,lerror,i1,i2,i3,i4,i5,i6,
  512. . parahot3,idimpara3,deps,esigmaf,rkappat,rkappac,
  513. . lnoconv,lcp,U)
  514. endif
  515.  
  516. if (lnoconv) then
  517. return
  518. endif
  519.  
  520. if (lerror) then
  521. write(2,*) 'lerror is returned from ALGOPLAS to
  522. . DPRAN3D. Ntot =',ntot
  523. return
  524. endif
  525.  
  526. Dlambdat = x(i1)
  527. Dlambdac = x(i2)
  528. if (n.gt.1) then
  529. c we are currently in a step subdivision (for algorithmic reason)
  530. c we take the value saved at the end of the previous step-subdivision
  531. rkappait = parahot3(idimpara3-35)
  532. rkappaic = parahot3(idimpara3-34)
  533. else
  534. c we take the value at the previous (converged) time step
  535. rkappait= parahot3(idimpara3-37)
  536. rkappaic= parahot3(idimpara3-36)
  537. endif
  538.  
  539. c update of the hardening parameters
  540. rkappat = rkappait+Dlambdat
  541. rkappac = rkappaic+Dlambdac
  542.  
  543. endif
  544.  
  545. c store the stress at the end of this iteration
  546. do iloc=i1,i6
  547. ESIGMA(iloc) = ESIGMAF(iloc)
  548. end do
  549. do iloc=i1,i6
  550. parahot3(idimpara3-i7+iloc) = ESIGMAF(iloc)
  551. end do
  552.  
  553. c store the accumulated plastic strains at the end of this iteration
  554. parahot3(idimpara3-35) = rkappat
  555. parahot3(idimpara3-34) = rkappac
  556.  
  557. c THE EFFECTIVE STRESSES HAVE BEEN FOUND
  558. c **************************************
  559.  
  560. c ================================================================
  561. c ! TANGENT MATRIX. (effective!) !
  562. c ! ---------------------------- !
  563. c ! Here, we calculate the effective elastoplastic tangent !
  564. c ! modulus Dt. As this is a damage-plastic model, we have !
  565. c ! then to correct this modulus Dt with the damage terms to !
  566. c ! obtain the (nominal) consistent algorithmic tangent modulus !
  567. c ! Ct. This will be done in DAMTG3D. Here, we calculate Dt. !
  568. c ================================================================
  569.  
  570. c APEX of Drucker-Prager surface
  571. alpha = (fb-fc)/(2.*fb-fc)
  572. f_compr = DPfunc(esigmae,alpha)
  573. if ((ABS(f_compr).le.100.).and.(rkc.eq.1)) then
  574. do jloc=i1,i6
  575. do iloc=i1,i6
  576. H66t(iloc,jloc) = 0.d0
  577. end do
  578. enddo
  579. do iloc=i1,i6
  580. parahot3(idimpara3-19+iloc) = EPSMEC(iloc) +
  581. . parahot3(idimpara3-25+iloc)
  582. end do
  583. RETURN
  584. endif
  585.  
  586. if (lcp) then
  587. c recalculate the full 3D matrix
  588. do jloc=i1,i6
  589. do iloc=i1,i6
  590. H66(iloc,jloc) = r0
  591. end do
  592. end do
  593. rloc = Eo/((r1+poison)*(r1-2.0d0*poison))
  594. H66(i1,i1) = rloc*(r1-poison)
  595. H66(i2,i2) = rloc*(r1-poison)
  596. H66(i3,i3) = rloc*(r1-poison)
  597. H66(i1,i2) = rloc*poison
  598. H66(i1,i3) = rloc*poison
  599. H66(i2,i3) = rloc*poison
  600. H66(i2,i1) = rloc*poison
  601. H66(i3,i1) = rloc*poison
  602. H66(i3,i2) = rloc*poison
  603. rloc2 = Eo/(2.0d0*(r1+poison))
  604. H66(i4,i4) =rloc2
  605. H66(i5,i5) =rloc2
  606. H66(i6,i6) =rloc2
  607. endif
  608. if (lcp) then
  609. tau_tension = fhardct(rkappat,ft)
  610. psi1 = tau_tension - (esigmaf(1)+esigmaf(2))/2
  611. if (abs(psi1).le.fc*1.e-5) then
  612. do jloc=1,6
  613. do iloc=1,6
  614. H66t(iloc,jloc) = 0
  615. enddo
  616. enddo
  617. do iloc=i1,i6
  618. parahot3(idimpara3-19+iloc) = EPSMEC(iloc) +
  619. . parahot3(idimpara3-25+iloc)
  620. end do
  621. RETURN
  622. endif
  623. endif
  624.  
  625. c First, calculation of Da, the effective algorithmic modulus
  626. c **************************************************
  627.  
  628. call ddsigrank(esigmaf,d2gt,i1,i2,i3,i4,i5,i6,lcp,
  629. . parahot3,idimpara3)
  630.  
  631. do jloc=i1,i6
  632. do iloc=i1,i6
  633. TRAV66b(iloc,jloc) = Dlambdat*d2gt(iloc,jloc)
  634. c [6;6]
  635. end do
  636. end do
  637. call mulAB(H66,TRAV66b,TRAV66b,i6,i6,i6,i6)
  638. c [6;6] = [6;6] x [6;6]
  639.  
  640. tau_c = fhardcc(rkappac,fc,fc0,rkappa1,
  641. . ac,bc,lerror)
  642. if (lerror) then
  643. write(2,*) '>Error returned from compressive hardening
  644. . to DPRAN3d'
  645. return
  646. endif
  647. call mulATB(pi2,esigmaf,vloc1,i6,i1,i6,i1)
  648. c [1] = [1;6] x [6;1]
  649. psi2 = (1d0-alpha)*tau_c - alpha*vloc1(i1)
  650.  
  651. c First come the second order terms.
  652. call mulab(P2,ESIGMAF,TRAV6,i6,i6,i6,i1)
  653. c [6;1] = [6;6] x [6;1]
  654.  
  655. c t t t
  656. c T66 = T6 T6 = P2 s s P2
  657. call mulABT(trav6,trav6,trav66,i6,i1,i6,i1)
  658. c [6;6] = [6;1] x [1;6]
  659. test = 1.E+6
  660. psi2lim=max(1.E-7*test,1.d-3)
  661. if (abs(psi2).le.psi2lim) then
  662. do jloc=i1,i6
  663. do iloc=i1,i6
  664. TRAV66(iloc,jloc) = r0
  665. end do
  666. end do
  667. else
  668. do jloc=i1,i6
  669. do iloc=i1,i6
  670. TRAV66(iloc,jloc) = (P2(iloc,jloc)-TRAV66(iloc,jloc)
  671. . /(r2*psi2*psi2))/(r2*psi2)
  672. c [6;6]
  673. end do
  674. end do
  675. endif
  676.  
  677. do jloc=i1,i6
  678. do iloc=i1,i6
  679. TRAV66(iloc,jloc) = Dlambdac*TRAV66(iloc,jloc)
  680. end do
  681. end do
  682. call mulAB(H66,TRAV66,TRAV66,i6,i6,i6,i6)
  683. c **** ***** [6;6]=[6;6]x[6;6]
  684. do jloc=i1,i6
  685. do iloc=i1,i6
  686. TRAV66(iloc,jloc) = U(iloc,jloc)+TRAV66(iloc,jloc)
  687. c TRAV66(iloc,jloc)=U(iloc,jloc)+TRAV66(iloc,jloc)
  688. . +TRAV66b(iloc,jloc)
  689. c TRAV66b requires the calculation of d2gt - we prefer not to take it
  690. end do
  691. end do
  692. c t t -1
  693. c Da = ( I + Dlambdat H66 d2gt + Dlambdac H66 (P2/2psi2)-(P2 s s P2)/(4psi2³) ) H66
  694.  
  695. call invert9(TRAV66,i6,i6)
  696. c **** ******* [6;6]
  697. call mulab(TRAV66,H66,Da,i6,i6,i6,i6)
  698. c [6;6]=[6;6]x[6;6]
  699.  
  700. c The effective algorithmic modulus Da has been obtained
  701. c Now, calculation of the effective elastoplastic tangent modulus Dt
  702. c ******************************************************************
  703. c Now come the first order terms
  704.  
  705. dtau1 = dfhardct(rkappat)
  706. dtau2 = dfhardcc(rkappac,fc,fc0,rkappa1,ac,bc)
  707.  
  708. c Calculation of n1, n2, n3, n4
  709.  
  710. call dsigRank(esigmaf,dgt,i3,i6,parahot3,idimpara3,
  711. c **** ******** [6;1]
  712. . rkappait,lcp)
  713. rloc = 0.d0
  714. do iloc=1,6
  715. rloc = rloc + (dgt(iloc))**2
  716. enddo
  717. if (rloc.eq.0.d0) rkt = 0.d0
  718. do iloc=i1,i6
  719. dft(iloc) = rkt*dgt(iloc)
  720. c [6;1] = [1] x [6;1]
  721. end do
  722.  
  723. e1 = r1-rkt*(r1+dtau1)
  724. c [1]
  725. e2 = r1-rkc*(r1+dtau2)
  726. c [1]
  727.  
  728. fc = parahot3(i3)
  729. c Compressive strength
  730. fb = parahot3(i7)
  731. c Biaxial compressive strength
  732.  
  733. call dsigDrucker(fc,fb,esigmaf,dfc,alpha,P2,
  734. c **** **************** [6;1]
  735. . pi2,i1,i6,lcp)
  736. do iloc=i1,i6
  737. dfc(iloc) = rkc*dfc(iloc)
  738. c [6;1] = [1]x[6;1]
  739. end do
  740. ag = parahot3(i8)
  741. call dsiggc(esigmaf,dgc,ag,P2,pi2,i1,i6,fc,fb,lcp)
  742. c **** **************** [6;1]
  743.  
  744. call mulATB(dft,Da,trav16,i6,i1,i6,i6)
  745. c [1;6] = [1;6]x[6;6]
  746. do iloc=i1,i6
  747. dftD(i1,iloc) = trav16(i1,iloc)
  748. end do
  749.  
  750. call mulab(Da,dgt,Ddgt,i6,i6,i6,i1)
  751. c [6;1] = [6;6] x [6;1]
  752.  
  753. call mulATB(dfc,Da,trav16,i6,i1,i6,i6)
  754. c [1;6] = [1;6]x[6;6]
  755. do iloc=i1,i6
  756. dfcD(i1,iloc) = trav16(i1,iloc)
  757. end do
  758.  
  759. call mulab(Da,dgc,Ddgc,i6,i6,i6,i1)
  760. c [6;1] = [6;6]x[6;1]
  761.  
  762. c n1
  763. call mulATB(dft,Ddgc,trav1,i6,i1,i6,i1)
  764. c **** ****** [1]=[1;6]x[6;1]
  765. dftDdgc = trav1(i1)
  766. c [1]
  767. do iloc=i1,i6
  768. rn1(i1,iloc) = -dftDdgc*dfcD(i1,iloc)
  769. end do
  770. c [1;6]=[1]x[1;6]
  771.  
  772. c n2 [1]=[1;6]x[6;1]
  773. call mulATB(dfc,Ddgc,trav1,i6,i1,i6,i1)
  774. c **** ******
  775. dfcDdgc = trav1(i1)
  776. do iloc=i1,i6
  777. rn2(i1,iloc) = (dfcDdgc-e2)*dftD(i1,iloc)
  778. end do
  779.  
  780. c n3 [1]=[1;6]x[6;1]
  781. call mulATB(dft,Ddgt,trav1,i6,i1,i6,i1)
  782. c **** ******
  783. dftDdgt = trav1(i1)
  784. do iloc=i1,i6
  785. rn3(i1,iloc) = (dftDdgt-e1)*dfcD(i1,iloc)
  786. end do
  787.  
  788. c n4 [1]=[1;6]x[6;1]
  789. call mulATB(dfc,Ddgt,trav1,i6,i1,i6,i1)
  790. c **** ******
  791. dfcDdgt = trav1(i1)
  792. do iloc=i1,i6
  793. rn4(i1,iloc) = -(dfcDdgt)*dftD(i1,iloc)
  794. end do
  795.  
  796. c denominator [1]
  797. Denom = (dftDdgt-e1)*(dfcDdgc-e2) - (dfcDdgt*dftDdgc)
  798.  
  799. c Calculation of rt1 and rc1 that are used in Mater 2 for calculation of Ct
  800. do iloc=i1,i6
  801. rt1(i1,iloc) = (rn1(i1,iloc)+rn2(i1,iloc))/denom
  802. rc1(i1,iloc) = (rn3(i1,iloc)+rn4(i1,iloc))/denom
  803. end do
  804.  
  805. do iloc=i1,i6
  806. vloc(iloc) = rn1(i1,iloc)
  807. end do
  808. call mulABT(Ddgt,vloc,Ddgtn1,i6,i1,i6,i1)
  809. c **** ****** [6;6]=[6;1]x[1;6]
  810. do iloc=i1,i6
  811. vloc(iloc) = rn2(i1,iloc)
  812. end do
  813. call mulABT(Ddgt,vloc,Ddgtn2,i6,i1,i6,i1)
  814. c **** ****** [6;6]=[6;1]x[1;6]
  815. do iloc=i1,i6
  816. vloc(iloc) = rn3(i1,iloc)
  817. end do
  818. call mulABT(Ddgc,vloc,Ddgcn3,i6,i1,i6,i1)
  819. c **** ****** [6;6]=[6;1]x[1;6]
  820. do iloc=i1,i6
  821. vloc(iloc) = rn4(i1,iloc)
  822. end do
  823. call mulABT(Ddgc,vloc,Ddgcn4,i6,i1,i6,i1)
  824. c **** ****** [6;6]=[6;1]x[1;6]
  825.  
  826. do jloc=i1,i6
  827. do iloc=i1,i6
  828. TRAV66(iloc,jloc) = Ddgtn1(iloc,jloc) + Ddgtn2(iloc,jloc)
  829. . + Ddgcn3(iloc,jloc) + Ddgcn4(iloc,jloc)
  830. c [6;6]
  831. TRAV66(iloc,jloc) = TRAV66(iloc,jloc)/Denom
  832. c [6;6]
  833. end do
  834. end do
  835.  
  836. do jloc=i1,i6
  837. do iloc=i1,i6
  838. H66t(iloc,jloc) = Da(iloc,jloc) - TRAV66(iloc,jloc)
  839. end do
  840. end do
  841. c END OF CALCULATION OF EFFECTIVE TANGENT MATRIX
  842. c **************************************************
  843.  
  844. endif
  845. c END OF PLASTICITY
  846.  
  847. c store the mechanical strain at the end of this iteration
  848. c mechanical strain = epsmec + epspl
  849. do iloc=i1,i6
  850. parahot3(idimpara3-19+iloc) = EPSMEC(iloc) +
  851. . parahot3(idimpara3-25+iloc)
  852. end do
  853.  
  854. RETURN
  855. END
  856.  
  857.  
  858.  
  859.  
  860.  

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