Télécharger cflue2.eso

Retour à la liste

Numérotation des lignes :

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

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