%%%%單調加載程序

clear;
%%楊召煥?靜三軸試驗結果
%臨界狀態參數he初始應力狀態
%Mc=1.3;Me=1.3;laimuda=0.162;kaipa=0.0243;?
%miu=0.3;e0=0.974;
%st3=300;st1=st3;jmax=-500;?AAA=1;beta1=1;
%楊騰
Mc=1.3;Me=0.65;laimuda=0.162;kaipa=0.0243;?
miu=0.3;e0=0.974;
st3=50;st1=st3;jmax=500;?beta1=1;


%%%%Li?靜三軸
%Mc=0.772;?Me=0.614;laimuda=0.173;kaipa=0.034;miu=0.30;e0=1.021;st1=500;st3=500;jmax=350;?????
%AAA=0.75;beta1=4;

%%%Weald?clay?靜三軸
%Mc=0.95;Me=0.95;laimuda=0.093;kaipa=0.025;miu=0.3;e0=0.632;st1=207;st3=207;max=119;???????????????
%cy=0;?r=2.;?R=0.50;???beta=0;??L=2.;numb=0;?????uu=2.6;UU=120;

%Mc=1.1392;?Me=0.8256;laimuda=0.0508;kaipa=0.0108;miu=0.15;e0=0.62;
%st1=393;st3=393;max=362.5;?cy=0;?beta=-7.2;?r=2;?
%R=0.50;?HL=4.9;%numb=0;?uu=3.4;UU=60;


v0=1+e0;????????????????????%重塑土三軸等壓固結后比容
X0=beta1*v0/（laimuda-kaipa）;%臨時變量

dj=0.1;
djmax=0:dj:jmax;????????????%將jmax按照dj的增量步，將其劃分成增量的表達形式
dqj=dj+zeros（1jmax/dj+1）;??%靜力階段??有效應力形式?偏應力增量

dpzj=dqj/3;?????????????????%靜力階段??總應力形式?平均應力增量
qj=st1-st3+djmax;
pzj=（st1+djmax+2*st3）/3;????%靜力階段??總應力?平均應力全量

pj（1）=pzj（1）;???????????????%靜力階段??有效應力??平均應力初值
pc（1）=pj（1）;????????????????%初始固結壓力

%在加載面（Loading?surface）等向固結的參量pc
uwj（1）=0;???????%定義循環加載時的累計孔壓值該值等于0
stn1j（1）=0;?????%軸向應變初值
stn3j（1）=0;?????%剪切應變初值

pvstnj（1）=0;????%塑性體積應變
pqstnj（1）=0;????%塑性剪切應變
pqstnabsj（1）=0;?%用于計算塑性剪切應變增量的絕對值

evstnj（1）=0;????%彈性體積應變
eqstnj（1）=0;????%彈性剪切應變
tqstnj（1）=0;????%循環加載開始階段的剪切應變全量初始值

for?j=1:jmax/dj+1
???????
????????M=Mc;???????
????????
????????p1（j）=pj（j）;????
????????q1（j）=qj（j）;
???????????????????????
????????%pFpp（j）=2*p1（j）-pc（j）;???????????%定義邊界面F對邊界面上的p點的偏導數；
????????pFpp（j）=M^2*（2*p1（j）-pc（j））;
????????pFpq（j）=2*q1（j）;??????????????????????????????%定義邊界面F對邊界面上的q點的偏導數；
????????pFppc（j）=-M^2*p1（j）;
???????
????????%pFpvstn（j）=-X0*pc（j）.*p1（j）;?????????????????%定義邊界面F對塑性體積應變的偏導數????????????
????????%kp（j）=-（pFpp（j）.*pFpvstn（j））;????????????????%定義邊界面上的塑性模量
????????
????????kp（j）=-X0*（pFpp（j）.*pFppc（j）.*pc（j））;?????????%定義邊界面上的塑性模量??????????
????????K（j）=（（1+e0）/kaipa）*p1（j）;?