Разработка численной схемы стекловаренной печи

сти продукции, рост рентабельности производства и накоплений предприятия.

Более полное использование основных производственных фондов приводит также к сокращению необходимости ввода новых производств ен-ных мощностей при изменении объема производства, а следовательно, к рациональному использованию финансовых ресурсов предприятия, например, на увеличение доли отчислений от прибыли в фонд потребления, направление большей части фонда накопления на механизацию и автоматизацию технологических процессов и другие инновации [5]. Кроме того, ускоряется оборачиваемость основных производственных фондов, что в значительной мере способствует решению проблемы сокращения разрыва в сроках физического и морального износа, ускорения темпов обновления основных фондов. Наконец, эффективное использование основных фондов тесно связано и с другой ключевой задачей современного периода экономической реформы - повышением качества выпускаемой продукции, поскольку в условиях рыночной конкуренции высококачественная продукция пользуется наибольшим спросом и быстрее реализуется.

Применение системы Global-EAM показывает, что успешное функционирование основных производственных фондов зависит от того, насколько полно реализуются экстенсивные и интенсивные факторы улучшения их использования. Экстенсивное улучшение использования основных фондов предполагает, что, с одной стороны, увеличится время работы действующего оборудования в календарный период, а с другой - повысится удельный вес действующего оборудования в составе всего оборудования, имеющегося на предприятии.

Интенсивное улучшение использования основных фондов предполагает повышение степени загрузки оборудования в единицу времени. Более высокая загрузка технологического оборудования может быть достигнута при модернизации действующих машин и механизмов, установлении оптимального режима их работы. Оптимальный режим технологического процесса обеспечивает увеличение выпуска продукции без изменения состава основных производственных фондов, без роста численности работающих и при снижении расхода материальных ресурсов на единицу продукции.

Важно отметить, что процесс использования моделей управления формированием и использованием основных производственных фондов с развитием рыночных отношений и обострением конкуренции будет постоянно совершенствоваться и нарастать.

Связано это с тем, что без моделирования процессов руководители не смогут правильно оценить все возрастающие потоки информации и, следовательно, своевременно и качественно разрабатывать и принимать управленческие решения по вопросам управления формированием и использованием основных производственных фондов предприятия.

Литература

1. Сафронов Н. Экономика организации (предприятия). М.: ИНФРА-М, 2005. 618 с.

2. URL: www.asutp.ru (дата обращения: 05.07.2009).

3. URL: www.bitec.ru (дата обращения: 05.07.2009).

4. URL: http://global-eam.ru/index.php?id=211 (дата обращения: 12.07.2009).

5. Разинькова О.П. Комплексный анализ эффективности формирования и использования основных фондов предприятия / Тверь: Изд-во ТГТУ, 2008.

УДК 666.1.031.2

РАЗРАБОТКА ЧИСЛЕННОЙ СХЕМЫ СТЕКЛОВАРЕННОЙ ПЕЧИ

В.С. Швыдкий, д.т.н.; С.Е. Собянин

В статье приведен пример построения дискретного аналога трехмерной математической модели гидродинамики и теплообмена ванны стекловаренной печи. Выполнены необходимые преобразования уравнений неразрывности и балансов импульса, проведен переход к уравнениям дискретного аналога на разнесенной конечноразностной сетке. Также даны рекомендации по решению уравнений дискретного аналога с использованием итерационного метода и по нахождению граничных условий.

Ключевые слова: математическая модель, стекловаренная печь, численная схема, дискретный аналог, аппроксимация.

При разработке численной схемы стекловаренной печи необходимо обязательно учесть осо-

бенности стекломассы: большую плотность и вязкость, медленное движение (практически ползу-

щее). Учитывая их, можно уверенно отнести стекломассу к несжимаемой жидкости.

Конечно, плотность стекломассы зависит от температуры, а температура в бассейне изменяется, но, во-первых, эта зависимость слабая и, во-вторых, температура стекломассы в ванне изменяется от 1000 до 1580 оС, что приводит к относительно небольшим изменениям плотности. Иными словами, с достаточной точностью можно считать для стекломассы div v=0 и записать уравнения движения в виде:

» в проекции на ось х

г ди ди диЛ

и--+ V--+ w—

дх дУ дг

= -дР + 2-д

дх дх

Ц

+

д_

дУ

Ц

ди дv

— + —

дУ дх

+

д дг

Ц

ди дм

— + —

дг дх

да

дх

+

в проекции на ось у

дч дч дч

и--+ V--+ w —

дх ду дг

др д

= —- +—

ду дх

дч ди

— + —)

дх ду

+

+2-д ду

с я. л дч

ду

+

д_ дг

г дч дwЛ

