Télécharger labord.eso

Retour à la liste

Numérotation des lignes :

  1. C LABORD SOURCE CHAT 12/04/06 21:15:21 7348
  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=DMAX1(Y1,YLIM1)
  108. YLIM2=DMAX1(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.  

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