C HARTS SOURCE CB215821 25/04/08 21:15:19 12227 C C======================================================================= C= MODELE HYPERELASTIQUE HARTSMITH = C= EN GRANDES TRANSFORMATIONS = C= CONTRAINTES PLANES Modele INCOMPRESSIBLE (direction 3) = C= DEFORMATIONS PLANES Modele Quasi INCOMPRESSIBLE = C= TRIDIMENSIONNEL Modele Quasi INCOMPRESSIBLE = C= Contribution de Laurent Gornet - Ecole Centrale de Nantes (2009) = C======================================================================= C= = C= Exemple d'utilisation d'un modele UMAT en grandes transformations = C= = C= Pour plus d'informations, voir la presentation de L. Gornet lors = C= du Club Cast3m 2006, disponible sur le site Web de Cast3m. = C= = C= Note : Actuellement en grandes deformations dans PASAPAS, le modele = C= ne peut contenir que des modèles type UMAT. On ne peut pas = C= "melanger" les derivees objectives dans Cast3m. = C======================================================================= SUBROUTINE HARTS (STRESS, STATEV, DDSDDE, STRAN, DSTRAN, & TIME, DTIME, TEMP, DTEMP, PREDEF, DPRED, & NDI, NSHR, NTENS, NSTATV, PROPS, NPROPS, & DFGRD0, DFGRD1, KSTEP, KINC,SSE,spd,scd,rpl, & DDSDDT,DRPLDE,DRPLDT,cmname,coords,drot,PNEWDT, & celent,NOEL, NPT, LAYER, KSPT ) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) CHARACTER*16 CMNAME INTEGER NDI, NSHR, NTENS, NSTATV, NPROPS, KSTEP, KINC & NOEL, NPT, LAYER, KSPT REAL*8 STRESS(NTENS), STATEV(NSTATV),SSE, SPD, SCD, RPL, & DDSDDE(NTENS,NTENS), STRAN(NTENS), DSTRAN(NTENS), & DDSDDT(NTENS), DRPLDE(NTENS), DRPLDT, & TIME(2), DTIME, TEMP, DTEMP, PREDEF(*), DPRED(*), & PROPS(NPROPS), DFGRD0(3,3), DFGRD1(3,3), & COORDS(3),DROT(3,3),PNEWDT,CELENT REAL*8 CG11,CG22,CG12,CG33, CInv11,CInv22,CInv12, DLT, & I1,I2, dWdI1,dWdI2, phyd, Coe1,Coe2,Coe3 & F11,F21,F12,F22, S11,S12,S21,S22 & F13,F31,F23,F32,F33 & CGg1,CGg2,CGg3,CGg4,CGg5,CGg6 & C2Gg1,C2Gg2,C2Gg3,C2Gg4,C2Gg5,C2Gg6 & dWisodI1bar,dWisodI2bar,eg1,eg2,Gam1,Gam2 REAL*8 Cc11,Cc22,Cc12,Cc33 REAL*8 I1B,I2B,cinv33,detF,mdetF,Jm2s3,Jm1s3 REAL*8 S33,dwvdj,zz,C1,C2 C=============================================== C 1 DECLARER VOTRE MATERIAU HARTSMITH : REAL*8 G,K1,K2 C REAL*8 zero,one,two,three,four,six Parameter(zero=0.d0,one=1.d0,two=2.d0,three=3.d0, & four=4.d0,six=6.d0) PARAMETER (cm5s3=-1.66666666666666666666666666666666666D0) PARAMETER (cm7s3=-2.33333333333333333333333333333333333D0) PARAMETER (cm2s3=-0.66666666666666666666666666666666666D0) PARAMETER (cm1s3=-0.33333333333333333333333333333333333D0) PARAMETER (cm4s3=-1.33333333333333333333333333333333333D0) c print*,'GORNET', ndi if(ndi.ne.-2) go to 100 * formulation 2D contraintes planes C PARAMETRES MATERIAU DU MODELE HARTSMITH : C =============================================== C= Dans le cas du modele de GORNET DESMORAT, C= la densite d'energie de deformation est definie par : C= W = f(I1,I2) I3= 1 C =============================================== C 3 PARAMETRES MATERIAU DU MODELE HARTSMITH : C =============================================== G = PROPS(3) K1 = PROPS(4) K2 = PROPS(5) C Deformation totale stockage perso dans var 4 5 6 C ******************** GREEN LAGRANGE ************** C F11 = dfgrd1(1,1) F21 = dfgrd1(2,1) F12 = dfgrd1(1,2) F22 = dfgrd1(2,2) C******************* C print*,F11 C print*,F21 C print*,F12 C print*,F22 C C ******************** CC = FT * F C cc11 = F11**2 + F21**2 cc22 = F12**2 + F22**2 cc12 = (F11*F12)+(F21*F22) C C print*,cc11 C print*,cc22 C print*,cc22 C ******************** CC-1 DLT = CC11*CC22 - CC12*CC12 cinv11 = 1.0D0/DLT*CC22 cinv22 = 1.0D0/DLT*CC11 cinv12 = -1.0D0/DLT*CC12 C C*********************************************************** C det J =1 ==>c cc33=1.0D0/(cc11*cc22-cc12**2) C trace C i1= cc11+cc22+cc33 i2= .50D0*(i1**2-cc11**2-cc22**2-cc33**2-2*cc12**2) C=============================================================== C 2 INTEGRER VOTRE MATERIAU HARTSMITH : DWDI1 = ?, DWDI2 = ? C Constantes materielles C W = f(i1,i2) C dwdi1 = G *EXP(K1*(i1-3.0D0)*(i1-3.0D0)) C dwdi2 = G*K2/i2 IF (ABS(I2).ne.0.D0) then dwdi2 = G*K2/i2 ELSE dwdi2=1.0D15 END IF C=============================================================== C p : pression hydrostatique obtenu ici a partir de S33 = 0. p=2.0D0*((dwdi1+i1*dwdi2)-dwdi2*cc33)*cc33 C******************** CONTRAINTES PK2 ****************** C S11= (2.0D0*((dwdi1+i1*dwdi2)-dwdi2*cc11)-p*cinv11) S22= (2.0D0*((dwdi1+i1*dwdi2)-dwdi2*cc22)-p*cinv22) S12= (2.0D0*(-dwdi2*cc12)-p*cinv12) C C C******************** CONTRAINTES Cauchy ****************** C A FAIRE PK2 a transformer par cauchy = j-1 F S FT C ici incompressible j =1 a modifier par F et FT, C*************************************************************** C234567 stress(1)=(S11*F11**2)+(2.0D0*F11*F12*S12)+(F12**2*S22) stress(2)=(S11*F21**2)+(2.0D0*F21*F22*S12)+(F22**2*S22) stress(4)=(F21*F11*S11)+(F21*F12*S12)+(F22*F11*S12) . +(F22*F12*S22) RETURN 100 CONTINUE * formulation 3D massive if(ndi.ne.2) go to 200 C PARAMETRES MATERIAU : C ======================= C= Dans le cas du modele Mooney Rivlin compressible, C= l'energie de deformation est decomposee en deux termes decouples : C= - la partie isochorique ou incompressible Wiso, fonction des inva- C= riants du tenseur de Cauchy-Green modifie ; C= - la partie purement volumique Wvol, dependant de J=det(F). C= Pour le present modele, nous avons : C= - Wiso = Wiso(I1bar,I2bar) = Coe2 * (I1bar-3.) + Coe3 * (I2bar-3.) C= - Wvol = Wvol(J) = 1/Coe1 * (J-1)*(J-1) C= Coe2 et Coe3 : coefficients du materiau C= le module de cisaillement est egal a mu = 2*Coe2 C= bbar : tenseur de Cauchy-Green gauche modifie C= Par definition, bbar = J**(-2/3)*(F.Ft) = J**(-2/3)*CGg C= I1bar : 1er invariant de bbar (= trace(bbar)) C= Dans le cas de la quasi-incompressibilite, c.a.d. J proche de 1, C= Wvol peut etre interpretee comme une fonction de penalisation. C* Youn = PROPS(1) C* Pois = PROPS(2) c Coe1 = PROPS(5) c Coe2 = PROPS(3) c Coe3 = PROPS(4) C =============================================== C 2 PARAMETRES MATERIAU DU MODELE HARTSMITH : C =============================================== G = PROPS(3) K1 = PROPS(4) K2 = PROPS(5) Coe1 = PROPS(6) C GRADIENT DE LA TRANSFORMATION (FIN DU PAS) : C ============================================== F11 = DFGRD1(1,1) F12 = DFGRD1(1,2) F13 = DFGRD1(1,3) F21 = DFGRD1(2,1) F22 = DFGRD1(2,2) F23 = DFGRD1(2,3) F31 = DFGRD1(3,1) F32 = DFGRD1(3,2) F33 = DFGRD1(3,3) C JACOBIEN DE LA TRANSFORMATION (FIN DU PAS) : C ============================================== detF = F11*F22*F33 - F12*F21*F33 + F12*F23*F31 & + F13*F32*F21 - F13*F31*F22 - F23*F32*F11 C TENSEUR DES DEFORMATIONS DE CAUCHY-GREEN GAUCHE C ================================================= C= Tenseur de Cauchy-Green gauche CGg = F.Ft C G11 G22 G33 G12 G13 G23 CGg1 = F11*F11 + F12*F12 + F13*F13 CGg2 = F21*F21 + F22*F22 + F23*F23 CGg3 = F33*F33 + F31*F31 + F32*F32 CGg4 = F11*F21 + F12*F22 + F13*F23 CGg5 = F11*F31 + F12*F32 + F13*F33 CGg6 = F21*F31 + F22*F32 + F23*F33 C= Tenseur de Cauchy-Green gauche CGg2 = (F.Ft)*(F.Ft) C indices (1 a 6) =(11 22 33 12 13 23 C2Gg1 = CGg1*CGg1 + CGg4*CGg4 + CGg5*CGg5 C2Gg2 = CGg4*CGg4 + CGg2*CGg2 + CGg6*CGg6 C2Gg3 = CGg5*CGg5 + CGg6*CGg6 + CGg3*CGg3 C2Gg4 = CGg1*CGg4 + CGg4*CGg2 + CGg5*CGg6 C2Gg5 = CGg1*CGg5 + CGg4*CGg6 + CGg5*CGg3 C2Gg6 = CGg4*CGg5 + CGg2*CGg6 + CGg6*CGg3 C INVARIANTS : I1B= J**-1/3 I1, I2B= J**-4/3 I2 I1=(CGg1+CGg2+CGg3) i2= .5D0*(i1**2-C2Gg1-C2Gg2-C2Gg3-2.D0*(C2Gg4+C2Gg5+C2Gg6)) I1B = (detF**cm2s3)*I1 I2B = (detF**cm4s3)*I2 C I1B = (detF**cm2s3)*I1 I2B = (detF**cm4s3)*I2 C CONTRAINTES DE CAUCHY (FIN DU PAS) : C ====================================== C= Les contraintes de Cauchy STRESS se decomposent en deux termes : C= - STRESS = SCvol + SCiso C= - SCvol = phyd.Iden = (dWvol(J)/dJ).Iden C= avec Iden = tenseur identite d'ordre 2, C= phyd = pression hydrostatique, C= avec Gam1 = (dWiso/dI1bar+I1bar.dWiso/dI2bar) C= Gam2 = (dWiso/dI2bar) C= - SCiso = 2. J**(-5/3)Gam1 CGg - 2. J**(-7/3) Gam2 C2Gg C C ENTRER LE MATERIAU ICI : C=============================================================== C 2 INTEGRER VOTRE MATERIAU HARTSMITH : DWDI1 = ?, DWDI2 = ? C Constantes materielles C W = f(i1,i2) C dWisodI1bar = G *EXP(K1*(i1-3.0D0)*(i1-3.0D0)) C dwdi2 = G*K2/i2 IF (ABS(I2).ne.0.D0) then dWisodI2bar = G*K2/i2 ELSE dWisodI2bar=1.0D15 END IF c dWisodI1bar = Coe2 c dWisodI2bar = Coe3 Gam1 = (dWisodI1bar+I1B*dWisodI2bar) Gam2 = (dWisodI2bar) phyd = 2.0D0 * (detF-1.0D0) / Coe1 eg1 = 2.0D0 * (detF**cm5s3) * Gam1 eg2 = 2.0D0 * (detF**cm7s3) * Gam2 STRESS(1) = eg1 * CGg1 - eg2 * C2Gg1 + phyd STRESS(2) = eg1 * CGg2 - eg2 * C2Gg2 + phyd STRESS(3) = eg1 * CGg3 - eg2 * C2Gg3 + phyd STRESS(4) = eg1 * CGg4 - eg2 * C2Gg4 STRESS(5) = eg1 * CGg5 - eg2 * C2Gg5 STRESS(6) = eg1 * CGg6 - eg2 * C2Gg6 return 200 if(ndi.ne.-1) go to 300 C Deformations Planes C G = PROPS(3) K1 = PROPS(4) K2 = PROPS(5) Coe1 = PROPS(6) C Deformation totale stockage perso dans var 4 5 6 C ******************** GREEN LAGRANGE ************** C F11 = dfgrd1(1,1) F21 = dfgrd1(2,1) F12 = dfgrd1(1,2) F22 = dfgrd1(2,2) F33 = dfgrd1(3,3) detF =F11*F22-F12*F21 C ******************** CC = FT * F C cc11 = F11**2 + F21**2 cc22 = F12**2 + F22**2 cc12 = (F11*F12)+(F21*F22) C ******************** CC-1 DLT = CC11*CC22 - CC12*CC12 cinv11 = 1.0D0/DLT*CC22 cinv22 = 1.0D0/DLT*CC11 cinv12 = -1.0D0/DLT*CC12 cinv33 = 1.0D0/DLT C C*********************************************************** C det J =1 ==>c cc33=1.0D0 C trace C et CB i1= cc11+cc22+cc33 i2= .50D0*(i1**2-cc11**2-cc22**2-cc33**2-2*cc12**2) C I1B = (detF**cm2s3)*I1 I2B = (detF**cm4s3)*I2 C=============================================================== C 2 INTEGRER VOTRE MATERIAU : DWDI1B = ?, DWDI2B = ? C attention nom !!! dwdi1 = DWDI1B et dwdi2 = DWDI2B C Wiso = f(i1b,i2b) C C=============================================================== C 2 INTEGRER VOTRE MATERIAU HARTSMITH : DWDI1 = ?, DWDI2 = ? C Constantes materielles C W = f(i1,i2) C dwdi1 = G *EXP(K1*(i1-3.0D0)*(i1-3.0D0)) C dwdi2 = G*K2/i2 IF (ABS(I2).ne.0.D0) then dwdi2 = G*K2/i2 ELSE dwdi2=1.0D15 END IF C=============================================================== dwvdj= 2.0D0 * (detF-1.0D0) / Coe1 C=============================================================== Jm2s3=detF**cm2s3 Jm1s3=detF**cm1s3 phyd = Jm1s3*dwvdj C C******************** CONTRAINTES PK2 ****************** C C Epaisseur a ajouter dans les formules!!! C234567 S11=2.0D0*Jm2s3*((dwdi1+i1b*dwdi2)-dwdi2*cc11)+phyd*cinv11 S22=2.0D0*Jm2s3*((dwdi1+i1b*dwdi2)-dwdi2*cc22)+phyd*cinv22 S12= 2.0D0*Jm2s3*(-dwdi2*cc12)+phyd*cinv12 S33=2.0D0*Jm2s3*((dwdi1+i1b*dwdi2)-dwdi2*cc33)+phyd*cinv33 C C C******************** CONTRAINTES Cauchy ****************** C A FAIRE PK2 a transformer par cauchy = j-1 F S FT C ici incompressible j =1 a modifier par F et FT, C*************************************************************** mdetF = 1.0D0/detF C stress(1)=mdetF*((S11*F11**2)+(2.0D0*F11*F12*S12)+(F12**2*S22)) stress(2)=mdetF*((S11*F21**2)+(2.0D0*F21*F22*S12)+(F22**2*S22)) stress(4)=mdetF*((F21*F11*S11)+(F21*F12*S12)+(F22*F11*S12) . +(F22*F12*S22)) stress(3)= mdetF*S33 return 300 kinc=-2 return END