Télécharger la plaquette Cast3M


fr en

Télécharger dynar.eso

Retour à la liste

Numérotation des lignes : 

dynar

  1. C DYNAR     SOURCE    PV        21/01/14    21:15:00     10847          
  2.       SUBROUTINE DYNAR(WRK52,WRK53,WRK54,IECOU)
  3. c
  4. C=======================================================================
  5. C
  6. C        MODELE VISCOPLASTIQUE VISCOENDOMMAGEABLE POUR LA DYNAMIQUE
  7. C                    - PROGRAMME PRINCIPAL -
  8. C                          VERSION 2.0
  9. C
  10. C=======================================================================
  11. C     CREATION :F.GATUINGT
  12. C=======================================================================
  13. C=========== Hypopthese: au temps t=0 tout est nul   ===================
  14. C=======================================================================
  15.         IMPLICIT REAL*8(A-H,O-Z)
  16.  
  17. -INC PPARAM
  18. -INC CCOPTIO
  19. -INC DECHE
  20.         DIMENSION EPSPTDT(6),EPSTE(6),EPSPOPT(6),EPSPO(6),SIGPO(6),
  21.      &     sigt(6),sigtdt(6),EPST(6),EPSTDT(6),EPSPT(6),EPSTDTE(6),
  22.      &   S(3,3),EPSIPP(3),ARRAY(3,3),SIGTDTE(6),DFIDSIG(6),EPSPOPTDT(6)
  23. C
  24.         PARAMETER (XUnDemi=0.5D0,XZERO=0.D0,UN=1.D0,DEUX=2.D0)
  25.         PARAMETER (XPETIT=1.D-12)
  26.         PARAMETER (TROIS=3.D0)
  27. C
  28.       SEGMENT IECOU
  29. *      COMMON/IECOU/NYOG,NYNU,NYALFA,NYSMAX,NYN,NYM,NYKK,
  30.        INTEGER icow1,icow2,icow3,icow4,icow5,icow6,icow7,
  31. C      INTEGER NYOG, NYNU, NYALFA,NYSMAX,NYN, NYM,  NYKK,
  32.      1  icow8,icow9,icow10,icow11,icow12,icow13,icow14,icow15,icow16,
  33. C    .  NYALF1,NYBET1,NYR,  NYA,   NYRHO,NSIGY,  NNKX,  NYKX,  IND,
  34.      2 icow17,icow18,icow19,icow20,icow21,icow22,icow23,icow24,
  35. C    .  NSOM,  NINV, NINCMA,NCOMP, JELEM, LEGAUS,INAT,  NCXMAT,
  36.      3 icow25,icow26,icow27,icow28,icow29,icow30,ICARA,
  37. C    .  LTRAC, MFR,  IELE,  NHRM,   NBNN, NBELEM,ICARA,
  38.      4  icow32,icow33,NSTRS1,MFR1,icow36,icow37,icow38,
  39. C    .  LW2,   NDEF,  NSTRSS,MFR1,  NBGMAT,NELMAT,MSOUPA,
  40.      5  icow39,icow40,icow41,icow42,icow43,icow44
  41. C    .  NUMAT1,LENDO, NBBB,  NNVARI,KERR1, MELEME
  42.         INTEGER icow45,icow46,icow47,icow48,icow49,icow50,
  43.      .          icow51,icow52,icow53,icow54,icow55,icow56
  44.      .          icow57,icow58
  45.       ENDSEGMENT
  46. C
  47. C
  48. C
  49. C
  50. C - Recuperation des donnees materiau
  51. C
  52.         Eini=XMAT(1)
  53.         XNU=XMAT(2)
  54.         D0=XMAT(5)
  55.         Q1=XMAT(6)
  56.         Q2=XMAT(7)
  57.         Q3=XMAT(8)
  58.         SIGM0=XMAT(9)
  59.         XN=XMAT(10)
  60.         XKPLAS = XMAT(12)
  61.         XMPLAS = XMAT(11)
  62.         XK=XMAT(13)
  63.         XMDt=XMAT(14)
  64.         XNDt=XMAT(15)
  65.         At=XMAT(16)
  66.         Bt=XMAT(17)
  67.         XMDc=XMAT(18)
  68.         XNDc=XMAT(19)
  69.         Ac=XMAT(20)
  70.         Bc=XMAT(21)
  71.         ED0=XMAT(22)
  72. C - Pas de temps
  73.         DELTAT=TEMPF-TEMP0
  74.         deltat = max (1d-30,deltat)
  75. C - Contrainte au temps T
  76.         SIGT(1)=SIG0(1)
  77.         SIGT(2)=SIG0(2)
  78.         SIGT(3)=SIG0(3)
  79.         SIGT(4)=SIG0(4)
  80.         SIGT(5)=SIG0(5)
  81.         SIGT(6)=SIG0(6)
  82. C - Deformation au temps T
  83.         EPST(1)=EPST0(1)
  84.         EPST(2)=EPST0(2)
  85.         EPST(3)=EPST0(3)
  86.         EPST(4)=EPST0(4)/2.
  87.         EPST(5)=EPST0(5)/2.
  88.         EPST(6)=EPST0(6)/2.
  89. C - Increment de deformation
  90.         DEPST(1)=DEPST(1)
  91.         DEPST(2)=DEPST(2)
  92.         DEPST(3)=DEPST(3)
  93.         DEPST(4)=DEPST(4)/2.
  94.         DEPST(5)=DEPST(5)/2.
  95.         DEPST(6)=DEPST(6)/2.
  96. C======================================================================
  97. C - Subcycling si le pas de temps est trop grand=======================
  98. C======================================================================
  99.         IF (DELTAT.GT.1.D-3) THEN
  100.           coef=100
  101.         ELSE
  102.           coef=max (1.d0 , deltat/ 1.d-5)
  103.         ENDIF
  104.         nbiter=int(coef)
  105.         DO i=1,6
  106.         DEPST(i)=DEPST(i)/nbiter
  107.         ENDDO
  108.         DELTAT=DELTAT/nbiter
  109.         l=0
  110.         DO WHILE (l.LT.nbiter)
  111. C - Debut des boucles
  112. C=======================================================================
  113. C - Recuperation des variables historiques
  114. C=======================================================================
  115.         DT=VAR0(1)
  116.         fT=VAR0(2)+D0
  117.         EPSPT(1)=VAR0(3)
  118.         EPSPT(2)=VAR0(4)
  119.         EPSPT(3)=VAR0(5)
  120.         EPSPT(4)=VAR0(6)
  121.         EPSPT(5)=VAR0(7)
  122.         EPSPT(6)=VAR0(8)
  123.         SIGMT=VAR0(9)+SIGM0
  124.         dtt=VAR0(10)
  125.         dct=VAR0(11)
  126.         ES=VAR0(12)+ED0
  127. C======================================================================
  128. C - Deformation au temps t+dt
  129. C======================================================================
  130.         EPSTDT(1)=EPST(1)+DEPST(1)
  131.         EPSTDT(2)=EPST(2)+DEPST(2)
  132.         EPSTDT(3)=EPST(3)+DEPST(3)
  133.         EPSTDT(4)=EPST(4)+DEPST(4)
  134.         EPSTDT(5)=EPST(5)+DEPST(5)
  135.         EPSTDT(6)=EPST(6)+DEPST(6)
  136.         TRACEPSTDT=EPSTDT(1)+EPSTDT(2)+EPSTDT(3)
  137. C======================================================================
  138. C - Vitesse de Deformation
  139. C======================================================================
  140.         DO i=1,6
  141.           EPSPO(i)=DEPST(i)/DELTAT
  142.         ENDDO
  143.         TRACEPSPO=EPSPO(1)+EPSPO(2)+EPSPO(3)
  144. C=======================================================================
  145. C - Module d'elasticite tangent de la Matrice
  146. C=======================================================================
  147.         ETAN=(Eini)/XN*(SIGMt/SIGM0)**(1-XN)
  148. C=======================================================================
  149. C - Coefficients de compressibilite et cisaillement de la matrice sans les pores
  150. C=======================================================================
  151.         XKm=Eini/(TROIS*(UN-DEUX*XNU))
  152.         XGm=Eini/(DEUX*(UN+XNU))
  153. C=======================================================================
  154. C - Coeficients de compressibilite et cisaillement de la matrice avec les pores (Mori-Tanaka)
  155. C=======================================================================
  156.         XKporo=4.D0*XKm*XGm*(UN-fT)/(4.D0*XGm+3.D0*XKm*fT)
  157.         XGporo=XGm*(UN-fT)/(UN+fT*
  158.      &  (6.D0*XKm+12.D0*XGm)/(9.D0*XKm+8.D0*XGm))
  159. C ======================================================================
  160. C - Coefficients de compressibilite et cisaillement de la matrice endommagee
  161. C ======================================================================
  162.         XKendo=(UN-Dt)*XKporo
  163.         XGendo=(UN-Dt)*XGporo
  164. C ======================================================================
  165. C Calcul des invariants de la contrainte
  166. C ======================================================================
  167.         XI1t=SIGT(1)+SIGT(2)+SIGT(3)
  168.         XJ2t=0.5D0*((SIGT(1)-SIGT(2))**2
  169.      &  +(SIGT(2)-SIGT(3))**2+(SIGT(3)-SIGT(1))**2
  170.      &  +6.D0*(SIGT(4)**2+SIGT(5)**2+SIGT(6)**2))
  171. C======================================================================
  172. C - Deformation elastique au temps t+dt
  173. C======================================================================
  174. C Calcul du multiplicateur plastique au temps t
  175.         Fit=XJ2t/(SIGMT**2)+2*Q1*fT*COSH((Q2*XI1t)/(2*SIGMT))
  176.      &     -(1+(Q3*fT)**2)
  177.         IF (Fit.LT.0.) THEN
  178.           Fit=0.
  179.         ENDIF
  180.         xlambdapo=(fT/(1-fT))*(Fit/XKPLAS)**XMPLAS
  181.  
  182. C - Calcul des deformations plastique au temps t
  183. C Derivees partielles de fi
  184.         DFIDSIG(1)=0.5*(4*SIGT(1)-2*SIGT(2)-2*SIGT(3))/SIGMT**2
  185.      &  +Q1*Q2*ft*SINH(Q2*XI1t/(2*SIGMt))/SIGMt
  186.         DFIDSIG(2)=0.5*(-2*SIGT(1)+4*SIGT(2)-2*SIGT(3))/SIGMT**2
  187.      &  +Q1*Q2*ft*SINH(Q2*XI1t/(2*SIGMt))/SIGMt
  188.         DFIDSIG(3)=0.5*(-2*SIGT(1)-2*SIGT(2)+4*SIGT(3))/SIGMT**2
  189.      &  +Q1*Q2*ft*SINH(Q2*XI1t/(2*SIGMt))/SIGMt
  190.         DFIDSIG(4)=6/(SIGMt**2)*SIGT(4)
  191.         DFIDSIG(5)=6/(SIGMt**2)*SIGT(5)
  192.         DFIDSIG(6)=6/(SIGMt**2)*SIGT(6)
  193. C
  194.         DO i=1,6
  195.           EPSPOPT(i)=xlambdapo*DFIDSIG(i)
  196.         ENDDO
  197.         TRACEPSPOPT=EPSPOPT(1)+EPSPOPT(2)+EPSPOPT(3)
  198. C - Déformation viscoplastique
  199.         DO i=1,6
  200.           EPSPTDT(i)=EPSPT(i)+EPSPOPT(i)*DELTAT
  201.         ENDDO
  202.         TRACEPSPTDT=EPSPTDT(1)+EPSPTDT(2)+EPSPTDT(3)
  203. C - Deformation elastique
  204.         DO i=1,6
  205.           EPSTDTE(i)=EPSTDT(i)-EPSPTDT(i)
  206.         ENDDO
  207.         TRACEPSTDTE=EPSTDTE(1)+EPSTDTE(2)+EPSTDTE(3)
  208. C======================================================================
  209. C - Deformation equivalente elastique au sens de Mazars a t
  210.         ARRAY(1,1)=EPSTDTE(1)
  211.         ARRAY(2,2)=EPSTDTE(2)
  212.         ARRAY(3,3)=EPSTDTE(3)
  213.         ARRAY(1,2)=EPSTDTE(4)
  214.         ARRAY(2,3)=EPSTDTE(5)
  215.         ARRAY(1,3)=EPSTDTE(6)
  216.         ARRAY(3,1)=ARRAY(1,3)
  217.         ARRAY(3,2)=ARRAY(2,3)
  218.         ARRAY(2,1)=ARRAY(1,2)
  219. C
  220.         CALL JACOB3(ARRAY,3,EPSIPP,S)
  221.         EPSE=MAX( XZERO , EPSIPP(1) )**2 +
  222.      1       MAX( XZERO , EPSIPP(2) )**2 +
  223.      2       MAX( XZERO , EPSIPP(3) )**2
  224.         EPSE=SQRT(EPSE)
  225. C - Nouveau seuil de Mazars
  226.         ES=max (ES, EPSE)
  227. C ==============================================================
  228. C - Prediction élastique de la contrainte (pas d'évolution de d)
  229.         DO i=1,3
  230.           SIGPO(i)=XKendo*(TRACEPSPO-TRACEPSPOPT)+2*XGendo*
  231.      &    ((EPSPO(i)-1./3.*TRACEPSPO)-(EPSPOPT(i)-1./3.*TRACEPSPOPT))
  232.         ENDDO
  233.         DO i=4,6
  234.           SIGPO(i)=2*XGendo*((EPSPO(i))
  235.      &     -(EPSPOPT(i)))
  236.         ENDDO
  237. C
  238.         DO i=1,6
  239.           SIGTDTE(i)=SIGT(i)+SIGPO(i)*DELTAT
  240.         ENDDO
  241. C ============================================================
  242. C - Calcul des parametres avec la contrainte predite elastique
  243. C - Calcul des invariants à t
  244.         XI1tdte=SIGTDTE(1)+SIGTDTE(2)+SIGTDTE(3)
  245. C
  246.         XJ2tdte=0.5*((SIGTDTE(1)-SIGTDTE(2))**2
  247.      &  +(SIGTDTE(2)-SIGTDTE(3))**2+(SIGTDTE(3)-SIGTDTE(1))**2
  248.      &  +6*(SIGTDTE(4)**2+SIGTDTE(5)**2+SIGTDTE(6)**2))
  249. C
  250. C Derivees partielles de fi
  251.         DFIDSIG(1)=0.5*(4*SIGTDTE(1)
  252.      &  -2*SIGTDTE(2)-2*SIGTDTE(3))/SIGMT**2
  253.      &  +Q1*Q2*ft*SINH(Q2*XI1tdte/(2*SIGMt))/SIGMt
  254.         DFIDSIG(2)=0.5*(-2*SIGTDTE(1)+4*SIGTDTE(2)
  255.      &  -2*SIGTDTE(3))/SIGMT**2
  256.      &  +Q1*Q2*ft*SINH(Q2*XI1tdte/(2*SIGMt))/SIGMt
  257.         DFIDSIG(3)=0.5*(-2*SIGTDTE(1)-2*SIGTDTE(2)
  258.      &  +4*SIGTDTE(3))/SIGMT**2
  259.      &  +Q1*Q2*ft*SINH(Q2*XI1tdte/(2*SIGMt))/SIGMt
  260.         DFIDSIG(4)=6/(SIGMt**2)*SIGTDTE(4)
  261.         DFIDSIG(5)=6/(SIGMt**2)*SIGTDTE(5)
  262.         DFIDSIG(6)=6/(SIGMt**2)*SIGTDTE(6)
  263. C
  264.         TRACEDFIDSIG=DFIDSIG(1)+DFIDSIG(2)+DFIDSIG(3)
  265. C
  266.         DFIDSIGM=-2*XJ2tdte/(SIGMt**3)-Q1*Q2*ft*XI1tdte
  267.      &  *SINH(Q2*XI1tdte/(2*SIGMt))/(SIGMt**2)
  268. C
  269.         DFIDft=2*Q1*COSH(Q2*XI1tdte/(2*SIGMt))-2*ft*Q3*Q3
  270. C===================================================================
  271. C - Calcul de la surface seuil fi à l'aide de la contrainte prédite
  272.         Fi=XJ2tdte/(SIGMT**2)+2*Q1*fT*COSH((Q2*XI1tdte)/(2*SIGMT))
  273.      &     -(1+(Q3*fT)**2)
  274. C
  275. C======================================================================
  276. C========= Organigramme du calcul de la contrainte a t+Dt =============
  277. C======================================================================
  278. C
  279. C Test sur le signe de la deformation equivalente
  280.  
  281.  
  282.         IF(EPSE.GT.XZERO) THEN
  283.           IF(fi.LT.XZERO)THEN
  284.             IF(EPSE.LT.ES) THEN
  285. C =====================================================================
  286. C - Elasticite Seule
  287.               dtdt=dt
  288.               ftdt=ft
  289.               dttdt=dtt
  290.               dctdt=dct
  291.               sigmtdt=sigmt
  292.               DO i=1,6
  293.                 sigtdt(i)=sigtdte(i)
  294.               ENDDO
  295.             ELSE
  296. C =====================================================================
  297. C - Visco-Endommagement Seul
  298.               sigmtdt=sigmt
  299.               ftdt=ft
  300.               CALL VISCOE(XMAT,EPSIPP,EPSE,DTT,DCT,DT
  301.      &                    ,DELTAT,DTTDT,DCTDT,DTDT)
  302. C
  303.               Dpo=(Dtdt-Dt)/DELTAT
  304. C
  305. C - Raideur endommagee et ses derivees
  306.               XKendo=(1-Dtdt)*XKporo
  307.               XGendo=(1-Dtdt)*XGporo
  308.               XKpo=-Dpo*XKporo
  309.               XGpo=-Dpo*XGporo
  310. C
  311. C - Calcul de la contrainte
  312.               DO i=1,3
  313.                SIGPO(i)=XKendo*(TRACEPSPO-TRACEPSPOPT)+2*XGendo
  314.      &         *((EPSPO(i)-1./3.*TRACEPSPO)-(EPSPOPT(i)
  315.      &         -1./3.*TRACEPSPOPT))
  316.      &         +XKpo*(TRACEPSTDT-TRACEPSPTDT)+2*XGpo
  317.      &          *((EPSTDT(i)-1./3.*TRACEPSTDT)-(EPSPTDT(i)
  318.      &         -1./3.*TRACEPSPTDT))
  319.               ENDDO
  320.               DO i=4,6
  321.                 SIGPO(i)=2*XGendo*((EPSPO(i))
  322.      &          -(EPSPOPT(i)))+2*XGpo
  323.      &        *((EPSTDT(i))-(EPSPTDT(i)))
  324.               ENDDO
  325. C
  326.               DO i=1,6
  327.                 SIGTDT(i)=SIGT(i)+SIGPO(i)*DELTAT
  328.               ENDDO
  329.            ENDIF
  330.          ELSE
  331.            IF(EPSE.LT.ES) THEN
  332. C =====================================================================
  333. C - Plasticite sans endommagement
  334.              Dtdt=Dt
  335.              dttdt=dtt
  336.              dctdt=dct
  337. C - Evolution de la contrainte dans la matrice
  338.              EPSPOPM=-xlambdapo*DFIDSIGM
  339.              SIGPOM=(Etan*Eini/(Eini-Etan))*EPSPOPM
  340.              SIGMTDT=SIGMT+SIGPOM*DELTAT
  341. C - Evolution de la porosite
  342.              fpo=XK*((UN-ft)*ft*TRACEPSPOPT)
  343.              ftdt=ft+fpo*DELTAT
  344.              IF (ftdt.LT.XPETIT) THEN
  345.                 ftdt=XZERO
  346.              ENDIF
  347. C - Raideur endommagee et ses derives
  348.              XKporo=4*XKm*XGm*(1-fTDT)/(4*XGm+3*XKm*fTDT)
  349.              XGporo=XGm*(1-fTDT)/(1+fTDT*(6*XKm+12*XGm)/(9*XKm+8*XGm))
  350.              XKendo=(1-Dtdt)*XKporo
  351.              XGendo=(1-Dtdt)*XGporo
  352.              XKpo=(1-Dtdt)*(-4*XKm*XGm*(4*XGm+3*XKm)
  353.      &        /(4*XGm+3*XKm*fTDT)**2)*fpo
  354.              XGpo=(1-Dtdt)*(-5*XGm*(32*XGm**2+60*XGm*XKm+27*XKm**2)
  355.      &        /(4*(2+3*fTDT)*XGm+3*(3+2*fTDT)*XKm)**2)*fpo
  356. C =============================
  357. C - Correction de la contrainte
  358.              DO i=1,3
  359.                SIGPO(i)=XKendo*(TRACEPSPO-TRACEPSPOPT)+2*XGendo
  360.      &         *((EPSPO(i)-1./3.*TRACEPSPO)-(EPSPOPT(i)
  361.      &         -1./3.*TRACEPSPOPT))
  362.      &         +XKpo*(TRACEPSTDT-TRACEPSPTDT)+2*XGpo*
  363.      &         ((EPSTDT(i)-1./3.*TRACEPSTDT)-(EPSPTDT(i)
  364.      &          -1./3.*TRACEPSPTDT))
  365.              ENDDO
  366.              DO i=4,6
  367.                SIGPO(i)=2*XGendo*((EPSPO(i))
  368.      &          -(EPSPOPT(i)))+2*XGpo*
  369.      &         ((EPSTDT(i))-(EPSPTDT(i)))
  370.              ENDDO
  371. C
  372.              DO i=1,6
  373.                 SIGTDT(i)=SIGT(i)+SIGPO(i)*DELTAT
  374.              ENDDO
  375.            ELSE
  376. C =====================================================================
  377. C - Plasticite avec endommagement
  378. C - Evolution de l'endommagement
  379.               CALL VISCOE(XMAT,EPSIPP,EPSE,DTT,DCT,DT
  380.      &                    ,DELTAT,DTTDT,DCTDT,DTDT)
  381. C
  382.              Dpo=(Dtdt-Dt)/DELTAT
  383.              IF (Dpo.LE.0.) THEN
  384.                Dtdt=Dt
  385.                Dpo=0.
  386.              ENDIF
  387. C - Evolution de la contrainte dans la matrice
  388.              EPSPOPM=-xlambdapo*DFIDSIGM
  389.              SIGPOM=(Etan*Eini/(Eini-Etan))*EPSPOPM
  390.              SIGMTDT=SIGMT+SIGPOM*DELTAT
  391. C - Evolution de la porosite
  392.              fpo=XK*((1-ft)*ft*TRACEPSPOPT)
  393.              ftdt=ft+fpo*DELTAT
  394.              IF (ftdt.LT.0.) THEN
  395.                 ftdt=0.
  396.              ENDIF
  397. C - Raideur endommagee et ses derives
  398.              XKporo=4*XKm*XGm*(1-fTDT)/(4*XGm+3*XKm*fTDT)
  399.              XGporo=XGm*(1-fTDT)/(1+fTDT*(6*XKm+12*XGm)/(9*XKm+8*XGm))
  400.              XKendo=(1-dtdt)*XKporo
  401.              XGendo=(1-dtdt)*XGporo
  402.              XKpo=-Dpo*XKporo+(1-Dtdt)*
  403.      &        (-4*XKm*XGm*(4*XGm+3*XKm)/(4*XGm+3*XKm*fTDT)**2)*fpo
  404.              XGpo=-Dpo*XGporo
  405.      &        +(1-Dtdt)*(-5*XGm*(32*XGm**2+60*XGm*XKm+27*XKm**2)
  406.      &        /(4*(2+3*fTDT)*XGm+3*(3+2*fTDT)*XKm)**2)*fpo
  407. C =============================
  408. C - Correction de la contrainte
  409.              DO i=1,3
  410.                SIGPO(i)=XKendo*(TRACEPSPO-TRACEPSPOPT)+2*XGendo
  411.      &         *((EPSPO(i)-1./3.*TRACEPSPO)-(EPSPOPT(i)
  412.      &          -1./3.*TRACEPSPOPT))
  413.      &         +XKpo*(TRACEPSTDT-TRACEPSPTDT)+2*XGpo*
  414.      &        ((EPSTDT(i)-1./3.*TRACEPSTDT)-(EPSPTDT(i)
  415.      &         -1./3.*TRACEPSPTDT))
  416.              ENDDO
  417.              DO i=4,6
  418.                SIGPO(i)=2*XGendo*((EPSPO(i))
  419.      &           -(EPSPOPT(i)))+2*XGpo*
  420.      &        ((EPSTDT(i))-(EPSPTDT(i)))
  421.              ENDDO
  422. C
  423.              DO i=1,6
  424.                 SIGTDT(i)=SIGT(i)+SIGPO(i)*DELTAT
  425.              ENDDO
  426.            ENDIF
  427.          ENDIF
  428.         ELSE
  429.          IF(fi.LT.XZERO)THEN
  430. C =====================================================================
  431. C - Elasticite Seule
  432.            Dtdt=Dt
  433.            ftdt=ft
  434.            dttdt=dtt
  435.            dctdt=dct
  436.            sigmtdt=sigmt
  437.            DO i=1,6
  438.              sigtdt(i)=sigtdte(i)
  439.            ENDDO
  440.          ELSE
  441. C =====================================================================
  442. C - Plasticite sans endommagement
  443.            Dtdt=Dt
  444.            dttdt=dtt
  445.            dctdt=dct
  446. C - Evolution de la contrainte dans la matrice
  447.            EPSPOPM=-xlambdapo*DFIDSIGM
  448.            SIGPOM=(Etan*Eini/(Eini-Etan))*EPSPOPM
  449.            SIGMTDT=SIGMT+SIGPOM*DELTAT
  450. C - Evolution de la porosite
  451.            fpo=XK*((1-ft)*ft*TRACEPSPOPT)
  452.            ftdt=ft+fpo*DELTAT
  453.            IF (ftdt.LT.0.) THEN
  454.              ftdt=0.
  455.            ENDIF
  456. C - Raideur endommagee et ses derives
  457.            XKporo=4*XKm*XGm*(1-fTDT)/(4*XGm+3*XKm*fTDT)
  458.            XGporo=XGm*(1-fTDT)/(1+fTDT*(6*XKm+12*XGm)/(9*XKm+8*XGm))
  459.            XKendo=(1-Dtdt)*XKporo
  460.            XGendo=(1-Dtdt)*XGporo
  461.            XKpo=(1-Dtdt)*(-4*XKm*XGm*(4*XGm+3*XKm)
  462.      &      /(4*XGm+3*XKm*fTDT)**2)*fpo
  463.            XGpo=(1-Dtdt)*(-5*XGm*(32*XGm**2+60*XGm*XKm+27*XKm**2)
  464.      &      /(4*(2+3*fTDT)*XGm+3*(3+2*fTDT)*XKm)**2)*fpo
  465. C =============================
  466. C - Correction de la contrainte
  467.            DO i=1,3
  468.              SIGPO(i)=XKendo*(TRACEPSPO-TRACEPSPOPT)+2*XGendo
  469.      &       *((EPSPO(i)-1./3.*TRACEPSPO)-(EPSPOPT(i)
  470.      &       -1./3.*TRACEPSPOPT))
  471.      &       +XKpo*(TRACEPSTDT-TRACEPSPTDT)+2*XGpo*
  472.      &      ((EPSTDT(i)-1./3.*TRACEPSTDT)-(EPSPTDT(i)
  473.      &       -1./3.*TRACEPSPTDT))
  474.            ENDDO
  475.            DO i=4,6
  476.              SIGPO(i)=2*XGendo*((EPSPO(i))
  477.      &        -(EPSPOPT(i)))+2*XGpo*
  478.      &      ((EPSTDT(i))-(EPSPTDT(i)))
  479.            ENDDO
  480. C
  481.            DO i=1,6
  482.              SIGTDT(i)=SIGT(i)+SIGPO(i)*DELTAT
  483.            ENDDO
  484.          ENDIF
  485.         ENDIF
  486. C
  487. C======================================================================
  488. C========= Fin de l'organigramme du calcul de la contrainte a t+dt ====
  489. C======================================================================
  490. C
  491. C =====================================================================
  492. C - Evolution des variables historiques
  493.         VAR0(1)=DTDT
  494.         VAR0(2)=fTDT-D0
  495.         VAR0(3)=EPSPTDT(1)
  496.         VAR0(4)=EPSPTDT(2)
  497.         VAR0(5)=EPSPTDT(3)
  498.         VAR0(6)=EPSPTDT(4)
  499.         VAR0(7)=EPSPTDT(5)
  500.         VAR0(8)=EPSPTDT(6)
  501.         VAR0(9)=SIGMTDT-SIGM0
  502.         VAR0(10)=dttdt
  503.         VAR0(11)=dctdt
  504.         VAR0(12)=ES-ED0
  505. C =====================================================================
  506.         SIGT(1)=SIGTDT(1)
  507.         SIGT(2)=SIGTDT(2)
  508.         SIGT(3)=SIGTDT(3)
  509.         SIGT(4)=SIGTDT(4)
  510.         SIGT(5)=SIGTDT(5)
  511.         SIGT(6)=SIGTDT(6)
  512. C =====================================================================
  513.         EPST(1)=EPSTDT(1)
  514.         EPST(2)=EPSTDT(2)
  515.         EPST(3)=EPSTDT(3)
  516.         EPST(4)=EPSTDT(4)
  517.         EPST(5)=EPSTDT(5)
  518.         EPST(6)=EPSTDT(6)
  519. C =====================================================================
  520. C  - Compteur du Subcycling
  521.         l=l+1
  522.         ENDDO
  523. C =====================================================================
  524. C ========Fin du Subcycling=================
  525. C =====================================================================
  526. C - Evolution des variables historiques
  527.         VARF(1)=DTDT
  528.         VARF(2)=fTDT-D0
  529.         VARF(3)=EPSPTDT(1)
  530.         VARF(4)=EPSPTDT(2)
  531.         VARF(5)=EPSPTDT(3)
  532.         VARF(6)=EPSPTDT(4)
  533.         VARF(7)=EPSPTDT(5)
  534.         VARF(8)=EPSPTDT(6)
  535.         VARF(9)=SIGMTDT-SIGM0
  536.         VARF(10)=dttdt
  537.         VARF(11)=dctdt
  538.         VARF(12)=ES-ED0
  539. C======================================================================
  540. C - Mise a jour des contraintes
  541.         SIGF(1)=SIGTDT(1)
  542.         SIGF(2)=SIGTDT(2)
  543.         SIGF(3)=SIGTDT(3)
  544.         SIGF(4)=SIGTDT(4)
  545.         SIGF(5)=SIGTDT(5)
  546.         SIGF(6)=SIGTDT(6)
  547. C
  548.         END
  549.  
  550.  
  551.  
  552.  
  553.  
  554.  
  555.  
  556.  
  557.  
  558.  
  559.  
  560.  

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