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