Télécharger la plaquette Cast3M


fr en

Télécharger desmnl.eso

Retour à la liste

Numérotation des lignes : 

desmnl

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

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