Télécharger mooney.eso

Retour à la liste

Numérotation des lignes :

  1. C MOONEY SOURCE CHAT 11/03/28 21:15:02 6923
  2. C=======================================================================
  3. C= MODELE HYPERELASTIQUE MOONEY-RIVLIN INCOMPRESSIBLE =
  4. C= EN GRANDES TRANSFORMATIONS =
  5. C= =
  6. C= Modele actuellement disponible uniquement en CONTRAINTE PLANE =
  7. C= (2D PLAN CONT ou 3D avec contraintes hors plan (1,2) nulles) =
  8. C= =
  9. C= Exemple d'utilisation d'un modele UMAT en grandes transformations =
  10. C= Contribution de Laurent Gornet - Ecole Centrale de Nantes (2006) =
  11. C= =
  12. C= Pour plus d'informations, voir la presentation de L. Gornet lors =
  13. C= du Club Cast3m 2006, disponible sur le site Web de Cast3m. =
  14. C= =
  15. C= Note : Actuellement en grandes deformations dans PASAPAS, le modele =
  16. C= ne peut contenir que des modèles type UMAT. On ne peut pas =
  17. C= "melanger" les derivees objectives dans Cast3m. =
  18. C=======================================================================
  19.  
  20. SUBROUTINE MOONEY (STRESS, STATEV, DDSDDE, STRAN, DSTRAN,
  21. & TIME, DTIME, TEMP, DTEMP, PREDEF, DPRED,
  22. & NDI, NSHR, NTENS, NSTATV, PROPS, NPROPS,
  23. & DFGRD0, DFGRD1, KSTEP, KINC,SSE,spd,scd,rpl,
  24. & DDSDDT,DRPLDE,DRPLDT,cmname,coords,drot,PNEWDT,
  25. & celent,NOEL, NPT, LAYER, KSPT )
  26.  
  27. IMPLICIT INTEGER(I-N)
  28. IMPLICIT REAL*8(A-H,O-Z)
  29. CHARACTER*16 CMNAME
  30.  
  31. INTEGER NDI, NSHR, NTENS, NSTATV, NPROPS, KSTEP, KINC
  32. & NOEL, NPT, LAYER, KSPT
  33. REAL*8 STRESS(NTENS), STATEV(NSTATV),SSE, SPD, SCD, RPL,
  34. & DDSDDE(NTENS,NTENS), STRAN(NTENS), DSTRAN(NTENS),
  35. & DDSDDT(NTENS), DRPLDE(NTENS), DRPLDT,
  36. & TIME(2), DTIME, TEMP, DTEMP, PREDEF(*), DPRED(*),
  37. & PROPS(NPROPS), DFGRD0(3,3), DFGRD1(3,3),
  38. & COORDS(3),DROT(3,3),PNEWDT,CELENT
  39. REAL*8 CG11,CG22,CG12,CG33, CInv11,CInv22,CInv12, DLT,
  40. & I1,I2, dWdI1,dWdI2, phyd, Coe1,Coe2,
  41. & F11,F21,F12,F22, S11,S12,S21,S22
  42. REAL*8 F33,I1B,I2B,cinv33,detF,mdetF,Jm2s3,Jm1s3
  43. REAL*8 S33,dwvdj,C1,C2
  44.  
  45. Parameter(zero=0.d0,one=1.d0,two=2.d0,three=3.d0,four=4.d0,
  46. & six=6.d0)
  47. PARAMETER (cm5s3=-1.66666666666666666666666666666666666)
  48. PARAMETER (cm7s3=-2.33333333333333333333333333333333333)
  49. PARAMETER (cm2s3=-0.66666666666666666666666666666666666)
  50. PARAMETER (cm1s3=-0.33333333333333333333333333333333333)
  51. PARAMETER (cm4s3=-1.33333333333333333333333333333333333)
  52. if(ndi.ne.-2) go to 100
  53. * formulation 2D contraintes planes
  54. C PARAMETRES MATERIAU DU MODELE MOONEY-RIVLIN :
  55. C ===============================================
  56. C= Dans le cas du modele de Mooney-Rivlin,
  57. C= la densite d'energie de deformation est definie par :
  58. C= W = Coe1 * (I1-3.) + Coe2 * (I2-3.)
  59. C= et le module de cisaillement est donne par :
  60. C= mu = 2*(Coe1+Coe2)
  61. C= Coe1 et Coe2 : coefficients du materiau
  62. C= I1, I2 et I3 : trois invariants du tenseur de Cauchy-Green droit
  63. C= Incompressibilite ==> I3 = 1
  64. C= On a alors : dW/dI1 = Coe1 et dW/dI2 = Coe2
  65. C= mu = 2*(Coe1+Coe2) = YOUN/3. (car NU proche de 0.5)
  66. Coe1 = PROPS(3)
  67. Coe2 = PROPS(4)
  68. dWdI1 = Coe1
  69. dWdI2 = Coe2
  70.  
  71. C GRADIENT DE LA TRANSFORMATION (FIN DU PAS) :
  72. C ==============================================
  73. C= Les composantes F31, F32, F13, F23 sont nulles.
  74. C= La composante F33 est non nulle.
  75. F11 = DFGRD1(1,1)
  76. F12 = DFGRD1(1,2)
  77. F21 = DFGRD1(2,1)
  78. F22 = DFGRD1(2,2)
  79. * write(6,*) ' f11 f12 f21 f22',f11, f12, f21, f22
  80. C* F33 = DFGRD1(3,3)
  81.  
  82. C TENSEUR DE CAUCHY-GREEN DROIT (FIN DU PAS) :
  83. C ==============================================
  84. C= CG : tenseur des dilatations de Cauchy-Green droit et CG = Ft*F
  85. C= CG33 : composante dans l'epaisseur de la membrane, calculee par la
  86. C= conservation du volume (incompressibilite) : I3 = det(CG) = 1
  87. CG11 = F11*F11 + F21*F21
  88. CG22 = F12*F12 + F22*F22
  89. CG12 = F11*F12 + F21*F22
  90. C* CG33 = F33*F33
  91. DLT = CG11*CG22 - CG12*CG12
  92. CG33 = 1.D0 / DLT
  93.  
  94. C INVARIANTS DU TENSEUR DE CAUCHY-GREEN DROIT :
  95. C ===============================================
  96. C= I1 : 1er invariant et I1 = trace(CG)
  97. C= I2 : 2e invariant (inutilise ici)
  98. C= I3 : 3e invariant (I3 = det(CG) = 1)
  99. I1 = CG11 + CG22 + CG33
  100. C* I2 = 0.5D0*(I1*I1-CG11*CG11-CG22*CG22-CG33*GC33) - CG12*CG12
  101.  
  102. C INVERSE DU TENSEUR DE CAUCHY-GREEN DROIT :
  103. C ============================================
  104. CInv11 = CG22/DLT
  105. CInv22 = CG11/DLT
  106. CInv12 = -CG12/DLT
  107.  
  108. C CONTRAINTES DE PIOLA-KIRCHHOFF DE SECONDE ESPECE :
  109. C ====================================================
  110. C= S ou PK2 : second tenseur des contraintes de Piola-Kirchhoff
  111. C= S = PK2 = -phyd.Inv(CG) + 2.(dW/dI1 + I1.dW/dI2).Iden - 2.dW/dI2.CG
  112. C= Iden : tenseur identite d'ordre 2
  113. C= phyd : pression hydrostatique
  114. C= La pression hydrostatique est calculee a partir de la condition de
  115. C= contrainte planes soit S33 = 0.
  116. zz = dWdI1 + I1*dWdI2
  117. phyd = 2.D0*(zz-dWdI2*CG33) * CG33
  118. S11 = 2.D0*(zz-dWdI2*CG11) - phyd*CInv11
  119. S22 = 2.D0*(zz-dWdI2*CG22) - phyd*CInv22
  120. S12 = -2.D0*( dWdI2*CG12) - phyd*CInv12
  121.  
  122. C CONTRAINTES DE CAUCHY (FIN DU PAS) :
  123. C ======================================
  124. C= STRESS : Contraintes de Cauchy = (det(F)-1) F S Ft
  125. C= L'hypothese d'incompressibilite se traduit par det(F)=1.
  126. STRESS(1) = (S11*F11*F11) + (2.D0*F11*F12*S12) + (F12*F12*S22)
  127. STRESS(2) = (S11*F21*F21) + (2.D0*F21*F22*S12) + (F22*F22*S22)
  128. C* STRESS(3) = 0.
  129. STRESS(4) = (F21*F11*S11) + (F21*F12*S12) + (F22*F11*S12)
  130. & + (F22*F12*S22)
  131. * write(6,*) 'stres',STRESS(1),STRESS(2),STRESS(4)
  132. RETURN
  133. 100 CONTINUE
  134. * formulation 3D massive
  135. if(ndi.ne.2) go to 200
  136. C PARAMETRES MATERIAU :
  137. C =======================
  138. C= Dans le cas du modele Mooney Rivlin compressible,
  139. C= l'energie de deformation est decomposee en deux termes decouples :
  140. C= - la partie isochorique ou incompressible Wiso, fonction des inva-
  141. C= riants du tenseur de Cauchy-Green modifie ;
  142. C= - la partie purement volumique Wvol, dependant de J=det(F).
  143. C= Pour le present modele, nous avons :
  144. C= - Wiso = Wiso(I1bar,I2bar) = Coe2 * (I1bar-3.) + Coe3 * (I2bar-3.)
  145. C= - Wvol = Wvol(J) = 1/Coe1 * (J-1)*(J-1)
  146. C= Coe2 et Coe3 : coefficients du materiau
  147. C= le module de cisaillement est egal a mu = 2*Coe2
  148. C= bbar : tenseur de Cauchy-Green gauche modifie
  149. C= Par definition, bbar = J**(-2/3)*(F.Ft) = J**(-2/3)*CGg
  150. C= I1bar : 1er invariant de bbar (= trace(bbar))
  151. C= Dans le cas de la quasi-incompressibilite, c.a.d. J proche de 1,
  152. C= Wvol peut etre interpretee comme une fonction de penalisation.
  153. C* Youn = PROPS(1)
  154. C* Pois = PROPS(2)
  155. Coe1 = PROPS(3)
  156. Coe2 = PROPS(4)
  157. Coe3 = PROPS(5)
  158.  
  159. C GRADIENT DE LA TRANSFORMATION (FIN DU PAS) :
  160. C ==============================================
  161. F11 = DFGRD1(1,1)
  162. F12 = DFGRD1(1,2)
  163. F13 = DFGRD1(1,3)
  164. F21 = DFGRD1(2,1)
  165. F22 = DFGRD1(2,2)
  166. F23 = DFGRD1(2,3)
  167. F31 = DFGRD1(3,1)
  168. F32 = DFGRD1(3,2)
  169. F33 = DFGRD1(3,3)
  170.  
  171. C JACOBIEN DE LA TRANSFORMATION (FIN DU PAS) :
  172. C ==============================================
  173. detF = F11*F22*F33 - F12*F21*F33 + F12*F23*F31
  174. & + F13*F32*F21 - F13*F31*F22 - F23*F32*F11
  175.  
  176. C TENSEUR DES DEFORMATIONS DE CAUCHY-GREEN GAUCHE
  177. C =================================================
  178. C= Tenseur de Cauchy-Green gauche CGg = F.Ft
  179. C G11 G22 G33 G12 G13 G23
  180. CGg1 = F11*F11 + F12*F12 + F13*F13
  181. CGg2 = F21*F21 + F22*F22 + F23*F23
  182. CGg3 = F33*F33 + F31*F31 + F32*F32
  183. CGg4 = F11*F21 + F12*F22 + F13*F23
  184. CGg5 = F11*F31 + F12*F32 + F13*F33
  185. CGg6 = F21*F31 + F22*F32 + F23*F33
  186. C= Tenseur de Cauchy-Green gauche CGg2 = (F.Ft)*(F.Ft)
  187. C indices (1 a 6) =(11 22 33 12 13 23
  188. C2Gg1 = CGg1*CGg1 + CGg4*CGg4 + CGg5*CGg5
  189. C2Gg2 = CGg4*CGg4 + CGg2*CGg2 + CGg6*CGg6
  190. C2Gg3 = CGg5*CGg5 + CGg6*CGg6 + CGg3*CGg3
  191. C2Gg4 = CGg1*CGg4 + CGg4*CGg2 + CGg5*CGg6
  192. C2Gg5 = CGg1*CGg5 + CGg4*CGg6 + CGg5*CGg3
  193. C2Gg6 = CGg4*CGg5 + CGg2*CGg6 + CGg6*CGg3
  194. C INVARIANTS : I1B= J**-1/3 I1
  195. I1B = (detF**cm2s3)*(CGg1+CGg2+CGg3)
  196. C CONTRAINTES DE CAUCHY (FIN DU PAS) :
  197. C ======================================
  198. C= Les contraintes de Cauchy STRESS se decomposent en deux termes :
  199. C= - STRESS = SCvol + SCiso
  200. C= - SCvol = phyd.Iden = (dWvol(J)/dJ).Iden
  201. C= avec Iden = tenseur identite d'ordre 2,
  202. C= phyd = pression hydrostatique,
  203. C= avec Gam1 = (dWiso/dI1bar+I1bar.dWiso/dI2bar)
  204. C= Gam2 = (dWiso/dI2bar)
  205. C= - SCiso = 2. J**(-5/3)Gam1 CGg - 2. J**(-7/3) Gam2 C2Gg
  206. C
  207. C ENTRER LE MATERIAU ICI :
  208. dWisodI1bar = Coe2
  209. dWisodI2bar = Coe3
  210.  
  211. Gam1 = (dWisodI1bar+I1B*dWisodI2bar)
  212. Gam2 = (dWisodI2bar)
  213.  
  214. phyd = 2. * (detF-1.0D0) / Coe1
  215.  
  216. eg1 = 2. * (detF**cm5s3) * Gam1
  217. eg2 = 2. * (detF**cm7s3) * Gam2
  218.  
  219.  
  220. STRESS(1) = eg1 * CGg1 - eg2 * C2Gg1 + phyd
  221. STRESS(2) = eg1 * CGg2 - eg2 * C2Gg2 + phyd
  222. STRESS(3) = eg1 * CGg3 - eg2 * C2Gg3 + phyd
  223. STRESS(4) = eg1 * CGg4 - eg2 * C2Gg4
  224. STRESS(5) = eg1 * CGg5 - eg2 * C2Gg5
  225. STRESS(6) = eg1 * CGg6 - eg2 * C2Gg6
  226. return
  227. 200 if(ndi.ne.-1) go to 300
  228. C1 = PROPS(3)
  229. C2 = PROPS(4)
  230. Coe1 = PROPS(5)
  231. C Deformation totale stockage perso dans var 4 5 6
  232. C ******************** GREEN LAGRANGE **************
  233. C
  234. F11 = dfgrd1(1,1)
  235. F21 = dfgrd1(2,1)
  236. F12 = dfgrd1(1,2)
  237. F22 = dfgrd1(2,2)
  238. F33 = dfgrd1(3,3)
  239. detF =F11*F22-F12*F21
  240. C ******************** CC = FT * F
  241. C
  242. cc11 = F11**2 + F21**2
  243. cc22 = F12**2 + F22**2
  244. cc12 = (F11*F12)+(F21*F22)
  245. C ******************** CC-1
  246. DLT = CC11*CC22 - CC12*CC12
  247. cinv11 = 1.0D0/DLT*CC22
  248. cinv22 = 1.0D0/DLT*CC11
  249. cinv12 = -1.0D0/DLT*CC12
  250. cinv33 = 1.0D0/DLT
  251. C
  252. C***********************************************************
  253. C det J =1 ==>c
  254. cc33=1.0D0
  255. C trace C et CB
  256. i1= cc11+cc22+cc33
  257. i2= .50D0*(i1**2-cc11**2-cc22**2-cc33**2-2*cc12**2)
  258. C
  259. I1B = (detF**cm2s3)*I1
  260. I2B = (detF**cm4s3)*I2
  261. C===============================================================
  262. C 2 INTEGRER VOTRE MATERIAU MOONEY : DWDI1B = ?, DWDI2B = ?
  263. C attention nom !!! dwdi1 = DWDI1B et dwdi2 = DWDI2B
  264. C Wiso = f(i1b,i2b)
  265. C
  266. dwdi1 = C1
  267. dwdi2 = C2
  268. dwvdj= 2. * (detF-1.0D0) / Coe1
  269. C===============================================================
  270. Jm2s3=detF**cm2s3
  271. Jm1s3=detF**cm1s3
  272. phyd = Jm1s3*dwvdj
  273. C===============================================================
  274. C p : pression hydrostatique obtenu ici a partir de S33 = 0.
  275. p=2.0D0*((dwdi1+i1*dwdi2)-dwdi2*cc33)*cc33
  276. C******************** CONTRAINTES PK2 ******************
  277. C
  278. C Epaisseur a ajouter dans les formules!!!
  279. C234567
  280. S11=2.0D0*Jm2s3*((dwdi1+i1b*dwdi2)-dwdi2*cc11)+phyd*cinv11
  281. S22=2.0D0*Jm2s3*((dwdi1+i1b*dwdi2)-dwdi2*cc22)+phyd*cinv22
  282. S12= 2.0D0*Jm2s3*(-dwdi2*cc12)+phyd*cinv12
  283. S33 =2.0D0*Jm2s3*((dwdi1+i1b*dwdi2)-dwdi2*cc33)+phyd*cinv33
  284. C
  285. C
  286. C******************** CONTRAINTES Cauchy ******************
  287. C A FAIRE PK2 a transformer par cauchy = j-1 F S FT
  288. C ici incompressible j =1 a modifier par F et FT,
  289. C***************************************************************
  290. mdetF = 1.0D0/detF
  291. C
  292. stress(1)=mdetF*((S11*F11**2)+(2.0D0*F11*F12*S12)+(F12**2*S22))
  293. stress(2)=mdetF*((S11*F21**2)+(2.0D0*F21*F22*S12)+(F22**2*S22))
  294. stress(4)=mdetF*((F21*F11*S11)+(F21*F12*S12)+(F22*F11*S12)
  295. . +(F22*F12*S22))
  296. stress(3)= mdetF*S33
  297. return
  298. 300 kinc=-2
  299. return
  300. END
  301.  
  302.  
  303.  
  304.  
  305.  
  306.  
  307.  

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