Visco.zip_BP模型_border_hardening_应力轴_随动强化模型

function Visco % Visco is the main program that calculates the response of an % elastic viscoplastic metal for uniaxial stress global mat sigt dtt Z jp1 stres %------------------------------------ % Specify the material constants mat = zeros(1,8); E = 200.0; D0 = 1.0e+8; n = 1.0; Z0 = 10.0; Z1 = 15.0; m1 = 0.05e+3; Z3 = 5.0; m2 =0.15e+3; mat(1) = E; % Young's modulus (GPa) mat(2) = D0; % Controls maximum plastic strain rate (1/s) mat(3) = n; % Controls rate sensitivity mat(4) = Z0; % Initial value of isotropic hardening ZI (GPa) mat(5) = Z1; % Maximum value of ZI (GPa) mat(6) = m1; % Controls rate of isotropic hardening (1/GPa) mat(7) = Z3; % Maximum value of directional hardening ZD (GPa) mat(8) = m2; % Controls rate of directional hardening (1/GPa) %------------------------------------- % Specify the specify number of loading cycles nload = 1; % Number of loading cycles % Specify cycle 1 nstep(1) = 400; % Number of steps in the cycle e11f(1) = 5.0e-2; % Final value of the total strain for the cycle rate(1) = 1.0e-3; % Magnitude of the total strain rate for the cycle(1/s) % If rate = 0 then it is assumed that relaxation occurs dtrel(1) = 0.0; % Relaxation time for the cycle. If rate is positive then the magnitude of dtrel is ignored % Specify cycle 2 nstep(2) = 400; e11f(2) = -2.5e-2; rate(2) = 1.0e-3; dtrel(2) = 0.0e+0; % Specify cycle 3 nstep(3) = 400; e11f(3) = 2.5e-2; rate(3) = 1.0e-3; dtrel(3) = 0.0; % Specify cycle 4 nstep(4) = 400; e11f(4) = -2.5e-2; rate(4) = 1.0e-3; dtrel(4) = 0.0e+0; %-------------------------------------- % Calculate the total number of steps ntotal=0; for i=1:nload; ntotal=ntotal+nstep(i); end %---------------------------------------------- % Initialize the stres array % stres(:,1)=time (s); % stres(:,2)=e11=total strain % stres(:,3)=ep11=plastic strain; % stres(:,4)=ZI=isotropic hardening (GPa); % stres(:,5)=bet11=directional hardening parameter (GPa); % stres(:,6)=ZD=directional hardening (GPa); % stres(:,7)=sig11=axial stress (GPa); % stres(:,8)=dWp=rate of plastic work (GPa); % stres(:,9)=value of the function f; stres = zeros(ntotal+1,9); stres(1,4)=Z0; %-------------------------------------------------- % Calculate the time increments for each step dt=zeros(nload); rate(1)=sign(e11f(1))*rate(1); dt(1)=e11f(1)/rate(1)/nstep(1); for i=2:nload if rate(i)>0.0 de11=e11f(i)-e11f(i-1); rate(i)=sign(de11)*rate(i); dt(i)=de11/rate(i)/nstep(i); else dt(i)=dtrel(i)/nstep(i); end end %-------------------------------------------------- % Calculate the response istep=0; for i=1:nload; for jj=1:nstep(i); % Calculate time dtt=dt(i); j=istep+jj; jp1=j+1; stres(jp1,1)=stres(j,1)+dtt; % Calculate total strain stres(jp1,2) = stres(j,2)+dtt*rate(i); % Calculate elastic trial value of stress sigt=E*(stres(j+1,2)-stres(j,3)); % Calculate the scale factor lamda Z = stres(j,4) + stres(j,6); % Use old value of hardening Z lamda=1.0; % Initial guess for lamda lamda=fzero('fun',lamda); % Calculate stress stres(jp1,7)=lamda*sigt; % Calculate plastic strain stres(jp1,3)=stres(jp1,2)-stres(jp1,7)/E; % Calculate increment of plastic work dtWp = lamda*(1.0-lamda)*sigt^2/E; % Calculate rate of plastic work stres(jp1,8)=dtWp/dtt; % Calculate isotropic hardening stres(jp1,4)=Z1-(Z1-stres(j,4))*exp(-m1*dtWp); % Calculate directional hardening u11=sign(sigt); stres(jp1,5)=Z3*u11-(Z3*u11-stres(j,5))*exp(-m2*dtWp); stres(jp1,6)=stres(jp1,5)*u11; end istep=istep+nstep(i); end function f = fun(lamda) % The function fun determines the function % who's root gives the scalar value lambda global mat sigt dtt Z jp1 stres % Input material parameters E=mat(1); D0=mat(2); n=mat(3); % Calculate function f sig=abs(sigt); if sig>0.0 factor=2*dtt*D0/sqrt(3.0)*E/sig; f=1.0-lamda-factor*exp(-0.5*(Z/lamda/sig)^(2.0*n)); else f=0.0; end stres(jp1,9)=f;