Télécharger la plaquette Cast3M


fr en

Télécharger micro1.eso

Retour à la liste

Numérotation des lignes : 

micro1

  1. C MICRO1    SOURCE    CHAT      05/01/13    01:44:59     5004
  2.       SUBROUTINE MICRO1(WRK0,NMATT,WRK1,NSTRS,DSIGT,
  3.      1                          NVARI,IDECAL,KERRE,LIMPR)
  4.       IMPLICIT INTEGER(I-N)
  5.  
  6. -INC PPARAM
  7. -INC CCOPTIO
  8.       INTEGER NSTRS,NVARI,KERRE,ISTRS,NMATT,ITER,NSIG,IDECAL
  9.       LOGICAL CONV,LCP,LIMPR
  10.       REAL*8 SEQ,TRDFDS,YOUNG,XNU,DEUXMU,LAMB,R0,DR,DR0,RINF,RI
  11.       REAL*8 FPC,MP,EPSE,R,TEST,TRSIG,F0,DELTL,DELTLT
  12.       REAL*8 DEPS3,EPS3,DFDELTL,DFDEPS3,DS3DELTL,DS3DEPS3
  13.       REAL*8 A,B,C,A1,B1,C1,A2,B2,C2
  14.       CHARACTER*8 CMATE
  15.       REAL*8 DSIGT(6)
  16. *
  17.       SEGMENT WRK0
  18.        REAL*8 XMAT(NMATT)
  19.       ENDSEGMENT
  20. *
  21.       SEGMENT WRK1
  22.         REAL*8 DDHOOK(LHOOK,LHOOK),SIG0(NSTRS),DEPST(NSTRS)
  23.         REAL*8 SIGF(NSTRS),VAR0(NVARI),VARF(NVARI)
  24.         REAL*8 DEFP(NSTRS),XCAR(ICARA)
  25.       ENDSEGMENT
  26. *
  27.       SEGMENT WRKK1
  28.         REAL*8 SIGE0(NSTRS),DFDS(NSTRS)
  29.       END SEGMENT
  30. *
  31.       SEGINI WRKK1
  32.       YOUNG=XMAT(1)
  33.       XNU=XMAT(2)
  34.       FPC=XMAT(7)
  35.       MP=XMAT(8)
  36.       EPSE=VAR0(1)
  37.       RI=VAR0(2)
  38.       R0=0.4D0*FPC
  39.       RINF=0.6d0*FPC
  40.       LCP=IFOUR.EQ.-2
  41.       DEUXMU=YOUNG/(1.D0+XNU)
  42.       LAMB=XNU*DEUXMU/(1.D0-2.D0*XNU)
  43.       TEST=1.D-6*R0
  44.       DR0=0.D0
  45.       KERRE=0
  46. C
  47. C     calcul de la matrice elastique
  48. C
  49.       CMATE = 'ISOTROPE'
  50.       KCAS=1
  51.       CALL DOHMAS(XMAT,CMATE,IFOUR,NSTRS,KCAS,DDHOOK,IRTD)
  52.       IF ( IRTD .EQ. 0) THEN
  53.         KERRE=123
  54.         return
  55.       end if
  56. C
  57. C      calcul de l'increment de contrainte
  58. C
  59.       CALL MATVE1 (DDHOOK,DEPST,NSTRS,NSTRS,DSIGT,1)
  60. *      print*,'IFOUR=',IFOUR
  61. *      print*,'DDHOOK'
  62. *      do i=1,3
  63. *      print*,(ddhook(i,j),j=1,3)
  64. *      end do
  65. *      do i=1,3
  66. *      print*,'DEPST(',I,')=',DEPST(i)
  67. *      print*,'DSIGT(',I,')=',DSIGT(I)
  68. *      end do
  69. *
  70. *     calcul des parametres du critere de NADAI
  71. *
  72. C-
  73. C- A,B,C : PARAMETRES DU CRITERE DE PLASTICITE (NADAI)
  74. C- ILS SONT DIFFERENTS SUIVANT LE CHARGEMENT
  75. C-
  76.  
  77.       CONV=.FALSE.
  78.       CALL ZDANUL (DEFP,NSTRS)
  79.       ITER=0
  80.       EPS3=0.d0
  81.       DELTLT=0.D0
  82.       R=RI
  83.       A1=0.4d0
  84.       B1=(1.d0 - A1) / (1.D0 + A1)
  85.       C1=(2.D0 * A1) / (1.D0 + A1)
  86.       A2=1.16d0
  87.       B2=(A2 - 1.D0) / (2.D0 * A2 - 1.D0)
  88.       C2=A2 / (2.D0 * A2 - 1.D0)
  89. *
  90. *     NSIG est la dimension du tenseur des contraintes:
  91. *     NSIG=3 sauf 2 en contraintes planes
  92. *
  93.       IF (LCP)THEN
  94.         NSIG=2
  95.       ELSE
  96.         NSIG=3
  97.       END IF
  98. *
  99. *     on recupere les contraintes anelastiques due a l'endommagement
  100. *
  101.       DO ISTRS=1,NSTRS
  102.          SIGE0(ISTRS)=VAR0(ISTRS+IDECAL)
  103.       END DO
  104. *
  105. *     prediction elastique de la contrainte effective
  106. *     contrainte + contrainte anelastique + increment de contrainte elastique
  107. *
  108. *     l'increment de contrainte est decoupe:
  109. *     1) point de contact avec la surface seuil ???
  110. *     2) decoupage du reste de facon a ce que le pas en contrainte
  111. *        equivalente soit inferieur a disons pour le moment E*1.e-04
  112. *
  113. *     en deux coups les gros
  114. *
  115. *
  116. *      DS=0.d0
  117. *      do istrs=1,nstrs
  118. *        ds=ds + dsigt(istrs)**2
  119. *      end do
  120. *      nstep=int(1.e+4*sqrt(ds)/youn)
  121.       DO ISTRS=1,NSTRS
  122.          SIGE0(ISTRS)=SIG0(ISTRS)+SIGE0(ISTRS)+DSIGT(ISTRS)
  123.       END DO
  124.       if(limpr)then
  125.        DO ISTRS=1,NSTRS
  126.          PRINT*,'SIG0(',ISTRS,')=',SIGE0(ISTRS)
  127.          PRINT*,'DSIGT(',ISTRS,')=',DSIGT(ISTRS)
  128.        END DO
  129.       endif
  130.  
  131.       DO WHILE (ITER .LE. 100 .AND. .NOT. CONV)
  132.         ITER=ITER+1
  133. *
  134. *       on calcule la contrainte equivalente de Mises (SEQ)
  135. *               et la trace des contraintes (TRSIG)
  136. *
  137. *       on utilise la fonction modulo:
  138. *       MOD( 1,3)= 1
  139. *       MOD( 2,3)= 2
  140. *       MOD( 3,3)= 0
  141. *       MOD( 4,3)= 1
  142. *       MOD( 5,3)= 2
  143. *       MOD( 6,3)= 0
  144. *
  145.         SEQ=0.D0
  146.         TRSIG=0.d0
  147. * ici NSIG change en 3 le 6/11/95 (clb)
  148.         DO ISTRS=1,3
  149.           TRSIG=TRSIG + SIGE0(ISTRS)
  150.           SEQ=SEQ+SIGE0(ISTRS)*(SIGE0(ISTRS)-SIGE0(MOD(ISTRS,3)+1))
  151.         END DO
  152.         DO ISTRS=4,NSTRS
  153.           SEQ=SEQ + 3.D0*SIGE0(ISTRS)**2
  154.         END DO
  155.         if (limpr) then
  156.           print*,'VMIS**2=',SEQ,'TRSIG=',TRSIG
  157.         end if
  158.         SEQ=SQRT(SEQ)
  159. *
  160. *       en fonction de la valeur de la trace de sigma, on determine les
  161. *       coef b et c a applique a la surface seuil
  162. *
  163.         IF (MOD(ITER-1,3).EQ.0 )then
  164.           IF ( TRSIG .GT. (R+R0)*(C2-C1)/(B2-B1))THEN
  165.             B=B1
  166.             C=C1
  167.           ELSE
  168.             B=B2
  169.             C=C2
  170.           END IF
  171.         ENDIF
  172. *
  173. *       fonction seuil
  174. *
  175.         F0=(SEQ + B*TRSIG)/C - R0 -R
  176.         if (limpr) then
  177.           PRINT*,'B=',B,'C=',C,'F0=',F0
  178.           PRINT*,'SIGMAEQ=',(SEQ + B*TRSIG)/C
  179.         end if
  180. *
  181. *       le critere de convergence est verifie si la valeur absolue
  182. *       de la fonction seuil est plus petite que le critere ou si
  183. *       la fonction seuil est negative pour la premiere iteration
  184. *
  185.         CONV=(ABS(F0).LE.TEST .AND.
  186.      1       (.NOT. LCP .OR. ABS(SIGE0(3)) .LE. TEST ) )
  187.      2        .OR. (F0 .LE. 0.D0 .AND. ITER .EQ. 1)
  188. *
  189. *       dans le cas des contraintes planes on rajoute aussi
  190. *       que la valeur absolue de sigma33 est inferieure au critere
  191. *
  192.         IF (.NOT. CONV) THEN
  193. *
  194. *       calcul de dF/dSigma
  195. *
  196.           DO ISTRS=1,3
  197.             DFDS(ISTRS)=(2.D0*SIGE0(ISTRS)-SIGE0(MOD(ISTRS,3)+1)
  198.      &              -SIGE0(MOD(ISTRS+1,3)+1) )/(2.D0*C*SEQ) +B/C
  199.           END DO
  200.           DO ISTRS=4,NSTRS
  201.             DFDS(ISTRS)=3.D0*SIGE0(ISTRS)/(2.d0*C*SEQ)
  202.           END DO
  203.           TRDFDS=0.D0
  204.           DO ISTRS=1,NSIG
  205.             TRDFDS=TRDFDS+DFDS(ISTRS)
  206.           END DO
  207.           if (limpr) then
  208. *           on verifie que les valeur ppales sont bonnes
  209.             do istrs=1,nstrs
  210.               print*,'dfds(',istrs,')=',dfds(istrs)
  211.             end do
  212.             print*,'trdfds=',trdfds
  213.           end if
  214. *
  215. *        calcul de dsigma/deltl=C * depsp
  216. *        attention !! dans le cas des contraintes planes:epsp33 n'est
  217. *        pas une inconnue directe mais c'est depse33
  218. *
  219. *
  220. *         R=MP*p  -----> DR=MP * Dp  =MP*DELTL
  221. *         R=Rinf(1-exp(-bR))  ---> dR=b(Rinf-R)deltl
  222. *
  223. *          DR = (RINF - R)*MP
  224. *          DR=MP
  225. *          R=MP*(p)**0.5 ---> DR=0.5 * MP *P**(-05)dp
  226.           IF ((EPSE+DELTLT).GE. 1.D-20 )THEN
  227.              DR=0.5*MP*(EPSE+DELTLT)**(-0.5)
  228.           ELSE
  229.              DR=0.D0
  230.           END IF
  231.           DFDELTL=-DR
  232.           DO ISTRS=1,NSIG
  233.             DFDELTL=DFDELTL-DFDS(ISTRS)*(DFDS(ISTRS)*DEUXMU+TRDFDS*LAMB)
  234.           END DO
  235.           DO ISTRS=4,NSTRS
  236.             DFDELTL=DFDELTL - 2.d0*DEUXMU * DFDS(ISTRS)**2
  237.           END DO
  238.           IF (LCP) THEN
  239.             DFDELTL=DFDELTL-DFDS(3)*TRDFDS*LAMB
  240.             DFDEPS3=LAMB*TRDFDS + (LAMB+DEUXMU)*DFDS(3)
  241.             DS3DELTL=-LAMB*TRDFDS
  242.             DS3DEPS3=LAMB + DEUXMU
  243.             DELTL=(SIGE0(3)*DFDEPS3 - F0*DS3DEPS3)/
  244.      &            (DS3DEPS3*DFDELTL - DS3DELTL*DFDEPS3)
  245.             DEPS3=-(SIGE0(3)-DS3DELTL*DELTL)/DS3DEPS3
  246.             EPS3=EPS3+DEPS3
  247.             DELTLT=DELTLT+DELTL
  248.             DO ISTRS=1,NSIG
  249.                SIGE0(ISTRS)=SIGE0(ISTRS)+LAMB*DEPS3
  250.             END DO
  251.             SIGE0(3)=SIGE0(3)+(LAMB+DEUXMU)*DEPS3-LAMB*DELTL*TRDFDS
  252.           ELSE
  253.             DELTL=-F0/DFDELTL
  254.             DELTLT=DELTLT+DELTL
  255.           END IF
  256.           DO ISTRS=1,NSIG
  257.             SIGE0(ISTRS)=SIGE0(ISTRS)
  258.      &               -DELTL*(DEUXMU*DFDS(ISTRS)+LAMB*TRDFDS)
  259.           END DO
  260.           DO ISTRS=4,NSTRS
  261.             SIGE0(ISTRS)=SIGE0(ISTRS)-DELTL*DEUXMU*DFDS(ISTRS)
  262.           END DO
  263.         ELSE
  264.                 DELTL=0.D0
  265.         END IF
  266.         R=R + DR*DELTL
  267.         DO ISTRS=1,3
  268.           DEFP(ISTRS)=DEFP(ISTRS)+DELTL*DFDS(ISTRS)
  269.         END DO
  270.         DO ISTRS=4,NSTRS
  271.           DEFP(ISTRS)=DEFP(ISTRS)+2.D0*DELTL*DFDS(ISTRS)
  272.         END DO
  273.          if (limpr) then
  274.            print*,'ITER=',ITER,'F0=',F0,'DELTL=', DELTL
  275.            print*,'DELTLT=',DELTLT,'R=',R
  276.            do istrs=1,nstrs
  277.              print*,'SIGMA(',ISTRS,')=', SIGE0(ISTRS)
  278.            end do
  279.          end if
  280.       END DO
  281. *      EN CAS DE NON CONVERGENCE: ON MET KERRE a 1
  282.       IF (.NOT. CONV) THEN
  283.               KERRE=1
  284.       END IF
  285. *      PRINT*,'CONVERGENCE OBTENUE APRES ',ITER, 'ITERATIONS'
  286. *
  287. *  on remplit les variables finales
  288. *
  289.       VARF(1)=EPSE+DELTLT
  290.       VARF(2)=R
  291. *
  292. *  ON REND LES CONTRAINTES
  293. *
  294.       DO ISTRS=1,NSTRS
  295.         SIGF(ISTRS)=SIGE0(ISTRS)
  296. *        PRINT*,'SIGF(',ISTRS,')=',SIGF(ISTRS)
  297.       END DO
  298. *
  299. * ON REND L'INCREMENT DE DEFORMATIONS PLASTIQUES
  300. *
  301.  
  302.       SEGSUP WRKK1
  303.       RETURN
  304.       END
  305. *
  306.  
  307.  
  308.  
  309.  

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