— + —

дг ду

в проекции на ось z

дм дм дм

и--+ V--+ м —

дх ду дг

др д

= —- +—

дг дх

дм ди

— + —

дх дг

+

+

д ду

г дм дvл

— + —

ду дг

+

дг

дг

Р£.

Вязкость стекломассы существенно зависит от температуры. Однако уравнения (1)-(3) - это уравнения баланса импульса для бесконечно малого объема. При составлении этого баланса можно считать, что температура стекломассы, следовательно, и ее вязкость, в пределах элементарного объема одинаковы. Тогда ц можно вынести из-под знака производных, что приводит к упрощению уравнений.

Например, в уравнении (1) выделим слагаемые

д2и

д2 V

ц—2 +Ц-

дх дудх

+ Ц

д2

м

дх

дzдx л

ди дv дм

— + — +—

дх ду дг

= 0,

то же сделаем и для уравнений (2) и (3). Заметим, что такой подход не исключает учета зависимости

ц от температуры, поэтому окончательный вид уравнений движения стекломассы можно представить следующим образом:

Р

ди ди ди

и--+ V--+ w—

дх ду дг

др д

= —- +—

дх дх

Г Я. Л

ди

V дх У

+

ду

ди

ду

+

дг

дv ду дч

и--+ у--+ w—

дх ду дг

др д

= —- +—

ду дх

ду

V дх У

+

ду

дУ

+

дг

ди

V дг J

ду

V дг

дм дм дм

и--+ V--+ м—

дх ду дг

др д

— + —

дг дх

с

Ц

дм дх

+

+

д_

ду

Ц

дм

дУ

+

дг

дм дг

Уравнение неразрывности можно упростить. Хотя формально необходимо учитывать наличие внутренних источников (стоков) массы, фактически сделать это нельзя из-за неопределенности кинетики физико-химических превращений и отсутствия математического описания. Поэтому используем уравнение неразрывности в форме

— (ри ) + — (ру) + — (р™ ) = 0.

Здесь возникло некоторое противоречие со сделанным утверждением о несжимаемости стекломассы, однако, во-первых, введение плотности под знаки производных не слишком усложняет уравнение, а, во-вторых, это важно для дальнейшего учета теплообмена.

Построение дискретного аналога всегда нужно выполнять для безразмерных уравнений, поскольку в этом случае легче оценить порядок аппроксимации и устойчивость численной схемы. Введем соответствующие безразмерные параметры. Предварительно, однако, приведем уравнения движения к дивергентному виду, объединив с уравнением неразрывности:

др д

— + —

дх дх

рии —ц

ди

дх

+

+

д ду

риу — ц

ди

ду

+

др д

— +—

ду дх

руи — ц

ду дх

дz

+

pиw — ц

ди дг

= 0;

+

д_ дУ

руу—ц

ду дУ

д

+—

дг

ру^ —ц

ду дг

= 0;

д

д

Р

д

д

Р

Р

Р

д

Р

д

д

др д

— +—

дz дх

pwu — ц

дw

дХ

+

+

ду

pwv—ц

дw

д

+—

pww — ц

дw

~дк

+pg = 0.

Отсчитаем компоненты скорости в долях среднерасходной скорости стекломассы в протоке ^)=Р/(86,4р0^2Ь2) м/с, где ро - масштабное значение плотности стекломассы (кг/м3); Р - производительность печи, т/сутки; z2 - высота протока, м; Ь2 - ширина протока, м. Индекс 0 у тепло-физических параметров характеризует значение при масштабной температуре Т0=1000 оС. В качестве характерной длины выберем длину ванны К Тогда имеем

и = 1,У = ^ = = ^ У0 V V L

*0

*0

Y = = = = L L L L

z . Ь. х, р _ ц

Z. = , В. = ^, X, = , р = , ц = ^.

Ь Ь Ь р0 ц

Здесь h - глубина слоя стекломассы в ванне, м; Ь - ширина ванны, м; х1 - длина загрузочной части ванны, м. Подстановка этих соотношений в консервативную форму записи уравнения движения в проекции на ось х приводит к выражению

± +± -д [р. а+...=о,

l дх L дХ

L

дХ

или

дР д

— + —

дХ дХ

ри -Ади

V Re дХ J

+

д

+—

дY

_д_ дЬ

дР

+

й ди '

рЦУ - — — Re дУ)

¿1 ди

рUW -—-

Re дЬ

+

= 0,

д

— +-

дУ дХ

риу ду

Re дХ

+

д

+—

дУ

_д дЬ

дР

ру 2 -А дУ

V Re дУ у

+

+

рУW -Аду Re дЬ

= 0,

д

— +-

