新人求救！！！急啊！！

4everlove

这是一个精馏塔仿真程序，我怎么看不懂啊！我因为做毕设刚刚接触mathematica，以前制学过c语言，看这个有点吃力，那位大侠替我解释下，或者给我发个mathematica程序设计的简单知识，比如代码解释什么的！我急啊！！




*
This is for the operating condition:
The TOTAL NUMBER OF STGRES EQUALS: 20+reboiler = 7/6/7 +reboiler.
This is a static model for an ideal reactive distillation column.
The reversible reaction involved in the distillation
column is: A + B = C + D.
The production spec. is 0.95 for both overhead and bottom.
The total holdup on reactive stages is 6 kmol.

RDSIM---reactive distillation simulation;
propdata---property data base;
xtov---from liquid composition to vapor flow calculation;
NewRap---Newton-Rapsion iteration subroutine

Mbalance---material balance subroutine;
Reaction---reaction subroutine;
acf---activity energy of forward reaction, cal/mol;
acb---activity energy of backward reaction, cal/mol;
af---preexponential factor of backward reaction;
alpha---relative volatility;
Avp--- vapor pressure coefficient;
bf---preexponential factor of forward reaction;
Bvp---vapor pressure coefficient;
Bc---specified D composition in bottom product;
component 1---A;
component 2---B;
component 3---C;
component 4---D;
cov---Newton-Raphsion covergent criterion;
cov1, cov2, cov3---error of vapor-liquid equilibrium;
dx---step length of Newton-Raphsion iteration;
error---error of component material balance equations;
ext---adjustion of step length;
fcomp---feed composition;
feedflow---feed flow rate;
holdup---stage holdup;
hor---heat of reaction, cal/mol;
hov---heat of vaporization,cal/mol;
i---component index;
j---stage index;
Lc---subsitute of liquid;
liquid---liquid flow rate;
mstep---Newton-Raphsion iteration number;
nstep---number of material balance calculation;
nt---total number of stages;
nc---total number of components;
ns---number of stages in the stripping section;
nfeed1---number of first feeding stage;
nfeed2---number of second feeding stage;
nrb---number of the bottom reactive stage;
press--- pressure of reactive distillation column;
ptimes---frequency of printing calculation results;
qfeed----feed thermal condition;
Qcon---condenser duty;
Qrea---exothermic reaction heat;
Qreb---reboiler duty;
rmax---maximum error of material balance;
rate---reaction rate, kmol/s;
R1---ideal gas law constant;
RR---reflux rate;
T0---ambient temperature;
tc--- subsitute of temp;
Tc---specified C composition in top product ;
temp---stage temperature;
vapor--- vapor flow rate;
vc---- subsitute of vapor;
xc---subsitute of xcomp;
xcomp--- liquid phase composition;
xtry, xtry1--- working arrays for decison variables;
yc---subsitute of ycomp;
ycomp---vapor phase composition;
*)

Clear[RDSIM, NewRap, Mbalance, Reaction, Tc, Bc,
nfeed1, nfeed2, ns, nt, xtry,
feedflow, qfeed, fcomp, press, holdup,
xcomp, ycomp, liquid, vapor, temp, RR,
xc,yc,Lc,vc,
ratio, Avp, Bvp, error, rmax,
alpha, T0, mstep, nstep];

