Télécharger desmor.eso

Retour à la liste

Numérotation des lignes :

  1. C DESMOR SOURCE BP208322 17/03/01 21:17:01 9325
  2. C DESMORAT (V2)
  3. SUBROUTINE DESMOR(WRK52,WRK53,WRK54,NVARI,Iecou)
  4. C
  5. C variables en entree
  6. C
  7. C WRK52, WRK53, WRK54 pointeur actif a l'entree
  8. C NVARI nombre de variables internes (doit etre egal a ?)
  9. C
  10. C ISTEP flag utilise pour separer les etapes dans un calcul non local
  11. C ISTEP=0 -----> calcul local
  12. C ISTEP=1 -----> calcul non local etape 1 on calcule les seuils
  13. C ISTEP=2 -----> calcul non local etape 2 on continue le calcul
  14. C a partir des seuils moyennes
  15. C
  16. C variables en sortie
  17. C
  18. C VARINF variables internes finales
  19. C
  20. C SIGMAF contraintes finales
  21. C
  22. C
  23. IMPLICIT REAL*8(A-H,O-Z)
  24. C
  25. -INC CCOPTIO
  26. -INC DECHE
  27. C
  28. C
  29. DIMENSION XD0(3,3),DF(3,3),POSEPSI2(3,3),EPS33(3,3)
  30. DIMENSION SIGTILDE33(3,3),SIGDEV33(3,3),SIGD33(3,3),SIGTILDE(6)
  31. DIMENSION UNMOINSDVP(3),HM1(3,3),DFD(3,3),UNMOINSD(3,3),DFVP(3)
  32. DIMENSION ROT(3,3),ROTT(3,3),VECP33T(3,3),TRAV(3,3),VECP33(3,3)
  33. DIMENSION EPSIPOS(3),HM1DIAG(3),EPS(2,2),DELTAF22(2,2),EPSIPP(3)
  34. DIMENSION VPDF(2),VECPDF(2,2),EPSTF(6)
  35. C
  36. PARAMETER (XZERO=0.D0 , UN=1.D0 , DEUX=2.D0, XPETIT=1.D-12)
  37. C
  38. SEGMENT IECOU
  39. * COMMON/IECOU/NYOG,NYNU,NYALFA,NYSMAX,NYN,NYM,NYKK,
  40. INTEGER icow1,icow2,icow3,icow4,icow5,icow6,icow7,
  41. C INTEGER NYOG, NYNU, NYALFA,NYSMAX,NYN, NYM, NYKK,
  42. 1 icow8,icow9,icow10,icow11,icow12,icow13,icow14,icow15,icow16,
  43. C . NYALF1,NYBET1,NYR, NYA, NYRHO,NSIGY, NNKX, NYKX, IND,
  44. 2 icow17,icow18,icow19,icow20,icow21,icow22,icow23,icow24,
  45. C . NSOM, NINV, NINCMA,NCOMP, JELEM, LEGAUS,INAT, NCXMAT,
  46. 3 icow25,icow26,icow27,icow28,icow29,icow30,ICARA,
  47. C . LTRAC, MFR, IELE, NHRM, NBNN, NBELEM,ICARA,
  48. 4 icow32,icow33,NSTRS1,MFR1,icow36,icow37,icow38,
  49. C . LW2, NDEF, NSTRSS,MFR1, NBGMAT,NELMAT,MSOUPA,
  50. 5 icow39,icow40,icow41,icow42,icow43,icow44
  51. C . NUMAT1,LENDO, NBBB, NNVARI,KERR1, MELEME
  52. INTEGER icow45,icow46,icow47,icow48,icow49,icow50,
  53. . icow51,icow52,icow53,icow54,icow55,icow56
  54. . icow57,icow58
  55. ENDSEGMENT
  56. C
  57. C
  58. C Recuperation des variables initiales dans les tableaux
  59. C
  60. CMATE= 'ISOTROPE'
  61. YOUN = XMAT(1)
  62. XNU = XMAT(2)
  63. XK0 = XMAT(5)
  64. AA = XMAT(6)
  65. A = XMAT(7)
  66. etaC = XMAT(8)
  67. etaT = XMAT(9)
  68. Xk = XMAT(10)
  69. DC=XMAT(11)
  70. C
  71. C Recuperation de l'endommagement initial
  72. C
  73. XD0(1,1)=VAR0(2)
  74. XD0(2,2)=VAR0(3)
  75. XD0(3,3)=VAR0(4)
  76. XD0(1,2)=VAR0(5)
  77. XD0(1,3)=VAR0(6)
  78. XD0(2,3)=VAR0(7)
  79. XD0(2,1)=XD0(1,2)
  80. XD0(3,1)=XD0(1,3)
  81. XD0(3,2)=XD0(2,3)
  82. C
  83. C Recuperation de la base d'endommagement
  84. C
  85. ROT(1,1)=VAR0(8)
  86. ROT(1,2)=VAR0(9)
  87. ROT(1,3)=VAR0(10)
  88. ROT(2,1)=VAR0(11)
  89. ROT(2,2)=VAR0(12)
  90. ROT(2,3)=VAR0(13)
  91. ROT(3,1)=VAR0(14)
  92. ROT(3,2)=VAR0(15)
  93. ROT(3,3)=VAR0(16)
  94. C
  95. C Recuperation de l'etat d'endommagement critique
  96. C
  97. IFAILURE=nint(VAR0(17))
  98. C
  99. C Matrice de raideur elastique
  100. CALL ZERO(DDHOOK,6,6)
  101. AUX0=YOUN/((UN+XNU)*(UN-XNU-XNU))
  102. AUX=AUX0*(UN-XNU)
  103. AUX1=AUX0*XNU
  104. GEGE=0.5D0*YOUN/(UN+XNU)
  105. DDHOOK(1,1)=AUX
  106. DDHOOK(2,1)=AUX1
  107. DDHOOK(3,1)=AUX1
  108. DDHOOK(1,2)=AUX1
  109. DDHOOK(2,2)=AUX
  110. DDHOOK(3,2)=AUX1
  111. DDHOOK(1,3)=AUX1
  112. DDHOOK(2,3)=AUX1
  113. DDHOOK(3,3)=AUX
  114. DDHOOK(4,4)=GEGE
  115. DDHOOK(5,5)=GEGE
  116. DDHOOK(6,6)=GEGE
  117. C
  118. C
  119. C Deformation totale
  120. C
  121. DO ISTRS=1,NSTRS
  122. EPSTF(ISTRS)=VAR0(17+ISTRS)+DEPST(ISTRS)
  123. ENDDO
  124. C - Trace de la deformation
  125. TRAEPS=EPSTF(1)+EPSTF(2)+EPSTF(3)
  126. C
  127. C calcul des deformations principales
  128. CALL ENDOCB(EPSTF,EPS33,2,IFOUR)
  129. CALL JACOB3(EPS33,3,EPSIPP,VECP33)
  130. EPSIPOS(1)= MAX( XZERO , EPSIPP(1) )**2
  131. EPSIPOS(2)= MAX( XZERO , EPSIPP(2) )**2
  132. EPSIPOS(3)= MAX( XZERO , EPSIPP(3) )**2
  133. EPSTILNLOC=SQRT(EPSIPOS(1)+EPSIPOS(2)+EPSIPOS(3))
  134. C
  135. IF (ISTEP.EQ.0) THEN
  136. EPSTIL=EPSTILNLOC
  137. VARF(1)=EPSTIL
  138. ELSE IF (ISTEP.EQ.1) THEN
  139. VARF(1)=EPSTILNLOC
  140. DO I=2,23
  141. VARF(I)=VAR0(I)
  142. ENDDO
  143. GOTO 2000
  144. ELSE IF (ISTEP.EQ.2) THEN
  145. c Over nonlocal, coef 2
  146. C EPSTIL=DEUX*VAR0(1)-EPSTILNLOC
  147. EPSTIL=VAR0(1)
  148. ENDIF
  149. C
  150. C**********************************************************
  151. C Debut de l'intégration numerique du comportement
  152. C**********************************************************
  153. C
  154. C calcul des seuils
  155. C
  156. C partie positive de la matrice des deformations
  157. CALL PRODT2(POSEPSI2,EPSIPOS,VECP33,3)
  158. C
  159. C calcul de la trace de D au pas n
  160. TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3)
  161. C
  162. C calcul du seuil
  163. XKAPPA=A*TAN(TRAD0/(A*AA)+ATAN(XK0/A))
  164. FCRIT=EPSTIL+Xk*TRAEPS-XKAPPA
  165. C
  166. C TEST sur le seuil
  167. IF(FCRIT.LT.XZERO) THEN
  168. C on est elastique
  169. C endommagement et deformations permanentes inchanges
  170. DO I=1,3
  171. DO J=1,3
  172. DF(I,J)=XD0(I,J)
  173. ENDDO
  174. ENDDO
  175. DH=(DF(1,1)+DF(2,2)+DF(3,3))/3.D0
  176. ELSE
  177. C
  178. C debut de la correction d'endom.
  179. C
  180. C Calcul de Kappa^(-1)=Tr Dn+1
  181. C
  182. TRADF=A*AA*(ATAN((EPSTIL+Xk*TRAEPS)/A)-ATAN(XK0/A))
  183. C
  184. C test sur l'eventuelle presence d'une direction d'endommagement fixee
  185. C
  186. IF(IFAILURE.EQ.0)THEN
  187. C actualisation de l'endommagement, Dc pas encore atteint (cas "normal")
  188. IF (EPSTIL.EQ.XZERO) THEN
  189. DLAMBDA=XZERO
  190. ELSE
  191. DLAMBDA=(TRADF-TRAD0)/(EPSTIL**2)
  192. ENDIF
  193. IF(DLAMBDA.LT.XZERO) THEN
  194. DLAMBDA=XZERO
  195. ENDIF
  196. CALL ZERO(DF,3,3)
  197. DO I=1,3
  198. DO J=1,3
  199. DF(I,J)=XD0(I,J)+DLAMBDA*POSEPSI2(I,J)
  200. ENDDO
  201. ENDDO
  202. C On borne l'endommagement D(1) a Dc si D>Dc et Fixed-crack
  203. CALL JACOD3(DF,3,DFVP)
  204. IF(DFVP(1).GE.DC) THEN
  205. IFAILURE=1
  206. C mise a jour de XD0 pour traiter IFAILURE=1 dans le meme pas
  207. DO I=1,3
  208. DO J=1,3
  209. XD0(I,J)=DF(I,J)
  210. ENDDO
  211. ENDDO
  212. TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3)
  213. TRADF=XD0(1,1)+XD0(2,2)+XD0(3,3)
  214. CALL JACOB3(XD0,3,DFVP,ROT)
  215. ENDIF
  216. ENDIF
  217. C
  218. IF(IFAILURE.EQ.1)THEN
  219. C une direction d'endom est fixee a Dc
  220. C
  221. C On calcule le multiplicateur d'endommagement dans la bonne base
  222. C On met les deformations positives dans le repère principal
  223. C de l'endommagement
  224. CALL JACOD3(XD0,3,DFVP)
  225. CALL TRSPOD(ROT,3,3,ROTT)
  226. CALL MATMAT(POSEPSI2,ROT,3,3,3,TRAV)
  227. CALL MATMAT(ROTT,TRAV,3,3,3,POSEPSI2)
  228. C
  229. IF(ABS(POSEPSI2(2,2)+POSEPSI2(3,3)).GT.XPETIT) THEN
  230. DLAMBDA=(TRADF-TRAD0)/(POSEPSI2(2,2)+POSEPSI2(3,3))
  231. ELSE
  232. DLAMBDA=XZERO
  233. ENDIF
  234. IF(DLAMBDA.LT.XZERO) THEN
  235. DLAMBDA=XZERO
  236. ENDIF
  237. C
  238. C On corrige l'endommagement
  239. CALL ZERO(DF,3,3)
  240. DF(1,1)=Dc
  241. DF(2,2)=DFVP(2)+DLAMBDA*POSEPSI2(2,2)
  242. DF(3,3)=DFVP(3)+DLAMBDA*POSEPSI2(3,3)
  243. DF(2,3)=DLAMBDA*POSEPSI2(2,3)
  244. DF(3,2)=DF(2,3)
  245. C
  246. CALL MATMAT(DF,ROTT,3,3,3,TRAV)
  247. CALL MATMAT(ROT,TRAV,3,3,3,DF)
  248. C
  249. C On verifie la deuxieme direction d'endommagement
  250. CALL JACOD3(DF,3,DFVP)
  251. IF(DFVP(2).GE.Dc)THEN
  252. IFAILURE=2
  253. PRINT*, 'IFAILURE=', IFAILURE
  254. C mise a jour de XD0 pour traiter IFAILURE=2 dans le meme pas
  255. DO I=1,3
  256. DO J=1,3
  257. XD0(I,J)=DF(I,J)
  258. ENDDO
  259. ENDDO
  260. TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3)
  261. TRADF=XD0(1,1)+XD0(2,2)+XD0(3,3)
  262. CALL JACOB3(XD0,3,DFVP,ROT)
  263. ENDIF
  264. ENDIF
  265. C
  266. IF(IFAILURE.EQ.2) THEN
  267. C deux directions d'endom sont fixee a Dc
  268. CALL TRSPOD(ROT,3,3,ROTT)
  269. C
  270. DF33=TRADF-DEUX*Dc
  271. C print*, 'DF33=',DF33
  272. IF(DF33.LT.XZERO) THEN
  273. DF33=XPETIT
  274. ENDIF
  275. C On corrige l'endommagement
  276. CALL ZERO(DF,3,3)
  277. DF(1,1)=Dc
  278. DF(2,2)=Dc
  279. DF(3,3)=DF33
  280. C
  281. CALL MATMAT(DF,ROTT,3,3,3,TRAV)
  282. CALL MATMAT(ROT,TRAV,3,3,3,DF)
  283. C
  284. C On verifie la troisieme direction d'endommagement
  285. IF(DF33.GE.Dc)THEN
  286. IFAILURE=3
  287. C mise a jour de XD0 pour traiter IFAILURE=3 dans le même pas
  288. DO I=1,3
  289. DO J=1,3
  290. XD0(I,J)=DF(I,J)
  291. ENDDO
  292. ENDDO
  293. TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3)
  294. CALL JACOB3(XD0,3,DFVP,ROT)
  295. ENDIF
  296. ENDIF
  297. C
  298. IF(IFAILURE.EQ.3)THEN
  299. C On est totalement rompu, endommagement isotrope egal a DC
  300. CALL TRSPOD(ROT,3,3,ROTT)
  301. C On corrige l'endommagement
  302. CALL ZERO(DF,3,3)
  303. DF(1,1)=Dc
  304. DF(2,2)=Dc
  305. DF(3,3)=Dc
  306. C
  307. CALL MATMAT(DF,ROTT,3,3,3,TRAV)
  308. CALL MATMAT(ROT,TRAV,3,3,3,DF)
  309. ENDIF
  310. C
  311. C Fin de la correction d'endommagement
  312. C
  313. C
  314. DH=(DF(1,1)+DF(2,2)+DF(3,3))/3.D0
  315. C
  316. ENDIF
  317. C
  318. C Fin des differents cas (elasticite ou endommagement)
  319. C
  320. C Calcul de la contrainte finale
  321. C
  322. DO I=1,3
  323. DO J=1,3
  324. UNMOINSD(I,J)=-UN*DF(I,J)
  325. ENDDO
  326. ENDDO
  327. UNMOINSD(1,1)=UN-DF(1,1)
  328. UNMOINSD(2,2)=UN-DF(2,2)
  329. UNMOINSD(3,3)=UN-DF(3,3)
  330. CALL JACOB3(UNMOINSD,3,UNMOINSDVP,VECP33)
  331. HM1DIAG(1)=SQRT(UNMOINSDVP(1))
  332. HM1DIAG(2)=SQRT(UNMOINSDVP(2))
  333. HM1DIAG(3)=SQRT(UNMOINSDVP(3))
  334. CALL PRODT2(HM1,HM1DIAG,VECP33,3)
  335. C
  336. C contrainte effective
  337. CALL MATVE1(DDHOOK,EPSTF,NSTRS,NSTRS,SIGTILDE,2)
  338. CALL ENDOCB(SIGTILDE,SIGTILDE33,1,IFOUR)
  339. CALL MATMAT(HM1,SIGTILDE33,3,3,3,TRAV)
  340. CALL MATMAT(TRAV,HM1,3,3,3,SIGD33)
  341. XTRAV=XZERO
  342. DO I=1,3
  343. DO J=1,3
  344. XTRAV=XTRAV+UNMOINSD(I,J)*SIGTILDE33(I,J)
  345. ENDDO
  346. ENDDO
  347. DO I=1,3
  348. DO J=1,3
  349. SIGD33(I,J)=SIGD33(I,J)-(XTRAV/(3.D0*(UN-DH)))*UNMOINSD(I,J)
  350. ENDDO
  351. ENDDO
  352. TRASIGTILDE=SIGTILDE(1)+SIGTILDE(2)+SIGTILDE(3)
  353. C
  354. IF(TRASIGTILDE.GE.XZERO) THEN
  355. COEF=UN-etaT*DH
  356. IF(COEF.LE.XZERO) THEN
  357. COEF=UN-Dc
  358. ENDIF
  359. TRASIG=COEF*TRASIGTILDE
  360. ELSE
  361. TRASIG=(UN-etaC*DH)*TRASIGTILDE
  362. ENDIF
  363. C
  364. SIGF(1)=SIGD33(1,1)+TRASIG/3.d0
  365. SIGF(2)=SIGD33(2,2)+TRASIG/3.d0
  366. SIGF(3)=SIGD33(3,3)+TRASIG/3.d0
  367. SIGF(4)=SIGD33(1,2)
  368. SIGF(5)=SIGD33(1,3)
  369. SIGF(6)=SIGD33(2,3)
  370. C
  371. C
  372. C stockage des variables du probleme
  373. C
  374. VARF(2)=DF(1,1)
  375. VARF(3)=DF(2,2)
  376. VARF(4)=DF(3,3)
  377. VARF(5)=DF(1,2)
  378. VARF(6)=DF(1,3)
  379. VARF(7)=DF(2,3)
  380. VARF(8)=ROT(1,1)
  381. VARF(9)=ROT(1,2)
  382. VARF(10)=ROT(1,3)
  383. VARF(11)=ROT(2,1)
  384. VARF(12)=ROT(2,2)
  385. VARF(13)=ROT(2,3)
  386. VARF(14)=ROT(3,1)
  387. VARF(15)=ROT(3,2)
  388. VARF(16)=ROT(3,3)
  389. C
  390. C flag pour l'endommagement maximum
  391. VARF(17)=IFAILURE
  392. C
  393. C Deformations totales finales
  394. DO ISTRS=1,NSTRS
  395. VARF(17+ISTRS)=EPSTF(ISTRS)
  396. ENDDO
  397. C
  398. 2000 CONTINUE
  399. C
  400. C
  401. RETURN
  402. END
  403. C
  404.  
  405.  
  406.  
  407.  
  408.  
  409.  
  410.  
  411.  
  412.  
  413.  

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