Предложения по модификации кода MATlab - PullRequest
Новая
       47

Предложения по модификации кода MATlab

0 голосов
/ 28 апреля 2020

В настоящее время этот код дает значение t для n=1. Я хочу t значение для каждого n (скажем, 400). Итак, мне нужно 400 t значений. Какую модификацию я должен сделать? 

clear all
close all
clc


Ep_in=0;                                     %Volumetric strain

%Storage parameters
%Input values
a=0.02;
mu=0.3;                              
K=2000;                                     %Bulk modulus(in kPa)                                    
n_p=1;                                      %Peak stress ratio
M=1;                                      %Critical stress ratio
G=(3*K*(1-2*mu))/(2*(1+mu));
tol=0.01;

%Base state 
peff_in=100;                                 %Mean effective stress(in kPa)
q_in=0;                                      %Deviatoric stress 
ny_in=q_in/peff_in;                          %Initial yield stress ratio
Eq_in=0;                                     %Shear strain
v=400;
p=zeros(1,v);
q=zeros(1,v);

%Loop
n=1;
peff(1)=peff_in;
q(1)=q_in;
ny(1)=ny_in;
Eq(1)=Eq_in;                                %Deviatoric strain
Ep(1)=Ep_in;                                %Volumetric strain
m=1;
x=100;
r(m)=x;
dEq=0.005;

while (n<=v)
  dEp=x*dEq;
  dp_eff_t=(K*dEp)-(((3*G*K*(M-ny(n)).*dEq)+(-1*K*K*ny(n).*(M-ny(n))*dEp))/((3*G-(K*ny(n).*(M- 
  ny(n))))+(peff(n).*((n_p-ny(n)).^2)./(a*n_p))));
  dq_t=((3*G*dEq)-(((-3*G*K*ny(n).*dEp)+(9*G*G*dEq))/((3*G-(K*ny(n).*(M-ny(n))))+(peff(n).*((n_p- 
  ny(n)).^2)./(a*n_p)))));
  R(m)=dq_t/dp_eff_t
  if R(m)>3
    z=m
    x1=r(m-1)
    x2=r(m)
    y1=R(m-1)
    y2=R(m)
    t=x2;
    for t=x2:0.001:x1
      dEp=t*dEq;
      idp_eff_t=(K*dEp)-(((3*G*K*(M-ny(n)).*dEq)+(-1*K*K*ny(n).*(M-ny(n))*dEp))/((3*G-(K*ny(n).*(M- 
              ny(n))))+(peff(n).*((n_p-ny(n)).^2)./(a*n_p))));
      idq_t=((3*G*dEq)-(((-3*G*K*ny(n).*dEp)+(9*G*G*dEq))/((3*G-(K*ny(n).*(M-ny(n))))+(peff(n).*((n_p- 
          ny(n)).^2)./(a*n_p)))));
      iR=idq_t/idp_eff_t
      if abs(iR-3)<=tol
        sid=t
        while n<=v
          dEp=t*dEq;
          idp_eff_t=(K*dEp)-(((3*G*K*(M-ny(n)).*dEq)+(-1*K*K*ny(n).*(M-ny(n))*dEp))/((3*G-(K*ny(n).*(M-ny(n))))+(peff(n).*((n_p-ny(n)).^2)./(a*n_p))));
          idq_t=((3*G*dEq)-(((-3*G*K*ny(n).*dEp)+(9*G*G*dEq))/((3*G-(K*ny(n).*(M-ny(n))))+(peff(n).*((n_p-ny(n)).^2)./(a*n_p)))));
          dq(n)=idq_t;
          dp_eff(n)=idp_eff_t;
          S(n)= dq(n)/dp_eff(n);
          q(n+1)=q(n)+dq(n);
          peff(n+1)=peff(n)+dp_eff(n);
          ny(n+1)=(q(n+1))./(peff(n+1));
          Eq(n+1)=Eq(n)+dEq; 
          Ep(n+1)=Ep(n)+dEp; 
          n=n+1; 
        end
        break
      end
    end  
  end 

  m=m+1
  x=x/2; 

  r(m)=x
  n


  pause

end
R
n
z
dq;
dp_eff;
S
q;
ny;
figure(1)
plot(Eq,ny/n_p,'s--g')
xlabel('Deviatoric strain')
ylabel('Stress ratio')

figure(2)
plot(Eq,q,'s--r')
xlabel('deviatoric strain')
ylabel('deviatoric stress')
...