RDSIM[nfeed1_, nfeed2_, ns_, Tc_, Bc_, COPROD1_, COPROD2_]:=
Module[{NewRap, Mbalance, xtry,
feedflow, qfeed, fcomp,press,holdup,
xcomp, ycomp, liquid, vapor, temp, Qheat, RR,
xc,yc,Lc,vc,ratio,
Avp, Bvp, error, rmax, alpha, T0, mstep, nstep},

propdata[]:=Module[{},
(* Ideal gas law Constant, cal/mol/K *)
R1=1.987;
(* Vapor Pressure Constants: Ln P = Avp - Bvp/T *)
Avp[1]=12.3463; Avp[2]=11.6531;
Avp[3]=13.0394; Avp[4]=10.9600;
Bvp[1]=3862;Bvp[2]=3862;
Bvp[3]=3862;Bvp[4]=3862;
(* Relative Volatilities *)
alpha[1]=4;alpha[2]=2;
alpha[3]=8;alpha[4]=1;
(*  Activity energy, cal/mol *)
acf=30000;
acb=140047/4.184;
(*  Heat of reaction, cal/mol *)
hor=14527/4.184;
(*  Heat of vaporization, mol/cal *)
hov=6944;
(*  Preexponential factors, kf=af*Exp[-acf/(R*T)] and
kb=ab*Exp[-acb/(R*T)] *)
af=0.008*Exp[acf/(R1*366)];
ab=0.004*Exp[acb/(R1*366)];
(* Zero Absolute Temperature *)
T0=273.15;
];

xtov[]:=Module[{cov1, cov2, cov3, tguess1,
tguess2, tguess3, ysum, Vp},
Do[
{xc[4,j]=1-xc[1,j]-xc[2,j]-xc[3,j];
If[xc[4,j]<=0,
{xc[4,j]=0;
Do[xc[i,j]=xc[i,j]/(xc[1,j]+xc[2,j]+xc[3,j]+xc[4,j]),{i,1,nc}]}
]
},{j,1,nt}];
Do[
Do[
xtry[(i-1)*nt+k]=xc[i,k]
,{i,1,nc-1}];
,{k,1,nt}];
(* Stage temperature calculation *)
(* Guessing Temperature *)
Do[{
tguess1=Bvp[4]/(Avp[4]-Log[press])+ 0.1;
tguess2=Bvp[3]/(Avp[3]-Log[press])- 0.1;
Do[{
Vp[i]= Avp[i]-Bvp[i]/tguess1;
Vp[i]=Exp[Vp[i]]
},{i,1,nc}];
cov1=press;
Do[
cov1=cov1-xc[i,j]*Vp[i]
,{i,1,nc}];
Do[{
Vp[i]= Avp[i]-Bvp[i]/tguess2;
Vp[i]=Exp[Vp[i]]
},{i,1,nc}];
cov2=press;
Do[
cov2=cov2-xc[i,j]*Vp[i]
,{i,1,nc}];
cov3=press;
While[Abs[cov3]>=1.0*10^-8,
{tguess3=(tguess1+tguess2)/2;
Do[{
Vp[i]= Avp[i]-Bvp[i]/tguess3;
Vp[i]=Exp[Vp[i]]
},{i,1,nc}];
cov3=press;
Do[
cov3=cov3-xc[i,j]*Vp[i]
,{i,1,nc}];
If[cov3*cov1>=0,
{cov1=cov3;tguess1=tguess3},{cov2=cov3;tguess2=tguess3}];
}];
tc[j]=tguess3;
(* Vapor phase composition calculation *)
ysum=0;
Do[
ysum=ysum+alpha[i]*xc[i,j]
,{i,1,nc}];
Do[
yc[i,j]=alpha[i]*xc[i,j]/ysum,
{i,1,nc}];
},{j,1,nt}];
(* liquid flow rate calculation *)
Reaction[];
Lc[1]=xtry[(nc-1)*nt+1];vc[nt]=xtry[(nc-1)*nt+2]/hov;
Do[Lc[j]=Lc[1],{j, 2, nrt-1}];
Do[Lc[j]=Lc[j-1]-rate[j]*hor/hov,{j, nrt, nrb}];
Do[Lc[j]=Lc[j]+qfeed[nfeed1]*feedflow[nfeed1],{j, nfeed1, nrb}];
Do[Lc[j]=Lc[j]+qfeed[nfeed2]*feedflow[nfeed2],{j, nfeed2, nrb}];
Do[Lc[j]=Lc[nrb],{j, nrb+1, nt-1}];
Lc[nt]=Lc[nt-1]-vc[nt];
(* Vapor flow rate calculation *)
Do[vc[j]=vc[nt],{j, nt-1, nrb+1, -1}];
Do[vc[j]=vc[j+1]+rate[j]*hor/hov,{j, nrb, nrt, -1}];
Do[vc[j]=vc[j]+(1-qfeed[nfeed2])*feedflow[nfeed2],{j, nfeed2, nrt, -1}];
Do[vc[j]=vc[j]+(1-qfeed[nfeed2])*feedflow[nfeed2],{j, nfeed1, nrt, -1}];
Do[vc[j]=vc[nrt],{j, nrt-1, 2, -1}];
vc[1]=vc[2];
];

Reaction[]:=Module[{kf, kb},
Do[{rate[j]=0;
Qrea[j]=0;
},{j,1,nt}];
Do[{
kf=af*Exp[-acf/(R1*tc[j])];
kb=ab*Exp[-acb/(R1*tc[j])];
rate[j]=holdup[j](kf*xc[1,j]*xc[2,j] - kb*xc[3,j]*xc[4,j]);
Qrea[j]=rate[j]*hor;
},{j, nrt, nrb}]
];

NewRap[]:= Module[{kk, ptimes, xtry1,
k, m, x1, JAC, JAC1, com1, com2, ext,
nom, er, er1, cov, jvalue, dx},
nstep=0;mstep=1;ptimes=1;
(*
(* Initial Guess *)
xtry[(nc-1)*nt+1]=0.02;
xtry[(nc-1)*nt+2]=200;
xtry[2*nt+1]=Tc;
xtry[1]=xtry[nt+1]=(1-xtry[2*nt+1])/3;
xtry[nt]=xtry[2*nt]=xtry[3*nt]=(1-Bc)/3;
Do[
Do[
xtry[(i-1)*nt+k]=xtry[(i-1)*nt+1]-
(xtry[(i-1)*nt+1]-xtry[(i-1)*nt+nt])/(nt-1)*(k-1)
,{k,2,nt-1}]
,{i,1,nc-1}];
*)
init=ReadList["xstatic.nb", Number];
Do[xtry[k1]=Part[init , k1], {k1, 1, (nc-1)*nt+2}];
Do[
Do[
xc[i,k]=xtry[(i-1)*nt+k]
,{i,1,nc-1}]
,{k,1,nt}];
xtov[];
nstep=nstep+1;
Print["The Number of Jacobi Inversion Calcuation:mstep = ", mstep];
Print["The Number of XtoV Calcuation:nstep = ", nstep];
Mbalance[];
Do[er[k]=error[k],{k,1,(nc-1)*nt+2}];
cov=Max[Table[Abs[er[k]],{k,1,(nc-1)*nt+2}]];
Print["cov=", cov, "  ","Reflux rate =", xtry[(nc-1)*nt+1],
"  ","Qreb = ", xtry[(nc-1)*nt+2]];
While[cov>=1*10^-8,{
Do[{xtry1[m]=xtry[m];
If[m>=(nc-1)*nt+1,xtry[m]=xtry[m]*0.98,
xtry[m]=xtry[m]*0.98];
Do[
Do[
xc[i,k]=xtry[(i-1)*nt+k]
,{i,1,nc-1}]
,{k,1,nt}];
xtov[];
nstep=nstep+1;
Mbalance[];
Do[er1[k]=error[k],{k,1,(nc-1)*nt+2}];
Do[{com1[k,m]=com2[k,m]=
(er[k]-er1[k])/(xtry1[m]-xtry[m]);
},{k,1,(nc-1)*nt+2}];
xtry[m]=xtry1[m]
},{m,1,(nc-1)*nt+2}];
JAC1=Table[com1[k,m],{k,1,(nc-1)*nt+2},{m,1,(nc-1)*nt+2}];
JAC1=Inverse[JAC1];
jvalue=Det[JAC1];
Do[{
Do[com2[k,m]=-er[k],{k,1,(nc-1)*nt+2}];
nom=Table[com2[k,kk],{k,1,(nc-1)*nt+2},{kk,1,(nc-1)*nt+2}];
dx=jvalue*Det[nom];
ext=0.0001;
If[mstep>=200,ext=0.005];
If[mstep>=220,ext=0.01];If[mstep>=240,ext=0.02];
If[mstep>=260,ext=0.04];If[mstep>=280,ext=0.1];
xtry[m]=xtry[m]+ext*dx;
If[m<=(nc-1)*nt && xtry[m]<=0,xtry[m]=0.000001];
If[m<=(nc-1)*nt && xtry[m]>=1.0,xtry[m]=xtry[m]-ext*dx];
If[m>=(nc-1)*nt+1 && xtry[m]<=0,{}];
Do[com2[kk,m]=com1[kk,m],{kk,1,(nc-1)*nt+2}]
},{m,1,(nc-1)*nt+2}];
Do[
Do[
xc[i,k]=xtry[(i-1)*nt+k]
,{i,1,nc-1}];
,{k,1,nt}];
xtov[];
nstep=nstep+1;mstep=mstep+1;
Mbalance[];
Do[er[k]=error[k],{k,1,(nc-1)*nt+2}];
cov=Max[Table[Abs[er[k]],{k,1,(nc-1)*nt+2}]];
ptimes=ptimes+1;
Print["The Number of XtoV Calcuation:nstep=", nstep,"  The Number of Jacobi Inversion Calcuation: mstep = ", mstep];
If[ptimes>=30,{
Print["xc[1,j]      ", "xc[2,j]      ", "xc[3,j]      ", "xc[4,j]      ","Lc[j]      ", "tc[j]", "     rate[j]"];
Do[
Print[xc[1,k], "    ",xc[2,k],"    ", xc[3,k], "    ",xc[4,k], "     ",Lc[k],"    ",tc[k], "       ", rate[k]]
,{k,1,nt}];
Print[];
Print["yc[1,j]      ", "yc[2,j]      ", "yc[3,j]      ", "yc[4,j]      ", vc[j]];
Do[
Print[yc[1,k],"    ", yc[2,k], "    ",yc[3,k],"    ", yc[4,k],"     ",vc[k]]
,{k,1,nt}];
Print[];ptimes=0}];
Print["cov=", cov, "  ","Reflux Rate =", xtry[(nc-1)*nt+1],
"  ","Qreb = ", xtry[(nc-1)*nt+2]];
}];
rmax=cov;
Do[{vapor[k]=vc[k];
liquid[k]=Lc[k];
temp[k]=tc[k];
Do[{
xcomp[i,k]=xc[i,k];
ycomp[i,k]=yc[i,k];
},{i,1,nc}];
},{k,1,nt}];
RR=xtry[(nc-1)*nt+1];
Qcon=vapor[1]*hov;
Qreb=xtry[(nc-1)*nt+2];
];

Mbalance[]:=Module[{k    },
Reaction[];
Do[
{error[(i-1)*nt+1]=Lc[1]*xc[i,1]+(vc[1]-xtry[(nc-1)*nt+1])*yc[i,1]
-vc[2]*yc[i,2];

Do[error[(i-1)*nt+k]=Lc[k]*xc[i,k]+vc[k]*yc[i,k]
-Lc[k-1]*xc[i,k-1]-vc[k+1]*yc[i,k+1]
,{k, 2, nfeed1-1}];

error[(i-1)*nt+nfeed1]=Lc[nfeed1]*xc[i,nfeed1]+vc[nfeed1]*yc[i,nfeed1]
-Lc[nfeed1-1]*xc[i,nfeed1-1]-vc[nfeed1+1]*yc[i,nfeed1+1]-feedflow[nfeed1]*fcomp[i, nfeed1];

Do[error[(i-1)*nt+k]=Lc[k]*xc[i,k]+vc[k]*yc[i,k]
-Lc[k-1]*xc[i,k-1]-vc[k+1]*yc[i,k+1],{k,nfeed1+1,nfeed2-1}];

error[(i-1)*nt+nfeed2]=Lc[nfeed2]*xc[i,nfeed2]+vc[nfeed2]*yc[i,nfeed2]
-Lc[nfeed2-1]*xc[i,nfeed2-1]-vc[nfeed2+1]*yc[i,nfeed2+1]-feedflow[nfeed2]*fcomp[i, nfeed2];

Do[error[(i-1)*nt+k]=Lc[k]*xc[i,k]+vc[k]*yc[i,k]
-Lc[k-1]*xc[i,k-1]-vc[k+1]*yc[i,k+1],{k,nfeed2+1,nt-1}];

error[(i-1)*nt+nt]=Lc[nt]*xc[i,nt]+
vc[nt]*yc[i,nt]-Lc[nt-1]*xc[i,nt-1]
},{i,1,nc-1}];
Do[
Do[
If[i<=2, error[(i-1)*nt+k]=error[(i-1)*nt+k]+rate[k],error[(i-1)*nt+k]=error[(i-1)*nt+k]-rate[k]]
,{i, 1, nc-1}]
,{k, nrt, nrb}];
error[(nc-1)*nt+1]=Tc-yc[3,1];
error[(nc-1)*nt+2]=1-Bc-xc[1,nt]-xc[2,nt]-xc[3,nt];
];

(* Main Program *)
propdata[];
(* Number of components *)
nc=4;
nrt=nfeed1-COPROD1;nrb=nfeed2+COPROD2;
nt=nfeed2+ns+1;
Print["nfeed1 = ",nfeed1," nfeed2 = ",nfeed2,
" nrt = ",nrt, " nrb = ",nrb," nt = ",nt];
(* Operating Pressure, Ba *)
press=9;
(* Stage holdup, kmol *)
Do[holdup[j]=1.0*6/(nrb-nrt+1), {j, 1, nt}];
(* Feed information *)
Do[Do[fcomp[i,j]=0,{i, 1, nc}], {j, 1, nt}];
Do[feedflow[j]=qfeed[j]=0, {j, 1, nt}];
(* The Feed F1 information *)
fcomp[1,nfeed1]=0.0;fcomp[2,nfeed1]=1.0;
fcomp[3,nfeed1]=0.0;fcomp[4,nfeed1]=0.0;
qfeed[nfeed1]=1.0;
feedflow[nfeed1]=0.0126;
(* The Feed F2 information *)
fcomp[1,nfeed2]=1.0;fcomp[2,nfeed2]=0.0;
fcomp[3,nfeed2]=0.0;fcomp[4,nfeed2]=0.0;
qfeed[nfeed2]=1.0;
feedflow[nfeed2]=0.0126;
Print["The beginning of calculation"];
NewRap[];
Print["SUCCESS of NEWTON-RAPHSION ITERATIVE SOLUTION"];
Print["The Number of Jacobi Inversion Calcuation: mstep = ", mstep];
Print["The Maximum Error of Component Balance:rmax=",rmax];
Print[" Reflux Rate = ", RR ];
Print["Qcon = ",Qcon,"   Qreb = ",Qreb];
Print[];
Print["xcomp[1,j]        ", "xcomp[2,j]        ",
"xcomp[3,j]        ", "xcomp[4,j]"];
Do[
Print[xcomp[1,j], "        ", xcomp[2,j], "        ",
xcomp[3,j], "        ", xcomp[4,j]],
{j,1,nt}];
Print[];
Print["ycomp[1,j]        ", "ycomp[2,j]        ",
"ycomp[3,j]        ", "ycomp[4,j]"];
Do[
Print[ycomp[1,j], "        ", ycomp[2,j], "        ",
ycomp[3,j], "        ", ycomp[4,j]],
{j, 1, nt}];
Print[];
Do[
{Print["temp[",j,"]=",temp[j],"   Qrea[",j,"]=",Qrea[j]];
},{j, 1, nt}];
Print[];
Do[Print["rate[",j,"]=", rate[j]],{j, 1, nt}];
Print[];
sumR=0;
Do[sumR = sumR + rate[j],{j, 1, nt}];
Print["sumR = ", sumR];
Do[Print["vapor[",j,"]=",vapor[j]],{j, 1, nt}];
Print[];
Do[Print["liquid[",j,"]=",liquid[j]],{j, 1, nt}];
Print[];
Do[
Do[{xcomp[i,j]>>>xfile;ycomp[i,j]>>>yfile},{i,1,nc}]
,{j, 1, nt}];
Do[{temp[j]>>>tfile;Qrea[j]>>>qrfile;
vapor[j]>>>vfile;liquid[j]>>>lfile
},{j, 1, nt}];
Do[xtry[k]>>>xstatic,{k, 1, (nc-1)*nt+2}];
];
RDSIM[9, 12, 7, 0.95, 0.95, 2, 5]

