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    JK148537  23/08/21    21:15:12     11723          
  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. C
  30. C
  31.       DIMENSION XD0(3,3),DF(3,3),POSEPSI2(3,3),EPS33(3,3)
  32.       DIMENSION SIGTILDE33(3,3),SIGDEV33(3,3),SIGD33(3,3),SIGTILDE(6)
  33.       DIMENSION UNMOINSDVP(3),HM1(3,3),DFD(3,3),UNMOINSD(3,3),DFVP(3)
  34.       DIMENSION ROT(3,3),ROTT(3,3),VECP33T(3,3),TRAV(3,3),VECP33(3,3)
  35.       DIMENSION EPSIPOS(3),HM1DIAG(3),EPS(2,2),DELTAF22(2,2),EPSIPP(3)
  36.       DIMENSION VPDF(2),VECPDF(2,2),EPSTF(6)
  37. C
  38.       PARAMETER (XZERO=0.D0 , UN=1.D0 , DEUX=2.D0, XPETIT=1.D-12)
  39. C
  40.       SEGMENT IECOU
  41. *      COMMON/IECOU/NYOG,NYNU,NYALFA,NYSMAX,NYN,NYM,NYKK,
  42.        INTEGER icow1,icow2,icow3,icow4,icow5,icow6,icow7,
  43. C      INTEGER NYOG, NYNU, NYALFA,NYSMAX,NYN, NYM,  NYKK,
  44.      1  icow8,icow9,icow10,icow11,icow12,icow13,icow14,icow15,icow16,
  45. C    .  NYALF1,NYBET1,NYR,  NYA,   NYRHO,NSIGY,  NNKX,  NYKX,  IND,
  46.      2 icow17,icow18,icow19,icow20,icow21,icow22,icow23,icow24,
  47. C    .  NSOM,  NINV, NINCMA,NCOMP, JELEM, LEGAUS,INAT,  NCXMAT,
  48.      3 icow25,icow26,icow27,icow28,icow29,icow30,ICARA,
  49. C    .  LTRAC, MFR,  IELE,  NHRM,   NBNN, NBELEM,ICARA,
  50.      4  icow32,icow33,NSTRS1,MFR1,icow36,icow37,icow38,
  51. C    .  LW2,   NDEF,  NSTRSS,MFR1,  NBGMAT,NELMAT,MSOUPA,
  52.      5  icow39,icow40,icow41,icow42,icow43,icow44
  53. C    .  NUMAT1,LENDO, NBBB,  NNVARI,KERR1, MELEME
  54.         INTEGER icow45,icow46,icow47,icow48,icow49,icow50,
  55.      .          icow51,icow52,icow53,icow54,icow55,icow56
  56.      .          icow57,icow58
  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.       DC=XMAT(10)
  71.       Xk = XMAT(14)
  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 Debut de l'intégration numerique 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.             PRINT*, 'IFAILURE=', IFAILURE
  256. C           mise a jour de XD0 pour traiter IFAILURE=2 dans le meme pas
  257.             DO I=1,3
  258.              DO J=1,3
  259.               XD0(I,J)=DF(I,J)
  260.              ENDDO
  261.             ENDDO
  262.             TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3)
  263.             TRADF=XD0(1,1)+XD0(2,2)+XD0(3,3)
  264.             CALL JACOB3(XD0,3,DFVP,ROT)
  265.           ENDIF
  266.         ENDIF
  267. C
  268.         IF(IFAILURE.EQ.2) THEN
  269. C         deux directions d'endom sont fixee a Dc
  270.           CALL TRSPOD(ROT,3,3,ROTT)
  271. C
  272.           DF33=TRADF-DEUX*Dc
  273. C         print*, 'DF33=',DF33
  274.           IF(DF33.LT.XZERO) THEN
  275.             DF33=XPETIT
  276.           ENDIF
  277. C         On corrige l'endommagement
  278.           CALL ZERO(DF,3,3)
  279.           DF(1,1)=Dc
  280.           DF(2,2)=Dc
  281.           DF(3,3)=DF33
  282. C
  283.           CALL MATMAT(DF,ROTT,3,3,3,TRAV)
  284.           CALL MATMAT(ROT,TRAV,3,3,3,DF)
  285. C
  286. C         On verifie la troisieme direction d'endommagement
  287.           IF(DF33.GE.Dc)THEN
  288.             IFAILURE=3
  289. C           mise a jour de XD0 pour traiter IFAILURE=3 dans le même pas
  290.             DO I=1,3
  291.              DO J=1,3
  292.               XD0(I,J)=DF(I,J)
  293.              ENDDO
  294.             ENDDO
  295.             TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3)
  296.             CALL JACOB3(XD0,3,DFVP,ROT)
  297.           ENDIF
  298.         ENDIF
  299. C
  300.         IF(IFAILURE.EQ.3)THEN
  301. C         On est totalement rompu, endommagement isotrope egal a DC
  302.           CALL TRSPOD(ROT,3,3,ROTT)
  303. C         On corrige l'endommagement
  304.           CALL ZERO(DF,3,3)
  305.           DF(1,1)=Dc
  306.           DF(2,2)=Dc
  307.           DF(3,3)=Dc
  308. C
  309.           CALL MATMAT(DF,ROTT,3,3,3,TRAV)
  310.           CALL MATMAT(ROT,TRAV,3,3,3,DF)
  311.         ENDIF
  312. C
  313. C   Fin de la correction d'endommagement
  314. C
  315. C
  316.         DH=(DF(1,1)+DF(2,2)+DF(3,3))/3.D0
  317. C
  318.       ENDIF
  319. C
  320. C  Fin des differents cas (elasticite ou endommagement)
  321. C
  322. C     Calcul de la contrainte finale
  323. C
  324.       DO I=1,3
  325.         DO J=1,3
  326.           UNMOINSD(I,J)=-UN*DF(I,J)
  327.         ENDDO
  328.       ENDDO
  329.       UNMOINSD(1,1)=UN-DF(1,1)
  330.       UNMOINSD(2,2)=UN-DF(2,2)
  331.       UNMOINSD(3,3)=UN-DF(3,3)
  332.       CALL JACOB3(UNMOINSD,3,UNMOINSDVP,VECP33)
  333.       HM1DIAG(1)=SQRT(UNMOINSDVP(1))
  334.       HM1DIAG(2)=SQRT(UNMOINSDVP(2))
  335.       HM1DIAG(3)=SQRT(UNMOINSDVP(3))
  336.       CALL PRODT2(HM1,HM1DIAG,VECP33,3)
  337. C
  338. C    contrainte effective
  339.       CALL MATVE1(DDHOOK,EPSTF,NSTRS,NSTRS,SIGTILDE,2)
  340.       CALL ENDOCB(SIGTILDE,SIGTILDE33,1,IFOUR)
  341.       CALL MATMAT(HM1,SIGTILDE33,3,3,3,TRAV)
  342.       CALL MATMAT(TRAV,HM1,3,3,3,SIGD33)
  343.       XTRAV=XZERO
  344.       DO I=1,3
  345.         DO J=1,3
  346.           XTRAV=XTRAV+UNMOINSD(I,J)*SIGTILDE33(I,J)
  347.         ENDDO
  348.       ENDDO
  349.       DO I=1,3
  350.         DO J=1,3
  351.           SIGD33(I,J)=SIGD33(I,J)-(XTRAV/(3.D0*(UN-DH)))*UNMOINSD(I,J)
  352.         ENDDO
  353.       ENDDO
  354.       TRASIGTILDE=SIGTILDE(1)+SIGTILDE(2)+SIGTILDE(3)
  355. C
  356.       IF(TRASIGTILDE.GE.XZERO) THEN
  357.         COEF=UN-etaT*DH
  358.         IF(COEF.LE.XZERO) THEN
  359.           COEF=UN-Dc
  360.         ENDIF
  361.         TRASIG=COEF*TRASIGTILDE
  362.       ELSE
  363.         TRASIG=(UN-etaC*DH)*TRASIGTILDE
  364.       ENDIF
  365. C
  366.       SIGF(1)=SIGD33(1,1)+TRASIG/3.d0
  367.       SIGF(2)=SIGD33(2,2)+TRASIG/3.d0
  368.       SIGF(3)=SIGD33(3,3)+TRASIG/3.d0
  369.       SIGF(4)=SIGD33(1,2)
  370.       SIGF(5)=SIGD33(1,3)
  371.       SIGF(6)=SIGD33(2,3)
  372. C
  373. C
  374. C stockage des variables du probleme
  375. C
  376.       VARF(2)=DF(1,1)
  377.       VARF(3)=DF(2,2)
  378.       VARF(4)=DF(3,3)
  379.       VARF(5)=DF(1,2)
  380.       VARF(6)=DF(1,3)
  381.       VARF(7)=DF(2,3)
  382.       VARF(8)=ROT(1,1)
  383.       VARF(9)=ROT(1,2)
  384.       VARF(10)=ROT(1,3)
  385.       VARF(11)=ROT(2,1)
  386.       VARF(12)=ROT(2,2)
  387.       VARF(13)=ROT(2,3)
  388.       VARF(14)=ROT(3,1)
  389.       VARF(15)=ROT(3,2)
  390.       VARF(16)=ROT(3,3)
  391. C
  392. C flag pour l'endommagement maximum
  393.       VARF(17)=IFAILURE
  394. C
  395. C  Deformations totales finales
  396.       DO ISTRS=1,NSTRS
  397.        VARF(17+ISTRS)=EPSTF(ISTRS)
  398.       ENDDO
  399. C
  400. 2000  CONTINUE
  401. C
  402. C
  403.       RETURN
  404.       END
  405. C
  406.  
  407.  
  408.  
  409.  
  410.  
  411.  
  412.  
  413.  
  414.  
  415.  
  416.  
  417.  

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