Программа оптимизации теплообменного аппарата
Для разработки предложенной программы была использована среда PascalABC совместно с офисным пакетом Microsoft Excel 2010.
Текст программы
Programsimone;
Constdx=0.0001;
Q=10000000; //Производительность
Tv1=64; //Начальная температура воды
Tv2=138; //Конечная температура воды
Tn=175.4; //Температура пара
R_zagr=0.00015; //Степень загрязнения
teta=0.98; //Коэффициент тепловых потерь
Lam_st=43; //Коэффициент теплопроводности
Delta_st=0.0025; //Толщина стенки
//Ограничения
dmax=10;
dmin=0.00000001;
Wmax=10;
Wmin=0.0000001;
//Параметры стенки
Pr_c=1.278;
Kw=0.0001;
Tc=138.2;
Lam_c=0.685;
mu_c=204.106/1000000;
L_tr=3;
//Параметры воды
Ro=957.66;
Cp=4221;
Lam_v=0.6832;
v=0.293/1000000;
Pr=1.735;
//Параметры пара
r1=2030.4;
mu_s=157.5/1000000;
Lambda_s=0.676;
Pr_s=1.023;
Ro_s=4.66;
w11=0.162;
//Цена
Cena_F=2000;
Cena_ElEn=5;
v_s=0.173/1000000;
dT=37.7;
H1=3;
typearray2D = array of array ofreal;
Varx,x0,df,g :array ofreal; d2f:array2D;
eps,h,grad,s,r,fc,P,c :real;
i,j,n,n1,k,l :integer;
Alfa_vod,K_tp,F_ta,dP,N_nasosa,Zatr_F,
Zatr_ElEn,Delta_t,G_vod,Lambda,Alfa_para,a,an,eps1,Q_ta,HdTkr: real;
Functionbarier(g:array ofreal): real;
labelstop;
Vari:integer;
s:real;
Begin
s:=0;
barier:=1e30;
Fori:=1 ton1 do
ifg[i]<0 thens:=s+1/g[i] else gotostop;
barier:=-r*s;
stop: end;
Functionf(x: array ofreal):real;
Begin
eps1:=Power((Power(Lam_c/Lambda_s ,3)*(mu_c/mu_s)),0.125);
an:=0.725*Power( ( Power(Lam_v,3)*9.81*(Ro-Ro_s)*r1*1000)/(v*x[1]*(Tn-Tc)) ,0.25 )*eps1;
a:=25.7*Power((Ro_s*w11*w11)/(9.81*Ro*x[1]),0.08)*Power((an*x[1]/Lam_v),(-0.5))*an;
Alfa_para:=(a*0.84)/Power(10,0.07);
Delta_t:=((Tn-Tv1)-(Tn-Tv2))/ln((Tn-Tv1)/(Tn-Tv2));
Alfa_vod:=(0.021*Power((x[2]*x[1]/v),0.8)*Power(Pr,0.43)*Power((Pr/Pr_c),0.25)*Lam_v)/x[1];
K_tp:=1/(1/Alfa_para+Delta_st/Lam_st+1/Alfa_vod+R_zagr);
Q_ta:=Q*teta;
F_ta:=Q_ta/(K_tp*Delta_t);
Lambda:=0.11*Power((Kw/x[1])+(v*68/x[2]*x[1]),0.25);
dP:=Lambda*(L_tr/x[1])*(x[2]*x[2]/2)*Ro+7.5*(Ro*x[2]*x[2]/2);
G_vod:=Q_ta/(Cp*(Tv2-Tv1));
N_nasosa:=(G_vod*dP)/(Ro*0.85);
Zatr_F:=F_ta*Cena_F;
Zatr_ElEn:=(N_nasosa*8400/1000)*Cena_ElEn;
fc:=Zatr_F+Zatr_ElEn;
//Ограничения
g[1]:=-dmax+x[1];
g[2]:=-x[1]+dmin;
g[3]:=-Wmax+x[2];
g[4]:=-x[2]+Wmin;
P:=barier(g);
F:=fc+P;
End;
//Процедураобращенияматрицы
Procedureinvert(n,q: integer; matr1:array2D; varmatr: array2D);
labelM1,M2;
Var
a:array of array ofreal;
i,j,k,m : integer;
t :real;
Begin
SetLength(a, n + 1);
Fori := low(a) tohigh(a) do
SetLength(a[i], 2*(n+1));
m:=2*n; q:=0;
Fori:=1 ton do
Forj:= 1 tom do
Ifj<=n thena[i,j]:=matr1[i,j] else
Ifj=n+i thena[i,j]:=1.0 elsea[i,j]:=0;
Fori:=1 ton do
Begink:=i;
M1: ifa[k,i]=0 then
Beginq:= 1;
ifk<n thenk:=k+1 else gotoM2;
gotoM1;
End;
Ifq=1 then
Forj:=1 tom do
Begin
t:=a[k,j]; a[k,j]:=a[i,j]; a[i,j]:=t
End;
Forj:=m downtoi doa[i,j]:=a[i,j]/a[i,i];
Fork:= 1 ton do
Ifk<>i then
Forj:=m downto1 do
a[k,j]:=a[k,j]-a[i,j]*a[k,i];
End; { i }
q:=0;
Fori:= 1 ton do
Forj:= 1 ton domatr[i,j]:=a[i,j+n];
M2: End; {invert}
// Процедура вычисления координат вектора градиента и формирования матрицы вторых производных
ProcedureFor_Mat_d2f(n:integer; vargrad:real);
Var
i,j :integer;
s,f0: real;
Begin
//Расчет первых производных
f0:=f(x); s:=0;
Fori:= 1 ton do begin
x[i]:=x[i]+dx;
df[i]:=(f(x)-f0)/dx;
s:=sqr(df[i]);
x[i]:=x[i]-dx;
End;
grad:=sqrt(s);
//Расчет вторых производных
Fori:=1 ton do
Begin
s:=-2*f(x);
x[i]:=x[i]+dx;
s:=s+f(x);
x[i]:=x[i]-2*dx;
s:=s+f(x);
x[i]:=x[i]+dx;
d2f[i,i]:=s/sqr(dx);
End;
//Расчет смешанных производных
Fori:= 1 ton-1 do
Forj:=i+1 ton do
Begin
s:=f(x); // 1
x[i]:=x[i]-dx; x[j]:=x[j]-dx;
s:=s+f(x); //4
x[j]:=x[j]+dx;
s:=s-f(x); //2
x[i]:=x[i]+dx; x[j]:=x[j]-dx;
s:=s-f(x); //3
x[j]:=x[j]+dx;
d2f[i,j]:=s/sqr(dx);
d2f[j,i]:=d2f[i,j];
End;
End; // For_Mat_d2f
Proceduretek_koord;
Vari,j :integer;
s: real;
Begin
Repeat
For_Mat_d2f(n,grad);
invert(n,q,d2f,d2f);
Ifq = 1 then
Begin
writeln('Определительравеннулю');
exit{goto stop}
End;
Fori:= 1 ton do
Begin
s:=0;
Forj:= 1 ton do
s:=s+d2f[i,j]*df[j];
x[i]:=x[i]-s;
End;
Untilgrad<eps;
End; // tek_koord
Procedureprint;
Vari : integer;
Begin
Writeln('Итерация ',k, ' Параметрштрафаr=',r: 12:8);
Fori:=1 ton do
Begin
Writeln('x0[',i,'] = ',x0[i], 'x[',i,'] = ',x[i]);
x0[i]:=x[i]
End;
writeln('Функцияштрафа = ',P);
writeln
End;
Begin
Writeln('Исходныеданные');
Writeln;
Write('Введите размерность задачи оптимизации n= ');
Readln(n);
Write('Введите точность вычислений eps= ');
Readln(eps);
Writeln('Введите начальные значения переменных');
SetLength(x,n + 1);
SetLength(x0,n+1);
SetLength(df,n+1);
SetLength(g,n + 3);
SetLength(d2f, n + 1);
Fori := low(d2f) tohigh(d2f) do
SetLength(d2f[i], n+1);
Fori:= 1 ton do
Begin
Write('x[',i,']=');
Readln(x[i]);
End;
Write('Введите число ограничений = ');
Readln(n1);
Write('Введите начальное значение параметра штрафа r=');
Readln(r);
Write('Введите число уменьшения параметра штрафа с=');
Readln(c);
Writeln;
k:=0;
Repeat
k:=k+1;
tek_koord;
print;
r:=r/c;
Untilabs(P)<0.01;
Writeln;
Writeln('Peзyльтaтыоптимизации');
Writeln('c использованием обратной штрафной функции:');
Writeln;
Fori:= 1 ton dowriteln('x[',i,'] = ',x[i]);
Writeln;
Writeln('Значение функции цели = ',f(x));
Writeln;
Writeln('Значение функции цели = ',f(x));
Writeln;
Writeln('Значение функции цели = ',f(x));
Writeln;
Writeln('Оптимальные значения ');
Writeln;
Writeln('диаметр d = ',x[1]:2:3,' м');
Writeln;
Writeln('скоростьw = ',x[2]:2:3,' м3/с');
Writeln;
Writeln('ЗатратыЗ= ',f(x):8:3,' руб');
Writeln;
Writeln('ПлощадьтеплообменногоаппаратаF= ',F_ta:5:3,' м2');
Readln;
End.
Результаты оптимизации данной программы быть представлены на рисунке 1.2.
Рис.1.2. - Результаты оптимизации теплообменного аппарата
З=f(dвн,w)
Рис.1.3. – Изменение затрат
В программе Surfer 13 был построен график зависимости скорости в теплообменном аппарате от диаметра его труб.
Рис 1.4. – Зависимость диаметра от скорости теплоносителя с ограничениями
Из результатов оптимизации функции затрат для рассматриваемого аппарата видно, что минимум (З=140607.015 руб) достигается при dвн=17.915 мм, F=63.477 
и скорости теплоносителя в трубах w=1.185 м/с.
По каталогу был выбран стандартный теплообменный аппарат, тип ихарактеристики данного аппарата представлены в таблице 1.1:
Таблица 1.1. - Характеристики теплообменного аппарата
| Тип аппарата | Диаметр кожуха, мм | Условное давление, МПа | Число ходов по трубам | Наружный диаметр труб, мм | Площадь поверхности теплообменника, м2 | Площадь проходного сечения одного хода по трубам, м2 | Площадь проходного сечения по межтрубному пространству, м2 | ||||
Наружный  | Внутренний  | В кожухе | В трубах | При длине прямого участка 3000 мм | При толщине стенки труб, мм | В вырезе перегородки | Между перегородками | ||||
| 1,8 | |||||||||||
| ТН | 1,6 | 1,6 | 0,039 | 0,037 | 0,037 | 0,048 |