Полак_и_др__Вычислительные_методы_в_химической_кинетике (972296), страница 36
Текст из файла (страница 36)
Если все Х; отрицательны, то Яр(ехр(Аг)) падает с ростом г, в случае наличия положительного корня наблюдается рост этой величины. По изменению величины Яр(ехр(АгЦ можно оценивать и величину положительного корня, и разброс отрицательных корней, т.е. разброс характерных времен процесса. В описанном алгоритме не производится оценка положительного корня, так как в методе зта оценка делается по величине глобальной ошибки, однако такая оценка может быть полезной в процедурах дифференциального спуска и нахождения производных по параметрам. В работе [1) предложен аналогичный описанному метод решения линейной задачи. Некоторью отличия состоят в вычислении матричной экс- 145 в противном случае последующие действия определяются установленным режимом.
В случае режима *'выбор шага" полагаем И = Ы2, переходим к выполнению п, 3, в случае режима "удвоения" устанавливаем режим "сетка" и переходим к выполнению и, 9, в случае режима "сетка" перь ходим к выполнению п. В. 5. Полагаем г = 1+и, у(г) = у(г) + ц (И). Используя формулу (5.53), вычисляем б у'(г) . Если установлен режим "сетка", переходим к выполнению п. 7, При других режимах переходим к выполнению и. 6. 6. По функциональному уравнению (5.43) вычисляем значение С (2И), полагаем И = 2И, устанавливаем режим "удвоения".
Переходим к выполнению п. 4. 7. Если сделано и шагов в режиме "сетке"', где и — размерность систе. мы уравнений (это означает, что на новом шаге интегрирования в уравнении (5.40) используется якобиан, полученный на предыдущих шагах. Таким образом, при решении уравнения (5.40) в режиме "сетка" не пересчитывается оператор С(И)), или величина глобальной ошибки превышает заданное значение, то переходим к выполнению и. 2, в противном случае — к выполнению п. 4. В.
Если не исчерпан запас запоминаемых матриц С(И) (в режиме "удвоения" запоминается 10 последних матриц С(И) ), то выполняем и. 9, в противном случае — п. 2. 9. Полагаем И = И72, С(И) = С(И(2). Переходим к выполнению п. 4. С помощью программы, построенной по такому алгоритму, было проведено моделирование реальных химических процессов [123, 141). В работе [123) исследовался процесс окисления углеводородов в присутствии двух ингибиторов, в работе [141) — реакция фторирования дифторметана, причем в этой реакции наблюдались большие времена задержки и наличие резко растущего решения по отдельным компонентам, что свидетельствует о наличии положительного собственного значения в спектре якобиана. Заметим, что если спектр якобиана принадлежит положительной области, то, как видно из уравнения (5.52), полная ошибка вычислений будет быстро расти, так как она умножается на экспоненту в положительной степени.
Таким образом, на участках, где якобиан имеет положительный корень, необходимо ограничивать шаг интегрирования. Так как в алгоритме вычисляется матричная экспонента, то можно непосредственно выявить положительный корень якобиана. Действительно, пусть Хг — собственные значения якобиана. Ызвестно, что [37) поненты на начальном шаге интегрирования. Здесь используется не разложение в ряд Тейлора, а приближение в виде (Š— АИ) '.
В случае действительного отрицательного спектра матрицы А зто позволяет выбирать достаточно большой начальный шаг интегрирования, однако при наличии положительного собственного значения такая аппроксимация в области М>1 оказывается неверной. Поэтому при выборе шага необходимо анализировать положительную часть спектра матрицы А, что, по-видимому, сведет на нет выигрыш, полученный от начального увеличения шага.
Оценка полохачтельного корня по методу минимакса, как правило, невозможна, так как положительный корень по модулю может быть меньше модуля больших отрицательных корней. В работе (161] предлагается вести решение линейной системы по формуле ух+~ = ух +ехр(АИ) ух. (5.57) где у„— решение на предыдущем шаге интегрирования. Преимущество этого способа решения линейной системы по сравне- нию с расчетом интеграла от матричной экспоненты с использованием рекурренты (5.461 состоит в выигрыше необходимого числа операций.
Однако в вычислении по формуле (5.57), на наш взгляд, есть два суще. ственых недостатка. Во-первых, вычисления матричной экспоненты по формуле Пойя ме- нее устойчивы, чем по рек урренте ехр АИ = Е + В, ехр 2АИ = Е + 2В + Вз, В [2И) = 2В(И] 4 В(И) В(И), (5.58) если в матричной экспоненте присутствуют недиагональные элементы, существенно отличные от диагональных. Но в этом случае нет выигрыша в числе операций. Во-вторых, законы накопления ошибок в уравнениях (5.571 и (5.461 различны.
Запишем линейное приближение для погрешности решения линейной задачи с использованием формулы (5.571: Ь у„, = Бух + ехр(АИ) Ье„. (5.59) Видно, что даже в случае действительного отрицательного спектра матрицы А ошибка может накапливаться. Аналогичное уравнение для ошибки при решении линейной системы по формуле (5.46) имеет вид ЬС„+, = ЬС„ехр(АИ) + ехр(АИ) ЬС„. (5.60) В случае устойчивой матрицы А ошибка будет нивелироваться эа счет члена ехр(АИ], т.е. такое решение линейной системы обладает свойством устойчивости. Основу всех методов локальной линеаризации составляют методы интегрирования систем линейных дифференциальных уравнений, т.е. они так или иначе связаны с приближенным вычислением матричной экспоненты.
В работе (95] предложено однополюсное дробно-рациональное приближение экспоненты в комплексной области. Известно, что неявные методы Рунге-Кутта при интегрировании линейной системы дифференциальных уравнений приводят к дробно-рациональной аппроксимации Падэ и, следовательно, трудоемки, так как фактически требуют обращения матричных многочленов.
Неявные линейные многошаговые методы дают аппроксимацию ехр(Аг) главным корнем р(Аг) характеристического 144 полиноме данного метода. Если Р— порядок метода, то р(Ас) = ехр(Ас) + + 0(сл+'). Неявные линейные многошаговые методы менее трудоемки по сравнению с неявными методами Рунге-Кутта, но у них малый порядок точности, как у обычной неявной схемы, либо неверная асимптотика в левой полуплоскости, что свойственно более точным схемам. Препложенное дробно. рациональное приближение экспоненты частичными суммами ряда экспоненты по полиномам Лагерра в значительной мере лишено недостатков, свойственных раэностным методам. Следуя работе [95), представим матричную экспоненту в виде ряда хз ехр(Ас) =ехр — = [1+с)~+ Х Еи(х) ( — г)и, )г) (1,а> — 1, (5.61) 1+с и"-О глез (Ас/х) (1 — Ас/хГс; х > О;а — целое: Еи (х) — многочлен Лагерра.
Для приближения экспоненты используются частичные суммы этого ряда: ги — 1 Ри г (х, а, Ас) = (1 — Ас/х) ' Х Е„(х) [(-Ас/х) (1 — Ас/х) ')". и=э (5.62) В работе получены оценки для оптимальных значений параметров х, а, т изусловия минимума модуля ошибки )б,„. с (х,а, Ас) )= ) ехр(Ас)— Р,„с (х, а, Ас) ) на от(хицательной полуоси при фиксированных вычисли. тельных затратах применительно к решению жестких систем.
Показано, что оптимальные в этом смысле значения параметров х иасвт. Лолучены также асимптотические оценки ошибки частичных сумм (5.62) при оптимальных значениях параметров как в отрицательной полуплоскости л/с, так и в положительной окружности с центром в нуле и радиусом М/ /т. Метод является жестко устойчивым, с большим радиусом в правой полуплоскости -) М/т (. В работе [95) также исследована эффективность вычисления матричной экспоненты с помощью ее аппроксимации частичными суммами вида (5.62).
Показано, что этот алгоритм является эффективным при работе с ленточными матрицами. Необходимо записать обратную матрицу (Е— — Ас/х) ' в факторизованном виде, а затем использовать процедуры перемножения векторов и матриц на факторизованную матрицу (Š— Ас/ /х) г. Вычисление полиномов Лагерра нужно производить с помощью известных рекуррентных формул. В заключение отметим, что в настоящее время метсщы локальной линеаризации становятся все более попупярными для решения обыкновенных дифференциальных уравнений.
Особенно зто касается решения жестких систем, в которых линейная задача во многом является определяюшвй. Большое распространение этих методов связано с тем, что они используют хорошо разработанный аппарат линейной алгебры. Это, в свою очередь, облегчает апгоритмиэацию метода для программирования на ЭВМ. В случаях, когда решаемая численно система жестких обыкновенных дифференциальных уравнений существенно нелинейна, наиболее оправданным, по-видимому, является применение разностных методов, из которых в настоящее время наиболее эффективным является метод Гира.
В случаях, когда спектр якобиана содержит большие положительные собственные значения, целесообразно использовать методы локальной линеаризации. 14У 4. пРимеРы Решений пРямОЙ кинетической 3АдАчи Ниже будет рассмотрено несколько примеров решения прямой кинетической задачи. На зтих примерах мы попытаемся проиллюстрировать, какая физико-химическая информация может быть получена при решении прямой кинетической задачи, а также какие вычислительные трудности возникают в рассматриваемых конкретных случаях.
Моделирование кинетики окисления метане Рассмотрим решение прямой кинетической задачи дпя реакции окисления метана, следуя в основном результатам работы [23) . Механизм процесса приведен в табл. бЛ. Он основан на предложенной Н.Н. Семеновым радикально-цепной схеме процесса окисления метана. Стадия (1) — зарождение цепи, стадия (6) — разветвление цепи, стадии (2) — (б) и (7)— (12) — продолжение цепи, стадии (13) — (17) — обрывы цепи. Исходные концентрации метана и кислорода задавались в соотношении [СН4): [Оз) = 0,29: 0,71. Интегрирование системы обыкновенных дифференциальных уравнений химической кинетики проводилось методова Рунге — Кутта с автоматическим выбором шага с относительной погрешностью 10 -10 ~, одмако в соответствии с предложенным в [22, 23] алгоритьюм интегрирования систем жестких дифференциальных уравнений (см. раздел 2) полная система обыкновенных дифференциальных уравнений заменялась укороченной, совместно с которой решалась система алгебраических уравнений для концентраций "быстрых" компонент СН700, ОН, НСО.
В данном случае расчеты упрощались тем, что алгебраические уравнения оказапись независимыми. За счет применения принципа квазистационарно- Тв бп и На 6.1 Мвааннзм пкиепения метена Е, нкап)мопь реакции Номер реакции ' Размерность А реакций — п. моль ' с ~. 148 (1) (2) (31 (4) 60 (6) (7) (8) (9) ( го) (11) (12) (13) (14) (16) (16) (17) СН, +О, -СН, +НО. СН„+ О, — СН,ОО СН 00 СН О+ОН СН, +ОН СН +Н,О СН,О+ ОН -Н,О+ НСО СН10+ О, НСО+ НО. НСО+ Оь СО+ НО, СН +НО, Н~О, +СН, СН,О+ НО, Н О, + НСО СО+ ОН СО, + Н СН +Н СН +Н, СН О+Н НСО+Н, СНь +СН С,Н„ НО, +Нов Н,О, +О, ОН + НО Н,О'+ О н+но, н +О, СНь + НО, СН„+ О, дпп реакции (3)- с ',дпп остальных 1),О 8,0 13,0 11,0 12,0 Н,о 8,0 11,0 11,0 1 1,0 10,1 8,6 10,0 10,0 10.0 10,0 10,0 55,0 0 20,0 8,5 6.3 32.0 0 19,76 8,6 7,0 тцг 2,0 0 0 0 0 0 су, а6.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
