Программа оптимизации теплообменного аппарата

Для разработки предложенной программы была использована среда 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