Список литературы

1. Теплофизические свойства криопродуктов / Л.А Акулов., Е.И. Борзенко, В.Н. Новотельнов, А.В. Зайцев. – СПб.: Политехника, 2001. – 242 с.

2. Архаров А.М., Беляков В.П.и др. Криогенные системы. Основы проектирования аппаратов, установок и систем. – М.: Машиностроение, 1999. – 719 с.

3. Борзенко Е.И., Зайцев А.В.Установки и системы низкотемпературной техники. Автоматизированный расчет и моделирование процессов криогенных установок и систем: Учеб. пособие. – СПб.: СПбГУНиПТ, 2006. – 232 с.

4. Алексеев В.П., Вайнштейн Т.Е., Герасимов П.В.Расчет и моделирование аппаратов криогенных установок. – Л.: Энергоатомиздат, 1987. – 277 с.

ПриложениЕ 1 Текст подпрограммы теплового и гидравлического расчета пластинчато-ребристого конденсатора-испарителя PlRkiCalc

!Подпрограмма расчета пластинчато-ребристого конденсатора-испарителя

subroutine PlrkiCalc(trez)

!1. ––––––––––––––––––––––––––––––––––––––– Исходные данные ––––––––––––––––––––––––––––––––––––––––––––––––––––––

!свойства смеси компонентов

common /sv4/TS,AMs,AMps,rs,rsk,ros,rops,ss,amus,amups,als,cps

!исходные данные

common/dan/Q,pi,pk,xi1,xi2,yi1,yi2,yk1,yk2,Rsl,A,w0, &

hk,hst,co,hp,b1,b2,b3,hli,hlk,si,sk,di,dk,dst,db

!результаты

common/rez/Ti,Re,pis,TSi,dTi,ali,qei,Tk,TSk,dTk,alk,qek,dTst,QQ,dTs,Fke,nk,nk1,np, &

w0vp,w0v,wop,fop,hop,pdv,dpe,dp,dpd,dpu,dpv,dps,dppod

character(80)trez

DIMENSION xxx(20),fff(20) кусочно-линейная аппроксимация (рис.6)

DATA xxx/.2,.3,.4,.5,.6,.7,.8,1.,1.5,2.1,3.,4.,5.,6.,7.,8.,9.,10.,20.,30./

DATA fff/14.2,9.1,7.,6.08,5.27,4.62,4.06,3.66,3.2,2.85, &

2.6,2.24,2.04,2.,1.85,1.75,1.67,1.58,1.5,1.48/

call smes !константы для расчета свойств смеси

alst=145. !теплопроводность материала стенки

eps=1.e-4 !погрешность расчетов

!2. ––––––––––––––––––––––––––––––-–– Геометрические параметры ––––––––––––––––––––––––––––––––––––––––––––

dei=2*(si-di)*(hli-di)/(si+hli-2*di) !эквивалентный диаметр каналов испарителя

hr=hp-hk-b3 !рабочая высота пакета (в т.ч. парогенерирующего канала)

Fi=2*hr*(b1-2*b3)*(si+hli-2*di)/si !площадь раб. поверхности канала испарителя

Fk=2*hr*(b1-2*b3)*(sk+hlk-2*dk)/sk ! - " - канала конденсатора

Fri=2*hr*(b1-2*b3)*(hli-di)/si ! - " - ребер канала испарителя

Frk=2*hr*(b1-2*b3)*(hlk-dk)/sk ! - " - ребер канала конденсатора

Fsi=(b1-2*b3)*(si-di)*(hli-di)/si !площадь поперечного сечения канала испарителя

!3. ––––––––––––––––––––––––––––––––––––- Расчет канала конденсации ––––––––––––––––––––––––––––––––––––––––

! Объемные доли компонентов в конденсате на входе

xk1v=yk1; xk2v=yk2 !начальное приближение xк1в и xк2в

x1=0.; x2=0. !начальное xк1в и xк2в с предыдущей итерации

do while (abs(x1-xk1v)>=eps.or.abs(x2-xk2v)>=eps) !заголовок цикла, погрешн.

call olet(pk,xk1v,xk2v,a1v,a2v) !относительные летучестеи a1в, a2в

