Télécharger la plaquette Cast3M


fr en

Télécharger acti3.eso

Retour à la liste

Numérotation des lignes : 

acti3

  1. C ACTI3     SOURCE    CB215821  26/08/24    21:15:07     12622          
  2. C ACTI3
  3.       SUBROUTINE ACTI3(SIG0,SIGF,D,NSTRSS,BETINSA)
  4. C
  5. C ==================================================================
  6. C         Deux criteres de traction compression: DRUCKER PRAGER
  7. C                          3D
  8. C ==================================================================
  9. C CE SOUS-PROGRAMME EST APPELE DANS "BONE3D".
  10. C
  11.       IMPLICIT INTEGER(I-N)
  12.       IMPLICIT REAL*8(A-H,O-Z)
  13.       DIMENSION SIG0(NSTRSS),SIGF(NSTRSS)
  14.       DIMENSION DFSIG1(6),DFSIG2(6),VEC1(6),VEC2(6),SIGE(6)
  15.       DIMENSION AC1(6),AC2(6)
  16.       DIMENSION DJAC0(2,2),DJI(4,4),A(6,6),AI(6,6),AIM(6,6),D(6,6)
  17.       DIMENSION FCRI0(2),FCRI1(2),DEPSI(6),DLAM1(2),DLAM0(2)
  18. C
  19.       SEGMENT BETINSA
  20.        REAL*8 RT,RC,YOUN,XNU,GFT,GFC,CAR
  21.        REAL*8 DKT,DKC,SEQT,SEQC,ENDT,ENDC
  22.        INTEGER IFIS,IPLA,IBB,IGAU
  23.       ENDSEGMENT
  24. C
  25.       IAPEX=0
  26.       PRB=1.D-5
  27.       PRB2=1.D-2
  28.       ITER=1
  29.       ITR = 1500
  30.       IET=0
  31.       IEC=0
  32.       IBROY=0
  33.       ITANG=0
  34.       CRIMAX=0.D0
  35.       SEQ = 0.D0
  36.       SEQ1 = 0.D0
  37.       SEQ2 = 0.D0
  38.       CALL ZERO(SIGE,6,1)
  39. C
  40.       DO 10 I=1,NSTRSS
  41.       SIGE(I)=SIGF(I)
  42.    10 CONTINUE
  43. C
  44. C ************ Le point est fissure pour 1ere fois ***********
  45. C
  46.    12 CONTINUE
  47.       CALL DRUTRA(SIGF,SEQTE,BETINSA)
  48.       CALL DRUCOM(SIGF,SEQCE,BETINSA)
  49.       FCRI0(1) = SEQTE - SEQT
  50.       FCRI0(2) = SEQCE - SEQC
  51.       CRIMAX=ABS(100.D0*FCRI0(1))
  52.       IF (FCRI0(1).LT.0.D0.AND.FCRI0(2).LT.0.D0) THEN
  53.         IET=1
  54.         IEC=1
  55.       ENDIF
  56.       IF (FCRI0(1).LT.0.D0.AND.FCRI0(2).GE.0.D0) THEN
  57.         IET=1
  58.         IEC=0
  59.       ENDIF
  60.       IF (FCRI0(1).GE.0.D0.AND.FCRI0(2).LT.0.D0) THEN
  61.         IET=0
  62.         IEC=1
  63.       ENDIF
  64. C
  65.    15 CONTINUE
  66. C
  67.       DO 16 I=1,NSTRSS
  68.       SIGF(I)=SIGE(I)
  69.    16 CONTINUE
  70. C
  71.       IF (IET.EQ.1.AND.IEC.EQ.1) THEN
  72.         GOTO 100
  73.       ENDIF
  74. C
  75.       IF (IET.EQ.1.AND.IEC.EQ.0) THEN
  76.         CALL ACTI2(SIG0,SIGF,D,NSTRSS,BETINSA)
  77.         GOTO 100
  78.       ENDIF
  79. C
  80.       IF (IET.EQ.0.AND.IEC.EQ.1) THEN
  81.         CALL ACTI1(SIG0,SIGF,D,NSTRSS,BETINSA)
  82.         GOTO 100
  83.       ENDIF
  84. C
  85.       DK1=DKT
  86.       DK2=DKC
  87.       DLAM0(1)=0.D0
  88.       DLAM0(2)=0.D0
  89. C
  90.    18 CONTINUE
  91. C
  92. C ************ CALCUL DU JACOBIEN INITIAL  ********************
  93. C
  94. C ---------------- direction de traction ------------------
  95. C
  96.       CALL DRUTR1(SIGF,DFSIG1,BETINSA)
  97. C
  98.       DO 101 I=1,NSTRSS
  99.       AC1(I)=0.D0
  100.       DO 20 J=1,NSTRSS
  101.       AC1(I)=AC1(I)+D(I,J)*DFSIG1(J)
  102.    20 CONTINUE
  103.   101 CONTINUE
  104. C
  105. C ---------------- direction de compression ----------------
  106. C
  107.       CALL DRUCO1(SIGF,DFSIG2,BETINSA)
  108. C
  109.       DO 102 I=1,NSTRSS
  110.       AC2(I)=0.D0
  111.       DO 25 J=1,NSTRSS
  112.       AC2(I)=AC2(I)+D(I,J)*DFSIG2(J)
  113.    25 CONTINUE
  114.   102 CONTINUE
  115. C------------------------------------------------------------
  116.       F1DF1=0.D0
  117.       DO 30 J=1,NSTRSS
  118.       F1DF1=F1DF1+DFSIG1(J)*AC1(J)
  119.    30 CONTINUE
  120.       F1DF2=0.D0
  121.       DO 31 J=1,NSTRSS
  122.       F1DF2=F1DF2+DFSIG1(J)*AC2(J)
  123.    31 CONTINUE
  124.       F2DF2=0.D0
  125.       DO 32 J=1,NSTRSS
  126.       F2DF2=F2DF2+DFSIG2(J)*AC2(J)
  127.    32 CONTINUE
  128.       F2DF1=0.D0
  129.       DO 33 J=1,NSTRSS
  130.       F2DF1=F2DF1+DFSIG2(J)*AC1(J)
  131.    33 CONTINUE
  132. C
  133.       CALL ENDAME(1,BETINSA)
  134.       CALL FORECR(DK1,PAECT,1,SEQ,BETINSA)
  135.       CALL ENDAME(2,BETINSA)
  136.       CALL FORECR(DK2,PAECC,2,SEQ,BETINSA)
  137. C
  138.       DJAC0(1,1)=-(F1DF1+PAECT)
  139.       DJAC0(1,2)=-(F1DF2)
  140.       DJAC0(2,2)=-(F2DF2+PAECC)
  141.       DJAC0(2,1)=-(F2DF1)
  142. C
  143.       DO 103 I=1,2
  144.       DO 43 J=1,2
  145.       DJI(I,J)=DJAC0(I,J)
  146.    43 CONTINUE
  147.   103 CONTINUE
  148. C
  149.       CALL INVMA2(DJI,2,ISING)
  150.       IF (ISING.EQ.1) THEN
  151.        WRITE(*,*)'MATRICE DJI singuliere ds ACTI3'
  152.       ENDIF
  153. C
  154. C ************ DEBUT ITERATION INTERNES *******************
  155. C
  156.    40 CONTINUE
  157. C
  158. C *************** Determination de DK et DLAM ******************
  159. C
  160. C
  161. C
  162.       DLAM1(1)=DLAM0(1)-DJI(1,1)*FCRI0(1)-DJI(1,2)*FCRI0(2)
  163.       DLAM1(2)=DLAM0(2)-DJI(2,1)*FCRI0(1)-DJI(2,2)*FCRI0(2)
  164. C
  165.       IF (DLAM1(1).LE.0.D0) IET=1
  166.       IF (DLAM1(2).LE.0.D0) IEC=1
  167.       IF (IET.EQ.1.OR.IEC.EQ.1) THEN
  168. C       WRITE(*,*)'Dans ACTI3, DLAMDA1 est negatif:',DLAM1(1)
  169. C       WRITE(*,*)'Dans ACTI3, DLAMDA2 est negatif:',DLAM1(2)
  170. C       WRITE(*,*)'A l iteration :',ITER
  171.        GOTO 15
  172.       ENDIF
  173. C
  174.       DK1=DKT+DLAM1(1)
  175.       DK2=DKC+DLAM1(2)
  176. C
  177.       CALL ENDAME(1,BETINSA)
  178.       CALL FORECR(DK1,PAECT,1,SEQ1,BETINSA)
  179.       CALL ENDAME(2,BETINSA)
  180.       CALL FORECR(DK2,PAECC,2,SEQ2,BETINSA)
  181. C
  182. C ************** Determination de DPHI1 et 2 ******************
  183. C
  184.       CALL DRUTR2(SIGE,SEQ1,DPHI1,DLAM1(1),VEC1,BETINSA)
  185.       CALL DRUCO2(SIGE,SEQ2,DPHI2,DLAM1(2),VEC2,BETINSA)
  186. C
  187. C--------------- Cas de l'apex -----------------------------
  188. C
  189.       IF (ABS(DPHI1).LE.10E-10.AND.
  190.      *ABS(DPHI2).GT.10E-10) THEN
  191.       IAPEX=1
  192. C      WRITE(*,*)'IAPEX ds ACTI3 =',IAPEX
  193. C      WRITE(*,*)'Dans l element',IBB
  194. C      WRITE(*,*)'et au point d intégration',IGAU
  195.       DO 104 I=1,NSTRSS
  196.       DO 50 J=1,NSTRSS
  197.       AI(I,J)=0.D0
  198.    50 CONTINUE
  199.   104 CONTINUE
  200.         AI(1,1)=1./3.
  201.         AI(1,2)=AI(1,1)
  202.         AI(1,2)=AI(1,1)
  203.         AI(1,3)=AI(1,1)
  204.         AI(2,1)=AI(1,1)
  205.         AI(2,2)=AI(1,1)
  206.         AI(2,3)=AI(1,1)
  207.         AI(3,1)=AI(1,1)
  208.         AI(3,2)=AI(1,1)
  209.         AI(3,3)=AI(1,1)
  210.         GOTO 75
  211.       ENDIF
  212. C
  213. C---------- Cas du critere reduit a un point ----------------
  214. C
  215.       IF (ABS(DPHI2).LE.10E-10) THEN
  216.       IAPEX=2
  217. C      WRITE(*,*)'IAPEX ds ACTI3 =',IAPEX
  218. C      WRITE(*,*)'Dans l element',IBB
  219. C      WRITE(*,*)'et au point d intégration',IGAU
  220.       DO 105 I=1,NSTRSS
  221.       DO 52 J=1,NSTRSS
  222.       AI(I,J)=0.D0
  223.    52 CONTINUE
  224.   105 CONTINUE
  225.         AI(1,1)=1./3.
  226.         AI(1,2)=AI(1,1)
  227.         AI(1,2)=AI(1,1)
  228.         AI(1,3)=AI(1,1)
  229.         AI(2,1)=AI(1,1)
  230.         AI(2,2)=AI(1,1)
  231.         AI(2,3)=AI(1,1)
  232.         AI(3,1)=AI(1,1)
  233.         AI(3,2)=AI(1,1)
  234.         AI(3,3)=AI(1,1)
  235.         GOTO 75
  236.       ENDIF
  237. C
  238. C ************** Mise a jour des contraintes ***************
  239. C
  240. C ---------------- calcul de la matrice A ------------------
  241. C
  242.       DO 106 I=1,NSTRSS
  243.       DO 60 J=1,NSTRSS
  244.       A(I,J)=0.D0
  245.    60 CONTINUE
  246.   106 CONTINUE
  247. C
  248.       DG=YOUN/(1.D0+XNU)
  249. C
  250.       A(1,1)=1.D0+2.*(DLAM1(1)*DG)/2.D0/DPHI1
  251.       A(1,1)=A(1,1)+2.*(DLAM1(2)*DG)/2.D0/DPHI2
  252.       A(2,2)=A(1,1)
  253.       A(3,3)=A(1,1)
  254.       A(1,2)=-(DLAM1(1)*DG)/2.D0/DPHI1-(DLAM1(2)*DG)/2.D0/DPHI2
  255.       A(1,3)=A(1,2)
  256.       A(2,1)=A(1,2)
  257.       A(2,3)=A(1,2)
  258.       A(3,1)=A(1,2)
  259.       A(3,2)=A(1,2)
  260.       A(4,4)=1.D0+3.*(DLAM1(1)*DG)/2.D0/DPHI1
  261.       A(4,4)=A(4,4)+3.*(DLAM1(2)*DG)/2.D0/DPHI2
  262.       A(5,5)=A(4,4)
  263.       A(6,6)=A(4,4)
  264. C
  265. C -------------- invertion de la matrice A -----------------
  266. C
  267.       CALL ZERO(AI,6,6)
  268.       CALL ZERO(AIM,6,6)
  269. C
  270.       DO 107 I=1,3
  271.       DO 70 J=1,3
  272.       AIM(I,J)=A(I,J)
  273.    70 CONTINUE
  274.   107 CONTINUE
  275.       CALL INVMA2(AIM,3,ISING)
  276.       IF (ISING.EQ.1) THEN
  277.        WRITE(*,*)'MATRICE AIM singuliere ds ACTI3'
  278.       ENDIF
  279.       DO 108 I=1,3
  280.       DO 72 J=1,3
  281.       AI(I,J)=AIM(I,J)
  282.    72 CONTINUE
  283.   108 CONTINUE
  284.       AI(4,4) = 1./A(4,4)
  285.       AI(5,5) = 1./A(5,5)
  286.       AI(6,6) = 1./A(6,6)
  287. C
  288. C -------------- mise a jour des contraintes  ------------
  289. C
  290.    75 CONTINUE
  291. C
  292.       DO 80 I=1,NSTRSS
  293.       DEPSI(I)=SIGE(I)-DLAM1(1)*VEC1(I)
  294.       DEPSI(I)=DEPSI(I)-DLAM1(2)*VEC2(I)
  295.    80 CONTINUE
  296. C
  297.       DO 109 I=1,NSTRSS
  298.       SIGF(I)=0.0D+00
  299.       DO 90 J=1,NSTRSS
  300.       SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J)
  301.    90 CONTINUE
  302.   109 CONTINUE
  303. C
  304. C ******** Verification des criteres ****************
  305. C
  306.       CALL DRUTRA(SIGF,SEQTT,BETINSA)
  307.       FCRI1(1) = SEQTT - SEQ1
  308.       CALL DRUCOM(SIGF,SEQCC,BETINSA)
  309.       FCRI1(2) = SEQCC - SEQ2
  310. C
  311.       IF (IAPEX.EQ.2) FCRI1(1)=0.D0
  312. C
  313.       IF (IBROY.EQ.0.AND.(ABS(FCRI1(1)).GE.CRIMAX.OR
  314.      *      .ABS(FCRI1(2)).GE.CRIMAX)) THEN
  315. C      WRITE(*,*)'****************************************'
  316. C      WRITE(*,*)'LE RESIDU DIVERGE AVEC BROYDEN'
  317. C      WRITE(*,*)'on passe donc a la secante'
  318. C      WRITE(*,*)'Dans l element',IBB
  319. C      WRITE(*,*)'et au point d intégration',IGAU
  320. C      WRITE(*,*)'CRIMAX=',CRIMAX
  321. C      WRITE(*,*)'****************************************'
  322.       ITER=ITR
  323.       ENDIF
  324. C
  325. C *******  Compteur sur la methode de resolution ****
  326. C
  327.       IF (IBROY.EQ.0.AND.ITER.EQ.ITR) THEN
  328.         IBROY=1
  329.         ITANG=1
  330.         ITER=1
  331.         IAPEX=0
  332.         IET=0
  333.         IEC=0
  334.         GOTO 12
  335.       ENDIF
  336. C
  337. C ******* non convergence **************************
  338. C
  339.       IF ((ABS(FCRI1(1)).GT.PRB.OR.ABS(FCRI1(2)).GT.PRB)
  340.      *.AND.ITER.LT.ITR) THEN
  341.       IF (IBROY.EQ.0) THEN
  342.        CALL BROYDI(DJI,DLAM0,DLAM1,FCRI0,FCRI1)
  343.        DLAM0(1)=DLAM1(1)
  344.        DLAM0(2)=DLAM1(2)
  345.        FCRI0(1)=FCRI1(1)
  346.        FCRI0(2)=FCRI1(2)
  347.        IF (ITER.GE.(ITR-1)) THEN
  348. C        WRITE(*,*)'***********************'
  349. C        WRITE(*,*)'BROYDEN n a pas aboutit'
  350. C        WRITE(*,*)'ITER=',ITER
  351. C        WRITE(*,*)'FCRIT=',FCRI0(1)
  352. C        WRITE(*,*)'FCRIC=',FCRI0(2)
  353. C        WRITE(*,*)'Dans l element',IBB
  354. C        WRITE(*,*)'et au point d intégration',IGAU
  355. C        WRITE(*,*)'***********************'
  356.        ENDIF
  357.        ITER=ITER+1
  358.        GOTO 40
  359.       ENDIF
  360.       IF (IBROY.EQ.1.AND.ITANG.EQ.1) THEN
  361.        DLAM0(1)=DLAM1(1)
  362.        DLAM0(2)=DLAM1(2)
  363.        FCRI0(1)=FCRI1(1)
  364.        FCRI0(2)=FCRI1(2)
  365.        IF (ITER.GE.(ITR-5)) THEN
  366. C        WRITE(*,*)'ITER=',ITER
  367. C        WRITE(*,*)'FCRIT=',FCRI0(1)
  368. C        WRITE(*,*)'FCRIC=',FCRI0(2)
  369.        ENDIF
  370.        ITER=ITER+1
  371.        GOTO 18
  372.       ENDIF
  373.       ENDIF
  374.       IF (ITER.GE.ITR.AND.(ABS(FCRI1(1)).GT.PRB2.OR.
  375.      *      ABS(FCRI1(2)).GT.PRB2)) THEN
  376.         WRITE(*,*)'NON CONVERGENCE INTERNE dans BEHAV3'
  377.         WRITE(*,*)'Dans l element',IBB
  378.         WRITE(*,*)'et au point d intégration',IGAU
  379.         WRITE(*,*)'FCRIT=',FCRI0(1)
  380.         WRITE(*,*)'FCRIC=',FCRI0(2)
  381.         WRITE(*,*)'IPLA=',IPLA
  382.         WRITE(*,*)'IFIS=',IFIS
  383. C        STOP
  384.       ENDIF
  385. C
  386. C **************** Fin des iterations internes *******************
  387. C
  388.       DKT=DK1
  389.       DKC=DK2
  390.       SEQT=SEQ1
  391.       SEQC=SEQ2
  392. C
  393. C ********************************************************************
  394.   100 CONTINUE
  395.       IET=0
  396.       IEC=0
  397.       CRIMAX=0.D0
  398. C
  399.       RETURN
  400.       END
  401.  
  402.  
  403.  
  404.  
  405.  
  406.  
  407.  
  408.  
  409.  

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