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  16/04/21    21:15:04     8920
  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 20 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. C
  104. C ---------------- direction de compression ----------------
  105. C
  106.       CALL DRUCO1(SIGF,DFSIG2,BETINSA)
  107. C
  108.       DO 25 I=1,NSTRSS
  109.       AC2(I)=0.D0
  110.       DO 25 J=1,NSTRSS
  111.       AC2(I)=AC2(I)+D(I,J)*DFSIG2(J)
  112.    25 CONTINUE
  113. C------------------------------------------------------------
  114.       F1DF1=0.D0
  115.       DO 30 J=1,NSTRSS
  116.       F1DF1=F1DF1+DFSIG1(J)*AC1(J)
  117.    30 CONTINUE
  118.       F1DF2=0.D0
  119.       DO 31 J=1,NSTRSS
  120.       F1DF2=F1DF2+DFSIG1(J)*AC2(J)
  121.    31 CONTINUE
  122.       F2DF2=0.D0
  123.       DO 32 J=1,NSTRSS
  124.       F2DF2=F2DF2+DFSIG2(J)*AC2(J)
  125.    32 CONTINUE
  126.       F2DF1=0.D0
  127.       DO 33 J=1,NSTRSS
  128.       F2DF1=F2DF1+DFSIG2(J)*AC1(J)
  129.    33 CONTINUE
  130. C
  131.       CALL ENDAME(1,BETINSA)
  132.       CALL FORECR(DK1,PAECT,1,SEQ,BETINSA)
  133.       CALL ENDAME(2,BETINSA)
  134.       CALL FORECR(DK2,PAECC,2,SEQ,BETINSA)
  135. C
  136.       DJAC0(1,1)=-(F1DF1+PAECT)
  137.       DJAC0(1,2)=-(F1DF2)
  138.       DJAC0(2,2)=-(F2DF2+PAECC)
  139.       DJAC0(2,1)=-(F2DF1)
  140. C
  141.       DO 43 I=1,2
  142.       DO 43 J=1,2
  143.       DJI(I,J)=DJAC0(I,J)
  144.    43 CONTINUE
  145. C
  146.       CALL INVMA2(DJI,2,ISING)
  147.       IF (ISING.EQ.1) THEN
  148.        WRITE(*,*)'MATRICE DJI singuliere ds ACTI3'
  149.       ENDIF
  150. C
  151. C ************ DEBUT ITERATION INTERNES *******************
  152. C
  153.    40 CONTINUE
  154. C
  155. C *************** Determination de DK et DLAM ******************
  156. C
  157. C
  158. C
  159.       DLAM1(1)=DLAM0(1)-DJI(1,1)*FCRI0(1)-DJI(1,2)*FCRI0(2)
  160.       DLAM1(2)=DLAM0(2)-DJI(2,1)*FCRI0(1)-DJI(2,2)*FCRI0(2)
  161. C
  162.       IF (DLAM1(1).LE.0.D0) IET=1
  163.       IF (DLAM1(2).LE.0.D0) IEC=1
  164.       IF (IET.EQ.1.OR.IEC.EQ.1) THEN
  165. C       WRITE(*,*)'Dans ACTI3, DLAMDA1 est negatif:',DLAM1(1)
  166. C       WRITE(*,*)'Dans ACTI3, DLAMDA2 est negatif:',DLAM1(2)
  167. C       WRITE(*,*)'A l iteration :',ITER
  168.        GOTO 15
  169.       ENDIF
  170. C
  171.       DK1=DKT+DLAM1(1)
  172.       DK2=DKC+DLAM1(2)
  173. C
  174.       CALL ENDAME(1,BETINSA)
  175.       CALL FORECR(DK1,PAECT,1,SEQ1,BETINSA)
  176.       CALL ENDAME(2,BETINSA)
  177.       CALL FORECR(DK2,PAECC,2,SEQ2,BETINSA)
  178. C
  179. C ************** Determination de DPHI1 et 2 ******************
  180. C
  181.       CALL DRUTR2(SIGE,SEQ1,DPHI1,DLAM1(1),VEC1,BETINSA)
  182.       CALL DRUCO2(SIGE,SEQ2,DPHI2,DLAM1(2),VEC2,BETINSA)
  183. C
  184. C--------------- Cas de l'apex -----------------------------
  185. C
  186.       IF (ABS(DPHI1).LE.10E-10.AND.
  187.      *ABS(DPHI2).GT.10E-10) THEN
  188.       IAPEX=1
  189. C      WRITE(*,*)'IAPEX ds ACTI3 =',IAPEX
  190. C      WRITE(*,*)'Dans l element',IBB
  191. C      WRITE(*,*)'et au point d intégration',IGAU
  192.       DO 50 I=1,NSTRSS
  193.       DO 50 J=1,NSTRSS
  194.       AI(I,J)=0.D0
  195.    50 CONTINUE
  196.         AI(1,1)=1./3.
  197.         AI(1,2)=AI(1,1)
  198.         AI(1,2)=AI(1,1)
  199.         AI(1,3)=AI(1,1)
  200.         AI(2,1)=AI(1,1)
  201.         AI(2,2)=AI(1,1)
  202.         AI(2,3)=AI(1,1)
  203.         AI(3,1)=AI(1,1)
  204.         AI(3,2)=AI(1,1)
  205.         AI(3,3)=AI(1,1)
  206.         GOTO 75
  207.       ENDIF
  208. C
  209. C---------- Cas du critere reduit a un point ----------------
  210. C
  211.       IF (ABS(DPHI2).LE.10E-10) THEN
  212.       IAPEX=2
  213. C      WRITE(*,*)'IAPEX ds ACTI3 =',IAPEX
  214. C      WRITE(*,*)'Dans l element',IBB
  215. C      WRITE(*,*)'et au point d intégration',IGAU
  216.       DO 52 I=1,NSTRSS
  217.       DO 52 J=1,NSTRSS
  218.       AI(I,J)=0.D0
  219.    52 CONTINUE
  220.         AI(1,1)=1./3.
  221.         AI(1,2)=AI(1,1)
  222.         AI(1,2)=AI(1,1)
  223.         AI(1,3)=AI(1,1)
  224.         AI(2,1)=AI(1,1)
  225.         AI(2,2)=AI(1,1)
  226.         AI(2,3)=AI(1,1)
  227.         AI(3,1)=AI(1,1)
  228.         AI(3,2)=AI(1,1)
  229.         AI(3,3)=AI(1,1)
  230.         GOTO 75
  231.       ENDIF
  232. C
  233. C ************** Mise a jour des contraintes ***************
  234. C
  235. C ---------------- calcul de la matrice A ------------------
  236. C
  237.       DO 60 I=1,NSTRSS
  238.       DO 60 J=1,NSTRSS
  239.       A(I,J)=0.D0
  240.    60 CONTINUE
  241. C
  242.       DG=YOUN/(1.D0+XNU)
  243. C
  244.       A(1,1)=1.D0+2.*(DLAM1(1)*DG)/2.D0/DPHI1
  245.       A(1,1)=A(1,1)+2.*(DLAM1(2)*DG)/2.D0/DPHI2
  246.       A(2,2)=A(1,1)
  247.       A(3,3)=A(1,1)
  248.       A(1,2)=-(DLAM1(1)*DG)/2.D0/DPHI1-(DLAM1(2)*DG)/2.D0/DPHI2
  249.       A(1,3)=A(1,2)
  250.       A(2,1)=A(1,2)
  251.       A(2,3)=A(1,2)
  252.       A(3,1)=A(1,2)
  253.       A(3,2)=A(1,2)
  254.       A(4,4)=1.D0+3.*(DLAM1(1)*DG)/2.D0/DPHI1
  255.       A(4,4)=A(4,4)+3.*(DLAM1(2)*DG)/2.D0/DPHI2
  256.       A(5,5)=A(4,4)
  257.       A(6,6)=A(4,4)
  258. C
  259. C -------------- invertion de la matrice A -----------------
  260. C
  261.       CALL ZERO(AI,6,6)
  262.       CALL ZERO(AIM,6,6)
  263. C
  264.       DO 70 I=1,3
  265.       DO 70 J=1,3
  266.       AIM(I,J)=A(I,J)
  267.    70 CONTINUE
  268.       CALL INVMA2(AIM,3,ISING)
  269.       IF (ISING.EQ.1) THEN
  270.        WRITE(*,*)'MATRICE AIM singuliere ds ACTI3'
  271.       ENDIF
  272.       DO 72 I=1,3
  273.       DO 72 J=1,3
  274.       AI(I,J)=AIM(I,J)
  275.    72 CONTINUE
  276.       AI(4,4) = 1./A(4,4)
  277.       AI(5,5) = 1./A(5,5)
  278.       AI(6,6) = 1./A(6,6)
  279. C
  280. C -------------- mise a jour des contraintes  ------------
  281. C
  282.    75 CONTINUE
  283. C
  284.       DO 80 I=1,NSTRSS
  285.       DEPSI(I)=SIGE(I)-DLAM1(1)*VEC1(I)
  286.       DEPSI(I)=DEPSI(I)-DLAM1(2)*VEC2(I)
  287.    80 CONTINUE
  288. C
  289.       DO 90 I=1,NSTRSS
  290.       SIGF(I)=0.0D+00
  291.       DO 90 J=1,NSTRSS
  292.       SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J)
  293.    90 CONTINUE
  294. C
  295. C ******** Verification des criteres ****************
  296. C
  297.       CALL DRUTRA(SIGF,SEQTT,BETINSA)
  298.       FCRI1(1) = SEQTT - SEQ1
  299.       CALL DRUCOM(SIGF,SEQCC,BETINSA)
  300.       FCRI1(2) = SEQCC - SEQ2
  301. C
  302.       IF (IAPEX.EQ.2) FCRI1(1)=0.D0
  303. C
  304.       IF (IBROY.EQ.0.AND.(ABS(FCRI1(1)).GE.CRIMAX.OR
  305.      *      .ABS(FCRI1(2)).GE.CRIMAX)) THEN
  306. C      WRITE(*,*)'****************************************'
  307. C      WRITE(*,*)'LE RESIDU DIVERGE AVEC BROYDEN'
  308. C      WRITE(*,*)'on passe donc a la secante'
  309. C      WRITE(*,*)'Dans l element',IBB
  310. C      WRITE(*,*)'et au point d intégration',IGAU
  311. C      WRITE(*,*)'CRIMAX=',CRIMAX
  312. C      WRITE(*,*)'****************************************'
  313.       ITER=ITR
  314.       ENDIF
  315. C
  316. C *******  Compteur sur la methode de resolution ****
  317. C
  318.       IF (IBROY.EQ.0.AND.ITER.EQ.ITR) THEN
  319.         IBROY=1
  320.         ITANG=1
  321.         ITER=1
  322.         IAPEX=0
  323.         IET=0
  324.         IEC=0
  325.         GOTO 12
  326.       ENDIF
  327. C
  328. C ******* non convergence **************************
  329. C
  330.       IF ((ABS(FCRI1(1)).GT.PRB.OR.ABS(FCRI1(2)).GT.PRB)
  331.      *.AND.ITER.LT.ITR) THEN
  332.       IF (IBROY.EQ.0) THEN
  333.        CALL BROYDI(DJI,DLAM0,DLAM1,FCRI0,FCRI1)
  334.        DLAM0(1)=DLAM1(1)
  335.        DLAM0(2)=DLAM1(2)
  336.        FCRI0(1)=FCRI1(1)
  337.        FCRI0(2)=FCRI1(2)
  338.        IF (ITER.GE.(ITR-1)) THEN
  339. C        WRITE(*,*)'***********************'
  340. C        WRITE(*,*)'BROYDEN n a pas aboutit'
  341. C        WRITE(*,*)'ITER=',ITER
  342. C        WRITE(*,*)'FCRIT=',FCRI0(1)
  343. C        WRITE(*,*)'FCRIC=',FCRI0(2)
  344. C        WRITE(*,*)'Dans l element',IBB
  345. C        WRITE(*,*)'et au point d intégration',IGAU
  346. C        WRITE(*,*)'***********************'
  347.        ENDIF
  348.        ITER=ITER+1
  349.        GOTO 40
  350.       ENDIF
  351.       IF (IBROY.EQ.1.AND.ITANG.EQ.1) THEN
  352.        DLAM0(1)=DLAM1(1)
  353.        DLAM0(2)=DLAM1(2)
  354.        FCRI0(1)=FCRI1(1)
  355.        FCRI0(2)=FCRI1(2)
  356.        IF (ITER.GE.(ITR-5)) THEN
  357. C        WRITE(*,*)'ITER=',ITER
  358. C        WRITE(*,*)'FCRIT=',FCRI0(1)
  359. C        WRITE(*,*)'FCRIC=',FCRI0(2)
  360.        ENDIF
  361.        ITER=ITER+1
  362.        GOTO 18
  363.       ENDIF
  364.       ENDIF
  365.       IF (ITER.GE.ITR.AND.(ABS(FCRI1(1)).GT.PRB2.OR.
  366.      *      ABS(FCRI1(2)).GT.PRB2)) THEN
  367.         WRITE(*,*)'NON CONVERGENCE INTERNE dans BEHAV3'
  368.         WRITE(*,*)'Dans l element',IBB
  369.         WRITE(*,*)'et au point d intégration',IGAU
  370.         WRITE(*,*)'FCRIT=',FCRI0(1)
  371.         WRITE(*,*)'FCRIC=',FCRI0(2)
  372.         WRITE(*,*)'IPLA=',IPLA
  373.         WRITE(*,*)'IFIS=',IFIS
  374. C        STOP
  375.       ENDIF
  376. C
  377. C **************** Fin des iterations internes *******************
  378. C
  379.       DKT=DK1
  380.       DKC=DK2
  381.       SEQT=SEQ1
  382.       SEQC=SEQ2
  383. C
  384. C ********************************************************************
  385.   100 CONTINUE
  386.       IET=0
  387.       IEC=0
  388.       CRIMAX=0.D0
  389. C
  390.       RETURN
  391.       END
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  

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