Télécharger gdgd.eso

Retour à la liste

Numérotation des lignes :

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

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