x1=xk1v; x2=xk2v !старые значения xк1в и xк2в

cc=(yk1/a1v+yk2/a2v+1-yk1-yk2)

xk1v=(yk1/a1v)/cc; xk2v=(yk2/a2v)/cc !новые значения xк1в и xк2в

enddo !конец цикла

!4. Промежуточный анализ результатов

call dolm(xk1v,xk2v,xxk1,xxk2) !массовые доли комп. в конденсате на входе

call temp(pk,Tkv,xk1v,xk2v) !температура конденсации в верхнем сечении Tкв

call ssmes(pk,xk1v,xk2v) !свойства смеси в конденсате на входе

!5. Объемные доли компонентов в конденсате на выходе

xk1n=yk1; xk2n=yk2 !начальное приближение xк1н и xк2н

x1=0.; x2=0. !начальное xк1н и xк2н с предыдущей итерации

do while (abs(x1-xk1n)>=eps.or.abs(x2-xk2n)>=eps) !заголовок цикла, погрешн.

do while (abs(x1-xk1v)>=eps.or.abs(x2-xk2v)>=eps) !заголовок цикла, погрешн.

call olet(pk,xk1n,xk2n,a1n,a2n) !относительные летучести a1н, a2н

x1=xk1n; x2=xk2n !старые значения xк1н и xк2н

cc=a1n*x1+a2n*x2+1-x1-x2

c1=a1n*cc; c2=a2n*cc

cc=Q+rsk*A

xk1n=yk1*cc/(Q+c1*rsk*A) !новое значение xк1н

xk2n=yk2*cc/(Q+c2*rsk*A) !новое значение xк2н

enddo

!6. Средний состав пленки конденсата

xk1=(xk1v+xk1n)/2.; xk2=(xk2v+xk2n)/2. !объемные доли

call dolm(xk1,xk2,xxk1,xxk2) !массовые доли

!7. Средняя температура и свойства пленки конденсата

call temp(pk,Tk,xk1,xk2)

call ssmes(pk,xk1,xk2)

TSk=TS; AMsk=AMs; AMpsk=AMps; rosk=ros; ropsk=rops; ssk=ss;amusk=amus

cdTk=(1./als)*(amus/ros/9.81)**.333*(hr/(rsk*amus))**.282 !комплекс для расчета dTк

!8. Цикл по эффективному q при конденсации

qeki=800.; qeka=4500. !мин и макс q

do while(abs(qeka-qeki)>.001) ) !––––––––––––-начало цикла по qek–––––––––-

qek=.5*(qeki+qeka)

dTk=(qek/.95)**1.282*cdTk !темп. напор на стороне конденсации

alk=qek/dTk !коэф. теплоотдачи конд.

amk=sqrt(2*alk/(alst*dk)) !параметр ребра кан.конд.

pdrk=.5*amk*hlk; pdrk=tanh(pdrk)/pdrk !КПД ребра канала конденсации

pdk=1-(1-pdrk)*Frk/Fk !КПД поверхности конденсации

QQ=qek*Fk*pdk !кол-во выд.теплоты в кан.конд.

dTst=QQ*dst/(2*(b1-2*b3)*hr*alst) !температурный напор в проставочном листе

!9. ––––––––––––––––––––––––––––––––––––––– Расчет канала кипения ––––––––––––––––––––––––––––––––––––––––––––

eee=1.

qei=QQ/Fi !эфф.плотность Q в канале кип.

do while (eee>.001) !––––––––––––-начало цикла по qei–––––––––-

!массовые доли компонентов жидкости на входе

call dolm(xi1,xi2,xxi1,xxi2)

!10. Температура и свойства кипящей смеси на входе

call temp(pi,Ti,xi1,xi2)

call ssmes(pi,xi1,xi2)

Re=w0*dei*ros/amus !число Рейнольдса в испарителе

xxi=4.*hr*qei/(ros*rs*dei*w0) !массовая доля пара на выходе из испарителя

ckr=1./xxi !кратность циркуляции

Go=Rsl*rs*1.e3/AMs/Q ! отн. расх. ж. на вых. из кан.кип.

! Объемные доли компонентов на выходе и входе в каналы кипения

