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 :
© Cast3M 2003 - All rights reserved.
Disclaimer