Télécharger algop1cp.eso

Retour à la liste

Numérotation des lignes :

  1. C ALGOP1CP SOURCE FANDEUR 14/03/25 21:15:04 7993
  2. Subroutine algoplas1cp(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 of Broyden (Newton-Raphson secante)
  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.  
  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. dimension esigmae(6)
  57. c trial stress vector, = elastic predictor
  58. **** dimension trav1(i1),trav2(i2),trav2b(i2),trav2c(i2),
  59. **** . trav6(i6),trav22(i2,i2)
  60. dimension trav1(1),trav2(2),trav2b(2),trav2c(2),
  61. . trav6(6),trav22(2,2)
  62. **** dimension P2(i6,i6),U(i6,i6),pi2(i6),v2loc(i2,i2)
  63. dimension P2(6,6),U(6,6),pi2(6),v2loc(2,2)
  64. **** dimension deps(i6),esigmain(6),dsigma(6),esigmaf(6)
  65. dimension deps(6),esigmain(6),dsigma(6),esigmaf(6)
  66. dimension parahot3(idimpara3)
  67.  
  68. DATA itermax /90/, tolx /1.e-4/
  69.  
  70. *****
  71. * MESSAGES D'ERREUR ( SUPPRESSION DES AUTOMATIC OBJECTS)
  72. IF(I1.GT.1) PRINT *, ' ALGOPLAS1CP - ERREUR I1 = ', I1, ' > 1 '
  73. IF(I2.GT.2) PRINT *, ' ALGOPLAS1CP - ERREUR I2 = ', I2, ' > 2 '
  74. IF(I6.GT.6) PRINT *, ' ALGOPLAS1CP - ERREUR I6 = ', I6, ' > 6 '
  75. *****
  76.  
  77. i7 = 7
  78. i8 = 8
  79. i9 = 9
  80. i10=10
  81. i11=11
  82. i12=12
  83. i13=13
  84. i14=14
  85. i15=15
  86. i16=16
  87. i17=17
  88. i18=18
  89. i19=19
  90. i20=20
  91. i21=21
  92. i22=22
  93. i23=23
  94. i24=24
  95. i25=25
  96. i26=26
  97. i27=27
  98. i28=28
  99. i29=29
  100. i30=30
  101. i31=31
  102. i32=32
  103. i33=33
  104. i34=34
  105. i35=35
  106. i36=36
  107. i37=37
  108. i38=38
  109. i39=39
  110.  
  111. r1 = 1.
  112. r0 = 0.
  113. r2 = 2.
  114. r10 = 10.
  115.  
  116. lerror = .false.
  117. lproblem = .false.
  118. G = Eo/(r2*(r1+poison))
  119. lfirstiter = .true.
  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 initialization counter loop plane stress
  132. iderin = 0
  133. iter = 0
  134.  
  135. c initialization
  136. do iloc=1,6
  137. esigmain(iloc) = esigmae(iloc)
  138. esigmaf(iloc) = esigmae(iloc)
  139. enddo
  140. rkappat = rkappait
  141. rkappac = rkappaic
  142. test = 0.d0
  143. dxprevt = 0.0d0
  144. dxprevc = 0.0d0
  145.  
  146. c calculate a first f[x(1),x(2),esigmae,rkappait,rkappaic]
  147. call find_f3d(f,esigmae,rkappait,rkappaic,x(i1),x(i2),alpha,
  148. . pi2,fc,ft,ag,fc0,rkappa1,ac,bc,H66,P2,lg,ntot,
  149. . i2,i6,fb,parahot3,idimpara3,lerror,lcp,esigmae,U)
  150. if (lerror) then
  151. return
  152. endif
  153.  
  154. if ((f(i1).lt.10.).and.(f(i2).lt.10.)) then
  155. return
  156. endif
  157. if (f(1).lt.10.) rkt = 0.d0
  158.  
  159. f(i1) = rkt*f(i1)
  160. f(i2) = rkc*f(i2)
  161.  
  162. if (f(i1).gt.r0) lftact = .true.
  163. if (f(i2).gt.r0) lfcact = .true.
  164.  
  165. fn(i1) = f(i1)
  166. fn(i2) = f(i2)
  167.  
  168. c =================================================================
  169. c ! !
  170. c ! loop "Newton-Raphson !
  171. c ! ==================== !
  172. c =================================================================
  173.  
  174. 11 iter=iter+i1
  175. 10 continue
  176.  
  177. if (iter.lt.2) then
  178.  
  179. c ===================================================================
  180. c ! !
  181. c ! initialization of JACOBIAN B !
  182. c ! ============================ !
  183. c ===================================================================
  184.  
  185. c initialization of Jacobian B in x = (0;0) by the following equation
  186. c B(11) = - (dft/dsigma H66 dft/dsigma + dtaut)
  187. c B(12) = - (dft/dsigma H66 dfc/dsigma)
  188. c B(21) = - (dfc/dsigma H66 dft/dsigma)
  189. c B(22) = - (dfc/dsigma H66 dfc/dsigma + dtauc)
  190.  
  191. call dsigRank(esigmaf,dft,i3,i6,parahot3,idimpara3,
  192. . rkappait,lcp)
  193. c ********
  194. do iloc=i1,i6
  195. DFT(iloc) = rkt * DFT(iloc)
  196. enddo
  197.  
  198. call dsigDrucker(fc,fb,esigmaf,dfc,alpha,P2,pi2,i1,i6,lcp)
  199. c ***********
  200. do iloc=i1,i6
  201. DFC(iloc) = rkc * DFC(iloc)
  202. enddo
  203.  
  204. call mulATB(dft,H66,trav6,i6,i1,i6,i6)
  205. c ******
  206. call mulAB(trav6,dft,trav1,i1,i6,i6,i1)
  207. c *****
  208. df11 = trav1(i1)
  209.  
  210. call mulAB(trav6,dfc,trav1,i1,i6,i6,i1)
  211. c *****
  212. df12 = trav1(i1)
  213. df21 = df12
  214.  
  215. call mulATB(dfc,H66,trav6,i6,i1,i6,i6)
  216. c ******
  217. call mulAB(trav6,dfc,trav1,i1,i6,i6,i1)
  218. c *****
  219. df22 = trav1(i1)
  220.  
  221. c derivatives of the hardening functions
  222. dtaut = dfhardct(rkappat)
  223. dtauc = dfhardcc(rkappac,fc,fc0,rkappa1,ac,bc)
  224.  
  225. B(i1,i1) = -(df11+rkt*dtaut) +(r1-rkt)
  226. B(i1,i2) = -df12
  227. B(i2,i1) = -df21
  228. B(i2,i2) = -(df22+rkc*dtauc) +(r1-rkc)
  229.  
  230. c ===================================================================
  231. c ! !
  232. c ! invert the initial JACOBIAN !
  233. c ! =========================== !
  234. c ===================================================================
  235.  
  236. dtmB = B(i1,i1)*B(i2,i2)-B(i1,i2)*B(i2,i1)
  237. if (dtmB.eq.r0) then
  238. lerror = .true.
  239. return
  240. else
  241. Binv(i1,i1) = B(i2,i2)/dtmB
  242. Binv(i1,i2) =-B(i1,i2)/dtmB
  243. Binv(i2,i1) =-B(i2,i1)/dtmB
  244. Binv(i2,i2) = B(i1,i1)/dtmB
  245. endif
  246.  
  247. else
  248.  
  249. c **********************************************************
  250. c ! UPDATE THE INVERT JACOBIAN BY PAR SHERMAN-MORRISON !
  251. c **********************************************************
  252.  
  253. c ! YIELD FUNCTIONS !
  254. c
  255. c calculate f[0.d0,0.d0,esigmaf,rkappat,rkappac]
  256. call find_f3d(fn,esigmaf,rkappat,rkappac,0.d0,0.d0,
  257. . alpha,pi2,fc,ft,ag,fc0,rkappa1,ac,bc,H66,P2,lg,ntot,
  258. . i2,i6,fb,parahot3,idimpara3,lerror,lcp,esigmae,U)
  259. if (lerror) then
  260. return
  261. endif
  262.  
  263. c Reactivation of yield surface
  264. if (fn(i1).gt.r10) rkt=1.
  265. if (fn(i2).gt.r10) rkc=1.
  266.  
  267. fn(i1) = rkt*fn(i1)
  268. fn(i2) = rkc*fn(i2)
  269.  
  270. rloc = xnew(i1)
  271. xnew(i1) = rkt*rloc
  272. c xnew(1) = dlambdat
  273. rloc = xnew(i2)
  274. xnew(i2) = rkc*rloc
  275. c xnew(2) = dlambdac
  276.  
  277. DF(i1) = FN(i1) - F(i1)
  278. DF(i2) = FN(i2) - F(i2)
  279.  
  280. c calculation of the new Jacobian Binv_i+1
  281. c Now, we apply the formula of Sherman-Morrison
  282. c Binv.df
  283. call mulab(Binv,DF,TRAV2,i2,i2,i2,i1)
  284. c dx - Binv.df
  285. TRAV2b(i1) = DX(i1)-TRAV2(i1)
  286. TRAV2b(i2) = DX(i2)-TRAV2(i2)
  287. c dx.Binv
  288. call mulATB(dx,Binv,trav2c,i2,i1,i2,i2)
  289. c ******
  290. c (dx-Binv.df).(dx.Binv)
  291. call mulAB(trav2b,trav2c,trav22,i2,i1,i1,i2)
  292. c *****
  293. c dx.Binv.df
  294. call mulATB(dx,trav2,trav1,i2,i1,i2,i1)
  295. c ******
  296. rloc = r1/trav1(i1)
  297. do jloc=i1,i2
  298. do iloc=i1,i2
  299. TRAV22(iloc,jloc) = rloc*TRAV22(iloc,jloc)
  300. v2loc(iloc,jloc) = BINV(iloc,jloc)
  301. BINV(iloc,jloc) = v2loc(iloc,jloc) + TRAV22(iloc,jloc)
  302. end do
  303. end do
  304.  
  305. X(i1) = XNEW(i1)
  306. X(i2) = XNEW(i2)
  307.  
  308. F(i1) = FN(i1)
  309. F(i2) = FN(i2)
  310.  
  311. endif
  312.  
  313. c ******************************
  314. c ! CALCULATION OF DX_i !
  315. c ******************************
  316. c -1
  317. c DX_i = - J_i F(X_i)
  318. c
  319. TRAV22(i1,i1) = -BINV(i1,i1)
  320. TRAV22(i1,i2) = -BINV(i1,i2)
  321. TRAV22(i2,i1) = -BINV(i2,i1)
  322. TRAV22(i2,i2) = -BINV(i2,i2)
  323. call mulab(trav22,F,DX,i2,i2,i2,i1)
  324. if (abs(dx(i1)) .lt. 1.E-16) dx(i1) = r0
  325. if (abs(dx(i2)) .lt. 1.E-16) dx(i2) = r0
  326. if (((abs(dx(i1)) .gt. 1.).and.(rkt.eq.1.)).or.
  327. . ((abs(dx(i2)) .gt. 1.).and.(rkc.eq.1.))) then
  328. lnoconv = .true.
  329. lerror = .true.
  330. return
  331. endif
  332. if ( (abs(dx(i1)) .gt. 1.).and.(rkt.eq.0.) )
  333. . dx(1) = 1. * dx(1)/abs(dx(i1))
  334. if ( (abs(dx(i2)) .gt. 1.).and.(rkc.eq.0.) )
  335. . dx(2) = 1. * dx(2)/abs(dx(i2))
  336.  
  337. if (iter.gt.2) then
  338. r_z = 1.50d0 * dxprevt
  339. if (abs(dx(1)).gt.abs(r_z)) then
  340. rloc = dx(1)
  341. dx(1) = r_z * (rloc/abs(rloc))
  342. endif
  343. r_z = 1.50d0 * dxprevc
  344. if (abs(dx(2)).gt.abs(r_z)) then
  345. rloc = dx(2)
  346. dx(2) = r_z * (rloc/abs(rloc))
  347. endif
  348. endif
  349.  
  350. dxprevt = abs(dx(1))
  351. dxprevc = abs(dx(2))
  352.  
  353. c *****************************
  354. c ! CALCULATION OF X_i+1 !
  355. c *****************************
  356. c
  357. c X_i+1 = X_i + DX_i
  358. c
  359. XNEW(i1) = X(i1) + DX(i1)
  360. XNEW(i2) = X(i2) + DX(i2)
  361.  
  362. c The strain hardening parameters can only grow, they cannot decrease.
  363.  
  364. if ((xnew(i1).le.r0).and.(xnew(i2).le.r0)) then
  365. lnoconv = .true.
  366. lerror = .true.
  367. return
  368. endif
  369.  
  370. if (xnew(i1).eq.r0) then
  371. rkt = r0
  372. endif
  373. if (xnew(i2).eq.r0) then
  374. rkc = r0
  375. endif
  376. c new calculation of dx if we have reactivated a yield surface
  377. if ((xnew(i1).lt.r0).and.(rkc.eq.1.)) then
  378. rkt = r0
  379. if (lfirstiter) then
  380. x(1) = 0.d0
  381. x(2) = 0.d0
  382. lfirstiter = .false.
  383. go to 87
  384. else if (iter.lt.70) then
  385. go to 11
  386. else
  387. lerror = .true.
  388. lnoconv = .true.
  389. return
  390. endif
  391. endif
  392. if ((xnew(i2).lt.r0).and.(rkt.eq.1.)) then
  393. rkc = r0
  394. if (lfirstiter) then
  395. x(1) = 0.d0
  396. x(2) = 0.d0
  397. lfirstiter = .false.
  398. go to 87
  399. else if (iter.lt.70) then
  400. go to 11
  401. else
  402. lerror = .true.
  403. lnoconv = .true.
  404. return
  405. endif
  406. endif
  407.  
  408. if (xnew(i1).gt.r0) rkt=1.
  409. if (xnew(i2).gt.r0) rkc=1.
  410.  
  411. rloc = xnew(i1)
  412. xnew(i1) = rkt*rloc
  413. c xnew(1) = dlambdat
  414. rloc = xnew(i2)
  415. xnew(i2) = rkc*rloc
  416. c xnew(2) = dlambdac
  417.  
  418. c dxcor(rected) for calculation of dsigma
  419. c dx = xnew-x because if xnew<0 => xnew=0 and the dx must be corrected
  420. dxcor1 = xnew(1) - x(1)
  421. dxcor2 = xnew(2) - x(2)
  422.  
  423. c **************************
  424. c ! CALCULATION OF DSIGMA !
  425. c **************************
  426. c
  427. call sigf3D_implicit(esigmain,dxcor1,dxcor2,esigmaf,
  428. . H66,P2,pi2,ag,i6,
  429. . fc,fb,parahot3,idimpara3,lerror,lcp,esigmae,U)
  430. if (lerror) then
  431. return
  432. endif
  433. do i=1,6
  434. dsigma(i) = esigmaf(i) - esigmain(i)
  435. enddo
  436.  
  437. c ***************************************
  438. c ! PLANE STRESS ITERATIONS !
  439. c ***************************************
  440. c
  441. bb = abs(dsigma(3))
  442. esigmaf(3) = esigmain(3)-dsigma(3)
  443. esigmaf(1) = esigmain(1)-dsigma(3)*poison/(1.d0-poison)
  444. esigmaf(2) = esigmain(2)-dsigma(3)*poison/(1.d0-poison)
  445. if (bb.gt.(10.d0)) then
  446. if (iderin.eq.0) then
  447. iter = iter-1
  448. endif
  449. iderin = iderin+1
  450. if (iderin.gt.50) then
  451. write(6,*) ' problem internal iterations plane stress'
  452. return
  453. endif
  454. go to 10
  455. else
  456. continue
  457. endif
  458.  
  459. c *************************************
  460. c ! UPDATE OF THE STRESSES !
  461. c *************************************
  462. c
  463. do i=1,6
  464. esigmaf(i) = esigmaf(i) + dsigma(i)
  465. end do
  466. esigmaf(3) = 0.d0
  467. do i=1,6
  468. esigmain(i) = esigmaf(i)
  469. end do
  470.  
  471. c *************************************************
  472. c ! UPDATE OF THE HARDENING VARIABLES !
  473. c *************************************************
  474. X(i1) = XNEW(i1)
  475. X(i2) = XNEW(i2)
  476.  
  477. c update of the hardening parameters
  478. rkappait= parahot3(idimpara3-37)
  479. rkappaic= parahot3(idimpara3-36)
  480. rkappat = rkappait+x(1)
  481. rkappac = rkappaic+x(2)
  482.  
  483. c *************************
  484. c ! YIELD FUNCTIONS !
  485. c *************************
  486. c
  487. call find_f3d(fn,esigmaf,rkappat,rkappac,0.d0,0.d0,
  488. . alpha,pi2,fc,ft,ag,fc0,rkappa1,ac,bc,H66,P2,lg,ntot,
  489. . i2,i6,fb,parahot3,idimpara3,lerror,lcp,esigmae,U)
  490. if (lerror) then
  491. return
  492. endif
  493. if (fn(i1).gt.r10) rkt=1.
  494. if (fn(i2).gt.r10) rkc=1.
  495.  
  496. fn(i1) = rkt*fn(i1)
  497. fn(i2) = rkc*fn(i2)
  498.  
  499. c **********************************
  500. c ! TEST TO END THE LOOP !
  501. c **********************************
  502. if ((abs(dx(i1)).lt.1.E-10).and.(abs(dx(i2)).lt.1.E-10).and.
  503. & (abs(fn(i1)).le.10.).and.(abs(fn(i2)).le.10.)) then
  504. return
  505. endif
  506. test = r0
  507. if (iter.ne.i1) then
  508. c The plasticity functions must be satisfied with a precision of 10 Pa, i.e. 0.000010 N/m²
  509. if ((abs(fn(i1)).le.10.).and.(abs(fn(i2)).le.10.)) then
  510. c convergence has been obtained
  511. return
  512. endif
  513.  
  514. do k=1,2
  515. if (dx(k).ne.r0) then
  516. temp = abs(dx(k)/(1.E-10))
  517. if (temp.gt.test) test = temp
  518. endif
  519. end do
  520. if (test.lt.tolx) then
  521. return
  522. endif
  523. endif
  524. c END OF THE LOOP EITHER BECAUSE THE EQUATION IS SOLVED
  525. c OR BECAUSE THE X VARIABLE DOES NOT VARY ANY MORE
  526. c OTHERWISE, GO ON WITH THE ITERATIONS
  527. c ************************************
  528.  
  529. if (iter.lt.70) then
  530. go to 11
  531. else
  532. lerror = .true.
  533. lnoconv = .true.
  534. endif
  535.  
  536. RETURN
  537. END
  538.  
  539.  
  540.  

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