Télécharger la plaquette Cast3M


fr en

Télécharger codjf2.eso

Retour à la liste

Numérotation des lignes : 

codjf2

  1. C CODJF2    SOURCE    PV        17/12/08    21:16:15     9660           
  2. C MODJF2    SOURCE    AM        00/12/13    21:41:10     4045
  3.       SUBROUTINE CODJF2(WRK52,WRK53,EPSP1,EPSP2,YLIM1,YLIM2,Z1,Z2,Y1,Y2,
  4.      1                  JSIGN,SIG1,SIG2,CONV,DT1)
  5. C
  6. C     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  7. C     %                                                               %
  8. C     %  INTEGRATION DU MODJFD2 UNILATERAL AVEC FERMETURE LINEAIRE     %
  9. C     %          DE FISSURES EN CONTRAINTES PLANES                    %
  10. C     %                (MODCLB NON ASSOCIE)                           %
  11. C     %      A L"AIDE DE LA METHODE DE DICHOTOMIE                     %
  12. C     %    AJOUT VICOENDOMMAGEMENT J.F. DUBE                          %
  13. C     %                                                               %
  14. C     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  15. C
  16. C
  17. C
  18.       IMPLICIT INTEGER(I-N)
  19. -INC DECHE
  20. C        VARIABLES D"ENTREES
  21. C
  22. C
  23.       REAL*8 EPSP1,EPSP2
  24. C
  25. C        VARIABLES DE SORTIE
  26. C
  27. C
  28.       REAL*8 SIG1,SIG2,Z1,Z2,Y1,Y2,DT1
  29.       INTEGER JSIGN(3)
  30.       LOGICAL CONV
  31. C
  32. C
  33. C        PARAMETRES D"EVOLUTION
  34. C
  35. C
  36.       REAL*8 E,NU,A1,A2,B1,B2,Y01,Y02,SIGF1
  37.       REAL*8 BETA1,BETA2,YLIM1,YLIM2
  38. C
  39. C
  40. C        VARIABLES INTERNES
  41. C
  42. C
  43.       REAL*8 Y001,Y002,DY01,DY02,hdp1,hdp2
  44.       REAL*8 DELTASIG
  45.       INTEGER NDIM
  46.       PARAMETER (NDIM=6)
  47.       REAL*8 JACOB(NDIM,NDIM),J15,J16,J25,J26
  48.       REAL*8 DELAFON(NDIM)
  49.       REAL*8 X(NDIM),DX(NDIM),X0(NDIM),XX,XIG1,XIG2
  50.       REAL*8 ERREUR,TEST
  51.       REAL*8 H1,H2,HT,HT2,trf
  52.       REAL*8 ALFA1,ALFA2
  53.       REAL*8 SIGY
  54.       REAL*8 UN,DEUX
  55.       REAL*8 ZMAX
  56.       INTEGER ITER ,I,DIMPB
  57.       LOGICAL LDOM1,LDOM2,evol1,evol2
  58.       REAL*8 M1,M2,ALPHA1,ALPHA2
  59.       REAL*8 x3,x4,x5,x6,dx5,dx6
  60.       INTEGER IMAX
  61.       REAL*8 pas1,erreurd1,xd1,xg1,pas2,erreurd2,xd2,xg2
  62.       REAL*8 erreur_max,fctd1,fctd2
  63.       logical passage1,passage2
  64.       UN=1.D0
  65.       DEUX=2.D0
  66.       ZMAX=1.d+04
  67.       E     = XMAT(1)
  68.       NU    = XMAT(2)
  69.       Y001  = XMAT(5)
  70.       Y002  = XMAT(6)
  71.       A1    = XMAT(7)
  72.       A2    = XMAT(8)
  73.       B1    = XMAT(9)
  74.       B2    = XMAT(10)
  75.       BETA1 = XMAT(11)
  76.       BETA2 = XMAT(12)
  77.       SIGF1  = XMAT(13)
  78.       M1    = XMAT(14)
  79.       M2    = XMAT(15)
  80.       ALPHA1= XMAT(16)
  81.       ALPHA2= XMAT(17)
  82.       SIGY  =SQRT(2.D0*E*Y001)
  83. *      DELTASIG=SIGY*1.D-4
  84.       DELTASIG=0.d0
  85. *      print*,'E,NU,Y01,Y02,A1,A2,B1,B2,BETA1,BETA2'
  86. *      print*,E,NU,Y01,Y02,A1,A2,B1,B2,BETA1,BETA2
  87.       TEST=1.D-12
  88.       ERREUR=1.D0
  89.       ITER=0
  90.       X(1)=SIG1
  91.       X(2)=SIG2
  92.       X(3)=Y1
  93.       X(4)=Y2
  94.       X(5)=Z1
  95.       X(6)=Z2
  96.       DO I=1, NDIM
  97.          X0(I)=X(I)
  98.       END DO
  99.       x5=Z1
  100.       x6=z2
  101. C
  102. C
  103. C
  104.       JACOB(1,1) = -1.d0
  105.       JACOB(1,2) = 0.D0
  106.       JACOB(1,3) = 0.D0
  107.       JACOB(1,4) = 0.D0
  108.       JACOB(1,5) = 0.D0
  109.       JACOB(1,6) = 0.D0
  110.       JACOB(2,1) = 0.D0
  111.       JACOB(2,2) = -1.d0
  112.       JACOB(2,3) = 0.D0
  113.       JACOB(2,4) = 0.D0
  114.       JACOB(2,5) = 0.D0
  115.       JACOB(2,6) = 0.D0
  116.       JACOB(3,1) = 0.D0
  117.       JACOB(3,2) = 0.D0
  118.       JACOB(3,3) = -1.d0
  119.       JACOB(3,4) = 0.D0
  120.       JACOB(3,5) = 0.D0
  121.       JACOB(3,6) = 0.D0
  122.       JACOB(4,1) = 0.D0
  123.       JACOB(4,2) = 0.D0
  124.       JACOB(4,3) = 0.D0
  125.       JACOB(4,4) = -1.d0
  126.       JACOB(4,5) = 0.D0
  127.       JACOB(4,6) = 0.D0
  128.       JACOB(5,1) = 0.D0
  129.       JACOB(5,2) = 0.D0
  130.       JACOB(5,3) = 0.D0
  131.       JACOB(5,4) = 0.D0
  132.       JACOB(5,5) = -1.d0
  133.       JACOB(5,6) = 0.D0
  134.       JACOB(6,1) = 0.D0
  135.       JACOB(6,2) = 0.D0
  136.       JACOB(6,3) = 0.D0
  137.       JACOB(6,4) = 0.D0
  138.       JACOB(6,5) = 0.D0
  139.       JACOB(6,6) = -1.d0
  140.       Erreur_max=1.0d-10
  141.       pas1=1.0
  142.       Xd1=z1
  143.       Xg1=z1
  144.       imax=0
  145.       Erreurd1=1.0
  146.       passage1=.true.
  147.       Do while ( (Erreurd1 .GT. Erreur_max) .and. (pas1 .ge. 1.d-6) )
  148.          x(5)=xd1+pas1
  149.  
  150. c         print*,'*** Z1 ***'
  151. c         print*,x(5)
  152.  
  153. **      Calcul de Y01 modifie par l'effet de vitesse
  154. c        dp1=(DIM(x(5),z1)/(x(5)*z1*dt1))
  155. c         y01=y001+m1*dp1**alpha1
  156.  
  157.          pas2=1.0
  158.          Xd2=z2
  159.          Xg2=x(6)
  160.          Erreurd2=1.0
  161.          passage2=.true.
  162.          X(4)=Y2
  163.          Do while ((Erreurd2 .GT. Erreur_max) .and. (pas2 .ge. 1.d-6))
  165.             if (imax .ge. 10000) then
  166.               print*,'pb dans jfd2'
  167.               conv=.false.
  168.               return
  169.             end if
  170.             x(6)=xd2+pas2
  171.  
  172. c            print*,'*** Z2 ***'
  173. c            print*,x(6)
  174.  
  175. **      Calcul de Y02 modifie par l'effet de vitesse
  176. c            dp2=(DIM(x(6),z2)/(x(6)*z2*dt1))
  177. c            y02=y002+m2*dp2**alpha2
  178. **      Calcul de sigma et Y1 et Y2
  179.             ERREUR=1.D0
  180.             ITER=0
  181.             X(1)=Sig1
  182.             X(2)=Sig2
  183.             DO WHILE (ERREUR .GT. TEST .AND. ITER .LE. 100)
  184.                ITER=ITER+1
  185.                IF (X(1) .GE. DELTASIG) THEN
  186.                   H1=UN
  187.                ELSE IF (X(1) .LE. -DELTASIG) THEN
  188.                   H1=0.D0
  189.                ELSE
  190.                   H1=(X(1)+DELTASIG)/(DEUX*DELTASIG)
  191.                END IF
  192.                IF (X(2) .GE. DELTASIG) THEN
  193.                   H2=UN
  194.                ELSE IF (X(2) .LE. -DELTASIG) THEN
  195.                   H2=0.D0
  196.                ELSE
  197.                   H2=(X(2)+DELTASIG)/(DEUX*DELTASIG)
  198.                END IF
  199. c
  200. c     smooth transition on ht ht2
  201. c
  202.                trf = x(1) + x(2)
  203.                IF(trf.GE.deltasig)THEN
  204.                   HT=UN
  205.                   HT2=0.D0
  206.                else if(trf.ge.(-deltasig)) then
  207.                   ht = (trf+deltasig)/DEUX/deltasig
  208.                   ht2 = UN-ht
  209.                else if(trf.ge.(-sigf+deltasig)) then
  210.                   ht = 0.D+00
  211.                   ht2 = UN
  212.                else if(trf.ge.(-sigf-deltasig)) then
  213.                   ht = 0.D+00
  214. c           ht2 =UN- (trf-(sigf-deltasig))/(DEUX*sigf)
  215.                   ht2 = (trf+sigf+deltasig)/DEUX/deltasig
  216.                else
  217.                   ht = 0.D+00
  218.                   ht2 = 0.D+00
  219.                endif
  220. c         print*,'DTjfd',dt1
  221. **      Calcul de Y01 modifie par l'effet de vitesse
  222.                dy01=(DIM(x(5),z1)/(x(5)*z1*dt1))
  223.                hdp1=1.0
  224.                if (dy01 .gt. 0.0) hdp1=0.0
  225.                y01=y001+m1*dy01**alpha1
  226. **      Calcul de Y02 modifie par l'effet de vitesse
  227.                dy02=(DIM(x(6),z2)/(x(6)*z2*dt1))
  228.                hdp2=1.0
  229.                if (dy02 .gt. 0.0) hdp2=0.0
  230.                y02=y002+m2*dy02**alpha2
  231. c         print*,'dans clbeco iter=',iter,' y01=',y01,'y02=',y02
  232.  
  233.                ldom1=.false.
  234.                evol1=.false.
  235.                ldom2=.false.
  236.                evol2=.false.
  237. c         print*,'y1',x(3),'ylim1',ylim1
  238. c         print*,'z1',x(5),'x5',x5
  239. c         print*,'ldom1',ldom1,'ldom2',ldom2
  240. C        =====> PROBLEME ELASTIQUE
  241. C
  242. C
  243. C        CALCUL DU JACOBIEN
  244. C
  245. C
  246.          JACOB(3,1) = H1*X(1)*X(5)*X(5)/E
  247. C
  248.          JACOB(3,2) = H2*X(2)*X(5)*X(5)/E
  249. C
  250.          JACOB(4,1) = (UN-H1)*X(1)*X(6)*X(6)/E
  251. C
  252.          JACOB(4,2) = (UN-H2)*X(2)*X(6)*X(6)/E
  253. C
  254. C        CALCUL DU SECOND MEMBRE
  255. C
  256.       XIG1 =((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP2+(X(6)-UN)
  257.      1  *BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)
  258.      2  *HT2+BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT+
  259.      3  BETA2*(X(6)*(X(5)*H2+H2-UN)-X(5)*H2+X(6)**2*(UN-H2))+E*EPSP1*
  260.      4  (X(6)*(H2-UN)-X(5)*H2))*SIGF1+(BETA1*E*(X(5)*EPSP2-EPSP2)+
  261.      5  (UN-X(5))*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-UN)-
  262.      6 X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF1+
  263.      7  (DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*
  264.      8  (-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))
  265.      9  *BETA1*H1)*HT2)
  266. C
  267.       XIG2 =((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP1+(X(6)-UN)
  268.      1  *BETA2)*NU+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)*
  269.      2  X(6)*BETA1)*HT2+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+
  270.      3  (X(5)-UN)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2-X(5)*EPSP2)+(-X(6)**2+
  271.      4  (X(5)+UN)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPSP2+(X(6)**2-X(6))*
  272.      5 BETA2)*SIGF1+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)-UN)*BETA1*E*EPSP1)
  273.      6  *HT2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*
  274.      7  (UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*
  275.      8  BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*
  276.      9  (H2-DEUX)+DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)
  277.  
  278.          DELAFON(1) =X(1)-XIG1
  279. C
  280.          DELAFON(2) =X(2)-XIG2
  281. C
  282.          XX= (H1*X(1)*X(1)+H2*X(2)*X(2))*X(5)*X(5)/(DEUX*E)
  283.          DELAFON(3) = X(3)-XX
  284. C
  285.          XX=((UN-H1)*X(1)*X(1)+(UN-H2)*X(2)*X(2) )
  286.      1      *X(6)*X(6)/(DEUX*E)
  287.          DELAFON(4) = X(4)-XX
  288. C
  289.          IF (LDOM1) THEN
  290. C
  291. C        =====> ENDOMMAGEMENT EN TRACTION
  292. C
  293.       J15 = (((BETA1*HT2+BETA1*HT)*NU+BETA1*(DEUX*X(5)*H2-H2+X(6)*(UN-H2
  294.      &   ))*HT2+BETA1*(DEUX*X(5)*H2-H2+X(6)*(UN-H2))*HT+BETA2*(X(6)*H2-H
  295.      &  2)-E*EPSP1*H2)*SIGF1+(BETA1*E*EPSP2-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(6
  297.      &  )*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA1
  298.      &   *(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-UN
  299.      &   )*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X(6)*(DEUX*H2-UN)-DEU
  300.      &  X*X(5)*H2)-X(6)*H2)*SIGF1-DEUX*BETA1*HT2*NU+(BETA1*(-DEUX*X(5)*H
  301.      &   2+H2+X(6)*(H2-DEUX))+(X(6)-DEUX*X(5)+UN)*BETA1*H1)*HT2)*((((X(5
  302.      &   )-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP2+(X(6)-UN)*BETA2)*NU+
  303.      &   BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT2+BETA1*
  304.      &   (X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT+BETA2*(X(6)*(
  305.      &   X(5)*H2+H2-UN)-X(5)*H2+X(6)**2*(UN-H2))+E*EPSP1*(X(6)*(H2-UN)-X
  306.      &  (5)*H2))*SIGF1+(BETA1*E*(X(5)*EPSP2-EPSP2)+(UN-X(5))*BETA1*E*EPS
  307.      &   P1)*HT2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*
  308.      &  (UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*BE
  309.      &   TA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)
  310.      &   +DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  311. C
  312. C
  313.       J25 = (((BETA1*HT2+BETA1*HT)*NU+((-X(6)+DEUX*X(5)-UN)*BETA1*H1+X(6
  314.      &   )*BETA1)*HT2+((-X(6)+DEUX*X(5)-UN)*BETA1*H1+X(6)*BETA1)*HT+((X(
  315.      &  6)-UN)*BETA2-E*EPSP2)*H1)*SIGF1+(BETA1*E*EPSP1-BETA1*E*EPSP2)*HT
  316.      &   2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2
  317.      &  ))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*BETA1*HT
  318.      &   2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX)
  319.      &   )+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X(6)*(DEUX
  320.      &  *H2-UN)-DEUX*X(5)*H2)-X(6)*H2)*SIGF1-DEUX*BETA1*HT2*NU+(BETA1*(-
  321.      &   DEUX*X(5)*H2+H2+X(6)*(H2-DEUX))+(X(6)-DEUX*X(5)+UN)*BETA1*H1)*H
  322.      &   T2)*((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP1+(X(6)-UN)
  323.      &   *BETA2)*NU+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)*X(
  324.      &   6)*BETA1)*HT2+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)
  325.      &   *X(6)*BETA1)*HT+(E*(X(6)*EPSP2-X(5)*EPSP2)+(-X(6)**2+(X(5)+UN)*
  326.      &  X(6)-X(5))*BETA2)*H1-X(6)*E*EPSP2+(X(6)**2-X(6))*BETA2)*SIGF1+(B
  327.      &   ETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)-UN)*BETA1*E*EPSP1)*HT2)/((NU**2
  328.      &   +H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(
  329.      &  6)*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA
  330.      &   1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-U
  331.      &   N)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  332. C
  333.             JACOB(1,5) = J15
  334.             JACOB(2,5) = J25
  335.             JACOB(3,5) = (H1*X(1)*X(1)+H2*X(2)*X(2))*X(5)/E
  336.             JACOB(4,5) = 0.D0
  337.             JACOB(5,5) = -1.D0
  338.      &     -A1*B1*(A1*(X(3)-Y01))**(B1-UN)
  339.      &     *alpha1*m1*dy01**alpha1/(dy01+hdp1)/(x(5)**2*dt1)
  340.             JACOB(5,1) = 0.D0
  341.             JACOB(5,2) = 0.D0
  342.             JACOB(5,3) = A1*B1*(A1*(X(3)-Y01))**(B1-UN)
  343.             JACOB(5,4) = 0.D0
  344.             DELAFON(5) = X(5)-((A1*(X(3)-Y01))**B1+UN)
  345.  
  346.             IF (LDOM2) THEN
  347. C
  348. C        =====> ENDOMMAGEMENT EN COMPRESSION
  349. C
  350.       J16 = (BETA2*NU+BETA1*(H2+X(5)*(1-H2)-1)*HT2+BETA1*(H2+X(5)*(1-H2)
  351.      1  -1)*HT+BETA2*(X(5)*H2+H2+2*X(6)*(1-H2)-1)+E*EPSP1*(H2-1))*SIGF1/
  352.      2   ((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6)**2*(1-H2))-X(5)*
  353.      3  X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))*BETA1*HT2*NU+(BETA1*(-X
  354.      4   (5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+((X(5)-1)*X(6)-X(5)
  355.      5   **2+X(5))*BETA1*H1)*HT2)-((H1*(X(5)*(2*H2-1)+2*X(6)*(1-H2))-X(5
  356.      6  )*H2+2*X(6)*(H2-1))*SIGF1+(BETA1*(-H2+X(5)*(H2-2)+2)+(X(5)-1)*BE
  357.      7   TA1*H1)*HT2)*((((X(5)-1)*BETA1*HT2+(X(5)-1)*BETA1*HT-E*EPSP2+(X
  358.      8   (6)-1)*BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)
  359.      9   *H2)*HT2+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)*H2)*HT+
  360.      :   BETA2*(X(6)*(X(5)*H2+H2-1)-X(5)*H2+X(6)**2*(1-H2))+E*EPSP1*(X(6
  361.      ;  )*(H2-1)-X(5)*H2))*SIGF1+(BETA1*E*(X(5)*EPSP2-EPSP2)+(1-X(5))*BE
  362.      <   TA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6
  363.      =  )**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))*BETA1
  364.      >   *HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+(
  365.      ?   (X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  366. C
  367.       J26 = (BETA2*NU+((1-X(5))*BETA1*H1+(X(5)-1)*BETA1)*HT2+((1-X(5))*B
  368.      1   ETA1*H1+(X(5)-1)*BETA1)*HT+(E*EPSP2+(-2*X(6)+X(5)+1)*BETA2)*H1-
  369.      2  E*EPSP2+(2*X(6)-1)*BETA2)*SIGF1/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X
  370.      3  (5)**2*H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-
  371.      4   2*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)
  372.      5   *(H2-2)+2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X
  373.      6  (5)*(2*H2-1)+2*X(6)*(1-H2))-X(5)*H2+2*X(6)*(H2-1))*SIGF1+(BETA1*
  374.      7   (-H2+X(5)*(H2-2)+2)+(X(5)-1)*BETA1*H1)*HT2)*((((X(5)-1)*BETA1*H
  375.      8   T2+(X(5)-1)*BETA1*HT-E*EPSP1+(X(6)-1)*BETA2)*NU+(((1-X(5))*X(6)
  376.      9   +X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT2+(((1-X(5))*X(6
  377.      :   )+X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2
  378.      ;   -X(5)*EPSP2)+(-X(6)**2+(X(5)+1)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPS
  379.      <  P2+(X(6)**2-X(6))*BETA2)*SIGF1+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)
  380.      =   -1)*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*
  381.      >  H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))
  382.      ?   *BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)
  383.      @   +2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  384. C
  385. C
  386.                JACOB(1,6) = J16
  387.                JACOB(2,6) = J26
  388.                JACOB(3,6) = 0.D0
  389.                JACOB(4,6) = ((1.D0-H1)*X(1)*X(1)+(1.D0-H2)
  390.      &                      *X(2)*X(2))*X(6)/E
  391.                JACOB(5,6) = 0.D0
  392.                JACOB(6,6) = -1.D0
  393.      &     -A2*B2*(A2*(X(4)-Y02))**(B2-UN)
  394.      &     *alpha2*m2*dy02**alpha2/(dy02+hdp2)/(x(6)**2*dt1)
  395.                JACOB(6,1) = 0.D0
  396.                JACOB(6,2) = 0.D0
  397.                JACOB(6,3) = 0.D0
  398. *               JACOB(6,4) = (A2*B2*(X(4)-Y02)**(B2-1.0))
  399.                JACOB(6,4) = A2*B2*(A2*(X(4)-Y02))**(B2-1.0)
  400.                JACOB(6,5) = 0.D0
  401. *               DELAFON(6) = X(6)-(A2*(X(4)-Y02)**B2+1.0)
  402.                DELAFON(6) = X(6)-((A2*(X(4)-Y02))**B2+1.0)
  403.                DIMPB=6
  404.             ELSE
  405.                DIMPB=5
  406.             END IF
  407.          ELSE
  408. C
  409. C        =====> PAS D'ENDOMMAGEMENT EN TRACTION
  410. C
  411.             IF (LDOM2) THEN
  412. C
  413. C        =====> ENDOMMAGEMENT EN COMPRESSION
  414. C
  415.       J16 = (BETA2*NU+BETA1*(H2+X(5)*(1-H2)-1)*HT2+BETA1*(H2+X(5)*(1-H2)
  416.      1  -1)*HT+BETA2*(X(5)*H2+H2+2*X(6)*(1-H2)-1)+E*EPSP1*(H2-1))*SIGF1/
  417.      2   ((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6)**2*(1-H2))-X(5)*
  418.      3  X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))*BETA1*HT2*NU+(BETA1*(-X
  419.      4   (5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+((X(5)-1)*X(6)-X(5)
  420.      5   **2+X(5))*BETA1*H1)*HT2)-((H1*(X(5)*(2*H2-1)+2*X(6)*(1-H2))-X(5
  421.      6  )*H2+2*X(6)*(H2-1))*SIGF1+(BETA1*(-H2+X(5)*(H2-2)+2)+(X(5)-1)*BE
  422.      7   TA1*H1)*HT2)*((((X(5)-1)*BETA1*HT2+(X(5)-1)*BETA1*HT-E*EPSP2+(X
  423.      8   (6)-1)*BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)
  424.      9   *H2)*HT2+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)*H2)*HT+
  425.      :   BETA2*(X(6)*(X(5)*H2+H2-1)-X(5)*H2+X(6)**2*(1-H2))+E*EPSP1*(X(6
  426.      ;  )*(H2-1)-X(5)*H2))*SIGF1+(BETA1*E*(X(5)*EPSP2-EPSP2)+(1-X(5))*BE
  427.      <   TA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6
  428.      =  )**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))*BETA1
  429.      >   *HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+(
  430.      ?   (X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  431. C
  432.       J26 = (BETA2*NU+((1-X(5))*BETA1*H1+(X(5)-1)*BETA1)*HT2+((1-X(5))*B
  433.      1   ETA1*H1+(X(5)-1)*BETA1)*HT+(E*EPSP2+(-2*X(6)+X(5)+1)*BETA2)*H1-
  434.      2  E*EPSP2+(2*X(6)-1)*BETA2)*SIGF1/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X
  435.      3  (5)**2*H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-
  436.      4   2*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)
  437.      5   *(H2-2)+2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X
  438.      6  (5)*(2*H2-1)+2*X(6)*(1-H2))-X(5)*H2+2*X(6)*(H2-1))*SIGF1+(BETA1*
  439.      7   (-H2+X(5)*(H2-2)+2)+(X(5)-1)*BETA1*H1)*HT2)*((((X(5)-1)*BETA1*H
  440.      8   T2+(X(5)-1)*BETA1*HT-E*EPSP1+(X(6)-1)*BETA2)*NU+(((1-X(5))*X(6)
  441.      9   +X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT2+(((1-X(5))*X(6
  442.      :   )+X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2
  443.      ;   -X(5)*EPSP2)+(-X(6)**2+(X(5)+1)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPS
  444.      <  P2+(X(6)**2-X(6))*BETA2)*SIGF1+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)
  445.      =   -1)*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*
  446.      >  H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))
  447.      ?   *BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)
  448.      @   +2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  449. C
  450. C
  451.  
  452.                JACOB(1,5) = J16
  453.                JACOB(2,5) = J26
  454.                JACOB(3,5) = 0.D0
  455.                JACOB(4,5) = ((1.D0-H1)*X(1)*X(1)+(1.D0-H2)
  456.      &                      *X(2)*X(2))*X(6)/E
  457.                JACOB(5,5) = -1.D0
  458.      &     -A2*B2*(A2*(X(4)-Y02))**(B2-UN)
  459.      &     *alpha2*m2*dy02**alpha2/(dy02+hdp2)/(x(6)**2*dt1)
  460.                JACOB(5,1) = 0.D0
  461.                JACOB(5,2) = 0.D0
  462.                JACOB(5,3) = 0.D0
  463. *               JACOB(5,4) = (A2*B2*(X(4)-Y02)**(B2-1.0))
  464.                JACOB(5,4) = A2*B2*(A2*(X(4)-Y02))**(B2-1.0)
  465. *               DELAFON(5) = X(6)-(A2*(X(4)-Y02)**B2+1.0)
  466.                DELAFON(5) = X(6)-((A2*(X(4)-Y02))**B2+1.0)
  467.                DIMPB=5
  468.             ELSE
  469.                DIMPB=4
  470.             END IF
  471.          END IF
  472. C
  473. C        RESOLUTION
  474. C
  475.          CALL GAUSCL(JACOB,DX,DELAFON,NDIM,DIMPB,0)
  476. C
  477. C        CALCUL DE L"ERREUR ET DE LA NOUVELLE SOLUTION
  478. C
  479.                ERREUR=0.D0
  480.                do i=1,ndim
  481.                   erreur=erreur+(delafon(i)*delafon(i))
  482.                enddo
  483.                ERREUR=0.D0
  484.                DO I=1,2
  485.                   ERREUR=ERREUR+DX(I)*DX(I)/(SIGY*SIGY)
  486.                   X(I)=X(I)+DX(I)
  487. c            if (x(i) .ge. 0.) then
  488. c               x(i)=MIN(1.d+15,x(i))
  489. c            else
  490. c               x(i)= MAX(-1.d+15,x(i))
  491. c            endif
  492.                END DO
  493.                X(3)=X(3)+DX(3)
  494.                X3= MAX(X(3),Y01)
  495. c         X3=MIN(ABS(X3),1.d+07)
  496.                ERREUR=ERREUR+DX(3)*DX(3)/(x3*x3)
  497.                X(4)=X(4)+DX(4)
  498.                X4= MAX(X(4),Y02)
  499. c         X4=MIN(ABS(X4),1.d+07)
  500.                ERREUR=ERREUR+DX(4)*DX(4)/(x4*x4)
  501. c         if (.not. ldom1) then
  502. c            dx5=(x(5)-x5)*1.0/4.0
  503. c            x(5)=dx5+x5
  504. c            x(5)= MAX(z1,x(5))
  505. c            ERREUR=ERREUR+DX5*DX5/(z1*z1)
  506. c         endif
  507. c         if (.not. ldom2) then
  508. c            dx6=(x(6)-x6)/4.
  509. c            x(6)=dx6+x6
  510. c            x(6)= MAX(z2,x(6))
  511. c            ERREUR=ERREUR+DX6*DX6/(z2*z2)
  512. c         endif
  513.                IF (LDOM1) THEN
  514.                  IF ((X(5)+DX(5)).LT.Z1) THEN
  515.                     DX(5)=(Z1-X(5))
  516.                  END IF
  517.                  IF ( X(5) .GT. ZMAX)THEN
  518.                    DX(5)=ZMAX + 0.1D0 - X(5)
  519.                  END IF
  520.                  X(5)=X(5)+DX(5)
  521.                  ERREUR=ERREUR+DX(5)*DX(5)/(z1*z1)
  522.                  IF (LDOM2) THEN
  523.                     IF ((X(6)+DX(6)).LT.Z2) THEN
  524.                        DX(6)=(Z2-X(6))
  525.                     END IF
  526.                     IF ( X(6) .GT. ZMAX)THEN
  527.                       DX(6)=ZMAX + 0.1D0 - X(6)
  528.                     END IF
  529.                     X(6)=X(6)+DX(6)
  530.                     ERREUR=ERREUR+DX(6)*DX(6)/(z2*z2)
  531.                  END IF
  532.                ELSE IF (LDOM2) THEN
  533.                  IF ((X(6)+DX(5)).LT.Z2) THEN
  534.                     DX(5)=(Z2-X(6))
  535.                  END IF
  536.                  IF ( X(6) .GT. ZMAX)THEN
  537.                    DX(5)=ZMAX + 0.1D0 - X(6)
  538.                  END IF
  539.                  X(6)=X(6)+DX(5)
  540.                  ERREUR=ERREUR+DX(5)*DX(5)/(z2*z2)
  541.                END IF
  542.  
  543.  
  544.          IF (ERREUR .GT. 1.D+100)THEN
  545.             CONV=.FALSE.
  546.            print*,'erreur dans modele',erreur
  547.            print*,'iteration num',iter
  548.            print*,'delafon',delafon
  549.             RETURN
  550.          END IF
  551. C
  552. C        QUELQUES IMPRESSIONS EN CAS D'ERREUR
  553. C
  554.          if (.not. conv )then
  555.             PRINT*,'ITERATION  ',ITER,'ERREUR',erreur
  556.             print*,'LDOM1  :',LDOM1,'   LDOM2  :',LDOM2
  557.             PRINT*,'EPS1= ',EPSP1,'  EPS2= ',EPSP2
  558.             PRINT*,'SIG1:',X(1) ,'   DSIG1:',DX(1),'ERR :',DELAFON(1)
  559.             PRINT*,'SIG2:',X(2) ,'   DSIG2:',DX(2),'ERR :',DELAFON(2)
  560.             PRINT*,'Y1  :',X(3) ,'   DY1  :',DX(3),'ERR :',DELAFON(3)
  561.             PRINT*,'Y2  :',X(4) ,'   DY2  :',DX(4),'ERR :',DELAFON(4)
  562.             PRINT*,'Z1  :',X(5) ,'   DZ1  :',DX(5),'ERR :',DELAFON(5)
  563.             PRINT*,'Z2  :',X(6) ,'   DZ2  :',DX(6),'ERR :',DELAFON(6)
  564.             PRINT*,'H1',H1,'H2',H2,ht,ht2
  565.             PRINT*
  566.             write(*,200)jacob
  567.           END IF
  568.  
  569.             END DO
  570.  
  571. ***     pour z2   ****
  572. ********************
  573.             fctd2=0.
  574.             If ( x(4) .Ge. Y02 ) Then
  575. **      si z converge est > a x(6) alors fctd est < 0
  576.                Fctd2 =-( 1.0+(A2*(-Y02 + X(4)))**B2-X(6) )
  577.                if (fctd2.ge.0.0 .and. x(6).eq.z2) fctd2 =0.0
  578.             Else
  579.                if (x(4) .le. y002 ) then
  580.                   fctd2=x(6)-z2
  581. c                  x(6)=z2
  582.                else
  583.                   Fctd2 =x(6)-z2
  584.                endif
  585.             Endif
  586.             If ( Fctd2 .lt. 0.0) then
  587.                Xg2=Xd2
  588.                if (passage2) pas2=pas2*2.0
  589.                Xd2=Xd2+pas2
  590.             Else
  591.                Xd2=Xg2
  592.                Pas2=Pas2/2.0
  593.                Xd2=Xd2+pas2
  594.                passage2=.false.
  595.             Endif
  596.             Erreurd2=(fctd2*fctd2)/(z2*z2)
  597.             if (x(6) .gt. zmax) then
  598.                x(6)=zmax
  599.                Erreurd2=0.0
  600.             endif
  601.  
  602. c            print*,'Erreurd2'
  603. c            print*,Erreurd2,x(4),y02,y002
  604.  
  605.          ENDDO
  606.  
  607. ***     pour z1   ****
  608. ********************
  609.          fctd1=0.0
  610.          If ( x(3) .Ge. Y01 ) Then
  611. **      si z converge est > a x(5) alors fctd1 est < 0
  612.             Fctd1 =-( 1.0+(A1*(-Y01 + X(3)))**B1-X(5) )
  613.             if (fctd1.ge.0.0 .and. x(5).eq.z1) fctd1=0.0
  614.          Else
  615.             if (x(3) .le. y001) then
  616.                fctd1=x(5)-z2
  617.                x(5)=z1
  618.             else
  619.                Fctd1 =x(5)-z2
  620.             endif
  621.          Endif
  622.          If ( Fctd1 .lt. 0.0) then
  623.             Xg1=Xd1
  624.             if (passage1) pas1=pas1*2.0
  625.             Xd1=Xd1+pas1
  626.          Else
  627.             Xd1=Xg1
  628.             Pas1=Pas1/2.0
  629.             Xd1=Xd1+pas1
  630.             passage1=.false.
  631.          Endif
  632.          Erreurd1=(Fctd1*Fctd1)/(z1*z1)
  633.          if (x(5) .gt. zmax) then
  634.             x(5)=zmax
  635.             Erreurd1=0.0
  636.          endif
  637.  
  638. c         print*,'Erreurd1'
  639. c         print*,Erreurd1,x(3),y02,Y002
  640.  
  641.       ENDDO
  642.  
  643.       IF (ERREUR .GT. TEST ) THEN
  644.          CONV=.FALSE.
  645.       ELSE
  646.          CONV=.TRUE.
  647.          SIG1= X(1)
  648.          SIG2= X(2)
  649.          Y1    = X(3)
  650.          Y2    = X(4)
  651.          Z1    = X(5)
  652.          Z2    = X(6)
  653.          IF (SIG1 .GE. 0.D0) THEN
  654.             JSIGN(1)=1
  655.          ELSE
  656.             JSIGN(1)=0
  657.          END IF
  658.          IF (SIG2.GE. 0.D0)THEN
  659.             JSIGN(2)=1
  660.          ELSE
  661.             JSIGN(2)=0
  662.          END IF
  663.          JSIGN(3)=1
  664.       END IF
  665.       RETURN
  666. 200   format(6(6(2x,e11.4),/))
  667.       END
  668.  
  669.  
  670.  
  671.  
  672.  
  673.  
  674.  
  675.  
  676.  
  677.  
  678.  
  679.  
  680.  
  681.  
  682.  
  683.  

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