Télécharger algop2cp.eso

Retour à la liste

Numérotation des lignes :

  1. C ALGOP2CP SOURCE CB215821 16/04/21 21:15:07 8920
  2. Subroutine Algoplas2cp(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,esigmaf,rkappat,rkappac,
  6. . lnoconv,lcp,U)
  7.  
  8. c SUBROUTINE added by THG
  9.  
  10. c THIS SUBROUTINE IS CALLED BY DPRAN3D
  11. c the system is solved by the method Newton
  12.  
  13. c rkt x function_t(x(1),x(2)) + (1-rkt) x(1) = 0
  14. c rkc x function_c(x(1),x(2)) + (1-rkc) x(2) = 0
  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(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),esigmain(6),dsigma(6),esigmaf(6)
  64. dimension parahot3(idimpara3),rk(2),xtest(2),dsiglt(6),dsiglc(6)
  65. dimension B11(1),B12(1),B21(1),B22(1)
  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 *, ' ALGOPLAS2CP - ERREUR I1 = ', I1, ' > 1 '
  73. IF(I2.GT.2) PRINT *, ' ALGOPLAS2CP - ERREUR I2 = ', I2, ' > 2 '
  74. IF(I6.GT.6) PRINT *, ' ALGOPLAS2CP - 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. r0 = 0.
  114. r2 = 2.
  115. r10 = 10.
  116.  
  117. lerror = .false.
  118. lproblem = .false.
  119. lonechance = .true.
  120. G = Eo/(r2*(r1+poison))
  121. rk(1) = rkt
  122. rk(2) = rkc
  123.  
  124. c initialization of fn
  125. FN(i1) = r0
  126. FN(i2) = r0
  127.  
  128. c initialization counter loop plane stress
  129. iderin = 0
  130. iter = 0
  131.  
  132. c initialization
  133. do iloc=1,6
  134. esigmain(iloc) = esigmae(iloc)
  135. esigmaf(iloc) = esigmae(iloc)
  136. enddo
  137. rkappat = rkappait
  138. rkappac = rkappaic
  139. test = 0.d0
  140. xnew(1) = 0.0d0
  141. xnew(2) = 0.0d0
  142.  
  143. c calculate a first f[x(1),x(2),esigmae,rkappait,rkappaic]
  144. call find_f3d(f,esigmae,rkappait,rkappaic,x(i1),x(i2),alpha,
  145. . pi2,fc,ft,ag,fc0,rkappa1,ac,bc,H66,P2,lg,ntot,
  146. . i2,i6,fb,parahot3,idimpara3,lerror,lcp,esigmae,U)
  147. if (lerror) then
  148. return
  149. endif
  150.  
  151. if ((f(i1).lt.10.).and.(f(i2).lt.10.)) then
  152. return
  153. endif
  154. if (f(1).lt.10.) rk(1) = 0.d0
  155.  
  156. fn(i1) = f(i1)
  157. fn(i2) = f(i2)
  158.  
  159. c =================================================================
  160. c ! !
  161. c ! loop "Newton-Raphson !
  162. c ! ==================== !
  163. c =================================================================
  164.  
  165. 11 iter=iter+i1
  166. 10 continue
  167.  
  168.  
  169. c *************************************
  170. c ! CALCULATION OF THE JACOBIAN !
  171. c *************************************
  172.  
  173. call dsigRank(esigmaf,dft,i3,i6,parahot3,idimpara3,
  174. . rkappait,lcp)
  175. call dsigDrucker(fc,fb,esigmaf,dfc,alpha,P2,pi2,i1,i6,lcp)
  176.  
  177. call dsiglambda(xnew(1),xnew(2),P2,pi2,ag,fc,fb,lcp,rkappait,
  178. . parahot3,idimpara3,U,H66,esigmaf,dsiglt,dsiglc)
  179.  
  180. call mulATB(dft,dsiglt,B11,6,1,6,1)
  181. call mulATB(dft,dsiglc,B12,6,1,6,1)
  182. call mulATB(dfc,dsiglt,B21,6,1,6,1)
  183. call mulATB(dfc,dsiglc,B22,6,1,6,1)
  184.  
  185. dtaut = dfhardct(rkappat)
  186. dtauc = dfhardcc(rkappac,fc,fc0,rkappa1,ac,bc)
  187.  
  188. B(1,1) = B11(1) - dtaut
  189. B(1,2) = rk(1) * rk(2) * B12(1)
  190. B(2,1) = rk(2) * rk(1) * B21(1)
  191. B(2,2) = B22(1) - (1.-alpha) * dtauc
  192.  
  193. dtmB = B(i1,i1)*B(i2,i2)-B(i1,i2)*B(i2,i1)
  194. if (dtmB.eq.r0) then
  195. lerror = .true.
  196. return
  197. else
  198. Binv(i1,i1) = B(i2,i2)/dtmB
  199. Binv(i1,i2) =-B(i1,i2)/dtmB
  200. Binv(i2,i1) =-B(i2,i1)/dtmB
  201. Binv(i2,i2) = B(i1,i1)/dtmB
  202. endif
  203.  
  204. c *************************
  205. c ! CALCULATION OF DX_i !
  206. c *************************
  207. c -1
  208. c DX_i = - J_i F(X_i)
  209. c
  210. TRAV22(i1,i1) = -BINV(i1,i1)
  211. TRAV22(i1,i2) = -BINV(i1,i2)
  212. TRAV22(i2,i1) = -BINV(i2,i1)
  213. TRAV22(i2,i2) = -BINV(i2,i2)
  214. call mulab(trav22,F,DX,i2,i2,i2,i1)
  215. if (ABS(dx(i1)) .lt. 1.E-18) dx(i1) = r0
  216. if (ABS(dx(i2)) .lt. 1.E-18) dx(i2) = r0
  217. if (((ABS(dx(i1)) .gt. 1.).and.(rk(1).eq.1.)).or.
  218. . ((ABS(dx(i2)) .gt. 1.).and.(rk(2).eq.1.))) then
  219. lnoconv = .true.
  220. lerror = .true.
  221. return
  222. endif
  223. if ( (ABS(dx(i1)) .gt. 1.).and.(rk(1).eq.0.) )
  224. . dx(1) = 1. * dx(1)/ABS(dx(i1))
  225. if ( (ABS(dx(i2)) .gt. 1.).and.(rk(2).eq.0.) )
  226. . dx(2) = 1. * dx(2)/ABS(dx(i2))
  227.  
  228. c *************************
  229. c ! CALCULATION OF X_i+1 !
  230. c *************************
  231. c
  232. c X_i+1 = X_i + DX_i
  233. c
  234. do iloc=1,2
  235. Xtest(iloc) = X(iloc) + DX(iloc)
  236. if ((Xtest(iloc).le.(0.)).and.(rk(iloc).eq.1)) then
  237. do jloc=1,2
  238. if ((rk(jloc).eq.1).and.(jloc.ne.iloc))
  239. . DX(jloc) = - f(jloc)*BINV(jloc,jloc)
  240. enddo
  241. endif
  242. enddo
  243.  
  244. XNEW(i1) = X(i1) + DX(i1)
  245. XNEW(i2) = X(i2) + DX(i2)
  246.  
  247. XNEW(i1) = max(0.D0,XNEW(i1))
  248. XNEW(i2) = max(0.D0,XNEW(i2))
  249.  
  250. do iloc=1,2
  251. if (XNEW(iloc).eq.(0.)) then
  252. rk(iloc) = 0.
  253. else
  254. rk(iloc) = 1.
  255. endif
  256. enddo
  257.  
  258. c The strain hardening parameters can only grow, they cannot decrease.
  259. if ((xnew(i1).le.r0).and.(xnew(i2).le.r0)) then
  260. lnoconv = .true.
  261. lerror = .true.
  262. return
  263. endif
  264.  
  265. dxcor1 = xnew(1) - x(1)
  266. dxcor2 = xnew(2) - x(2)
  267.  
  268. c **************************
  269. c ! CALCULATION OF DSIGMA !
  270. c **************************
  271. c
  272. call sigf3D_implicit(esigmain,dxcor1,dxcor2,esigmaf,H66,
  273. . P2,pi2,ag,i6,fc,fb,parahot3,idimpara3,lerror,lcp,esigmae,U)
  274. if (lerror) then
  275. return
  276. endif
  277. do i=1,6
  278. dsigma(i) = esigmaf(i) - esigmain(i)
  279. enddo
  280.  
  281. c ***********************************
  282. c ! PLANE STRESS ITERATIONS !
  283. c ***********************************
  284. c
  285. bb = abs(dsigma(3)+esigmaf(3))
  286. esigmaf(3) = -dsigma(3)
  287. esigmaf(1) = esigmain(1)-dsigma(3)*poison/(1.d0-poison)
  288. esigmaf(2) = esigmain(2)-dsigma(3)*poison/(1.d0-poison)
  289. if (bb.gt.(10.d0)) then
  290. if (iderin.eq.0) then
  291. iter = iter-1
  292. endif
  293. iderin = iderin+1
  294. if (iderin.gt.50) then
  295. write(6,*) ' problem in internal iterations plane stress'
  296. return
  297. endif
  298. go to 10
  299. else
  300. continue
  301. endif
  302.  
  303. c *********************************
  304. c ! UPDATE OF THE STRESSES !
  305. c *********************************
  306. c
  307. do i=1,6
  308. esigmaf(i) = esigmaf(i) + dsigma(i)
  309. end do
  310. esigmaf(3) = 0.d0
  311. do i=1,6
  312. esigmain(i) = esigmaf(i)
  313. end do
  314.  
  315. c *********************************************
  316. c ! UPDATE OF THE HARDENING PARAMETERS !
  317. c *********************************************
  318. X(i1) = XNEW(i1)
  319. X(i2) = XNEW(i2)
  320.  
  321. c update of the hardening parameters
  322. rkappait= parahot3(idimpara3-37)
  323. rkappaic= parahot3(idimpara3-36)
  324. rkappat = rkappait+x(1)
  325. rkappac = rkappaic+x(2)
  326.  
  327. c ****************************
  328. c ! YIELD FUNCTIONS !
  329. c ****************************
  330. c
  331. c calculate f[0.d0,0.d0,esigmaf,rkappat,rkappac]
  332. call find_f3d(fn,esigmaf,rkappat,rkappac,0.d0,0.d0,
  333. . alpha,pi2,fc,ft,ag,fc0,rkappa1,ac,bc,H66,P2,lg,ntot,
  334. . i2,i6,fb,parahot3,idimpara3,lerror,lcp,esigmae,U)
  335. if (lerror) then
  336. return
  337. endif
  338. if (fn(i1).gt.r10) rk(1)=1.
  339. if (fn(i2).gt.r10) rk(2)=1.
  340.  
  341. rkt = rk(1)
  342. rkc = rk(2)
  343.  
  344. X(i1) = XNEW(i1)
  345. X(i2) = XNEW(i2)
  346.  
  347. F(i1) = FN(i1)
  348. F(i2) = FN(i2)
  349.  
  350. c *******************************
  351. c ! TEST TO END THE LOOP !
  352. c *******************************
  353.  
  354. test = r0
  355. if (iter.ne.i1) then
  356.  
  357. c The plasticity functions must be satisfied with a precision of 10 Pa, i.e. 0.000010 N/m²
  358. if ((ABS(rk(1)*fn(i1)).le.10.).and.
  359. . (ABS(rk(2)*fn(i2)).le.10.)) then
  360. c convergence has been obtained
  361. return
  362. endif
  363. c
  364. do k=1,2
  365. if (dx(k).ne.r0) then
  366. temp = abs(rk(k)*dx(k)/(1.E-10))
  367. if (temp.gt.test) test = temp
  368. endif
  369. end do
  370. if (test.lt.tolx) then
  371. if ((ABS(rk(1)*fn(i1)).le.100.).and.
  372. . (ABS(rk(2)*fn(i2)).le.100.)) then
  373. return
  374. else
  375. lnoconv = .true.
  376. return
  377. endif
  378. endif
  379. endif
  380. c END OF THE LOOP EITHER BECAUSE THE EQUATION IS SOLVED
  381. c OR BECAUSE THE X VARIABLE DOES NOT VARY ANY MORE
  382. c OTHERWISE, GO ON WITH THE ITERATIONS
  383. c ************************************
  384.  
  385. if (iter.lt.70) then
  386. go to 11
  387. else
  388. if ((ABS(rk(1)*fn(i1)).le.10000.).and.
  389. . (ABS(rk(2)*fn(i2)).le.10000.)
  390. . .and.(lonechance)) then
  391. iter = 2
  392. lonechance = .false.
  393. go to 11
  394. else
  395. lerror = .true.
  396. lnoconv = .true.
  397. endif
  398. endif
  399.  
  400. RETURN
  401. END
  402.  
  403.  
  404.  
  405.  
  406.  
  407.  

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