Télécharger algop1.eso

Retour à la liste

Numérotation des lignes :

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

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