function main

clear all

clc

global n G dl M M1 b1 b2 C4 C6 C7 C8 C9 a

T=273.2+112.13;                                 %反应器内平均温度,K

RMI=1.1;                                        %氯化氢乙炔分子比

u=0.062*3.5;                                    %气体线速度,m/s

C_A=exp(1.1835+1.055e-3*T)*4.184e3;             %乙炔热容,J/kg/K

C_H=0.192*4.184e3;                              %氯化氢热容,J/Kg/K

C_V=exp(2.238+1.256e-3*T)*4.184*1000/62.5;      %氯乙烯热容,J/Kg/K

L_A=exp(9.6084+1.413*log(T))*0.1;               %乙炔导热系数,W/m

L_H=exp(-7.574+0.9769*log(T))*0.1;              %氯化氢导热系数,W/m/K

L_V=0.018;                                      %氯乙烯导热系数,W/m/K

r_A=0.0115*(T/353.2)^0.97*0.001;                %乙炔黏度,Pa.s

r_H=0.0165*(T/353.2)^0.95*0.001;                %氯化氢黏度,pa.s

r_V=8.601e-3*(T/250.1)^0.97*0.001;              %氯乙烯黏度,Pa.s

rou_A=0.9130;                                   %乙炔密度,Kg/m3

rou_H=1.2817;                                   %氯化氢密度,Kg/m3

rou_V=2.1946;                                   %氯乙烯密度,Kg/m3

x_A=0.5;

C_av=(C_A*(1.0-x_A)+C_H*(RMI-x_A)+C_V*x_A)/(RMI+1.0-x_A);           %平均热容,J/Kg/K

L_av=(L_A*(1.0-x_A)+L_H*(RMI-x_A)+L_V*x_A)/(RMI+1.0-x_A);           %平均导热系数,W/m/K

r_av=(r_A*(1.0-x_A)+r_H*(RMI-x_A)+r_V*x_A)/(RMI+1.0-x_A);           %平均黏度,Pa.s

rou_av=(rou_A*(1.0-x_A)+rou_H*(RMI-x_A)+rou_V*x_A)/(RMI+1.0-x_A);   %平均密度,Kg/m3

re=0.044*u*rou_av/r_av;                                             %雷诺数

pr=C_av*r_av/L_av;                                                  %普朗特数

l_er=L_av*(10.3+0.037*re*pr);                                       %床层有效导热系数,W/(mK)

dp=0.004327;                                                        %催化剂颗粒直径,m

hw=(2.58*re^(1/3)*pr^(1/3)+0.094*re^0.8*pr^0.4)*L_av/dp;            %床层与器壁间导热系数,W/(m2*K)

%网格和步长

n=5+1;                  %径向有6个点,包括圆心

m=268;                  %？？？m是什么

dr=0.0044;              %径向微元长度

dl=0.0112;              %轴向微元长度

%参数

Cp=1.206e3;             %反应混合物比热容,J/(kg*K)

dH=-109e3;              %反应热,J/mol

rou_b=610;              %催化剂堆积密度,kg/m^3

G=u*rou_av;             %总质量速率,kg/(m^2.s)

v=u*3600/3/2.1;         %乙炔空速,hr-1

C_A0=20.7;              %进口乙炔浓度,mol/m3

Er=dp*u/10;             %径向混合扩散系数,m2/s

Tw=371.2;               %管壁温度,K

%方程的系数

a1=l_er/(G*Cp);         %(3-1)方程中的a1,b1

b1=-dH*rou_b/G/Cp;

a2=Er/u;                %(3-2)方程中的a1,b1

b2=rou_b/u/C_A0;

M=a1*dl/dr^2;

M1=a2*dl/dr^2;

C1=4*dr*hw/(-l_er);     %(3-25)(3-26)中的系数c1-c9 ？？？c1计算公式

C2=(1+C1*M*(1-1/2/n)/4/(1+M))/(1-C1*M*(1+1/2/n)/4/(1+M));

C3=C1/(1+M)/(1-C1*M*(1+1/2/n)/4/(1+M));

C4=b1*C1*dl/4/(1+M)/(1-C1*M*(1+1/2/n)/4/(1+M));

C5=C1*Tw/(1-C1*M*(1+1/2/n)/4/(1+M));

C6=M*(1-1/2/n+C2*(1+1/2/n))/2/(1+M);

C7=(M*C3*(1+1/2/n))/2/(1+M)+(1-M)/(1+M);

C8=(M*C4*(1+1/2/n)/2/(1+M)+b1*dl/2)/(1+M);

C9=M*C5*(1+1/2/n)/2/(1+M);

%初始变量,反应器入口

a=0.209;                %？？？a=0.209

T(1:n-1)=362.8;         %进口温度

T(n)=371.2;             %管壁温度

x(1:n)=0;               %入口氯乙烯浓度

XO=[T;x];

Tresult(1,:)=XO(1,:);

xresult(1,:)=XO(2,:);

%i：径向，j：轴向

%计算轴向起点处的F(i,j), G(i,j) 和 rc(i,j)

%计算对应于j=1时的F(i,1), G(i,1)

FGrc(T,x);

% 沿轴向移动

options=optimset('tolx',1e-20);

options=optimset(options,'tolfun',1e-20);

for j=1:m-1

j=j+1;

X=fsolve(@TxEquations,XO,options);  %求解非线性方程组，代入初值为X0

T=X(1,:);

x=X(2,:);

Tresult(j,:)=T;

xresult(j,:)=x;

%轴向平均温度，转化率用梯形面积法计算？？？没有中点温度，且壁面温度多了一份，是否认为T(1)=T(6)

xa(j,:)=sum(2*xresult(j,2)+4*xresult(j,3)+6*xresult(j,4)+8*xresult(j,5)+5*xresult(j,6))/25;%径向均温

Ta(j,:)=sum(2*Tresult(j,2)+4*Tresult(j,3)+6*Tresult(j,4)+8*Tresult(j,5)+5*Tresult(j,6))/25;%径向平均转化率

%计算F(i,j)和G(i,j)以及反应速度re(i,j)供下次迭代使用

FGrc(T,x)

end

l=dl.*[0:j-1];

r=dr.*[0:n-1];

l=l';

%结果作图

surf(r,l,Tresult)   %反应管轴径向温度分布

xlabel('r(m)')

ylabel('l(m)')

zlabel('T(K)')

figure

plot(l,Ta)          %径向均温轴向分布(模拟出的均温有很大波动)

xlabel('l(m)')

ylabel('T_a_v(K)')

figure

plot(l,xa)          %径向转化率轴向分布

xlabel('l(m)')

ylabel('j_a_v')

figure

surf(r,l,xresult)   %反应管轴径向转化率分布

xlabel('r(m)')

ylabel('l(m)')

zlabel('x')