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

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