Télécharger la plaquette Cast3M


fr en

Télécharger behav2.eso

Retour à la liste

Numérotation des lignes : 

behav2

  1. C BEHAV2    SOURCE    CB215821  16/04/21    21:15:18     8920
  2.       SUBROUTINE BEHAV2(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      Un seul critere de compression: Von Mises ou Drucker-Prager
  10. C
  11. C ==================================================================
  12. C CE SOUS-PROGRAMME EST APPELE DANS "BONE".
  13. C
  14.       IMPLICIT INTEGER(I-N)
  15.       IMPLICIT REAL*8(A-H,O-Z)
  16.       DIMENSION DFSIG(4),SIGF(4),DSIGT(4),VEC1(4),DEPSI(4),DSTRN(4)
  17.       DIMENSION SIGR(4),AC(4),SIGE(4),SIGRV(4),SIGP(4)
  18.       DIMENSION D(4,4),DEP(4,4),D2FSIG(4,4),P2(4,4),DP(4,4)
  19.       DIMENSION A(4,4),AI(4,4),AH(4,4)
  20. *
  21.       SEGMENT BETJEF
  22.        REAL*8 AA,BETA,COLI,PALF,YOUN,XNU,GFC,GFT,CAR,ETA,TDEF,
  23.      &          TCON,DPSTF1,DPSTF2,TETA,PDT,TP0
  24.        INTEGER ICT,ICC,IMOD,IVIS,ITR,
  25.      &                ISIM,IBB,IGAU,IZON
  26.       ENDSEGMENT
  27. *
  28.       SEGMENT VISCO
  29.        REAL*8 DPSTV1,DPSTV2,SIGV1,SIGV2,ENDV
  30.       ENDSEGMENT
  31.       SEGMENT NECH0
  32.        REAL*8 DT,DC,ALFG,S0,ENDO
  33.       ENDSEGMENT
  34.       SEGMENT NECH1
  35.        REAL*8 ENDL
  36.       ENDSEGMENT
  37. *
  38. *
  39. C
  40. *     COMMON /DBETJEF/AA,BETA,COLI,PALF,YOUN,XNU,GFC,GFT,CAR,ETA,TDEF,
  41. *    &              TCON,DPSTF1,DPSTF2,TETA,PDT,ICT,ICC,IMOD,IVIS,ITR,
  42. *    &                ISIM,IBB,IGAU,IZON
  43. C
  44.       IAPEX=0
  45.       PRB=1.D-10
  46.       PRB2=1.D-6
  47.       ITER=1
  48.       IBROY=0
  49.       ITANG=0
  50.       Rb=COLI
  51.       Ft=PALF*COLI
  52.       IEC2=0
  53.       CRIMAX=0.D0
  54.       SEQ = 0.D0
  55.       SEQ2 = 0.D0
  56.       CALL ZERO(SIGE,4,1)
  57. ***** CALL ZERO(V1,4,1)
  58.       CALL ZERO(P2,4,4)
  59.       CALL ZERO(D,4,4)
  60. C
  61.       DO 10 I=1,NSTRS
  62.       SIGE(I)=SIGF(I)
  63.    10 CONTINUE
  64. C
  65.       IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  66.       AD=YOUN/(1.D0-XNU*XNU)
  67.       D(1,1)=AD
  68.       D(2,2)=D(1,1)
  69.       D(3,3)=AD*(1.D0-XNU)/2.D0
  70.       D(1,2)=AD*XNU
  71.       D(2,1)=D(1,2)
  72.       ENDIF
  73. C
  74.       IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  75.       ADD=YOUN/((1.D0+XNU)*(1.D0-2.D0*XNU))
  76.       D(1,1)=ADD*(1.D0-XNU)
  77.       D(2,2)=D(1,1)
  78.       D(3,3)=D(1,1)
  79.       D(1,2)=ADD*XNU
  80.       D(2,1)=D(1,2)
  81.       D(1,3)=D(1,2)
  82.       D(2,3)=D(1,2)
  83.       D(3,1)=D(1,2)
  84.       D(3,2)=D(1,2)
  85.       D(4,4)=0.5*ADD*(1.D0-2.D0*XNU)
  86.       ENDIF
  87. C
  88. C ************ Le point est comprime pour 1ere fois ***********
  89. C
  90.       IF (IPLA.EQ.0) THEN
  91.        IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  92.         CALL CRIVON(SIGF,SEQ,NSTRS,BETJEF)
  93.        ENDIF
  94.        IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN
  95.         CALL CRIDRU(SIGF,SEQ,NSTRS,BETJEF)
  96.        ENDIF
  97.         DEPS2=0.D0
  98.         SIG2=Rb*AA
  99.         IF (SEQ.GT.SIG2) THEN
  100.          IPLA=1
  101.         ENDIF
  102.       ENDIF
  103.       IF (IFIS.EQ.0) THEN
  104.         DEPS1=0.D0
  105.         SIG1=Ft
  106.       ENDIF
  107. C
  108.     8 CONTINUE
  109. C
  110.        IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  111.         CALL CRIVON(SIGE,SEQ,NSTRS,BETJEF)
  112.        ENDIF
  113.        IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN
  114.         CALL CRIDRU(SIGE,SEQ,NSTRS,BETJEF)
  115.        ENDIF
  116.       FCRIC0=SEQ-SIG2
  117.       CRIMAX=ABS(100.D0*FCRIC0)
  118.       IF (FCRIC0.LE.0.D0) THEN
  119.         IEC2=1
  120.       ENDIF
  121.       DK=DEPS2
  122.       DLAM0=0.D0
  123. C
  124.    12 CONTINUE
  125. C
  126.       DO 15 I=1,NSTRS
  127.       SIGF(I)=SIGE(I)
  128.    15 CONTINUE
  129. C
  130.       IF (IEC2.EQ.1) THEN
  131.       DO 16 I=1,NSTRS
  132.       DO 16 J=1,NSTRS
  133.       DEP(I,J)=D(I,J)
  134.    16 CONTINUE
  135.        GOTO 100
  136.       ENDIF
  137. C
  138.     9 CONTINUE
  139. C
  140.       IF(IVIS.LE.2) THEN
  141.        CALL UNICOU(DK,PAEC,2,SEQ2,BETJEF)
  142.       ELSE
  143.        CALL UNICO1(DK,PAEC,2,SEQ2,BETJEF,NECH0,NECH1)
  144.       ENDIF
  145. C
  146. C ************ Preparation des criteres ********************
  147. C
  148. C ------------------ Von Mises -----------------------------
  149.       IF (IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  150.         A1=0.D0
  151.         A2=1.D0
  152.       ENDIF
  153. C ----------------- Drucker Prager -------------------------
  154.       IF (IMOD.EQ.3.OR.IMOD.EQ.4) THEN
  155.         A1=(BETA-1.D0)/(2*BETA-1.D0)
  156.         A2=BETA/(2*BETA-1.D0)
  157.       ENDIF
  158. C ----------------------------------------------------------
  159.       IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  160.       SX=SIGF(1)
  161.       SY=SIGF(2)
  162.       SXY=SIGF(3)
  163.         R1=(SX-SY)*(SX-SY)+SX*SX+SY*SY+6*SXY*SXY
  164.         R=SQRT(0.5*R1)
  165.       IF (ABS(R).LE.1.D-10) THEN
  166.       WRITE(*,*)'INDETERMINATION de DFSIG ds BEHAV2'
  167.       STOP
  168.       ENDIF
  169.        DFSIG(1)=(2.D0*SX-SY)/2.D0/R/A2+A1/A2
  170.        DFSIG(2)=(2.D0*SY-SX)/2.D0/R/A2+A1/A2
  171.        DFSIG(3)=3.D0*SXY/R/A2
  172.       P2(1,1)=2.D0/A2/A2
  173.       P2(1,2)=-1.D0/A2/A2
  174.       P2(1,3)=0.D0/A2/A2
  175.       P2(2,1)=-1.D0/A2/A2
  176.       P2(2,2)=2.D0/A2/A2
  177.       P2(2,3)=0.D0/A2/A2
  178.       P2(3,1)=0.D0/A2/A2
  179.       P2(3,2)=0.D0/A2/A2
  180.       P2(3,3)=6.D0/A2/A2
  181.       ENDIF
  182. C
  183.       IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  184.       SX=SIGF(1)
  185.       SY=SIGF(2)
  186.       SZ=SIGF(3)
  187.       SXY=SIGF(4)
  188.       TX=(SX-SY)*(SX-SY)
  189.       TY=(SX-SZ)*(SX-SZ)
  190.       TZ=(SY-SZ)*(SY-SZ)
  191.          R1=TX+TY+TZ+6*SXY*SXY
  192.          R=SQRT(0.5*R1)
  193.       IF (ABS(R).LE.1.D-10) THEN
  194.       WRITE(*,*)'INDETERMINATION de DFSIG ds BEHAV2'
  195.       ENDIF
  196.         DFSIG(1)=(2.D0*SX-SY-SZ)/2.D0/R/A2+A1/A2
  197.         DFSIG(2)=(2.D0*SY-SX-SZ)/2.D0/R/A2+A1/A2
  198.         DFSIG(3)=(2.D0*SZ-SX-SY)/2.D0/R/A2+A1/A2
  199.         DFSIG(4)=3.D0*SXY/R/A2
  200.       P2(1,1)=2.D0/A2/A2
  201.       P2(1,2)=-1.D0/A2/A2
  202.       P2(1,3)=-1.D0/A2/A2
  203.       P2(1,4)=0.D0/A2/A2
  204.       P2(2,1)=-1.D0/A2/A2
  205.       P2(2,2)=2.D0/A2/A2
  206.       P2(2,3)=-1.D0/A2/A2
  207.       P2(2,4)=0.D0/A2/A2
  208.       P2(3,1)=-1.D0/A2/A2
  209.       P2(3,2)=-1.D0/A2/A2
  210.       P2(3,3)=2.D0/A2/A2
  211.       P2(3,4)=0.D0/A2/A2
  212.       P2(4,1)=0.D0/A2/A2
  213.       P2(4,2)=0.D0/A2/A2
  214.       P2(4,3)=0.D0/A2/A2
  215.       P2(4,4)=6.D0/A2/A2
  216.       ENDIF
  217. C
  218.       DO 25 I=1,NSTRS
  219.       AC(I)=0.D0
  220.       DO 25 J=1,NSTRS
  221.       AC(I)=AC(I)+D(I,J)*DFSIG(J)
  222.    25 CONTINUE
  223. C
  224. C ************ Preparation du jacobien  ***********************
  225. C
  226.       FDF=0.D0
  227.       DO 30 J=1,NSTRS
  228.       FDF=FDF+DFSIG(J)*AC(J)
  229.    30 CONTINUE
  230. C
  231.       DJAC0=-(FDF+PAEC)
  232. C
  233. C ************ Debut des iterations internes *******************
  234. C
  235.    40 CONTINUE
  236. C
  237. C *************** Determination de DK et DLAM ******************
  238. C
  239. C
  240.       DLAM1=DLAM0-FCRIC0/DJAC0
  241. C
  242.       IF (DLAM1.LT.0.D0) IEC2=1
  243.       IF (IEC2.EQ.1) THEN
  244. C       WRITE(*,*)'Dans BEHAV2, DLAMDA est negatif:',DLAM1
  245. C       WRITE(*,*)'A l iteration :',ITER
  246.         GOTO 12
  247.       ENDIF
  248. C
  249.       DK=DEPS2+DLAM1
  250. C
  251.       IF(IVIS.LE.2) THEN
  252.        CALL UNICOU(DK,PAEC,2,SEQ2,BETJEF)
  253.       ELSE
  254.        CALL UNICO1(DK,PAEC,2,SEQ2,BETJEF,NECH0,NECH1)
  255.       ENDIF
  256. C
  257. C ************** Determination de DPHI2 *********************
  258. C
  259.       IF (IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  260.        DPHI2=SEQ2
  261.       ENDIF
  262.       IF (IMOD.EQ.3) THEN
  263.        SOM=SIGF(1)+SIGF(2)+SIGF(3)
  264.        DPHI2=SEQ2-A1*SOM/A2
  265.       ENDIF
  266.       IF (IMOD.EQ.4) THEN
  267.       SOM=SIGE(1)+SIGE(2)+SIGE(3)
  268.       DPHI2=A1*SOM/A2-A1*A1*DLAM1*3.D0*ADD*(1.D0+XNU)/A2/A2
  269.       DPHI2=SEQ2-DPHI2
  270.       ENDIF
  271. C      IF (DPHI2.LT.0.D0) THEN
  272. C       WRITE(*,*)'ATTENTION DPHI2 ds BEHAV2 NEGATIF'
  273. C       WRITE(*,*)'DPHI2=',DPHI2
  274. C      ENDIF
  275. C
  276. C---------- Cas du critere reduit a un point ----------------
  277. C
  278.       IF (ABS(DPHI2).LE.10E-10) THEN
  279.       IAPEX=1
  280. C      WRITE(*,*)'IAPEX ds BEHAV2 =',IAPEX
  281.       DO 50 I=1,NSTRS
  282.       DO 50 J=1,NSTRS
  283.       AI(I,J)=0.D0
  284.    50 CONTINUE
  285.       IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  286.       DO 53 I=1,3
  287.       DO 53 J=1,3
  288.       AI(I,J)=1.D0/3.D0
  289.    53 CONTINUE
  290.       ENDIF
  291.         GOTO 75
  292.       ENDIF
  293. C
  294. C ************** Mise a jour des contraintes ***************
  295. C
  296. C ---------------- calcul de la matrice A ------------------
  297. C
  298.       DO 60 I=1,NSTRS
  299.       DO 60 J=1,NSTRS
  300.       A(I,J)=0.D0
  301.    60 CONTINUE
  302. C
  303.       DG=YOUN/2.D0/(1.D0+XNU)
  304. C
  305.       IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  306.       DW1=AD*(2.D0-XNU)
  307.       DW2=AD*(2.D0*XNU-1.D0)
  308.       A(1,1)=1.D0+(DLAM1*DW1)/2.D0/DPHI2/A2/A2
  309.       A(2,2)=A(1,1)
  310.       A(1,2)=(DLAM1*DW2)/2.D0/DPHI2/A2/A2
  311.       A(2,1)=A(1,2)
  312.       A(3,3)=1.D0+(DLAM1*3.D0*DG)/DPHI2/A2/A2
  313.       ENDIF
  314.       IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  315.       A(1,1)=1.D0+(DLAM1*2.D0*DG)/DPHI2/A2/A2
  316.       A(2,2)=A(1,1)
  317.       A(3,3)=1.D0+(DLAM1*2.D0*DG)/DPHI2/A2/A2
  318.       A(1,2)=-(DLAM1*DG)/DPHI2/A2/A2
  319.       A(2,1)=A(1,2)
  320.       A(1,3)=-(DLAM1*DG)/DPHI2/A2/A2
  321.       A(3,1)=A(1,3)
  322.       A(3,2)=A(1,3)
  323.       A(2,3)=A(1,3)
  324.       A(4,4)=1.D0+(DLAM1*DG*3.D0)/DPHI2/A2/A2
  325.       ENDIF
  326. C
  327. C -------------- invertion de la matrice A -----------------
  328. C
  329.       DO 70 I=1,NSTRS
  330.       DO 70 J=1,NSTRS
  331.       AI(I,J)=A(I,J)
  332.    70 CONTINUE
  333.       CALL INVMA2(AI,NSTRS,ISING)
  334.       IF (ISING.EQ.1) THEN
  335.        WRITE(*,*)'MATRICE AI singuliere ds BEHAV2'
  336.        WRITE(*,*)'NSTRS=',NSTRS
  337.       ENDIF
  338. C
  339. C -------------- mise a jour des contraintes  --------------
  340. C
  341.    75 CONTINUE
  342.       IF(IMOD.EQ.3) THEN
  343.        VEC1(1)=AD*(1.D0+XNU)
  344.        VEC1(2)=AD*(1.D0+XNU)
  345.        VEC1(3)=0.D0
  346.        VEC1(4)=0.D0
  347.       ENDIF
  348.       IF(IMOD.EQ.4) THEN
  349.        VEC1(1)=ADD*(1.D0+XNU)
  350.        VEC1(2)=ADD*(1.D0+XNU)
  351.        VEC1(3)=ADD*(1.D0+XNU)*2.D0*XNU
  352.        VEC1(4)=0.D0
  353.       ENDIF
  354.       IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  355.        VEC1(1)=0.D0
  356.        VEC1(2)=0.D0
  357.        VEC1(3)=0.D0
  358.        VEC1(4)=0.D0
  359.       ENDIF
  360. C
  361.       DO 80 I=1,NSTRS
  362.       DEPSI(I)=SIGE(I)-A1*DLAM1*VEC1(I)/A2
  363.    80 CONTINUE
  364. C
  365.       DO 90 I=1,NSTRS
  366.       SIGF(I)=0.0D+00
  367.       DO 90 J=1,NSTRS
  368.       SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J)
  369.    90 CONTINUE
  370. C
  371. C ************ Verification du critere ********************
  372. C
  373.       IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  374.        CALL CRIVON(SIGF,SEQ,NSTRS,BETJEF)
  375.       ENDIF
  376.       IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN
  377.        CALL CRIDRU(SIGF,SEQ,NSTRS,BETJEF)
  378.       ENDIF
  379.       FCRIC1=SEQ-SEQ2
  380. C
  381.       IF (IBROY.EQ.0.AND.ABS(FCRIC1).GE.CRIMAX) THEN
  382. C      WRITE(*,*)'****************************************'
  383. C      WRITE(*,*)'LE RESIDU DIVERGE AVEC BROYDEN'
  384. C      WRITE(*,*)'on passe donc a la secante'
  385. C      WRITE(*,*)'Dans l element',IBB
  386. C      WRITE(*,*)'et au point d intégration',IGAU
  387. C      WRITE(*,*)'CRIMAX=',CRIMAX
  388. C      WRITE(*,*)'****************************************'
  389.       ITER=ITR
  390.       ENDIF
  391. C
  392. C *******  Compteur sur la methode de resolution ****
  393. C
  394.       IF (IBROY.EQ.0.AND.ITER.EQ.ITR) THEN
  395.         IBROY=1
  396.         ITANG=1
  397.         IAPEX=0
  398.         IEC2=0
  399.         ITER=1
  400.         GOTO 8
  401.       ENDIF
  402. C
  403. C ******* non convergence **************************
  404. C
  405.       IF (ABS(FCRIC1).GT.PRB.AND.ITER.LT.ITR) THEN
  406.        IF (IBROY.EQ.0) THEN
  407.        DJAC0=(FCRIC0-FCRIC1)/(DLAM0-DLAM1)
  408.        DLAM0=DLAM1
  409.        FCRIC0=FCRIC1
  410.        ITER=ITER+1
  411.        IF (ITER.GE.(ITR-1)) THEN
  412. C        WRITE(*,*)'***********************'
  413. C        WRITE(*,*)'BROYDEN n a pas aboutit'
  414. C        WRITE(*,*)'Dans l element',IBB
  415. C        WRITE(*,*)'et au point d intégration',IGAU
  416. C        WRITE(*,*)'ITER=',ITER
  417. C        WRITE(*,*)'FCRI=',FCRIC1
  418. C        WRITE(*,*)'***********************'
  419.        ENDIF
  420.        GOTO 40
  421.       ENDIF
  422.       IF (IBROY.EQ.1.AND.ITANG.EQ.1) THEN
  423.        DLAM0=DLAM1
  424.        FCRIC0=FCRIC1
  425.        ITER=ITER+1
  426.        IF (ITER.GE.(ITR-5)) THEN
  427. C        WRITE(*,*)'ITER=',ITER
  428. C        WRITE(*,*)'FCRI=',FCRIC1
  429.        ENDIF
  430.        GOTO 9
  431.       ENDIF
  432.       ENDIF
  433.       IF (ITER.GE.ITR.AND.ABS(FCRIC1).GT.PRB2) THEN
  434.         WRITE(*,*)'NON CONVERGENCE INTERNE dans BEHAV2'
  435.         WRITE(*,*)'FCRI=',FCRIC1
  436.         WRITE(*,*)'Dans l element',IBB
  437.         WRITE(*,*)'et au point d intégration',IGAU
  438.         WRITE(*,*)'IPLA=',IPLA
  439.         WRITE(*,*)'IFIS=',IFIS
  440. C        STOP
  441.       ENDIF
  442. C
  443. C ************* Calcul de la DEP ****************
  444. C
  445.       IF (IAPEX.EQ.0) THEN
  446.       CALL DERI2(D2FSIG,DPHI2,P2,SIGF,NSTRS,BETJEF)
  447.       CALL LADEP(SIGF,DEP,PAEC,DPHI,NSTRS,IFOU,
  448.      & DFSIG,D2FSIG,DLAM1,D,DP,BETJEF)
  449.       ELSE IF (IAPEX.EQ.1) THEN
  450. C      WRITE(*,*)'ELAS2',IAPEX
  451.       DO 96 I=1,NSTRS
  452.       DO 96 J=1,NSTRS
  453.          DEP(I,J)=0.D0
  454.       DO 96 K=1,NSTRS
  455.          DEP(I,J)=DEP(I,J)+AI(I,K)*D(K,J)
  456.    96 CONTINUE
  457.       ENDIF
  458.       IAPEX=0
  459.       CRIMAX=0.D0
  460. *      DO 97 I=1,NSTRS
  461. *      DO 97 J=1,NSTRS
  462. *         DEP(I,J)=D(I,J)
  463. *   97 CONTINUE
  464. C
  465. C ***********************************************
  466. C
  467.       DEPS2=DK
  468.       SIG2=SEQ2
  469.       IF(IVIS.LE.2) THEN
  470.        CALL INDICA(0.D0,DK,IFIS,IPLA,2,BETJEF)
  471.       ELSE
  472.        CALL INDIC1(0.D0,DK,IFIS,IPLA,2,BETJEF,NECH0,NECH1)
  473.       ENDIF
  474. C ********************** CALCUL VISCOPLASTIQUE ***********************
  475.       IF (IVIS.EQ.1) THEN
  476.         CALL VISPLA(SIGR,SIGF,DSIGT,NSTRS,DEPS1,DEPS2,
  477.      &     SIGP,SIGRV,DSTRN,D,BETJEF,VISCO)
  478.       ENDIF
  479. C ********************************************************************
  480.   100 CONTINUE
  481.       IEC2=0
  482. C
  483.       RETURN
  484.       END
  485.  
  486.  
  487.  
  488.  
  489.  
  490.  
  491.  
  492.  
  493.  

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