книги / Моделирование физико- химических процессов нефтепереработки и нефтехимии

Глава V

Методы проведения расчетов и исследования

математических описаний

Выше показано, что математические описания химико-техноло­ гических процессов представляют собой системы алгебраических или дифференциальных уравнений. Здесь приведем описание не­ которых численных методов, позволяющих выполнять расчеты таких систем. Далее рассмотрим существенные для математиче­ ского моделирования методы исследования таких систем: опреде­ ление чувствительности решения к величинам параметров и, если число возможных решений больше одного, — определение устой­ чивого решения и на его основе — устойчивого режима работы химико-технологического процесса.

1- ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ СИСТЕМ АЛГЕБРАИЧЕСКИХ УРАВНЕНИИ

Укажем вначале, что решением системы алгебраических урав­ нений

/х (*1» ••• •» xk) —О

(V .4 )

f p ( xi i •••> x k) —О

является такой набор (наборы) чисел oslt ..., ak, который в случае замены ими х г, ..., хк обращает систему (V.1) в тождество.

Наиболее обстоятельно разработаны методы решения системы

aplxl"t* арЗх2~Ь * •* Н“ apkxk— dp

когда р > к, обычно применяют описанный в главе I метод наи­ меньших квадратов, который позволяет свести ее к нормальной системе (при р = к). Если р = к и система (V.2) является опре­ деленной, т. е. имеет одно решение, его можно найти описанным в главе I методом, по которому х,- = Аг/А. Здесь А — определи­

ли

тель, составленный из коэффициентов левой части системы (V.2)', a получается из А заменой г-того столбца на столбец d :

Пример такого расчета дан в главе I. Укажем, что по извест­ ной теореме Крамера, система (V.2) является определенной, если А Ф 0. Другой метод точного решения системы линейных уравне­ ний (Гаусса) приведен ниже (стр. 201).

Для сложных систем (V.2) при больших к нахождение точ­ ного решения потребует выполнения большого числа расчетов; поэтому часто ищут не точное, а приближенное решение этой системы, используя различные итерационные методы. Как пра­ вило, программы для ЭВМ при использовании итерационных мето­ дов значительно компактнее и время вычислений гораздо меньше. Известен [1] ряд итерационных методов решения системы (V.2), однако каждый из них применим лишь в ограниченной области условий, позволяющих быстро свести итерационным процессом «плохое» решение к хорошему. Вне этой области сходимость решения будет медленной.

Проиллюстрируем часто используемый метод простой итера­

диагональные элементы (подчеркнуты) по модулю больше суммы модулей недиагональных элементов строки. Запишем систему (V.2) в матричной форме: А х = d- Решая ее относительно х , перейдем' к эквивалентной системе вида:

Для перехода к этой системе г-тое уравнение системы (V.2) нужно решить относительно xt, так чтоx

Выберем начальное приближение — вектор-столбец х 0:

