Télécharger la plaquette Cast3M


fr en

Télécharger labord.eso

Retour à la liste

Numérotation des lignes : 

labord

  1. C LABORD    SOURCE    CB215821  17/11/30    21:16:43     9639           
  2. C     ==================================================================
  3. C     =       SUBROUTINE INTEGRANT LE MODELE D'ENDOMMAGEMENT           =
  4. C     =     UNILATERAL EN UNIAXIAL AVEC PILOTAGE EN DEFORMATION        =
  5. C     =              FG + FR (LMT-Cachan, fevrier 2001)                =
  6. C     ==================================================================
  7.       SUBROUTINE LABORD(XMAT,DEPS,SIG0,VAR0,SIGF,VARF)
  8.  
  9.       IMPLICIT INTEGER(I-N)
  10.       IMPLICIT REAL*8 (a-h,o-z)
  11.  
  12.       REAL*8 XMAT(13),DEPS(3),SIG0(3),VAR0(8),SIGF(3),VARF(8)
  13.  
  14. C---------------------------------------------------------------
  15. C     MODELE DE CLB (CLB_UNI) SPECIALISE POUR LE MODELE A FIBRE
  16. C---------------------------------------------------------------
  17. C
  18. C     variables en entree
  19. C
  20. C     XMAT : Caracteristique du materiaux
  21. C
  22. C     EPS : Deformation totale
  23. C
  24. C     DEPS : Increment de deformations
  25. C
  26. C     SIG0 : Contraintes initiales
  27. C
  28. C     VAR0 : Variables internes initiales
  29. C
  30. C
  31. C     variables en sortie
  32. C
  33. C     SIGF : Contraintes finales
  34. C
  35. C     VARF : Variables internes finales
  36. C
  37. C     declaration des variables
  38. C
  39.       logical change
  40.  
  41. C     ==================================================================
  42. C     =     Initialisation des Constantes                              =
  43. C     ==================================================================
  44.       E=XMAT(1)
  45.       XNU=XMAT(2)
  46.       Y01=XMAT(5)
  47.       Y02=XMAT(6)
  48.       A1=XMAT(7)
  49.       A2=XMAT(8)
  50.       B1=XMAT(9)
  51.       B2=XMAT(10)
  52.       BETA1=XMAT(11)
  53.       BETA2=XMAT(12)
  54.       SIGF1=XMAT(13)
  55. C     ==================================================================
  56. C     =     Initialisation des VARIABLES HISTORIQUES                   =
  57. C     ==================================================================
  58.       ESEC = VAR0(1)
  59.       D1 = VAR0(2)
  60.       D2 = VAR0(3)
  61.       YLIM1 = VAR0(4)
  62.       YLIM2 = VAR0(5)
  63.       EPS10 = VAR0(6)
  64. C
  65.       EPS = EPS10 + DEPS(1)
  66.       y1=YLIM1-1.D0
  67.       y2=YLIM2-1.D0
  68. C     ==================================================================
  69. ***
  70.       change=.true.
  71.       do while (change)
  72. ***
  73.       X1=(BETA1*D1/(1.D0-D1)+BETA2*D2/(1.D0-D2))/E
  74.       X2=(BETA2*D2-SIGF1)/(1.D0-D2)/E
  75.       if ( X1 .LE. Eps ) then
  76.         EPSF=Eps-BETA2*D2/(1.0D0-D2)/E
  77.         CALL CAS13(EPSF,SIG,ESEC,E,D1,BETA1,A1,B1,Y1,YLIM1,Y01)
  78.         icas1=1
  79.       else if ( Eps .LE. X2 ) then
  80.         CALL CAS13(Eps,SIG,ESEC,E,D2,BETA2,A2,B2,Y2,YLIM2,Y02)
  81.         icas1=3
  82.       else
  83.         call CAS2(Eps,BETA1,E,D1,SIGF1,BETA2,D2,SIG,ESEC)
  84.         icas1=2
  85.       endif
  86.       X1=(BETA1*D1/(1.D0-D1)+BETA2*D2/(1.D0-D2))/E
  87.       X2=(BETA2*D2-SIGF1)/(1.D0-D2)/E
  88.       if ( X1 .LE. Eps ) then
  89.         icas2=1
  90.       else if ( Eps .LE. X2 ) then
  91.         icas2=3
  92.       else
  93.         icas2=2
  94.       endif
  95.       if (icas1 .eq. icas2) then
  96.          change=.false.
  97.       endif
  98. ***
  99.       enddo
  100. ***
  101. c     Si on garde le module du beton constant
  102.       ESEC=E
  103.  
  104. C     ==================================================================
  105. C     =     reactualisation des seuils                                 =
  106. C     ==================================================================
  107.       YLIM1= MAX(Y1,YLIM1)
  108.       YLIM2= MAX(Y2,YLIM2)
  109.  
  110. C     ==================================================================
  111. C     =     reactualisation des VARIABLES HISTORIQUES                  =
  112. C     ==================================================================
  113.       VARF(1) = ESEC
  114.       VARF(2) = D1
  115.       VARF(3) = D2
  116.       VARF(4) = YLIM1
  117.       VARF(5) = YLIM2
  118.       VARF(6) = EPS
  119.       VARF(7) = VAR0(7) + DEPS(2)
  120.       VARF(8) = VAR0(8) + DEPS(3)
  121. C     ==================================================================
  122. C     =     Calcul des contraintes finales                             =
  123. C     ==================================================================
  124.       SIGF(1) = SIG
  125.       r_z = 0.5D0 * E / (1.D0+XNU)
  126.       SIGF(2) = r_z * VARF(7)
  127.       SIGF(3) = r_z * VARF(8)
  128.  
  129.       RETURN
  130.       END
  131.  
  132. C     ==================================================================
  133. C     =               CALCUL DE SIGMA DANS LE CAS 1 & 3                =
  134. C     =             ( le calcul se fait par dichotomie )               =
  135. C     ==================================================================
  136.       SUBROUTINE CAS13(EPS,SIG,ESEC,E,D,BETA,A,B,Y,YLIM,Y00)
  137.  
  138.       IMPLICIT INTEGER(I-N)
  139.       IMPLICIT REAL*8(A-H,O-Z)
  140.  
  141.       REAL*8 m0,m
  142.  
  143.       dt = 1.0d0
  144.       m0 = 0.0d0
  145.       alpha = 0.0d0
  146.  
  147. c     Erreur relative par rapport a D
  148.       Erreur_max=1.0D-6
  149.       Nbpas=50
  150.  
  151. c     Dans le cas No 1 ou 3, on recherche D13 dans [0,1]
  152.       X1= d-0.01D0
  153.       X2= 1.0D0
  154.  
  155.       pas=(X2-X1)/Nbpas
  156.       m=m0
  157.       if ( EPS .ge. 0.D0) then
  158.          m=m0/650.D0
  159.       endif
  160.  
  161. C     ------------------------------------------------------------------
  162. C     -   Avant de faire les calculs, on verifie que l'on depasse le   -
  163. C     -   seuil YLIM                                                   -
  164. C     ------------------------------------------------------------------
  165.       D13=D
  166.       y0=y00
  167.       YY = (E*EPS+BETA)*(E*EPS+BETA)-BETA*BETA/(1.0D0-D13)/(1.0D0-D13)
  168.       YY = YY/(2.0D0*E)
  169.  
  170. ***   Prise en compte de la viscosite
  171.  
  172. ****  coefficient de viscosite
  173. c      m=3.d4
  174. c      alpha=0.18d0
  175. ****
  176. **      Si yy est inferieur a y0 on sort sans evolution de d
  177.       If ( YY .LE. Y0 ) Then
  178.         y=yy
  179. c       Module Secant
  180.         ESEC=E*(1.0D0-D)
  181. c       Valeur de la contrainte
  182.         SIG=EPS*ESEC-BETA*D
  183.         return
  184.       endif
  185.  
  186. ***********
  187. **      Calcul du seuil
  188. c      y=(E*EPS+BETA)*(E*EPS+BETA)/(2.D0*E)
  189.  
  190. **      Calcul de l endommagement
  191.       Erreur=1.D0
  192.       XG=d
  193.       Xd=X1
  194.       Do while ( (Erreur .GT. Erreur_max) .and. (pas .ge. 1.d-8) )
  195.  
  196. c       Valeur de la fonction a la borne de droite D13=XD
  197.         D13=XD
  198.         YY = (E*EPS+BETA)*(E*EPS+BETA)/(2.0D0*E)
  199.      &      -(BETA*BETA/(1.D0-D13)/(1.D0-D13))/(2.D0*E)
  200. c        if (yy .lt. ylim) then
  201. c          yy=ylim
  202. c        endif
  203.         y0=y00+m*( (abs(d13-d)+(d13-d))/(2.D0*dt) )**alpha
  204.         If ( YY .GT. Y0 ) Then
  205. **      si Dconverge est > a D13 alors fctd est < 0
  206.           Fctd =  1.0D0 - (1.0D0+(A*(YY-Y0))**B)*(1.0D0-D13)
  207.         Else
  208.           Fctd =  1.d0
  209.         Endif
  210.  
  211.         If ( Fctd .lt. 0.D0) then
  212.           Xg=Xd
  213.           Xd=Xd+pas
  214.         Else
  215.           Xd=Xg
  216.           Pas=Pas * 0.5D0
  217.           Xd=Xd+pas
  218.         Endif
  219.  
  220.         If ( Xg .gt. X2 ) Then
  221.           Stop ' . DANS CAS13 . Pas de Zero.'
  222.         Endif
  223.  
  224.         Erreur=abs(Fctd)
  225.       Enddo
  226.  
  227. c     On verifie que d13 > d . Si ce n'est pas le cas rien n'evolue.
  228. c      print*, (abs(d13-d)+(d13-d))/(2.D0*dt)
  229. ***
  230.       y=y00+1.D0/a*(d13/(1.D0-d13))**(1.D0/b)
  231. ***
  232.       If ( d13 .Ge. d ) then
  233.         D = D13
  234.       Endif
  235.       if (d .gt. 0.9999999D0) then
  236.          d=0.9999999D0
  237.       endif
  238.  
  239. c     Module Secant
  240.       ESEC=E*(1.0D0-D)
  241. c     Valeur de la contrainte
  242.       SIG=EPS*ESEC-BETA*D
  243.  
  244.       Return
  245.       End
  246.  
  247. C     ==================================================================
  248. C     =                 CALCUL DE SIGMA DANS LE CAS 2                  =
  249. C     =             ( le calcul se fait par dichotomie )               =
  250. C     ==================================================================
  251.       SUBROUTINE CAS2(EPS,BETA1,E,D1,SIGF1,BETA2,D2,SIG,ESEC)
  252.  
  253.       IMPLICIT INTEGER(I-N)
  254.       IMPLICIT REAL*8(A-H,O-Z)
  255.  
  256.       fprim=  (e*eps*(1.D0-d2)-beta2*d2+sigf1)
  257.      &      / (sigf1+beta1*d1*(1.D0-d2)/(1.D0-d1))
  258.  
  259. c     La solution est la borne de droite.
  260.       SIG=e*eps*(1.D0-d2)-beta2*d2-beta1*d1*(1.D0-d2)*fprim/(1.D0-d1)
  261.  
  262. c     ESEC est le module du beton de compression
  263.       ESEC= E*(1.0D0-D2)
  264.  
  265.       Return
  266.       End
  267.  
  268. C=======================================================================
  269. C REMARQUE : LES SOUS-PROGRAMMES RAIDTANi (avec i=1,2,3) NE SONT
  270. C            PAS UTILISES !
  271. C=======================================================================
  272.       SUBROUTINE RAIDTAN1 (EPS,SIG,ETAN,E,D1,BETA1,D2,BETA2,A,B,Y,Y0)
  273.  
  274.       IMPLICIT INTEGER(I-N)
  275.       IMPLICIT REAL*8(A-H,O-Z)
  276.  
  277. **      derive de d par rapport a y
  278.         dddy= b*(a**b)*(y-y0)**(b-1.D0)/(1.D0+(a*(y-y0))**b)**2
  279.  
  280. **      derive de d par rapport a epsilon
  281.         dddeps=dddy*( e*eps+beta1-beta2*d2/(1.D0-d2) )
  282.         dddeps=dddeps/( 1.D0+dddy*beta1**2/(e*(1.D0-d1)**3) )
  283.  
  284. **      derive de sigma par rapport a epsilon -- calcul de E tangent
  285.         etan=e*(1.D0-d1)-e*eps*dddeps-beta1*dddeps+beta2*d2*dddeps/
  286.      $  (1.D0-d2)
  287.  
  288.         return
  289.         end
  290. C---------
  291.       SUBROUTINE RAIDTAN2 (EPS,BETA1,E,D1,SIGF1,BETA2,D2,SIG,ETAN)
  292.  
  293.       IMPLICIT INTEGER(I-N)
  294.       IMPLICIT REAL*8(A-H,O-Z)
  295.  
  296.         etan=E*(1.D0-d2)/( sigf1+beta1*d1*(1.D0-d2)/(1.D0-d1) )
  297.         etan=etan*beta1*d1*(1.D0-d2)/(1.D0-d1)
  298.         etan=e*(1.D0-d2)-etan
  299.  
  300.         return
  301.         end
  302. C---------
  303.       SUBROUTINE RAIDTAN3 (EPS,SIG,ETAN,E,D1,BETA1,D2,BETA2,A,B,Y,Y0)
  304.  
  305.       IMPLICIT INTEGER(I-N)
  306.       IMPLICIT REAL*8(A-H,O-Z)
  307.  
  308. **      derive de d par rapport a y
  309.         dddy= b*(a**b)*(y-y0)**(b-1.D0)/(1.D0+(a*(y-y0))**b)**2
  310.  
  311. **      derive de d par rapport a epsilon
  312.         dddeps=dddy*(e*eps+beta2)
  313.         dddeps=dddeps/(1.D0+dddy*beta2**2/(e*(1.D0-d2)**3))
  314.  
  315. **      derive de sigma par rapport a epsilon -- calcul de E tangent
  316.         etan=e*(1.D0-d2)-e*eps*dddeps-beta2*dddeps
  317.  
  318.         return
  319.         end
  320.  
  321. C     ==================================================================
  322. C     =                 CALCUL DE SIGMA DANS LE CAS 2                  =
  323. C     =             ( le calcul se fait par dichotomie )               =
  324. C     ==================================================================
  325.       SUBROUTINE CAS21(EPS,BETA1,E,D1,SIGF1,BETA2,D2,SIG,ESEC)
  326.  
  327.       IMPLICIT INTEGER(I-N)
  328.       IMPLICIT REAL*8(A-H,O-Z)
  329.  
  330. c     Erreur relative par rapport a Sigma
  331.       Erreur_max=1.0D-6
  332.       Nbpas=50
  333.       itermax=1000
  334.       iter=0
  335. c     Dans le cas No 2, on recherche Y dans [-1,0]   Y=sig/sigf1
  336.       X1= -1.0D0
  337.       X2=  0.0D0
  338.  
  339.       pas=(X2-X1)/Nbpas
  340.  
  341. c     Valeur initiale des bornes de gauche et droite
  342.       XG=X1
  343.       XD=XG+PAS
  344.  
  345. c     Valeur de la fonction a la borne de gauche y=XG
  346.       Y=XG
  347.       Fctg    =  EPS*E - BETA1*D1*(1.0D0-3.0D0*Y*Y-2.0D0*Y*Y*Y)/
  348.      $  (1.0D0-D1)
  349.      &        - BETA2*D2/(1.0D0-D2) - Y*SIGF1/(1.0D0-D1)
  350.      &        - SIGF1*(1.0D0/(1.0D0-D1)-1.0D0/(1.0D0-D2))*(2.0D0+Y)*Y*Y
  351.       If (Fctg .GE. 0.0D0) then
  352.         isigg=1
  353.       else
  354.         isigg=0
  355.       endif
  356.       Erreur=abs(Fctg)
  357. C
  358.       Do while ( (Erreur .GT. Erreur_max) .and. (iter .le. itermax) )
  359.          iter=iter+1
  360. c       Valeur de la fonction a la borne de droite y=XD
  361.         Y=XD
  362.         Fctd    =  EPS*E - BETA1*D1*(1.0D0-3.0D0*Y*Y-2.0D0*Y*Y*Y)/
  363.      &   (1.0D0-D1)    - BETA2*D2/(1.0D0-D2) - Y*SIGF1/(1.0D0-D1)
  364.      &       - SIGF1*(1.0D0/(1.0D0-D1)-1.0D0/(1.0D0-D2))*(2.0D0+Y)*Y*Y
  365.         If (Fctd .GE. 0.0D0) then
  366.           isigd=1
  367.         else
  368.           isigd=0
  369.         endif
  370.         Erreur=abs(Fctd)
  371.         If ( isigg .EQ. isigd) then
  372.           Fctg=Fctd
  373.           isigg=isigd
  374.           Xg=Xd
  375.           Xd=Xd+pas
  376.         Else
  377.           Pas=Pas * 0.5D0
  378.           Xd=Xg+pas
  379.         Endif
  380.         If ( Xg .gt. X2 ) Then
  381.           Stop ' . DANS CAS2 . Pas de Zero.'
  382.         Endif
  383.       Enddo
  384. c     La solution est la borne de droite.
  385.       SIG=Y*Sigf1
  386. c     Si on fait varier lineairement ESEC en fonction de Sig
  387.       ESEC= E*(1.0D0 - D1 + Y*(D2-D1))
  388.       Return
  389.       End
  390.  
  391.  
  392.  
  393.  
  394.  

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