call olet(pi,xi1,xi2,a1,a2)

cc=a1*xi1+a2*xi2+1.-xi1-xi2

c1=a1/cc; c2=a2/cc

cc=ckr*(1.+Go)

xi1v=xi1*cc/(ckr*c1+Go*(ckr+c1-1.)) !объемная доля на выходе из испарителя

xi2v=xi2*cc/(ckr*c2+Go*(ckr+c2-1.)) !- " -

xi1n=xi1v*(ckr+c1-1.)/ckr !объемная доля на входе в испаритель

xi2n=xi2v*(ckr+c2-1.)/ckr !- " -

!11. Усредненные параметры смеси в каналах испарителя

h=6100. !

pis=pi+.5e-4*h*ros*hr*9.81 !среднее давление в каналах кипения

xi1s=.5*(xi1n+xi1v); xi2s=.5*(xi2n+xi2v) !средние объемные доли комп. в жидкости

call temp(pis,Tis,xi1s,xi2s)

call ssmes(pis,xi1s,xi2s)

RM=w0*Fsi*ros !массовый расход жидкости на входе в парогенератор

Rem=qei*dei*ros/(rs*amus*rops) !модифиц. Re потока массы от теплообм.пов-ти

!12. Температурный напор на стороне кипения

dTi=.15*QQ*Re**1.1*(dei/hr)**1.1/(RM*cps*Rem**.75) !

ali=qei/dTi !коэффициент теплоотдачи при кипении

ami=sqrt(2*ali/(alst*di)) !параметр ребра канала кипения

pdri=.5*ami*hli; pdri=tanh(pdri)/pdri !КПД ребра канала кипения

pdi=1-(1-pdri)*Fri/Fi !КПД поверхности кипения

!13. Новое значение плотности теплового потока

qein=QQ/Fi/pdi

eee=(qein-qei)/qein

qei=qein

enddo !––––––––––––––––––окончание цикла по qei–––––––––––-

! свойства тройной смеси в испарителе

TSi=TS; AMsi=AMs; AMpsi=AMps; rosi=ros; ropsi=rops

rsi=rs; ssi=ss; amusi=amus; amupsi=amups

!14. ––––––––––––––––––––––––––––––––––––––––– Баланс температур ––––––––––––––––––––––––––––––––

call temp(pi,Tir,xi1v,xi2v)

if(xi1<.965) Tir=Tis

dTr=Tk-Tir !располагаемый температурный напор

dTs=dTk+dTst+dTi !полученный суммарный температурный напор

if(dTs>dTr)then; qeka=qek; else; qeki=qek; endif !метод половинного деления

enddo !–––––––––––––––––-окончание цикла по qek–––––––

!15. ––––––––––––––––––––––––––––––––––-Параметры конструкции ––––––––––––––––––––––––––––––––––––––––––––––

Fke=Q/qek !площадь эффективной поверхности конденсатора

nk=Fke/(Fk*pdk)+1 !число каналов конденсации

nk1=(b2-(hli+2.*db))/(hli+hlk+2.*dst) !число каналов конденсации в пакете

np=nk/nk1+1 !число пакетов в конденсаторе–испарителе

!16. ––––––––––––––––––––––––––– Гидравлика, опускная система –––––––––––––––––––––––––––––––––––––––––––

w0vp=4*hr*qei/(dei*rsi*ropsi) !приведеная скорость пара на вых. из испарителя

w0v=w0*(1.-(w0vp/w0)*(ropsi/rosi)) !приведеная скорость смеси на вых. испарит.

beta=w0vp/(w0v+w0vp) !объемное расходное паросодержание на вых. испарит.

Vop=Fsi*w0v*(nk1+1.)*3600. !объемный расход жидкости в опускной системе

cmb=10000.; if((cmb*hst**1.5)<10.)then

cmb=6400.; if((cmb*hst**1.5)>5.)then

cmb=2800./(1.-720*hst**1.5); endif; endif !коэффициент расхода перелива

b=Vop/(cmb*hst**1.5) !параметр перелива

u0=.35*sqrt(9.81*(rosi-ropsi)*dei/rosi) !скорость всплытия пузырьков в исп.(fi<=.7)

