clear?all?
syms?o?oo??a11?a12?a13?a21?a22?a23?a31?a32?a33?b11?b12?b13?...
????b21?b22?b23?b31?b32?b33?Nss?Nsl?Nsg?Nls?Nll?Nlg?Ngs?Ngl?Ngg?es0?el0??...
????eg0?nl?ese?ng?sil?sig?ns0?nl0?ng0?ks?kp2?rs?ms?rse?mse?w?kp?krse?kse?kipe?...
????krpe?Ks?Kl?Kg?thl?thg?A1?A2?A3?A4?A5?B1?B2?B3?B4?B5?Sr0?n0?k?krl?krg?m?w?kr?m11...
????m12?m13?m14?m15?m16?m21?m22?m23?m24?m25?m26?m31?m32?m33?m34?m35?m36?m41?m42?m43...
????m44?m45?m46?m51?m52?m53?m54?m55?m56?m61?m62?m63?m64?m65?m66?arp?brs?ap1?ap2?ap3?bts?q1?q2?q3?q4?q5?q6?...
????ktp?kts?kr

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w=1000*2*3.141593????;?%??
Sr0=0.85?????????????;?%
H=100????????????????;?%
k=3*10^（-13）?????????;?%?
%%%%%%%%%%%%%%%非飽和孔隙介質(zhì)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m=0.5?;???????????????%
n=1/（1-m）;
avg=0.0001;
es0=2700???????????;??%?
el0=1000???????????;??%?
eg0=1.2????????????;??%?
n0=0.3?????????????;??%?
ns0=1-n0???????????;??%?
nl0=n0*Sr0?????????;??%?水的的初始體積分數(shù)
ng0=n0*（1-Sr0）?????;??%?空氣的初始體積分數(shù)
Ks=35*10^9?????????;??%?土顆粒的體積模量
Kl=2.2*10^9????????;??%?水的體積模量
Kg=0.1*10^6????????;??%?空氣的體積模量
nl=0.001???????????;??%?水的粘滯系數(shù)
ng=18*10^（-6）??????;??%?空氣的粘滯系數(shù)
rs=9*10^9??????????;??%?土骨架的彈性常數(shù)
ms=4*10^9??????????;??%?土骨架的彈性常數(shù)
%%%%%%%%%%%%%%%均質(zhì)彈性體%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rse=12*10^9?????????;??%?彈性體的彈性常數(shù)?
mse=5*10^9?????????;??%?彈性體的彈性常數(shù)??
ese=2460???????????;??%?彈性體的密度
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
krl=（Sr0）^（1/2）*（1-（1-Sr0^（1/m））^m）^2;
krg=（1-Sr0）^（1/2）*（1-Sr0^（1/m））^（2*m）;
sil=nl0^2*nl/（k*krl）;
sig=ng0^2*ng/（k*krg）??;
thl=（（1-m）/（avg*m））*（Sr0^（-1/m）-1）^（-m）*Sr0^（-1/m-1）*Sr0*（1-Sr0）;
thg=（（1-m）/（avg*m））*（Sr0^（-1/m）-1）^（-m）*Sr0^（-1/m-1）*Sr0*（1-Sr0）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A1=-thl-Kl;
A2=thg+Kg;
A4=-thl;
A5=thg;
B1=Sr0*（Kl+thl）+Ks*nl0/ns0;
B2=（1-Sr0）*（Kg+thg）+Ks*ng0/ns0;
B3=Ks;
B4=Sr0*thl+Ks*nl0/ns0;
B5=（1-Sr0）*thg+Ks*ng0/ns0;
Nss=Ks*（ns0+（-B3*A2*nl0）/（A2*B1-A1*B2）+（A1*B3*ng0）/（A2*B1-A1*B2））;%?Nss廇惀Rss椆???
Nsl=（A2*B3*nl0*Kl）/（A2*B1-A1*B2）;
Nls=（A2*B3*nl0*Kl）/（A2*B1-A1*B2）;
Nsg=（-A1*B3*ng0*Kg）/（A2*B1-A1*B2）;
Ngs=（-A1*B3*ng0*Kg）/（A2*B1-A1*B2）;
Nll=nl0*Kl*（A2*B4-A4*B2）/（A2*B1-A1*B2）;
Nlg=ng0*Kg*（A4*B1-A1*B4）/（A2*B1-A1*B2）;
Ngl=ng0*Kg*（A4*B1-A1*B4）/（A2*B1-A1*B2）;
Ngg=ng0*Kg*（A5*B1-A1*B5）/（A2*B1-A1*B2）;
M=[ns0*es000;0nl0*el00;00ng0*eg0];
C=[sil+sig-sil-sig;-silsil0;-sig0sig];
Kp=[Nss+ns0*rs+2*ns0*msNslNsg;NlsNllNlg;NgsNglNgg]???;????????
Ks=[ns0*ms00;000;000];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%