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  16/04/21    21:15:19     8920
  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 17 I=1,NSTRS
  157.       DO 17 J=1,NSTRS
  158.       DEP(I,J)=D(I,J)
  159.    17 CONTINUE
  160.         GOTO 100
  161.       ENDIF
  162. C
  163.       IF (IEC1.EQ.1.AND.IEC2.EQ.0) THEN
  164.         CALL BEHAV2(SIGR,DSTRN,DEPS1,DEPS2,IPLA
  165.      & ,SIG1,SIG2,IFIS,SIGF,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP,
  166.      &  BETJEF,VISCO,NECH0,NECH1)
  167.         GOTO 100
  168.       ENDIF
  169. C
  170.       IF (IEC1.EQ.0.AND.IEC2.EQ.1) THEN
  171.         CALL BEHAV1(SIGR,DSTRN,DEPS1,DEPS2,IPLA
  172.      & ,SIG1,SIG2,IFIS,SIGF,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP,
  173.      &  BETJEF,VISCO,NECH0,NECH1)
  174.         GOTO 100
  175.       ENDIF
  176. C
  177.       DK1=DEPS1
  178.       DK2=DEPS2
  179.       DLAM0(1)=0.D0
  180.       DLAM0(2)=0.D0
  181. C
  182.    18 CONTINUE
  183. C
  184. C ************ Preparation des criteres  ********************
  185. C
  186. C ------------------ Von Mises -----------------------------
  187.       IF (IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  188.         A1=0.D0
  189.         A2=1.D0
  190.       ENDIF
  191. C ----------------- Drucker Prager -------------------------
  192.       IF (IMOD.EQ.3.OR.IMOD.EQ.4) THEN
  193.         A1=(BETA-1.D0)/(2*BETA-1.D0)
  194.         A2=BETA/(2*BETA-1.D0)
  195.       ENDIF
  196. C ----------------------------------------------------------
  197. C
  198. C ---------------- direction de traction ------------------
  199. C
  200. C--------------- Determination du parametre d'ecrouissage ----------
  201. C
  202.       IF(IVIS.LE.2) THEN
  203.        CALL UNICOU(DK1,PAECT,1,SEQ1,BETJEF)
  204.       ELSE
  205.        CALL UNICO1(DK1,PAECT,1,SEQ1,BETJEF,NECH0,NECH1)
  206.       ENDIF
  207. C
  208. C------------------------------------------------------------
  209.       PI=4.D0*ATAN(1.D0)
  210.       PHIC=V1(4)*(PI/180.D0)
  211.       COSA=COS(PHIC)
  212.       SINA=SIN(PHIC)
  213. C
  214.       IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  215.           DFSIG1(1)=COSA*COSA
  216.           DFSIG1(2)=SINA*SINA
  217.           DFSIG1(3)=2.D0*SINA*COSA
  218.       DO 19 I=1,NSTRS
  219.       DO 19 J=1,NSTRS
  220.        P1(I,J)=0.D0
  221.    19 CONTINUE
  222.       P1(1,1)=1.D0/2.D0
  223.       P1(1,2)=-1.D0/2.D0
  224.       P1(2,1)=-1.D0/2.D0
  225.       P1(2,2)=1.D0/2.D0
  226.       P1(3,3)=2.D0
  227.       ENDIF
  228. C
  229.       IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  230.           DFSIG1(1)=COSA*COSA
  231.           DFSIG1(2)=SINA*SINA
  232.           DFSIG1(3)=0.D0
  233.           DFSIG1(4)=2.D0*SINA*COSA
  234.       DO 21 I=1,NSTRS
  235.       DO 21 J=1,NSTRS
  236.        P1(I,J)=0.D0
  237.    21 CONTINUE
  238.       P1(1,1)=1.D0/2.D0
  239.       P1(1,2)=-1.D0/2.D0
  240.       P1(2,1)=-1.D0/2.D0
  241.       P1(2,2)=1.D0/2.D0
  242.       P1(4,4)=2.D0
  243.       ENDIF
  244. C
  245.       DO 20 I=1,NSTRS
  246.       AC1(I)=0.D0
  247.       DO 20 J=1,NSTRS
  248.       AC1(I)=AC1(I)+D(I,J)*DFSIG1(J)
  249.    20 CONTINUE
  250. C
  251. C ---------------- direction de compression ----------------
  252. C
  253.       IF(IVIS.LE.2) THEN
  254.        CALL UNICOU(DK2,PAEC,2,SEQ2,BETJEF)
  255.       ELSE
  256.        CALL UNICO1(DK2,PAEC,2,SEQ2,BETJEF,NECH0,NECH1)
  257.       ENDIF
  258. C
  259.       IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  260.       SX=SIGF(1)
  261.       SY=SIGF(2)
  262.       SXY=SIGF(3)
  263.       R1=(SX-SY)*(SX-SY)+SX*SX+SY*SY+6*SXY*SXY
  264.       R=SQRT(0.5*R1)
  265.       IF (ABS(R).LE.1.D-10) THEN
  266.       WRITE(*,*)'INDETERMINATION de DFSIG ds BEHAV3'
  267.       STOP
  268.       ENDIF
  269.           DFSIG2(1)=(2.D0*SX-SY)/2.D0/R/A2+A1/A2
  270.           DFSIG2(2)=(2.D0*SY-SX)/2.D0/R/A2+A1/A2
  271.           DFSIG2(3)=3.D0*SXY/R/A2
  272.       P2(1,1)=2.D0/A2/A2
  273.       P2(1,2)=-1.D0/A2/A2
  274.       P2(1,3)=0.D0/A2/A2
  275.       P2(2,1)=-1.D0/A2/A2
  276.       P2(2,2)=2.D0/A2/A2
  277.       P2(2,3)=0.D0/A2/A2
  278.       P2(3,1)=0.D0/A2/A2
  279.       P2(3,2)=0.D0/A2/A2
  280.       P2(3,3)=6.D0/A2/A2
  281.       ENDIF
  282. C
  283.       IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  284.       SX=SIGF(1)
  285.       SY=SIGF(2)
  286.       SZ=SIGF(3)
  287.       SXY=SIGF(4)
  288.       TX=(SX-SY)*(SX-SY)
  289.       TY=(SX-SZ)*(SX-SZ)
  290.       TZ=(SY-SZ)*(SY-SZ)
  291.       R1=TX+TY+TZ+6*SXY*SXY
  292.       R=SQRT(0.5*R1)
  293.       IF (ABS(R).LE.1.D-10) THEN
  294.       WRITE(*,*)'INDETERMINATION de DFSIG ds BEHAV3'
  295.       STOP
  296.       ENDIF
  297.           DFSIG2(1)=(2.D0*SX-SY-SZ)/2.D0/R/A2+A1/A2
  298.           DFSIG2(2)=(2.D0*SY-SX-SZ)/2.D0/R/A2+A1/A2
  299.           DFSIG2(3)=(2.D0*SZ-SX-SY)/2.D0/R/A2+A1/A2
  300.           DFSIG2(4)=3.D0*SXY/R/A2
  301.       P2(1,1)=2.D0/A2/A2
  302.       P2(1,2)=-1.D0/A2/A2
  303.       P2(1,3)=-1.D0/A2/A2
  304.       P2(1,4)=0.D0/A2/A2
  305.       P2(2,1)=-1.D0/A2/A2
  306.       P2(2,2)=2.D0/A2/A2
  307.       P2(2,3)=-1.D0/A2/A2
  308.       P2(2,4)=0.D0/A2/A2
  309.       P2(3,1)=-1.D0/A2/A2
  310.       P2(3,2)=-1.D0/A2/A2
  311.       P2(3,3)=2.D0/A2/A2
  312.       P2(3,4)=0.D0/A2/A2
  313.       P2(4,1)=0.D0/A2/A2
  314.       P2(4,2)=0.D0/A2/A2
  315.       P2(4,3)=0.D0/A2/A2
  316.       P2(4,4)=6.D0/A2/A2
  317.       ENDIF
  318. C
  319.       DO 25 I=1,NSTRS
  320.       AC2(I)=0.D0
  321.       DO 25 J=1,NSTRS
  322.       AC2(I)=AC2(I)+D(I,J)*DFSIG2(J)
  323.    25 CONTINUE
  324. C
  325. C ************ Preparation du jacobien  ********************
  326. C
  327.       F1DF1=0.D0
  328.       DO 30 J=1,NSTRS
  329.       F1DF1=F1DF1+DFSIG1(J)*AC1(J)
  330.    30 CONTINUE
  331.       F1DF2=0.D0
  332.       DO 31 J=1,NSTRS
  333.       F1DF2=F1DF2+DFSIG1(J)*AC2(J)
  334.    31 CONTINUE
  335.       F2DF2=0.D0
  336.       DO 32 J=1,NSTRS
  337.       F2DF2=F2DF2+DFSIG2(J)*AC2(J)
  338.    32 CONTINUE
  339.       F2DF1=0.D0
  340.       DO 33 J=1,NSTRS
  341.       F2DF1=F2DF1+DFSIG2(J)*AC1(J)
  342.    33 CONTINUE
  343. C
  344.       DJAC0(1,1)=-(F1DF1+PAECT)
  345.       DJAC0(1,2)=-(F1DF2)
  346.       DJAC0(2,2)=-(F2DF2+PAEC)
  347.       DJAC0(2,1)=-(F2DF1)
  348. C
  349.       DO 43 I=1,2
  350.       DO 43 J=1,2
  351.       DJI(I,J)=DJAC0(I,J)
  352.    43 CONTINUE
  353. C
  354.       CALL INVMA2(DJI,2,ISING)
  355.       IF (ISING.EQ.1) THEN
  356.        WRITE(*,*)'MATRICE DJI singuliere ds BEHAV3'
  357.       ENDIF
  358. C
  359. C ************ Debut des iterations internes *******************
  360. C
  361.    40 CONTINUE
  362. C
  363. C *************** Determination de DK et DLAM ******************
  364. C
  365. C
  366. C
  367.       DLAM1(1)=DLAM0(1)-DJI(1,1)*FCRI0(1)-DJI(1,2)*FCRI0(2)
  368.       DLAM1(2)=DLAM0(2)-DJI(2,1)*FCRI0(1)-DJI(2,2)*FCRI0(2)
  369. C
  370.       IF (DLAM1(1).LE.0.D0) IEC1=1
  371.       IF (DLAM1(2).LE.0.D0) IEC2=1
  372.       IF (IEC1.EQ.1.OR.IEC2.EQ.1) THEN
  373. C       WRITE(*,*)'Dans BEHAV3, DLAMDA1 est negatif:',DLAM1(1)
  374. C       WRITE(*,*)'Dans BEHAV3, DLAMDA2 est negatif:',DLAM1(2)
  375. C       WRITE(*,*)'A l iteration :',ITER
  376.        GOTO 15
  377.       ENDIF
  378. C
  379.       DK1=DEPS1+DLAM1(1)
  380.       DK2=DEPS2+DLAM1(2)
  381. C
  382.       IF(IVIS.LE.2) THEN
  383.        CALL UNICOU(DK1,PAECT,1,SEQ1,BETJEF)
  384.        CALL UNICOU(DK2,PAEC,2,SEQ2,BETJEF)
  385.       ELSE
  386.        CALL UNICO1(DK1,PAECT,1,SEQ1,BETJEF,NECH0,NECH1)
  387.        CALL UNICO1(DK2,PAEC,2,SEQ2,BETJEF,NECH0,NECH1)
  388.       ENDIF
  389. C
  390. C ************** Determination de DPHI *********************
  391. C
  392.       IF (IMOD.EQ.1) THEN
  393.       AD1=YOUN/(1.D0-XNU)
  394.       coef=2.D0*SEQ2/(2.D0*SEQ2+DLAM1(2)*AD1)
  395.       DPHI1=SEQ1-0.5*coef*((SIGE(1)+SIGE(2))-AD1*DLAM1(1))
  396.       TO=SIGF(3)
  397.       DPHI2=SEQ2
  398.       ENDIF
  399.       IF (IMOD.EQ.2) THEN
  400.       DG=YOUN/2.D0/(1.D0+XNU)
  401.       AD1=YOUN/(1.D0+XNU)/(1.D0-2.D0*XNU)
  402.       coef1=(SEQ2+2.D0*DLAM1(2)*DG)/(SEQ2+3.D0*DLAM1(2)*DG)
  403.       coef2=(2.D0*DLAM1(2)*DG)/(SEQ2+3.D0*DLAM1(2)*DG)
  404.       CAS1=(SIGE(1)+SIGE(2)-0.5*DLAM1(1)*2.D0*AD1)
  405.       CAS2=SIGE(3)-DLAM1(1)*XNU*AD1
  406.       DPHI1=SEQ1-0.5*(coef1*CAS1+coef2*CAS2)
  407.       TO=SIGF(4)
  408.       DPHI2=SEQ2
  409.       ENDIF
  410.       IF (IMOD.EQ.3) THEN
  411.       DPHI2=SEQ2-A1/A2*(SIGF(1)+SIGF(2))
  412.       AD1=YOUN/(1.D0-XNU)
  413.       coef=2.D0*DPHI2/(2.D0*DPHI2+DLAM1(2)*AD1)
  414.       SOM=A1/A2*DLAM1(2)*2.D0*AD1
  415.       DPHI1=SEQ1-0.5*coef*((SIGE(1)+SIGE(2))-AD1*DLAM1(1)-SOM)
  416.       TO=SIGF(3)
  417.       ENDIF
  418.       IF (IMOD.EQ.4) THEN
  419.       SOM=SIGE(1)+SIGE(2)+SIGE(3)
  420.       DPHI2=A1*SOM/A2
  421.       DPHI2=DPHI2-A1*DLAM1(1)*ADD*(1.D0+XNU)/A2
  422.       DPHI2=DPHI2-A1*A1*DLAM1(2)*3.D0*ADD*(1.D0+XNU)/A2/A2
  423.       DPHI2=SEQ2-DPHI2
  424.       DG=YOUN/2.D0/(1.D0+XNU)/A2/A2
  425.       AD1=YOUN/(1.D0+XNU)/(1.D0-2.D0*XNU)
  426.       AD2=YOUN/(1.D0-2.D0*XNU)
  427.       coef1=(DPHI2+2.D0*DLAM1(2)*DG)/(DPHI2+3.D0*DLAM1(2)*DG)
  428.       coef2=(2.D0*DLAM1(2)*DG)/(DPHI2+3.D0*DLAM1(2)*DG)
  429.       CAS1=(SIGE(1)+SIGE(2)-DLAM1(1)*AD1-2.D0*A1/A2*DLAM1(2)*AD2)
  430.       CAS2=SIGE(3)-DLAM1(1)*XNU*AD1-A1/A2*DLAM1(2)*AD2
  431.       DPHI1=SEQ1-0.5*(coef1*CAS1+coef2*CAS2)
  432.       TO=SIGF(4)
  433.       ENDIF
  434. C      IF (DPHI1.LT.0.D0) THEN
  435. C       WRITE(*,*)'ATTENTION DPHI1 NEGATIF'
  436. C       WRITE(*,*)'DPHI1=',DPHI1
  437. C      ENDIF
  438. C      IF (DPHI2.LT.0.D0) THEN
  439. C       WRITE(*,*)'ATTENTION DPHI2 ds BEHAV3 NEGATIF'
  440. C       WRITE(*,*)'DPHI2=',DPHI2
  441. C      ENDIF
  442. C
  443. C--------------- Cas de l'apex -----------------------------
  444. C
  445.       IF (ABS(DPHI1).LE.10E-10.AND.
  446.      *ABS(DPHI2).GT.10E-10) THEN
  447.       IAPEX=1
  448. C      WRITE(*,*)'IAPEX ds BEHAV3 =',IAPEX
  449. C      WRITE(*,*)'Dans l element',IBB
  450. C      WRITE(*,*)'et au point d intégration',IGAU
  451.       DO 50 I=1,NSTRS
  452.       DO 50 J=1,NSTRS
  453.       AI(I,J)=0.D0
  454.    50 CONTINUE
  455.         AI(1,1)=0.5
  456.         AI(1,2)=0.5
  457.         AI(2,1)=AI(1,2)
  458.         AI(2,2)=AI(1,1)
  459.         IF (IMOD.EQ.1.OR.IMOD.EQ.3) AI(3,3)=0.D0
  460.         IF (IMOD.EQ.2.OR.IMOD.EQ.4) AI(3,3)=1.D0
  461.         GOTO 75
  462.       ENDIF
  463. C
  464. C---------- Cas du critere reduit a un point ----------------
  465. C
  466.       IF (ABS(DPHI2).LE.10E-10) THEN
  467.       IAPEX=2
  468. C      WRITE(*,*)'IAPEX ds BEHAV3 =',IAPEX
  469. C      WRITE(*,*)'Dans l element',IBB
  470. C      WRITE(*,*)'et au point d intégration',IGAU
  471.       DO 52 I=1,NSTRS
  472.       DO 52 J=1,NSTRS
  473.       AI(I,J)=0.D0
  474.    52 CONTINUE
  475.       IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  476.       DO 53 I=1,3
  477.       DO 53 J=1,3
  478.       AI(I,J)=1.D0/3.D0
  479.    53 CONTINUE
  480.       ENDIF
  481.         GOTO 75
  482.       ENDIF
  483. C
  484. C ************** Mise a jour des contraintes ***************
  485. C
  486. C ---------------- calcul de la matrice A ------------------
  487. C
  488.       DO 60 I=1,NSTRS
  489.       DO 60 J=1,NSTRS
  490.       A(I,J)=0.D0
  491.    60 CONTINUE
  492. C
  493.       DG=YOUN/2.D0/(1.D0+XNU)
  494. C
  495.       IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  496.       DW1=AD*(2.D0-XNU)
  497.       DW2=AD*(2.D0*XNU-1.D0)
  498.       A(1,1)=1.D0+(DLAM1(1)*DG)/2.D0/DPHI1
  499.            A(1,1)=A(1,1)+(DLAM1(2)*DW1)/2.D0/DPHI2/A2/A2
  500.       A(2,2)=A(1,1)
  501.       A(1,2)=-(DLAM1(1)*DG)/2.D0/DPHI1
  502.            A(1,2)=A(1,2)+(DLAM1(2)*DW2)/2.D0/DPHI2/A2/A2
  503.       A(2,1)=A(1,2)
  504.       A(3,3)=(DLAM1(1)*DG)/DPHI1
  505.       A(3,3)=A(3,3)+1.D0+(DLAM1(2)*3.D0*DG)/DPHI2/A2/A2
  506.       ENDIF
  507. C
  508.       IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  509.       A(1,1)=1.D0+(DLAM1(1)*DG)/2.D0/DPHI1
  510.            A(1,1)=A(1,1)+(DLAM1(2)*2.D0*DG)/DPHI2/A2/A2
  511.       A(2,2)=A(1,1)
  512.       A(3,3)=1.D0+(DLAM1(2)*2.D0*DG)/DPHI2/A2/A2
  513.       A(1,2)=-(DLAM1(1)*DG)/2.D0/DPHI1
  514.            A(1,2)=A(1,2)-(DLAM1(2)*DG)/DPHI2/A2/A2
  515.       A(2,1)=A(1,2)
  516.       A(1,3)=-(DLAM1(2)*DG)/DPHI2/A2/A2
  517.       A(3,1)=A(1,3)
  518.       A(3,2)=A(1,3)
  519.       A(2,3)=A(1,3)
  520.       A(4,4)=1.D0+(DLAM1(1)*DG)/DPHI1
  521.            A(4,4)=A(4,4)+(DLAM1(2)*DG*3.D0)/DPHI2/A2/A2
  522.       ENDIF
  523. C
  524. C -------------- invertion de la matrice A -----------------
  525. C
  526.       DO 70 I=1,NSTRS
  527.       DO 70 J=1,NSTRS
  528.       AI(I,J)=A(I,J)
  529.    70 CONTINUE
  530. C
  531.       CALL INVMA2(AI,NSTRS,ISING)
  532.       IF (ISING.EQ.1) THEN
  533.        WRITE(*,*)'MATRICE AI singuliere ds BEHAV3'
  534.       ENDIF
  535. C
  536. C -------------- mise a jour des contraintes  ------------
  537. C
  538.    75 CONTINUE
  539.       VEC1(1)=AD1
  540.       VEC1(2)=AD1
  541.       VEC1(3)=0.D0
  542.       VEC1(4)=0.D0
  543.       IF(IMOD.EQ.3) THEN
  544.        VEC2(1)=AD1*(1.D0+XNU)
  545.        VEC2(2)=AD1*(1.D0+XNU)
  546.        VEC2(3)=0.D0
  547.        VEC2(4)=0.D0
  548.       ENDIF
  549.       IF(IMOD.EQ.4) THEN
  550.        VEC2(1)=AD1*(1.D0+XNU)
  551.        VEC2(2)=AD1*(1.D0+XNU)
  552.        VEC2(3)=AD1*(1.D0+XNU)
  553.        VEC2(4)=0.D0
  554.       ENDIF
  555.       IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  556.        VEC2(1)=0.D0
  557.        VEC2(2)=0.D0
  558.        VEC2(3)=0.D0
  559.        VEC2(4)=0.D0
  560.       ENDIF
  561. C
  562.       DO 80 I=1,NSTRS
  563.       DEPSI(I)=SIGE(I)-0.5*DLAM1(1)*VEC1(I)
  564.       DEPSI(I)=DEPSI(I)-A1*DLAM1(2)*VEC2(I)/A2
  565.    80 CONTINUE
  566. C
  567.       DO 90 I=1,NSTRS
  568.       SIGF(I)=0.0D+00
  569.       DO 90 J=1,NSTRS
  570.       SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J)
  571.    90 CONTINUE
  572. C
  573. C ******** Verification des criteres ****************
  574. C
  575.       CALL PRINC(SIGF,V1,NSTRS)
  576.       IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  577.        CALL CRIVON(SIGF,SEQ,NSTRS,BETJEF)
  578.       ENDIF
  579.       IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN
  580.        CALL CRIDRU(SIGF,SEQ,NSTRS,BETJEF)
  581.       ENDIF
  582.       FCRI1(1)=V1(1)-SEQ1
  583.       IF (IAPEX.EQ.2) FCRI1(1)=0.D0
  584.       FCRI1(2)=SEQ-SEQ2
  585. C
  586.       IF (IBROY.EQ.0.AND.(ABS(FCRI1(1)).GE.CRIMAX.OR
  587.      *      .ABS(FCRI1(2)).GE.CRIMAX)) THEN
  588. C      WRITE(*,*)'****************************************'
  589. C      WRITE(*,*)'LE RESIDU DIVERGE AVEC BROYDEN'
  590. C      WRITE(*,*)'on passe donc a la secante'
  591. C      WRITE(*,*)'Dans l element',IBB
  592. C      WRITE(*,*)'et au point d intégration',IGAU
  593. C      WRITE(*,*)'CRIMAX=',CRIMAX
  594. C      WRITE(*,*)'****************************************'
  595.       ITER=ITR
  596.       ENDIF
  597. C
  598. C *******  Compteur sur la methode de resolution ****
  599. C
  600.       IF (IBROY.EQ.0.AND.ITER.EQ.ITR) THEN
  601.         IBROY=1
  602.         ITANG=1
  603.         ITER=1
  604.         IAPEX=0
  605.         IEC1=0
  606.         IEC2=0
  607.         GOTO 12
  608.       ENDIF
  609. C
  610. C ******* non convergence **************************
  611. C
  612.       IF ((ABS(FCRI1(1)).GT.PRB.OR.ABS(FCRI1(2)).GT.PRB)
  613.      *.AND.ITER.LT.ITR) THEN
  614.       IF (IBROY.EQ.0) THEN
  615.        CALL BROYDI(DJI,DLAM0,DLAM1,FCRI0,FCRI1)
  616.        DLAM0(1)=DLAM1(1)
  617.        DLAM0(2)=DLAM1(2)
  618.        FCRI0(1)=FCRI1(1)
  619.        FCRI0(2)=FCRI1(2)
  620.        IF (ITER.GE.(ITR-1)) THEN
  621. C        WRITE(*,*)'***********************'
  622. C        WRITE(*,*)'BROYDEN n a pas aboutit'
  623. C        WRITE(*,*)'ITER=',ITER
  624. C        WRITE(*,*)'FCRIT=',FCRI0(1)
  625. C        WRITE(*,*)'FCRIC=',FCRI0(2)
  626. C        WRITE(*,*)'Dans l element',IBB
  627. C        WRITE(*,*)'et au point d intégration',IGAU
  628. C        WRITE(*,*)'***********************'
  629.        ENDIF
  630.        ITER=ITER+1
  631.        GOTO 40
  632.       ENDIF
  633.       IF (IBROY.EQ.1.AND.ITANG.EQ.1) THEN
  634.        DLAM0(1)=DLAM1(1)
  635.        DLAM0(2)=DLAM1(2)
  636.        FCRI0(1)=FCRI1(1)
  637.        FCRI0(2)=FCRI1(2)
  638.        CALL PRINC(SIGF,V1,NSTRS)
  639.        IF (ITER.GE.(ITR-5)) THEN
  640. C        WRITE(*,*)'ITER=',ITER
  641. C        WRITE(*,*)'FCRIT=',FCRI0(1)
  642. C        WRITE(*,*)'FCRIC=',FCRI0(2)
  643.        ENDIF
  644.        ITER=ITER+1
  645.        GOTO 18
  646.       ENDIF
  647.       ENDIF
  648.       IF (ITER.GE.ITR.AND.(ABS(FCRI1(1)).GT.PRB2.OR.
  649.      *      ABS(FCRI1(2)).GT.PRB2)) THEN
  650.         WRITE(*,*)'NON CONVERGENCE INTERNE dans BEHAV3'
  651.         WRITE(*,*)'Dans l element',IBB
  652.         WRITE(*,*)'et au point d intégration',IGAU
  653.         WRITE(*,*)'FCRIT=',FCRI0(1)
  654.         WRITE(*,*)'FCRIC=',FCRI0(2)
  655.         WRITE(*,*)'IPLA=',IPLA
  656.         WRITE(*,*)'IFIS=',IFIS
  657. C        STOP
  658.       ENDIF
  659. C
  660. C **************** Fin des iterations internes *******************
  661. C
  662. C
  663. C ************* Calcul de la DEP ****************
  664. C
  665.       IF (IAPEX.EQ.0) THEN
  666.       CALL DERI2(D2FSI1,DPHI1,P1,SIGF,NSTRS,BETJEF)
  667.       CALL DERI2(D2FSI2,DPHI2,P2,SIGF,NSTRS,BETJEF)
  668.       CALL LADEP1(SIGF,DEP,PAECT,PAEC,NSTRS,IFOU,
  669.      & DFSIG1,DFSIG2,D2FSI1,D2FSI2,DLAM1,BETJEF)
  670.       ELSE IF (IAPEX.EQ.1) THEN
  671.       DO 95 I=1,NSTRS
  672.       DO 95 J=1,NSTRS
  673.          AH(I,J)=0.D0
  674.       DO 95 K=1,NSTRS
  675.          AH(I,J)=AH(I,J)+AI(I,K)*D(K,J)
  676.    95 CONTINUE
  677.        DFSIG1(1)=0.5*SQRT(2.D0)
  678.        DFSIG1(2)=0.5*SQRT(2.D0)
  679.        DFSIG1(3)=0.D0
  680.        IF (IMOD.EQ.2) DFSIG1(4)=0.D0
  681.        CALL CREDEP(AH,DFSIG1,PAEC,NSTRS,DEP,BETJEF)
  682.       ELSE IF (IAPEX.EQ.2) THEN
  683.       DO 96 I=1,NSTRS
  684.       DO 96 J=1,NSTRS
  685.          DEP(I,J)=0.D0
  686.       DO 96 K=1,NSTRS
  687.          DEP(I,J)=DEP(I,J)+AI(I,K)*D(K,J)
  688.    96 CONTINUE
  689.       ENDIF
  690.       IAPEX=0
  691. C      DO 97 I=1,NSTRS
  692. C      DO 97 J=1,NSTRS
  693. C         DEP(I,J)=D(I,J)
  694. C   97 CONTINUE
  695. C
  696. C ***********************************************
  697. C
  698.       DEPS1=DK1
  699.       DEPS2=DK2
  700.       SIG1=SEQ1
  701.       SIG2=SEQ2
  702.       IF(IVIS.LE.2) THEN
  703.        CALL INDICA(DK1,DK2,IFIS,IPLA,3,BETJEF)
  704.       ELSE
  705.        CALL INDIC1(DK1,DK2,IFIS,IPLA,3,BETJEF,NECH0,NECH1)
  706.       ENDIF
  707. C ********************** CALCUL VISCOPLASTIQUE ***********************
  708.       IF (IVIS.EQ.1) THEN
  709.         CALL VISPLA(SIGR,SIGF,DSIGT,NSTRS,DEPS1,DEPS2,
  710.      &                  SIGP,SIGRV,DSTRN,D,BETJEF,VISCO)
  711.       ENDIF
  712. C ********************************************************************
  713.   100 CONTINUE
  714.       IEC1=0
  715.       IEC2=0
  716.       CRIMAX=0.D0
  717. C
  718.       RETURN
  719.       END
  720.  
  721.  
  722.  
  723.  
  724.  
  725.  
  726.  

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