化工原理里边对精馏塔编程

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）         ; % 
%%%%%%%%%%%%%%%非饱和孔隙介质%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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         ;  % 水的的初始体积分数
ng0=n0*（1-Sr0）     ;  % 空气的初始体积分数
Ks=35*10^9         ;  % 土颗粒的体积模量
Kl=2.2*10^9        ;  % 水的体积模量
Kg=0.1*10^6        ;  % 空气的体积模量
nl=0.001           ;  % 水的粘滞系数
ng=18*10^（-6）      ;  % 空气的粘滞系数
rs=9*10^9          ;  % 土骨架的弹性常数
ms=4*10^9          ;  % 土骨架的弹性常数
%%%%%%%%%%%%%%%均质弹性体%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rse=12*10^9         ;  % 弹性体的弹性常数 
mse=5*10^9         ;  % 弹性体的弹性常数  
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];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%