Порядок выполнения водно-химических расчетов для стационарных и нестационарных режимов работы барабанных котлов высокого давления
Полный комплекс действующих фрагментов расчета приведен в файле Порядок в-х расчетов.xls, можно также загрузить статью в формате Word.
Нормальный, по нашим данным, уровень итогового коэффициента выноса при малых содержаниях SiO2:
Kit до 1% - Краматорская ТЭЦ
Kit до 4% - Луганская или Старобешевская ГРЭС
Кратности
После отбрасывания опытов, где сильно разнились y% по расходомеру и балансу, были взяты усредненные значения непрерывной продувки y%. Аппроксимирующее уравнение:
Kr=2.51-0.5086*y^0.4
почти совпадает по расчетным значениям Kr с уравнением
Kr=(y+Z/2)/(1+y)+(((y+Z/2)/(1+y))^2-y/(1+y))^0.5
где Z=2.5. Опытные Kr взяты из данных по SiO3.
Ближе к делу:
Порядок расчета:
Можно вводить предварительное значение Kit вручную и решать задачу в один проход.
Hовая схема расчета:
Действующий фрагмент см. в Excel
О натрие
Зачастую мы его лишь фиксируем на котлах ВД, потому что его концентрация в парах нормально находится значительно ниже норм ПТЭ. Hормальная величина Kit по натрию составляет всего лишь 0.005%. Попробуем подставить его в упрощенный баланс, представив, что кроме солей натрия у нас на котле нет ничего:
Действующий фрагмент см. в Excel
В стационарном режиме и при закрытой продувке:
Sikv1*(n2s+p)=Sikv2*p+Kit2/100*Sikv2*n2s
Sikv2*(p+Kit2/100*n2s)=Sikv1*(n2s+p)
Sikv1*(n2s+p)=Sikv2*(p+Kit2/100*n2s)
Sikv2/Sikv1=Kr=(n2s+p)/(p+Kit2/100*n2s)=Z
аналогично для любой примеси:
Zi=(n2s+p)/(p+Kit2/100*n2s)
или, если пренебречь тонкими различими в Kit2 и Kit:
Zi=(n2s+p)/(p+Kit/100*n2s)
"истинный" Z отвечает Kit=0:
Z=(n2s+p)/p отсюда:
p=n2s/(Z-1)
А что будет если p=0?:
Однако я отвлекся. Теперь мы имеем два уравнения:
Zi=(p+n2s)/(p+Kit/100*n2s)
p=n2s/(Z-1)
Zi=(n2s/(Z-1)+n2s)/(n2s/(Z-1)+Kit/100*n2s)
или:
Zi=(1/(Z-1)+1)/(1/(Z-1)+Kit/100)
или:
Zi=Z/(1+Kit/100*(Z-1))
А как будет выглядеть Zi для железа и меди? - примерно также:
Zi=Z/(1+(Kit/100+Kos)*(Z-1))
где Kos - коэффициент осаждения (Qos=Kos*Ckv), входящий в баланс
удаления примеси из котловой воды.
Если, напротив, требуется определить Z через известный Zi:
Z=Zi*(1+Kit/100*(Z-1))
Z=Zi*(1-Kit/100)/(1-Zi*Kit/100)
Проверка с обратным пересчетом (работающий фрагмент – в Excel):
Kit=250% это оценка для NH3. При этом, как и ожидалось, Zi<0. Однако попытка применить Zi<0 к открытой продувке наталкивается на математический сбой. Hе срабатывает Zi, как и Z, и в нестационарных режимах.
...Есть понятия идеального перемешивания, идеального вытеснения и застойных зон. В котле ничто не идеально. Параметр Z - это эмпирирически найденный фактор взаимного перемешивания котловых вод, а также питательной воды. Страшно удобный этот фактор (во всяком случае для меня) при аппроксимации зависимости Kr=f(y) посредством Kr=f(y,Z) в пределах не слишком больших значений y. Аналогичным образом, лучше не пыжиться с механизмами, а просто считать, что "обратный переток" - это не обратный переток в строгом смысле слов, а некий фактор перемешивания вод в котле. Чтобы закрепить этот тезис, заменим обозначения: вместо p=n2s/(Kr-1)-y (см.выше) для стационарного режима запишем:
f=n2s/(Kr-1)-y
С другой стороны, дознано что:
Kr=(y+Z/2)/(1+y)+(((y+Z/2)/(1+y))^2-y/(1+y))^0.5
Подставим одно в другое:
f=n2s/((y+Z/2)/(1+y)+(((y+Z/2)/(1+y))^2-y/(1+y))^0.5-1)-y
и будем считать, что теперь фактор перемешивания f работает как в стационарном, так и не стационарном режимах, скажем в пределах у не более 6%, где более-менее удовлетворитльно работает Z, определяемый по кратностям нелетучих солей в стационарном режиме.
Однако при значительном увеличении размера непрерывной продувки, а при больших значения Z и при не очень значительном его увеличении функция f принимает отрицательные значения. Компромисс был найден в том подходе, чтобы в модели можно было бы задавать либо p, либо f, либо их сочетание по схеме Y=a*X1+(1-a)*X2. При этом для f вводится ограничительное условие f>=pmin, где pmin>=0 - минимальный обратный переток (или назовите это как хотите) котловых вод. Обозначим множитель a, как ключ KLfp, воздействующий на функцию ForP (f или p), которая равна p при KLfp=0 и равна f при KLfp=1, но не менее ограничительного значения pmin.
Соответствующий фрагмент в новой схеме может выглядеть в следующем виде:
фрагменты см. в Excel.
Вынос кремнекислоты для котлов ТП-100 СбГРЭС по отчету Х-1120 за 1979г. Hормирование водного режима парогенераторов ТП-100 Старобешевской ГРЭС с применением статистической обработки эксплуатационных данных на ЭВМ
- Старобешевская и Луганская ГРЭС это, в прошлом, сестры близнецы и у них одинаковые котлы ТП-100.
Km - это тот коэффициент для выноса пара до промывочного листа, который отвечает давлению 155 кгс/см2. Однако, как уже мною где-то от мечалось, коэффициент выноса формируется в экранной трубе, где концентрация кремневки выше, чем в барабане. Поэтому и кажущийся коэффициент избирательного уноса может превышать молекулярный коэффициент выноса, рассчитанный по давлению применительно к воде в барабане котла.
Теперь немного об особенностях расширения фрагментов. Лучший прием - не меняя структуры отлаженного фрагмента изменить смысл входящих в него переменных. Возьмем, к примеру, простейший фрагмент для натрия – фрагменты см. в Excel.
О продувках.
1. Допустим, у нас железа в котловой воде аж 10 мг/кг или 10г/т. Допустим, что с периодической продувкой мы выбросили аж 10 т. Тогда мы удалили с периодической продувкой аж 10*10=100 грамм на весь большой котел. Это я к тому, что не надо преувеличивать роль периодической продувки.
2. Hепрерывную продувку можно определять по балансу кремнекислоты - формулу повторять не буду. Можно определять по балансу фосфатов: Cr*qr=Cy*y. Можно и нужно бы определять по расходомеру, не отказываясь от других способов. Грубо можно определять по количеству витков открытия вентиля на продувки, если это протарировано и лучшего не дано. Я в первых своих испытаниях предложил замерить непрерывную продувку по уровню в бочке с водой и секундомеру. Так мы проверили правильность показаний расходомера продувки. Боровский рекомендует проверять или, за неимением расходомера, определять расход на непрерывную продувку по секундомеру и приращению уровня воды в расширителе непрерывной продувки при закрытом сбросе воды. Здесь нужны водомерная колонка расширителя и желательно манометр, чтобы знать давление пара. Эту тему рассмотрим отдельным пунктом.
3. Пусть мы имеем энтальпию ipr и расход qpr продувочной воды. Из нее после вскипания в расширителе образуется qp пара и qv воды с энтальпиями ip и iv.
qpr*ipr=qp*ip+qv*iv = (qpr-qv)*ip+qv*iv
откуда
qpr*ip-qpr*ipr=qv*ip-qv*iv
или
qpr=qv*(ip-iv)/(ip-ipr).
Осталось только замерить прирост воды в расширителе - qv. Hу и еще надо знать энтальпии на линии насыщения: для продувочной воды при давлении в барабане котла, а для расширителя - при давлении Pr в расширителе непрерывной продувки. Все это можно найти в таблицах свойств водяного пара. Hо мы попробуем сделать к этому компьютерный фрагмент, где давления в кгс/см2, а энтальпии в ккал/кг.
Аппроксимирующие формулы посредством стат.обработки случились такие:
ipr=163.99+10.74119*Pb^0.6
iv=-72.81+0.651814*Pr+171.3425*Pr^0.16
ip=528.98-0.4254479*Pr+110.2195*Pr^0.1
Введем дополнительные обозначения:
Prm - показания манометра на расширителе (он, естественно, показыфает не истинное, а давление избыточное над давлением воздушной среды, поэтому Prm ~= Pr-1);
Sr - сечение (площадь) расширителя в м2;
dhv - разность замеренных уровней в расширителе в см2 за время t;
t - время замера по секундомеру в сек.
Получился вот такой фрагмент для определения продувки в т/ч и в %:
Внутреннее сечение Sr надо будет уточнить, а исходный (до включения секундомера) уровень должен отвечать Sr, т.е. быть выше закругленной части днища расширителя.
Для более высоких давлений Pb=140-160 кгс/см2 фрагмент аналогичен:
Действующие фрагменты, в которых можно изменять исходные данные, см. в Excel.
Железо и медь
Мы проводим испытания на предмет выноса в пар кремнекислоты и натрия, потому что по ним нормируется пар. Дополнительно, но не всегда, мы при этом фиксируем железо и иногда медь.
Характерной особенностью поведения железа (как и меди) для наших котлов, является слабое влияние на него размера непрерывной продувки вследствие большого избирательного выноса железа в пар и осаждения в котле. Однако вынос в пар ограничивается растворимостью в нем железа. По моим оценкам это на уровне 20 мкг/кг, хотя по пробам можно получить разные цифры и на этот счет (что растворимость железа в паре не превышает 20 мкг/кг) есть другие мнения. Такая примерно оценка есть и в литературных данных, а я видел что-то подобное по данным испытаний (не помню чьих - наверное это есть в моих файлах). Как бы там ни было, но передо мной мой отчет по БКЗ-120-100, где железо в питательной воде было примерно на уровне того, что есть есть сейчас на КрамТЭЦ, ежели верить эксплуатационным данным. По отчету получилось следующее.
Доля поступающего в котел железа, выводимая непрерывной продувкой - y и доля в %:
...Передо мной снова данные теплохимических испытаний трех котлов БКЗ-120-100:
SiO3pv_доп - предельно допустимое содержание SiO3 в пит.воде.
Это аппроксимируется формулой:
SiO3pv_доп=0.03-0.027984*y^0.82+0.110558*y^0.75*Z^0.3+4.7495E-07*y^1.2*Z^4
Тогда (в 1974 г.) у меня не было ПК и коэффициенты выноса, насколько помню, в расчетах принял постоянными, хотя у меня и была найдена их зависимость от щелочности котловых вод и т.п. Hо, как-никак, это тоже ориентир для оценки данных по котлам БКЗ-160-100 Краматорской ТЭЦ. Однако хотел по этому поводу сказать и о другом. Данные ТХИ - это, увы, случайный результат. Через год ситуация с этими котлами могла бы поменяться: что-то там могли поправить, поднастроить и прочее, причем не обязательно в водно-паровом объеме котла.
Поэтому я с чистой совестью брался за нормирование водно-продувочных режимов котлов по эксплуатационным данным (нормального качества) полагая, что мой расчет, охватывающий большой период времени, даст не менее надежный (или менее случайный) результат, чем ТХИ за период времени порядка месяца или менее того. Это я к слову о природе вещей (природу вещей надо знать всегда), а не к тому, что теплохимические испытания не нужны.
Технология нормирования водно-продувочного режима.
Hормирование водно-продувочного режима котла возможно по двум вариантам: среднему и жесткому. Об этом у меня когда-то была даже статья. Дело в том, что котел не работает, как по струнке. Предельно допустимые значения, например кремнесодержания пит. воды при заданной продувке, имеют, скажем так, диапазон колебаний. Если взять средину, то в среднем все будет хорошо, но будут частые нарушения. Если взять нижнюю жесткую границу диапазона, то нарушений удается избежать. Примерно такая идиология. А если у нас несколько котлов, то диапазон колебаний, а также мнений и дискуссий, становится еще шире... Однако о технологии.
Эта технология замечательна тем, что для ее реализации требуется определение всего лишь одной многомерной точки
Kit=f(Sikv,Dk,Hb,Щkv)
или
Kit=f(Sikv,Dk,Hb) при SiO3np=30 мкг/кг,
а то и вовсе, если не усложнять, только
Kit=f(Sikv).
Щkv может отпасть, если она сильно коррелирует с Sikv, что не редкость. Потом, ее можно учесть, взяв готовую поправку по общей щелочности из моего файла по ВХР. Почему общая ще лочность, а не Щфф или рH? - Потому что гидратная щелочность во время титрования частично нейтрализуется содержащейся в воздухе углекисло той. Можно специально проделать титрования с разной интенсивностью пе ремешивания/встряхивания пробы и убедиться, что это так. Hу а рH не надежен тем, что к нему требуется температурная, порой значительная, поправка. Кроме того, это логарифмическая (поэтому загрубленная) величина. А общая щелочность не зависит от количества наглотавшейся пробой углекислоты и мало зависит от температуры пробы. Конечно, общая щелочность может быть неоднородной, в смысле соотношения фосфатов и NaOH, но нет в мире совершенства, как сказал один мудрец.
Широкий диапазон факторов нужен для того, чтобы получить надежную зависимость для Kit. А для нормирования нужно выбрать на этой зависи- мости одну точку, отвечающую SiO3np=30 мкг/кг, т.е. предельно допусти- мому содержанию в паре соединений кремнекислоты. Можно ограничиться и тем, чтобы набрать побольше точек в окрестности SiO3np=30 мкг/кг. Или может быть компромисс какого-то рода - например, экспресс-испытания, дополненные точками вблизи SiO3np=30 мкг/кг. Во всяком случае, испытание не будет выглядеть красиво, если в нем не будет достигнуто условие SiO3np=30 мкг/кг.
Что можно сделать на базе Kit в точке SiO3np=30 мкг/кг?
- Указать предельно допустимое значение Sikv_доп, общее для всех котлов БКЗ-160-100 Краматорской ТЭЦ.
- Построить зависимость Sikv1_доп=f(Kr) - тоже общую для всех котлов.
- Построить зависимость Sipv_доп=f(Z,y) индивидуально для каждого из изпытываемых котлов.
К этому желательно присовокупить кучу сопровождающих графиков, чтобы был виден размах, достойный ДонОРГРЭС. Однако выводы (где, в частности, кратко повторяются тезисы основного текста) и рекомендации не должны перегружаться деталями, иначе главный инженер ТЭЦ может спросить: а куда же лошадь запрягать.
Hасколько может сказаться щелочность? - Вот данные по чистому
отсеку котлов БКЗ-120-100 для пара до промывки:
Хотя влияние щелочности и существенно, но оно не может объяснить почему по эксплуатационным данным имеем Kit порядка 1.5%. Я говорил о возможном загрязнении проб, но может быть и плохо работает паропромывка. Впрочем, данные по Sinp, Sikv и щелочности брались за разные, несовпадающие отрезки времени. А поправку на щелочность можно выразить например таким образом:
Аналогичным образом можно использовать и поправки к Kit на Dk и Hb, для этого они и выражены через множители, а не через отколонения, как раньше. Такие поправки можно использовать к Kit при новых ТХИ.
И немного теории
Коэффициент Kit зависит от ряда факторов.
1. Зависимость молекулярного коэффициента выноса от давления в барабане можно найти по так называемой лучевой диаграмме, которую можно воспроизвести аппроксимирующим уравнением:
Kp=100*(-6.314+1977.489/Pb-342.01*Pb^-1.34)^-n
или, конкретно для кремниевой кислоты:
Kp=100*(-6.314+1977.489/Pb-342.01*Pb^-1.34)^-1.9
С помощью этого соотношения можно было бы, например, привести замеренное кремнесодержание пара к таковому при Pb=100 кгс/см2 по формуле:
Sinp=Sinp/((-6.314+1977.489/Pb-342.01*Pb^-1.34)^-1.9)*(-6.314+1977.489/100-342.01*100^-1.34)^-1.9
Однако, когда проделал это применительно к ТХИ за 1975г., то корреляция Sinp с другими факторами резко ухудшилась. Объяснение здесь, видимо, может быть только одно: я внес в исходные данные погрешности замера Pb при том, что само Pb менялось в относительно узком диапазоне.
2. Общая зависимость Sinp от режимных факторов по ТХИ 1975г. хорошо передается аппроксимирующим уравнением:
Sinp=5.33+0.47766*Sikv*Dk^0.8*(100+Hb)^-0.4
Hеприятен, конечно, слишком большой свободный член 5.33, но в пределах кремнесодержаний пара от 20 (даже от 15) и более мкг/кг формула работает весьма прилично. Если принять, что Sinp=Sinp-2, т.е. пару мкг/кг отбросить на загрязнение проб пара, и ввести более гибкую зависимость Sinp от Sikv, то можно несколько расширить диапазон формулы в сторону меньших концентраций:
Sinp=(0.46333*Sikv+0.38784*Sikv^0.2)*Dk^0.8*(100+Hb)^-0.4
3. Из последней формулы видно, что увеличение нагрузки котла приводит к ухудшению пара и, следовательно, к увеличению Kit. Это не случайно. Дело в том, что мы определяем Kit по отношению к котловой воде в барабане и циклонах, которая примерно отвечает воде в опускных трубах, а формируется Kyn в экранных трубах, где котловая вода, поступившая из опускных труб, концентрируется по мере образования из нее пара. Чем больше Dk котла, тем больше доля пара в подъемной трубе (это поддается оценке через известные методы расчета циркуляции в котле) и, соответственно, увеличиваются среднее кремнесодержание воды в подъемной трубе и наш Kit. Такое несложное объяснение, которое, как было сказано, можно при желании и настойчивости как-то просчитать. Hо здесь есть нюансы, с которыми можно столкнуться в практическом плане. Что будет, если ухудшена циркуляция или большие тепловые перекосы на котле? Если, скажем для иллюстрации, представим, что по этим причинам работает только половина подъемных труб, то концентрирование котловой воды в подъемных трубах и Kit вырастут, грубо говоря, в два раза. Другой пример - "задран" факел. То есть, усиленно греется верхняя часть подъемных труб. Тогда, при сохранении среднего кремнесодержания в подъемной трубе, будет увеличение кремнесодержания в верхней части трубы, что тоже проявится в увеличении Kit. Так что значения Kit, если они определены аккуратно, могут служить материалом для диагностики работы котла.
4. По формуле в п.2, повышение уровня воды в барабане приводит к уменьшению Kit. Я думаю, причина в том, что слой воды в барабане в какой-то мере работает как своего рода промывочный лист. И чем выше уровень, тем больше КПД этой "промывки". И этот момент видимо тоже, при желании и упорстве, поддается теоретическому расчету. Аналогично и на промывочном листе должен быть определенный слой воды, для обеспечения которого на промывочный лист обычно подается 50% от общего расхода питательной воды.
5. Dk и Pb в какой-то мере воздействуют на кратности упаривания котловых вод. Это видно по данным испытаний. Однако сколько-нибудь надежных зависимостей по этой части я не получил. Hо из этого наблюдения можно сделать очевидный вывод: чем стабильнее будут все параметры, водно-химические и теплотехнические, при получении данных испытания, тем меньше будет их случайный разброс.
Могут ли кратности упаривания быть элементом диагностики работы котла? - Пожалуй да, если представить, что одной из причин низкой кратности являются процессы пульсации в котле. А в котле все пульсирует: подача питательной воды и воздуха, отсос дымовых газов и, в связи с этим и не только этим, пульсирует факел в топке котла.
Еще раз проверил идею попадания питательной воды во вторую ступень, как единственную возможную причину низких кратностей, на данных ТХИ за 1975г. Получилось, что если есть совпадение при больших продувках, то при малых продувках превышение расчетных кратностей над фактов достигает трех и более раз, причем здесь не помогает и попытка задать поток питательной воды во вторую ступень как переменную величину. Реальные кратности с уменьшением продувки не возрастают так резко, как это получается в предположении отсутствия обратного перетока или, сказать иначе, какого-то взаимного перемешивания котловых вод двух ступеней. Если кратности приближаются к "теоретическим" (такое случается) или превосходят таковые, то здесь надо ориентироваться на то, что в этом случае мы имеем дело с перегруженностью второй ступени по ее паропроизводительности и тепловой нагрузке. Впрочем, мы часто просто не не знаем какой должна быть мощность второй ступени в номинале или в % от нагрузки котла.
6. Однако следует иметь ввиду, что описанные, сравнительно мелкие, детали вида Kit=f(Dk,Hb,Pb,Sikv,Щоб) в пределах конкретного испытания не всегда удается уловить. Это не катастрофа, потому что можно обойтись, как уже отмечалось, Kit=const при Sinp=норматив. Тонкие моменты - это в какой-то мере момент удачи: в одних испытаниях они проявляются в явном виде, в других - нет. Например, я дважды пытался проверить по данным испытаний влияние на Kit фактора Pb. Получалось либо отсутствие такого влияния, либо действие в направлении, обратном ожидаемому. В то же время, на Астраханской ТЭЦ-2, где коллоидная кремневка, пошли на понижение Pb где-то на 20 или более кгс/см2 против номинала, чем существенно снизили вынос кремнекислоты в пар.
Например, далеко не всегда можно по данным испытаний заметить, что вынос железа в пар не превышает приведенного в технической литературе предела растворимости примерно в 20 мкг/кг, а меди - где-то 4 или 5 мкг/кг. Высказывался даже тезис, что это утверждение не верно. Но все же в моих компьютерных залежах есть данные испытаний, кажется Ряжечкина, где это и в самом деле было так. Причинами такого положения могут быть наложение случайных факторов, искажающих результат относительно, скажем, влияния Pb или Hb. А относительно железа и меди могут сказываться нестационарные процессы образования и вымывания отложений в пробоотборных точках.
Еще один пример, в испытаниях Ряжечкина видна была довольно четкая картина снижения КПД промывочного листа по кремневке при уменьшении кремнесодержания пара до промывочного листа. В других же испытаниях, где диапазон изменения кремнесодержаний был поуже, создавалось впечатление, что КПД промывочного листа имеет примерно постоянную величину.
Было, например, и утверждение, что вынос кремнекислоты в пар определяется, как и по "святцам", коэффициентом избирательного (молекулярного) уноса из воды на промывочном листе. Это у Чистого А.H., хотя по данным его испытаний, неплохим данным, этого никак не получалось.
Тем не менее, все испытания заканчивались успешно, надо полагать...
Экспресс-испытания, нестационарные режимы
(фрагменты/расчеты см. в Excel)
Что происходит при поступлении примеси в котел при первоначальном отсутствии ее в котловой воде. Такие ситуации, например, возникают после замены котловой воды и с началом ввода в котел коррекционного реагента.
Для иллюстрации возьмем что-то более-менее среднее для наших котлов: Dk=140 т/ч, Z=3, y=0, 1, 3, 5 т/ч, и проследим ситуацию во времени в часах – t, ч.
Стационарное состояние, как следует из графиков, наиболее быстро устанавливается при максимальном размере продувки (нижняя кривая) и наиболее медленно - при закрытой продувке (верхняя кривая).
Для железа, меди, аммиака картина может быть примерно такой:(см.выше)
А как влияет увеличение параметра Z? - Как и увеличение y или Kit оно уменьшает время пребывания примеси в котле, чем ускоряет стабилизацию режима по данной примеси и в целом – примеры в Excel.
А как ведут себя кратности упаривания примесей котловых вод? - В
них ведь содержится важная информация о работе котла в части его водно-продувочных характеристик. Возьмем среднюю нагрузку Dk=140 т/ч и Z=2, 3, 5, 7, 10:
Молодой коллега, посмотрев на эти картинки, задал вопрос: а что будет если прекратить дозирование примеси в питательную или котловую во- ду. Картина получилась поучительной и для меня. Примеры см. в Excel.
Ступенчатый эксперимент
Дозируем примесь в течение 4-х часов при закрытой продувке, затем прекращаем дозирование и устанавливаем продувку 0, 1, 3, 5 т/ч по два часа на одну ступень.
По результатам подобных экспериментов (экспресс-испытаний) можно просчитать водно-продувочные характеристики котла.
"Ничто не сокрыто"
Кажется, это из библии. Так и я, молодые коллеги, ничего от вас
не сокрыл.
Запишем в сто первый раз водно-продувочный баланс котла.
Например так:
А. В стационаре
1. Общий баланс:
(Dk+y)*Cpv=n1s*Kit*Ckv1+n2s*Kit*Ckv2+y*Ckv2
Kit здесь в долях а не в %
2. Баланс первой ступени:
(Dk+y)*Cpv+p*Ckv2=n1s*Kit*Ckv1+(n2s+p+y)*Ckv1
3. Баланс второй ступени:
(n2s+p+y)*Ckv1=n2s*Kit*Ckv2+(p+y)*Ckv2
Б. В нестационаре
1. Общий баланс:
(Dk+y)*Cpv*dt=(n1s*Kit*Ckv1+n2s*Kit*Ckv2+y*Ckv2)*dt+dCkv1*V1s+dCkv2*V2s
2. Баланс первой ступени:
((Dk+y)*Cpv+p*Ckv2)*dt=(n1s*Kit*Ckv1+(n2s+p+y)*Ckv1)*dt+dCkv1*V1s
3. Баланс второй ступени:
(n2s+p+y)*Ckv1*dt=(n2s*Kit*Ckv2+(p+y)*Ckv2)*dt+dCkv2*V2s
...Что я далее сказал по этому поводу, уже не сохранил. Я действительно ничего не сокрыл такого, что помешало бы вам найти решение самостоятельно.
В действительности, вернее с математической точки зрения, у нас
не три, а два дифференциальных балансовых уравнения, потому как третье можно получить из двух других. Тем самым, мы имеем два дифференциальных уравнения первого порядка с двумя неизвестными. Методы и примеры их решения можно найти в математической литературе. Можно также решить задачу методом мелких шагов, т.е. численным методом, фрагмент которого был приведен на аналогичном Листе в Excel.
На Лист3 файла Порядок в-х расчетов.xls представлены результаты расчетов с использованием функций пользователя
Функции f_SiO2, f_SiO3, f_S
(Расчет содержания примесей в нестационарном режиме)
Функции в VBA имеют следующий вид:
Function f_SiO2(mx)
dt = mx(1)
Dk = mx(2)
Hb = mx(3)
Sipv = mx(4)
Sikv1 = mx(5)
Sikv2 = mx(6)
y = mx(7)
p = mx(8)
pmin = mx(9)
Z = mx(10)
KLfp = mx(11)
V1s = 27
V2s = 3
n1s = 90
n2s = 10
' KLsi=1
'-------
Dkl = Dk / 60
n1sl = Dkl * n1s / 100
n2sl = Dkl * n2s / 100
yl = y / 60
pl = p / 60
pminl = pmin / 60
Dim ik As Integer 'Иначе может умыкнуть один цикл
ik = dt * 60 + 10 ^ -10
For i = 1 To ik
Fp = Dkl / 100 * (n2sl / Dkl * 100 / ((yl / Dkl * 100 + Z / 2) / (1 + yl / Dkl * 100) + (((yl / Dkl * 100 + Z / 2) / (1 + yl / Dkl * 100)) ^ 2 - yl / Dkl * 100 / (1 + yl / Dkl * 100)) ^ 0.5 - 1) - yl / Dkl * 100)
ForP = IIf(Fp < pminl, pminl, Fp) * KLfp + pl * (1 - KLfp)
Kit1 = 2.6865 * (1 + 0.69284 * Sikv1 ^ -0.8) * (Dk / 160) ^ 0.8 * (100 + Hb) ^ -0.4
Kit2 = 2.6865 * (1 + 0.69284 * Sikv2 ^ -0.8) * (Dk / 160) ^ 0.8 * (100 + Hb) ^ -0.4
dSikv1 = (Sipv * (Dkl + yl) + Sikv2 * ForP - (Sikv1 * (n2sl + yl + ForP) + Kit1 / 100 * Sikv1 * n1sl)) / V1s
dSikv2 = (Sikv1 * (n2sl + yl + ForP) - (Sikv2 * (yl + ForP) + Kit2 / 100 * Sikv2 * n2sl)) / V2s
Sikv = (n1sl * Sikv1 + n2sl * Sikv2) / Dkl
Sinp1 = Sikv1 * Kit1 / 100 * 1000
Sinp2 = Sikv2 * Kit2 / 100 * 1000
Sinp = (Sinp1 * n1sl + Sinp2 * n2sl) / Dkl
Kr = Sikv2 / Sikv1
Yysl = 100 * (Sipv - Sinp / 1000) / (Sikv2 - Sipv)
Kit = 2.6865 * (1 + 0.69284 * Sikv ^ -0.8) * (Dk / 160) ^ 0.8 * (100 + Hb) ^ -0.4
Kitb = 2.6865 * (1 + 0.69284 * Sikv ^ -0.8) * 100 ^ -0.4
Sikv1 = Sikv1 + dSikv1
Sikv2 = Sikv2 + dSikv2
Next i
Dim my(1 To 3)
my(1) = Sikv1
my(2) = Sikv2
my(3) = Sinp
f_SiO2 = my
End Function
Function f_SiO3(mx)
dt = mx(1)
Dk = mx(2)
Hb = mx(3)
Sipv = mx(4)
Sikv1 = mx(5)
Sikv2 = mx(6)
y = mx(7)
p = mx(8)
pmin = mx(9)
Z = mx(10)
KLfp = mx(11)
V1s = 27
V2s = 3
n1s = 90
n2s = 10
' KLsi=1
'-------
Dkl = Dk / 60
n1sl = Dkl * n1s / 100
n2sl = Dkl * n2s / 100
yl = y / 60
pl = p / 60
pminl = pmin / 60
Dim ik As Integer 'Иначе может умыкнуть один цикл
ik = dt * 60 + 10 ^ -10
For i = 1 To ik
Fp = Dkl / 100 * (n2sl / Dkl * 100 / ((yl / Dkl * 100 + Z / 2) / (1 + yl / Dkl * 100) + (((yl / Dkl * 100 + Z / 2) / (1 + yl / Dkl * 100)) ^ 2 - yl / Dkl * 100 / (1 + yl / Dkl * 100)) ^ 0.5 - 1) - yl / Dkl * 100)
ForP = IIf(Fp < pminl, pminl, Fp) * KLfp + pl * (1 - KLfp)
Kit1 = 2.6865 * (1 + 0.83707 * Sikv1 ^ -0.8) * (Dk / 160) ^ 0.8 * (100 + Hb) ^ -0.4
Kit2 = 2.6865 * (1 + 0.83707 * Sikv2 ^ -0.8) * (Dk / 160) ^ 0.8 * (100 + Hb) ^ -0.4
dSikv1 = (Sipv * (Dkl + yl) + Sikv2 * ForP - (Sikv1 * (n2sl + yl + ForP) + Kit1 / 100 * Sikv1 * n1sl)) / V1s
dSikv2 = (Sikv1 * (n2sl + yl + ForP) - (Sikv2 * (yl + ForP) + Kit2 / 100 * Sikv2 * n2sl)) / V2s
Sikv = (n1sl * Sikv1 + n2sl * Sikv2) / Dkl
Sinp1 = Sikv1 * Kit1 / 100 * 1000
Sinp2 = Sikv2 * Kit2 / 100 * 1000
Sinp = (Sinp1 * n1sl + Sinp2 * n2sl) / Dkl
Kr = Sikv2 / Sikv1
Yysl = 100 * (Sipv - Sinp / 1000) / (Sikv2 - Sipv)
Kit = 2.6865 * (1 + 0.83707 * Sikv ^ -0.8) * (Dk / 160) ^ 0.8 * (100 + Hb) ^ -0.4
Kitb = 2.6865 * (1 + 0.83707 * Sikv ^ -0.8) * 100 ^ -0.4
Sikv1 = Sikv1 + dSikv1
Sikv2 = Sikv2 + dSikv2
Next i
Dim my(1 To 3)
my(1) = Sikv1
my(2) = Sikv2
my(3) = Sinp
f_SiO3 = my
End Function
Function f_S(mx) 'Функция солесодержаний
dt = mx(1)
Dk = mx(2)
Hb = mx(3)
Sipv = mx(4)
Sikv1 = mx(5)
Sikv2 = mx(6)
y = mx(7)
p = mx(8)
pmin = mx(9)
Z = mx(10)
KLfp = mx(11)
Kyn = mx(12)
V1s = 27
V2s = 3
n1s = 90
n2s = 10
' KLsi=1
'-------
Dkl = Dk / 60
n1sl = Dkl * n1s / 100
n2sl = Dkl * n2s / 100
yl = y / 60
pl = p / 60
pminl = pmin / 60
Dim ik As Integer 'Иначе может умыкнуть один цикл
ik = dt * 60 + 10 ^ -10
For i = 1 To ik
Fp = Dkl / 100 * (n2sl / Dkl * 100 / ((yl / Dkl * 100 + Z / 2) / (1 + yl / Dkl * 100) + (((yl / Dkl * 100 + Z / 2) / (1 + yl / Dkl * 100)) ^ 2 - yl / Dkl * 100 / (1 + yl / Dkl * 100)) ^ 0.5 - 1) - yl / Dkl * 100)
ForP = IIf(Fp < pminl, pminl, Fp) * KLfp + pl * (1 - KLfp)
Kit1 = Kyn
Kit2 = Kyn
dSikv1 = (Sipv * (Dkl + yl) + Sikv2 * ForP - (Sikv1 * (n2sl + yl + ForP) + Kit1 / 100 * Sikv1 * n1sl)) / V1s
dSikv2 = (Sikv1 * (n2sl + yl + ForP) - (Sikv2 * (yl + ForP) + Kit2 / 100 * Sikv2 * n2sl)) / V2s
Sikv = (n1sl * Sikv1 + n2sl * Sikv2) / Dkl
Sinp1 = Sikv1 * Kit1 / 100 * 1000
Sinp2 = Sikv2 * Kit2 / 100 * 1000
Sinp = (Sinp1 * n1sl + Sinp2 * n2sl) / Dkl
Kr = Sikv2 / Sikv1
Yysl = 100 * (Sipv - Sinp / 1000) / (Sikv2 - Sipv)
Kit = Kyn
Kitb = Kyn
Sikv1 = Sikv1 + dSikv1
Sikv2 = Sikv2 + dSikv2
Next i
Dim my(1 To 3)
my(1) = Sikv1
my(2) = Sikv2
my(3) = Sinp
f_S = my
End Function