Hi Debasis

I recommend you the following:

a) Make a function in Matlab that accepts stress tensor. Calculate the Eigenvectors and Eigenvalues using eig() function.

b) Use the analytical definitions of Tresca, Von Mises, M-C and DP to set up the analytical functions with respect to sig1, sig2 and sig3. Use symbolic sig1, sig2 and sig3 variables.

c) Make a loop over values of sig1, sig2 and sig3 to form different triads of (sig1, sig2, sig3) using each analytical definition of Tresca, VM, MC and DP separately.

d) plot each triad of (sig1, sig2, sig3) using plotting functions depending on your desired method of representation.

Try to use meshgrid(), and  ezmesh() functions to form triads and plot 'em.

Hope it helps.

Hi Debasis

I recommend you the following:

a) Make a function in Matlab that accepts stress tensor. Calculate the Eigenvectors and Eigenvalues using eig() function.

b) Use the analytical definitions of Tresca, Von Mises, M-C and DP to set up the analytical functions with respect to sig1, sig2 and sig3. Use symbolic sig1, sig2 and sig3 variables.

c) Make a loop over values of sig1, sig2 and sig3 to form different triads of (sig1, sig2, sig3) using each analytical definition of Tresca, VM, MC and DP separately.

d) plot each triad of (sig1, sig2, sig3) using plotting functions depending on your desired method of representation.

Try to use meshgrid(), and  ezmesh() functions to form triads and plot 'em.

Hope it helps.

This is what I would do.. probably there are some things to fix... e.g. when traction is on probably axes must be shifted... but try it your self... cheers

subplot(2,2,1),MatsuokaNakai(0,30,135,35)

subplot(2,2,2),MatsuokaNakai(0,30,135,125)

subplot(2,2,3),MatsuokaNakai(0,30,-90,90)

subplot(2,2,4),MatsuokaNakai(0,30,45,45)

%********************************************************

%

%   MatsuokaNakai:    function for the 3D representation of the 

%               Matsuoka Nakai failure cone envelop 

%               MatNak(pt,div,azim,elev);

%               pt -> isotropic traction resistance "isotropic traction"

%               div -> number of divisions

%               azim -> horizontal angle

%               elev -> vertical angle

%

%               by Matteo Ciantia v.01 

%

%               All rights reserved

%

%********************************************************

function [xx,yy,zz] = MatsuokaNakai(pt,div,azim,elev);   %(0,30,135,35)

if nargin < 4, elev = 0;end

if nargin < 3, azim = 0;end

if nargin < 2, div = 40;end

if nargin < 1, pt = 0;end

ps=1000;

km=13;

pz=ps+km*pt;

mu=10;

j=1;

shift=0;

dpz=pz/(div);

dpteta=(360/(div));

for pm=1:dpz:pz

  i=1;

  for teta=0:dpteta:360

    rad=teta*pi/180;

    rad3=3*rad;   

     t=(cos(rad3));

    %%%%%%%%soluzione equazione di 3? grado in eta con teta parametro!!!

     if t-0.0001;

    eta(i,1)=0.0;eta(i,2)=-sqrt((3*(mu-9))/(mu-3));eta(i,3)=sqrt((3*(mu-9))/(mu-3));

    etatrue(i)=eta(i,3);

    else

    a2=-(9/2)*((mu-3)/(mu*t));

    a1=0;

    a0=(27/2)*((mu-9)/(mu*t));

        qpol=((1/3)*a1)-(1/9)*(a2^2);

    rpol=((1/6)*(a1*a2-3*a0))-((1/27)*(a2^3));

    d(i)=(qpol^3)+(rpol^2);

    pol = [1  a2  a1  a0];

    sol = roots(pol);

    eta(i,1)=sol(1);eta(i,2)=sol(2);eta(i,3)=sol(3);   

        if  t>0.00001

         etatrue(i)=eta(i,2);

        elseif t

Thank u so much Dr. Matteo  and Dr. Hooshiar...

Dear Debasis

In the use of 3D plot options in MATLAB, you can draw your diagrams. A 3D plot programming of my research was appended in the bellow.  It is hoped that it will be useful.

best regard

--------------------------------------------------------------------

% 3D diagram

global X CS NNN L AT XI

global U UI EQ1 EQ NSTEP

%

% 3d diagram for time (1) or freqency (2)

coding = 1;

if coding == 1

tdr = U;

tdrint = UI;

for i = 1:NSTEP

%break

tdr(2,i) = (tdr(4,i)+tdrint(L/2,i))/2;

tdr(1,i) = (tdr(2,i)+tdrint(L/2,i))/2;

tdr(3,i) = (tdr(2,i)+tdr(4,i))/2;

tdr(NNN,i) = tdr(1,i);

tdr(NNN-1,i) = tdr(2,i);

tdr(NNN-2,i) = tdr(3,i);

%tdr(round(NNN/2),i) = tdr(round(NNN/2)-1,i)+(tdr(round(NNN/2)-1,i)-tdr(round(NNN/2)-2,i))/4;

end

tdr_tot = [tdrint(1:L/2,:);tdr(:,:);tdrint(L/2+1:end,:)];

eq_tot = [EQ(:,1:L/2)';EQ1(:,:)';EQ(:,L/2+1:end)'];

xnew = [XI(1:L/2),X,XI(L/2+1:end)];

time = AT:AT:NSTEP*AT;

figure(1)

hold on

for i=1:1:length(xnew)

%break

plot3(time,xnew(i)*ones(1,length(time))/200,tdr_tot(i,:)/0.001,'k')

end

for i=61:1:154

break

plot3(time,xnew(i)*ones(1,length(time))/200,eq_tot(i,:)/0.001,'k-')

end

view(0,70)

%[tt,xxx] = meshgrid(time,xnew/X(1));

%mesh(tt,xxx,tdr_tot./0.001)

return

end

nfft = 2^nextpow2(3000);

fftsignal = zeros(NNN,nfft);

fftinter = zeros(L,nfft);

ffteq1 = zeros(NNN,nfft);

ffteq = zeros(L,nfft);

for i = 1 : NNN

fftsignal(i,:) = fft(U(i,:),nfft);

ffteq1(i,:) = fft(EQ1(:,i),nfft);

end

for i = 1 : L

fftinter(i,:) = fft(UI(i,:),nfft);

ffteq(i,:) = fft(EQ(:,i),nfft);

end

% correction of sem-sine shaped valley at corner points

for i = 1:nfft

%break

fftsignal(2,i) = (fftsignal(4,i)+fftinter(L/2,i))/2;

fftsignal(1,i) = (fftsignal(2,i)+fftinter(L/2,i))/2;

fftsignal(3,i) = (fftsignal(2,i)+fftsignal(4,i))/2;

fftsignal(NNN,i) = fftsignal(1,i);

fftsignal(NNN-1,i) = fftsignal(2,i);

fftsignal(NNN-2,i) = fftsignal(3,i);

fftsignal(round(NNN/2),i) = fftsignal(round(NNN/2)-1,i)+(fftsignal(round(NNN/2)-1,i)-fftsignal(round(NNN/2)-2,i))/4;

end

%Next, calculate the frequency axis, which is defined by the sampling rate

T = AT;

fs = 1/T;

f = fs/2*linspace(0,1,nfft/2+1);

% dimensionless frequency

df = 2*X(1)*f./CS;

% select dimensionless frequency " max "

sdf_max = 4.0;

sdf_min = 0.125;

vec_max = find(df>=sdf_max);

vec_min = find(df