C234567890123456789012345678901234567890123456789012345678901234567890 SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, 1 RPL,DDSDDT,DRPLDE,DRPLDT, 2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME, 3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT, 4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC) INCLUDE 'ABA_PARAM.INC' CHARACTER*80 cmname DIMENSION STRESS(NTENS),STATEV(NSTATV), 1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS), 2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1), 3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3) REAL::SSE,SPD,SCD,RPL,DRPLDE,DRPLDT,TEMP,DTEMP,CELENT,x(6,6) c Integer:: KurvenWoelbung,KurvenWoelbung_alt,Lastaenderung,Kinc, c 1 NDI,NSHR,NTENS,NSTATV,PNEWDT,NPROPS,NOEL,NPT,LAYER,KSPT,KSTEP c 2 ,K1 Real :: Nue_m,Nue_f,Nue12_UD,Nue21_UD,Nue23_UD,E1_UD 1 ,E2_UD,G12_UD,G23_UD,eps_sp(6),eps_v(6),eps_v_ges(6),eps_p(6) 2 ,x1(6,6),x2(6,6),S_ges(6,6),S_f(6,6),S_m(6,6),Zeit,Zeit_neu, 3 DDSDDE_alt(6,6),DDSDDE_0(6,6),DDSDDE_inv(6,6),DDSDDE_neu(6,6), 4 DSTRAN_alt(6),deps_sp(6),deps_v(6),deps_p(6),sigma_sp(6), 5 sigma_v(6),sigma_p(6),eps_v_0(6),Deps_v_ges(6),DDSDDE_c(6,6), 6 wahre_Spannung(6),ZERO,ONE,A11,A22,A33,A44,A55,A66,eigenwert1, 7 eigenwert2,eigenwert3,S1111,S2222,S2233,S2211,S1122,S2323,S1212, 8 Faservolumenanteil,fak1(6,6),fak2(6,6),Delta,y1(6),y2(6), 9 E0,E_f,G_f,E,E_c,G_m,E_mm,E_m ,stress_eq PARAMETER(ZERO = 0.D0, ONE = 1.D0) !fitfunktion spezifische Deklarationen.. c--------------- Hauptprogrammteil ---------------- A11=Props(1) A22=Props(2) A33=Props(3) A44=Props(4) A55=Props(5) A66=Props(6) eigenwert1=Props(7) eigenwert2=Props(8) eigenwert3=Props(9) c Komponenten des Eshelby Tensor - vorab berechnet S1111=0.025062 S2222=0.706938 S2233=0.125446 S2211=0.322594 S1122=0.000009 S2323=0.290321 S1212=0.245530 c Materialeigenschaften E0=1100.000000 Nue_m=0.400000 E_f=75000.000000 Nue_f=0.220000 G_f=30737.704918 FaserVolumenAnteil=0.180000 Delta=1.831035 open(unit=11,File='C:\uni\Abaqus\Umatttest\ausg.out', &status='unknown') Zeit=Dtime*Kinc stress_eq=sqrt(stress(1)**2+stress(2)**2+stress(3)**2- & stress(1)*stress(2)-stress(2)*stress(3)-stress(1)*stress(3)+ & 3*(stress(4)**2+stress(5)**2+stress(6)**2)) call fitfunktion(Zeit,Stress_eq,E_m) G_m=E_m/(2.0+2.0*Nue_m) c Aufspaltung der Dehnungsanteile: write(11,*)'E_m=',E_m,G_m call GetUDProps(G_f,E_f,E_m,Nue_f,Nue_m, & FaserVolumenAnteil,S1111,S2222,S2211,S1122,S2233,S2323,S1212, & G_m,E1_UD,E2_UD,G12_UD,G23_UD,Nue12_UD,Nue23_UD, & Nue21_UD) if (noel.EQ.1) THEN write(11,*)'2' end if call GetDDSDDE (eigenwert1,eigenwert2,eigenwert3,A11,A22,A33, 1 A44,A55,A66,E1_UD,E2_UD,G12_UD,Nue12_UD,Nue23_UD, 2 Nue21_UD,Delta,x) DDSDDE=x dstress=Zero DO 70 K1=1,NTENS DO 60 K2=1,NTENS STRESS(K1)=STRESS(K1)+DDSDDE(K2,K1)*DSTRAN(K2) 60 CONTINUE 70 CONTINUE c dstress=matmul(DDSDDE,DSTRAN) c STRESS=STRESS+dstress END !(Hauptprogramm) c ****************************************************************** c ************ SUBROUTINEN DIE VERWENDET WERDEN ******************** c ****************************************************************** SUBROUTINE GetUdProps(G_f,E_f,E_m,Nue_f,Nue_m, 1 FaserVolumenAnteil,S1111,S2222,S2211,S1122,S2233,S2323,S1212, 2 G_m,E1_UD,E2_UD,G12_UD,G23_UD,Nue12_UD,Nue23_UD, 3 Nue21_UD) !INCLUDE 'ABA_PARAM.INC' c Bestimmung alle Verbundkennwerte fuer den Fall, dass alle Fasern c unidirektional orientiert sind, nach Tandon Weng c (aufbauend auf 1. Mean Field Therorie(Mori-Tanaka) 2. Einschluss- c problem nach Eshelby) C C CHARACTER*8 cmname Real::Nue_f,Nue_m,Nue12_UD,Nue23_UD,Nue21_UD,lambda_f,lambda_m !Tandon-Weng UD-Eigenschaften, S=Eshelby Tensor lambda_f=Nue_f*E_f/((1+Nue_f)*(1-2*Nue_f)) lambda_m=Nue_m*E_m/((1+Nue_m)*(1-2*Nue_m)) D1=1+(2*G_f-2*G_m)/(lambda_f-lambda_m) D2=(lambda_m+2*G_m)/(lambda_f-lambda_m) D3=lambda_m/(lambda_f-lambda_m) B1=FaserVolumenAnteil*D1+D2+(1-FaserVolumenAnteil)* & (D1*S1111+2*S2211) B2=FaserVolumenAnteil+D3+(1-FaserVolumenAnteil)* & (D1*S1122+S2222+S2233) B3=FaserVolumenAnteil+D3+(1-FaserVolumenAnteil)* & (S1111+(1+D1)*S2211) B4=FaserVolumenAnteil*D1+D2+(1-FaserVolumenAnteil)* & (S1122+D1*S2222+S2233) B5=FaserVolumenAnteil+D3+(1-FaserVolumenAnteil)* & (S1122+S2222+D1*S2233) A=2*B2*B3-B1*(B4+B5) A1=D1*(B4+B5)-2*B2 A2=(1+D1)*B2-(B4+B5) A3=B1-D1*B3 A4=(1+D1)*B1-2*B3 A5=(1-D1)/(B4-B5) E1_UD=E_m/(1+FaserVolumenAnteil*(A1+2*Nue_m*A2)/A) E2_UD=E_m/(1+FaserVolumenAnteil*(-2*Nue_m*A3+(1-Nue_m)*A4 &+(1+Nue_m)*A*A5)/(2*A)) G12_UD=G_m*(1+FaserVolumenAnteil/((G_m/(G_f-G_m))+ &(2*S1212-2*S1212*FaserVolumenAnteil))) G23_UD=G_m*(1+FaserVolumenAnteil/((G_m/(G_f-G_m))+ &(2*S2323-2*S2323*FaserVolumenAnteil))) Nue12_UD=(Nue_m*A-FaserVolumenAnteil*(A3-Nue_m*A4))/ &(A+FaserVolumenAnteil*(A1+2*Nue_m*A2)) Nue21_UD=Nue12_UD*E2_UD/E1_UD Nue23_UD=0.5*E2_UD/G23_UD-1 write(11,*)'Nues',E1_UD,G12_UD,Nue12_UD RETURN END SUBROUTINE GetDDSDDE(eigenwert1,eigenwert2,eigenwert3,A11,A22,A33, 1 A44,A55,A66,E1_UD,E2_UD,G12_UD,Nue12_UD,Nue23_UD, 2 Nue21_UD,Delta,DDSDDE) !INCLUDE 'ABA_PARAM.INC' c Bestimmung des Elastizitaetstensors, Orientierungsmittelung nach c Tucker (orthotropic fitted closure -ORF) um Faserorientierungs- c Verteilung zu beruecksichtigen C Dimension:: DDSDDE(6,6) real ::Nue12_UD,Nue23_UD,Nue21_UD,ONE real::C11_UD,C22_UD,C33_UD,C44_UD,C12_UD,C13_UD,C23_UD, 1 C55_UD,C66_UD,Delta PARAMETER(ZERO = 0.D0, ONE = 1.D0) !write(11,*)'1' !Delta=1.0*((1+Nue23_UD)*(1-Nue23_UD-2*Nue21_UD*Nue12_UD)) DDSDDE=ZERO DO 20 K1=1,NTENS DO 10 K2=1,NTENS DDSDDE(K2,K1)=0.D0 10 CONTINUE 20 CONTINUE C C11_UD=(1-Nue23_UD*Nue23_UD)*E1_UD*Delta C22_UD=(1-Nue12_UD*Nue21_UD)*E2_UD*Delta C33_UD=C22_UD C12_UD=(Nue21_UD+Nue23_UD)*E2_UD*Delta C13_UD=C12_UD !write(11,*)'2' C23_UD=(Nue23_UD+Nue12_UD*Nue21_UD)*E2_UD*Delta C44_UD=G12_UD C55_UD=G12_UD !C66_UD=(1-Nue23_UD-2*Nue21_UD*Nue12_UD)*E2_UD*(2*Delta) !write(11,*)'Cs:',C11_UD,C22_UD,C12_UD,C23_UD,C44_UD B1=C11_UD+C22_UD-2*C12_UD-4*C44_UD B2=C12_UD-C23_UD B3=C44_UD+0.5*(C23_UD-C22_UD) B4=C23_UD B5=0.5*(C22_UD-C23_UD) !write(11,*)'4',B1,B2,B3,B4,B5 DDSDDE(1,1)=B1*A11+2*B2*eigenwert1+4*B3+B4+2*B5 DDSDDE(2,2)=B1*A22+2*B2*eigenwert2+4*B3+B4+2*B5 DDSDDE(3,3)=B1*A33+2*B2*eigenwert3+4*B3+B4+2*B5 DDSDDE(1,2)=B1*A44+B2*(eigenwert1+eigenwert2)+B4 DDSDDE(1,3)=B1*A55+B2*(eigenwert1+eigenwert3)+B4 DDSDDE(2,3)=B1*A66+B2*(eigenwert2+eigenwert3)+B4 DDSDDE(2,1)=DDSDDE(1,2) DDSDDE(3,1)=DDSDDE(1,3) DDSDDE(3,2)=DDSDDE(2,3) DDSDDE(4,4)=B1*A44+B3*(eigenwert1+eigenwert2)+B5 DDSDDE(5,5)=B1*A55+B3*(eigenwert3+eigenwert1)+B5 DDSDDE(6,6)=B1*A66+B3*(eigenwert2+eigenwert3)+B5 write(11,*)'5',DDSDDE RETURN END Subroutine fitfunktion(Time,Stress,E_m) real :: E_m,stress,time,e e=2.718281828 Time=abs(Time/100000.000000) Stress=abs(Stress/16.000000) write(11,*)'FF3',time,stress E_m=1100.000000/(1081299.36+ &1457611.42+ &34659.272)*(1081299.36*e**(-0.113138151*time*stress**1.56572474)+ &1457611.42*e**(-13.8167772*time**0.217218944*stress)+ &34659.272*e**(-90.8730882*(time*stress)**0.00882264013)) &-0.0990403305*(1-e**(-3.01439369*time**0.5*stress**0.5)) c E_m=1100.0 write(11,*)'FF4',E_m Return End