example.m

%This script demonstrates the WTS model with typical experimental values found for glassy PMMA.

%typical input variables for PMMA
T=293;                %[K]      temperature
G=990;                %[MPa]    initial elastic shear modulus at strain->0
PR=0.3;               %         initial Poisson's ratio at strain->0
GnH=10;               %[MPa]    neo-Hookean shear modulus in strain hardening
A0=8.6785 * 10^(-23); %[s]      factor involving the fundamental vibration energy
H0=172590;            %[J/mole] energy barrier for plastic flow at zero stress
a=1.8;                %[nm^3]   parameter determining the strain dependence of the activation volume Vact
b=0.5;                %[nm^3]   parameter determining the strain dependence of the activation volume Vact
tol=0.1;              % absolute error tolerance ODE solver

%estimate the bulk modulus K
K=G*(2/3)*(1 + PR)/(1-2*PR);

%define constant true uniaxial compressive strain-rates
strainrate(1)=-0.001; %[1/s]
strainrate(2)=-0.003; %[1/s]
strainrate(3)=-0.01;  %[1/s]

%define range of true uniaxial compressive strain
strainint=[0 -0.8];

%determine and plot the total true axial stress (including the viscous and neo-Hookean contributions),
%the relative volume deformation J, for the true axial strain and true strain-rates defined above
figure;
for i=1:length(strainrate)
    [T33,truestrain,J,S33,R33]=WTS(K,G,GnH,a,b,A0,H0,strainint,strainrate(i),T,tol);
    
    % plot true axial stress - true axial strain
    subplot(2,1,1)
    plot(-truestrain,-T33,'k');hold on;
    plot(-truestrain,-S33,'r');hold on;
    plot(-truestrain,-R33,'b');hold on;
    xlabel('true axial strain \epsilon_z_z [-]')
    ylabel('total true axial stress \sigma [-MPa]')
    legend('total true axial stress', ...
        'viscous stress contribution','elastic stress contribution')
    
    % plot true axial stress - neo-hookean axial strain
    subplot(2,1,2)
    neohookstrain=exp(2*truestrain)-exp(-truestrain);
    plot(-neohookstrain,-T33,'k');hold on;
    plot(-neohookstrain,-S33,'r');hold on;
    plot(-neohookstrain,-R33,'b');hold on;
    xlabel('\lambda^2 - \lambda^-^1 [-]')
    ylabel('total true stress -\sigma [-MPa]')
    legend('total true axial stress', ...
        'viscous stress contribution','elastic stress contribution')
    
end