Télécharger cmicr1.eso

Retour à la liste

Numérotation des lignes :

  1. C CMICR1 SOURCE BP208322 17/03/01 21:15:50 9325
  2. C MICRO1 SOURCE AM 00/12/13 21:40:25 4045
  3. SUBROUTINE CMICR1(WRK52,WRK53,WRK54,NSTRS1,NVARI,IDECAL
  4. & ,LIMPR,DEFPT,EPSE1,EPSE)
  5. IMPLICIT INTEGER(I-N)
  6. -INC CCOPTIO
  7. -INC DECHE
  8. INTEGER NSTRS1,NVARI,ISTRS,ITER,NSIG,IDECAL
  9. LOGICAL CONV,LCP,LIMPR
  10. REAL*8 SEQ,TRDFDS,YOUNG,XNU,DEUXMU,LAMB,R0,DR,DR0,RINF,RI,FTR
  11. REAL*8 FPC,MP,EPSE,R,TEST,TRSIG,F0,DELTL,DELTLT,EPSE1,EPSD0
  12. REAL*8 DEPS3,EPS3,DFDELTL,DFDEPS3,DS3DELTL,DS3DEPS3
  13. REAL*8 A,B,C,A1,B1,C1,A2,B2,C2,DEFPT(6)
  14. *
  15. *
  16. SEGMENT WRKK1
  17. REAL*8 SIGE0(NSTRS1),DFDS(NSTRS1)
  18. END SEGMENT
  19. *
  20. SEGINI WRKK1
  21. YOUNG=XMAT(1)
  22. XNU=XMAT(2)
  23. FPC=XMAT(7)
  24. MP=XMAT(8)
  25. EPSD0=XMAT(5)
  26. EPSE=VAR0(1)
  27. RI=VAR0(2)
  28. R0=EPSD0*YOUNG
  29. * R0=0.4D0*FPC
  30. RINF=0.6d0*FPC
  31. LCP=IFOUR.EQ.-2
  32. DEUXMU=YOUNG/(1.D0+XNU)
  33. LAMB=XNU*DEUXMU/(1.D0-2.D0*XNU)
  34. TEST=1.D-6*R0
  35. DR=0.D0
  36. KERRE=0
  37.  
  38.  
  39. C
  40. C calcul de la matrice elastique
  41. C
  42. CMATE = 'ISOTROPE'
  43. KCAS=1
  44. CALL DOHMAS(XMAT,CMATE,IFOUR,NSTRS1,KCAS,DDHOOK,IRTD)
  45. IF ( IRTD .EQ. 0) THEN
  46. KERRE=123
  47. return
  48. end if
  49. C
  50. C calcul de l'increment de contrainte
  51. C
  52. CALL MATVE1 (DDHOOK,DEPST,NSTRS1,NSTRS1,DSIGT,1)
  53. * print*,'IFOUR=',IFOUR
  54. * print*,'DDHOOK'
  55. * do i=1,3
  56. * print*,(ddhook(i,j),j=1,3)
  57. * end do
  58. * do i=1,3
  59. * print*,'DEPST(',I,')=',DEPST(i)
  60. * print*,'DSIGT(',I,')=',DSIGT(I)
  61. * end do
  62. *
  63. * calcul des parametres du critere de NADAI
  64. *
  65. C-
  66. C- A,B,C : PARAMETRES DU CRITERE DE PLASTICITE (NADAI)
  67. C- ILS SONT DIFFERENTS SUIVANT LE CHARGEMENT
  68. C-
  69.  
  70. CONV=.FALSE.
  71. * CONV=.TRUE.
  72. CALL ZDANUL (DEFP,NSTRS1)
  73. ITER=0
  74. EPS3=0.d0
  75. DELTLT=0.D0
  76. R=RI
  77. A1=0.4d0
  78. B1=(1.d0 - A1) / (1.D0 + A1)
  79. C1=(2.D0 * A1) / (1.D0 + A1)
  80. A2=1.16d0
  81. * B2=(A2 - 1.D0) / (2.D0 * A2 - 1.D0)
  82. B2=0.2
  83. * C2=A2 / (2.D0 * A2 - 1.D0)
  84. C2=0.8
  85. *
  86. * NSIG est la dimension du tenseur des contraintes:
  87. * NSIG=3 sauf 2 en contraintes planes
  88. *
  89. IF (LCP)THEN
  90. NSIG=2
  91. ELSE
  92. NSIG=3
  93. END IF
  94. *
  95. * on recupere les contraintes anelastiques due a l'endommagement
  96. *
  97. DO ISTRS=1,NSTRS1
  98. SIGE0(ISTRS)=VAR0(ISTRS+IDECAL)
  99. END DO
  100. *
  101. * prediction elastique de la contrainte effective
  102. * contrainte + contrainte anelastique + increment de contrainte elastique
  103. *
  104. * l'increment de contrainte est decoupe:
  105. * 1) point de contact avec la surface seuil ???
  106. * 2) decoupage du reste de facon a ce que le pas en contrainte
  107. * equivalente soit inferieur a disons pour le moment E*1.e-04
  108. *
  109. * en deux coups les gros
  110. *
  111. *
  112. * DS=0.d0
  113. * do istrs=1,nstrs
  114. * ds=ds + dsigt(istrs)**2
  115. * end do
  116. * nstep=int(1.e+4*sqrt(ds)/youn)
  117. DO ISTRS=1,NSTRS1
  118. SIGE0(ISTRS)=SIG0(ISTRS)+SIGE0(ISTRS)+DSIGT(ISTRS)
  119. END DO
  120. if(limpr)then
  121. DO ISTRS=1,NSTRS1
  122. PRINT*,'SIG0(',ISTRS,')=',SIGE0(ISTRS)
  123. PRINT*,'DSIGT(',ISTRS,')=',DSIGT(ISTRS)
  124. END DO
  125. endif
  126.  
  127. DO WHILE (ITER .LE. 100 .AND. .NOT. CONV)
  128. ITER=ITER+1
  129. *
  130. * on calcule la contrainte equivalente de Mises (SEQ)
  131. * et la trace des contraintes (TRSIG)
  132. *
  133. * on utilise la fonction modulo:
  134. * MOD( 1,3)= 1
  135. * MOD( 2,3)= 2
  136. * MOD( 3,3)= 0
  137. * MOD( 4,3)= 1
  138. * MOD( 5,3)= 2
  139. * MOD( 6,3)= 0
  140. *
  141. SEQ=0.D0
  142. TRSIG=0.d0
  143. * ici NSIG change en 3 le 6/11/95 (clb)
  144. DO ISTRS=1,3
  145. TRSIG=TRSIG + SIGE0(ISTRS)
  146. SEQ=SEQ+SIGE0(ISTRS)*(SIGE0(ISTRS)-SIGE0(MOD(ISTRS,3)+1))
  147. END DO
  148. DO ISTRS=4,NSTRS1
  149. SEQ=SEQ + 3.D0*SIGE0(ISTRS)**2
  150. END DO
  151. if (limpr) then
  152. print*,'VMIS**2=',SEQ,'TRSIG=',TRSIG
  153. end if
  154. SEQ=SQRT(SEQ)
  155. *
  156. * en fonction de la valeur de la trace de sigma, on determine les
  157. * coef b et c a applique a la surface seuil
  158. *
  159. IF (MOD(ITER-1,3).EQ.0 )then
  160. IF ( TRSIG .GT. (4.D0*(R+R0)*(C2-C1)/(B2-B1)))THEN
  161. * PRINT*,'OUI TRSIG .GT. (R+R0)'
  162. * PRINT*,'(R+R0)*(C2-C1)/(B2-B1)',4.D0*(R+R0)*(C2-C1)/(B2-B1)
  163. * B=B1
  164. * C=C1
  165. B=1.D0
  166. C=2.D0
  167. * C=1.D0
  168. * B=((((0.4D0*FPC)/FTR)*C)-1)
  169. * PRINT*,'B',B
  170. * PRINT*,'C',C
  171. ELSE
  172. * PRINT*,'NON TRSIG .GT. (R+R0)'
  173. * PRINT*,'(R+R0)*(C2-C1)/(B2-B1)',4*(R+R0)*(C2-C1)/(B2-B1)
  174. * B=B2
  175. * C=C2
  176. * B=-1.D0
  177. * C=8.0D0
  178. B=0.25D0
  179. C=3.0D0
  180. * PRINT*,'B',B
  181. * PRINT*,'C',C
  182. ENDIF
  183. ENDIF
  184. * B=1.D0
  185. * C=2.D0
  186. *
  187. * fonction seuil
  188. *
  189. F0=(SEQ + B*TRSIG)/C - (R0 + R)
  190. if (limpr) then
  191. PRINT*,'B=',B,'C=',C,'F0=',F0
  192. PRINT*,'SIGMAEQ=',(SEQ + B*TRSIG)/C
  193. PRINT*,'R0',R0
  194. PRINT*,'(1.D0+B)/C',(1.D0+B)/C
  195. PRINT*,'B', B
  196. PRINT*,'C', C
  197. PRINT*,'TRSIG',TRSIG,'SEQ',SEQ
  198. PRINT*,'F0',F0
  199. end if
  200. *
  201. * le critere de convergence est verifie si la valeur absolue
  202. * de la fonction seuil est plus petite que le critere ou si
  203. * la fonction seuil est negative pour la premiere iteration
  204. *
  205. CONV=(ABS(F0).LE.TEST .AND.
  206. 1 (.NOT. LCP .OR. ABS(SIGE0(3)) .LE. TEST ) )
  207. 2 .OR. (F0 .LE. 0.D0 .AND. ITER .EQ. 1)
  208. *
  209. * dans le cas des contraintes planes on rajoute aussi
  210. * que la valeur absolue de sigma33 est inferieure au critere
  211. *
  212. * PRINT*,'CONV',CONV
  213. * PRINT*,'F(0)',F0
  214. IF (.NOT. CONV) THEN
  215. * PRINT*,'On ecoule'
  216. *
  217. * calcul de dF/dSigma
  218. *
  219. DO ISTRS=1,3
  220. DFDS(ISTRS)=(2.D0*SIGE0(ISTRS)-SIGE0(MOD(ISTRS,3)+1)
  221. & -SIGE0(MOD(ISTRS+1,3)+1) )/(2.D0*C*SEQ) +B/C
  222. END DO
  223. * PRINT*,'DFDS(11)',DFDS(1)
  224. * PRINT*,'DFDS(22)',DFDS(2)
  225. * PRINT*,'DFDS(33)',DFDS(3)
  226. DO ISTRS=4,NSTRS1
  227. DFDS(ISTRS)=3.D0*SIGE0(ISTRS)/(2.d0*C*SEQ)
  228. END DO
  229. TRDFDS=0.D0
  230. DO ISTRS=1,NSIG
  231. TRDFDS=TRDFDS+DFDS(ISTRS)
  232. END DO
  233. if (limpr) then
  234. * on verifie que les valeur ppales sont bonnes
  235. do istrs=1,nstrs
  236. print*,'dfds(',istrs,')=',dfds(istrs)
  237. end do
  238. print*,'trdfds=',trdfds
  239. end if
  240. *
  241. * calcul de dsigma/deltl=C * depsp
  242. * attention !! dans le cas des contraintes planes:epsp33 n'est
  243. * pas une inconnue directe mais c'est depse33
  244. *
  245. *
  246. * R=MP*p -----> DR=MP * Dp =MP*DELTL
  247. * R=Rinf(1-exp(-bR)) ---> dR=b(Rinf-R)deltl
  248. *
  249. * DR = (RINF - R)*MP
  250. * DR=MP
  251. * R=MP*(p)**0.5 ---> DR=0.5 * MP *P**(-05)dp
  252. * Si On prend un écrouissage linéaire
  253. * R=MP*(p)**1.D0 ---> DR= MP dp
  254. *
  255.  
  256. * IF ((EPSE+DELTLT).GE. 1.D-20 )THEN
  257. * DR=0.5*MP*(EPSE+DELTLT)**(-0.5)
  258. DR=MP
  259. * PRINT*,'OUI (EPSE+DELTLT).GE. 1.D-20 )',DR
  260. * ELSE
  261. * DR=0.D0
  262. * END IF
  263. DFDELTL=-DR
  264. DO ISTRS=1,NSIG
  265. DFDELTL=DFDELTL-DFDS(ISTRS)*(DFDS(ISTRS)*DEUXMU+TRDFDS*LAMB)
  266. END DO
  267. DO ISTRS=4,NSTRS1
  268. DFDELTL=DFDELTL - 2.d0*DEUXMU * DFDS(ISTRS)**2
  269. END DO
  270. IF (LCP) THEN
  271. DFDELTL=DFDELTL-DFDS(3)*TRDFDS*LAMB
  272. DFDEPS3=LAMB*TRDFDS + (LAMB+DEUXMU)*DFDS(3)
  273. DS3DELTL=-LAMB*TRDFDS
  274. DS3DEPS3=LAMB + DEUXMU
  275. DELTL=(SIGE0(3)*DFDEPS3 - F0*DS3DEPS3)/
  276. & (DS3DEPS3*DFDELTL - DS3DELTL*DFDEPS3)
  277. DEPS3=-(SIGE0(3)-DS3DELTL*DELTL)/DS3DEPS3
  278. EPS3=EPS3+DEPS3
  279. DELTLT=DELTLT+DELTL
  280. DO ISTRS=1,NSIG
  281. SIGE0(ISTRS)=SIGE0(ISTRS)+LAMB*DEPS3
  282. END DO
  283. SIGE0(3)=SIGE0(3)+(LAMB+DEUXMU)*DEPS3-LAMB*DELTL*TRDFDS
  284. ELSE
  285. DELTL=-F0/DFDELTL
  286. DELTLT=DELTLT+DELTL
  287. END IF
  288. DO ISTRS=1,NSIG
  289. SIGE0(ISTRS)=SIGE0(ISTRS)
  290. & -DELTL*(DEUXMU*DFDS(ISTRS)+LAMB*TRDFDS)
  291. END DO
  292. DO ISTRS=4,NSTRS1
  293. SIGE0(ISTRS)=SIGE0(ISTRS)-DELTL*DEUXMU*DFDS(ISTRS)
  294. END DO
  295. ELSE
  296. DELTL=0.D0
  297. END IF
  298. R=R + DR*DELTL
  299. DO ISTRS=1,3
  300. DEFP(ISTRS)=DEFP(ISTRS)+DELTL*DFDS(ISTRS)
  301. END DO
  302. DO ISTRS=4,NSTRS1
  303. DEFP(ISTRS)=DEFP(ISTRS)+2.D0*DELTL*DFDS(ISTRS)
  304. END DO
  305. if (limpr) then
  306. print*,'ITER=',ITER,'F0=',F0,'DELTL=', DELTL
  307. print*,'DELTLT=',DELTLT,'DR=',DR,'R',R
  308. do istrs=1,nstrs
  309. print*,'SIGMA(',ISTRS,')=', SIGE0(ISTRS)
  310. print*,'DEFP(',ISTRS,')=', DEFP(ISTRS)
  311. end do
  312. end if
  313. * PRINT*,'EPSE+DELTLT',EPSE+DELTLT
  314. END DO
  315. * EN CAS DE NON CONVERGENCE: ON MET KERRE a 1
  316. IF (.NOT. CONV) THEN
  317. KERRE=1
  318. END IF
  319. * PRINT*,'CONVERGENCE OBTENUE APRES ',ITER, 'ITERATIONS'
  320. *
  321. * on remplit les variables finales
  322.  
  323. **
  324. * Je rajoute ça pour les déformations plastiques
  325. *
  326. DO ISTRS=1,NSTRS1
  327. DEFPT(ISTRS)=DEFP(ISTRS)
  328. ENDDO
  329. VARF(1)=EPSE+DELTLT
  330. VARF(2)=R
  331. EPSE1=EPSE*((1.D0+B)/C)
  332. * EPSE=EPSE1*(C/(1.D0+B))
  333. *
  334. * ON REND LES CONTRAINTES
  335. *
  336. DO ISTRS=1,NSTRS1
  337. SIGF(ISTRS)=SIGE0(ISTRS)
  338. * PRINT*,'SIGF(',ISTRS,')=',SIGF(ISTRS)
  339. END DO
  340.  
  341. * ON REND L'INCREMENT DE DEFORMATIONS PLASTIQUES
  342. *
  343.  
  344. SEGSUP WRKK1
  345. RETURN
  346. END
  347. *
  348.  
  349.  
  350.  
  351.  
  352.  
  353.  
  354.  
  355.  
  356.  
  357.  

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