Полак__Применение_вычислительной_математики_в_химической_и_физической_кинетике (972294), страница 31
Текст из файла (страница 31)
С целью автоматизации этой процедуры нами составлена стандартная программа для ЭВМ типа М-20, реализующая метод оврагов (см. !47, 62[, а также стр. 103 настоящей книги). Согласно этой программе, отыскиваются такие значения констант скоростей, которые обеспечивают достижение минимума суммы квадратов отклонений расчетных величин концентраций от опытных. Коли ошибки последних подчиняются нормальному закону, то определяемые таким образом константы скоростей являются наиболее вероятными для данного механизма.
Мы считаем необходимым еще раз обратить внимание на то, что метод оврагов является нелокальным методом нахождения минимумов функций многих переменных. Другими словами, в районе начальных значений констант определяются все (если их несколько) минимумы суммы квадратов отклонений, так что имеется возможность определить координаты наиболее глубокого минимума.
Как отмечалось на стр. 102, метод оврагов наиболее эффективен в случае отыскания минимумов так называемых «хорошо организованных» функций. Под «хорошей организацией» в данном случае «Значение Ьг ~( О,! часа обусловлено максимально(! величпнои шага интегрирования (ОД часа); выбор последнего мотивирован ниже. следует понимать то, что все константы скоростей можно разбить на две группы, причем в одну из них войдут те константы, изменение которых существенно меняет сумму квадратов отклонений, а в другую — константы, изменение которых мало влияет на изменение суммы квадратов отклонений. Чем меньше меняет сумму квадратов отклонений изменение константы на некоторую величину, тем больше диапазон неопределенности этой константы, т.
е. тем больше ошибка, с которой она определяется. В настоящей работе за ошибку констант принималась такая величина о (й), прн которой изменение )с от )г — о ()с) до Й + о ()с) оставляло расчетные концентрации в пределах 5%-ного коридора ошибок (см. стр. 111). Перейдем теперь к методике расчета концентраций компонентов реакции. В кннетике обычноприннмается, что концентрации активных частиц (в нашем случае Н„)), и 80) стацнонарны практически в течение всего времени реакции.
Однако для рассматриваемой системы такая стационарность имеет место лишь прибольших временах реакции и, напротив, не существует при г (20 час. Поскольку абсолютные значения скоростей изменения [(Р),! и [(Н),! невелики, то при выборе достаточно малого шага * можно записать условия квазистационарностии для г (20 час в виде системы нелинейных алгебраических уравнений '( [(Н) ! 0 ~([(")т! 0 '( (об) 0 (98) ш и ' лт Подчеркнем, что уравнения (98) справедливы только в пределах выбранного шага интегрирования. Эта система решалась итерационным методом Ньютона [105), причем необходимые для решения производные соответствующих правых частей по концентрациям активных частиц находились численно [195), что позволило сделать программу стандартной.
Система уравнений квазистационарности решалась при тех временах облучения, для которых имелись опытные данные по составу газовой фааы. Необходимые для расчета величины [ИОН[ и [ЯО()! определяли по уравнениям баланса. Полученные значения [(Н),), [П),! и [ЯО! аппроксимировали полиномом по времени [В) = Ь, + Ьг(+ Ьтг', (99) коэффициенты которого находили методом наименьших квадратов.
Из формулы (99) следует, что при ( = 0 [В) = Ье, в то время как строго при ( = 0 [В! = О. Однако введение члена Ь, вполином практически не сказывается на результатах расчета, поскольку установление квазистационарного режима происходит за столь малые времена, что изменением концентрации стабильных веществ при этом можно пренебречь. Между тем, если положить Ь,:= О, то е Шага, достаточно малого для соблюдения уранненнй (98) с заданной точностью. для удовлетворительного описания зависимости К = у (1) требуется уже другая аппроксимация. Заметим, что величины квази- стационарных концентраций (Н)„(П), и (80) зависят от коккретных значений констант скоростей.
Поэтому по мере уточнения последних уточнялись и концентрации активных частиц. При такой методике расчета систему дифференциальных уравнений кинетики удалось свести к пяти дифференциальным уравнениям (для скоростей изменения концентраций Н, П, НП, ЗОН и ЗОВ). На каждом шаге интегрирования этой системы концентрации (Н)„(П), и 80 вычислялись по соответствующим полиномам второй степени.
Обычно численное решение задачи Коши проводится методом Рунге — Кутта. В процессе нашей работы выяснилось, что можно почти в 2 раза сократить объем вычислений, не теряя при этом точности, если для решения задачи Коши использовать алгоритм, предложенный недавно Фелбергом (196]. Для ускорения расчетов выбирали шаг, который обеспечивал выполнение балансовых уравнений с точностью до 2% .
Значение времени, соответствующее такому шагу, зависело от конкретных значений констант и участка кинетической кривой. Естественно, что минимальнын шаг был принят на начальном участке, но и при больших временах реакции максимальная величина шага не превышала 0,1 часа. Рассмотрим теперь результаты расчета механизмов № 1 — 4. Механизм № 1 был отвергнут, поскольку при удовлетворительном описании экспериментальных данных константы скоростей внутри рассматриваемых групп реакций отличались на 2 — 3 порядка (нарушение требования 9), что, конечно, не имеет физического смысла.
Заметим, что в работе (192) механизм, основанный на реакциях группыВ, также был отвергнут, но ко другой причине: в связи с тем, что радиационная деструкция гидроксильных групп происходит с недостаточным для объяснения наблюдаемой скорости обмена радиационно-химическим выходом (6 0,1; тогда как для исследуемого процесса С 2 — 5). Механизм № 2 был отвергнут, поскольку при его расчете не удается получить результатов, соответствующих эксперименту (см. требование 4). Существенно наличие качественных расхождений — не удается получить максимума, который имеется в экспериментальной зависимости концентрации НП от времени. От механизма № 3 пришлось отказаться в связи с невыполнением требований 7, 8 и 10.
Наиболее подходящим (удовлетворяющим всем перечисленным выше требованиям) оказался механизм № 4. В связи с этим рассмотрим результаты его расчета более подробно. Как было отмечено выше, механизм № 4 включает реакции групп А, Р, Е, Е и 6. Следовательно, для скоростей соответствующих реакций мы имеем систему из четырнадцати дифференциальных уравнений, которая для решения сводится (по числу реагирующих компонент) к системе пяти дифференциальных и трех алгебраических уравнений вида 146 (00). В принципе, используя уравнения баланса, можно свести зту систему и к двум дифференциальным уравнениям. Однако мы предпочли оставить все уравнения баланса дчя контроля за результатами расчета.
Чтобы больше не возвращаться к атому вопросу, подчеркнем, что все контрольные уравнения баланса выполнялись с ошибкой не больше 2 %. Выполнение требования 4 видно из данных расчета, графически представленных на рис. 34. Выполнение требования 6 обеспечивалось, как это видно из изложенного выше, о, и акино ооо Рнс. 35. Зависимость количества атомов Н(б) и [) (+) на поверхности силикаголя от количества водорода и дейтерия в газовой фазе аоо А — величины 21Н,) + [НМ к 2[П,) + +[НШ,  — величины [Н) к [П]. Тачки соотззготвуют временам облучеякя, указаккым ка рмс.
Зо зо 6 часы (нй [ш, ц чяоы [нп вял 1,5 5 10 15 0,9 1,4 2,0 2,4 2,3 1,9 1,5 11 7,9 6,8 5,6 4,9 4,7 4,5 4,7 г 2,6 2,6 17 20 25 1,0 0,9 0,9 Подчеркнем, что для проверки выполнения требования 8 важны не абсолютные, а относительные значения концентраций о. ч Заметим, что абсолютное значение концентрации хомосорбнрованных атомов водорода) ((Н), + ([))л) = (3,4 + 0,1) ° 10п оло з, тогда как, согласно работе (194), койцонтрация радиационно-нндуцпрозанных центров адсорбцпи водорода на позорхностп сплякаголя составляет 3,3 ° 10п сло ' 147 о ою ~оо Л и ззяолыул самой методикой расчета. Выполнение требования 7 иллюстрируется рис.
35. Кстати, тангенс угла наклона прямой на рис. 35 приблизительно равен 0,01, что подтверждает выполнение требования 6. Выполнение требования 8 видно из приводимой ниже табл. 11. Таблица И Зависимость концентрации промежуточных продуктов (в 10 " олч з) от времени облучения Из приведенных виже значений констант скоростей, полученных на ЭВМ при расчете механизма М 4, видно, что требование 9 удовлетворительно выполняется.
Значения констант приведены по группам (размерность всех констант выражена в еек '), так как расчет был проведен для безраамерных концентраций: Группа Л (Ул =- (0,95~ 0,06),10-4. «ц (1 4 х 0 4) 10-4. Усп1 = (1,3~0,1) 10 ' Группа УЭ (кх =- (6, 7 х 4,5) 10 '"; Ус и =- (9,9 зс 2,4) 1 О '-', йхп = = (4,2+ 2,5) 10 "; Уехтп = (6,2 +1,4) 10 "). Группа Е(Усдч — — (2,4+ 0 7) 10-ез Угхч = (1,0 ~0,5) 10 "). Группа Р(Усхч,=(4,3~2,1) 10 "; кхчп =-(6,2+1,4) 10 "' Уехчгп = (2,4 ~ 0,2) 10 ") Группа С(Усх~х = (2,1~0 4) 10 "' кхх = (1,1*0,4) 10 "). — — == (Уех1 + йхтч) (5ПН)е [(В)..)~-о.