дЬ дХ

й дW

pUW---

Re дХ

+

+

+

_д_ дУ _д_ Ж

_ж„,г й дW

pVW-—--

Re дУ

pw ^ -

V Re дЬ у

+

+ Р = 0,

Fr

где P=p/(р0V02) - безразмерное давление, или число Эйлера; Re=р0V0L/ц0 - число Рейнольдса; Fr=V02/(gL) - число Фруда.

Введем обозначения:

Fl= ри -Ади, С1= риУ-А ди, Re дХ Re дУ

Н1= р UW -А

Re дЬ

F2= риУ-А дУ, 02= рУ2 -А дУ,

Re дХ Re дУ

„ .,„„ Д

Н2= р VW —---,

Re дЬ

Fз= pUW -А Оз= -А

Re дХ Re сУ

Н3= р Ш2 —---.

Re дЬ

Тогда уравнения движения можно переписать в компактном виде:

дР сН1

- + -

- + -

- + -

= 0,

дХ дХ дУ дЬ

дР дF, да, дН,

-+ —2 +-2 +-2 = 0,

дУ дХ дУ дЬ

дР дК дН

+

+

+

+ = 0.

дЬ дХ дУ дЬ ^

Для аппроксимации уравнения неразрывности и уравнений системы (8) воспользуемся шахматной, или разнесенной, конечноразностной сеткой. Осуществим переход к разностным соотношениям на различных контрольных объемах (КО). Нужно отметить, что сетка на самом деле едина (в геометрическом плане) и индексация в различных КО изменяется лишь с целью более четкого выделения искомых параметров. В ряде случаев такая условность затрудняет построение дискретных аналогов (возникает неопределенность индексации), поэтому в данной работе следует связать индексацию с геометрическим центром КО, на котором аппроксимируется уравнение неразрывности. В этом центре определим все параметры состояния среды: давление, температуру и концентрацию.

Используя контрольный объем, соответствующий данным рисунка 1, запишем уравнение неразрывности в виде

[(ри),+1/2,к,1 - (Ри)з-1/2,кД ]AYAZ +

+ [

+ ),,к, 1+1/2 - 0™),, к,1-1/2 ]AXAY = 0-

Но по заданному условию теплофизические свойства берутся в центре КО и далее остаются неизменными во всем объеме параллелепипеда. Тогда во всех слагаемых плотность одинакова и, следовательно, исключается из (9). В результате имеем

(иж/2,к,г - иЬ1/2,к,г) AYAZ + + Кк+т, , - Vj,k-1/1, г) AXAZ + + Кк,М/2 - /-1/2) AXAY = 0.

Можно видеть, что данная форма аппроксимации сохраняет физический смысл уравнения неразрывности как разности между массовыми расходами втекающей в КО и вытекающей из него стекломассы. По сравнению с первоначальной трактовкой изменена индексация компонентов скорости. Согласно физическому смыслу баланса массы определяемыми величинами являются компоненты скорости в геометрическом центре граней, то есть и^^,/, Ц|+1/2,к,/, |-1/2,/, |+1/2,/, и ^|,к,/+1/2.

Перейдем к уравнениям баланса импульсов. В проекции на ось Х имеем дН,

дР дЕ дG,

-+—L +-1- +-

дХ дХ дY &Ъ

± = 0.

В новой индексации не все суммарные индексы попадают в целочисленные узлы. В соответствии с рисунком 2 координаты КО для оси Х изменяются в пределах: по Х от Х| до Х|+1; по Y от Yk-l/2 до Yk+l/2; по Z от Z/-l/2 до Z/+l/2. Применяя метод контрольного объема к уравнению х-компо-ненты импульса (8) с КО, показанным на рисунке 2, можно получить дискретное уравнение

(^+1,кд - ) AYAZ +

+ (Си+1/2, к+1/2,1 - °1л+1/2,к-1/2,1 ) AXAZ +

+ (Н^+1/2,к,г+1/2 - H1j+1/2,k,г-1/2 ) +

+ (р^ - Рца ) AYAZ = 0.

Центральным параметром данного КО является компонент Ц|+1/2,к/, поэтому все аппроксимации выполняются относительно этой величины. После небольших математических преобразований получим следующее уравнение, записанное в каноническом виде:

Я}+112,к,1 из+1/2,к,г + 2 апЪ ип> + пЬ

+ayaz (pj+l,k,г - р,к,г) = 0,

где ^ аЦ) ииЬ означает все конвективные и диф-

пЬ

фузионные вклады из соседних узлов.

ин ,к,1

и] ,к ,1+1

и] ,к ,1-1

Рис. 2. Контрольный объем для аппроксимации уравнения импульса по оси х