smarten

这么长的程序，应该把nb文件压缩成zip或者rar之后放上来比较好。你这程序不能运行，不知道是不是我没拷贝好。
另外就是这个程序没用Mathematica的知识，大概应该叫C程序吧。我觉得你可以学习一下Mathematica，然后把这些程序重写大概总共20行。 有空再学习一下FreddyMusic的demo， 做一个像样的仿真，有计算有图形。

4everlove

2# smarten

4everlove

是吗？？我不懂啊！！导师知识给了我们这个程序说是mathematica进行仿真

smarten

没看到你说什么

4everlove

这次看到了吗？？我们是要用mathematica啊！您觉得我这个程序不是mathematica语言吗？？反正我觉得不是c语言！

smarten

我是说写这个程序的人没怎么学过Mathematica。程序不难，你大概花一个上午就可以看明白了吧。不然你们导师要不高兴了，他花了半个月写的程序，你花一分钟就搞通了，太没面子了吧。

4everlove

呵呵，那倒不至于，这个程序知识我们毕设的一部分，下周三就要去找老师了，我想看明白或者带着点什么问题就找他！可问题是我看的不太明白啊！编程我没有什么基础啊，就大一学过c语言，还马马虎虎！你觉得看懂这个程序一个没什么基础的人现在要看什么！

