Télécharger la plaquette Cast3M


fr en

Télécharger behav3.eso

Retour à la liste

Numérotation des lignes : 

behav3

  1. C BEHAV3    SOURCE    CB215821  26/08/24    21:15:15     12622          
  2.       SUBROUTINE BEHAV3(SIGR,DSTRN,DEPS1,DEPS2,IPLA,
  3.      & SIG1,SIG2,IFIS,SIGF,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP,
  4.      & BETJEF,VISCO,NECH0,NECH1)
  5. C
  6. C ==================================================================
  7. C
  8. C          MODELE DE PLASTICITE EN TRAITEMENT PRE ET POST PIC
  9. C            Deux criteres de traction compression:
  10. C                   Rankine et Von Mises
  11. C
  12. C ==================================================================
  13. C CE SOUS-PROGRAMME EST APPELE DANS "BONE".
  14. C
  15.       IMPLICIT INTEGER(I-N)
  16.       IMPLICIT REAL*8(A-H,O-Z)
  17.       DIMENSION DFSIG1(4),DFSIG2(4),SIGF(4),DSIGT(4),VEC1(4),VEC2(4)
  18.       DIMENSION V1(4),AC1(4),AC2(4),SIGE(4),SIGP(4)
  19.       DIMENSION DLAM1(2),DLAM0(2),SIGR(4),IS1(2)
  20.       DIMENSION FCRI0(2),FCRI1(2),DEPSI(4),DSTRN(4),SIGRV(4)
  21.       DIMENSION D2FSI1(4,4),D2FSI2(4,4),DEP(4,4),P1(4,4),P2(4,4)
  22.       DIMENSION D(4,4),A(4,4),AI(4,4),DJAC0(4,4),DJI(4,4),AH(4,4)
  23. C
  24. *
  25.       SEGMENT BETJEF
  26.        REAL*8 AA,BETA,COLI,PALF,YOUN,XNU,GFC,GFT,CAR,ETA,TDEF,
  27.      &          TCON,DPSTF1,DPSTF2,TETA,PDT,TP0
  28.        INTEGER ICT,ICC,IMOD,IVIS,ITR,
  29.      &                ISIM,IBB,IGAU,IZON
  30.       ENDSEGMENT
  31.       SEGMENT VISCO
  32.        REAL*8 DPSTV1,DPSTV2,SIGV1,SIGV2,ENDV
  33.       ENDSEGMENT
  34.       SEGMENT NECH0
  35.        REAL*8 DT,DC,ALFG,S0,ENDO
  36.       ENDSEGMENT
  37.       SEGMENT NECH1
  38.        REAL*8 ENDL
  39.       ENDSEGMENT
  40. *
  41. *
  42. *     COMMON /DBETJEF/AA,BETA,COLI,PALF,YOUN,XNU,GFC,GFT,CAR,ETA,TDEF,
  43. *    &              TCON,DPSTF1,DPSTF2,TETA,PDT,ICT,ICC,IMOD,IVIS,ITR,
  44. *    &                ISIM,IBB,IGAU,IZON
  45. C
  46.       MNSTRS=NSTRS
  47. C
  48.       IAPEX=0
  49.       Ft=PALF*COLI
  50.       Rb=COLI
  51.       PRB=1.D-10
  52.       PRB2=1.D-6
  53.       ITER=1
  54.       IEC1=0
  55.       IEC2=0
  56.       IBROY=0
  57.       ITANG=0
  58.       CRIMAX=0.D0
  59.       SEQ = 0.D0
  60.       SEQ1 = 0.D0
  61.       SEQ2 = 0.D0
  62.       CALL ZERO(SIGE,4,1)
  63.       CALL ZERO(V1,4,1)
  64.       CALL ZERO(P1,4,4)
  65.       CALL ZERO(P2,4,4)
  66.       CALL ZERO(D,4,4)
  67. C
  68.       DO 10 I=1,NSTRS
  69.       SIGE(I)=SIGF(I)
  70.    10 CONTINUE
  71. C
  72.       IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  73.       AD=YOUN/(1.D0-XNU*XNU)
  74.       D(1,1)=AD
  75.       D(2,2)=D(1,1)
  76.       D(3,3)=AD*(1.D0-XNU)/2.D0
  77.       D(1,2)=AD*XNU
  78.       D(2,1)=D(1,2)
  79.       ENDIF
  80. C
  81.       IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  82.       ADD=YOUN/((1.D0+XNU)*(1.D0-2.D0*XNU))
  83.       D(1,1)=ADD*(1.D0-XNU)
  84.       D(2,2)=D(1,1)
  85.       D(3,3)=D(1,1)
  86.       D(1,2)=ADD*XNU
  87.       D(2,1)=D(1,2)
  88.       D(1,3)=D(1,2)
  89.       D(2,3)=D(1,2)
  90.       D(3,1)=D(1,2)
  91.       D(3,2)=D(1,2)
  92.       D(4,4)=0.5*ADD*(1.D0-2.D0*XNU)
  93.       ENDIF
  94. C
  95. C ************ Le point est fissure pour 1ere fois ***********
  96. C
  97.       IF (IFIS.EQ.0) THEN
  98.         CALL PRINC(SIGF,V1,NSTRS)
  99.         DEPS1=0.D0
  100.         SIG1=Ft
  101.         IF (V1(1).GT.Ft) THEN
  102.          IFIS=1
  103.         ENDIF
  104.       ENDIF
  105.       IF (IPLA.EQ.0) THEN
  106.        IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  107.         CALL CRIVON(SIGF,SEQ,NSTRS,BETJEF)
  108.        ENDIF
  109.        IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN
  110.         CALL CRIDRU(SIGF,SEQ,NSTRS,BETJEF)
  111.        ENDIF
  112.         DEPS2=0.D0
  113.         SIG2=Rb*AA
  114.         IF (SEQ.GT.SIG2) THEN
  115.          IPLA=1
  116.         ENDIF
  117.       ENDIF
  118. C
  119.    12 CONTINUE
  120. C
  121.       CALL PRINC(SIGE,V1,NSTRS)
  122.       TETA=V1(4)
  123.        IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  124.         CALL CRIVON(SIGE,SEQ,NSTRS,BETJEF)
  125.        ENDIF
  126.        IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN
  127.         CALL CRIDRU(SIGE,SEQ,NSTRS,BETJEF)
  128.        ENDIF
  129.       FCRI0(1)=V1(1)-SIG1
  130.       FCRI0(2)=SEQ-SIG2
  131.       IF (ABS(FCRI0(1)).GT.ABS(FCRI0(2))) THEN
  132.        CRIMAX=ABS(100.D0*FCRI0(1))
  133.       ELSE
  134.        CRIMAX=ABS(100.D0*FCRI0(2))
  135.       ENDIF
  136.       IF (FCRI0(1).LT.0.D0.AND.FCRI0(2).LT.0.D0) THEN
  137.         IEC1=1
  138.         IEC2=1
  139.       ENDIF
  140.       IF (FCRI0(1).LT.0.D0.AND.FCRI0(2).GE.0.D0) THEN
  141.         IEC1=1
  142.         IEC2=0
  143.       ENDIF
  144.       IF (FCRI0(1).GE.0.D0.AND.FCRI0(2).LT.0.D0) THEN
  145.         IEC1=0
  146.         IEC2=1
  147.       ENDIF
  148. C
  149.    15 CONTINUE
  150. C
  151.       DO 16 I=1,NSTRS
  152.       SIGF(I)=SIGE(I)
  153.    16 CONTINUE
  154. C
  155.       IF (IEC1.EQ.1.AND.IEC2.EQ.1) THEN
  156.       DO 101 I=1,NSTRS
  157.       DO 17 J=1,NSTRS
  158.       DEP(I,J)=D(I,J)
  159.    17 CONTINUE
  160.   101 CONTINUE
  161.         GOTO 100
  162.       ENDIF
  163. C
  164.       IF (IEC1.EQ.1.AND.IEC2.EQ.0) THEN
  165.         CALL BEHAV2(SIGR,DSTRN,DEPS1,DEPS2,IPLA
  166.      & ,SIG1,SIG2,IFIS,SIGF,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP,
  167.      &  BETJEF,VISCO,NECH0,NECH1)
  168.         GOTO 100
  169.       ENDIF
  170. C
  171.       IF (IEC1.EQ.0.AND.IEC2.EQ.1) THEN
  172.         CALL BEHAV1(SIGR,DSTRN,DEPS1,DEPS2,IPLA
  173.      & ,SIG1,SIG2,IFIS,SIGF,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP,
  174.      &  BETJEF,VISCO,NECH0,NECH1)
  175.         GOTO 100
  176.       ENDIF
  177. C
  178.       DK1=DEPS1
  179.       DK2=DEPS2
  180.       DLAM0(1)=0.D0
  181.       DLAM0(2)=0.D0
  182. C
  183.    18 CONTINUE
  184. C
  185. C ************ Preparation des criteres  ********************
  186. C
  187. C ------------------ Von Mises -----------------------------
  188.       IF (IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  189.         A1=0.D0
  190.         A2=1.D0
  191.       ENDIF
  192. C ----------------- Drucker Prager -------------------------
  193.       IF (IMOD.EQ.3.OR.IMOD.EQ.4) THEN
  194.         A1=(BETA-1.D0)/(2*BETA-1.D0)
  195.         A2=BETA/(2*BETA-1.D0)
  196.       ENDIF
  197. C ----------------------------------------------------------
  198. C
  199. C ---------------- direction de traction ------------------
  200. C
  201. C--------------- Determination du parametre d'ecrouissage ----------
  202. C
  203.       IF(IVIS.LE.2) THEN
  204.        CALL UNICOU(DK1,PAECT,1,SEQ1,BETJEF)
  205.       ELSE
  206.        CALL UNICO1(DK1,PAECT,1,SEQ1,BETJEF,NECH0,NECH1)
  207.       ENDIF
  208. C
  209. C------------------------------------------------------------
  210.       PI=4.D0*ATAN(1.D0)
  211.       PHIC=V1(4)*(PI/180.D0)
  212.       COSA=COS(PHIC)
  213.       SINA=SIN(PHIC)
  214. C
  215.       IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  216.           DFSIG1(1)=COSA*COSA
  217.           DFSIG1(2)=SINA*SINA
  218.           DFSIG1(3)=2.D0*SINA*COSA
  219.       DO 102 I=1,NSTRS
  220.       DO 19 J=1,NSTRS
  221.        P1(I,J)=0.D0
  222.    19 CONTINUE
  223.   102 CONTINUE
  224.       P1(1,1)=1.D0/2.D0
  225.       P1(1,2)=-1.D0/2.D0
  226.       P1(2,1)=-1.D0/2.D0
  227.       P1(2,2)=1.D0/2.D0
  228.       P1(3,3)=2.D0
  229.       ENDIF
  230. C
  231.       IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  232.           DFSIG1(1)=COSA*COSA
  233.           DFSIG1(2)=SINA*SINA
  234.           DFSIG1(3)=0.D0
  235.           DFSIG1(4)=2.D0*SINA*COSA
  236.       DO 103 I=1,NSTRS
  237.       DO 21 J=1,NSTRS
  238.        P1(I,J)=0.D0
  239.    21 CONTINUE
  240.   103 CONTINUE
  241.       P1(1,1)=1.D0/2.D0
  242.       P1(1,2)=-1.D0/2.D0
  243.       P1(2,1)=-1.D0/2.D0
  244.       P1(2,2)=1.D0/2.D0
  245.       P1(4,4)=2.D0
  246.       ENDIF
  247. C
  248.       DO 104 I=1,NSTRS
  249.       AC1(I)=0.D0
  250.       DO 20 J=1,NSTRS
  251.       AC1(I)=AC1(I)+D(I,J)*DFSIG1(J)
  252.    20 CONTINUE
  253.   104 CONTINUE
  254. C
  255. C ---------------- direction de compression ----------------
  256. C
  257.       IF(IVIS.LE.2) THEN
  258.        CALL UNICOU(DK2,PAEC,2,SEQ2,BETJEF)
  259.       ELSE
  260.        CALL UNICO1(DK2,PAEC,2,SEQ2,BETJEF,NECH0,NECH1)
  261.       ENDIF
  262. C
  263.       IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  264.       SX=SIGF(1)
  265.       SY=SIGF(2)
  266.       SXY=SIGF(3)
  267.       R1=(SX-SY)*(SX-SY)+SX*SX+SY*SY+6*SXY*SXY
  268.       R=SQRT(0.5*R1)
  269.       IF (ABS(R).LE.1.D-10) THEN
  270.       WRITE(*,*)'INDETERMINATION de DFSIG ds BEHAV3'
  271.       STOP
  272.       ENDIF
  273.           DFSIG2(1)=(2.D0*SX-SY)/2.D0/R/A2+A1/A2
  274.           DFSIG2(2)=(2.D0*SY-SX)/2.D0/R/A2+A1/A2
  275.           DFSIG2(3)=3.D0*SXY/R/A2
  276.       P2(1,1)=2.D0/A2/A2
  277.       P2(1,2)=-1.D0/A2/A2
  278.       P2(1,3)=0.D0/A2/A2
  279.       P2(2,1)=-1.D0/A2/A2
  280.       P2(2,2)=2.D0/A2/A2
  281.       P2(2,3)=0.D0/A2/A2
  282.       P2(3,1)=0.D0/A2/A2
  283.       P2(3,2)=0.D0/A2/A2
  284.       P2(3,3)=6.D0/A2/A2
  285.       ENDIF
  286. C
  287.       IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  288.       SX=SIGF(1)
  289.       SY=SIGF(2)
  290.       SZ=SIGF(3)
  291.       SXY=SIGF(4)
  292.       TX=(SX-SY)*(SX-SY)
  293.       TY=(SX-SZ)*(SX-SZ)
  294.       TZ=(SY-SZ)*(SY-SZ)
  295.       R1=TX+TY+TZ+6*SXY*SXY
  296.       R=SQRT(0.5*R1)
  297.       IF (ABS(R).LE.1.D-10) THEN
  298.       WRITE(*,*)'INDETERMINATION de DFSIG ds BEHAV3'
  299.       STOP
  300.       ENDIF
  301.           DFSIG2(1)=(2.D0*SX-SY-SZ)/2.D0/R/A2+A1/A2
  302.           DFSIG2(2)=(2.D0*SY-SX-SZ)/2.D0/R/A2+A1/A2
  303.           DFSIG2(3)=(2.D0*SZ-SX-SY)/2.D0/R/A2+A1/A2
  304.           DFSIG2(4)=3.D0*SXY/R/A2
  305.       P2(1,1)=2.D0/A2/A2
  306.       P2(1,2)=-1.D0/A2/A2
  307.       P2(1,3)=-1.D0/A2/A2
  308.       P2(1,4)=0.D0/A2/A2
  309.       P2(2,1)=-1.D0/A2/A2
  310.       P2(2,2)=2.D0/A2/A2
  311.       P2(2,3)=-1.D0/A2/A2
  312.       P2(2,4)=0.D0/A2/A2
  313.       P2(3,1)=-1.D0/A2/A2
  314.       P2(3,2)=-1.D0/A2/A2
  315.       P2(3,3)=2.D0/A2/A2
  316.       P2(3,4)=0.D0/A2/A2
  317.       P2(4,1)=0.D0/A2/A2
  318.       P2(4,2)=0.D0/A2/A2
  319.       P2(4,3)=0.D0/A2/A2
  320.       P2(4,4)=6.D0/A2/A2
  321.       ENDIF
  322. C
  323.       DO 105 I=1,NSTRS
  324.       AC2(I)=0.D0
  325.       DO 25 J=1,NSTRS
  326.       AC2(I)=AC2(I)+D(I,J)*DFSIG2(J)
  327.    25 CONTINUE
  328.   105 CONTINUE
  329. C
  330. C ************ Preparation du jacobien  ********************
  331. C
  332.       F1DF1=0.D0
  333.       DO 30 J=1,NSTRS
  334.       F1DF1=F1DF1+DFSIG1(J)*AC1(J)
  335.    30 CONTINUE
  336.       F1DF2=0.D0
  337.       DO 31 J=1,NSTRS
  338.       F1DF2=F1DF2+DFSIG1(J)*AC2(J)
  339.    31 CONTINUE
  340.       F2DF2=0.D0
  341.       DO 32 J=1,NSTRS
  342.       F2DF2=F2DF2+DFSIG2(J)*AC2(J)
  343.    32 CONTINUE
  344.       F2DF1=0.D0
  345.       DO 33 J=1,NSTRS
  346.       F2DF1=F2DF1+DFSIG2(J)*AC1(J)
  347.    33 CONTINUE
  348. C
  349.       DJAC0(1,1)=-(F1DF1+PAECT)
  350.       DJAC0(1,2)=-(F1DF2)
  351.       DJAC0(2,2)=-(F2DF2+PAEC)
  352.       DJAC0(2,1)=-(F2DF1)
  353. C
  354.       DO 106 I=1,2
  355.       DO 43 J=1,2
  356.       DJI(I,J)=DJAC0(I,J)
  357.    43 CONTINUE
  358.   106 CONTINUE
  359. C
  360.       CALL INVMA2(DJI,2,ISING)
  361.       IF (ISING.EQ.1) THEN
  362.        WRITE(*,*)'MATRICE DJI singuliere ds BEHAV3'
  363.       ENDIF
  364. C
  365. C ************ Debut des iterations internes *******************
  366. C
  367.    40 CONTINUE
  368. C
  369. C *************** Determination de DK et DLAM ******************
  370. C
  371. C
  372. C
  373.       DLAM1(1)=DLAM0(1)-DJI(1,1)*FCRI0(1)-DJI(1,2)*FCRI0(2)
  374.       DLAM1(2)=DLAM0(2)-DJI(2,1)*FCRI0(1)-DJI(2,2)*FCRI0(2)
  375. C
  376.       IF (DLAM1(1).LE.0.D0) IEC1=1
  377.       IF (DLAM1(2).LE.0.D0) IEC2=1
  378.       IF (IEC1.EQ.1.OR.IEC2.EQ.1) THEN
  379. C       WRITE(*,*)'Dans BEHAV3, DLAMDA1 est negatif:',DLAM1(1)
  380. C       WRITE(*,*)'Dans BEHAV3, DLAMDA2 est negatif:',DLAM1(2)
  381. C       WRITE(*,*)'A l iteration :',ITER
  382.        GOTO 15
  383.       ENDIF
  384. C
  385.       DK1=DEPS1+DLAM1(1)
  386.       DK2=DEPS2+DLAM1(2)
  387. C
  388.       IF(IVIS.LE.2) THEN
  389.        CALL UNICOU(DK1,PAECT,1,SEQ1,BETJEF)
  390.        CALL UNICOU(DK2,PAEC,2,SEQ2,BETJEF)
  391.       ELSE
  392.        CALL UNICO1(DK1,PAECT,1,SEQ1,BETJEF,NECH0,NECH1)
  393.        CALL UNICO1(DK2,PAEC,2,SEQ2,BETJEF,NECH0,NECH1)
  394.       ENDIF
  395. C
  396. C ************** Determination de DPHI *********************
  397. C
  398.       IF (IMOD.EQ.1) THEN
  399.       AD1=YOUN/(1.D0-XNU)
  400.       coef=2.D0*SEQ2/(2.D0*SEQ2+DLAM1(2)*AD1)
  401.       DPHI1=SEQ1-0.5*coef*((SIGE(1)+SIGE(2))-AD1*DLAM1(1))
  402.       TO=SIGF(3)
  403.       DPHI2=SEQ2
  404.       ENDIF
  405.       IF (IMOD.EQ.2) THEN
  406.       DG=YOUN/2.D0/(1.D0+XNU)
  407.       AD1=YOUN/(1.D0+XNU)/(1.D0-2.D0*XNU)
  408.       coef1=(SEQ2+2.D0*DLAM1(2)*DG)/(SEQ2+3.D0*DLAM1(2)*DG)
  409.       coef2=(2.D0*DLAM1(2)*DG)/(SEQ2+3.D0*DLAM1(2)*DG)
  410.       CAS1=(SIGE(1)+SIGE(2)-0.5*DLAM1(1)*2.D0*AD1)
  411.       CAS2=SIGE(3)-DLAM1(1)*XNU*AD1
  412.       DPHI1=SEQ1-0.5*(coef1*CAS1+coef2*CAS2)
  413.       TO=SIGF(4)
  414.       DPHI2=SEQ2
  415.       ENDIF
  416.       IF (IMOD.EQ.3) THEN
  417.       DPHI2=SEQ2-A1/A2*(SIGF(1)+SIGF(2))
  418.       AD1=YOUN/(1.D0-XNU)
  419.       coef=2.D0*DPHI2/(2.D0*DPHI2+DLAM1(2)*AD1)
  420.       SOM=A1/A2*DLAM1(2)*2.D0*AD1
  421.       DPHI1=SEQ1-0.5*coef*((SIGE(1)+SIGE(2))-AD1*DLAM1(1)-SOM)
  422.       TO=SIGF(3)
  423.       ENDIF
  424.       IF (IMOD.EQ.4) THEN
  425.       SOM=SIGE(1)+SIGE(2)+SIGE(3)
  426.       DPHI2=A1*SOM/A2
  427.       DPHI2=DPHI2-A1*DLAM1(1)*ADD*(1.D0+XNU)/A2
  428.       DPHI2=DPHI2-A1*A1*DLAM1(2)*3.D0*ADD*(1.D0+XNU)/A2/A2
  429.       DPHI2=SEQ2-DPHI2
  430.       DG=YOUN/2.D0/(1.D0+XNU)/A2/A2
  431.       AD1=YOUN/(1.D0+XNU)/(1.D0-2.D0*XNU)
  432.       AD2=YOUN/(1.D0-2.D0*XNU)
  433.       coef1=(DPHI2+2.D0*DLAM1(2)*DG)/(DPHI2+3.D0*DLAM1(2)*DG)
  434.       coef2=(2.D0*DLAM1(2)*DG)/(DPHI2+3.D0*DLAM1(2)*DG)
  435.       CAS1=(SIGE(1)+SIGE(2)-DLAM1(1)*AD1-2.D0*A1/A2*DLAM1(2)*AD2)
  436.       CAS2=SIGE(3)-DLAM1(1)*XNU*AD1-A1/A2*DLAM1(2)*AD2
  437.       DPHI1=SEQ1-0.5*(coef1*CAS1+coef2*CAS2)
  438.       TO=SIGF(4)
  439.       ENDIF
  440. C      IF (DPHI1.LT.0.D0) THEN
  441. C       WRITE(*,*)'ATTENTION DPHI1 NEGATIF'
  442. C       WRITE(*,*)'DPHI1=',DPHI1
  443. C      ENDIF
  444. C      IF (DPHI2.LT.0.D0) THEN
  445. C       WRITE(*,*)'ATTENTION DPHI2 ds BEHAV3 NEGATIF'
  446. C       WRITE(*,*)'DPHI2=',DPHI2
  447. C      ENDIF
  448. C
  449. C--------------- Cas de l'apex -----------------------------
  450. C
  451.       IF (ABS(DPHI1).LE.10E-10.AND.
  452.      *ABS(DPHI2).GT.10E-10) THEN
  453.       IAPEX=1
  454. C      WRITE(*,*)'IAPEX ds BEHAV3 =',IAPEX
  455. C      WRITE(*,*)'Dans l element',IBB
  456. C      WRITE(*,*)'et au point d intégration',IGAU
  457.       DO 107 I=1,NSTRS
  458.       DO 50 J=1,NSTRS
  459.       AI(I,J)=0.D0
  460.    50 CONTINUE
  461.   107 CONTINUE
  462.         AI(1,1)=0.5
  463.         AI(1,2)=0.5
  464.         AI(2,1)=AI(1,2)
  465.         AI(2,2)=AI(1,1)
  466.         IF (IMOD.EQ.1.OR.IMOD.EQ.3) AI(3,3)=0.D0
  467.         IF (IMOD.EQ.2.OR.IMOD.EQ.4) AI(3,3)=1.D0
  468.         GOTO 75
  469.       ENDIF
  470. C
  471. C---------- Cas du critere reduit a un point ----------------
  472. C
  473.       IF (ABS(DPHI2).LE.10E-10) THEN
  474.       IAPEX=2
  475. C      WRITE(*,*)'IAPEX ds BEHAV3 =',IAPEX
  476. C      WRITE(*,*)'Dans l element',IBB
  477. C      WRITE(*,*)'et au point d intégration',IGAU
  478.       DO 108 I=1,NSTRS
  479.       DO 52 J=1,NSTRS
  480.       AI(I,J)=0.D0
  481.    52 CONTINUE
  482.   108 CONTINUE
  483.       IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  484.       DO 109 I=1,3
  485.       DO 53 J=1,3
  486.       AI(I,J)=1.D0/3.D0
  487.    53 CONTINUE
  488.   109 CONTINUE
  489.       ENDIF
  490.         GOTO 75
  491.       ENDIF
  492. C
  493. C ************** Mise a jour des contraintes ***************
  494. C
  495. C ---------------- calcul de la matrice A ------------------
  496. C
  497.       DO 110 I=1,NSTRS
  498.       DO 60 J=1,NSTRS
  499.       A(I,J)=0.D0
  500.    60 CONTINUE
  501.   110 CONTINUE
  502. C
  503.       DG=YOUN/2.D0/(1.D0+XNU)
  504. C
  505.       IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  506.       DW1=AD*(2.D0-XNU)
  507.       DW2=AD*(2.D0*XNU-1.D0)
  508.       A(1,1)=1.D0+(DLAM1(1)*DG)/2.D0/DPHI1
  509.            A(1,1)=A(1,1)+(DLAM1(2)*DW1)/2.D0/DPHI2/A2/A2
  510.       A(2,2)=A(1,1)
  511.       A(1,2)=-(DLAM1(1)*DG)/2.D0/DPHI1
  512.            A(1,2)=A(1,2)+(DLAM1(2)*DW2)/2.D0/DPHI2/A2/A2
  513.       A(2,1)=A(1,2)
  514.       A(3,3)=(DLAM1(1)*DG)/DPHI1
  515.       A(3,3)=A(3,3)+1.D0+(DLAM1(2)*3.D0*DG)/DPHI2/A2/A2
  516.       ENDIF
  517. C
  518.       IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  519.       A(1,1)=1.D0+(DLAM1(1)*DG)/2.D0/DPHI1
  520.            A(1,1)=A(1,1)+(DLAM1(2)*2.D0*DG)/DPHI2/A2/A2
  521.       A(2,2)=A(1,1)
  522.       A(3,3)=1.D0+(DLAM1(2)*2.D0*DG)/DPHI2/A2/A2
  523.       A(1,2)=-(DLAM1(1)*DG)/2.D0/DPHI1
  524.            A(1,2)=A(1,2)-(DLAM1(2)*DG)/DPHI2/A2/A2
  525.       A(2,1)=A(1,2)
  526.       A(1,3)=-(DLAM1(2)*DG)/DPHI2/A2/A2
  527.       A(3,1)=A(1,3)
  528.       A(3,2)=A(1,3)
  529.       A(2,3)=A(1,3)
  530.       A(4,4)=1.D0+(DLAM1(1)*DG)/DPHI1
  531.            A(4,4)=A(4,4)+(DLAM1(2)*DG*3.D0)/DPHI2/A2/A2
  532.       ENDIF
  533. C
  534. C -------------- invertion de la matrice A -----------------
  535. C
  536.       DO 111 I=1,NSTRS
  537.       DO 70 J=1,NSTRS
  538.       AI(I,J)=A(I,J)
  539.    70 CONTINUE
  540.   111 CONTINUE
  541. C
  542.       CALL INVMA2(AI,NSTRS,ISING)
  543.       IF (ISING.EQ.1) THEN
  544.        WRITE(*,*)'MATRICE AI singuliere ds BEHAV3'
  545.       ENDIF
  546. C
  547. C -------------- mise a jour des contraintes  ------------
  548. C
  549.    75 CONTINUE
  550.       VEC1(1)=AD1
  551.       VEC1(2)=AD1
  552.       VEC1(3)=0.D0
  553.       VEC1(4)=0.D0
  554.       IF(IMOD.EQ.3) THEN
  555.        VEC2(1)=AD1*(1.D0+XNU)
  556.        VEC2(2)=AD1*(1.D0+XNU)
  557.        VEC2(3)=0.D0
  558.        VEC2(4)=0.D0
  559.       ENDIF
  560.       IF(IMOD.EQ.4) THEN
  561.        VEC2(1)=AD1*(1.D0+XNU)
  562.        VEC2(2)=AD1*(1.D0+XNU)
  563.        VEC2(3)=AD1*(1.D0+XNU)
  564.        VEC2(4)=0.D0
  565.       ENDIF
  566.       IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  567.        VEC2(1)=0.D0
  568.        VEC2(2)=0.D0
  569.        VEC2(3)=0.D0
  570.        VEC2(4)=0.D0
  571.       ENDIF
  572. C
  573.       DO 80 I=1,NSTRS
  574.       DEPSI(I)=SIGE(I)-0.5*DLAM1(1)*VEC1(I)
  575.       DEPSI(I)=DEPSI(I)-A1*DLAM1(2)*VEC2(I)/A2
  576.    80 CONTINUE
  577. C
  578.       DO 112 I=1,NSTRS
  579.       SIGF(I)=0.0D+00
  580.       DO 90 J=1,NSTRS
  581.       SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J)
  582.    90 CONTINUE
  583.   112 CONTINUE
  584. C
  585. C ******** Verification des criteres ****************
  586. C
  587.       CALL PRINC(SIGF,V1,NSTRS)
  588.       IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  589.        CALL CRIVON(SIGF,SEQ,NSTRS,BETJEF)
  590.       ENDIF
  591.       IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN
  592.        CALL CRIDRU(SIGF,SEQ,NSTRS,BETJEF)
  593.       ENDIF
  594.       FCRI1(1)=V1(1)-SEQ1
  595.       IF (IAPEX.EQ.2) FCRI1(1)=0.D0
  596.       FCRI1(2)=SEQ-SEQ2
  597. C
  598.       IF (IBROY.EQ.0.AND.(ABS(FCRI1(1)).GE.CRIMAX.OR
  599.      *      .ABS(FCRI1(2)).GE.CRIMAX)) THEN
  600. C      WRITE(*,*)'****************************************'
  601. C      WRITE(*,*)'LE RESIDU DIVERGE AVEC BROYDEN'
  602. C      WRITE(*,*)'on passe donc a la secante'
  603. C      WRITE(*,*)'Dans l element',IBB
  604. C      WRITE(*,*)'et au point d intégration',IGAU
  605. C      WRITE(*,*)'CRIMAX=',CRIMAX
  606. C      WRITE(*,*)'****************************************'
  607.       ITER=ITR
  608.       ENDIF
  609. C
  610. C *******  Compteur sur la methode de resolution ****
  611. C
  612.       IF (IBROY.EQ.0.AND.ITER.EQ.ITR) THEN
  613.         IBROY=1
  614.         ITANG=1
  615.         ITER=1
  616.         IAPEX=0
  617.         IEC1=0
  618.         IEC2=0
  619.         GOTO 12
  620.       ENDIF
  621. C
  622. C ******* non convergence **************************
  623. C
  624.       IF ((ABS(FCRI1(1)).GT.PRB.OR.ABS(FCRI1(2)).GT.PRB)
  625.      *.AND.ITER.LT.ITR) THEN
  626.       IF (IBROY.EQ.0) THEN
  627.        CALL BROYDI(DJI,DLAM0,DLAM1,FCRI0,FCRI1)
  628.        DLAM0(1)=DLAM1(1)
  629.        DLAM0(2)=DLAM1(2)
  630.        FCRI0(1)=FCRI1(1)
  631.        FCRI0(2)=FCRI1(2)
  632.        IF (ITER.GE.(ITR-1)) THEN
  633. C        WRITE(*,*)'***********************'
  634. C        WRITE(*,*)'BROYDEN n a pas aboutit'
  635. C        WRITE(*,*)'ITER=',ITER
  636. C        WRITE(*,*)'FCRIT=',FCRI0(1)
  637. C        WRITE(*,*)'FCRIC=',FCRI0(2)
  638. C        WRITE(*,*)'Dans l element',IBB
  639. C        WRITE(*,*)'et au point d intégration',IGAU
  640. C        WRITE(*,*)'***********************'
  641.        ENDIF
  642.        ITER=ITER+1
  643.        GOTO 40
  644.       ENDIF
  645.       IF (IBROY.EQ.1.AND.ITANG.EQ.1) THEN
  646.        DLAM0(1)=DLAM1(1)
  647.        DLAM0(2)=DLAM1(2)
  648.        FCRI0(1)=FCRI1(1)
  649.        FCRI0(2)=FCRI1(2)
  650.        CALL PRINC(SIGF,V1,NSTRS)
  651.        IF (ITER.GE.(ITR-5)) THEN
  652. C        WRITE(*,*)'ITER=',ITER
  653. C        WRITE(*,*)'FCRIT=',FCRI0(1)
  654. C        WRITE(*,*)'FCRIC=',FCRI0(2)
  655.        ENDIF
  656.        ITER=ITER+1
  657.        GOTO 18
  658.       ENDIF
  659.       ENDIF
  660.       IF (ITER.GE.ITR.AND.(ABS(FCRI1(1)).GT.PRB2.OR.
  661.      *      ABS(FCRI1(2)).GT.PRB2)) THEN
  662.         WRITE(*,*)'NON CONVERGENCE INTERNE dans BEHAV3'
  663.         WRITE(*,*)'Dans l element',IBB
  664.         WRITE(*,*)'et au point d intégration',IGAU
  665.         WRITE(*,*)'FCRIT=',FCRI0(1)
  666.         WRITE(*,*)'FCRIC=',FCRI0(2)
  667.         WRITE(*,*)'IPLA=',IPLA
  668.         WRITE(*,*)'IFIS=',IFIS
  669. C        STOP
  670.       ENDIF
  671. C
  672. C **************** Fin des iterations internes *******************
  673. C
  674. C
  675. C ************* Calcul de la DEP ****************
  676. C
  677.       IF (IAPEX.EQ.0) THEN
  678.       CALL DERI2(D2FSI1,DPHI1,P1,SIGF,NSTRS,BETJEF)
  679.       CALL DERI2(D2FSI2,DPHI2,P2,SIGF,NSTRS,BETJEF)
  680.       CALL LADEP1(SIGF,DEP,PAECT,PAEC,NSTRS,IFOU,
  681.      & DFSIG1,DFSIG2,D2FSI1,D2FSI2,DLAM1,BETJEF)
  682.       ELSE IF (IAPEX.EQ.1) THEN
  683.       DO 114 I=1,NSTRS
  684.       DO 113 J=1,NSTRS
  685.          AH(I,J)=0.D0
  686.       DO 95 K=1,NSTRS
  687.          AH(I,J)=AH(I,J)+AI(I,K)*D(K,J)
  688.    95 CONTINUE
  689.   113 CONTINUE
  690.   114 CONTINUE
  691.        DFSIG1(1)=0.5*SQRT(2.D0)
  692.        DFSIG1(2)=0.5*SQRT(2.D0)
  693.        DFSIG1(3)=0.D0
  694.        IF (IMOD.EQ.2) DFSIG1(4)=0.D0
  695.        CALL CREDEP(AH,DFSIG1,PAEC,NSTRS,DEP,BETJEF)
  696.       ELSE IF (IAPEX.EQ.2) THEN
  697.       DO 116 I=1,NSTRS
  698.       DO 115 J=1,NSTRS
  699.          DEP(I,J)=0.D0
  700.       DO 96 K=1,NSTRS
  701.          DEP(I,J)=DEP(I,J)+AI(I,K)*D(K,J)
  702.    96 CONTINUE
  703.   115 CONTINUE
  704.   116 CONTINUE
  705.       ENDIF
  706.       IAPEX=0
  707. C      DO 97 I=1,NSTRS
  708. C      DO 97 J=1,NSTRS
  709. C         DEP(I,J)=D(I,J)
  710. C   97 CONTINUE
  711. C
  712. C ***********************************************
  713. C
  714.       DEPS1=DK1
  715.       DEPS2=DK2
  716.       SIG1=SEQ1
  717.       SIG2=SEQ2
  718.       IF(IVIS.LE.2) THEN
  719.        CALL INDICA(DK1,DK2,IFIS,IPLA,3,BETJEF)
  720.       ELSE
  721.        CALL INDIC1(DK1,DK2,IFIS,IPLA,3,BETJEF,NECH0,NECH1)
  722.       ENDIF
  723. C ********************** CALCUL VISCOPLASTIQUE ***********************
  724.       IF (IVIS.EQ.1) THEN
  725.         CALL VISPLA(SIGR,SIGF,DSIGT,NSTRS,DEPS1,DEPS2,
  726.      &                  SIGP,SIGRV,DSTRN,D,BETJEF,VISCO)
  727.       ENDIF
  728. C ********************************************************************
  729.   100 CONTINUE
  730.       IEC1=0
  731.       IEC2=0
  732.       CRIMAX=0.D0
  733. C
  734.       RETURN
  735.       END
  736.  
  737.  
  738.  
  739.  
  740.  
  741.  
  742.  
  743.  

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