Télécharger la plaquette Cast3M


fr en

Télécharger modjf2.eso

Retour à la liste

Numérotation des lignes : 

modjf2

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

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