smarten

你大概是看看书上的原理，自己重写。这个程序本身不全，不能运行。 xc 没有初值。或者你让你们老师给你们运行一次。

4everlove

您能给我些mathematic语言的一些基本语法的资料嘛？？最好有些例子，我在网上找了但似乎没有我想要的那种！！谢谢谢谢！！！

smarten

你要找到像你们老师那么用Mathematica的还真不容易。如果你要学mathematica的话，这个论坛上教程什么的不少，或者参考http://forum.simwe.com/thread-779926-1-1.html

4everlove

嗯，好的！谢谢你了。等我以后再有问题了再问你吧！
这个论坛高手不少啊！

4everlove

你一般什么时候在线啊？？

marveloustau

没什么难的，基本就和看C语言一样看就行了……你看不懂的原因也可能是因为程序太乱所以被吓到。如果把一开始的那一坨注释抄在每个变量在程序里第一次出现的位置就会好懂多了。

4everlove

有可能吧！但是我c语言那会学的就不太好。后来也没怎么接触过，就一点点基础了！比如下面这段程序，这个循环我就有点乱，您能替我解释一下嘛

Do[
{xc[4,j]=1-xc[1,j]-xc[2,j]-xc[3,j];
If[xc[4,j]<=0,
{xc[4,j]=0;
Do[xc[i,j]=xc[i,j]/(xc[1,j]+xc[2,j]+xc[3,j]+xc[4,j]),{i,1,nc}]}
]
},{j,1,nt}];
Do[
Do[
xtry[(i-1)*nt+k]=xc[i,k]
,{i,1,nc-1}];
,{k,1,nt}];