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  26/08/24    21:15:14     12622          
  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 101 I=1,NSTRS
  132.       DO 16 J=1,NSTRS
  133.       DEP(I,J)=D(I,J)
  134.    16 CONTINUE
  135.   101 CONTINUE
  136.        GOTO 100
  137.       ENDIF
  138. C
  139.     9 CONTINUE
  140. C
  141.       IF(IVIS.LE.2) THEN
  142.        CALL UNICOU(DK,PAEC,2,SEQ2,BETJEF)
  143.       ELSE
  144.        CALL UNICO1(DK,PAEC,2,SEQ2,BETJEF,NECH0,NECH1)
  145.       ENDIF
  146. C
  147. C ************ Preparation des criteres ********************
  148. C
  149. C ------------------ Von Mises -----------------------------
  150.       IF (IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  151.         A1=0.D0
  152.         A2=1.D0
  153.       ENDIF
  154. C ----------------- Drucker Prager -------------------------
  155.       IF (IMOD.EQ.3.OR.IMOD.EQ.4) THEN
  156.         A1=(BETA-1.D0)/(2*BETA-1.D0)
  157.         A2=BETA/(2*BETA-1.D0)
  158.       ENDIF
  159. C ----------------------------------------------------------
  160.       IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  161.       SX=SIGF(1)
  162.       SY=SIGF(2)
  163.       SXY=SIGF(3)
  164.         R1=(SX-SY)*(SX-SY)+SX*SX+SY*SY+6*SXY*SXY
  165.         R=SQRT(0.5*R1)
  166.       IF (ABS(R).LE.1.D-10) THEN
  167.       WRITE(*,*)'INDETERMINATION de DFSIG ds BEHAV2'
  168.       STOP
  169.       ENDIF
  170.        DFSIG(1)=(2.D0*SX-SY)/2.D0/R/A2+A1/A2
  171.        DFSIG(2)=(2.D0*SY-SX)/2.D0/R/A2+A1/A2
  172.        DFSIG(3)=3.D0*SXY/R/A2
  173.       P2(1,1)=2.D0/A2/A2
  174.       P2(1,2)=-1.D0/A2/A2
  175.       P2(1,3)=0.D0/A2/A2
  176.       P2(2,1)=-1.D0/A2/A2
  177.       P2(2,2)=2.D0/A2/A2
  178.       P2(2,3)=0.D0/A2/A2
  179.       P2(3,1)=0.D0/A2/A2
  180.       P2(3,2)=0.D0/A2/A2
  181.       P2(3,3)=6.D0/A2/A2
  182.       ENDIF
  183. C
  184.       IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  185.       SX=SIGF(1)
  186.       SY=SIGF(2)
  187.       SZ=SIGF(3)
  188.       SXY=SIGF(4)
  189.       TX=(SX-SY)*(SX-SY)
  190.       TY=(SX-SZ)*(SX-SZ)
  191.       TZ=(SY-SZ)*(SY-SZ)
  192.          R1=TX+TY+TZ+6*SXY*SXY
  193.          R=SQRT(0.5*R1)
  194.       IF (ABS(R).LE.1.D-10) THEN
  195.       WRITE(*,*)'INDETERMINATION de DFSIG ds BEHAV2'
  196.       ENDIF
  197.         DFSIG(1)=(2.D0*SX-SY-SZ)/2.D0/R/A2+A1/A2
  198.         DFSIG(2)=(2.D0*SY-SX-SZ)/2.D0/R/A2+A1/A2
  199.         DFSIG(3)=(2.D0*SZ-SX-SY)/2.D0/R/A2+A1/A2
  200.         DFSIG(4)=3.D0*SXY/R/A2
  201.       P2(1,1)=2.D0/A2/A2
  202.       P2(1,2)=-1.D0/A2/A2
  203.       P2(1,3)=-1.D0/A2/A2
  204.       P2(1,4)=0.D0/A2/A2
  205.       P2(2,1)=-1.D0/A2/A2
  206.       P2(2,2)=2.D0/A2/A2
  207.       P2(2,3)=-1.D0/A2/A2
  208.       P2(2,4)=0.D0/A2/A2
  209.       P2(3,1)=-1.D0/A2/A2
  210.       P2(3,2)=-1.D0/A2/A2
  211.       P2(3,3)=2.D0/A2/A2
  212.       P2(3,4)=0.D0/A2/A2
  213.       P2(4,1)=0.D0/A2/A2
  214.       P2(4,2)=0.D0/A2/A2
  215.       P2(4,3)=0.D0/A2/A2
  216.       P2(4,4)=6.D0/A2/A2
  217.       ENDIF
  218. C
  219.       DO 102 I=1,NSTRS
  220.       AC(I)=0.D0
  221.       DO 25 J=1,NSTRS
  222.       AC(I)=AC(I)+D(I,J)*DFSIG(J)
  223.    25 CONTINUE
  224.   102 CONTINUE
  225. C
  226. C ************ Preparation du jacobien  ***********************
  227. C
  228.       FDF=0.D0
  229.       DO 30 J=1,NSTRS
  230.       FDF=FDF+DFSIG(J)*AC(J)
  231.    30 CONTINUE
  232. C
  233.       DJAC0=-(FDF+PAEC)
  234. C
  235. C ************ Debut des iterations internes *******************
  236. C
  237.    40 CONTINUE
  238. C
  239. C *************** Determination de DK et DLAM ******************
  240. C
  241. C
  242.       DLAM1=DLAM0-FCRIC0/DJAC0
  243. C
  244.       IF (DLAM1.LT.0.D0) IEC2=1
  245.       IF (IEC2.EQ.1) THEN
  246. C       WRITE(*,*)'Dans BEHAV2, DLAMDA est negatif:',DLAM1
  247. C       WRITE(*,*)'A l iteration :',ITER
  248.         GOTO 12
  249.       ENDIF
  250. C
  251.       DK=DEPS2+DLAM1
  252. C
  253.       IF(IVIS.LE.2) THEN
  254.        CALL UNICOU(DK,PAEC,2,SEQ2,BETJEF)
  255.       ELSE
  256.        CALL UNICO1(DK,PAEC,2,SEQ2,BETJEF,NECH0,NECH1)
  257.       ENDIF
  258. C
  259. C ************** Determination de DPHI2 *********************
  260. C
  261.       IF (IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  262.        DPHI2=SEQ2
  263.       ENDIF
  264.       IF (IMOD.EQ.3) THEN
  265.        SOM=SIGF(1)+SIGF(2)+SIGF(3)
  266.        DPHI2=SEQ2-A1*SOM/A2
  267.       ENDIF
  268.       IF (IMOD.EQ.4) THEN
  269.       SOM=SIGE(1)+SIGE(2)+SIGE(3)
  270.       DPHI2=A1*SOM/A2-A1*A1*DLAM1*3.D0*ADD*(1.D0+XNU)/A2/A2
  271.       DPHI2=SEQ2-DPHI2
  272.       ENDIF
  273. C      IF (DPHI2.LT.0.D0) THEN
  274. C       WRITE(*,*)'ATTENTION DPHI2 ds BEHAV2 NEGATIF'
  275. C       WRITE(*,*)'DPHI2=',DPHI2
  276. C      ENDIF
  277. C
  278. C---------- Cas du critere reduit a un point ----------------
  279. C
  280.       IF (ABS(DPHI2).LE.10E-10) THEN
  281.       IAPEX=1
  282. C      WRITE(*,*)'IAPEX ds BEHAV2 =',IAPEX
  283.       DO 103 I=1,NSTRS
  284.       DO 50 J=1,NSTRS
  285.       AI(I,J)=0.D0
  286.    50 CONTINUE
  287.   103 CONTINUE
  288.       IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  289.       DO 104 I=1,3
  290.       DO 53 J=1,3
  291.       AI(I,J)=1.D0/3.D0
  292.    53 CONTINUE
  293.   104 CONTINUE
  294.       ENDIF
  295.         GOTO 75
  296.       ENDIF
  297. C
  298. C ************** Mise a jour des contraintes ***************
  299. C
  300. C ---------------- calcul de la matrice A ------------------
  301. C
  302.       DO 105 I=1,NSTRS
  303.       DO 60 J=1,NSTRS
  304.       A(I,J)=0.D0
  305.    60 CONTINUE
  306.   105 CONTINUE
  307. C
  308.       DG=YOUN/2.D0/(1.D0+XNU)
  309. C
  310.       IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  311.       DW1=AD*(2.D0-XNU)
  312.       DW2=AD*(2.D0*XNU-1.D0)
  313.       A(1,1)=1.D0+(DLAM1*DW1)/2.D0/DPHI2/A2/A2
  314.       A(2,2)=A(1,1)
  315.       A(1,2)=(DLAM1*DW2)/2.D0/DPHI2/A2/A2
  316.       A(2,1)=A(1,2)
  317.       A(3,3)=1.D0+(DLAM1*3.D0*DG)/DPHI2/A2/A2
  318.       ENDIF
  319.       IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  320.       A(1,1)=1.D0+(DLAM1*2.D0*DG)/DPHI2/A2/A2
  321.       A(2,2)=A(1,1)
  322.       A(3,3)=1.D0+(DLAM1*2.D0*DG)/DPHI2/A2/A2
  323.       A(1,2)=-(DLAM1*DG)/DPHI2/A2/A2
  324.       A(2,1)=A(1,2)
  325.       A(1,3)=-(DLAM1*DG)/DPHI2/A2/A2
  326.       A(3,1)=A(1,3)
  327.       A(3,2)=A(1,3)
  328.       A(2,3)=A(1,3)
  329.       A(4,4)=1.D0+(DLAM1*DG*3.D0)/DPHI2/A2/A2
  330.       ENDIF
  331. C
  332. C -------------- invertion de la matrice A -----------------
  333. C
  334.       DO 106 I=1,NSTRS
  335.       DO 70 J=1,NSTRS
  336.       AI(I,J)=A(I,J)
  337.    70 CONTINUE
  338.   106 CONTINUE
  339.       CALL INVMA2(AI,NSTRS,ISING)
  340.       IF (ISING.EQ.1) THEN
  341.        WRITE(*,*)'MATRICE AI singuliere ds BEHAV2'
  342.        WRITE(*,*)'NSTRS=',NSTRS
  343.       ENDIF
  344. C
  345. C -------------- mise a jour des contraintes  --------------
  346. C
  347.    75 CONTINUE
  348.       IF(IMOD.EQ.3) THEN
  349.        VEC1(1)=AD*(1.D0+XNU)
  350.        VEC1(2)=AD*(1.D0+XNU)
  351.        VEC1(3)=0.D0
  352.        VEC1(4)=0.D0
  353.       ENDIF
  354.       IF(IMOD.EQ.4) THEN
  355.        VEC1(1)=ADD*(1.D0+XNU)
  356.        VEC1(2)=ADD*(1.D0+XNU)
  357.        VEC1(3)=ADD*(1.D0+XNU)*2.D0*XNU
  358.        VEC1(4)=0.D0
  359.       ENDIF
  360.       IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  361.        VEC1(1)=0.D0
  362.        VEC1(2)=0.D0
  363.        VEC1(3)=0.D0
  364.        VEC1(4)=0.D0
  365.       ENDIF
  366. C
  367.       DO 80 I=1,NSTRS
  368.       DEPSI(I)=SIGE(I)-A1*DLAM1*VEC1(I)/A2
  369.    80 CONTINUE
  370. C
  371.       DO 107 I=1,NSTRS
  372.       SIGF(I)=0.0D+00
  373.       DO 90 J=1,NSTRS
  374.       SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J)
  375.    90 CONTINUE
  376.   107 CONTINUE
  377. C
  378. C ************ Verification du critere ********************
  379. C
  380.       IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  381.        CALL CRIVON(SIGF,SEQ,NSTRS,BETJEF)
  382.       ENDIF
  383.       IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN
  384.        CALL CRIDRU(SIGF,SEQ,NSTRS,BETJEF)
  385.       ENDIF
  386.       FCRIC1=SEQ-SEQ2
  387. C
  388.       IF (IBROY.EQ.0.AND.ABS(FCRIC1).GE.CRIMAX) THEN
  389. C      WRITE(*,*)'****************************************'
  390. C      WRITE(*,*)'LE RESIDU DIVERGE AVEC BROYDEN'
  391. C      WRITE(*,*)'on passe donc a la secante'
  392. C      WRITE(*,*)'Dans l element',IBB
  393. C      WRITE(*,*)'et au point d intégration',IGAU
  394. C      WRITE(*,*)'CRIMAX=',CRIMAX
  395. C      WRITE(*,*)'****************************************'
  396.       ITER=ITR
  397.       ENDIF
  398. C
  399. C *******  Compteur sur la methode de resolution ****
  400. C
  401.       IF (IBROY.EQ.0.AND.ITER.EQ.ITR) THEN
  402.         IBROY=1
  403.         ITANG=1
  404.         IAPEX=0
  405.         IEC2=0
  406.         ITER=1
  407.         GOTO 8
  408.       ENDIF
  409. C
  410. C ******* non convergence **************************
  411. C
  412.       IF (ABS(FCRIC1).GT.PRB.AND.ITER.LT.ITR) THEN
  413.        IF (IBROY.EQ.0) THEN
  414.        DJAC0=(FCRIC0-FCRIC1)/(DLAM0-DLAM1)
  415.        DLAM0=DLAM1
  416.        FCRIC0=FCRIC1
  417.        ITER=ITER+1
  418.        IF (ITER.GE.(ITR-1)) THEN
  419. C        WRITE(*,*)'***********************'
  420. C        WRITE(*,*)'BROYDEN n a pas aboutit'
  421. C        WRITE(*,*)'Dans l element',IBB
  422. C        WRITE(*,*)'et au point d intégration',IGAU
  423. C        WRITE(*,*)'ITER=',ITER
  424. C        WRITE(*,*)'FCRI=',FCRIC1
  425. C        WRITE(*,*)'***********************'
  426.        ENDIF
  427.        GOTO 40
  428.       ENDIF
  429.       IF (IBROY.EQ.1.AND.ITANG.EQ.1) THEN
  430.        DLAM0=DLAM1
  431.        FCRIC0=FCRIC1
  432.        ITER=ITER+1
  433.        IF (ITER.GE.(ITR-5)) THEN
  434. C        WRITE(*,*)'ITER=',ITER
  435. C        WRITE(*,*)'FCRI=',FCRIC1
  436.        ENDIF
  437.        GOTO 9
  438.       ENDIF
  439.       ENDIF
  440.       IF (ITER.GE.ITR.AND.ABS(FCRIC1).GT.PRB2) THEN
  441.         WRITE(*,*)'NON CONVERGENCE INTERNE dans BEHAV2'
  442.         WRITE(*,*)'FCRI=',FCRIC1
  443.         WRITE(*,*)'Dans l element',IBB
  444.         WRITE(*,*)'et au point d intégration',IGAU
  445.         WRITE(*,*)'IPLA=',IPLA
  446.         WRITE(*,*)'IFIS=',IFIS
  447. C        STOP
  448.       ENDIF
  449. C
  450. C ************* Calcul de la DEP ****************
  451. C
  452.       IF (IAPEX.EQ.0) THEN
  453.       CALL DERI2(D2FSIG,DPHI2,P2,SIGF,NSTRS,BETJEF)
  454.       CALL LADEP(SIGF,DEP,PAEC,DPHI,NSTRS,IFOU,
  455.      & DFSIG,D2FSIG,DLAM1,D,DP,BETJEF)
  456.       ELSE IF (IAPEX.EQ.1) THEN
  457. C      WRITE(*,*)'ELAS2',IAPEX
  458.       DO 109 I=1,NSTRS
  459.       DO 108 J=1,NSTRS
  460.          DEP(I,J)=0.D0
  461.       DO 96 K=1,NSTRS
  462.          DEP(I,J)=DEP(I,J)+AI(I,K)*D(K,J)
  463.    96 CONTINUE
  464.   108 CONTINUE
  465.   109 CONTINUE
  466.       ENDIF
  467.       IAPEX=0
  468.       CRIMAX=0.D0
  469. *      DO 97 I=1,NSTRS
  470. *      DO 97 J=1,NSTRS
  471. *         DEP(I,J)=D(I,J)
  472. *   97 CONTINUE
  473. C
  474. C ***********************************************
  475. C
  476.       DEPS2=DK
  477.       SIG2=SEQ2
  478.       IF(IVIS.LE.2) THEN
  479.        CALL INDICA(0.D0,DK,IFIS,IPLA,2,BETJEF)
  480.       ELSE
  481.        CALL INDIC1(0.D0,DK,IFIS,IPLA,2,BETJEF,NECH0,NECH1)
  482.       ENDIF
  483. C ********************** CALCUL VISCOPLASTIQUE ***********************
  484.       IF (IVIS.EQ.1) THEN
  485.         CALL VISPLA(SIGR,SIGF,DSIGT,NSTRS,DEPS1,DEPS2,
  486.      &     SIGP,SIGRV,DSTRN,D,BETJEF,VISCO)
  487.       ENDIF
  488. C ********************************************************************
  489.   100 CONTINUE
  490.       IEC2=0
  491. C
  492.       RETURN
  493.       END
  494.  
  495.  
  496.  
  497.  
  498.  
  499.  
  500.  
  501.  
  502.  
  503.  

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