Télécharger la plaquette Cast3M


fr en

Télécharger xxcre1.eso

Retour à la liste

Numérotation des lignes : 

xxcre1

  1. C XXCRE1    SOURCE    CB215821  17/11/30    21:17:47     9639           
  2.       SUBROUTINE XXCRE1(TAU,SIG,EPSV,VAR,SIGD,EPSVD,VARD,EPSTHD,
  3.      &DSPT,XMAT,XMAINF,XMASUP,NSTRS,NVARI,NCOMAT,DD,DDV,DDINV,
  4.      &DDINVp,YOUNG,NYOUNG,XNU,NXNU,MFR,XCAR,ICARA,IFOURB,INDIC,TI,
  5.      &TPOINT,TINF,TSUP,ITEST,ITHHER,VARF,IB,IGAU,KERRE)
  6. C-----------------------------------------------------------------------
  7. C Objet: Cette subroutine calcule les derivees des variables internes
  8. C        pour un materiau viscoplastique a endommagement et ecrouissage
  9. C        isotrope en regime anisotherme  de 2 manieres au choix suivant
  10. C        la valeur de ITEST,ITHHER
  11. C-----------------------------------------------------------------------
  12. C
  13. C-----------------------------------------------------------------------
  14. C Entree: TAU pas de temps
  15. C         TI temperature au debut du pas
  16. C         TPOINT derivee de la temperature
  17. C         TINF et TSUP bornes entre lesquelles TI est comprise
  18. C         SIG(NSTRS) tenseur des contraintes
  19. C         EPSV(NSTRS) tenseur des deformations
  20. C         VAR(NVARI) tableau contenant les variables internes
  21. C           pilotant les equations
  22. C           VAR(1)=p ; VAR(2)=D11 ; VAR(3)=D22;  VAR(4)=D33;
  23. C           VAR(5)=D12;  VAR(6)=D13;  VAR(7)=D23; VAR(8)=DMAX
  24. C         EPSTHD  vitesse de dilatation thermique au debut du pas
  25. C         DSPT(NSTRS,NSTRS) vitesse de deformation totale
  26. C         XMAT(NCOMAT) tableau des parametres scalaires du materiau
  27. C           a une temperature T donnee
  28. C         XMAINF(NCOMAT) tableau des parametres scalaires du materiau
  29. C           a la temperature TINF
  30. C         XMASUP(NCOMAT) tableau des parametres scalaires du materiau
  31. C           a la temperature TSUP
  32. C           XMAT(1)=YOUNG ;XMAT(2)=XNU ;XMAT(3)=N
  33. C           XMAT(4)=M ;XMAT(5)=KK; XMAT(6)=ALF2;
  34. C           XMAT(7)=BET2 ;XMAT(8)=r; XMAT(9)=A ;
  35. C           XMAT(10)=q ;XMAT(11)=ALPHAT
  36. C           XMAT(12)=RHO; XMAT(13)=SIGY
  37. C         MFR indice de la formulation mecanique(seulement massif
  38. C           ou coque pour les materiaux endommageables
  39. C         ICARA nombre de caracteristiques geometriques des elements
  40. C           finis
  41. C         XCAR(ICARA) tableau des caracteristiques geometriques des
  42. C           elements finis
  43. C         IFOURB= -2 EN CONTR. PLANES
  44. C                 -1 EN DEFORM. PLANES
  45. C                  0 EN AXISYMETRIE
  46. C                  1 EN SERIE DE FOURIER
  47. C                  2 EN TRIDIM
  48. C         INDIC=0, 1 OU -1 pour plasticite avec endommagement
  49. C              =2 OU -2 pour viscoplasticite avec endommagement
  50. C         ITEST = 0 pas uniformite des listes de temperatures pour tous les
  51. C           coefficients non lineaires
  52. C           d'ou interpolation sur les coefficients du materiau
  53. C               = 1 uniformite des listes de temperatures pour tous les
  54. C           coefficients non lineaires
  55. C           d'ou moyenne ponderee sur la variable elle-meme
  56. C           ex: D(T)=A(T)*K**B(T)
  57. C               T=TETA*TINF+(1-TETA)*TSUP
  58. C           si ITEST=0 on remplace A(T) par (TETA*A(TINF)+(1-TETA))*A(TSUP)
  59. C                                et B(T) par (TETA*B(TINF)+(1-TETA))*B(TSUP)
  60. C           si ITEST=1 D(T)=TETA*D(TINF)+(1-TETA)*D(TSUP)
  61. C         ITHHER = 0 pas de chargement thermique et materiau constant
  62. C                = 1 chargement thermique et materiau constant
  63. C                = 2 chargement thermique et materiau(T)
  64. C------------------------------------------------------------------------
  65. C
  66. C------------------------------------------------------------------------
  67. C Sortie: EPSVD(NSTRS) derivee du tenseur des deformations
  68. C         VARD(NVARI) tableau contenant les derivees des variables
  69. C           internes
  70. C         VARF(NVARI) état des variables internes à la fin du sous programme
  71. C         SIGD(NSTRS) derivee des contraintes
  72. C         DD(NSTRS,NSTRS) matrice de Hooke  au debut du pas
  73. C         DDV(NSTRS,NSTRS) derivee de DD
  74. C         DDINV(NSTRS,NSTRS) inverse de DD
  75. C------------------------------------------------------------------------
  76.       IMPLICIT INTEGER(I-N)
  77.       IMPLICIT REAL*8 (A-H,O-Z)
  78.       DIMENSION XMAINF(*),XMASUP(*),XCAR(*)
  79.       DIMENSION SIG(*),EPSV(*),VAR(*)
  80.       DIMENSION SIGD(*),EPSVD(*),VARD(*),VARF(NVARI)
  81.       DIMENSION XMAT(*),DSPT(*),YOUNG(*),XNU(*),EPSTHD(*)
  82.       DIMENSION DDV(NSTRS,NSTRS),DDINVp(NSTRS,NSTRS)
  83.       DIMENSION DD(NSTRS,NSTRS),DDINV(NSTRS,NSTRS)
  84.       DIMENSION VARD1(3),EPSVD1(6),SIGD1(6),VARD2(3),EPSVD2(6)
  85.       DIMENSION EPSED(6),XX(6),Di(3),VaP(3),AA(3),BB(6)
  86.       DIMENSION ROT(3,3),ROTT(3,3),CO(3,3),COM(3,3)
  87.       DIMENSION D1(3,3),D2(3,3),DV1(3,3),DI1(3,3),DI2(3,3)
  88.       DIMENSION AAA(3,3),BBB(3,3),CCC(3,3),VcP(3,3),FTH(3,3)
  89.       PARAMETER (UN=1.D0)
  90.       LOGICAL endonul
  91. C
  92. C----- le materiau est constant ou pas, on derive normalement les contraintes
  93. C
  94.       IFLAG=0
  95.       IF ( (ITEST.EQ.1).AND.(ITHHER.EQ.2) ) THEN
  96. C----- si le materiau depend de la temperature et si les listes sont identiques
  97.            IF (TINF.NE.TSUP) THEN
  98. C-------------- Moyenne ponderee sur les vitesses des variables
  99.                 IFLAG=1
  100.                 TETA=(TI-TSUP)/(TINF-TSUP)
  101.                 CALL DERIV1(TAU,SIG,EPSV,VAR,SIGD1,EPSVD1,VARD1,EPSTHD,
  102.      &           DSPT,XMAINF,NSTRS,NVARI,NCOMAT,DD,DDV,DDINV,DDINVp,
  103.      &           YOUNG,NYOUNG,XNU,NXNU,MFR,XCAR,ICARA,IFOURB,INDIC,
  104.      &           TINF,TPOINT,ITHHER,VARF,IB,IGAU,KERRRE)
  105.                 IF(KERRE.NE.0) GO TO 100
  106.  
  107.                 DO 10 I=1,NVARI
  108. 10                     VARD2(I)=TETA*VARD1(I)
  109.                      DO 20 I=1,NSTRS
  110. 20                     EPSVD2(I)=TETA*EPSVD1(I)
  111.                 CALL DERIV1(TAU,SIG,EPSV,VAR,SIGD1,EPSVD1,VARD1,EPSTHD,
  112.      &           DSPT,XMASUP,NSTRS,NVARI,NCOMAT,DD,DDV,DDINV,DDINVp,
  113.      &           YOUNG,NYOUNG,XNU,NXNU,MFR,XCAR,ICARA,IFOURB,INDIC,
  114.      &           TSUP,TPOINT,ITHHER,VARF,IB,IGAU,KERRRE)
  115.                  IF(KERRE.NE.0) GO TO 100
  116.  
  117.                      DO 30 I=1,NVARI
  118. 30                     VARD2(I)=VARD2(I)+(1.D0-TETA)*VARD1(I)
  119.                      DO 40 I=1,NSTRS
  120. 40                     EPSVD2(I)=EPSVD2(I)+(1.D0-TETA)*EPSVD1(I)
  121.            ENDIF
  122.       ENDIF
  123. C
  124.       CALL DERIV1(TAU,SIG,EPSV,VAR,SIGD,EPSVD,VARD,EPSTHD,
  125.      &DSPT,XMAT,NSTRS,NVARI,NCOMAT,DD,DDV,DDINV,DDINVp,YOUNG,
  126.      &NYOUNG,XNU,NXNU,MFR,XCAR,ICARA,IFOURB,INDIC,TI,TPOINT,
  127.      &ITHHER,VARF,IB,IGAU,KERRE)
  128.       IF(KERRE.NE.0) GO TO 100
  129. C
  130.       IF (IFLAG.EQ.0) GOTO 100
  131. C
  132. C --- Cas ou ITEST = 1 et ITHHER = 2
  133. C
  134.                 DO I=1,NVARI
  135.                      VARD(I)=VARD2(I)
  136.                 ENDDO
  137.                 DO I=1,NSTRS
  138.                      EPSVD(I)=EPSVD2(I)
  139.                 ENDDO
  140. C
  141.        endonul = .FALSE.
  142.        IF ((VAR(2).EQ.0).AND.(VAR(3).EQ.0).AND.(VAR(4).EQ.0).AND.
  143.      &     (VAR(5).EQ.0).AND.(VAR(6).EQ.0).AND.(VAR(7).EQ.0).AND.
  144.      &     (VAR(8).EQ.0)) THEN
  145.          endonul = .TRUE.
  146.        ENDIF
  147. C
  148. C----------------------------------------------
  149. C TENSEUR DES CONTRAINTES (dans repere general)
  150. C----------------------------------------------
  151. C
  152.        CALL ZERO(CO,3,3)
  153.        IF (mfr.eq.1) THEN
  154.        DO I=1,3
  155.          CO(I,I) = SIG(I)
  156.        ENDDO
  157.        CO(1,2) = SIG(4)
  158.        IF(NSTRS.NE.6) GOTO 50
  159.        CO(1,3) = SIG(5)
  160.        CO(2,3) = SIG(6)
  161. 50     CONTINUE
  162.        CO(2,1) = CO(1,2)
  163.        CO(3,1) = CO(1,3)
  164.        CO(3,2) = CO(2,3)
  165.        ENDIF
  166. C
  167. C------------------------------------
  168. C TENSEUR D ENDOMMAGEMENT ET
  169. C TENSEUR D ENDOMMAGEMENT DIAGONALISE
  170. C------------------------------------
  171. C
  172.       XD11=VAR(2)
  173.       XD22=VAR(3)
  174.       XD33=VAR(4)
  175.       XD12=VAR(5)
  176.       XD13=VAR(6)
  177.       XD23=VAR(7)
  178.        CALL ZERO(D1,3,3)
  179.        D1(1,1) = XD11
  180.        D1(2,2) = XD22
  181.        D1(3,3) = XD33
  182.        D1(1,2) = XD12
  183.        D1(1,3) = XD13
  184.        D1(2,3) = XD23
  185.        D1(2,1) = D1(1,2)
  186.        D1(3,1) = D1(1,3)
  187.        D1(3,2) = D1(2,3)
  188. C
  189.        CALL JACOB3(D1,3,Di,ROT)
  190. C
  191.        CALL ZERO(D2,3,3)
  192.        D2(1,1) = Di(1)
  193.        D2(2,2) = Di(2)
  194.        D2(3,3) = Di(3)
  195.        w1 = 1.D0 - Di(1)
  196.        w2 = 1.D0 - Di(2)
  197.        w3 = 1.D0 - Di(3)
  198.        IF (D2(1,1).EQ.1.D0) w1=1.D-6
  199. C
  200.        ADMAX =  MAX(Di(1),Di(2),Di(3))
  201. C
  202. C-------------------------------------------------------
  203. C TENSEUR DES CONTRAINTES MOTRICES (dans repere general)
  204. C-------------------------------------------------------
  205. C
  206.        EQ=XMAT(10)
  207.        DO I=1,3
  208.          DO J=1,3
  209.            ROTT(I,J)=ROT(I,J)
  210.          ENDDO
  211.        ENDDO
  212.        CALL INVALM (ROTT, 3 , 3 , KER, 1.D-12)
  213. C
  214.       IF (endonul) THEN
  215.         DO I=1,3
  216.           DO J=1,3
  217.             COM(I,J) = CO(I,J)
  218.           ENDDO
  219.         ENDDO
  220.       ELSE
  221. C --- Matrice de [(I-D)/det(I-D)]**(q/2) dans le repère principal de D ---
  222.         CALL ZERO(DI1,3,3)
  223.         DO I=1,3
  224.           XNB=(1.D0-Di(I))/(w1*w2*w3)
  225.           IF (XNB.GE.0.D0) THEN
  226.             DI1(I,I)=(XNB)**(EQ/2.D0)
  227.           ELSE
  228.             DI1(I,I)=0.D0
  229.             WRITE(6,*) 'erreur'
  230.           ENDIF
  231.         ENDDO
  232. C --- passage dans le repère de général ---
  233.         CALL MULMAT(AAA,ROT,DI1,3,3,3)
  234.         CALL MULMAT(DI2,AAA,ROTT,3,3,3)
  235. C
  236.         CALL MULMAT(AAA,DI2,CO,3,3,3)
  237.         CALL MULMAT(COM,AAA,DI2,3,3,3)
  238.       ENDIF
  239. C
  240. C----------------------------------------------------------------------
  241. C
  242. C CALCUL DE  J1,J2
  243. C   J1: CONTRAINTE MOYENNE
  244. C   J2: CONTRAINTE DE VON MISES
  245. C
  246. C----------------------------------------------------------------------
  247. C
  248.        AJ1 = (SIG(1)+SIG(2)+SIG(3))/(3.D0)
  249.        DO I=1,3
  250.          XX(I) = SIG(I)-AJ1
  251.        ENDDO
  252.        DO I=4,NSTRS
  253.          XX(I) = SIG(I)
  254.        ENDDO
  255.        AJ2 =  PROCON(XX,XX,NSTRS)
  256.        AJ2 =  SQRT(1.5D0*AJ2)
  257. C
  258. C----------------------
  259. C FORCE THERMODYNAMIQUE
  260. C----------------------
  261. C
  262. C --- Etat de triaxialité de la contrainte ---
  263.        tri = (AJ2)/(SIG(1)+SIG(2)+SIG(3))
  264. C --- Valeurs propres et vecteurs propres du tenseur des contraintes motrices ---
  265.        CALL JACOB3(COM,3,VaP,VcP)
  266.        ER=XMAT(8)
  267.        EA=XMAT(9)
  268.        DO I=1,3
  269.          AA(I) = ((tri*VaP(I)/(EA))+ABS(tri*VaP(I)/(EA)))/(2.D0)
  270.          AA(I) = (AA(I))**(ER)
  271.        ENDDO
  272. C -- Dans le repère général ---
  273.        CALL ZERO(FTH,3,3)
  274.        FTH(1,1) = AA(1)*(VcP(1,1))**2+AA(2)*(VcP(1,2))**2
  275.        FTH(1,1) = FTH(1,1)+AA(3)*(VcP(1,3))**2
  276.        FTH(2,2) = AA(1)*(VcP(2,1))**2+AA(2)*(VcP(2,2))**2
  277.        FTH(2,2) = FTH(2,2)+AA(3)*(VcP(2,3))**2
  278.        FTH(3,3) = AA(1)*(VcP(3,1))**2+AA(2)*(VcP(3,2))**2
  279.        FTH(3,3) = FTH(3,3)+AA(3)*(VcP(3,3))**2
  280.        FTH(1,2) = AA(1)*VcP(1,1)*VcP(2,1)+AA(2)*VcP(1,2)*VcP(2,2)
  281.        FTH(1,2) = FTH(1,2)+AA(3)*VcP(1,3)*VcP(2,3)
  282.        FTH(1,3) = AA(1)*VcP(1,1)*VcP(3,1)+AA(2)*VcP(1,2)*VcP(3,2)
  283.        FTH(1,3) = FTH(1,3)+AA(3)*VcP(1,3)*VcP(3,3)
  284.        FTH(2,3) = AA(1)*VcP(2,1)*VcP(3,1)+AA(2)*VcP(2,2)*VcP(3,2)
  285.        FTH(2,3) = FTH(1,3)+AA(3)*VcP(2,3)*VcP(3,3)
  286.        FTH(3,1) = FTH(1,3)
  287.        FTH(2,1) = FTH(1,2)
  288.        FTH(3,2) = FTH(2,3)
  289. C
  290.        BET2=XMAT(7)
  291.        DV1(1,1) = (BET2*FTH(1,1))+FTH(2,2)+FTH(3,3)
  292.        DV1(2,2) = (BET2*FTH(2,2))+FTH(1,1)+FTH(3,3)
  293.        DV1(3,3) = (BET2*FTH(3,3))+FTH(2,2)+FTH(1,1)
  294.        DV1(1,2) = (BET2-1.D0)*FTH(1,2)
  295.        DV1(1,3) = (BET2-1.D0)*FTH(1,3)
  296.        DV1(2,3) = (BET2-1.D0)*FTH(2,3)
  297.        DV1(2,1) = DV1(1,2)
  298.        DV1(3,1) = DV1(1,3)
  299.        DV1(3,2) = DV1(2,3)
  300. C
  301. C----------------------------------------------------------------------
  302. C
  303. C CALCUL DE LA VARIATION DES CONTRAINTES
  304. C
  305. C----------------------------------------------------------------------
  306. C
  307.       CALL DERTRA(NYOUNG,YOUNG,TI,YUNG,YUNGV,TINF,TSUP)
  308.       CALL DERTRA(NXNU,XNU,TI,ENU,ENUV,TINF,TSUP)
  309. C
  310. C --- Calcul de l inverse de la matrice de Hooke au debut du pas ---
  311.       CALL ELAST4(2,IFOURB,VAR,NVARI,XMAT,NCOMAT,YUNGV,ENUV,
  312.      &XCAR,ICARA,MFR,NSTRS,DDINV,DDV,KERRE,INDIC,ITHHER)
  313.       IF(KERRE.NE.0) GO TO 100
  314. C
  315. C --- Calcul de la matrice de Hooke et de sa derivee au debut du pas ---
  316.       CALL ELAST4(1,IFOURB,VAR,NVARI,XMAT,NCOMAT,YUNGV,ENUV,
  317.      &XCAR,ICARA,MFR,NSTRS,DD,DDV,KERRE,INDIC,ITHHER)
  318.       IF(KERRE.NE.0) GO TO 100
  319. C
  320. C --- Détermination de la variation de la matrive de Hooke inversée ---
  321.       CALL HOOKP(VAR,VARD,NVARI,XMAT,NCOMAT,MFR,NSTRS,DDINVp)
  322. C
  323. C --- Calcul de la variation de contrainte ---
  324.       CALL ZDANUL(SIGD,NSTRS)
  325.       DO I=1,3
  326.         XX(I) = CO(I,I)
  327.       ENDDO
  328.       XX(4) = CO(1,2)
  329.       XX(5) = CO(1,3)
  330.       XX(6) = CO(2,3)
  331. C
  332.       CALL ZDANUL(EPSED,NSTRS)
  333.       DO I=1,NSTRS
  334.         EPSED(I)=DSPT(I)-EPSVD(I)-EPSTHD(I)
  335.       ENDDO
  336. C
  337.       CALL MULMAT(BB,DDINVp,XX,NSTRS,1,NSTRS)
  338. C
  339.       DO I=1,NSTRS
  340.         XX(I)=EPSED(I)-BB(I)
  341.       ENDDO
  342. C
  343.       CALL MULMAT(SIGD,DD,XX,NSTRS,1,NSTRS)
  344. C
  345.  100  CONTINUE
  346.       RETURN
  347.       END
  348.  
  349.  
  350.  
  351.  
  352.  
  353.  
  354.  
  355.  
  356.  
  357.  

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