Télécharger algop2.eso

Retour à la liste

Numérotation des lignes :

  1. C ALGOP2 SOURCE CB215821 16/04/21 21:15:07 8920
  2. Subroutine algoplas2(rkt,rkc,esigmae,H66,P2,pi2,alpha,
  3. . rkappait,rkappaic,x,fc,fb,ft,Eo,poison,lg,ntot,
  4. . fc0,rkappa1,ac,bc,ag,lerror,i1,i2,i3,i4,i5,i6,
  5. . parahot3,idimpara3,deps,lnoconv,lcp,U)
  6.  
  7. c SUBROUTINE added by THG
  8.  
  9. c THIS SUBROUTINE IS CALLED BY DPRAN3D
  10. c the system is solved by the method Newton
  11.  
  12. c rkt x function_t(x(1),x(2)) + (1-rkt) x(1) = 0
  13. c rkc x function_c(x(1),x(2)) + (1-rkc) x(2) = 0
  14.  
  15. c with
  16. c rkt and rkc = 0 or 1
  17. c x(1)=Dlambdat
  18. c x(2)=Dlambdac
  19.  
  20. IMPLICIT REAL*8 (A-B,D-H,O-Z)
  21. implicit integer (I-K,M,N)
  22. implicit logical (L)
  23. implicit character*10 (C)
  24.  
  25. **** dimension B(i2,i2)
  26. dimension B(2,2)
  27. c jacobian matrix for the system of equations.
  28. **** dimension Binv(i2,i2)
  29. dimension Binv(2,2)
  30. c invert of the jacobian matrix for the system of equations.
  31. **** dimension df(i2)
  32. dimension df(2)
  33. **** dimension dfc(i6)
  34. dimension dfc(6)
  35. c derivative of the compression plasticity function
  36. **** dimension dft(i6)
  37. dimension dft(6)
  38. c derivative of the tension plasticity function
  39. **** dimension dx(i2)
  40. dimension dx(2)
  41. c increment of the solution
  42. **** dimension f(i2)
  43. dimension f(2)
  44. **** dimension fn(i2)
  45. dimension fn(2)
  46. **** dimension H66(i6,i6)
  47. dimension H66(6,6)
  48. c matrix of elasticity (sometimes denoted D)
  49. **** dimension x(i2)
  50. dimension x(2)
  51. c solution to be found (and initial guess at the entry of the subr.)
  52. **** dimension xnew(i2)
  53. dimension xnew(2)
  54. **** dimension esigmae(i6)
  55. dimension esigmae(6)
  56. c trial stress vector, = elastic predictor
  57. **** dimension trav1(i1),trav2(i2),trav2b(i2),trav2c(i2),
  58. **** trav6(i6),trav22(i2,i2)
  59. dimension trav1(6),trav2(2),trav2b(2),trav2c(2),
  60. . trav6(6),trav22(2,2)
  61. **** dimension P2(i6,i6),U(i6,i6),pi2(i6),v2loc(i2,i2)
  62. dimension P2(6,6),U(6,6),pi2(6),v2loc(2,2)
  63. dimension deps(i6)
  64. dimension parahot3(idimpara3),rk(2),esigmaf(6),dsiglc(6)
  65. dimension B11(1),B12(1),B21(1),B22(1),xtest(2),dsiglt(6)
  66.  
  67. DATA itermax /90/, tolx /1.e-4/
  68.  
  69.  
  70. *****
  71. * MESSAGES D'ERREUR ( SUPPRESSION DES AUTOMATIC OBJECTS)
  72. IF(I1.GT.1) PRINT *, ' ALGOPLAS2 - ERREUR I1 = ', I1, ' > 1 '
  73. IF(I2.GT.2) PRINT *, ' ALGOPLAS2 - ERREUR I2 = ', I2, ' > 2 '
  74. IF(I6.GT.6) PRINT *, ' ALGOPLAS2 - ERREUR I6 = ', I6, ' > 6 '
  75. *****
  76.  
  77.  
  78. i7 = 7
  79. i8 = 8
  80. i9 = 9
  81. i10=10
  82. i11=11
  83. i12=12
  84. i13=13
  85. i14=14
  86. i15=15
  87. i16=16
  88. i17=17
  89. i18=18
  90. i19=19
  91. i20=20
  92. i21=21
  93. i22=22
  94. i23=23
  95. i24=24
  96. i25=25
  97. i26=26
  98. i27=27
  99. i28=28
  100. i29=29
  101. i30=30
  102. i31=31
  103. i32=32
  104. i33=33
  105. i34=34
  106. i35=35
  107. i36=36
  108. i37=37
  109. i38=38
  110. i39=39
  111.  
  112. r1 = 1.
  113. r2 = 2.
  114. r0 = 0.
  115. r10 = 10.
  116. i0 = 0
  117. i1 = 1
  118.  
  119. lerror = .false.
  120. lproblem = .false.
  121. G = Eo/(r2*(r1+poison))
  122. rk(1) = rkt
  123. rk(2) = rkc
  124.  
  125. 87 continue
  126.  
  127. c initialization of fn
  128. FN(i1) = r0
  129. FN(i2) = r0
  130.  
  131. c initialization
  132. do iloc=1,6
  133. esigmaf(iloc) = esigmae(iloc)
  134. enddo
  135. rkappat = rkappait
  136. rkappac = rkappaic
  137. test = 0.d0
  138. xnew(1) = 0.0d0
  139. xnew(2) = 0.0d0
  140.  
  141. c calculate a first f[x(1),x(2),esigmae,rkappait,rkappaic]
  142. call find_f3d(f,esigmae,rkappait,rkappaic,x(i1),x(i2),alpha,
  143. . pi2,fc,ft,ag,fc0,rkappa1,ac,bc,H66,P2,lg,ntot,
  144. . i2,i6,fb,parahot3,idimpara3,lerror,lcp,esigmae,U)
  145. if (lerror) then
  146. return
  147. endif
  148.  
  149. if ((f(i1).lt.10.).and.(f(i2).lt.10.)) then
  150. return
  151. endif
  152.  
  153. if (f(1).lt.10.) rk(1) = 0.d0
  154.  
  155. fn(i1) = f(i1)
  156. fn(i2) = f(i2)
  157.  
  158. c =============================================================
  159. c ! !
  160. c ! initialization of JACOBIAN B !
  161. c ! ============================ !
  162. c =============================================================
  163.  
  164. c initialization of the Jacobian B by the following equation
  165. c B(11) = - (dft/dsigma dsigma/dlambdat - dtaut)
  166. c B(12) = - (rk(1)*rk(2)*dft/dsigma dsigma/dlambdac)
  167. c B(21) = - (rk(2)*rk(1)*dfc/dsigma dsigmaLdlambdat)
  168. c B(22) = - (dfc/dsigma dsigmaLdlambdac - (1-alpha)*dtauc)
  169.  
  170. call dsigRank(esigmaf,dft,i3,i6,parahot3,idimpara3,
  171. . rkappait,lcp)
  172. call dsigDrucker(fc,fb,esigmaf,dfc,alpha,P2,pi2,i1,i6,lcp)
  173.  
  174. call dsiglambda(xnew(1),xnew(2),P2,pi2,ag,fc,fb,lcp,rkappait,
  175. . parahot3,idimpara3,U,H66,esigmaf,dsiglt,dsiglc)
  176.  
  177. call mulATB(dft,dsiglt,B11,6,1,6,1)
  178. call mulATB(dft,dsiglc,B12,6,1,6,1)
  179. call mulATB(dfc,dsiglt,B21,6,1,6,1)
  180. call mulATB(dfc,dsiglc,B22,6,1,6,1)
  181.  
  182. dtaut = dfhardct(rkappat)
  183. dtauc = dfhardcc(rkappac,fc,fc0,rkappa1,ac,bc)
  184.  
  185. B(1,1) = B11(1) - dtaut
  186. B(1,2) = rk(1) * rk(2) * B12(1)
  187. B(2,1) = rk(2) * rk(1) * B21(1)
  188. B(2,2) = B22(1) - (1.-alpha) * dtauc
  189.  
  190. dtmB = B(i1,i1)*B(i2,i2)-B(i1,i2)*B(i2,i1)
  191. if (dtmB.eq.r0) then
  192. lerror = .true.
  193. print *,'ERROR BROYDEN dtmB=0'
  194. return
  195. else
  196. Binv(i1,i1) = B(i2,i2)/dtmB
  197. Binv(i1,i2) =-B(i1,i2)/dtmB
  198. Binv(i2,i1) =-B(i2,i1)/dtmB
  199. Binv(i2,i2) = B(i1,i1)/dtmB
  200. endif
  201.  
  202. c =================================================================
  203. c ! !
  204. c ! loop "Newton-Raphson !
  205. c ! ==================== !
  206. c =================================================================
  207.  
  208. do iter=i1,itermax
  209.  
  210. c **************************
  211. c ! CALCULATION OF DE DX_i !
  212. c **************************
  213. c -1
  214. c DX_i = - J_i F(X_i)
  215. c
  216. TRAV22(i1,i1) = -BINV(i1,i1)
  217. TRAV22(i1,i2) = -BINV(i1,i2)
  218. TRAV22(i2,i1) = -BINV(i2,i1)
  219. TRAV22(i2,i2) = -BINV(i2,i2)
  220. call mulab(trav22,F,DX,i2,i2,i2,i1)
  221. if (ABS(dx(i1)) .lt. 1.E-18) dx(i1) = r0
  222. if (ABS(dx(i2)) .lt. 1.E-18) dx(i2) = r0
  223. if (((ABS(dx(i1)) .gt. 1.).and.(rkt.eq.1.)).or.
  224. . ((ABS(dx(i2)) .gt. 1.).and.(rkc.eq.1.))) then
  225. lnoconv = .true.
  226. lerror = .true.
  227. return
  228. endif
  229. if ((ABS(dx(i1)).gt.1.).and.(rkt.eq.0.))
  230. . dx(1) = 1. * dx(1)/ABS(dx(i1))
  231. if ((ABS(dx(i2)).gt.1.).and.(rkc.eq.0.))
  232. . dx(2) = 1. * dx(2)/ABS(dx(i2))
  233.  
  234. c *************************
  235. c ! CALCULATION OF X_i+1 !
  236. c *************************
  237. c
  238. c X_i+1 = X_i + DX_i
  239. c
  240. do iloc=1,2
  241. Xtest(iloc) = X(iloc) + DX(iloc)
  242. if ((Xtest(iloc).le.(0.)).and.(rk(iloc).eq.1)) then
  243. do jloc=1,2
  244. if ((rk(jloc).eq.1).and.(jloc.ne.iloc))
  245. . DX(jloc) = - f(jloc)*BINV(jloc,jloc)
  246. enddo
  247. endif
  248. enddo
  249.  
  250. XNEW(i1) = X(i1) + DX(i1)
  251. XNEW(i2) = X(i2) + DX(i2)
  252. XNEW(i1) = max(0.D0,XNEW(i1))
  253. XNEW(i2) = max(0.D0,XNEW(i2))
  254.  
  255. do iloc=1,2
  256. if (XNEW(iloc).eq.(0.)) then
  257. rk(iloc) = 0.
  258. else
  259. rk(iloc) = 1.
  260. endif
  261. enddo
  262.  
  263. c The strain hardening parameters can only grow, they cannot decrease.
  264. if ((xnew(i1).le.r0).and.(xnew(i2).le.r0)) then
  265. lnoconv = .true.
  266. lerror = .true.
  267. return
  268. endif
  269.  
  270. c **********************************************
  271. c ! UPDATE OF THE HARDENING PARAMETERS !
  272. c **********************************************
  273.  
  274. X(i1) = XNEW(i1)
  275. X(i2) = XNEW(i2)
  276. c update of the hardening parameters
  277. rkappait= parahot3(idimpara3-37)
  278. rkappaic= parahot3(idimpara3-36)
  279. rkappat = rkappait+x(1)
  280. rkappac = rkappaic+x(2)
  281.  
  282. c ******************************
  283. c ! UPDATE OF THE STRESSES !
  284. c ******************************
  285.  
  286. call sigf3D_implicit(esigmae,xnew(1),xnew(2),esigmaf,H66,P2,
  287. . pi2,ag,i6,fc,fb,parahot3,idimpara3,lerror,lcp,esigmae,U)
  288. if (lerror) then
  289. return
  290. endif
  291.  
  292. c **************************
  293. c ! YIELD FUNCTIONS !
  294. c **************************
  295. c
  296. c calculate f[0.d0,0.d0,esigmaf,rkappat,rkappac]
  297. call find_f3d(fn,esigmaf,rkappat,rkappac,0.d0,0.d0,
  298. . alpha,pi2,fc,ft,ag,fc0,rkappa1,ac,bc,H66,P2,lg,ntot,
  299. . i2,i6,fb,parahot3,idimpara3,lerror,lcp,esigmae,U)
  300. if (lerror) then
  301. return
  302. endif
  303.  
  304. if (fn(i1).gt.r10) rk(1)=1.
  305. if (fn(i2).gt.r10) rk(2)=1.
  306.  
  307. rkt = rk(1)
  308. rkc = rk(2)
  309.  
  310. F(i1) = FN(i1)
  311. F(i2) = FN(i2)
  312.  
  313. c *************************************
  314. c ! UPDATE OF THE INVERT JACOBIAN !
  315. c *************************************
  316.  
  317. call dsigRank(esigmaf,dft,i3,i6,parahot3,idimpara3,
  318. . rkappait,lcp)
  319. call dsigDrucker(fc,fb,esigmaf,dfc,alpha,P2,pi2,i1,i6,lcp)
  320.  
  321. call dsiglambda(xnew(1),xnew(2),P2,pi2,ag,fc,fb,lcp,rkappait,
  322. . parahot3,idimpara3,U,H66,esigmaf,dsiglt,dsiglc)
  323.  
  324. call mulATB(dft,dsiglt,B11,6,1,6,1)
  325. call mulATB(dft,dsiglc,B12,6,1,6,1)
  326. call mulATB(dfc,dsiglt,B21,6,1,6,1)
  327. call mulATB(dfc,dsiglc,B22,6,1,6,1)
  328.  
  329. dtaut = dfhardct(rkappat)
  330. dtauc = dfhardcc(rkappac,fc,fc0,rkappa1,ac,bc)
  331.  
  332. B(1,1) = B11(1) - dtaut
  333. B(1,2) = rk(1) * rk(2) * B12(1)
  334. B(2,1) = rk(2) * rk(1) * B21(1)
  335. B(2,2) = B22(1) - (1.-alpha) * dtauc
  336.  
  337. dtmB = B(i1,i1)*B(i2,i2)-B(i1,i2)*B(i2,i1)
  338. if (dtmB.eq.r0) then
  339. lerror = .true.
  340. print *,'ERROR BROYDEN dtmB=0'
  341. return
  342. else
  343. Binv(i1,i1) = B(i2,i2)/dtmB
  344. Binv(i1,i2) =-B(i1,i2)/dtmB
  345. Binv(i2,i1) =-B(i2,i1)/dtmB
  346. Binv(i2,i2) = B(i1,i1)/dtmB
  347. endif
  348.  
  349. c ******************************
  350. c ! TEST END OF THE LOOP !
  351. c ******************************
  352.  
  353. test = r0
  354. if (iter.ne.i1) then
  355. c The plasticity functions must be satisfied with a precision of 10 Pa, i.e. 0.000010 N/m²
  356. if ((ABS(rk(1)*fn(i1)).le.10.).and.
  357. . (ABS(rk(2)*fn(i2)).le.10.)) then
  358. c convergence has been obtained
  359. return
  360. endif
  361.  
  362. do k=1,2
  363. if (dx(k).ne.r0) then
  364. temp = abs(rk(k)*dx(k)/(1.E-10))
  365. if (temp.gt.test) test = temp
  366. endif
  367. end do
  368. if (test.lt.tolx) then
  369. if ((ABS(rk(1)*fn(i1)).le.100.).and.
  370. . (ABS(rk(2)*fn(i2)).le.100.)) then
  371. return
  372. else
  373. lnoconv = .true.
  374. return
  375. endif
  376. endif
  377. endif
  378. c END OF THE LOOP EITHER BECAUSE THE EQUATION IS SOLVED
  379. c OR BECAUSE THE X VARIABLE DOES NOT VARY ANY MORE
  380. c OTHERWISE, GO ON WITH THE ITERATIONS
  381. c ************************************
  382.  
  383. F(i1) = FN(i1)
  384. F(i2) = FN(i2)
  385.  
  386. end do
  387. write(2,1) test,f(i1),f(i2),ntot
  388. 1 format( ' itermax exceeded in Broyden, test:',e7.2,
  389. . ', ft:',E8.2,' N/m², fc:',e8.2,' N/m²'/
  390. . ' ntot = ',i6)
  391. lerror = .true.
  392. lnoconv = .true.
  393.  
  394. RETURN
  395. END
  396.  
  397.  
  398.  
  399.  
  400.  
  401.  

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