Télécharger la plaquette Cast3M


fr en

Télécharger behav1.eso

Retour à la liste

Numérotation des lignes : 

behav1

  1. C BEHAV1    SOURCE    CB215821  16/04/21    21:15:17     8920
  2.       SUBROUTINE BEHAV1(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 POST PIC
  9. C                Un seul critere de traction: Rankine
  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),SIGRV(4)
  17.       DIMENSION V1(4),AC(4),SIGE(4),SIGP(4),SIGR(4),DSTRN(4)
  18.       DIMENSION D2FSIG(4,4),DEP(4,4),P1(4,4),D(4,4),DP(4,4),A(4,4)
  19.       DIMENSION 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. C
  38. *     COMMON /DBETJEF/AA,BETA,COLI,PALF,YOUN,XNU,GFC,GFT,CAR,ETA,TDEF,
  39. *    &              TCON,DPSTF1,DPSTF2,TETA,PDT,ICT,ICC,IMOD,IVIS,ITR,
  40. *    &                ISIM,IBB,IGAU,IZON
  41. C
  42. C
  43.       IAPEX=0
  44.       PRB=1.D-10
  45.       PRB2=1.D-6
  46.       ITER=1
  47.       ITANG=0
  48.       IBROY=0
  49.       Ft=PALF*COLI
  50.       CRIMAX=0.D0
  51.       SEQ = 0.D0
  52.       SEQ1 = 0.D0
  53.       CALL ZERO(SIGE,4,1)
  54.       CALL ZERO(V1,4,1)
  55.       CALL ZERO(P1,4,4)
  56.       CALL ZERO(D,4,4)
  57. C
  58.       DO 10 I=1,NSTRS
  59.       SIGE(I)=SIGF(I)
  60.    10 CONTINUE
  61. C
  62. C      DO 11 I=1,NSTRS
  63. C      DO 11 J=1,NSTRS
  64. C      D(I,J)=0.D0
  65. C   11 CONTINUE
  66. C
  67.       IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  68.       AD=YOUN/(1.D0-XNU*XNU)
  69.       D(1,1)=AD
  70.       D(2,2)=D(1,1)
  71.       D(3,3)=AD*(1.D0-XNU)/2.D0
  72.       D(1,2)=AD*XNU
  73.       D(2,1)=D(1,2)
  74.       ENDIF
  75. C
  76.       IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  77.       ADD=YOUN/((1.D0+XNU)*(1.D0-2.D0*XNU))
  78.       D(1,1)=ADD*(1.D0-XNU)
  79.       D(2,2)=D(1,1)
  80.       D(3,3)=D(1,1)
  81.       D(1,2)=ADD*XNU
  82.       D(2,1)=D(1,2)
  83.       D(1,3)=D(1,2)
  84.       D(2,3)=D(1,2)
  85.       D(3,1)=D(1,2)
  86.       D(3,2)=D(1,2)
  87.       D(4,4)=0.5*ADD*(1.D0-2.D0*XNU)
  88.       ENDIF
  89. C
  90. C ************ Le point est fissure pour 1ere fois ***********
  91. C
  92.       IF (IFIS.EQ.0) THEN
  93.         CALL PRINC(SIGF,V1,NSTRS)
  94.         DEPS1=0.D0
  95.         SIG1=Ft
  96.         IF (V1(1).GT.Ft) THEN
  97.          IFIS=1
  98.         ENDIF
  99.       ENDIF
  100.       IF (IPLA.EQ.0) THEN
  101.           DEPS2=0.D0
  102.           SIG2=COLI*AA
  103.       ENDIF
  104. C
  105.     8 CONTINUE
  106. C
  107.       DO 4 I=1,NSTRS
  108.        SIGF(I)=SIGE(I)
  109.     4 CONTINUE
  110. C
  111.       CALL PRINC(SIGF,V1,NSTRS)
  112.       TETA=V1(4)
  113.       DK=DEPS1
  114.       DLAM0=0.D0
  115.       FCRI0=V1(1)-SIG1
  116.       CRIMAX=ABS(100.D0*FCRI0)
  117. C
  118.       IF (FCRI0.LT.0.D0) THEN
  119.       DO 45 I=1,NSTRS
  120.       DO 45 J=1,NSTRS
  121.       DEP(I,J)=D(I,J)
  122.    45 CONTINUE
  123.         GOTO 100
  124.       ENDIF
  125. C
  126. C ************ Traitement du point deja fissure  ********************
  127. C
  128.     9 CONTINUE
  129.       PI=4.D0*ATAN(1.D0)
  130.       PHIC=V1(4)*(PI/180.D0)
  131.       COSA=COS(PHIC)
  132.       SINA=SIN(PHIC)
  133. C-------------------------------------------------------------------
  134.       IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  135.           DFSIG(1)=COSA*COSA
  136.           DFSIG(2)=SINA*SINA
  137.           DFSIG(3)=2.D0*SINA*COSA
  138.       DO 12 I=1,NSTRS
  139.       DO 12 J=1,NSTRS
  140.        P1(I,J)=0.D0
  141.    12 CONTINUE
  142.       P1(1,1)=1.D0/2.D0
  143.       P1(1,2)=-1.D0/2.D0
  144.       P1(2,1)=-1.D0/2.D0
  145.       P1(2,2)=1.D0/2.D0
  146.       P1(3,3)=2.D0
  147.       ENDIF
  148. C-------------------------------------------------------------------
  149.       IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  150.           DFSIG(1)=COSA*COSA
  151.           DFSIG(2)=SINA*SINA
  152.           DFSIG(3)=0.D0
  153.           DFSIG(4)=2.D0*SINA*COSA
  154.       DO 13 I=1,NSTRS
  155.       DO 13 J=1,NSTRS
  156.        P1(I,J)=0.D0
  157.    13 CONTINUE
  158.       P1(1,1)=1.D0/2.D0
  159.       P1(1,2)=-1.D0/2.D0
  160.       P1(2,1)=-1.D0/2.D0
  161.       P1(2,2)=1.D0/2.D0
  162.       P1(4,4)=2.D0
  163.       ENDIF
  164. C-------------------------------------------------------------------
  165.       DO 20 I=1,NSTRS
  166.       AC(I)=0.D0
  167.       DO 20 J=1,NSTRS
  168.       AC(I)=AC(I)+D(I,J)*DFSIG(I)
  169.    20 CONTINUE
  170. C
  171.       FDF=0.D0
  172.       DO 30 J=1,NSTRS
  173.       FDF=FDF+DFSIG(J)*AC(J)
  174.    30 CONTINUE
  175. C--------------- Determination du parametre d'ecrouissage ----------
  176. C
  177.       IF(IVIS.LE.2) THEN
  178.        CALL UNICOU(DK,PAEC,1,SEQ,BETJEF)
  179.       ELSE
  180.        CALL UNICO1(DK,PAEC,1,SEQ,BETJEF,NECH0,NECH1)
  181.       ENDIF
  182. C
  183. C-------------------------------------------------------------------
  184.       DJAC0=-(PAEC+FDF)
  185. C
  186. C ************ Debut des iterations internes ***************
  187. C
  188.    40 CONTINUE
  189. C
  190. C *************** Determination de DK  *********************
  191. C
  192.       DLAM1=-FCRI0/DJAC0+DLAM0
  193.       DK=DEPS1+DLAM1
  194. C      IF (DLAM1.LT.0.D0) THEN
  195. C       WRITE(*,*)'Dans BEHAV1, DLAMDA est negatif:',DLAM1
  196. C       WRITE(*,*)'A l iteration :',ITER
  197. C      ENDIF
  198. C
  199. C--------------- Estimation contrainte quivalente ----------
  200. C
  201.       IF(IVIS.LE.2) THEN
  202.        CALL UNICOU(DK,PAEC,1,SEQ1,BETJEF)
  203.       ELSE
  204.        CALL UNICO1(DK,PAEC,1,SEQ1,BETJEF,NECH0,NECH1)
  205.       ENDIF
  206. C
  207. C ************** Determination de DPHI *********************
  208. C
  209.       IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  210.       AD1=YOUN/(1.D0-XNU)
  211.       TO=SIGF(3)
  212.       ENDIF
  213.       IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  214.       AD1=YOUN/((1.D0+XNU)*(1.D0-2.D0*XNU))
  215.       TO=SIGF(4)
  216.       ENDIF
  217. C
  218.       DPHI=SEQ1-0.5*(SIGE(1)+SIGE(2))+0.5*AD1*DLAM1
  219. C      IF (ITER.EQ.1) THEN
  220. C      DPHI=0.25*(SIGF(1)-SIGF(2))*(SIGF(1)-SIGF(2))
  221. C      DPHI=DPHI+TO*TO
  222. C      DPHI=SQRT(DPHI)
  223. C      ENDIF
  224.       IF (DPHI.LT.0.D0) THEN
  225. C       WRITE(*,*)'ATTENTION DPHI NEGATIF'
  226. C       WRITE(*,*)'DPHI=',DPHI
  227.       ENDIF
  228. C
  229. C--------------- Cas de l'apex -----------------------------
  230. C
  231.       IF (ABS(DPHI).LE.10E-10) THEN
  232.       IAPEX=1
  233. C      WRITE(*,*)'IAPEX ds BEHAV1 =',IAPEX
  234. C      WRITE(*,*)'Dans l element',IBB
  235. C      WRITE(*,*)'et au point d intégration',IGAU
  236.       DO 50 I=1,NSTRS
  237.       DO 50 J=1,NSTRS
  238.       AI(I,J)=0.D0
  239.    50 CONTINUE
  240.         AI(1,1)=0.5
  241.         AI(1,2)=0.5
  242.         AI(2,1)=AI(1,2)
  243.         AI(2,2)=AI(1,1)
  244.         IF (IMOD.EQ.1.OR.IMOD.EQ.3) AI(3,3)=0.D0
  245.         IF (IMOD.EQ.2.OR.IMOD.EQ.4) AI(3,3)=1.D0
  246.         GOTO 75
  247.       ENDIF
  248. C
  249. C ************** Mise a jour des contraintes ***************
  250. C
  251. C ---------------- calcul de la matrice A ------------------
  252. C
  253.       DO 60 I=1,NSTRS
  254.       DO 60 J=1,NSTRS
  255.       A(I,J)=0.D0
  256.    60 CONTINUE
  257. C
  258.       DG=YOUN/2.D0/(1.D0+XNU)
  259. C
  260.       IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  261.       A(1,1)=1.D0+(DLAM1*DG)/2.D0/DPHI
  262.       A(2,2)=A(1,1)
  263.       A(1,2)=-(A(1,1)-1.D0)
  264.       A(2,1)=A(1,2)
  265.       A(3,3)=1.D0+(DLAM1*DG)/DPHI
  266.       ELSE
  267.       A(1,1)=1.D0+(DLAM1*DG)/2.D0/DPHI
  268.       A(2,2)=A(1,1)
  269.       A(1,2)=-(A(1,1)-1.D0)
  270.       A(2,1)=A(1,2)
  271.       A(3,3)=1.D0
  272.       A(4,4)=1.D0+(DLAM1*DG)/DPHI
  273.       ENDIF
  274. C
  275. C -------------- invertion de la matrice A -----------------
  276. C
  277.       DO 70 I=1,NSTRS
  278.       DO 70 J=1,NSTRS
  279.       AI(I,J)=A(I,J)
  280.    70 CONTINUE
  281.       CALL INVMA2(AI,NSTRS,ISING)
  282.       IF (ISING.EQ.1) THEN
  283.        WRITE(*,*)'MATRICE AI singuliere ds BEHAV1'
  284.       ENDIF
  285. C
  286. C -------------- mise a jour des contraintes  ------------
  287. C
  288.    75 CONTINUE
  289.       IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  290.       VEC1(1)=AD1
  291.       VEC1(2)=AD1
  292.       VEC1(3)=0.D0
  293.       VEC1(4)=0.D0
  294.       ENDIF
  295.       IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  296.       VEC1(1)=AD1
  297.       VEC1(2)=AD1
  298.       VEC1(3)=AD1*2.D0*XNU
  299.       VEC1(4)=0.D0
  300.       ENDIF
  301. C
  302.       DO 80 I=1,NSTRS
  303.       DEPSI(I)=SIGE(I)-0.5*DLAM1*VEC1(I)
  304.    80 CONTINUE
  305. C
  306.       DO 90 I=1,NSTRS
  307.       SIGF(I)=0.0D+00
  308.       DO 90 J=1,NSTRS
  309.       SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J)
  310.    90 CONTINUE
  311. C
  312. C ******** Verification du critere ****************
  313. C
  314.       CALL PRINC(SIGF,V1,NSTRS)
  315.       FCRI1=V1(1)-SEQ1
  316. C
  317.       IF (IBROY.EQ.0.AND.ABS(FCRI1).GE.CRIMAX) THEN
  318. C      WRITE(*,*)'****************************************'
  319. C      WRITE(*,*)'LE RESIDU DIVERGE AVEC BROYDEN'
  320. C      WRITE(*,*)'on passe donc a la secante'
  321. C      WRITE(*,*)'Dans l element',IBB
  322. C      WRITE(*,*)'et au point d intégration',IGAU
  323. C      WRITE(*,*)'CRIMAX=',CRIMAX
  324. C      WRITE(*,*)'****************************************'
  325.       ITER=ITR
  326.       ENDIF
  327. C
  328. C *******  Compteur sur la methode de resolution ****
  329. C
  330.       IF (IBROY.EQ.0.AND.ITER.EQ.ITR) THEN
  331.         IBROY=1
  332.         ITANG=1
  333.         ITER=1
  334.         IAPEX=0
  335.         GOTO 8
  336.       ENDIF
  337. C
  338. C ******* non convergence **************************
  339. C
  340.       IF (ABS(FCRI1).GT.PRB.AND.ITER.LT.ITR) THEN
  341.        IF (IBROY.EQ.0) THEN
  342.        DJAC1=(FCRI0-FCRI1)/(DLAM0-DLAM1)
  343.        DLAM0=DLAM1
  344.        DJAC0=DJAC1
  345.        FCRI0=FCRI1
  346.        ITER=ITER+1
  347.        IF (ITER.GE.(ITR-1)) THEN
  348. C        WRITE(*,*)'***********************'
  349. C        WRITE(*,*)'BROYDEN n a pas aboutit'
  350. C        WRITE(*,*)'Dans l element',IBB
  351. C        WRITE(*,*)'et au point d intégration',IGAU
  352. C        WRITE(*,*)'ITER=',ITER
  353. C        WRITE(*,*)'FCRI=',FCRI1
  354. C        WRITE(*,*)'***********************'
  355.        ENDIF
  356.        GOTO 40
  357.       ENDIF
  358.       IF (IBROY.EQ.1.AND.ITANG.EQ.1) THEN
  359.        DLAM0=DLAM1
  360.        FCRI0=FCRI1
  361.        CALL PRINC(SIGF,V1,NSTRS)
  362.        ITER=ITER+1
  363.        IF (ITER.GE.(ITR-5)) THEN
  364.         WRITE(*,*)'ITER=',ITER
  365.         WRITE(*,*)'FCRI=',FCRI1
  366.        ENDIF
  367.        GOTO 9
  368.       ENDIF
  369.       ENDIF
  370.       IF (ITER.GE.ITR.AND.ABS(FCRI1).GT.PRB2) THEN
  371.         WRITE(*,*)'NON CONVERGENCE INTERNE dans BEHAV1'
  372.         WRITE(*,*)'FCRI=',FCRI1
  373.         WRITE(*,*)'Dans l element',IBB
  374.         WRITE(*,*)'et au point d intégration',IGAU
  375.         WRITE(*,*)'IPLA=',IPLA
  376.         WRITE(*,*)'IFIS=',IFIS
  377. C        STOP
  378.       ENDIF
  379. C
  380. C ************* Calcul de la DEP ****************
  381. C
  382.       IF (IAPEX.EQ.0) THEN
  383.       CALL DERI2(D2FSIG,DPHI,P1,SIGF,NSTRS,BETJEF)
  384.       CALL LADEP(SIGF,DEP,PAEC,DPHI,NSTRS,IFOU,
  385.      & DFSIG,D2FSIG,DLAM1,D,DP,BETJEF)
  386.       ELSE
  387. C      WRITE(*,*)'ELAS',IAPEX
  388.       DO 95 I=1,NSTRS
  389.       DO 95 J=1,NSTRS
  390.          AH(I,J)=0.D0
  391.       DO 95 K=1,NSTRS
  392.          AH(I,J)=AH(I,J)+AI(I,K)*D(K,J)
  393.    95 CONTINUE
  394.        DFSIG(1)=0.5*SQRT(2.D0)
  395.        DFSIG(2)=0.5*SQRT(2.D0)
  396.        DFSIG(3)=0.D0
  397.        DFSIG(4)=0.D0
  398.        CALL CREDEP(AH,DFSIG,PAEC,NSTRS,DEP,BETJEF)
  399.       ENDIF
  400.       IAPEX=0
  401.       CRIMAX=0.D0
  402. C
  403. C ***********************************************
  404. C
  405.       IF (DK.LT.0.D0) THEN
  406.        WRITE(*,*)'Dans BEHAV1, DLAMDA est negatif:',DLAM1
  407.        WRITE(*,*)'A l iteration :',ITER
  408. C       STOP
  409.       ENDIF
  410.       DEPS1=DK
  411.       SIG1=V1(1)
  412.       IF (V1(1).LT.0.D0) V1(1)=0.D0
  413.       CALL INDICA(DK,0.0D0,IFIS,IPLA,1,BETJEF)
  414.       IF(IVIS.LE.2) THEN
  415.        CALL INDICA(DK,0.0D0,IFIS,IPLA,1,BETJEF)
  416.       ELSE
  417.        CALL INDIC1(DK,0.0D0,IFIS,IPLA,1,BETJEF,NECH0,NECH1)
  418.       ENDIF
  419. C ********************** CALCUL VISCOPLASTIQUE ***********************
  420.       IF (IVIS.EQ.1) THEN
  421.         CALL VISPLA(SIGR,SIGF,DSIGT,NSTRS,DEPS1,DEPS2,
  422.      &                   SIGP,SIGRV,DSTRN,D,BETJEF,VISCO)
  423.       ENDIF
  424. C ********************************************************************
  425.   100 CONTINUE
  426. C
  427.       RETURN
  428.       END
  429.  
  430.  
  431.  
  432.  
  433.  
  434.  
  435.  
  436.  
  437.  

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