Télécharger dynar.eso

Retour à la liste

Numérotation des lignes :

  1. C DYNAR SOURCE PV 17/12/08 21:17:22 9660
  2. SUBROUTINE DYNAR(WRK52,WRK53,WRK54,IECOU)
  3. c
  4. C=======================================================================
  5. C
  6. C MODELE VISCOPLASTIQUE VISCOENDOMMAGEABLE POUR LA DYNAMIQUE
  7. C - PROGRAMME PRINCIPAL -
  8. C VERSION 2.0
  9. C
  10. C=======================================================================
  11. C CREATION :F.GATUINGT
  12. C=======================================================================
  13. C=========== Hypopthese: au temps t=0 tout est nul ===================
  14. C=======================================================================
  15. IMPLICIT REAL*8(A-H,O-Z)
  16. -INC CCOPTIO
  17. -INC DECHE
  18. DIMENSION EPSPTDT(6),EPSTE(6),EPSPOPT(6),EPSPO(6),SIGPO(6),
  19. & sigt(6),sigtdt(6),EPST(6),EPSTDT(6),EPSPT(6),EPSTDTE(6),
  20. & S(3,3),EPSIPP(3),ARRAY(3,3),SIGTDTE(6),DFIDSIG(6),EPSPOPTDT(6)
  21. C
  22. PARAMETER (XUnDemi=0.5D0,XZERO=0.D0,UN=1.D0,DEUX=2.D0)
  23. PARAMETER (XPETIT=1.D-12)
  24. PARAMETER (TROIS=3.D0)
  25. C
  26. SEGMENT IECOU
  27. * COMMON/IECOU/NYOG,NYNU,NYALFA,NYSMAX,NYN,NYM,NYKK,
  28. INTEGER icow1,icow2,icow3,icow4,icow5,icow6,icow7,
  29. C INTEGER NYOG, NYNU, NYALFA,NYSMAX,NYN, NYM, NYKK,
  30. 1 icow8,icow9,icow10,icow11,icow12,icow13,icow14,icow15,icow16,
  31. C . NYALF1,NYBET1,NYR, NYA, NYRHO,NSIGY, NNKX, NYKX, IND,
  32. 2 icow17,icow18,icow19,icow20,icow21,icow22,icow23,icow24,
  33. C . NSOM, NINV, NINCMA,NCOMP, JELEM, LEGAUS,INAT, NCXMAT,
  34. 3 icow25,icow26,icow27,icow28,icow29,icow30,ICARA,
  35. C . LTRAC, MFR, IELE, NHRM, NBNN, NBELEM,ICARA,
  36. 4 icow32,icow33,NSTRS1,MFR1,icow36,icow37,icow38,
  37. C . LW2, NDEF, NSTRSS,MFR1, NBGMAT,NELMAT,MSOUPA,
  38. 5 icow39,icow40,icow41,icow42,icow43,icow44
  39. C . NUMAT1,LENDO, NBBB, NNVARI,KERR1, MELEME
  40. INTEGER icow45,icow46,icow47,icow48,icow49,icow50,
  41. . icow51,icow52,icow53,icow54,icow55,icow56
  42. . icow57,icow58
  43. ENDSEGMENT
  44. C
  45. C
  46. C
  47. C
  48. C - Recuperation des donnees materiau
  49. C
  50. Eini=XMAT(1)
  51. XNU=XMAT(2)
  52. D0=XMAT(5)
  53. Q1=XMAT(6)
  54. Q2=XMAT(7)
  55. Q3=XMAT(8)
  56. SIGM0=XMAT(9)
  57. XN=XMAT(10)
  58. XKPLAS = XMAT(12)
  59. XMPLAS = XMAT(11)
  60. XK=XMAT(13)
  61. XMDt=XMAT(14)
  62. XNDt=XMAT(15)
  63. At=XMAT(16)
  64. Bt=XMAT(17)
  65. XMDc=XMAT(18)
  66. XNDc=XMAT(19)
  67. Ac=XMAT(20)
  68. Bc=XMAT(21)
  69. ED0=XMAT(22)
  70. C - Pas de temps
  71. DELTAT=TEMPF-TEMP0
  72. C - Contrainte au temps T
  73. SIGT(1)=SIG0(1)
  74. SIGT(2)=SIG0(2)
  75. SIGT(3)=SIG0(3)
  76. SIGT(4)=SIG0(4)
  77. SIGT(5)=SIG0(5)
  78. SIGT(6)=SIG0(6)
  79. C - Deformation au temps T
  80. EPST(1)=EPST0(1)
  81. EPST(2)=EPST0(2)
  82. EPST(3)=EPST0(3)
  83. EPST(4)=EPST0(4)/2.
  84. EPST(5)=EPST0(5)/2.
  85. EPST(6)=EPST0(6)/2.
  86. C - Increment de deformation
  87. DEPST(1)=DEPST(1)
  88. DEPST(2)=DEPST(2)
  89. DEPST(3)=DEPST(3)
  90. DEPST(4)=DEPST(4)/2.
  91. DEPST(5)=DEPST(5)/2.
  92. DEPST(6)=DEPST(6)/2.
  93. C======================================================================
  94. C - Subcycling si le pas de temps est trop grand=======================
  95. C======================================================================
  96. IF (DELTAT.GT.1.D-3) THEN
  97. coef=100
  98. ELSE
  99. coef=max (1.d0 , deltat/ 1.d-5)
  100. ENDIF
  101. nbiter=int(coef)
  102. DO i=1,6
  103. DEPST(i)=DEPST(i)/nbiter
  104. ENDDO
  105. DELTAT=DELTAT/nbiter
  106. l=0
  107. DO WHILE (l.LT.nbiter)
  108. C - Debut des boucles
  109. C=======================================================================
  110. C - Recuperation des variables historiques
  111. C=======================================================================
  112. DT=VAR0(1)
  113. fT=VAR0(2)+D0
  114. EPSPT(1)=VAR0(3)
  115. EPSPT(2)=VAR0(4)
  116. EPSPT(3)=VAR0(5)
  117. EPSPT(4)=VAR0(6)
  118. EPSPT(5)=VAR0(7)
  119. EPSPT(6)=VAR0(8)
  120. SIGMT=VAR0(9)+SIGM0
  121. dtt=VAR0(10)
  122. dct=VAR0(11)
  123. ES=VAR0(12)+ED0
  124. C======================================================================
  125. C - Deformation au temps t+dt
  126. C======================================================================
  127. EPSTDT(1)=EPST(1)+DEPST(1)
  128. EPSTDT(2)=EPST(2)+DEPST(2)
  129. EPSTDT(3)=EPST(3)+DEPST(3)
  130. EPSTDT(4)=EPST(4)+DEPST(4)
  131. EPSTDT(5)=EPST(5)+DEPST(5)
  132. EPSTDT(6)=EPST(6)+DEPST(6)
  133. TRACEPSTDT=EPSTDT(1)+EPSTDT(2)+EPSTDT(3)
  134. C======================================================================
  135. C - Vitesse de Deformation
  136. C======================================================================
  137. DO i=1,6
  138. EPSPO(i)=DEPST(i)/DELTAT
  139. ENDDO
  140. TRACEPSPO=EPSPO(1)+EPSPO(2)+EPSPO(3)
  141. C=======================================================================
  142. C - Module d'elasticite tangent de la Matrice
  143. C=======================================================================
  144. ETAN=(Eini)/XN*(SIGMt/SIGM0)**(1-XN)
  145. C=======================================================================
  146. C - Coefficients de compressibilite et cisaillement de la matrice sans les pores
  147. C=======================================================================
  148. XKm=Eini/(TROIS*(UN-DEUX*XNU))
  149. XGm=Eini/(DEUX*(UN+XNU))
  150. C=======================================================================
  151. C - Coeficients de compressibilite et cisaillement de la matrice avec les pores (Mori-Tanaka)
  152. C=======================================================================
  153. XKporo=4.D0*XKm*XGm*(UN-fT)/(4.D0*XGm+3.D0*XKm*fT)
  154. XGporo=XGm*(UN-fT)/(UN+fT*
  155. & (6.D0*XKm+12.D0*XGm)/(9.D0*XKm+8.D0*XGm))
  156. C ======================================================================
  157. C - Coefficients de compressibilite et cisaillement de la matrice endommagee
  158. C ======================================================================
  159. XKendo=(UN-Dt)*XKporo
  160. XGendo=(UN-Dt)*XGporo
  161. C ======================================================================
  162. C Calcul des invariants de la contrainte
  163. C ======================================================================
  164. XI1t=SIGT(1)+SIGT(2)+SIGT(3)
  165. XJ2t=0.5D0*((SIGT(1)-SIGT(2))**2
  166. & +(SIGT(2)-SIGT(3))**2+(SIGT(3)-SIGT(1))**2
  167. & +6.D0*(SIGT(4)**2+SIGT(5)**2+SIGT(6)**2))
  168. C======================================================================
  169. C - Deformation elastique au temps t+dt
  170. C======================================================================
  171. C Calcul du multiplicateur plastique au temps t
  172. Fit=XJ2t/(SIGMT**2)+2*Q1*fT*COSH((Q2*XI1t)/(2*SIGMT))
  173. & -(1+(Q3*fT)**2)
  174. IF (Fit.LT.0.) THEN
  175. Fit=0.
  176. ENDIF
  177. xlambdapo=(fT/(1-fT))*(Fit/XKPLAS)**XMPLAS
  178.  
  179. C - Calcul des deformations plastique au temps t
  180. C Derivees partielles de fi
  181. DFIDSIG(1)=0.5*(4*SIGT(1)-2*SIGT(2)-2*SIGT(3))/SIGMT**2
  182. & +Q1*Q2*ft*SINH(Q2*XI1t/(2*SIGMt))/SIGMt
  183. DFIDSIG(2)=0.5*(-2*SIGT(1)+4*SIGT(2)-2*SIGT(3))/SIGMT**2
  184. & +Q1*Q2*ft*SINH(Q2*XI1t/(2*SIGMt))/SIGMt
  185. DFIDSIG(3)=0.5*(-2*SIGT(1)-2*SIGT(2)+4*SIGT(3))/SIGMT**2
  186. & +Q1*Q2*ft*SINH(Q2*XI1t/(2*SIGMt))/SIGMt
  187. DFIDSIG(4)=6/(SIGMt**2)*SIGT(4)
  188. DFIDSIG(5)=6/(SIGMt**2)*SIGT(5)
  189. DFIDSIG(6)=6/(SIGMt**2)*SIGT(6)
  190. C
  191. DO i=1,6
  192. EPSPOPT(i)=xlambdapo*DFIDSIG(i)
  193. ENDDO
  194. TRACEPSPOPT=EPSPOPT(1)+EPSPOPT(2)+EPSPOPT(3)
  195. C - Déformation viscoplastique
  196. DO i=1,6
  197. EPSPTDT(i)=EPSPT(i)+EPSPOPT(i)*DELTAT
  198. ENDDO
  199. TRACEPSPTDT=EPSPTDT(1)+EPSPTDT(2)+EPSPTDT(3)
  200. C - Deformation elastique
  201. DO i=1,6
  202. EPSTDTE(i)=EPSTDT(i)-EPSPTDT(i)
  203. ENDDO
  204. TRACEPSTDTE=EPSTDTE(1)+EPSTDTE(2)+EPSTDTE(3)
  205. C======================================================================
  206. C - Deformation equivalente elastique au sens de Mazars a t
  207. ARRAY(1,1)=EPSTDTE(1)
  208. ARRAY(2,2)=EPSTDTE(2)
  209. ARRAY(3,3)=EPSTDTE(3)
  210. ARRAY(1,2)=EPSTDTE(4)
  211. ARRAY(2,3)=EPSTDTE(5)
  212. ARRAY(1,3)=EPSTDTE(6)
  213. ARRAY(3,1)=ARRAY(1,3)
  214. ARRAY(3,2)=ARRAY(2,3)
  215. ARRAY(2,1)=ARRAY(1,2)
  216. C
  217. CALL JACOB3(ARRAY,3,EPSIPP,S)
  218. EPSE=MAX( XZERO , EPSIPP(1) )**2 +
  219. 1 MAX( XZERO , EPSIPP(2) )**2 +
  220. 2 MAX( XZERO , EPSIPP(3) )**2
  221. EPSE=SQRT(EPSE)
  222. C - Nouveau seuil de Mazars
  223. ES=max (ES, EPSE)
  224. C ==============================================================
  225. C - Prediction élastique de la contrainte (pas d'évolution de d)
  226. DO i=1,3
  227. SIGPO(i)=XKendo*(TRACEPSPO-TRACEPSPOPT)+2*XGendo*
  228. & ((EPSPO(i)-1./3.*TRACEPSPO)-(EPSPOPT(i)-1./3.*TRACEPSPOPT))
  229. ENDDO
  230. DO i=4,6
  231. SIGPO(i)=2*XGendo*((EPSPO(i))
  232. & -(EPSPOPT(i)))
  233. ENDDO
  234. C
  235. DO i=1,6
  236. SIGTDTE(i)=SIGT(i)+SIGPO(i)*DELTAT
  237. ENDDO
  238. C ============================================================
  239. C - Calcul des parametres avec la contrainte predite elastique
  240. C - Calcul des invariants à t
  241. XI1tdte=SIGTDTE(1)+SIGTDTE(2)+SIGTDTE(3)
  242. C
  243. XJ2tdte=0.5*((SIGTDTE(1)-SIGTDTE(2))**2
  244. & +(SIGTDTE(2)-SIGTDTE(3))**2+(SIGTDTE(3)-SIGTDTE(1))**2
  245. & +6*(SIGTDTE(4)**2+SIGTDTE(5)**2+SIGTDTE(6)**2))
  246. C
  247. C Derivees partielles de fi
  248. DFIDSIG(1)=0.5*(4*SIGTDTE(1)
  249. & -2*SIGTDTE(2)-2*SIGTDTE(3))/SIGMT**2
  250. & +Q1*Q2*ft*SINH(Q2*XI1tdte/(2*SIGMt))/SIGMt
  251. DFIDSIG(2)=0.5*(-2*SIGTDTE(1)+4*SIGTDTE(2)
  252. & -2*SIGTDTE(3))/SIGMT**2
  253. & +Q1*Q2*ft*SINH(Q2*XI1tdte/(2*SIGMt))/SIGMt
  254. DFIDSIG(3)=0.5*(-2*SIGTDTE(1)-2*SIGTDTE(2)
  255. & +4*SIGTDTE(3))/SIGMT**2
  256. & +Q1*Q2*ft*SINH(Q2*XI1tdte/(2*SIGMt))/SIGMt
  257. DFIDSIG(4)=6/(SIGMt**2)*SIGTDTE(4)
  258. DFIDSIG(5)=6/(SIGMt**2)*SIGTDTE(5)
  259. DFIDSIG(6)=6/(SIGMt**2)*SIGTDTE(6)
  260. C
  261. TRACEDFIDSIG=DFIDSIG(1)+DFIDSIG(2)+DFIDSIG(3)
  262. C
  263. DFIDSIGM=-2*XJ2tdte/(SIGMt**3)-Q1*Q2*ft*XI1tdte
  264. & *SINH(Q2*XI1tdte/(2*SIGMt))/(SIGMt**2)
  265. C
  266. DFIDft=2*Q1*COSH(Q2*XI1tdte/(2*SIGMt))-2*ft*Q3*Q3
  267. C===================================================================
  268. C - Calcul de la surface seuil fi à l'aide de la contrainte prédite
  269. Fi=XJ2tdte/(SIGMT**2)+2*Q1*fT*COSH((Q2*XI1tdte)/(2*SIGMT))
  270. & -(1+(Q3*fT)**2)
  271. C
  272. C======================================================================
  273. C========= Organigramme du calcul de la contrainte a t+Dt =============
  274. C======================================================================
  275. C
  276. C Test sur le signe de la deformation equivalente
  277.  
  278.  
  279. IF(EPSE.GT.XZERO) THEN
  280. IF(fi.LT.XZERO)THEN
  281. IF(EPSE.LT.ES) THEN
  282. C =====================================================================
  283. C - Elasticite Seule
  284. dtdt=dt
  285. ftdt=ft
  286. dttdt=dtt
  287. dctdt=dct
  288. sigmtdt=sigmt
  289. DO i=1,6
  290. sigtdt(i)=sigtdte(i)
  291. ENDDO
  292. ELSE
  293. C =====================================================================
  294. C - Visco-Endommagement Seul
  295. sigmtdt=sigmt
  296. ftdt=ft
  297. CALL VISCOE(XMAT,EPSIPP,EPSE,DTT,DCT,DT
  298. & ,DELTAT,DTTDT,DCTDT,DTDT)
  299. C
  300. Dpo=(Dtdt-Dt)/DELTAT
  301. C
  302. C - Raideur endommagee et ses derivees
  303. XKendo=(1-Dtdt)*XKporo
  304. XGendo=(1-Dtdt)*XGporo
  305. XKpo=-Dpo*XKporo
  306. XGpo=-Dpo*XGporo
  307. C
  308. C - Calcul de la contrainte
  309. DO i=1,3
  310. SIGPO(i)=XKendo*(TRACEPSPO-TRACEPSPOPT)+2*XGendo
  311. & *((EPSPO(i)-1./3.*TRACEPSPO)-(EPSPOPT(i)
  312. & -1./3.*TRACEPSPOPT))
  313. & +XKpo*(TRACEPSTDT-TRACEPSPTDT)+2*XGpo
  314. & *((EPSTDT(i)-1./3.*TRACEPSTDT)-(EPSPTDT(i)
  315. & -1./3.*TRACEPSPTDT))
  316. ENDDO
  317. DO i=4,6
  318. SIGPO(i)=2*XGendo*((EPSPO(i))
  319. & -(EPSPOPT(i)))+2*XGpo
  320. & *((EPSTDT(i))-(EPSPTDT(i)))
  321. ENDDO
  322. C
  323. DO i=1,6
  324. SIGTDT(i)=SIGT(i)+SIGPO(i)*DELTAT
  325. ENDDO
  326. ENDIF
  327. ELSE
  328. IF(EPSE.LT.ES) THEN
  329. C =====================================================================
  330. C - Plasticite sans endommagement
  331. Dtdt=Dt
  332. dttdt=dtt
  333. dctdt=dct
  334. C - Evolution de la contrainte dans la matrice
  335. EPSPOPM=-xlambdapo*DFIDSIGM
  336. SIGPOM=(Etan*Eini/(Eini-Etan))*EPSPOPM
  337. SIGMTDT=SIGMT+SIGPOM*DELTAT
  338. C - Evolution de la porosite
  339. fpo=XK*((UN-ft)*ft*TRACEPSPOPT)
  340. ftdt=ft+fpo*DELTAT
  341. IF (ftdt.LT.XPETIT) THEN
  342. ftdt=XZERO
  343. ENDIF
  344. C - Raideur endommagee et ses derives
  345. XKporo=4*XKm*XGm*(1-fTDT)/(4*XGm+3*XKm*fTDT)
  346. XGporo=XGm*(1-fTDT)/(1+fTDT*(6*XKm+12*XGm)/(9*XKm+8*XGm))
  347. XKendo=(1-Dtdt)*XKporo
  348. XGendo=(1-Dtdt)*XGporo
  349. XKpo=(1-Dtdt)*(-4*XKm*XGm*(4*XGm+3*XKm)
  350. & /(4*XGm+3*XKm*fTDT)**2)*fpo
  351. XGpo=(1-Dtdt)*(-5*XGm*(32*XGm**2+60*XGm*XKm+27*XKm**2)
  352. & /(4*(2+3*fTDT)*XGm+3*(3+2*fTDT)*XKm)**2)*fpo
  353. C =============================
  354. C - Correction de la contrainte
  355. DO i=1,3
  356. SIGPO(i)=XKendo*(TRACEPSPO-TRACEPSPOPT)+2*XGendo
  357. & *((EPSPO(i)-1./3.*TRACEPSPO)-(EPSPOPT(i)
  358. & -1./3.*TRACEPSPOPT))
  359. & +XKpo*(TRACEPSTDT-TRACEPSPTDT)+2*XGpo*
  360. & ((EPSTDT(i)-1./3.*TRACEPSTDT)-(EPSPTDT(i)
  361. & -1./3.*TRACEPSPTDT))
  362. ENDDO
  363. DO i=4,6
  364. SIGPO(i)=2*XGendo*((EPSPO(i))
  365. & -(EPSPOPT(i)))+2*XGpo*
  366. & ((EPSTDT(i))-(EPSPTDT(i)))
  367. ENDDO
  368. C
  369. DO i=1,6
  370. SIGTDT(i)=SIGT(i)+SIGPO(i)*DELTAT
  371. ENDDO
  372. ELSE
  373. C =====================================================================
  374. C - Plasticite avec endommagement
  375. C - Evolution de l'endommagement
  376. CALL VISCOE(XMAT,EPSIPP,EPSE,DTT,DCT,DT
  377. & ,DELTAT,DTTDT,DCTDT,DTDT)
  378. C
  379. Dpo=(Dtdt-Dt)/DELTAT
  380. IF (Dpo.LE.0.) THEN
  381. Dtdt=Dt
  382. Dpo=0.
  383. ENDIF
  384. C - Evolution de la contrainte dans la matrice
  385. EPSPOPM=-xlambdapo*DFIDSIGM
  386. SIGPOM=(Etan*Eini/(Eini-Etan))*EPSPOPM
  387. SIGMTDT=SIGMT+SIGPOM*DELTAT
  388. C - Evolution de la porosite
  389. fpo=XK*((1-ft)*ft*TRACEPSPOPT)
  390. ftdt=ft+fpo*DELTAT
  391. IF (ftdt.LT.0.) THEN
  392. ftdt=0.
  393. ENDIF
  394. C - Raideur endommagee et ses derives
  395. XKporo=4*XKm*XGm*(1-fTDT)/(4*XGm+3*XKm*fTDT)
  396. XGporo=XGm*(1-fTDT)/(1+fTDT*(6*XKm+12*XGm)/(9*XKm+8*XGm))
  397. XKendo=(1-dtdt)*XKporo
  398. XGendo=(1-dtdt)*XGporo
  399. XKpo=-Dpo*XKporo+(1-Dtdt)*
  400. & (-4*XKm*XGm*(4*XGm+3*XKm)/(4*XGm+3*XKm*fTDT)**2)*fpo
  401. XGpo=-Dpo*XGporo
  402. & +(1-Dtdt)*(-5*XGm*(32*XGm**2+60*XGm*XKm+27*XKm**2)
  403. & /(4*(2+3*fTDT)*XGm+3*(3+2*fTDT)*XKm)**2)*fpo
  404. C =============================
  405. C - Correction de la contrainte
  406. DO i=1,3
  407. SIGPO(i)=XKendo*(TRACEPSPO-TRACEPSPOPT)+2*XGendo
  408. & *((EPSPO(i)-1./3.*TRACEPSPO)-(EPSPOPT(i)
  409. & -1./3.*TRACEPSPOPT))
  410. & +XKpo*(TRACEPSTDT-TRACEPSPTDT)+2*XGpo*
  411. & ((EPSTDT(i)-1./3.*TRACEPSTDT)-(EPSPTDT(i)
  412. & -1./3.*TRACEPSPTDT))
  413. ENDDO
  414. DO i=4,6
  415. SIGPO(i)=2*XGendo*((EPSPO(i))
  416. & -(EPSPOPT(i)))+2*XGpo*
  417. & ((EPSTDT(i))-(EPSPTDT(i)))
  418. ENDDO
  419. C
  420. DO i=1,6
  421. SIGTDT(i)=SIGT(i)+SIGPO(i)*DELTAT
  422. ENDDO
  423. ENDIF
  424. ENDIF
  425. ELSE
  426. IF(fi.LT.XZERO)THEN
  427. C =====================================================================
  428. C - Elasticite Seule
  429. Dtdt=Dt
  430. ftdt=ft
  431. dttdt=dtt
  432. dctdt=dct
  433. sigmtdt=sigmt
  434. DO i=1,6
  435. sigtdt(i)=sigtdte(i)
  436. ENDDO
  437. ELSE
  438. C =====================================================================
  439. C - Plasticite sans endommagement
  440. Dtdt=Dt
  441. dttdt=dtt
  442. dctdt=dct
  443. C - Evolution de la contrainte dans la matrice
  444. EPSPOPM=-xlambdapo*DFIDSIGM
  445. SIGPOM=(Etan*Eini/(Eini-Etan))*EPSPOPM
  446. SIGMTDT=SIGMT+SIGPOM*DELTAT
  447. C - Evolution de la porosite
  448. fpo=XK*((1-ft)*ft*TRACEPSPOPT)
  449. ftdt=ft+fpo*DELTAT
  450. IF (ftdt.LT.0.) THEN
  451. ftdt=0.
  452. ENDIF
  453. C - Raideur endommagee et ses derives
  454. XKporo=4*XKm*XGm*(1-fTDT)/(4*XGm+3*XKm*fTDT)
  455. XGporo=XGm*(1-fTDT)/(1+fTDT*(6*XKm+12*XGm)/(9*XKm+8*XGm))
  456. XKendo=(1-Dtdt)*XKporo
  457. XGendo=(1-Dtdt)*XGporo
  458. XKpo=(1-Dtdt)*(-4*XKm*XGm*(4*XGm+3*XKm)
  459. & /(4*XGm+3*XKm*fTDT)**2)*fpo
  460. XGpo=(1-Dtdt)*(-5*XGm*(32*XGm**2+60*XGm*XKm+27*XKm**2)
  461. & /(4*(2+3*fTDT)*XGm+3*(3+2*fTDT)*XKm)**2)*fpo
  462. C =============================
  463. C - Correction de la contrainte
  464. DO i=1,3
  465. SIGPO(i)=XKendo*(TRACEPSPO-TRACEPSPOPT)+2*XGendo
  466. & *((EPSPO(i)-1./3.*TRACEPSPO)-(EPSPOPT(i)
  467. & -1./3.*TRACEPSPOPT))
  468. & +XKpo*(TRACEPSTDT-TRACEPSPTDT)+2*XGpo*
  469. & ((EPSTDT(i)-1./3.*TRACEPSTDT)-(EPSPTDT(i)
  470. & -1./3.*TRACEPSPTDT))
  471. ENDDO
  472. DO i=4,6
  473. SIGPO(i)=2*XGendo*((EPSPO(i))
  474. & -(EPSPOPT(i)))+2*XGpo*
  475. & ((EPSTDT(i))-(EPSPTDT(i)))
  476. ENDDO
  477. C
  478. DO i=1,6
  479. SIGTDT(i)=SIGT(i)+SIGPO(i)*DELTAT
  480. ENDDO
  481. ENDIF
  482. ENDIF
  483. C
  484. C======================================================================
  485. C========= Fin de l'organigramme du calcul de la contrainte a t+dt ====
  486. C======================================================================
  487. C
  488. C =====================================================================
  489. C - Evolution des variables historiques
  490. VAR0(1)=DTDT
  491. VAR0(2)=fTDT-D0
  492. VAR0(3)=EPSPTDT(1)
  493. VAR0(4)=EPSPTDT(2)
  494. VAR0(5)=EPSPTDT(3)
  495. VAR0(6)=EPSPTDT(4)
  496. VAR0(7)=EPSPTDT(5)
  497. VAR0(8)=EPSPTDT(6)
  498. VAR0(9)=SIGMTDT-SIGM0
  499. VAR0(10)=dttdt
  500. VAR0(11)=dctdt
  501. VAR0(12)=ES-ED0
  502. C =====================================================================
  503. SIGT(1)=SIGTDT(1)
  504. SIGT(2)=SIGTDT(2)
  505. SIGT(3)=SIGTDT(3)
  506. SIGT(4)=SIGTDT(4)
  507. SIGT(5)=SIGTDT(5)
  508. SIGT(6)=SIGTDT(6)
  509. C =====================================================================
  510. EPST(1)=EPSTDT(1)
  511. EPST(2)=EPSTDT(2)
  512. EPST(3)=EPSTDT(3)
  513. EPST(4)=EPSTDT(4)
  514. EPST(5)=EPSTDT(5)
  515. EPST(6)=EPSTDT(6)
  516. C =====================================================================
  517. C - Compteur du Subcycling
  518. l=l+1
  519. ENDDO
  520. C =====================================================================
  521. C ========Fin du Subcycling=================
  522. C =====================================================================
  523. C - Evolution des variables historiques
  524. VARF(1)=DTDT
  525. VARF(2)=fTDT-D0
  526. VARF(3)=EPSPTDT(1)
  527. VARF(4)=EPSPTDT(2)
  528. VARF(5)=EPSPTDT(3)
  529. VARF(6)=EPSPTDT(4)
  530. VARF(7)=EPSPTDT(5)
  531. VARF(8)=EPSPTDT(6)
  532. VARF(9)=SIGMTDT-SIGM0
  533. VARF(10)=dttdt
  534. VARF(11)=dctdt
  535. VARF(12)=ES-ED0
  536. C======================================================================
  537. C - Mise a jour des contraintes
  538. SIGF(1)=SIGTDT(1)
  539. SIGF(2)=SIGTDT(2)
  540. SIGF(3)=SIGTDT(3)
  541. SIGF(4)=SIGTDT(4)
  542. SIGF(5)=SIGTDT(5)
  543. SIGF(6)=SIGTDT(6)
  544. C
  545. END
  546.  
  547.  
  548.  
  549.  
  550.  
  551.  
  552.  
  553.  
  554.  
  555.  
  556.  

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