J2-Hardening plasticity

Below you may find the constitutive driver for the J2-Plasticity model with isotropic hardening, considering associated flow rule, which means that the flow potential is the same as the yield surface and therefore the normal to the yield surface is identical to the normal to the plastic potential surface. 


function [sig,eps,state_new,E_tang]=J2hardening(oldsig,oldeps,Deps,state_old,modpar)
% function [sig,eps,state_new,E_tang]=J2hardening(oldsig,oldeps,Deps,state_old,modpar)
%----------- Description ---------
% J2 elasto plasticity with isotropic hardening constitutive driver
% ------------ INPUT -------------
% oldsig    = stress vector from the last increment [9X1]
% oldeps    = strain vector from the last increment [9X1]
% Deps      = strain increment [9X1]
% state_old = structure array with the specified fields and values
% could contain eps_e, eps_p, kappa, mu, etc.   
% modpar    = contains material constants
% ----------- OUTPUT -------------
% sig       = stress vector from the current increment [9X1]
% eps       = strain vector from the current increment [9X1]
% state_new = updated structured values
% E_tang    = Tangent stiffness matrix
% Written by 
% Shah, 02 December 2011


E=modpar.E;
nu=modpar.nu;
H=modpar.H;
sy=modpar.sy;
oldeps_p=state_old.eps_p;
oldkap=state_old.kap;
K=E/3/(1-2*nu);
G=E/2/(1+nu);
v1=[1 1 1 0 0 0 0 0 0]';
v2=[0 0 0 1 1 1 1 1 1]';
I=diag(v1+0.5*v2);
Idev=I-(1/3)*(v1*v1');
Ee=2*G*Idev+K*(v1*v1');
I2=diag(v1+v2);
I2dev=I2-(1/3)*(v1*v1');
sigT=oldsig+Ee*Deps;
sigDT=I2dev*sigT;
sigET=sqrt(1.5)*norm(sigDT);
phiT=sigET-sy-oldkap;
sig_m1=(sigT(1)+sigT(2)+sigT(3))/3;
eps=oldeps+Deps;
if phiT<=0
    sig=sigT;
    mu=0;   
    kap=oldkap;
    eps_p=oldeps_p;
    E_tang=Ee;
else
    mu=phiT/(3*G+H);
    c=1-mu*3*G/sigET;
    sigD=sigDT*c;
    sig_m=(sigT(1)+sigT(2)+sigT(3))/3;
    sig=sigD+sig_m.*v1;
    nuT=(1.5)*(sigDT/sigET);
    kap=oldkap+H*mu;
    E_tang=Ee-2*G*((2*G/(3*G+H)).*((oldkap+sy)/sigET)*nuT*nuT'+mu*3*G*Idev/sigET);
    eps_p=oldeps_p+mu;
end
state_new.eps_p=eps_p;
state_new.J2=0;
state_new.I1=0;
state_new.kap=kap;
state_new.mu=mu;