Télécharger micro1.eso

Retour à la liste

Numérotation des lignes :

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

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