Télécharger cflue2.eso

Retour à la liste

Numérotation des lignes :

  1. C CFLUE2 SOURCE BP208322 17/03/01 21:15:31 9325
  2. SUBROUTINE CFLUE2(wrk52,wrk53,wrk54,wrk2,wrk3,
  3. & IB,IGAU,NBPGAU,NBGMAT,NELMAT,IFOURB)
  4. *
  5. * modele `lemaitre fluage` Lemaitre et Chaboche pp.289,411
  6. * voir aussi syntheses RUPTHER SEMT/LISN/RT/99-124/A
  7. * remarque ep et sigma' sont colineaires donc ep et (Sigma0 + K depst)' le sont
  8. * Runge_Kutta d ordre 2 : Sigma1 = 1/2 YK0 + 1/2 YK1
  9. * YK0 = seq0 * (Sigma0 + K depst) / seqtot
  10. * etablit increment de r et D
  11. * YK1 = seq1 * (YK0 + K depst) / seqtot
  12. * etablit increment de r et D
  13. * puis moyenne
  14. IMPLICIT INTEGER(I-N)
  15. IMPLICIT REAL*8 (A-H,O-Z)
  16. -INC CCOPTIO
  17. -INC DECHE
  18. *
  19. *
  20. *
  21. SEGMENT WRK2
  22. REAL*8 TRAC(LTRAC)
  23. ENDSEGMENT
  24. *
  25. SEGMENT WRK3
  26. REAL*8 WORK(LW),WORK2(LW2)
  27. ENDSEGMENT
  28. *
  29. dimension YK0(8),YK1(8),spri0(8),delep0(8),spri1(8),delep1(8),
  30. &DIV(8)
  31.  
  32.  
  33. d0 = var0(3)
  34. r0 = var0(2)
  35.  
  36. * cas isotrope
  37. ip1 = 4
  38. *
  39. youn0 = xmat0(1)
  40. sigy0 = xmat0(ip1+1)
  41. xn0 = xmat0(ip1+2)
  42. xm0 = xmat0(ip1+3)
  43. gk0 = xmat0(ip1+4)
  44. alp0 = xmat0(ip1+5)
  45. blp0 = xmat0(ip1+6)
  46. rini0 = xmat0(ip1+7)
  47. a0 = xmat0(ip1+8)
  48. pk0 = xmat0(ip1+9)
  49. x2mu0= xmat0(1)/(1.+xmat0(2))
  50.  
  51. if (d0.le.0.) d0 = 0.D0
  52. if (1. - d0.le.0.) then
  53. do ic = 1,nstrs
  54. sigf(ic) = 0.
  55. enddo
  56.  
  57. depse = 0.d0
  58. do is = 1, nstrs
  59. defp(is) = depst(is)
  60. depse = depse + defp(is)*defp(is)
  61. enddo
  62.  
  63. varf(1) = var0(1) + SQRT(depse*0.5)
  64. varf(2) = var0(2)
  65. varf(3) = 1.
  66. goto 1002
  67. endif
  68.  
  69. if (r0.lt.0.) then
  70. c write(6,*) 'erreur valeur r ', r0
  71. moterr(1:16) = conm
  72. moterr(17:24) = 'CFLUE2-3'
  73. call erreur(943)
  74. return
  75. endif
  76. if (r0.eq.0.) r0 = 1.e-10
  77.  
  78. if(ib.eq.1.and.igau.eq.1) then
  79. * write(6,*) 't0' , youn0, sigy0, xn0 ,xm0 ,gk0,pk0
  80. endif
  81.  
  82. youn1 = xmat(1)
  83. sigy1 = xmat(ip1+1)
  84. xn1 = xmat(ip1+2)
  85. xm1 = xmat(ip1+3)
  86. gk1 = xmat(ip1+4)
  87. alp1 = xmat(ip1+5)
  88. blp1 = xmat(ip1+6)
  89. rini1 = xmat(ip1+7)
  90. a1 = xmat(ip1+8)
  91. pk1 = xmat(ip1+9)
  92. x2mu1= xmat(1)/(1.+xmat(2))
  93.  
  94. if(ib.eq.1.and.igau.eq.1) then
  95. * write(6,*) 't1' , youn1, sigy1, xn1 ,xm1 ,gk1,pk1
  96. endif
  97. *
  98. delt = tempf - temp0
  99. * gerer la variable cachee avec les variations de temperature
  100.  
  101. C---------CARACTERISTIQUES GEOMETRIQUES---------------------------------
  102. C
  103. C COQUES
  104. C
  105. ALFAH=1.D0
  106. IF(MFR.EQ.3.OR.MFR.EQ.5.OR.MFR.EQ.9) THEN
  107. EP1=work(1)
  108. IF(MFR.NE.5) ALFAH=work(2)**2
  109. ENDIF
  110. C---------COQUES AVEC CT------------------------------------------------
  111. C ON NE TRAVAILLE QUE SUR LES 6 PREMIERES COMPOSANTES
  112.  
  113. IF(MFR.EQ.9) THEN
  114. NCONT=6
  115. ELSE
  116. NCONT=NSTRS
  117. ENDIF
  118.  
  119. * calcul YK0
  120. * determine direction et sens de eplas
  121. * calcul increments de contrainte / utilise xmatf
  122. * remarque on utilise les caracteristiques elastiques a la date finale
  123.  
  124. do ic = 1,valmat(/1)
  125. xmatf(ic) = valmat(ic)
  126. * completer pour le cas anisotrope
  127. if (ic.eq.1) xmatf(ic) = xmatf(ic) * (1. - d0)
  128. enddo
  129.  
  130. CALL CALSIG(depst,DDAUX,NSTAUX,CMATE,xmatf,VALCAR,
  131. & N2EL,N2PTEL,MFR,IFOURB,IB,IGAU,EPAIST,
  132. & NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,TXR,
  133. & XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI,dsigt,IRTD)
  134.  
  135. IF(IRTD.NE.1) THEN
  136. KERRE=69
  137. GOTO 1010
  138. ENDIF
  139.  
  140. DO I=1,NSTRS
  141. DSIGT(I)=SIG0(I) + dsigt(i)
  142. ENDDO
  143.  
  144. C---------CAS DES POUTRES-----------------------------------------------
  145.  
  146. IF(MFR.EQ.7) THEN
  147. DIV(1)=1.D0/work(4)
  148. DIV(2)=1.D0
  149. DIV(3)=1.D0
  150. DIV(4)=work(10)/work(1)
  151. DIV(5)=work(11)/work(2)
  152. DIV(6)=work(12)/work(3)
  153. IF(DIV(4).EQ.0.D0) DIV(4)=1.D-10/SQRT(work(1)*work(4))
  154. IF(DIV(5).EQ.0.D0) DIV(5)=1.D-10/SQRT(work(2)*work(4))
  155. IF(DIV(6).EQ.0.D0) DIV(6)=1.D-10/SQRT(work(3)*work(4))
  156. DO I=1,NCONT
  157. DSIGT(I)= DSIGT(I)*DIV(I)
  158. ENDDO
  159. ENDIF
  160.  
  161. C-----------------------------------------------------------------------
  162. C CALCUL DE LA CONTRAINTE EQUIVALENTE SEQ
  163. C-----------------------------------------------------------------------
  164. seqtot=VONMIS0(dsigt,NSTRS,MFR,IFOURB,EP1,ALFAH)
  165.  
  166. if (seqtot.gt.0.) then
  167. seq0 = seqtot
  168. else
  169. * contrainte spherique
  170. do ic = 1,nstrs
  171. YK0(ic) = dsigt(ic)
  172. enddo
  173.  
  174. do is = 1, nstrs
  175. delep0(is) = 0.
  176. enddo
  177.  
  178. delr0 = 0.
  179. deld0 = 0.
  180. goto 500
  181. endif
  182.  
  183. * point fixe pour determiner le multiplicateur de sigtot'/seqtot
  184. icaz = 1
  185. do ipfx = 1,10
  186. if(ib.eq.1.and.igau.eq.1) then
  187. c write(6,*) 'pt fixe' , ipfx, seq0,seqtot,icaz
  188. endif
  189. goto (70,80) icaz
  190. 70 continue
  191. * fonction
  192. if (gk1.gt.0..and.xm1.gt.0..and.r0.gt.0..and.(1.-d0).gt.0.) then
  193. delr0 = seq0/(1.-d0)/gk1/(r0 ** (1/xm1))
  194. else
  195. delr0 = 0.
  196. endif
  197. delr0 = max(delr0,0.d0)
  198. c cas delr0 nul ?
  199. if (xn1.gt.0..and.x2mu1.gt.0.) then
  200. delr0 = delr0 ** xn1
  201. else
  202. moterr(1:16) = conm
  203. moterr(17:24) = 'CFLUE2-1'
  204. call erreur(943)
  205. return
  206. endif
  207. xmult = 1.5 * delr0 * x2mu1 * delt / (1. - d0) / (1. - d0)
  208. if(ib.eq.1.and.igau.eq.1) then
  209. c write(6,*) 'delr0' , delr0
  210. endif
  211. seq01 = seqtot - xmult
  212. goto 90
  213.  
  214. * fonction associee
  215. 80 continue
  216. xmult = seqtot - seq0
  217. c en principe 1. - d0 > 0
  218. if (1. - d0.lt.0..or.delt.lt.0..or.x2mu1.le.0.) then
  219. * write(6,*) 'erreur parametres 1' , delt, x2mu1, d0
  220. moterr(1:16) = conm
  221. moterr(17:24) = 'CFLUE2-2'
  222. call erreur(943)
  223. return
  224. endif
  225. delr0 = xmult / delt * (1.- d0) * (1. - d0) /x2mu1 / 1.5
  226. if (xn1.gt.0) then
  227. delr0 = delr0 ** (1/xn1)
  228. * write(6,*) ' delr0 ', delr0
  229. else
  230. c write(6,*) 'erreur parametres 2' , xn1
  231. moterr(1:16) = conm
  232. moterr(17:24) = 'CFLUE2-4'
  233. call erreur(943)
  234. return
  235. endif
  236. if (gk1.gt.0..and.xm1.gt.0.) then
  237. seq01 = delr0 * (r0 ** (1/xm1)) * gk1 * (1. - d0)
  238. else
  239. c write(6,*) 'erreur parametres 3' , gk1,xm1
  240. moterr(1:16) = conm
  241. moterr(17:24) = 'CFLUE2-5'
  242. call erreur(943)
  243. return
  244. endif
  245.  
  246. 90 continue
  247. if (ipfx.eq.1) then
  248. if (seq01.lt.0) then
  249. icaz = 2
  250. seq01 = seqtot/2.
  251. else if (seq01.eq.0.) then
  252. icaz =1
  253. seq01 = seqtot/2.
  254. endif
  255. endif
  256. varseq = abs((seq01 - seq0) / seq0)
  257. if (ib.eq.1.and.igau.eq.1) then
  258. c write(6,*) 'variation relative', varseq
  259. endif
  260. seq0 = seq01
  261. if (seq0 .lt.0.) then
  262. c write(6,*) 'erreur point fixe', seqtot, seq0
  263. moterr(1:16) = conm
  264. moterr(17:24) = 'CFLUE2-6'
  265. call erreur(943)
  266. return
  267. endif
  268.  
  269. enddo
  270.  
  271. 100 if (seqtot - seq0 .lt.0.) then
  272. c write(6,*) 'erreur point fixe', seqtot, seq0
  273. moterr(1:16) = conm
  274. moterr(17:24) = 'CFLUE2-7'
  275. call erreur(943)
  276. return
  277. endif
  278.  
  279. *a revoir pour alp0 et beta0 .
  280. xki0 = seq0
  281. if ((xki0*a1).gt.0.and.(1. - d0).gt.0..
  282. &and.pk1.gt.0.and.a1.ne.0..and.rini1.gt.0.)then
  283. deld0 = ((xki0 /a1) ** rini1) * ((1. - d0) ** (pk1 * (-1.)))
  284. else
  285. deld0 = 0.
  286. endif
  287.  
  288. if(ib.eq.1.and.igau.eq.1) then
  289. * write(6,*) 'deld0' , deld0
  290. endif
  291.  
  292. trsig0 = (dsigt(1) + dsigt(2) + dsigt(3)) * 0.3333333333
  293. spri0(1) = dsigt(1) - trsig0
  294. spri0(2) = dsigt(2) - trsig0
  295. spri0(3) = dsigt(3) - trsig0
  296. do is = 4,nstrs
  297. spri0(is) = dsigt(is)
  298. enddo
  299.  
  300. if (deld0*delt.lt.1.) then
  301.  
  302. if (seq0.gt.0.and.(1. - d0).gt.0..and.seqtot - seq0.gt.0.) then
  303. coep0 = (seqtot - seq0) * (1. - d0) / x2mu1/ seqtot
  304. else
  305. coep0 = 0.d0
  306. endif
  307. do is = 1, nstrs
  308. delep0(is) = spri0(is) * coep0
  309. enddo
  310.  
  311. do ic = 1,nstrs
  312. YK0(ic) = spri0(ic) * seq0 / seqtot
  313. if (ic.le.3) YK0(ic) = YK0(ic) + trsig0 * seq0/seqtot
  314. enddo
  315.  
  316. else
  317.  
  318. do is = 1, nstrs
  319. delep0(is) = depst(is)
  320. enddo
  321.  
  322. do ic = 1,nstrs
  323. YK0(ic) = 0.
  324. enddo
  325. endif
  326.  
  327. 500 continue
  328. * calcul YK1
  329. * gerer la variable cachee avec les variations de temperature
  330. r1 = r0 + delr0 * delt
  331. if (r1.lt.0.) then
  332. c write(6,*) 'erreur valeur r ', r1
  333. moterr(1:16) = conm
  334. moterr(17:24) = 'CFLUE2-8'
  335. call erreur(943)
  336. return
  337. endif
  338. if (r1.eq.0.) r1 = 1.e-10
  339.  
  340. d1 = var0(3) + deld0*delt
  341. if (d1.le.0.) d1 = 0.D0
  342. if (d1.ge.1.) then
  343. do ic = 1,nstrs
  344. YK1(ic) = 0.
  345. enddo
  346. deld1 = 0.
  347. delr1 = 0.
  348. do is = 1, nstrs
  349. delep1(is)= depst(is)
  350. enddo
  351. goto 1000
  352. endif
  353.  
  354. do ic = 1,valmat(/1)
  355. xmatf(ic) = valmat(ic)
  356. * completer pour le cas anisotrope
  357. if (ic.eq.1) xmatf(ic) = xmatf(ic) * (1. - d1)
  358. enddo
  359.  
  360. CALL CALSIG(depst,DDAUX,NSTAUX,CMATE,xmatf,VALCAR,
  361. & N2EL,N2PTEL,MFR,IFOURB,IB,IGAU,EPAIST,
  362. & NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,TXR,
  363. & XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI,dsigt,IRTD)
  364.  
  365. IF(IRTD.NE.1) THEN
  366. KERRE=69
  367. GOTO 1010
  368. ENDIF
  369.  
  370. DO I=1,NSTRS
  371. DSIGT(I)=SIG0(I) + dsigt(i)
  372. ENDDO
  373. C---------CAS DES POUTRES-----------------------------------------------
  374.  
  375. IF(MFR.EQ.7) THEN
  376. DIV(1)=1.D0/work(4)
  377. DIV(2)=1.D0
  378. DIV(3)=1.D0
  379. DIV(4)=work(10)/work(1)
  380. DIV(5)=work(11)/work(2)
  381. DIV(6)=work(12)/work(3)
  382. IF(DIV(4).EQ.0.D0) DIV(4)=1.D-10/SQRT(work(1)*work(4))
  383. IF(DIV(5).EQ.0.D0) DIV(5)=1.D-10/SQRT(work(2)*work(4))
  384. IF(DIV(6).EQ.0.D0) DIV(6)=1.D-10/SQRT(work(3)*work(4))
  385. DO I=1,NCONT
  386. DSIGT(I)= DSIGT(I)*DIV(I)
  387. ENDDO
  388. ENDIF
  389.  
  390. C-----------------------------------------------------------------------
  391. C CALCUL DE LA CONTRAINTE EQUIVALENTE SEQ
  392. C-----------------------------------------------------------------------
  393. seqtot=VONMIS0(dsigt,NSTRS,MFR,IFOURB,EP1,ALFAH)
  394.  
  395.  
  396. if (seqtot.gt.0.) then
  397. seq1 = seqtot
  398. else
  399. * contrainte spherique
  400. do ic = 1,nstrs
  401. YK1(ic) = dsigt(ic)
  402. enddo
  403.  
  404. do is = 1, nstrs
  405. delep1(is) = 0.
  406. enddo
  407.  
  408. delr1 = 0.
  409. deld1 = 0.
  410. goto 1000
  411. endif
  412.  
  413. * point fixe pour determiner le multiplicateur de sigtot'/seqtot
  414. icaz = 1
  415. do ipfx = 1,10
  416. if(ib.eq.1.and.igau.eq.1) then
  417. c write(6,*) 'pt fixe' , ipfx, seq1,seqtot,icaz
  418. endif
  419. goto (570,580) icaz
  420.  
  421. 570 continue
  422. * fonction
  423. if (gk1.gt.0..and.xm1.gt.0..and.r1.gt.0..and.(1.-d1).gt.0.) then
  424. delr1 = seq1/(1.-d1)/gk1/(r1 ** (1/xm1))
  425. else
  426. delr1 = 0.
  427. endif
  428. delr1 = max(delr1,0.d0)
  429. c cas delr1 nul ?
  430. if (xn1.gt.0.) then
  431. delr1 = delr1 ** xn1
  432. else
  433. c write(6,*) 'erreur parametres 2' , xn1
  434. moterr(1:16) = conm
  435. moterr(17:24) = 'CFLUE2-9'
  436. call erreur(943)
  437. return
  438. endif
  439. if(ib.eq.1.and.igau.eq.1) then
  440. * write(6,*) 'delr1' , delr1
  441. endif
  442. xmult = 1.5 * delr1 * x2mu1 / (1.- d1) / (1. - d1) * delt
  443. seq11 = seqtot - xmult
  444. goto 590
  445.  
  446. * fonction associee
  447. 580 continue
  448. xmult = seqtot - seq1
  449. c en principe 1. - d1 > 0
  450. if (1. - d1.lt.0..or.delt.lt.0..or.x2mu1.le.0.) then
  451. c write(6,*) 'erreur parametres 1' , delt, x2mu1, d1
  452. moterr(1:16) = conm
  453. moterr(17:24) = 'CFLUE210'
  454. call erreur(943)
  455. return
  456. endif
  457. delr1 = xmult / delt * (1.- d1) * (1. - d1) /x2mu1 / 1.5
  458. if (xn1.gt.0) then
  459. delr1 = delr1 ** (1/xn1)
  460. * write(6,*) ' delr1 ', delr1
  461. else
  462. c write(6,*) 'erreur parametres 2' , xn1
  463. moterr(1:16) = conm
  464. moterr(17:24) = 'CFLUE211'
  465. call erreur(943)
  466. return
  467. endif
  468. if (gk1.gt.0..and.xm1.gt.0.) then
  469. seq11 = delr1 * (r1 ** (1/xm1)) * gk1 * (1. - d1)
  470. else
  471. c write(6,*) 'erreur parametres 3' , gk1,xm1
  472. moterr(1:16) = conm
  473. moterr(17:24) = 'CFLUE212'
  474. call erreur(943)
  475. return
  476. endif
  477.  
  478. 590 continue
  479. if (ipfx.eq.1) then
  480. if (seq11.lt.0) then
  481. icaz = 2
  482. seq11 = seqtot/2.
  483. else if (seq11.eq.0.) then
  484. icaz =1
  485. seq11 = seqtot/2.
  486. endif
  487. endif
  488. varseq = abs((seq11 - seq1) / seq1)
  489. if(ib.eq.1.and.igau.eq.1) then
  490. c write(6,*) 'variation relative', varseq
  491. endif
  492. seq1 = seq11
  493. if (seq1 .lt.0.) then
  494. c write(6,*) 'erreur point fixe', seqtot, seq1
  495. moterr(1:16) = conm
  496. moterr(17:24) = 'CFLUE213'
  497. call erreur(943)
  498. return
  499. endif
  500.  
  501.  
  502. enddo
  503.  
  504.  
  505. if (seqtot - seq1 .lt.0.) then
  506. c write(6,*) 'erreur point fixe', seqtot, seq1
  507. moterr(1:16) = conm
  508. moterr(17:24) = 'CFLUE214'
  509. call erreur(943)
  510. return
  511. endif
  512.  
  513.  
  514. *a revoir pour alp1 et bta1 .
  515. xki1 = seq1
  516. if ((xki1*a1).gt.0.and.(1. - d1).gt.0..
  517. &and.pk1.gt.0.) then
  518. deld1 = ((xki1 /a1) ** rini1) * ((1. - d1) ** (pk1 * (-1.)))
  519. else
  520. deld1 = 0.
  521. endif
  522. if(ib.eq.1.and.igau.eq.1) then
  523. * write(6,*) 'deld1' , deld1
  524. endif
  525.  
  526. trsig1 = (dsigt(1) + dsigt(2) + dsigt(3)) * 0.3333333333
  527. spri1(1) = dsigt(1) - trsig1
  528. spri1(2) = dsigt(2) - trsig1
  529. spri1(3) = dsigt(3) - trsig1
  530. do is = 4,nstrs
  531. spri1(is) = dsigt(is)
  532. enddo
  533.  
  534. if (deld1*delt.lt.1) then
  535. if (seq1.gt.0.and.(1. - d1).gt.0..and.seqtot - seq1.gt.0.) then
  536. coep1 = (seqtot - seq1) * (1. - d1) / x2mu1/ seqtot
  537. else
  538. coep1 = 0.
  539. endif
  540. do is = 1, nstrs
  541. delep1(is) = spri1(is) * coep1
  542. enddo
  543.  
  544. do ic = 1,nstrs
  545. YK1(ic) = spri1(ic) * seq1 / seqtot
  546. if (ic.le.3) YK1(ic) = YK1(ic) + trsig1 * seq1 / seqtot
  547. enddo
  548.  
  549. else
  550.  
  551. do is = 1, nstrs
  552. delep1(is) = depst(is)
  553. enddo
  554.  
  555. do ic = 1,nstrs
  556. YK1(ic) = 0.
  557. enddo
  558.  
  559. endif
  560.  
  561. c au final
  562. 1000 continue
  563. do ic = 1,nstrs
  564. sigf(ic) = 0.5*YK0(ic) + 0.5*YK1(ic)
  565. enddo
  566.  
  567. depse = 0.d0
  568. do is = 1, nstrs
  569. defp(is) = 0.5*(delep0(is) + delep1(is))
  570. depse = depse + defp(is)*defp(is)
  571. enddo
  572.  
  573. varf(1) = var0(1) + SQRT(depse*0.5)
  574. vtd = var0(3) + (deld0 + deld1)*0.5*delt
  575. if (vtd.le.1.) then
  576. varf(3) = vtd
  577. else
  578. varf(3) = 1.
  579. endif
  580. varf(2) = (r0 + r1)*0.5
  581.  
  582.  
  583. 1002 continue
  584. if(ib.eq.1.and.igau.eq.1) then
  585. * write(6,*) 't0' ,sigf(3),yk0(3),yk1(3),varf(3),depst(3)
  586. endif
  587.  
  588. 1010 continue
  589. RETURN
  590. END
  591.  
  592.  
  593.  
  594.  
  595.  
  596.  
  597.  
  598.  
  599.  
  600.  

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