%%%%單調(diào)加載程序

clear;
%%楊召煥?靜三軸試驗(yàn)結(jié)果
%臨界狀態(tài)參數(shù)he初始應(yīng)力狀態(tài)
%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;????????????????????%重塑土三軸等壓固結(jié)后比容
X0=beta1*v0/（laimuda-kaipa）;%臨時(shí)變量

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

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

pj（1）=pzj（1）;???????????????%靜力階段??有效應(yīng)力??平均應(yīng)力初值
pc（1）=pj（1）;????????????????%初始固結(jié)壓力

%在加載面（Loading?surface）等向固結(jié)的參量pc
uwj（1）=0;???????%定義循環(huán)加載時(shí)的累計(jì)孔壓值該值等于0
stn1j（1）=0;?????%軸向應(yīng)變初值
stn3j（1）=0;?????%剪切應(yīng)變初值

pvstnj（1）=0;????%塑性體積應(yīng)變
pqstnj（1）=0;????%塑性剪切應(yīng)變
pqstnabsj（1）=0;?%用于計(jì)算塑性剪切應(yīng)變?cè)隽康慕^對(duì)值

evstnj（1）=0;????%彈性體積應(yīng)變
eqstnj（1）=0;????%彈性剪切應(yīng)變
tqstnj（1）=0;????%循環(huán)加載開始階段的剪切應(yīng)變?nèi)砍跏贾?br/>
for?j=1:jmax/dj+1
???????
????????M=Mc;???????
????????
????????p1（j）=pj（j）;????
????????q1（j）=qj（j）;
???????????????????????
????????%pFpp（j）=2*p1（j）-pc（j）;???????????%定義邊界面F對(duì)邊界面上的p點(diǎn)的偏導(dǎo)數(shù)；
????????pFpp（j）=M^2*（2*p1（j）-pc（j））;
????????pFpq（j）=2*q1（j）;??????????????????????????????%定義邊界面F對(duì)邊界面上的q點(diǎn)的偏導(dǎo)數(shù)；
????????pFppc（j）=-M^2*p1（j）;
???????
????????%pFpvstn（j）=-X0*pc（j）.*p1（j）;?????????????????%定義邊界面F對(duì)塑性體積應(yīng)變的偏導(dǎo)數(shù)????????????
????????%kp（j）=-（pFpp（j）.*pFpvstn（j））;????????????????%定義邊界面上的塑性模量
????????
????????kp（j）=-X0*（pFpp（j）.*pFppc（j）.*pc（j））;?????????%定義邊界面上的塑性模量??????????
????????K（j）=（（1+e0）/kaipa）*p1（j）;?