Télécharger flufib.eso

Retour à la liste

Numérotation des lignes :

  1. C FLUFIB SOURCE FD218221 20/01/10 21:15:02 10494
  2. SUBROUTINE FLUFIB(IMOD,PROPS,DEPS,SIG0,VAR0,SIGF,VARF,TIME0,TIMEF,
  3. & NPROPS,NVAR)
  4. C-----------------------------------------------------------------------
  5. C RESOLUTION DES MODELES DE FLUAGE POUR LES POUTRES A FIBRE
  6. C - Implementation des lois de fluage pour les modeles de poutre
  7. C a fibre (modeles de section)
  8. C - Les modeles sont ecrits en dimension 1, dans la direction des
  9. C fibres de la poutre. Ainsi :
  10. C -- la contrainte equivalente est definie comme la valeur absolue
  11. C de la contrainte axiale (composante XX)
  12. C -- la vitesse de deformation equivalente de fluage est definie
  13. C comme la valeur absolue de la vitesse de deformation de
  14. C fluage axiale (composante XX)
  15. C - Le schema d'integration est de type Runge-Kutta a l'ordre 4
  16. C avec controle de l'erreur et adaptation du pas de temps
  17. C - Liste des modeles de fluage actuelement disponibles:
  18. C -- Norton (IMOD=15)
  19. C -- Polynomial (IMOD=16)
  20. C -- Blackburn (IMOD=17)
  21. C -- Blackburn 2 (IMOD=18)
  22. C -- Lemaitre (IMOD=19)
  23. C Francois Di Paola (nov 2019)
  24. C-----------------------------------------------------------------------
  25. C
  26. C Entrees
  27. C -------
  28. C IMOD : Numero du modele a fibre
  29. C PROPS : Caracteristiques du materiaux
  30. C NPROPS : Taille du tableau PROPS
  31. C DEPS : Increment de deformations totales sur le pas de temps (modele section)
  32. C SIG0 : Contraintes au debut du pas de temps (modele section)
  33. C VAR0 : Variables internes au debut du pas de temps
  34. C NVAR : Taille du tableau VAR0 (et aussi VARF)
  35. C TIME0 : Instant initial (debut du pas de temps)
  36. C TIMEF : Instant final (fin du pas de temps)
  37. C
  38. C Sorties
  39. C -------
  40. C SIGF : Contraintes a la fin du pas de temps (modele section)
  41. C VARF : Variables internes a la fin du pas de temps
  42. C
  43. C
  44. C Detail des composantes de materiau
  45. C ----------------------------------
  46. C PROPS(1) : 'YOUN' Module de Young
  47. C PROPS(2) : 'NU ' Coefficient de Poisson
  48. C PROPS(3) : 'RHO ' Masse volumique
  49. C PROPS(4) : 'ALPH' Coefficient de dilatation thermique
  50. C PROPS(.) : puis les autres composantes, selon le modele
  51. C
  52. C Detail des composantes des variables
  53. C ------------------------------------
  54. C DEPS(1) = D Epsilon_xx SIG0(1) SIGF(1) = Sig_xx
  55. C DEPS(2) = D Gamma_xy SIG0(2) SIGF(2) = Sig_xy
  56. C DEPS(3) = D Gamma_xz SIG0(3) SIGF(3) = Sig_xz
  57. C
  58. C
  59. C-----------------------------------------------------------------------
  60. C
  61. C Quelques declarations
  62. IMPLICIT INTEGER(I-N)
  63. IMPLICIT REAL*8(A-H,O-Z)
  64. PARAMETER(XZER=0.D0,XPETIT=1.D-6,
  65. & UNDEMI=0.5D0,UN=1.D0,DEUX=2.D0,TROIS=3.D0,SIX=6.D0)
  66. C
  67. C Declaration de quelques tableaux d'entree/sortie de la subroutine
  68. REAL*8 PROPS(NPROPS),DEPS(3),SIG0(3),VAR0(NVAR),SIGF(3),VARF(NVAR)
  69. C
  70. LOGICAL DTLIBR
  71. C
  72. C=======================================================================
  73. C 1 - Recuperation de l'etat mecanique au debut du pas
  74. C=======================================================================
  75. C Pas de temps
  76. DT=TIMEF-TIME0
  77. C Contrainte au debut du pas
  78. SIG=SIG0(1)
  79. C Variables internes au debut du pas
  80. IF (IMOD.EQ.19) THEN
  81. XXIN=VAR0(1)
  82. EPSIN=VAR0(2)
  83. ELSE
  84. EPSIN=VAR0(1)
  85. ENDIF
  86. C Increment de deformation totale sur le pas
  87. DEPST=DEPS(1)
  88. C Recuperation des parametres materiau
  89. YOU=PROPS(1)
  90. XNU=PROPS(2)
  91. G=YOU/(DEUX*(UN+XNU))
  92. C
  93. C=======================================================================
  94. C 2 - Initialisations pour le schema de Runge-Kutta 4
  95. C=======================================================================
  96. PRECIS=1.D-5
  97. DTLIBR=.TRUE.
  98. IF (DT.LT.0.0) THEN
  99. PRINT*,'ERREUR DT < 0 !'
  100. CALL ERREUR(414)
  101. RETURN
  102. ENDIF
  103. IF (DT.EQ.0.0) DT=1.D-20
  104. VSIG= YOU*DEPST/DT
  105. XMM1=XM-UN
  106. DTLEFT=DT
  107. BORNE=2.0
  108. RMAX=1.3
  109. RMIN=0.7
  110. DIV=7.0
  111. FAC=3.0
  112. MSOUPA=10000
  113. XMAX=YOU*1.D-3
  114. C
  115. C=======================================================================
  116. C 3 - Initialisations avant iterations en sous-increments
  117. C=======================================================================
  118. ITERO=0
  119. 10 CONTINUE
  120. ITERO=1+ITERO
  121. IF (ITERO.NE.1) THEN
  122. DTLIBR=.TRUE.
  123. PRECIS=PRECIS*7.D0
  124. IF (ITERO.GT.3) THEN
  125. PRINT*,'ERREUR : INCREMENT DE CHARGEMENT TROP GRAND !'
  126. CALL ERREUR(268)
  127. RETURN
  128. ENDIF
  129. ENDIF
  130. DTLEFT=DT
  131. TAU=DTLEFT
  132. ASIG=ABS(SIG)
  133. ERRABS=PRECIS*ASIG
  134. IF (XMAX.GT.ASIG) ERRABS=PRECIS*XMAX
  135. C
  136. C=======================================================================
  137. C 4 - Iterations en sous increments / fin si DTLEFT = 0
  138. C=======================================================================
  139. NSSINC=0
  140. NITERA=0
  141. 20 CONTINUE
  142. NSSINC=NSSINC+1
  143. IF (NSSINC.GT.MSOUPA) THEN
  144. DTLIBR=.FALSE.
  145. GOTO 10
  146. ENDIF
  147. C [1.1] Vitesses de contrainte et de deformation inelastique a partir de l'etat debut du pas
  148. SEQ=ABS(SIG)
  149. IF (SEQ.LT.XPETIT) THEN
  150. VEPSIN1=XZER
  151. IF (IMOD.EQ.19) VXX1=XZER
  152. ELSE
  153. C Modele de Norton
  154. IF(IMOD.EQ.15) THEN
  155. CALL NORTON(SEQ,EPSIN,TAU,PROPS(6),PROPS(7),PROPS(8),VEPSIN1)
  156. C Modele Polynomial
  157. ELSEIF(IMOD.EQ.16) THEN
  158. VEPSIN1=PROPS(6)+PROPS(7)*(SEQ**PROPS(8))+
  159. & PROPS(9)*(SEQ**PROPS(10)) +
  160. & PROPS(11)*(SEQ**PROPS(12))
  161. C Modele de Blackburn
  162. ELSEIF(IMOD.EQ.17) THEN
  163. CALL SOLUTN(20,SEQ,EPSIN,PROPS,DX,TIMEX)
  164. IF(TIMEX.LT.TAU) TIMEX=ABS(TAU)
  165. CALL EQUATN(20,SEQ,TIMEX,XZER,PROPS,VEPSINA,VEPSINB)
  166. VEPSIN1=VEPSINA+VEPSINB
  167. C Modele de Blackburn_2
  168. ELSEIF(IMOD.EQ.18) THEN
  169. CALL SOLUTN(61,SEQ,EPSIN,PROPS,DX,TIMEX)
  170. IF(TIMEX.LT.TAU) TIMEX=ABS(TAU)
  171. CALL EQUATN(61,SEQ,TIMEX,XZER,PROPS,VEPSINA,VEPSINB)
  172. VEPSIN1=VEPSINA+VEPSINB
  173. C Modele de Lemaitre
  174. ELSEIF(IMOD.EQ.19) THEN
  175. CALL LEMAIT(SEQ,XXIN,VXX1,TAU,PROPS,VEPSIN1)
  176. ENDIF
  177. ENDIF
  178. VSIG1=DSIGN(YOU*VEPSIN1,SIG)
  179. C Debut des iterations sur tau optimal / fin si RA
  180. NITERA=0
  181. 30 CONTINUE
  182. NITERA=NITERA+1
  183. TAU2=TAU*UNDEMI
  184. C [1.2] 1ere estimation a mi-pas, avec la vitesse debut de pas
  185. SIG1=SIG+(TAU2*(VSIG-VSIG1))
  186. EPSIN1=EPSIN+TAU2*VEPSIN1
  187. IF (IMOD.EQ.19) XX1=XXIN+TAU2*VXX1
  188. C [2.1] Vitesses de contrainte et de deformation inelastique a partir de l'etat a mi-pas [1.2]
  189. SEQ1=ABS(SIG1)
  190. IF (SEQ1.LT.XPETIT) THEN
  191. VEPSIN2=XZER
  192. IF (IMOD.EQ.19) VXX2=XZER
  193. ELSE
  194. C Modele de Norton
  195. IF(IMOD.EQ.15) THEN
  196. CALL NORTON(SEQ1,EPSIN1,TAU2,PROPS(6),PROPS(7),PROPS(8),
  197. & VEPSIN2)
  198. C Modele Polynomial
  199. ELSEIF(IMOD.EQ.16) THEN
  200. VEPSIN2=PROPS(6)+PROPS(7)*(SEQ1**PROPS(8))+
  201. & PROPS(9)*(SEQ1**PROPS(10)) +
  202. & PROPS(11)*(SEQ1**PROPS(12))
  203. C Modele de Blackburn
  204. ELSEIF(IMOD.EQ.17) THEN
  205. CALL SOLUTN(20,SEQ1,EPSIN1,PROPS,DX,TIMEX)
  206. IF(TIMEX.LT.TAU2) TIMEX=ABS(TAU2)
  207. CALL EQUATN(20,SEQ1,TIMEX,XZER,PROPS,VEPSINA,VEPSINB)
  208. VEPSIN2=VEPSINA+VEPSINB
  209. C Modele de Blackburn_2
  210. ELSEIF(IMOD.EQ.18) THEN
  211. CALL SOLUTN(61,SEQ1,EPSIN1,PROPS,DX,TIMEX)
  212. IF(TIMEX.LT.TAU2) TIMEX=ABS(TAU2)
  213. CALL EQUATN(61,SEQ1,TIMEX,XZER,PROPS,VEPSINA,VEPSINB)
  214. VEPSIN2=VEPSINA+VEPSINB
  215. C Modele de Lemaitre
  216. ELSEIF(IMOD.EQ.19) THEN
  217. CALL LEMAIT(SEQ1,XX1,VXX2,TAU2,PROPS,VEPSIN2)
  218. ENDIF
  219. ENDIF
  220. VSIG2=DSIGN(YOU*VEPSIN2,SIG1)
  221. C [2.2] 2nde estimation a mi-pas, avec la vitesse moyenne sur le mi-pas
  222. VEPSIN12=UNDEMI*(VEPSIN1+VEPSIN2)
  223. VSIG12=UNDEMI*(VSIG1+VSIG2)
  224. SIG2=SIG+(TAU2*(VSIG-VSIG12))
  225. EPSIN2=EPSIN+TAU2*VEPSIN12
  226. IF (IMOD.EQ.19) XX2=XXIN+TAU2*(UNDEMI*(VXX1+VXX2))
  227. C [3.1] Vitesses de contrainte et de deformation inelastique a partir de l'etat a mi-pas [2.2]
  228. SEQ2=ABS(SIG2)
  229. IF (SEQ2.LT.XPETIT) THEN
  230. VEPSIN3=XZER
  231. IF (IMOD.EQ.19) VXX3=XZER
  232. ELSE
  233. C Modele de Norton
  234. IF(IMOD.EQ.15) THEN
  235. CALL NORTON(SEQ2,EPSIN2,TAU2,PROPS(6),PROPS(7),PROPS(8),
  236. & VEPSIN3)
  237. C Modele Polynomial
  238. ELSEIF(IMOD.EQ.16) THEN
  239. VEPSIN3=PROPS(6)+PROPS(7)*(SEQ2**PROPS(8))+
  240. & PROPS(9)*(SEQ2**PROPS(10)) +
  241. & PROPS(11)*(SEQ2**PROPS(12))
  242. C Modele de Blackburn
  243. ELSEIF(IMOD.EQ.17) THEN
  244. CALL SOLUTN(20,SEQ2,EPSIN2,PROPS,DX,TIMEX)
  245. IF(TIMEX.LT.TAU2) TIMEX=ABS(TAU2)
  246. CALL EQUATN(20,SEQ2,TIMEX,XZER,PROPS,VEPSINA,VEPSINB)
  247. VEPSIN3=VEPSINA+VEPSINB
  248. C Modele de Blackburn_2
  249. ELSEIF(IMOD.EQ.18) THEN
  250. CALL SOLUTN(61,SEQ2,EPSIN2,PROPS,DX,TIMEX)
  251. IF(TIMEX.LT.TAU2) TIMEX=ABS(TAU2)
  252. CALL EQUATN(61,SEQ2,TIMEX,XZER,PROPS,VEPSINA,VEPSINB)
  253. VEPSIN3=VEPSINA+VEPSINB
  254. C Modele de Lemaitre
  255. ELSEIF(IMOD.EQ.19) THEN
  256. CALL LEMAIT(SEQ2,XX2,VXX3,TAU2,PROPS,VEPSIN3)
  257. ENDIF
  258. ENDIF
  259. VSIG3=DSIGN(YOU*VEPSIN3,SIG2)
  260. C [3.2] 1ere estimation a la fin du pas, avec la vitesse a mi-pas
  261. SIG3=SIG2+(TAU2*(VSIG-VSIG3))
  262. EPSIN3=EPSIN2+TAU2*VEPSIN3
  263. IF (IMOD.EQ.19) XX3=XXIN+TAU2*VXX3
  264. C [4.1] Vitesses de contrainte et de deformation inelastique a partir de l'etat a la fin du pas [3.2]
  265. SEQ3=ABS(SIG3)
  266. IF (SEQ3.LT.XPETIT) THEN
  267. VEPSIN4=XZER
  268. IF (IMOD.EQ.19) VXX4=XZER
  269. ELSE
  270. C Modele de Norton
  271. IF(IMOD.EQ.15) THEN
  272. CALL NORTON(SEQ3,EPSIN3,TAU2,PROPS(6),PROPS(7),PROPS(8),
  273. & VEPSIN4)
  274. C Modele Polynomial
  275. ELSEIF(IMOD.EQ.16) THEN
  276. VEPSIN4=PROPS(6)+PROPS(7)*(SEQ3**PROPS(8))+
  277. & PROPS(9)*(SEQ3**PROPS(10)) +
  278. & PROPS(11)*(SEQ3**PROPS(12))
  279. C Modele de Blackburn
  280. ELSEIF(IMOD.EQ.17) THEN
  281. CALL SOLUTN(20,SEQ3,EPSIN3,PROPS,DX,TIMEX)
  282. IF(TIMEX.LT.TAU2) TIMEX=ABS(TAU2)
  283. CALL EQUATN(20,SEQ3,TIMEX,XZER,PROPS,VEPSINA,VEPSINB)
  284. VEPSIN4=VEPSINA+VEPSINB
  285. C Modele de Blackburn_2
  286. ELSEIF(IMOD.EQ.18) THEN
  287. CALL SOLUTN(61,SEQ3,EPSIN3,PROPS,DX,TIMEX)
  288. IF(TIMEX.LT.TAU2) TIMEX=ABS(TAU2)
  289. CALL EQUATN(61,SEQ3,TIMEX,XZER,PROPS,VEPSINA,VEPSINB)
  290. VEPSIN4=VEPSINA+VEPSINB
  291. C Modele de Lemaitre
  292. ELSEIF(IMOD.EQ.19) THEN
  293. CALL LEMAIT(SEQ3,XX3,VXX4,TAU2,PROPS,VEPSIN4)
  294. ENDIF
  295. ENDIF
  296. VSIG4=DSIGN(YOU*VEPSIN4,SIG3)
  297. C [4.2] 2nde estimation a la fin du pas, avec la vitesse moyenne sur le mi-pas
  298. VEPSIN34=UNDEMI*(VEPSIN3+VEPSIN4)
  299. VSIG34=UNDEMI*(VSIG3+VSIG4)
  300. SIG4=SIG2+(TAU2*(VSIG-VSIG34))
  301. EPSIN4=EPSIN2+TAU2*VEPSIN34
  302. IF (IMOD.EQ.19) XX4=XXIN+TAU2*(UNDEMI*(VXX3+VXX4))
  303. C [5.1] Vitesse de contrainte et de deformation inelastique a partir de l'etat a la fin du pas [4.2]
  304. SEQ4=ABS(SIG4)
  305. IF (SEQ4.LT.XPETIT) THEN
  306. VEPSIN5=XZER
  307. IF (IMOD.EQ.19) VXX5=XZER
  308. ELSE
  309. C Modele de Norton
  310. IF(IMOD.EQ.15) THEN
  311. CALL NORTON(SEQ4,EPSIN4,TAU2,PROPS(6),PROPS(7),PROPS(8),
  312. & VEPSIN5)
  313. C Modele Polynomial
  314. ELSEIF(IMOD.EQ.16) THEN
  315. VEPSIN5=PROPS(6)+PROPS(7)*(SEQ4**PROPS(8))+
  316. & PROPS(9)*(SEQ4**PROPS(10)) +
  317. & PROPS(11)*(SEQ4**PROPS(12))
  318. C Modele de Blackburn
  319. ELSEIF(IMOD.EQ.17) THEN
  320. CALL SOLUTN(20,SEQ4,EPSIN4,PROPS,DX,TIMEX)
  321. IF(TIMEX.LT.TAU2) TIMEX=ABS(TAU2)
  322. CALL EQUATN(20,SEQ4,TIMEX,XZER,PROPS,VEPSINA,VEPSINB)
  323. VEPSIN5=VEPSINA+VEPSINB
  324. C Modele de Blackburn_2
  325. ELSEIF(IMOD.EQ.18) THEN
  326. CALL SOLUTN(61,SEQ4,EPSIN4,PROPS,DX,TIMEX)
  327. IF(TIMEX.LT.TAU2) TIMEX=ABS(TAU2)
  328. CALL EQUATN(61,SEQ4,TIMEX,XZER,PROPS,VEPSINA,VEPSINB)
  329. VEPSIN5=VEPSINA+VEPSINB
  330. C Modele de Lemaitre
  331. ELSEIF(IMOD.EQ.19) THEN
  332. CALL LEMAIT(SEQ4,XX4,VXX5,TAU2,PROPS,VEPSIN5)
  333. ENDIF
  334. ENDIF
  335. VSIG5=DSIGN(YOU*VEPSIN5,SIG4)
  336. C [5.2] 3eme estimation a la fin du pas avec les vitesses debut [1.1], mi [3.1] et fin de pas [5.1]
  337. VEPSI135=((VEPSIN1+VEPSIN5)/SIX)+(VEPSIN3*DEUX/TROIS)
  338. VSIG135=((VSIG1+VSIG5)/SIX)+(VSIG3*DEUX/TROIS)
  339. SIG5=SIG+(TAU*(VSIG-VSIG135))
  340. EPSIN5=EPSIN+TAU*VEPSI135
  341. IF (IMOD.EQ.19) VXX135=((VXX1+VXX5)/SIX)+(VXX3*DEUX/TROIS)
  342. IF (IMOD.EQ.19) XX5=XXIN+TAU*VXX135
  343. C Calcul du rapport : erreur calculee / erreur admise
  344. XX=SIG5-SIG4
  345. RA=ABS(XX)/ERRABS
  346. SQRA=SQRT(RA)
  347. C Test de fin d'iterations / mise a jour de tau
  348. C DIV =7 BORNE = 2
  349. C Si SQRA>7 alors TAU = TAU/7 puis nouvel essai
  350. C Si 2<RA<7*7 on vise RA = 1 puis nouvel essai
  351. IF (DTLIBR) THEN
  352. C Petite ruse pour dejouer l'optimisation
  353. RA1=RA*UN
  354. IF ((RA.GT.DIV*DIV).OR.(RA.NE.RA1)) THEN
  355. TAU=TAU/DIV
  356. IF(TAU.LT.1.D-30) THEN
  357. CALL ERREUR(5)
  358. ENDIF
  359. TTAUF=TTAU0+TAU
  360. GOTO 30
  361. ELSEIF (RA.GT.BORNE) THEN
  362. TAU = TAU/SQRA
  363. IF(TAU.LT.1.D-30) THEN
  364. CALL ERREUR(5)
  365. ENDIF
  366. TTAUF=TTAU0+TAU
  367. GOTO 30
  368. ENDIF
  369. ENDIF
  370. C
  371. C Fin des iterations sur tau optimal / mise a jour des variables
  372. C Ici RA < BORNE
  373. C On avance en temps
  374. SIG=SIG4
  375. EPSIN=EPSIN+TAU2*(VEPSIN12+VEPSIN34)
  376. IF (IMOD.EQ.19) XXIN=XX5
  377. C Test de fin des iterations en sous increments / mise a jour de tau
  378. C Si SQRA<1/3 TAU = TAU*3
  379. C Si 1/3<SQRA<RMIN on vise RA = 1
  380. C Si RMIN<SQRA<RMAX TAU inchange
  381. C Si SQRA>RMAX on vise RA = 1
  382. C Fin des boucles en sous increments si tau = DTLEFT
  383. IF (TAU.LT.DTLEFT) THEN
  384. TTAU0=TTAU0+TAU
  385. DTLEFT=DTLEFT-TAU
  386. IF (FAC*SQRA.LT.UN) THEN
  387. TAU=TAU*FAC
  388. ELSEIF ((SQRA.LT.RMIN).OR.(SQRA.GT.RMAX)) THEN
  389. TAU=TAU/SQRA
  390. ENDIF
  391. IF (TAU.GT.DTLEFT) TAU=DTLEFT
  392. TTAUF=TTAU0+TAU
  393. GOTO 20
  394. ENDIF
  395. C Ici, on a atteint une subdivision du pas de temps trop petite,
  396. C --> on sort en erreur
  397. IF (ABS(TAU-DTLEFT).GT.(TAU/1000.)) THEN
  398. PRINT*,'ERREUR NBR. MAX. DE SUBDIVISION DU PAS DE TEMPS ATTEINT'
  399. CALL ERREUR(996)
  400. RETURN
  401. ENDIF
  402. C
  403. C=======================================================================
  404. C 5 - Stockage des resultats et sortie
  405. C=======================================================================
  406. SIGF(1)=SIG
  407. SIGF(2)=SIG0(2)+(G*DEPS(2))
  408. SIGF(3)=SIG0(3)+(G*DEPS(3))
  409. IF (IMOD.EQ.19) THEN
  410. VARF(1)=XXIN
  411. VARF(2)=EPSIN
  412. ELSE
  413. VARF(1)=EPSIN
  414. ENDIF
  415. C Sortie du programme
  416. RETURN
  417. END
  418.  
  419.  
  420.  
  421.  

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