c0=1.061*w0**(-.065) !параметр распределения (fi<=.7)

cc=1./(c0*(1.-ropsi/rosi))

w0p=.7*(c0*w0-u0)/(1.-.7/cc) !приведенная скорость пара при fi=.7

u01=.35 !скорость всплытия пузырьков в испарителе (fi>.7)

c01=1.03 !параметр распределения (fi>.7)

cc1=u0+c0*w0; cc2=u01+c01*w0

fis=cc*(1.-alog10(1.+w0p/(cc*cc1))*cc*cc1/w0vp) !средн.истинное паросодерж. исп.

if(w0p<w0vp)then; fis=fis+cc*(w0p/w0vp-1.)+cc*(1.-w0p/w0vp- & ! alog10((1.+w0v/(cc*cc2))/(1.+w0p*c01/(c0*cc*cc2)))*cc*cc1/w0vp) !паросодерж.

endif

ross=rosi-(rosi-ropsi)*fis !средняя плотность двухфазной смеси в испарителе

wp=1.5*(ssi*(rosi-ropsi)*(9.81/ross)**2)**.25 !скор.вспл.пузырей в большом объеме

cw=.9032 !

wop=cw*wp !скорость жидкости в опускной системе

fop=Vop/(3600.*wop) !площадь сечения опускной системы

!17. ––––––––––––––––––––––––––- Цикл – поиск hop методом половинного деления ––––––-–––––––––

hopi=.2; hopa=1.

dpp=10.

do while(abs(dpp)>.1)

hop=.5*(hopi+hopa) !относит.уровень жидколсти в ОП задаем =0,2...1,0

!18. ––––––––––––––––––––––––––––––––––––––– Баланс давлений –––––––––––––––––––––––––––––––––––––––––––––––––––––

!ode=f(hop) !f(hop) – аппроксимация графиков рис. 5.

ode=.3087 !относительная длина экономайзерного участка

pdv=hr*9.81*(rosi*(hop-ode)-ross*(1.-ode)) !движущее давление циркуляции

cg=.5 !коэффициент гидравл.сопротивления на входе в испаритель

cgm=1. !коэффициент истинного гидравлического сопр.на выходе из испарителя

dp=cg*rosi*w0**2/2. !сопротивление жидкости на входе в испаритель

if(Re>2300.)then;cla=.3164/Re**.25

else;cla=64./Re; endif !коэффициент гидравлического сопр.

dpt=.5*cla*rosi*w0**2/dei !сопротивление трения жидкости на единицу длины

dpe=dpt*hr*ode !сопротивление на экономайзерном участке

im=1000*xxi

if(Re>2300.)then;ck=.25;else;ck=1.; endif !коэфф. для вычисления xx и cks

cc=sqrt(ropsi/rosi)*(amusi/amupsi)**ck/2.; cks=0.

do i=1,im

z=i/1000.

xx=cc*(1./z-1.)**((2.-ck)/2.) !параметр Мартинелли

do j=2,20;if (xxx(j-1)<=xx.and.xxx(j)>=xx)then !аппроксимация (fi ж) рис. 6

fig=fff(j-1)+(fff(j)-fff(j-1))*(xx-xxx(j-1))/(xxx(j)-xxx(j-1))

endif; enddo

cks=cks+(1.-z)**(2.-ck)*fig**2. !интеграл в ф-ле (55)

enddo

cks=1.e-3*cks/xxi !Кси

dpd=dpe*cks*hr*(1.-ode) !сопротивление трения двухфазному потоку

cc=w0+w0vp*(1.-ropsi/rosi)

if(w0p>w0vp)then; cfi=w0vp/(c01*cc+u01) !истинное объемн.паросодерж.вых.исп.

else; cfi=w0vp/(c0*cc+u0); endif паросодержание на вых. исп.

dpu=(rosi*w0v)**2*(xxi**2/(cfi*ropsi)+(1.-xxi)**2/((1.-cfi)*ropsi)-1./rosi) !потери на уск.

dpv=cgm*rosi*w0**2*(1.+w0vp/w0v*(1.-ropsi/rosi))/2. !сопр. двухф.потоку на вых.

