Télécharger la plaquette Cast3M


fr en

Télécharger modjfd.eso

Retour à la liste

Numérotation des lignes : 

modjfd

  1. C MODJFD    SOURCE    CB215821  17/11/30    21:16:52     9639           
  2.       SUBROUTINE MODJFD(WRK0,NMATT,EPSP1,EPSP2,YLIM1,YLIM2,Z1,Z2,Y1,Y2,
  3.      1                  JSIGN,SIG1,SIG2,CONV,DT)
  4. C
  5. C     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  6. C     %                                                               %
  7. C     %  INTEGRATION DU MODJFD UNILATERAL AVEC FERMETURE LINEAIRE     %
  8. C     %          DE FISSURES EN CONTRAINTES PLANES                    %
  9. C     %                (MODCLB NON ASSOCIE)                           %
  10. C     %      A L"AIDE DE LA METHODE DE NEWTON                         %
  11. C     %    AJOUT VICOENDOMMAGEMENT J.F. DUBE                          %
  12. C     %                                                               %
  13. C     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  14. C
  15. C
  16. C
  17. C        VARIABLES D"ENTREES
  18. C
  19. C
  20.       IMPLICIT INTEGER(I-N)
  21.       REAL*8 EPSP1,EPSP2
  22. C
  23. C        VARIABLES DE SORTIE
  24. C
  25. C
  26.       REAL*8 SIG1,SIG2,Z1,Z2,Y1,Y2,DT
  27.       INTEGER JSIGN(3)
  28.       LOGICAL CONV
  29. C
  30. C
  31. C        PARAMETRES D"EVOLUTION
  32. C
  33. C
  34.       REAL*8 E,NU,A1,A2,B1,B2,Y01,Y02,SIGF
  35.       REAL*8 BETA1,BETA2,YLIM1,YLIM2
  36. C
  37. C
  38. C        VARIABLES INTERNES
  39. C
  40. C
  41.       REAL*8 Y001,Y002,DY01,DY02,hdp1,hdp2
  42.       REAL*8 DELTASIG
  43.       INTEGER NDIM
  44.       PARAMETER (NDIM=6)
  45.       REAL*8 JACOB(NDIM,NDIM),J15,J16,J25,J26
  46.       REAL*8 DELAFON(NDIM)
  47.       REAL*8 X(NDIM),DX(NDIM),X0(NDIM),XX,XIG1,XIG2
  48.       REAL*8 ERREUR,TEST
  49.       REAL*8 H1,H2,HT,HT2,trf
  50.       REAL*8 ALFA1,ALFA2
  51.       REAL*8 SIGY
  52.       REAL*8 UN,DEUX
  53.       REAL*8 ZMAX
  54.       INTEGER ITER ,I,DIMPB
  55.       LOGICAL LDOM1,LDOM2,evol1,evol2
  56.       REAL*8 M1,M2,ALPHA1,ALPHA2
  57.       REAL*8 x3,x4,x5,x6,dx5,dx6
  58.       INTEGER NMATT
  59.       SEGMENT WRK0
  60.          REAL*8 XMAT(NMATT)
  61.       ENDSEGMENT
  62.       UN=1.D0
  63.       DEUX=2.D0
  64.       ZMAX=1.d+04
  65.       E     = XMAT(1)
  66.       NU    = XMAT(2)
  67.       Y001  = XMAT(5)
  68.       Y002  = XMAT(6)
  69.       A1    = XMAT(7)
  70.       A2    = XMAT(8)
  71.       B1    = XMAT(9)
  72.       B2    = XMAT(10)
  73.       BETA1 = XMAT(11)
  74.       BETA2 = XMAT(12)
  75.       SIGF  = XMAT(13)
  76.       M1    = XMAT(14)
  77.       M2    = XMAT(15)
  78.       ALPHA1= XMAT(16)
  79.       ALPHA2= XMAT(17)
  80.       SIGY  =SQRT(2.D0*E*Y001)
  81. *      DELTASIG=SIGY*1.D-4
  82.       DELTASIG=0.d0
  83. *      print*,'E,NU,Y01,Y02,A1,A2,B1,B2,BETA1,BETA2'
  84. *      print*,E,NU,Y01,Y02,A1,A2,B1,B2,BETA1,BETA2
  85.       TEST=1.D-12
  86.       ERREUR=1.D0
  87.       ITER=0
  88.       X(1)=SIG1
  89.       X(2)=SIG2
  90.       X(3)=Y1
  91.       X(4)=Y2
  92.       X(5)=Z1
  93.       X(6)=Z2
  94.       DO I=1, NDIM
  95.          X0(I)=X(I)
  96.       END DO
  97.       x5=Z1
  98.       x6=z2
  99. C
  100. C
  101. C
  102.       JACOB(1,1) = -1.d0
  103.       JACOB(1,2) = 0.D0
  104.       JACOB(1,3) = 0.D0
  105.       JACOB(1,4) = 0.D0
  106.       JACOB(1,5) = 0.D0
  107.       JACOB(1,6) = 0.D0
  108.       JACOB(2,1) = 0.D0
  109.       JACOB(2,2) = -1.d0
  110.       JACOB(2,3) = 0.D0
  111.       JACOB(2,4) = 0.D0
  112.       JACOB(2,5) = 0.D0
  113.       JACOB(2,6) = 0.D0
  114.       JACOB(3,1) = 0.D0
  115.       JACOB(3,2) = 0.D0
  116.       JACOB(3,3) = -1.d0
  117.       JACOB(3,4) = 0.D0
  118.       JACOB(3,5) = 0.D0
  119.       JACOB(3,6) = 0.D0
  120.       JACOB(4,1) = 0.D0
  121.       JACOB(4,2) = 0.D0
  122.       JACOB(4,3) = 0.D0
  123.       JACOB(4,4) = -1.d0
  124.       JACOB(4,5) = 0.D0
  125.       JACOB(4,6) = 0.D0
  126.       JACOB(5,1) = 0.D0
  127.       JACOB(5,2) = 0.D0
  128.       JACOB(5,3) = 0.D0
  129.       JACOB(5,4) = 0.D0
  130.       JACOB(5,5) = -1.d0
  131.       JACOB(5,6) = 0.D0
  132.       JACOB(6,1) = 0.D0
  133.       JACOB(6,2) = 0.D0
  134.       JACOB(6,3) = 0.D0
  135.       JACOB(6,4) = 0.D0
  136.       JACOB(6,5) = 0.D0
  137.       JACOB(6,6) = -1.d0
  138.  
  139.       DO WHILE (ERREUR .GT. TEST .AND. ITER .LE. 100)
  140.          ITER=ITER+1
  141.          IF (X(1) .GE. DELTASIG) THEN
  142.             H1=UN
  143.          ELSE IF (X(1) .LE. -DELTASIG) THEN
  144.             H1=0.D0
  145.          ELSE
  146.             H1=(X(1)+DELTASIG)/(DEUX*DELTASIG)
  147.          END IF
  148.          IF (X(2) .GE. DELTASIG) THEN
  149.             H2=UN
  150.          ELSE IF (X(2) .LE. -DELTASIG) THEN
  151.             H2=0.D0
  152.          ELSE
  153.             H2=(X(2)+DELTASIG)/(DEUX*DELTASIG)
  154.          END IF
  155. c
  156. c     smooth transition on ht ht2
  157. c
  158.          trf = x(1) + x(2)
  159.          IF(trf.GE.deltasig)THEN
  160.             HT=UN
  161.             HT2=0.D0
  162.          else if(trf.ge.(-deltasig)) then
  163.             ht = (trf+deltasig)/DEUX/deltasig
  164.             ht2 = UN-ht
  165.          else if(trf.ge.(-sigf+deltasig)) then
  166.             ht = 0.D+00
  167.             ht2 = UN
  168.          else if(trf.ge.(-sigf-deltasig)) then
  169.             ht = 0.D+00
  170. c           ht2 =UN- (trf-(sigf-deltasig))/(DEUX*sigf)
  171.             ht2 = (trf+sigf+deltasig)/DEUX/deltasig
  172.          else
  173.             ht = 0.D+00
  174.             ht2 = 0.D+00
  175.          endif
  176. c         print*,'DTjfd',dt
  177.          dy01=(dim(x(5),z1)/(x(5)*z1*dt))
  178.          hdp1=1.0
  179.          if (dy01 .gt. 0.0) hdp1=0.0
  180.          y01=y001+m1*dy01**alpha1
  181.          dy02=(dim(x(6),z2)/(x(6)*z2*dt))
  182.          hdp2=1.0
  183.          if (dy02 .gt. 0.0) hdp2=0.0
  184.          y02=y002+m2*dy02**alpha2
  185. c         print*,'iter=',iter,' dy01=',dy01,'dy02=',dy02
  186.  
  187.          ldom1=.false.
  188.          evol1=.false.
  189.          if (X(3).GT. Y01 .AND. X(5) .LT. ZMAX) then
  190.            if ( (1.+(A1*(-Y01 + X(3)))**B1-z1) .gt. 1.d-06) then
  191.               evol1=.true.
  192.               ldom1=.true.
  193.               x5=z1
  194.            endif
  195.          endif
  196.          ldom2=.false.
  197.          evol2=.false.
  198.          if (X(4).GT. Y02 .AND. X(6) .LT. ZMAX) then
  199.            if ( (1.+(A2*(-Y02 + X(4)))**B2-z2) .gt. 1.d-06) then
  200.               evol2=.true.
  201.               ldom2=.true.
  202.               x6=z2
  203.            endif
  204.          endif
  205. c         print*,'y1',x(3),'ylim1',ylim1
  206. c         print*,'z1',x(5),'x5',x5
  207. c         print*,'ldom1',ldom1,'ldom2',ldom2
  208. C        =====> PROBLEME ELASTIQUE
  209. C
  210. C
  211. C        CALCUL DU JACOBIEN
  212. C
  213. C
  214.          JACOB(3,1) = H1*X(1)*X(5)*X(5)/E
  215. C
  216.          JACOB(3,2) = H2*X(2)*X(5)*X(5)/E
  217. C
  218.          JACOB(4,1) = (UN-H1)*X(1)*X(6)*X(6)/E
  219. C
  220.          JACOB(4,2) = (UN-H2)*X(2)*X(6)*X(6)/E
  221. C
  222. C        CALCUL DU SECOND MEMBRE
  223. C
  224.       XIG1 =((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP2+(X(6)-UN)
  225.      1  *BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)
  226.      2  *HT2+BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT+
  227.      3  BETA2*(X(6)*(X(5)*H2+H2-UN)-X(5)*H2+X(6)**2*(UN-H2))+E*EPSP1*
  228.      4  (X(6)*(H2-UN)-X(5)*H2))*SIGF+(BETA1*E*(X(5)*EPSP2-EPSP2)+
  229.      5  (UN-X(5))*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-UN)-
  230.      6  X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF+
  231.      7  (DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*
  232.      8  (-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))
  233.      9  *BETA1*H1)*HT2)
  234. C
  235.       XIG2 =((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP1+(X(6)-UN)
  236.      1  *BETA2)*NU+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)*
  237.      2  X(6)*BETA1)*HT2+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+
  238.      3  (X(5)-UN)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2-X(5)*EPSP2)+(-X(6)**2+
  239.      4  (X(5)+UN)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPSP2+(X(6)**2-X(6))*
  240.      5  BETA2)*SIGF+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)-UN)*BETA1*E*EPSP1)
  241.      6  *HT2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*
  242.      7  (UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF+(DEUX-DEUX*X(5))*
  243.      8  BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*
  244.      9  (H2-DEUX)+DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)
  245.  
  246.          DELAFON(1) =X(1)-XIG1
  247. C
  248.          DELAFON(2) =X(2)-XIG2
  249. C
  250.          XX= (H1*X(1)*X(1)+H2*X(2)*X(2))*X(5)*X(5)/(DEUX*E)
  251.          DELAFON(3) = X(3)-XX
  252. C
  253.          XX=((UN-H1)*X(1)*X(1)+(UN-H2)*X(2)*X(2) )
  254.      1      *X(6)*X(6)/(DEUX*E)
  255.          DELAFON(4) = X(4)-XX
  256. C
  257.          IF (LDOM1) THEN
  258. C
  259. C        =====> ENDOMMAGEMENT EN TRACTION
  260. C
  261.       J15 = (((BETA1*HT2+BETA1*HT)*NU+BETA1*(DEUX*X(5)*H2-H2+X(6)*(UN-H2
  262.      &   ))*HT2+BETA1*(DEUX*X(5)*H2-H2+X(6)*(UN-H2))*HT+BETA2*(X(6)*H2-H
  263.      &   2)-E*EPSP1*H2)*SIGF+(BETA1*E*EPSP2-BETA1*E*EPSP1)*HT2)/((NU**2+
  264.      &   H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(6
  265.      &   )*H2+X(6)**2*(H2-UN))*SIGF+(DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA1
  266.      &   *(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-UN
  267.      &   )*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X(6)*(DEUX*H2-UN)-DEU
  268.      &   X*X(5)*H2)-X(6)*H2)*SIGF-DEUX*BETA1*HT2*NU+(BETA1*(-DEUX*X(5)*H
  269.      &   2+H2+X(6)*(H2-DEUX))+(X(6)-DEUX*X(5)+UN)*BETA1*H1)*HT2)*((((X(5
  270.      &   )-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP2+(X(6)-UN)*BETA2)*NU+
  271.      &   BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT2+BETA1*
  272.      &   (X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT+BETA2*(X(6)*(
  273.      &   X(5)*H2+H2-UN)-X(5)*H2+X(6)**2*(UN-H2))+E*EPSP1*(X(6)*(H2-UN)-X
  274.      &   (5)*H2))*SIGF+(BETA1*E*(X(5)*EPSP2-EPSP2)+(UN-X(5))*BETA1*E*EPS
  275.      &   P1)*HT2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*
  276.      &   (UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF+(DEUX-DEUX*X(5))*BE
  277.      &   TA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)
  278.      &   +DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  279. C
  280. C
  281.       J25 = (((BETA1*HT2+BETA1*HT)*NU+((-X(6)+DEUX*X(5)-UN)*BETA1*H1+X(6
  282.      &   )*BETA1)*HT2+((-X(6)+DEUX*X(5)-UN)*BETA1*H1+X(6)*BETA1)*HT+((X(
  283.      &   6)-UN)*BETA2-E*EPSP2)*H1)*SIGF+(BETA1*E*EPSP1-BETA1*E*EPSP2)*HT
  284.      &   2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2
  285.      &   ))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF+(DEUX-DEUX*X(5))*BETA1*HT
  286.      &   2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX)
  287.      &   )+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X(6)*(DEUX
  288.      &   *H2-UN)-DEUX*X(5)*H2)-X(6)*H2)*SIGF-DEUX*BETA1*HT2*NU+(BETA1*(-
  289.      &   DEUX*X(5)*H2+H2+X(6)*(H2-DEUX))+(X(6)-DEUX*X(5)+UN)*BETA1*H1)*H
  290.      &   T2)*((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP1+(X(6)-UN)
  291.      &   *BETA2)*NU+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)*X(
  292.      &   6)*BETA1)*HT2+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)
  293.      &   *X(6)*BETA1)*HT+(E*(X(6)*EPSP2-X(5)*EPSP2)+(-X(6)**2+(X(5)+UN)*
  294.      &   X(6)-X(5))*BETA2)*H1-X(6)*E*EPSP2+(X(6)**2-X(6))*BETA2)*SIGF+(B
  295.      &   ETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)-UN)*BETA1*E*EPSP1)*HT2)/((NU**2
  296.      &   +H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(
  297.      &   6)*H2+X(6)**2*(H2-UN))*SIGF+(DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA
  298.      &   1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-U
  299.      &   N)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  300. C
  301.             JACOB(1,5) = J15
  302.             JACOB(2,5) = J25
  303.             JACOB(3,5) = (H1*X(1)*X(1)+H2*X(2)*X(2))*X(5)/E
  304.             JACOB(4,5) = 0.D0
  305.             JACOB(5,5) = -1.D0
  306.      &     -A1*B1*(A1*(X(3)-Y01))**(B1-UN)
  307.      &     *alpha1*m1*dy01**alpha1/(dy01+hdp1)/(x(5)**2*dt)
  308.             JACOB(5,1) = 0.D0
  309.             JACOB(5,2) = 0.D0
  310.             JACOB(5,3) = A1*B1*(A1*(X(3)-Y01))**(B1-UN)
  311.             JACOB(5,4) = 0.D0
  312.             DELAFON(5) = X(5)-((A1*(X(3)-Y01))**B1+UN)
  313.  
  314.             IF (LDOM2) THEN
  315. C
  316. C        =====> ENDOMMAGEMENT EN COMPRESSION
  317. C
  318.       J16 = (BETA2*NU+BETA1*(H2+X(5)*(1-H2)-1)*HT2+BETA1*(H2+X(5)*(1-H2)
  319.      1   -1)*HT+BETA2*(X(5)*H2+H2+2*X(6)*(1-H2)-1)+E*EPSP1*(H2-1))*SIGF/
  320.      2   ((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6)**2*(1-H2))-X(5)*
  321.      3   X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-2*X(5))*BETA1*HT2*NU+(BETA1*(-X
  322.      4   (5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+((X(5)-1)*X(6)-X(5)
  323.      5   **2+X(5))*BETA1*H1)*HT2)-((H1*(X(5)*(2*H2-1)+2*X(6)*(1-H2))-X(5
  324.      6   )*H2+2*X(6)*(H2-1))*SIGF+(BETA1*(-H2+X(5)*(H2-2)+2)+(X(5)-1)*BE
  325.      7   TA1*H1)*HT2)*((((X(5)-1)*BETA1*HT2+(X(5)-1)*BETA1*HT-E*EPSP2+(X
  326.      8   (6)-1)*BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)
  327.      9   *H2)*HT2+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)*H2)*HT+
  328.      :   BETA2*(X(6)*(X(5)*H2+H2-1)-X(5)*H2+X(6)**2*(1-H2))+E*EPSP1*(X(6
  329.      ;   )*(H2-1)-X(5)*H2))*SIGF+(BETA1*E*(X(5)*EPSP2-EPSP2)+(1-X(5))*BE
  330.      <   TA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6
  331.      =   )**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-2*X(5))*BETA1
  332.      >   *HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+(
  333.      ?   (X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  334. C
  335.       J26 = (BETA2*NU+((1-X(5))*BETA1*H1+(X(5)-1)*BETA1)*HT2+((1-X(5))*B
  336.      1   ETA1*H1+(X(5)-1)*BETA1)*HT+(E*EPSP2+(-2*X(6)+X(5)+1)*BETA2)*H1-
  337.      2   E*EPSP2+(2*X(6)-1)*BETA2)*SIGF/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X
  338.      3   (5)**2*H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-
  339.      4   2*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)
  340.      5   *(H2-2)+2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X
  341.      6   (5)*(2*H2-1)+2*X(6)*(1-H2))-X(5)*H2+2*X(6)*(H2-1))*SIGF+(BETA1*
  342.      7   (-H2+X(5)*(H2-2)+2)+(X(5)-1)*BETA1*H1)*HT2)*((((X(5)-1)*BETA1*H
  343.      8   T2+(X(5)-1)*BETA1*HT-E*EPSP1+(X(6)-1)*BETA2)*NU+(((1-X(5))*X(6)
  344.      9   +X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT2+(((1-X(5))*X(6
  345.      :   )+X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2
  346.      ;   -X(5)*EPSP2)+(-X(6)**2+(X(5)+1)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPS
  347.      <   P2+(X(6)**2-X(6))*BETA2)*SIGF+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)
  348.      =   -1)*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*
  349.      >   H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-2*X(5))
  350.      ?   *BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)
  351.      @   +2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  352. C
  353. C
  354.                JACOB(1,6) = J16
  355.                JACOB(2,6) = J26
  356.                JACOB(3,6) = 0.D0
  357.                JACOB(4,6) = ((1.D0-H1)*X(1)*X(1)+(1.D0-H2)
  358.      &                      *X(2)*X(2))*X(6)/E
  359.                JACOB(5,6) = 0.D0
  360.                JACOB(6,6) = -1.D0
  361.      &     -A2*B2*(A2*(X(4)-Y02))**(B2-UN)
  362.      &     *alpha2*m2*dy02**alpha2/(dy02+hdp2)/(x(6)**2*dt)
  363.                JACOB(6,1) = 0.D0
  364.                JACOB(6,2) = 0.D0
  365.                JACOB(6,3) = 0.D0
  366. *               JACOB(6,4) = (A2*B2*(X(4)-Y02)**(B2-1.0))
  367.                JACOB(6,4) = A2*B2*(A2*(X(4)-Y02))**(B2-1.0)
  368.                JACOB(6,5) = 0.D0
  369. *               DELAFON(6) = X(6)-(A2*(X(4)-Y02)**B2+1.0)
  370.                DELAFON(6) = X(6)-((A2*(X(4)-Y02))**B2+1.0)
  371.                DIMPB=6
  372.             ELSE
  373.                DIMPB=5
  374.             END IF
  375.          ELSE
  376. C
  377. C        =====> PAS D'ENDOMMAGEMENT EN TRACTION
  378. C
  379.             IF (LDOM2) THEN
  380. C
  381. C        =====> ENDOMMAGEMENT EN COMPRESSION
  382. C
  383.       J16 = (BETA2*NU+BETA1*(H2+X(5)*(1-H2)-1)*HT2+BETA1*(H2+X(5)*(1-H2)
  384.      1   -1)*HT+BETA2*(X(5)*H2+H2+2*X(6)*(1-H2)-1)+E*EPSP1*(H2-1))*SIGF/
  385.      2   ((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6)**2*(1-H2))-X(5)*
  386.      3   X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-2*X(5))*BETA1*HT2*NU+(BETA1*(-X
  387.      4   (5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+((X(5)-1)*X(6)-X(5)
  388.      5   **2+X(5))*BETA1*H1)*HT2)-((H1*(X(5)*(2*H2-1)+2*X(6)*(1-H2))-X(5
  389.      6   )*H2+2*X(6)*(H2-1))*SIGF+(BETA1*(-H2+X(5)*(H2-2)+2)+(X(5)-1)*BE
  390.      7   TA1*H1)*HT2)*((((X(5)-1)*BETA1*HT2+(X(5)-1)*BETA1*HT-E*EPSP2+(X
  391.      8   (6)-1)*BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)
  392.      9   *H2)*HT2+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)*H2)*HT+
  393.      :   BETA2*(X(6)*(X(5)*H2+H2-1)-X(5)*H2+X(6)**2*(1-H2))+E*EPSP1*(X(6
  394.      ;   )*(H2-1)-X(5)*H2))*SIGF+(BETA1*E*(X(5)*EPSP2-EPSP2)+(1-X(5))*BE
  395.      <   TA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6
  396.      =   )**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-2*X(5))*BETA1
  397.      >   *HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+(
  398.      ?   (X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  399. C
  400.       J26 = (BETA2*NU+((1-X(5))*BETA1*H1+(X(5)-1)*BETA1)*HT2+((1-X(5))*B
  401.      1   ETA1*H1+(X(5)-1)*BETA1)*HT+(E*EPSP2+(-2*X(6)+X(5)+1)*BETA2)*H1-
  402.      2   E*EPSP2+(2*X(6)-1)*BETA2)*SIGF/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X
  403.      3   (5)**2*H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-
  404.      4   2*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)
  405.      5   *(H2-2)+2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X
  406.      6   (5)*(2*H2-1)+2*X(6)*(1-H2))-X(5)*H2+2*X(6)*(H2-1))*SIGF+(BETA1*
  407.      7   (-H2+X(5)*(H2-2)+2)+(X(5)-1)*BETA1*H1)*HT2)*((((X(5)-1)*BETA1*H
  408.      8   T2+(X(5)-1)*BETA1*HT-E*EPSP1+(X(6)-1)*BETA2)*NU+(((1-X(5))*X(6)
  409.      9   +X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT2+(((1-X(5))*X(6
  410.      :   )+X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2
  411.      ;   -X(5)*EPSP2)+(-X(6)**2+(X(5)+1)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPS
  412.      <   P2+(X(6)**2-X(6))*BETA2)*SIGF+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)
  413.      =   -1)*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*
  414.      >   H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-2*X(5))
  415.      ?   *BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)
  416.      @   +2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  417. C
  418. C
  419.  
  420.                JACOB(1,5) = J16
  421.                JACOB(2,5) = J26
  422.                JACOB(3,5) = 0.D0
  423.                JACOB(4,5) = ((1.D0-H1)*X(1)*X(1)+(1.D0-H2)
  424.      &                      *X(2)*X(2))*X(6)/E
  425.                JACOB(5,5) = -1.D0
  426.      &     -A2*B2*(A2*(X(4)-Y02))**(B2-UN)
  427.      &     *alpha2*m2*dy02**alpha2/(dy02+hdp2)/(x(6)**2*dt)
  428.                JACOB(5,1) = 0.D0
  429.                JACOB(5,2) = 0.D0
  430.                JACOB(5,3) = 0.D0
  431. *               JACOB(5,4) = (A2*B2*(X(4)-Y02)**(B2-1.0))
  432.                JACOB(5,4) = A2*B2*(A2*(X(4)-Y02))**(B2-1.0)
  433. *               DELAFON(5) = X(6)-(A2*(X(4)-Y02)**B2+1.0)
  434.                DELAFON(5) = X(6)-((A2*(X(4)-Y02))**B2+1.0)
  435.                DIMPB=5
  436.             ELSE
  437.                DIMPB=4
  438.             END IF
  439.          END IF
  440. C
  441. C        RESOLUTION
  442. C
  443.          CALL GAUSCL(JACOB,DX,DELAFON,NDIM,DIMPB,0)
  444. C
  445. C        CALCUL DE L"ERREUR ET DE LA NOUVELLE SOLUTION
  446. C
  447.          ERREUR=0.D0
  448.          DO I=1,2
  449.             ERREUR=ERREUR+DX(I)*DX(I)/(SIGY*SIGY)
  450.             X(I)=X(I)+DX(I)
  451.             if (x(i) .ge. 0.) then
  452.                x(i)=MIN(1.d+15,x(i))
  453.             else
  454.                x(i)=max(-1.d+15,x(i))
  455.             endif
  456.          END DO
  457.          X(3)=X(3)+DX(3)
  458.          X3=MAX(X(3),Y001)
  459.          X3=MIN(ABS(X3),1.d+07)
  460.          ERREUR=ERREUR+DX(3)*DX(3)/(x3*x3)
  461.          X(4)=X(4)+DX(4)
  462.          X4=MAX(X(4),Y002)
  463.          X4=MIN(ABS(X4),1.d+07)
  464.          ERREUR=ERREUR+DX(4)*DX(4)/(x4*x4)
  465.          if (.not. ldom1) then
  466.             dx5=(x(5)-x5)*1.0/4.0
  467.             x(5)=dx5+x5
  468.             x(5)=max(z1,x(5))
  469.             ERREUR=ERREUR+DX5*DX5/(z1*z1)
  470.          endif
  471.          if (.not. ldom2) then
  472.             dx6=(x(6)-x6)/4.
  473.             x(6)=dx6+x6
  474.             x(6)=max(z2,x(6))
  475.             ERREUR=ERREUR+DX6*DX6/(z2*z2)
  476.          endif
  477.          IF (LDOM1) THEN
  478.            IF ((X(5)+DX(5)).LT.Z1) THEN
  479.               DX(5)=(Z1-X(5))
  480.  
  481. c              if (evol1) dx(5)=dx(5)/2.
  482.  
  483.            END IF
  484.            IF ( X(5) .GT. ZMAX)THEN
  485.              DX(5)=ZMAX + 0.1D0 - X(5)
  486.            END IF
  487.            X(5)=X(5)+DX(5)
  488.            ERREUR=ERREUR+DX(5)*DX(5)/(z1*z1)
  489.            IF (LDOM2) THEN
  490.               IF ((X(6)+DX(6)).LT.Z2) THEN
  491.                  DX(6)=(Z2-X(6))
  492.  
  493. c              if (evol2) dx(6)=dx(6)/2.
  494.  
  495.               END IF
  496.               IF ( X(6) .GT. ZMAX)THEN
  497.                 DX(6)=ZMAX + 0.1D0 - X(6)
  498.               END IF
  499.               X(6)=X(6)+DX(6)
  500.               ERREUR=ERREUR+DX(6)*DX(6)/(z2*z2)
  501.            END IF
  502.          ELSE IF (LDOM2) THEN
  503.            IF ((X(6)+DX(5)).LT.Z2) THEN
  504.               DX(5)=(Z2-X(6))
  505.  
  506. c              if (evol2) dx(5)=dx(5)/2.
  507.  
  508.            END IF
  509.            IF ( X(6) .GT. ZMAX)THEN
  510.              DX(5)=ZMAX + 0.1D0 - X(6)
  511.            END IF
  512.            X(6)=X(6)+DX(5)
  513.            ERREUR=ERREUR+DX(5)*DX(5)/(z2*z2)
  514.          END IF
  515.  
  516.  
  517.          IF (ERREUR .GT. 1.D+100)THEN
  518.             CONV=.FALSE.
  519.            print*,'erreur dans modele',erreur
  520.            print*,'iteration num',iter
  521.            print*,'delafon',delafon
  522.            print*,z1,z2,x3,x4,sigy
  523.            print*,dx
  524.            print*,'y02=',y02,'y002=',y002
  525.             RETURN
  526.          END IF
  527. C
  528. C        QUELQUES IMPRESSIONS EN CAS D'ERREUR
  529. C
  530.          if (.not. conv )then
  531.             PRINT*,'ITERATION  ',ITER,'ERREUR',erreur
  532.             print*,'LDOM1  :',LDOM1,'   LDOM2  :',LDOM2
  533.             PRINT*,'EPS1= ',EPSP1,'  EPS2= ',EPSP2
  534.             PRINT*,'SIG1:',X(1) ,'   DSIG1:',DX(1),'ERR :',DELAFON(1)
  535.             PRINT*,'SIG2:',X(2) ,'   DSIG2:',DX(2),'ERR :',DELAFON(2)
  536.             PRINT*,'Y1  :',X(3) ,'   DY1  :',DX(3),'ERR :',DELAFON(3)
  537.             PRINT*,'Y2  :',X(4) ,'   DY2  :',DX(4),'ERR :',DELAFON(4)
  538.             PRINT*,'Z1  :',X(5) ,'   DZ1  :',DX(5),'ERR :',DELAFON(5)
  539.             PRINT*,'Z2  :',X(6) ,'   DZ2  :',DX(6),'ERR :',DELAFON(6)
  540.             PRINT*,'H1',H1,'H2',H2,ht,ht2
  541.             PRINT*
  542.             write(*,200)jacob
  543.          END IF
  544.  
  545.       END DO
  546.       IF (ERREUR .GT. TEST ) THEN
  547.          CONV=.FALSE.
  548.       ELSE
  549.          CONV=.TRUE.
  550.          SIG1= X(1)
  551.          SIG2= X(2)
  552.          Y1    = X(3)
  553.          Y2    = X(4)
  554.          Z1    = X(5)
  555.          Z2    = X(6)
  556.          IF (SIG1 .GE. 0.D0) THEN
  557.             JSIGN(1)=1
  558.          ELSE
  559.             JSIGN(1)=0
  560.          END IF
  561.          IF (SIG2.GE. 0.D0)THEN
  562.             JSIGN(2)=1
  563.          ELSE
  564.             JSIGN(2)=0
  565.          END IF
  566.          JSIGN(3)=1
  567.       END IF
  568.       RETURN
  569. 200   format(6(6(2x,e11.4),/))
  570.       END
  571.  
  572.  
  573.  
  574.  
  575.  
  576.  
  577.  

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