Télécharger fibmaz.eso

Retour à la liste

Numérotation des lignes :

  1. C FIBMAZ SOURCE CHAT 12/04/06 21:15:11 7348
  2. SUBROUTINE FIBMAZ (XMAT,DEPS,SIG0,VAR0,SIGF,VARF)
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8(A-H,O-Z)
  5. DIMENSION XMAT(*),DEPS(3),SIG0(3),VAR0(5),SIGF(3),VARF(5)
  6. C---------------------------------------------------------------
  7. C MODELE DE MAZARS (MAZARS_FIB) SPECIALISE POUR LE MODELE A FIBRE
  8. C Didier Combescure et Pierre Pegon Nov. 93
  9. C---------------------------------------------------------------
  10. C
  11. C
  12. C variables en entree
  13. C
  14. C XMAT : Caracteristique du materiaux
  15. C
  16. C EPS : Deformations initiales
  17. C
  18. C DEPS : Increment de deformations
  19. C
  20. C SIG0 : Contraintes initiales
  21. C
  22. C VAR0 : Variables internes initiales
  23. C
  24. C
  25. C variables en sortie
  26. C
  27. C SIGF : Contraintes finales
  28. C
  29. C VARF : Variables internes finales
  30. C
  31. C declaration des variables
  32. C
  33. C
  34. REAL*8 EPS33(3,3),EPSIPP(3),EPSILC(3),VALP33(3,3)
  35. REAL*8 SIGP(3),SIGPC(3)
  36. PARAMETER (XZERO=0.D0 , UN=1.D0 , DEUX=2.D0, XPETIT=1.D-12)
  37.  
  38. C
  39. C
  40. C recuperation des variables initiales dans les tableaux
  41. C
  42. C
  43. YOUN = XMAT(1)
  44. XNU = XMAT(2)
  45. EPSD0= XMAT(5)
  46. ACOM = XMAT(6)
  47. BCOM = XMAT(7)
  48. ATRA = XMAT(8)
  49. BTRA = XMAT(9)
  50. BETA = XMAT(10)
  51. DINI = VAR0(2)
  52. EPSX = VAR0(3)
  53. EPSY = VAR0(4)
  54. EPSZ = VAR0(5)
  55.  
  56. C
  57. C calcul des seuils
  58. C
  59. C
  60. C calcul de la deformation totale
  61. C
  62. DO 100 ISTRS=1,3
  63. VARF(ISTRS+2)=VAR0(ISTRS+2)+DEPS(ISTRS)
  64. 100 CONTINUE
  65. EPS33(1,1) = VARF(3)
  66. EPS33(1,2) = (VARF(4)/2.D0)
  67. EPS33(2,1) = EPS33(1,2)
  68. EPS33(1,3) = (VARF(5)/2.D0)
  69. EPS33(3,1) = EPS33(1,3)
  70. EPS33(2,2) = ((-1.D0*XNU)*VARF(3))
  71. EPS33(3,3) = EPS33(2,2)
  72. EPS33(3,2) = XZERO
  73. EPS33(2,3) = XZERO
  74.  
  75. C
  76. C calcul des deformations principales
  77. C
  78. C on diagonalise
  79. C
  80. CALL JACOB3 (EPS33,3,EPSIPP,VALP33)
  81.  
  82. C
  83. C on calcule les contraintes ppales
  84. C
  85. C
  86. SIGP(1)=((UN-XNU)*EPSIPP(1)+
  87. 1 (XNU*EPSIPP(2))+
  88. 2 (XNU*EPSIPP(3)))
  89. SIGP(2)=((UN-XNU)*EPSIPP(2)+
  90. 1 (XNU*EPSIPP(1))+
  91. 2 (XNU*EPSIPP(3)))
  92. SIGP(3)=((UN-XNU)*EPSIPP(3)+
  93. 1 (XNU*EPSIPP(2))+
  94. 2 (XNU*EPSIPP(1)))
  95. DO 200 ISTRS=1,3
  96. SIGP(ISTRS)=(SIGP(ISTRS)*YOUN/((UN+XNU)*(UN-2.D0*XNU)))
  97. 200 CONTINUE
  98.  
  99. C
  100. C on calcule le epsilontild
  101. C
  102. EPSTIL=MAX( XZERO , EPSIPP(1) )**2 +
  103. 1 MAX( XZERO , EPSIPP(2) )**2 +
  104. 2 MAX( XZERO , EPSIPP(3) )**2
  105. EPSTIL=SQRT (EPSTIL)
  106. C
  107. C on est toujours dans le cas istep=0
  108. C
  109. EPSTIM=EPSTIL
  110. VARF(1)=EPSTIL
  111. C
  112. C on calcule l'endommagement et les contraintes
  113. C
  114. C
  115. C on calcule ALFAT ALFAC DT et DC puis D
  116. C dans le cas ou le seuil initial est depasse
  117. C
  118. C on calcule le signe des contraintes elastiques
  119. C
  120. DO 300 ISTRS=1,3
  121. SIGPC(ISTRS)=MIN(XZERO,SIGP(ISTRS))
  122. 300 CONTINUE
  123. TRSIGC=SIGPC(1)+SIGPC(2)+SIGPC(3)
  124. C
  125. C ainsi que GAMMA et EPSTIM si BETA est different de zero
  126. C
  127. IF (BETA.GT.1.D-5) THEN
  128. IF (ABS(TRSIGC).GT.XPETIT) THEN
  129. GAMMA=SIGPC(1)*SIGPC(1)+SIGPC(2)*SIGPC(2)+
  130. 1 SIGPC(3)*SIGPC(3)
  131. GAMMA= SQRT(GAMMA)/ABS(TRSIGC)
  132. EPSTIM=EPSTIM*GAMMA
  133. END IF
  134. END IF
  135. C
  136. IF ( EPSTIM .GT. EPSD0) THEN
  137. C
  138. C calcul de l'endommagement
  139. C
  140. C on calcule les deformations dues aux contraintes negatives
  141. C
  142. DO 400 ISTRS=1,3
  143. EPSILC(ISTRS)=(SIGPC(ISTRS)*(UN+XNU)-TRSIGC*XNU)/YOUN
  144. 400 CONTINUE
  145. C
  146. C on en deduit ALFAT et ALFAC
  147. C
  148. ALFAC = MAX(XZERO,EPSIPP(1))*EPSILC(1) +
  149. 1 MAX(XZERO,EPSIPP(2))*EPSILC(2) +
  150. 2 MAX(XZERO,EPSIPP(3))*EPSILC(3)
  151. ALFAC = ALFAC/(EPSTIL*EPSTIL)
  152. ALFAT = UN - ALFAC
  153. C
  154. C amelioration de la reponse en cisaillement pour beta > 1.
  155. C
  156. IF (BETA .GT. UN) THEN
  157. IF ( ALFAT .GT. XPETIT ) THEN
  158. ALFAT=ALFAT**BETA
  159. END IF
  160. IF ( ALFAC .GT. XPETIT ) THEN
  161. ALFAC=ALFAC**BETA
  162. END IF
  163. END IF
  164.  
  165. C
  166. C on calcule DT et DC puis D
  167. C dans le cas ou le seuil initial est depasse
  168. C
  169. DT=UN - EPSD0*(UN-ATRA)/EPSTIM -
  170. 1 ATRA/EXP(BTRA*(EPSTIM-EPSD0))
  171. DC=UN - EPSD0*(UN-ACOM)/EPSTIM -
  172. 1 ACOM/EXP(BCOM*(EPSTIM-EPSD0))
  173. D = ALFAT*DT + ALFAC*DC
  174. C
  175. C on borne la valeur de D a 0.99999
  176. C
  177. D=MIN ( D , UN-1.D-05 )
  178. ELSE
  179. D=XZERO
  180. END IF
  181. C
  182. C on teste la croissance de D
  183. C
  184. D=MAX ( D , DINI )
  185. C
  186. C on le stocke dans les variables internes finales
  187. C
  188. VARF(2)= D
  189. C
  190. C on calcule les contraintes finales
  191. C
  192. SIGF(1) = YOUN*VARF(3)
  193. SIGF(2) = YOUN*VARF(4)/2.D0/(UN+XNU)
  194. SIGF(3) = YOUN*VARF(5)/2.D0/(UN+XNU)
  195. DO 500 ISTRS=1,3
  196. SIGF(ISTRS)=SIGF(ISTRS)*(UN-D)
  197. 500 CONTINUE
  198.  
  199. RETURN
  200. END
  201.  
  202.  
  203.  
  204.  

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