dps=9.81*rosi*hst !потери давления на преодоление столба жидкости

dppod=dp+dpe+dpd+dpu+dpv+dps

dpp=pdv-dppod !полезное давление циркуляции – сведение баланса

if(dpp>0.)then;hopa=hop; else;hopi=hop; endif

enddo !–––––––––––––––-конец цикла по hop (пока не достиенем dpp=0.)––––––––––-

19. ––––––––––––––––––––––––––––––––––––––––––––- Результаты расчета ––––––––––––––––––––––––––––––––––––––––––––

write(1,*);write(1,*)

write(1,*)' РЕЗУЛЬТАТЫ РАСЧЕТА'

write(1,*)

write(1,*)' КОНДЕНСАТОР'

write(1,*)' ==========='

write(1,*)'Объемн. доли компонент на входе: xk1v=',xk1v

write(1,*)' xk2v=',xk2v

write(1,*)'Относительные летучести на входе: a1v=',a1v

write(1,*)' a2v=',a2v

write(1,*)'Темп. конденсации в верхнем сеч.: Tkv=',Tkv

write(1,*)'Объемн. доли компонент на выходе: xk1n=',xk1n

write(1,*)' xk2n=',xk2n

write(1,*)'Относительные летучести на выходе: a1n=',a1n

write(1,*)' a2n=',a2n

write(1,*)'Средн.состав пленки - массовые доли: X1=',xxk1

write(1,*)' X2=',xxk2

write(1,*)'Средняя температура пленки: Tk=',Tk

write(1,*)' Свойства тройной смеси'

write(1,*)' ----------------------'

write(1,*)'Tемпература на линии насыщения: ',TSk

write(1,*)'Молярная масса в жидкой фазе: ',AMsk

write(1,*)'Молярная масса в паровoй фазе: ',AMpsk

write(1,*)'Плотность в жидкой фазе: ',rosk

write(1,*)'Плотность в паровoй фазе: ',ropsk

write(1,*)'Теплота фаз. перехода пар-жидк.: ',rsk

write(1,*)'Поверхостное натяжение: ',ssk

write(1,*)'Динамическая вязкость жидкости: ',amusk

write(1,*)' ---------------------'

write(1,*)'Температурный напор конденсации: dTk=',dTk

write(1,*)'Коэфф. теплоотдачи в канале конд.: ',alk

write(1,*)'Параметр ребра канала конденсации: ',amk

write(1,*)'КПД ребра: ',pdrk

write(1,*)'КПД пов. конд.: ',pdk

write(1,*)'Q конд (=Q исп): ',QQ

write(1,*)'Темп. напор в проставочном листе: dTst=',dTst

write(1,*)

write(1,*)' ИСПАРИТЕЛЬ'

write(1,*)' =========='

write(1,*)'Массовые доли жидкости на входе: Xi1=',xxi1

write(1,*)' Xi2=',xxi2

write(1,*)'Темп. смеси на входе в прогенератор: Ti=',Ti

write(1,*)'Re: ',Re

write(1,*)'Массовая доля пара на выходе: ',xxi

write(1,*)'Кратность циркуляции: ',ckr

write(1,*)'Относительный расход сливаемой жидкости: ',Go

write(1,*)'Объемн. доли на выходе испарителя: xi1v=',xi1v

write(1,*)' xi2v=',xi2v

write(1,*)'Объемн. доли на входе в испаритель: xi1n=',xi1n

write(1,*)' xi2n=',xi2n

write(1,*)'Среднее давление в испарителе: ',pis

write(1,*)'Средняя температура в испарителе: ',Tis

write(1,*)' Свойства тройной смеси'

write(1,*)' ----------------------'

write(1,*)'Температура на линии насыщения: ',TSi

write(1,*)'Молярная масса в жидкой фазе: ',AMsi

write(1,*)'Молярная масса в паровoй фазе: ',AMpsi

write(1,*)'Плотность в жидкой фазе: ',rosi

write(1,*)'Плотность в паровoй фазе: ',ropsi

write(1,*)'Теплота испарения: ',rsi

write(1,*)'Поверхостное натяжение: ',ssi

write(1,*)'Динамическая вязкость жидкости: ',amusi

