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;