Télécharger desmnl.eso

Retour à la liste

Numérotation des lignes :

  1. C DESMNL SOURCE CB215821 16/04/21 21:16:24 8920
  2. C DESMORAT
  3. SUBROUTINE DESMNL(WRK0,WRK1,WRK5,NSTRS,NVARI,LHOOK,ISTEP,
  4. & NMATT,ICARA,KERRE,MFR)
  5. C
  6. C variables en entree
  7. C
  8. C WRK0, WRK1 pointeur actif a l'entree
  9. C NVARI nombre de variables internes (doit etre egal a ?)
  10. C NMATT nombre de constantes du materiau
  11. C
  12. C ISTEP flag utilise pour separer les etapes dans un calcul non local
  13. C ISTEP=0 -----> calcul local
  14. C ISTEP=1 -----> calcul non local etape 1 on calcule les seuils
  15. C ISTEP=2 -----> calcul non local etape 2 on continue le calcul
  16. C a partir des seuils moyennes
  17. C
  18. C variables en sortie
  19. C
  20. C VARINF variables internes finales
  21. C
  22. C SIGMAF contraintes finales
  23. C
  24. C
  25. IMPLICIT REAL*8(A-H,O-Z)
  26. C
  27. -INC CCOPTIO
  28. -INC SMINTE
  29. -INC SMCOORD
  30. -INC SMELEME
  31. C
  32. C
  33. CHARACTER*8 CMATE
  34. REAL*8 XD0(3,3),DF(3,3),POSEPSI2(3,3),EPSIN(6),EPS33(3,3)
  35. REAL*8 SIGTILDE33(3,3),SIGDEV33(3,3),SIGD33(3,3),SIGTILDE(6)
  36. REAL*8 UNMOINSDVP(3),HM1(3,3),DFD(3,3),UNMOINSD(3,3),DFVP(3)
  37. REAL*8 ROT(3,3),ROTT(3,3),VECP33T(3,3),TRAV(3,3),VECP33(3,3)
  38. REAL*8 EPSIPOS(3),HM1DIAG(3),EPS(2,2),DELTAF22(2,2),EPSIPP(3)
  39. REAL*8 VPDF(2),VECPDF(2,2),EPSTF(6),DHOOK(6,6)
  40. C
  41. PARAMETER (XZERO=0.D0 , UN=1.D0 , DEUX=2.D0, XPETIT=1.D-12)
  42. C
  43. SEGMENT WRK0
  44. REAL*8 XMAT(NMATT)
  45. ENDSEGMENT
  46. C
  47. SEGMENT WRK1
  48. REAL*8 DDHOOK(LHOOK,LHOOK),SIG0(NSTRS),DEPST(NSTRS)
  49. REAL*8 SIGF(NSTRS),VAR0(NVARI),VARF(NVARI)
  50. REAL*8 DEFP(NSTRS),XCAR(ICARA)
  51. ENDSEGMENT
  52. C
  53. SEGMENT WRK5
  54. REAL*8 EPIN0(NSTRS),EPINF(NSTRS),EPST0(NSTRS)
  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 Début de l'intégration numérique 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. C mise a jour de XD0 pour traiter IFAILURE=2 dans le meme pas
  254. DO I=1,3
  255. DO J=1,3
  256. XD0(I,J)=DF(I,J)
  257. ENDDO
  258. ENDDO
  259. TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3)
  260. TRADF=XD0(1,1)+XD0(2,2)+XD0(3,3)
  261. CALL JACOB3(XD0,3,DFVP,ROT)
  262. ENDIF
  263. ENDIF
  264. C
  265. IF(IFAILURE.EQ.2) THEN
  266. C deux directions d'endom sont fixee a Dc
  267. CALL TRSPOD(ROT,3,3,ROTT)
  268. C
  269. DF33=TRADF-DEUX*Dc
  270. C print*, 'DF33=',DF33
  271. IF(DF33.LT.XZERO) THEN
  272. DF33=XPETIT
  273. ENDIF
  274. C On corrige l'endommagement
  275. CALL ZERO(DF,3,3)
  276. DF(1,1)=Dc
  277. DF(2,2)=Dc
  278. DF(3,3)=DF33
  279. C
  280. CALL MATMAT(DF,ROTT,3,3,3,TRAV)
  281. CALL MATMAT(ROT,TRAV,3,3,3,DF)
  282. C
  283. C On verifie la troisieme direction d'endommagement
  284. IF(DF33.GE.Dc)THEN
  285. IFAILURE=3
  286. C mise a jour de XD0 pour traiter IFAILURE=3 dans le même pas
  287. DO I=1,3
  288. DO J=1,3
  289. XD0(I,J)=DF(I,J)
  290. ENDDO
  291. ENDDO
  292. TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3)
  293. CALL JACOB3(XD0,3,DFVP,ROT)
  294. ENDIF
  295. ENDIF
  296. C
  297. IF(IFAILURE.EQ.3)THEN
  298. C On est totalement rompu, endommagement isotrope egal a DC
  299. CALL TRSPOD(ROT,3,3,ROTT)
  300. C On corrige l'endommagement
  301. CALL ZERO(DF,3,3)
  302. DF(1,1)=Dc
  303. DF(2,2)=Dc
  304. DF(3,3)=Dc
  305. C
  306. CALL MATMAT(DF,ROTT,3,3,3,TRAV)
  307. CALL MATMAT(ROT,TRAV,3,3,3,DF)
  308. ENDIF
  309. C
  310. C Fin de la correction d'endommagement
  311. C
  312. C
  313. DH=(DF(1,1)+DF(2,2)+DF(3,3))/3.D0
  314. C
  315. ENDIF
  316. C
  317. C Fin des differents cas (elasticite ou endommagement)
  318. C
  319. C Calcul de la contrainte finale
  320. C
  321. DO I=1,3
  322. DO J=1,3
  323. UNMOINSD(I,J)=-UN*DF(I,J)
  324. ENDDO
  325. ENDDO
  326. UNMOINSD(1,1)=UN-DF(1,1)
  327. UNMOINSD(2,2)=UN-DF(2,2)
  328. UNMOINSD(3,3)=UN-DF(3,3)
  329. CALL JACOB3(UNMOINSD,3,UNMOINSDVP,VECP33)
  330. HM1DIAG(1)=SQRT(UNMOINSDVP(1))
  331. HM1DIAG(2)=SQRT(UNMOINSDVP(2))
  332. HM1DIAG(3)=SQRT(UNMOINSDVP(3))
  333. CALL PRODT2(HM1,HM1DIAG,VECP33,3)
  334. C
  335. C contrainte effective
  336. CALL MATVE1(DDHOOK,EPSTF,NSTRS,NSTRS,SIGTILDE,2)
  337. CALL ENDOCB(SIGTILDE,SIGTILDE33,1,IFOUR)
  338. CALL MATMAT(HM1,SIGTILDE33,3,3,3,TRAV)
  339. CALL MATMAT(TRAV,HM1,3,3,3,SIGD33)
  340. XTRAV=XZERO
  341. DO I=1,3
  342. DO J=1,3
  343. XTRAV=XTRAV+UNMOINSD(I,J)*SIGTILDE33(I,J)
  344. ENDDO
  345. ENDDO
  346. DO I=1,3
  347. DO J=1,3
  348. SIGD33(I,J)=SIGD33(I,J)-(XTRAV/(3.D0*(UN-DH)))*UNMOINSD(I,J)
  349. ENDDO
  350. ENDDO
  351. TRASIGTILDE=SIGTILDE(1)+SIGTILDE(2)+SIGTILDE(3)
  352. C
  353. IF(TRASIGTILDE.GE.XZERO) THEN
  354. COEF=UN-etaT*DH
  355. IF(COEF.LE.XZERO) THEN
  356. COEF=UN-Dc
  357. ENDIF
  358. TRASIG=COEF*TRASIGTILDE
  359. ELSE
  360. TRASIG=(UN-etaC*DH)*TRASIGTILDE
  361. ENDIF
  362. C
  363. SIGF(1)=SIGD33(1,1)+TRASIG/3.d0
  364. SIGF(2)=SIGD33(2,2)+TRASIG/3.d0
  365. SIGF(3)=SIGD33(3,3)+TRASIG/3.d0
  366. SIGF(4)=SIGD33(1,2)
  367. SIGF(5)=SIGD33(1,3)
  368. SIGF(6)=SIGD33(2,3)
  369. C
  370. C
  371. C stockage des variables du probleme
  372. C
  373. VARF(2)=DF(1,1)
  374. VARF(3)=DF(2,2)
  375. VARF(4)=DF(3,3)
  376. VARF(5)=DF(1,2)
  377. VARF(6)=DF(1,3)
  378. VARF(7)=DF(2,3)
  379. VARF(8)=ROT(1,1)
  380. VARF(9)=ROT(1,2)
  381. VARF(10)=ROT(1,3)
  382. VARF(11)=ROT(2,1)
  383. VARF(12)=ROT(2,2)
  384. VARF(13)=ROT(2,3)
  385. VARF(14)=ROT(3,1)
  386. VARF(15)=ROT(3,2)
  387. VARF(16)=ROT(3,3)
  388. C
  389. C flag pour l'endommagement maximum
  390. VARF(17)=IFAILURE
  391. C
  392. C Deformations totales finales
  393. DO ISTRS=1,NSTRS
  394. VARF(17+ISTRS)=EPSTF(ISTRS)
  395. ENDDO
  396. C
  397. C
  398. 2000 CONTINUE
  399. C
  400. RETURN
  401. END
  402. C
  403.  
  404.  
  405.  
  406.  
  407.  
  408.  

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