Télécharger la plaquette Cast3M


fr en

Télécharger codjfd.eso

Retour à la liste

Numérotation des lignes : 

codjfd

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

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