Télécharger la plaquette Cast3M


fr en

Télécharger desmor.eso

Retour à la liste

Numérotation des lignes : 

desmor

  1. C DESMOR    SOURCE    OF166741  25/11/04    21:15:42     12349          
  2. C DESMORAT (V2)
  3.       SUBROUTINE DESMOR(WRK52,WRK53,WRK54,NVARI,iecou)
  4. C
  5. C     variables en entree
  6. C
  7. C     WRK52, WRK53, WRK54 pointeur actif a l'entree
  8. C     NVARI      nombre de variables internes (doit etre egal a ?)
  9. C
  10. C     ISTEP      flag utilise pour separer les etapes dans un calcul non local
  11. C                ISTEP=0 -----> calcul local
  12. C                ISTEP=1 -----> calcul non local etape 1 on calcule les seuils
  13. C                ISTEP=2 -----> calcul non local etape 2 on continue le calcul
  14. C                               a partir des seuils moyennes
  15. C
  16. C     variables en sortie
  17. C
  18. C     VARINF      variables internes finales
  19. C
  20. C     SIGMAF      contraintes finales
  21. C
  22. C
  23.       IMPLICIT REAL*8(A-H,O-Z)
  24. C
  25.  
  26. -INC PPARAM
  27. -INC CCOPTIO
  28. -INC DECHE
  29.  
  30. -INC TECOU
  31.  
  32.       DIMENSION XD0(3,3),DF(3,3),POSEPSI2(3,3),EPS33(3,3)
  33.       DIMENSION SIGTILDE33(3,3),SIGDEV33(3,3),SIGD33(3,3),SIGTILDE(6)
  34.       DIMENSION UNMOINSDVP(3),HM1(3,3),DFD(3,3),UNMOINSD(3,3),DFVP(3)
  35.       DIMENSION ROT(3,3),ROTT(3,3),VECP33T(3,3),TRAV(3,3),VECP33(3,3)
  36.       DIMENSION EPSIPOS(3),HM1DIAG(3),EPS(2,2),DELTAF22(2,2),EPSIPP(3)
  37.       DIMENSION VPDF(2),VECPDF(2,2),EPSTF(6)
  38.  
  39.       PARAMETER (XZERO=0.D0 , UN=1.D0 , DEUX=2.D0, XPETIT=1.D-12)
  40.  
  41. C
  42. C Recuperation des variables initiales dans les tableaux
  43. C
  44.       CMATE= 'ISOTROPE'
  45.       YOUN = XMAT(1)
  46.       XNU  = XMAT(2)
  47.       XK0  = XMAT(5)
  48.       AA   = XMAT(6)
  49.       A = XMAT(7)
  50.       etaC   = XMAT(8)
  51.       etaT   = XMAT(9)
  52.       DC=XMAT(10)
  53.       Xk = XMAT(14)
  54. C
  55. C Recuperation de l'endommagement initial
  56. C
  57.       XD0(1,1)=VAR0(2)
  58.       XD0(2,2)=VAR0(3)
  59.       XD0(3,3)=VAR0(4)
  60.       XD0(1,2)=VAR0(5)
  61.       XD0(1,3)=VAR0(6)
  62.       XD0(2,3)=VAR0(7)
  63.       XD0(2,1)=XD0(1,2)
  64.       XD0(3,1)=XD0(1,3)
  65.       XD0(3,2)=XD0(2,3)
  66. C
  67. C Recuperation de la base d'endommagement
  68. C
  69.       ROT(1,1)=VAR0(8)
  70.       ROT(1,2)=VAR0(9)
  71.       ROT(1,3)=VAR0(10)
  72.       ROT(2,1)=VAR0(11)
  73.       ROT(2,2)=VAR0(12)
  74.       ROT(2,3)=VAR0(13)
  75.       ROT(3,1)=VAR0(14)
  76.       ROT(3,2)=VAR0(15)
  77.       ROT(3,3)=VAR0(16)
  78. C
  79. C Recuperation de l'etat d'endommagement critique
  80. C
  81.       IFAILURE=nint(VAR0(17))
  82. C
  83. C  Matrice de raideur elastique
  84.       CALL ZERO(DDHOOK,6,6)
  85.       AUX0=YOUN/((UN+XNU)*(UN-XNU-XNU))
  86.       AUX=AUX0*(UN-XNU)
  87.       AUX1=AUX0*XNU
  88.       GEGE=0.5D0*YOUN/(UN+XNU)
  89.       DDHOOK(1,1)=AUX
  90.       DDHOOK(2,1)=AUX1
  91.       DDHOOK(3,1)=AUX1
  92.       DDHOOK(1,2)=AUX1
  93.       DDHOOK(2,2)=AUX
  94.       DDHOOK(3,2)=AUX1
  95.       DDHOOK(1,3)=AUX1
  96.       DDHOOK(2,3)=AUX1
  97.       DDHOOK(3,3)=AUX
  98.       DDHOOK(4,4)=GEGE
  99.       DDHOOK(5,5)=GEGE
  100.       DDHOOK(6,6)=GEGE
  101. C
  102. C Deformation totale
  103. C
  104.       DO ISTRS=1,NSTRS
  105.          EPSTF(ISTRS)=VAR0(17+ISTRS)+DEPST(ISTRS)
  106.       ENDDO
  107. C - Trace de la deformation
  108.       TRAEPS=EPSTF(1)+EPSTF(2)+EPSTF(3)
  109. C
  110. C  calcul des deformations principales
  111.       CALL ENDOCB(EPSTF,EPS33,2,IFOUR)
  112.       CALL JACOB3(EPS33,3,EPSIPP,VECP33)
  113.       EPSIPOS(1)= MAX( XZERO , EPSIPP(1) )**2
  114.       EPSIPOS(2)= MAX( XZERO , EPSIPP(2) )**2
  115.       EPSIPOS(3)= MAX( XZERO , EPSIPP(3) )**2
  116.       EPSTILNLOC=SQRT(EPSIPOS(1)+EPSIPOS(2)+EPSIPOS(3))
  117. C
  118.       IF (ISTEP.EQ.0) THEN
  119.          EPSTIL=EPSTILNLOC
  120.          VARF(1)=EPSTIL
  121.       ELSE IF (ISTEP.EQ.1) THEN
  122.          VARF(1)=EPSTILNLOC
  123.          DO I=2,23
  124.           VARF(I)=VAR0(I)
  125.          ENDDO
  126.          GOTO 2000
  127.       ELSE IF (ISTEP.EQ.2) THEN
  128. c Over nonlocal, coef 2
  129. C         EPSTIL=DEUX*VAR0(1)-EPSTILNLOC
  130.          EPSTIL=VAR0(1)
  131.       ENDIF
  132. C
  133. C**********************************************************
  134. C Debut de l'intégration numerique du comportement
  135. C**********************************************************
  136. C
  137. C     calcul des seuils
  138. C
  139. C  partie positive de la matrice des deformations
  140.       CALL PRODT2(POSEPSI2,EPSIPOS,VECP33,3)
  141. C
  142. C calcul de la trace de D au pas n
  143.       TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3)
  144. C
  145. C  calcul du seuil
  146.       XKAPPA=A*TAN(TRAD0/(A*AA)+ATAN(XK0/A))
  147.       FCRIT=EPSTIL+Xk*TRAEPS-XKAPPA
  148. C
  149. C  TEST sur le seuil
  150.       IF(FCRIT.LT.XZERO) THEN
  151. C       on est elastique
  152. C       endommagement et deformations permanentes inchanges
  153.         DO J=1,3
  154.          DO I=1,3
  155.           DF(I,J)=XD0(I,J)
  156.          ENDDO
  157.         ENDDO
  158.         DH=(DF(1,1)+DF(2,2)+DF(3,3))/3.D0
  159.       ELSE
  160. C
  161. C       debut de la correction d'endom.
  162. C
  163. C       Calcul de Kappa^(-1)=Tr Dn+1
  164. C
  165.         TRADF=A*AA*(ATAN((EPSTIL+Xk*TRAEPS)/A)-ATAN(XK0/A))
  166. C
  167. C       test sur l'eventuelle presence d'une direction d'endommagement fixee
  168. C
  169.         IF(IFAILURE.EQ.0)THEN
  170. C         actualisation de l'endommagement, Dc pas encore atteint (cas "normal")
  171.           IF (EPSTIL.EQ.XZERO) THEN
  172.            DLAMBDA=XZERO
  173.           ELSE
  174.            DLAMBDA=(TRADF-TRAD0)/(EPSTIL**2)
  175.           ENDIF
  176.           IF(DLAMBDA.LT.XZERO) THEN
  177.             DLAMBDA=XZERO
  178.           ENDIF
  179.           CALL ZERO(DF,3,3)
  180.           DO J=1,3
  181.            DO I=1,3
  182.             DF(I,J)=XD0(I,J)+DLAMBDA*POSEPSI2(I,J)
  183.            ENDDO
  184.           ENDDO
  185. C         On borne l'endommagement D(1) a Dc si D>Dc et Fixed-crack
  186.           CALL JACOD3(DF,3,DFVP)
  187.           IF(DFVP(1).GE.DC) THEN
  188.             IFAILURE=1
  189. C           mise a jour de XD0 pour traiter IFAILURE=1 dans le meme pas
  190.             DO J=1,3
  191.              DO I=1,3
  192.               XD0(I,J)=DF(I,J)
  193.              ENDDO
  194.             ENDDO
  195.             TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3)
  196.             TRADF=XD0(1,1)+XD0(2,2)+XD0(3,3)
  197.             CALL JACOB3(XD0,3,DFVP,ROT)
  198.            ENDIF
  199.         ENDIF
  200. C
  201.         IF(IFAILURE.EQ.1)THEN
  202. C       une direction d'endom est fixee a Dc
  203. C
  204. C       On calcule le multiplicateur d'endommagement dans la bonne base
  205. C       On met les deformations positives dans le repère principal
  206. C                 de l'endommagement
  207.           CALL JACOD3(XD0,3,DFVP)
  208.           CALL TRSPOD(ROT,3,3,ROTT)
  209.           CALL MATMAT(POSEPSI2,ROT,3,3,3,TRAV)
  210.           CALL MATMAT(ROTT,TRAV,3,3,3,POSEPSI2)
  211. C
  212.           IF(ABS(POSEPSI2(2,2)+POSEPSI2(3,3)).GT.XPETIT) THEN
  213.             DLAMBDA=(TRADF-TRAD0)/(POSEPSI2(2,2)+POSEPSI2(3,3))
  214.           ELSE
  215.             DLAMBDA=XZERO
  216.           ENDIF
  217.           IF(DLAMBDA.LT.XZERO) THEN
  218.             DLAMBDA=XZERO
  219.           ENDIF
  220. C
  221. C         On corrige l'endommagement
  222.           CALL ZERO(DF,3,3)
  223.           DF(1,1)=Dc
  224.           DF(2,2)=DFVP(2)+DLAMBDA*POSEPSI2(2,2)
  225.           DF(3,3)=DFVP(3)+DLAMBDA*POSEPSI2(3,3)
  226.           DF(2,3)=DLAMBDA*POSEPSI2(2,3)
  227.           DF(3,2)=DF(2,3)
  228. C
  229.           CALL MATMAT(DF,ROTT,3,3,3,TRAV)
  230.           CALL MATMAT(ROT,TRAV,3,3,3,DF)
  231. C
  232. C         On verifie la deuxieme direction d'endommagement
  233.           CALL JACOD3(DF,3,DFVP)
  234.           IF(DFVP(2).GE.Dc)THEN
  235.             IFAILURE=2
  236.             PRINT*, 'IFAILURE=', IFAILURE
  237. C           mise a jour de XD0 pour traiter IFAILURE=2 dans le meme pas
  238.             DO J=1,3
  239.               DO I=1,3
  240.                 XD0(I,J)=DF(I,J)
  241.              ENDDO
  242.             ENDDO
  243.             TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3)
  244.             TRADF=XD0(1,1)+XD0(2,2)+XD0(3,3)
  245.             CALL JACOB3(XD0,3,DFVP,ROT)
  246.           ENDIF
  247.         ENDIF
  248. C
  249.         IF(IFAILURE.EQ.2) THEN
  250. C         deux directions d'endom sont fixee a Dc
  251.           CALL TRSPOD(ROT,3,3,ROTT)
  252. C
  253.           DF33=TRADF-DEUX*Dc
  254. C         print*, 'DF33=',DF33
  255.           IF(DF33.LT.XZERO) THEN
  256.             DF33=XPETIT
  257.           ENDIF
  258. C         On corrige l'endommagement
  259.           CALL ZERO(DF,3,3)
  260.           DF(1,1)=Dc
  261.           DF(2,2)=Dc
  262.           DF(3,3)=DF33
  263. C
  264.           CALL MATMAT(DF,ROTT,3,3,3,TRAV)
  265.           CALL MATMAT(ROT,TRAV,3,3,3,DF)
  266. C
  267. C         On verifie la troisieme direction d'endommagement
  268.           IF(DF33.GE.Dc)THEN
  269.             IFAILURE=3
  270. C           mise a jour de XD0 pour traiter IFAILURE=3 dans le même pas
  271.             DO J=1,3
  272.               DO I=1,3
  273.                 XD0(I,J)=DF(I,J)
  274.              ENDDO
  275.             ENDDO
  276.             TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3)
  277.             CALL JACOB3(XD0,3,DFVP,ROT)
  278.           ENDIF
  279.         ENDIF
  280. C
  281.         IF(IFAILURE.EQ.3)THEN
  282. C         On est totalement rompu, endommagement isotrope egal a DC
  283.           CALL TRSPOD(ROT,3,3,ROTT)
  284. C         On corrige l'endommagement
  285.           CALL ZERO(DF,3,3)
  286.           DF(1,1)=Dc
  287.           DF(2,2)=Dc
  288.           DF(3,3)=Dc
  289. C
  290.           CALL MATMAT(DF,ROTT,3,3,3,TRAV)
  291.           CALL MATMAT(ROT,TRAV,3,3,3,DF)
  292.         ENDIF
  293. C
  294. C   Fin de la correction d'endommagement
  295. C
  296.         DH=(DF(1,1)+DF(2,2)+DF(3,3))/3.D0
  297. C
  298.       ENDIF
  299. C
  300. C  Fin des differents cas (elasticite ou endommagement)
  301. C
  302. C     Calcul de la contrainte finale
  303. C
  304.       DO J=1,3
  305.         DO I=1,3
  306.           UNMOINSD(I,J)=-UN*DF(I,J)
  307.         ENDDO
  308.       ENDDO
  309.       UNMOINSD(1,1)=UN-DF(1,1)
  310.       UNMOINSD(2,2)=UN-DF(2,2)
  311.       UNMOINSD(3,3)=UN-DF(3,3)
  312.       CALL JACOB3(UNMOINSD,3,UNMOINSDVP,VECP33)
  313.       HM1DIAG(1)=SQRT(UNMOINSDVP(1))
  314.       HM1DIAG(2)=SQRT(UNMOINSDVP(2))
  315.       HM1DIAG(3)=SQRT(UNMOINSDVP(3))
  316.       CALL PRODT2(HM1,HM1DIAG,VECP33,3)
  317. C
  318. C    contrainte effective
  319.       CALL MATVE1(DDHOOK,EPSTF,NSTRS,NSTRS,SIGTILDE,2)
  320.       CALL ENDOCB(SIGTILDE,SIGTILDE33,1,IFOUR)
  321.       CALL MATMAT(HM1,SIGTILDE33,3,3,3,TRAV)
  322.       CALL MATMAT(TRAV,HM1,3,3,3,SIGD33)
  323.       XTRAV=XZERO
  324.       DO J=1,3
  325.         DO I=1,3
  326.           XTRAV=XTRAV+UNMOINSD(I,J)*SIGTILDE33(I,J)
  327.         ENDDO
  328.       ENDDO
  329.       r_z = XTRAV/(3.D0*(UN-DH))
  330.       DO J=1,3
  331.         DO I=1,3
  332.           SIGD33(I,J)=SIGD33(I,J)-r_z*UNMOINSD(I,J)
  333.         ENDDO
  334.       ENDDO
  335.       TRASIGTILDE=SIGTILDE(1)+SIGTILDE(2)+SIGTILDE(3)
  336. C
  337.       IF(TRASIGTILDE.GE.XZERO) THEN
  338.         COEF=UN-etaT*DH
  339.         IF(COEF.LE.XZERO) THEN
  340.           COEF=UN-Dc
  341.         ENDIF
  342.         TRASIG=COEF*TRASIGTILDE/3.D0
  343.       ELSE
  344.         TRASIG=(UN-etaC*DH)*TRASIGTILDE/3.D0
  345.       ENDIF
  346. C
  347.       SIGF(1)=SIGD33(1,1)+TRASIG
  348.       SIGF(2)=SIGD33(2,2)+TRASIG
  349.       SIGF(3)=SIGD33(3,3)+TRASIG
  350.       SIGF(4)=SIGD33(1,2)
  351.       SIGF(5)=SIGD33(1,3)
  352.       SIGF(6)=SIGD33(2,3)
  353. C
  354. C stockage des variables du probleme
  355. C
  356.       VARF(2)=DF(1,1)
  357.       VARF(3)=DF(2,2)
  358.       VARF(4)=DF(3,3)
  359.       VARF(5)=DF(1,2)
  360.       VARF(6)=DF(1,3)
  361.       VARF(7)=DF(2,3)
  362.       VARF(8)=ROT(1,1)
  363.       VARF(9)=ROT(1,2)
  364.       VARF(10)=ROT(1,3)
  365.       VARF(11)=ROT(2,1)
  366.       VARF(12)=ROT(2,2)
  367.       VARF(13)=ROT(2,3)
  368.       VARF(14)=ROT(3,1)
  369.       VARF(15)=ROT(3,2)
  370.       VARF(16)=ROT(3,3)
  371. C
  372. C flag pour l'endommagement maximum
  373.       VARF(17)=IFAILURE
  374. C
  375. C  Deformations totales finales
  376.       DO ISTRS=1,NSTRS
  377.        VARF(17+ISTRS)=EPSTF(ISTRS)
  378.       ENDDO
  379. C
  380. 2000  CONTINUE
  381.  
  382.       RETURN
  383.       END
  384.  
  385.  
  386.  

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