Télécharger la plaquette Cast3M


fr en

Télécharger algop1cp.eso

Retour à la liste

Numérotation des lignes : 

algop1cp

  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