Можно заметить, что уравнение (11) нелинейно и, следовательно, может быть решено только методом итераций. Теплофизические параметры на гранях КО можно определить либо как среднегеометрическое из соседних узловых значений, либо как среднеарифметическое.

Аналогично рассматриваются и другие координатные направления. КО для оси у лежит в диапазоне: по Х - от Х| -1/2 до Х|+1/2; по Y - от Yk до Yk+1; по Z - от Z/_/ до Z/+1/2 (рис. 3). При помощи КО, изображенного на рисунке 3, дискретная форма уравнения у-компоненты импульса (второе уравнение системы (8)) может быть записана в виде

(^+1/2, к+1/2,1 — -р2;|-1/2,к+1/2,г) AYAZ + +(С2|,к+1,1 - С2|,к,1) л™ +

+(Н2з,к +1/2,1 +1/2 — Н2j,к+1/2,г-1/2 ) АХА'¥ +

+(Р1к+м - Р^к,г) АХ AZ = 0.

Здесь основной (искомой) переменной является Уцг+щ* После небольших математических преобразований получим уравнение для оси Y, записанное в канонической форме:

я;|,к+1/2,г Ч +1/2

V.....V + ^Х-. +

+AXAZ(Pj,k+u -Р|,к,г) = 0. Наконец, в соответствии с КО, представленным на рисунке 4, для z-компоненты импульса (третье уравнение системы (8)) можно записать

(F3j+1/2,k, 1+1/2 F3j-1/2,k,I+1/2 ) ayaz + +(G3j,k+1/2,1+1/2 — G3j,k-1/2,I+1/2 ) axaz +

+

+(pj,k,/+i - pj,k,/) axay + + pj,k,i + pj,k,l+1 axayaz = 0.

2fr

^-^-Vj.k.M^

Рис. 3. Контрольный объем для аппроксимации уравнения импульса по оси у

• Щ ,к ,1+1

<--->

Рис. 4. Контрольный объем уравнения импульсов по оси г

В данном случае определяющей величиной является компонента Wj,k,^+1/2, поэтому получаем следующее уравнение для оси Z, записанное в канонической форме:

aw W +V aw W +

Aj,k,I+1/2 v,j,k,i+1/2 T ^ Anb v,nb T

nb

+bw + AXAY(Pj,k,;+i -Pj,k,i) = 0. Подведем итоги. Гидродинамическая часть численной модели на данном этапе включает четыре уравнения: одно уравнение неразрывности (9а):

[uj+1/2, k,l - uj-1/2, k,i ]ayaz +

+ [vj, k+1/2, l - vj, k-1/2, i ]axaz +

+[wj,k,i+i/2 -wj,k,i-1/2]axay = 0

и три уравнения баланса импульсов: aj+1/2,k,l uj+1/2,k,l + 2 anb Unb +

nb

+ayaz (pj+i,k,i - pj,k,i) = 0,

aj,k+1/2,l Vj,k +112,1 + 2 3ri> Vnb + nb

+AXAZ(Pj,k+u -Pj,k,i) = 0, aIk,!+1/2 Wj,k,!+1/2 + Z aW Wnb +bW +

nb

+AXAY(Pj,k, i+i -Pj,k, i) = 0. Уравнения баланса импульсов дадут верные значения компонентов скорости, превращающие уравнение неразрывности в тождество, лишь в том случае, когда при их решении будет использоваться правильное поле давления. Но даже тогда для решения системы понадобятся граничные условия. Для решения уравнений дискретного аналога можно воспользоваться методом SIMPLE С. Патанкара. В плане граничных условий основную сложность представляет задание компонентов скорости на свободной поверхности. Для неизотермической задачи здесь можно попытаться воспользоваться эффектом Марангони. Однако для изотермического потока этот выход отсутствует, и приходится изыскивать какие-то другие способы.

Литература

1. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М.: Энергоатомиздат, 1984. 152 с.

2. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978. 592 с.

3. Дзюзер В.Я., Брылина Т.С., Собянин С.Е. Численное моделирование внешнего теплообмена в газоэлектрической стекловаренной печи с подковообразным направлением пламени // Стекло и керамика. 2008. № 5. С. 8-16.

К 40-летию НИИ «Центрпрограммсистем»

С благодарностью примем и разместим в юбилейном альбоме и на сайте института яркие воспоминания всех, кого судьба связала с «Центрпрограммсистем». Оцифруем фотографии, хронику, вырезки из периодики и другие материалы, приобретем старую вычислительную технику, приборы и тому подобное для экспозиции музея.

Контактные телефоны: (4822) 399-228, 399-288,

сайт в сети Интернет: ЦПС.РФ