Télécharger la plaquette Cast3M


fr en

Télécharger microp.eso

Retour à la liste

Numérotation des lignes : 

microp

  1. C MICROP    SOURCE    CB215821  17/11/30    21:16:49     9639           
  2.       SUBROUTINE MICROP (WRK0,WRK1,NSTRS,NVARI,NMATT,
  3.      &                                 ICARA,KERRE,MFR1,IFOURB)
  4. *
  5. *   modele d'endommagement microplan couple a la plasticite
  6. *   C. La Borderie + S. Fichant Oct. 95
  7. *   routines utilisees:
  8. *   micro1: plasticite nadai
  9. *            IDECAL=3 DANS LE CAS ISO IDECAL=8 DANS LE CAS ANISO
  10. *   jacob3: diagonalisation:
  11. *            attention jacob3 modifie la matrice a diagonaliser!!
  12. *   prodt et prodt2
  13. *   attention prodt2 ne fonctionne qu'avec la matrice des V. P. !!
  14. *
  15.       IMPLICIT INTEGER(I-N)
  16. C      IMPLICIT REAL*8(A-H,O-Z)
  17.       INTEGER NSTRS,NVARI,NMATT,ICARA,KERRE,MFR1,IFOURB
  18.       INTEGER ISTRS,I,J
  19.       REAL*8 YOUNG,XNU,EPSD0,BT,LAMB,DEUXMU,ALFA
  20.       REAL*8 DOM(6),SIGAN(6),TRSIG,DEF33(3,3),EPSIPP(3),VECP(3,3)
  21.       REAL*8 VECPT(3,3),DOM33(3,3),DEFRPDE(3,3),DOMRPDE(3,3)
  22.       REAL*8 D1,D2,D3,DOM3(3)
  23.       REAL*8 SIGPP(3),SIGPM(3),SIG33(3,3),SIG33P(3,3),SIG33M(3,3)
  24.       REAL*8 S33PRD(3,3),S33MRD(3,3)
  25.       REAL*8 LAMBDAP(6), LAMBDAM(6)
  26.       REAL*8 DSIGT(6)
  27.       LOGICAL COMP
  28.       INTEGER IDECAL
  29. *
  30.  
  31. -INC PPARAM
  32. -INC CCOPTIO
  33. *
  34.       SEGMENT WRK0
  35.        REAL*8 XMAT(NMATT)
  36.       ENDSEGMENT
  37. *
  38.       SEGMENT WRK1
  39.         REAL*8 DDHOOK(LHOOK,LHOOK),SIG0(NSTRS),DEPST(NSTRS)
  40.         REAL*8 SIGF(NSTRS),VAR0(NVARI),VARF(NVARI)
  41.         REAL*8 DEFP(NSTRS),XCAR(ICARA)
  42.       ENDSEGMENT
  43.       SEGMENT WRKK1
  44.         REAL*8 DEFELA(NSTRS)
  45.         REAL*8 DDEFEL(NSTRS)
  46.       ENDSEGMENT
  47. *
  48. *   on recupere les variable materielles
  49. *
  50.       YOUNG=XMAT(1)
  51.       XNU=XMAT(2)
  52.       EPSD0=XMAT(5)
  53.       BT=XMAT(6)
  54.       ALFA=XMAT(9)
  55.       DEUXMU=YOUNG/(1.D0+XNU)
  56.       LAMB=XNU*DEUXMU/(1.D0-2.D0*XNU)
  57. *      print*,'dans microp'
  58. *
  59. *     recuperation des variables internes d'endommagement
  60. *
  61.       DO ISTRS=1,NSTRS
  62.         DOM(ISTRS)=VAR0(2+ISTRS)
  63.       END DO
  64.       SEGINI WRKK1
  65. *
  66. *     on ecoule plastiquement sur la contrainte effective
  67. *
  68.       IDECAL=8
  69. *      print*,'appel a micro1'
  70.       CALL MICRO1(WRK0,NMATT,WRK1,NSTRS,DSIGT,
  71.      1                       NVARI,IDECAL,KERRE,.false.)
  72.       IF (KERRE .NE. 0) THEN
  73.            print*,'on n''a pas converge dans micro1'
  74.         CALL MICRO1(WRK0,NMATT,WRK1,NSTRS,DSIGT,
  75.      1                       NVARI,IDECAL,KERRE,.true.)
  76.          SEGSUP WRKK1
  77.          RETURN
  78.       ENDIF
  79. *
  80. *     on ecoule en endommagement sur les deformations elastiques
  81. *
  82. *      print*,'apres micro1'
  83. *
  84. *     calcul de l'increment deformations elastiques DDEFEL
  85. *     1) on calcule l'increment de deformations totales avec
  86. *        l'increment de contraintes elastique DSIGT
  87. *     2) on retranche l'increment de deformations plastiques DEFP
  88. *
  89.         TRSIG=0.D0
  90.         DO ISTRS=1,3
  91.           TRSIG=TRSIG+DSIGT(ISTRS)
  92.         END DO
  93.         DO ISTRS=1,3
  94.           DDEFEL(ISTRS)=( (1.D0+XNU)*DSIGT(ISTRS)-XNU*TRSIG)/YOUNG
  95.      1                 - DEFP(ISTRS)
  96.         END DO
  97.         DO ISTRS=4,NSTRS
  98.           DDEFEL(ISTRS)= (1.D0+XNU)*DSIGT(ISTRS)/YOUNG
  99.      1                 - 0.5d0*DEFP(ISTRS)
  100.         END DO
  101. *
  102. *      on diagonalise l'increment de deformations elastiques:
  103. *      1) on met sous forme 3x3 avec endoca
  104. *      2) on diagonalise avec jacob3
  105. *
  106. *       print*,'incr de def el dans rpg'
  107. *       print*,ddefel
  108. *       print*,'avant endoca ddefl'
  109.          CALL ENDOCA (DDEFEL,DEF33,1)
  110. *       print*,'apres endoca ddefl'
  111.          CALL JACOB3 (DEF33,IDIM,EPSIPP,VECP)
  112. *       print*,'apres jacob3 def33'
  113. *       print*,'incr deformations principales'
  114. *       print*,epsipp
  115. *
  116. *     calcul des deformations elastiques
  117. *
  118.          TRSIG=0.D0
  119. *      print*,'-----sigf------'
  120. *      print*,sigf(1),sigf(2),sigf(3)
  121.          DO ISTRS=1,3
  122.             TRSIG=TRSIG+SIGF(ISTRS)
  123.          END DO
  124.          DO ISTRS=1,3
  125.            DEFELA(ISTRS)=( (1.D0+XNU)*SIGF(ISTRS)-XNU*TRSIG)/YOUNG
  126.          END DO
  127. *      print*,'-----defela------'
  128. *      print*,defela(1),defela(2),defela(3)
  129.          DO ISTRS=4,NSTRS
  130.            DEFELA(ISTRS)= (1.D0+XNU)*SIGF(ISTRS)/YOUNG
  131.          END DO
  132. *         print*,'deformations elastiques dans rpg'
  133. *         print*,defela
  134. *
  135. *     on met les deformations sous forme de matrice 3x3
  136. *     puis on ecrit la matrice dans le repere de depsilon:DEFRPDE
  137. *
  138.          CALL ENDOCA (DEFELA,DEF33,1)
  139. *         print*,'deformations elastiques dans rpg(3x3)'
  140. *         print*,def33
  141.          CALL PRODT (DEFRPDE,DEF33,VECP,3,3)
  142. *          print*,'def elast dans rpddeps'
  143. *         print*,defrpde
  144. *
  145. *     on calcule l'endommagement resultant
  146. *
  147.          IF (DEFRPDE(1,1) .GT. EPSD0) THEN
  148.            D1=1.D0-EPSD0/DEFRPDE(1,1)*EXP(BT*( EPSD0 - DEFRPDE(1,1)))
  149.          ELSE
  150.            D1=0.D0
  151.          END IF
  152.          IF (DEFRPDE(2,2) .GT. EPSD0) THEN
  153.            D2=1.D0-EPSD0/DEFRPDE(2,2)*EXP(BT*( EPSD0 - DEFRPDE(2,2)))
  154.          ELSE
  155.            D2=0.D0
  156.          END IF
  157.          IF (DEFRPDE(3,3) .GT. EPSD0) THEN
  158.            D3=1.D0-EPSD0/DEFRPDE(3,3)*EXP(BT*( EPSD0 - DEFRPDE(3,3)))
  159.          ELSE
  160.            D3=0.D0
  161.          END IF
  162. *      print*,defrpde(1,1),epsd0,D1
  163. *      print*,defrpde(2,2),epsd0,D2
  164. *      print*,defrpde(3,3),epsd0,D3
  165. *
  166. *     on met l'endommagement initial dans le meme repere
  167. *
  168.          CALL ENDOCA(DOM,DOM33,1)
  169. *      print*,'prodt DOMRPDE'
  170.          CALL PRODT (DOMRPDE,DOM33,VECP,3,3)
  171. *       print*,'endommagement initial dans rpddeps'
  172. *       print*,domrpde
  173. *
  174. *     et on en deduit l'increment d'endommagement dans RPDE
  175. *
  176.          IF(d1.gt.domrpde(1,1))then
  177.             domrpde(1,1)=d1
  178.          endif
  179.          IF(d2.gt.domrpde(2,2))then
  180.             domrpde(2,2)=d2
  181.          endif
  182.          IF(d3.gt.domrpde(3,3))then
  183.             domrpde(3,3)=d3
  184.          endif
  185. *       print*,'endommagement final dans rpddeps'
  186. *       print*,domrpde
  187. *
  188. *     on remet D dans RPG
  189. *
  190.       DO I=1,3
  191.         DO J=1,3
  192.            VECPT(I,J)=VECP(J,I)
  193.         END DO
  194.       END DO
  195.       call prodt(DOM33,DOMRPDE,VECPT,3,3)
  196. *       print*,'endommagement final dans rpg'
  197. *       print*,dom33
  198.  
  199. *
  200. *     on separe les contraintes effectives en + et - dans rpsigma
  201. *
  202.          CALL ENDOCA (SIGF,SIG33,1)
  203.          CALL JACOB3 (SIG33,3,SIGPP,VECP)
  204. *         print*,'contraintes ppales'
  205. *         print*,sigpp
  206.          DO I=1,3
  207.            IF (SIGPP(I) .LT. 0.D0)THEN
  208.               SIGPM(I)=SIGPP(I)
  209.               SIGPP(I)=0.D0
  210.            ELSE
  211.               SIGPM(I)=0.D0
  212.            END IF
  213.          END DO
  214.          CALL PRODT2(SIG33P,SIGPP,VECP,3)
  215.          CALL PRODT2(SIG33M,SIGPM,VECP,3)
  216. *         print*,'contraintes dans rpg'
  217. *         print*,sig33p
  218. *         print*,sig33m
  219. *
  220. *      on met le tout dans le repere ppal d'endo
  221. *      attention jacob3 modifie la matrice fournie
  222. *       --> on passe une copie
  223.          DO I=1,3
  224.            DO J=1,3
  225.              SIG33(I,J)=DOM33(I,J)
  226.            END DO
  227.          END DO
  228.          CALL JACOB3(SIG33,3,DOM3,VECP)
  229. *         print*,'endom dans rpd'
  230. *         print*,dom3
  231. **************
  232. *  A REVOIR  EN ATENDANT MIEUX ON BORNE LES VALEURS PROPRES DE D ENTRE 0. ET 1.*
  233. **************
  234.          DO I=1,3
  235.            DOM3(I)=MAX(DOM3(I),0.D0)
  236.            DOM3(I)=MIN(DOM3(I),1.D0-1.d-6)
  237.          END DO
  238.          CALL PRODT2 (DOM33,DOM3,VECP,3)
  239. **************
  240. *  FIN A REVOIR
  241. **************
  242. *      print*,'dom3=',dom3
  243. *      print*,'dom33=',dom33
  244.          CALL PRODT (S33PRD,SIG33P,VECP,3,3)
  245.          CALL PRODT (S33MRD,SIG33M,VECP,3,3)
  246.  
  247. *      print*,'s33prd',s33prd(3,3)
  248. *      print*,'s33mrd',s33mrd(3,3)
  249. *      print*,'callambdap'
  250.          COMP=.FALSE.
  251.          CALL CLMBDA(DOM3(1),DOM3(2),DOM3(3),LAMBDAP,YOUNG,XNU,ALFA,
  252.      1                                                            COMP)
  253.          COMP=.TRUE.
  254.          CALL CLMBDA(DOM3(1),DOM3(2),DOM3(3),LAMBDAM,YOUNG,XNU,ALFA,
  255.      1                                                            COMP)
  256. *       print*,'apres callambdap'
  257.        SIG33P(1,1)=LAMBDAP(1)*S33PRD(1,1)+LAMBDAM(1)*S33MRD(1,1)
  258.        SIG33P(1,2)=LAMBDAP(4)*S33PRD(1,2)+LAMBDAM(4)*S33MRD(1,2)
  259.        SIG33P(1,3)=LAMBDAP(5)*S33PRD(1,3)+LAMBDAM(5)*S33MRD(1,3)
  260.        SIG33P(2,1)=SIG33P(1,2)
  261.        SIG33P(2,2)=LAMBDAP(2)*S33PRD(2,2)+LAMBDAM(2)*S33MRD(2,2)
  262.        SIG33P(2,3)=LAMBDAP(6)*S33PRD(2,3)+LAMBDAM(6)*S33MRD(2,3)
  263.        SIG33P(3,1)=SIG33P(1,3)
  264.        SIG33P(3,2)=SIG33P(2,3)
  265.        SIG33P(3,3)=LAMBDAP(3)*S33PRD(3,3)+LAMBDAM(3)*S33MRD(3,3)
  266. *       print*,'contraintes dans RPD'
  267. *       print*,sig33p
  268. *       print*,'sig33p(1,2)',sig33p(1,2)
  269. *      print*,'dsigt',dsigt(1),dsigt(2),dsigt(3)
  270. *      print*,'dsigt',dsigt(4),dsigt(5),dsigt(6)
  271. *      print*,'sigf',sigf(1),sigf(2),sigf(3)
  272. *      print*,'sigf',sigf(4),sigf(5),sigf(6)
  273. *
  274. *     on remet le tout dans le repere global
  275. *
  276. *      print*,'matrice des vect proprs'
  277. *      print*,vecp
  278.       DO I=1,3
  279.         DO J=1,3
  280.            VECPT(I,J)=VECP(J,I)
  281.         END DO
  282.       END DO
  283.       CALL PRODT (SIG33,SIG33P,VECPT,3,3)
  284. *      print*,'contraintes dans rpg'
  285. *      print*,sig33
  286. *
  287. *     on rend les contraintes et les variables internes finales
  288. *
  289.       SIGAN(1)=SIGF(1)-SIG33(1,1)
  290.       SIGF(1)=SIG33(1,1)
  291.       VARF(3)=MAX(DOM33(1,1),0.d0)
  292.       SIGAN(2)=SIGF(2)-SIG33(2,2)
  293.       SIGF(2)=SIG33(2,2)
  294.       VARF(4)=MAX(DOM33(2,2),0.d0)
  295.       SIGAN(3)=SIGF(3)-SIG33(3,3)
  296.       SIGF(3)=SIG33(3,3)
  297.       VARF(5)=MAX(DOM33(3,3),0.d0)
  298.       SIGAN(4)=SIGF(4)-SIG33(1,2)
  299.       SIGF(4)=SIG33(1,2)
  300.       VARF(6)=MAX(DOM33(1,2),0.d0)
  301.       VARF(7)=MAX(DOM33(1,3),0.d0)
  302.       VARF(8)=MAX(DOM33(2,3),0.d0)
  303.       IF(IFOUR.GE.1.OR.IFOUR.LE.-3) THEN
  304.         SIGAN(5)=SIGF(5)-SIG33(1,3)
  305.         SIGF(5)=SIG33(1,3)
  306.         SIGAN(6)=SIGF(6)-SIG33(2,3)
  307.         SIGF(6)=SIG33(2,3)
  308.       ELSE
  309.         SIGAN(5)=0.D0
  310.         SIGAN(6)=0.D0
  311.       END IF
  312.       DO ISTRS=1,6
  313.          VARF(ISTRS+8)=SIGAN(ISTRS)
  314.       END DO
  315. *      print*,'sigf',sigf(1),sigf(2),sigf(3)
  316. *      print*,'sigf',sigf(4),sigf(5),sigf(6)
  317. *      print*,'sigan',sigan(1),sigan(2),sigan(3)
  318. *      print*,'sigan',sigan(4),sigan(5),sigan(6)
  319.       SEGSUP WRKK1
  320.       RETURN
  321.       END
  322.  
  323.  
  324.  
  325.  
  326.  
  327.  

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