write(1,*)' ---------------------'

write(1,*)'Массовый расход на входе: ',RM

write(1,*)'Модифицированное Re потока от пов-ти: ',Rem

write(1,*)'Температурный напор кипения: dTi=',dTi

write(1,*)'Коэффициент теплоотдачи в канале кипения:',ali

write(1,*)'Параметр ребра канала кипения: ',ami

write(1,*)'КПД ребра: ',pdri

write(1,*)'КПД поверхности кипения: ',pdi

write(1,*)'Эффективная плотность q в канале кипения:',qei

write(1,*)

write(1,*)' БАЛАНС ТЕМПЕРАТУР'

write(1,*)' ================='

write(1,*)'Суммарный перепад температур: dTs=',dTs

write(1,*)'Действит. эффект. плотность q конд.: qek=',qek

write(1,*)

write(1,*)' ГЕОМЕТРИЧЕСКИЕ ПАРАМЕТРЫ'

write(1,*)' ========================'

write(1,*)'Эффективная поверхность конденсации: ',Fke

write(1,*)'Число каналов конденсатора: ',nk

write(1,*)'Число каналов в 1 пакете: ',nk1

write(1,*)'Число пакетов: ',np

write(1,*)

write(1,*)' ГИДРАВЛИКА, ОПУСКНАЯ СИСТЕМА'

write(1,*)' ============================='

write(1,*)'Привед. скорость пара на вых. испарителя:',w0vp

write(1,*)'Привед. скорость смеси на там же: ',w0v

write(1,*)'Объемн. расходное паросодержание там же: ',beta

write(1,*)'Объемн. расход жидкости в опускной сист.:',Vop

write(1,*)'Средн. истиное объемн. паросодерж в исп.:',fis

write(1,*)'Скорость жидкости в опускной системе: ',wop

write(1,*)'Сечение опускной системы: ',fop

write(1,*)'Относит. уровень жидкости в опускн.сист.:',hop

write(1,*)'Движущее давление циркуляции: ',pdv

write(1,*)'Сопротивление на экономайзерном участке: ',dpe

write(1,*)'Сопротивление жидкости на входе в исп.: ',dp

write(1,*)'Интеграл psi: ',cks

write(1,*)'Сопротивление трения двухфазному потоку: ',dpd

write(1,*)'Истинное объемн. паросод. на вых. исп.: ',cfi

write(1,*)'Потери давления на ускорение: ',dpu

write(1,*)'Сопротивление на выходе из испарителя: ',dpv

write(1,*)'Сопротивление на столб жидкости: ',dps

write(1,*)'Гидравлическое сопротивл. подъемн. части:',dppod

return

end

================================================================

!Основная подпрограмма для вычисления температуры

subroutine temp(p,T,x1,x2)

x3=1-x1-x2

B=-364.65*x3-345.48*x2-302.82*x1+x3*x2*(13.44-7.4*(x2-x3)+2.*(x2-x3)**2)+x3*x1* &

(93.24-58.*(x1-x3)+27.*(x1-x3)**2)+x2*x1*(54.92-26.1*(x1-x2)+8.*(x1-x2)**2)+5.*x1*x2*x3

G=2.4024*x3+2.5086*x2+2.8919*x1+x3*x2*(.084-.0215*(x2-x3)+.0103*(x2-x3)**2)+x3*x1* &

(.2864-.157*(x1-x3)+.061*(x1-x3)**2)+x2*x1*(.194-.086*(x1-x2)+.041*(x1-x2)**2)+.004*x1*x2*x3

T=tmp(p,B,B,G,0.)

if(T<90.19)then

B1=B-9.3*x3; DG=.1028*x3

T=tmp(p,B,B1,G,DG)

if(T<87.29)then

B1=B1-7.3*x2; DG=DG+.0839*x2

T=tmp(p,B,B1,G,DG)

if(T>120)print*,'PROBLEM: T>120 !!!'; pause ''

endif; endif

end subroutine

!================================================================

!Подпрограмма функция для вычисления температуры

function tmp(p,B,B1,G,DG)

tmp=B1/(log10(98.07*pp)-G-DG+.01*B)

end function