(*01

#0= ' x0kd

142

Можно, например, принять ж0[- = 0 или x Qi = dt. Далее по уравнению (V.3) определим улучшенное решение

x i — В х о ~(~ ф

и т. д., так что в результате г итераций получим:

а>+1= Я®г+ф

Если точное решение се может отличаться от найденного на величину меньшую или равную е, то необходимое число итераций

п рд + и> r - In (i- p )

где р, q, w — максимальные из величин |ац!аи [, |ж01|, \di!aii\. Ряд других итерационных методов описан в монографии [1]. Отметим, что при небольших к итерационный процесс теряет

смысл.

Рассмотрим теперь возможные методы решения системы (V. 1) для случая, когда среди функций / есть нелинейные. При этом возможно несколько решений этой системы, причем определить их можно в общем случае только численными методами.

Пусть, например, решается одно уравнение с одним неизвест­ ным и изменение х допустимо в интервале от а до Ъ. Проверяем /(a ), f (а-\-Ы 2), f (Ь) и выбираем для дальнейшего поиска ту половину начального интервала, на концах которой / (я) имеет противоположные знаки. Повторяя такой поиск, сузим интервал неопределенности до сколь угодно малой величины е. Число ите­ раций при этом

Укажем, что решение нелинейного уравнения с одним неизвест­ ным / (ж) — 0 можно рассматривать как задачу поиска минимума функции F (я), для которой / (я) = dF (x)ldx. Такая задача ре­ шается поисковыми методами (половинного деления, золотого сечения, стохастической аппроксимации), рассмотренными в главе VI.

Понятно, что решение системы (V.1) значительно сложнее решения одного нелинейного уравнения, но и для нее можно пред­ ложить итерационную процедуру. Прибавив я, к правой и левой части каждого i-того уравнения системы (V.1), перейдем к системе вида (р = к):

143

Система (V.4) похожа на аналогичную (V .3), использованную

вметоде простой итерации для линейных уравнений. Применим

ив этом случае ту же итерационную процедуру. Пусть после г

Для нелинейных функций фх, ..., срк уже не приходится выби­ рать, быстро или медленно сходится итерационный процесс. Важно, чтобы он вообще сходился. Показано [2], что сходимость будет выполняться, если начальное приближение достаточно близко к точному решению (а х ... ак) и если имеет место хотя бы одна из систем неравенств

Поскольку, однако, осх ... ak заранее неизвестны, то лишь реализация итерационного процесса может^выявить его эффектив­ ность.

Другой возможный метод нахождения решения нелинейной систем ы '(V .1) заключается в переходе к поиску экстремума функ­

ции многих переменных. Если ввести такую функцию F , что

T O F = F {x± ... Xk), и ее минимум достигается при таком наборе

величин otj, ..., <хк, заменяющих x v ... xk, который одновременно удовлетворяет решению системы (V.1). Методы поиска экстре­ мума функции многих переменных рассмотрены в главе VI.

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

A(GcT) = -qnvW (Ci)V + lcTF (Тва~ Т ) )

144

можно перейти к решению системы дифференциальных уравнений:

dCi

причем при достаточно больших т получим решение стационарной задачи.

2. ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

Рассмотрим вначале метод интегрирования одного дифферен­ циального уравнения

dy/dx = f ( x , у) у (х о) = i/о

с использованием разложения в ряд Тейлора. Задав небольшое

Производные порядка выше первого можно вычислить по

дУ

пли

V l - f x l k + fkf’y.k

Однако такое определение высших производных и нахождение по ним итеративным методом у при к =*= 1, 2 ,... неудобно, так как численное дифференцирование связано со значительными труд­ ностями. Поэтому ряд Тейлора используют лишь для оценки точности других методов, в которых ограничиваются вычислением производной только первого порядка, т. е. / (ж, у).

При численном решении систем дифференциальных уравнений наиболее часто используют методы Эйлера и Рунге — Кутта. С другими методами можно ознакомиться в книгах по вычисли­ тельной математике [2, 3]. Оба эти метода удобны при програм­ мировании решения на ЭВМ для тех случаев, когда все граничные (начальные) условия заданы при одном и том же значении аргу­ мента. Охарактеризуем кратко эти методы.

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

dy/dx —f (х, у) У Ш = Уо

145

Перейдем к уравнению, где dx заменяется конечным постоян­ ным шагом к. Тогда на шаге к

причем yk = у (®0 + кк).

По этому соотношению, начиная от у 0 и придавая к последо­ вательно значения 0, 1, 2, ..., можно найти у для любой заданной величины х.

Приведенный метод может быть без каких-либо изменений использован и для систем дифференциальных уравнений; в этом случае в приведенных соотношениях величины у и / следует рас­ сматривать как векторы или использовать их для любого из у(, образующих систему:

малых к дает ошибку, пропорциональную к2, т. е. равна const•к2. Если интервал х 0—х разбит на п частей и к — (х —х 0)/п, то такая ошибка совершается п раз, и суммарная ошибка будет пропорцио­ нальна ( а я 0)2/гс. Таким образом, увеличение точности в п раз требует увеличения в то же число раз точек деления. Именно этот недостаток ограничивает применение метода Эйлера. Если, од­ нако, зависимость у (х) близка к линейной (что довольно часто имеет место в прикладных расчетах), то коэффициент пропорцио­ нальности const мал, и метод Эйлера даже при небольших п даст точное решение.

Метод Рунге,— Кутта позволяет получить более высокую точ­

a*=/|j*fc-1—l " ’ Ук+ Н хЬ

т. е. производная определяется не в начальной, а в средней точке интервала (при х = xk -f- hi2)..

Для такого расчета найдем, пренебрегая-членами с к3, к* и т. д.,

что

№

t ^ . У k , _ yk*1—У к = У ф - \ - ^ 2Т ^ 8

т. е. суммарная (а не на одном шаге) ошибка пропорциональна А2 или обратно пропорциональна /г2. Таким образом, если при расчете по методу Эйлера с увеличением числа шагов в п раз ошибка уменьшается также в п раз, то по методу Рунге — Кутта она уменьшается в п 2 раз.

146

Точность расчета можно увеличить, если перейти к интерполя­ ционной формуле

1

Ук*1 — Ук +*б" [/ fob У*Ч-2об* + 2р*-|-Y*)1А

Интерполяционная формула дает суммарную ошибку, пропор­ циональную № или обратно пропорциональную п4 [3]. Как и прежде, соотношения для систем дифференциальных уравнений останутся теми же; в интерполяционных формулах у и / можно рассматривать как векторы или записывать их для каждого из уп образующих систему:

В расчетной практике необходимая точность достигается про­ ведением расчета при h t и h z = (1/2) h v Если результаты расчетов

Поскольку при описании процессов дифференциальными урав­ нениями второго и более высоких порядков * граничные условия могут быть заданы в разных точках (так называемая краевая задача), численные методы для этих случаев должны быть моди­ фицированы. Например, химический процесс в зерне пористого катализатора радиусом R , описываемый уравнением D (diCldxi) =

— f (С), обычно характеризуют краевыми условиями для концен­

заменить эквивалентной системой к уравнений первого порядка [например, приведенное уравнение второго порядка можно заме­

* Напомним, что порядком дифференциального уравнения называется максимальный порядок производной; для системы дифференциальных урав­ нений он соответствует числу уравнений первого порядка, содержащих про­ изводные различных функций. Число граничных (начальных) условий должно быть равно порядку.

147

Будем считать, что эта система имеет решение, притом един­ ственное. Наиболее часто такое решение находят численными методами, которые сводят краевую задачу к задаче с граничными условиями на одном конце (задача Коши). Если, например, к—р граничных условий заданы при х = а, а р условий — при х = b (фиксированные условия), то, выбрав р произвольных условий при х = а, будем решать задачу с к условиями при х = а (при этом р условий при х = Ъ не используются). Произвольные условия при х = а меняют таким образом, чтобы рассчитываемые У (b) удовлетворяли отброшенным фиксированным условиям.

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

сдвухточечными краевыми условиями у (а) = а , у (Ь) = р.

Заменим производные у " и уг конечными разностями для участка

и система уравнений (V.5), (V.6) — краевая задача — заменяется

148

выше (етр. 143), для ее решения созданы специальные методы, например метод прогонки. В этом методе уравнение типа (V.8) преобразуется в уравнение

yj = Mjyj+i+ Nj

где, как можно убедиться (подробнее см. [3]), Mj и Nj равны:

причем М 0 = О и iV0 = а .

Теперь легко осуществить итерационную программу вычисле­ ний Mj, Nj, а затем по ним и у/ в точках 1, 2, ..., п—1. Создан ряд вариантов сочетания методов конечных разностей и прогонки, когда краевые условия заданы не только в виде чисел, но и в виде функций.

3. МЕТОДЫ ЧИСЛЕННОГО ИНТЕГРИРОВАНИЯ УРАВНЕНИЙ В ЧАСТНЫХ ПРОИЗВОДНЫХ

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

Возможно и задание более сложных граничных условий, например, вида

Фх |Ук, у к |ж=0 = 0 ф2 |ук,у'к = 0 для любых т

Уравнения типа (V.9) относят к параболическим и, как и лю­ бые уравнения в частных производных, решают методом сеток. По этому методу всю область изменения х и х делят сеткой (рис. V-1). В узлах сетки рассматриваются функции дискретного аргумента (сеточные функции); на сетке производные заменяют отношением конечных разностей. Точность метода зависит от выбора сетки и способа аппроксимации производных. Координаты узла сетки в точке if, очевидно, следующие:

149

Функцию в точке ij обозначим у\. Достаточно ясно, что для

расчета первой производной dyidx в точке ij можно использовать две точки линии сетки при неизменном ту:

Рис. V-l. Сетка для численного интегрирования функции двух пе­ ременных: т и х.

При расчете суммы различных производных нужно использо­ вать несколько точек сетки, находящихся на разных линиях. Например, для расчета вблизи точки ij величины |ду/дх—д2у1дх21;/ можно использовать у в четырех близлежащих точках сетки: (£—1» j) (*> j)i (2 + 1> j) и (г, / + 1) — так называемый четырех­ точечный шаблон:

Производные, входящие в уравнение (V .9), рекомендуется вычислять по следующим разностным уравнениям четырехточеч ного шаблона:

150