Полак_и_др__Вычислительные_методы_в_химической_кинетике (972296), страница 49
Текст из файла (страница 49)
Е) =о>)" к(е)л(етие' — ь>] к(е')л(е,г)гте' — lг(е)л(е,г). (820) Эг времени, и в пределе соответственно стационарная функция распределения, которая может быть получена также и решением задачи на собственные значения для этой системы дифференциальных уравнений.
Следуя работам [138, 139), перейдем от функции распределения л(е, Г) к функции Г(е, г), являющейся поправочным фактором для равновесной больцмановской функции распределения: п(е,г) г(е,г) ' К(е) ° (8.26] Подставим это выражение в основное кинетическое уравнение (8.13) и, используя принцип детального равновесия и условие нормировки, получим для Г(е, Г) следующее уравнение: дК(е, ) — = ш))с (е; е) [ Це', г) — 7 (е, г ) ) сне' — Л (е)7(е, г ) . (8.27) дг о Для того чтобы провести дискретизацию задачи, необходимо ограничить интервал рассматриваемых энергий.
Из физических соображений ясно,что искомая функция распределения будет пренебрежимо мала при энергиях порядка нескольких энергий диссоциации. Для примера укажем, что в модели сильных столкновений квазистационарная функция распределения становится пренебрежимо мала при энергиях, превышающих критическую энергию Ео на величину порядка [Зеб) (чТ„что для больших молекул при не слишком высоких температурах составляет десятые доли от энергии диссоциации. Исходя из этих соображений, заменим бесконечный предел интегрирования в уравнении (8.27) на конечный и обозначим его Еь. Разбив рассматриваемый энергетический интервал на (У частей — (е,, ез, .... ен ), получимвточкахдискретизациисистемудифференциальных уравнений: 87(еп г) = ш) (г (е', е~) (7(е', г) — К(ео г!) оу' — д (ег) .
г (еь г), дг о (8.28) г' = 1, 2, ..., /У. Если интегралы в правой части системы вычислять с помощью тех ипи иных квадратурных формул, то пЬлучится конечная линейная система дифференциальных уравнений с постоянными коэффициентами: Ф =А 1, Ф(0) =(о, (829) где ( = ( у (еы г), Г (ез, г),..., 1(ек, г)); А — матрица коэффициентов. Решение этой системы уравнений имеет вид ( = ехр(Аг) (о. [8ЗО) Для того чтобы решение системы (8.28) удовлетворяло физическим требованиям, накладываемым на функцию распределения, необходимо.
чтобы при дискретизации задачи, т.е. при переходе от интегродифференциального уравнения к конечной системе линейных дифференциальных обыкновенных уравнений„матрица А сохранила свойства интегрального оператора уравнения (8.13) . Эта матрица должна быть симметризуемой и отрицательно определенной, т.е. обладать действительным отрицательным спектром.
Используя функцию 7 (е, г), можно получить константу скорости моно- молекулярной реакции по кинетической кривой методом наименьших квадратов. Константа скорости мономолекулярной реакции и поправочный фактор 7 (с, г) связаны с квазиравновесной функцией распределения соотношением лП3 = К(е) 7(е, г) ехр((гг). Решение нестационарного уравнения позволяет определить не только стационарную функцию распределения, но и время существования квазистационарного режима. В математическом смысле этот период определяется вторым по величине модуля собственным значением оператора интегродифференциального уравнения (8.13). Константа скорости мономолекулярной реакции и квазистационарная функция распределения могут быть найдены также и с помощью решения задачи на собственные значения. В случае дискретной задачи вычисление этих величин сводится к определению минимального по модулю собственного значения матрица А и соответствующего ему собственного вектора.
Для поиска собственного значения н соответствующего собственного вектора используется метод обратной итерации с релеевским сдвигом [142): (Ах„, ха) (хь хл) у = (А — )глЕ) хп ° х„,, = у/()у(! . (8.31) Процедура обратной итерации эффективна при решении частной проблемы собственных значений, т.е. при выделении одного собственного значения. В процессе вычислений приходится иметь дело с плохо обусловленной мат. Рицей, однако процедура быстро сходится, если удачно выбрано начальное приближение. Оператор(А — )г„Е) ' увеличивает проекцию вектора, на который он действует, в направлении собственного вектора и подавляет все остальные проекции.
В качестве начального приближения целесообразно вы. бирать равновесную функцию распределения. Вычислительная практика показывает, что такое начальное приближение обеспечивает сходимость про. цедуры обратных итераций к искомому минимальному по модулю собственному значению, т.е.
к константе скорости. Ядро уравнения (8.13) и матрица дискретного уравнения (8.29) не явля. ются симметричными, поэтому, прежде чем применять процедуру обратной итерации, необходимо симметризовать матрицу А. Это всегда возможно. так как ядро кинетического уравнения удовлетворяет принципу детального равновесия. Преобразование, симметризующее оператор, имеет вид л(е, г) = Р(е,г) К (е) . (8.32) Таким образом, вычисление нестационарной функции распределения сводится к решению линейной системы обыкновенных дифференциальных уравнений с постоянными коэффициентами, а вычисление квазистационарной функции распределения — к решению неполной проблемы собственных значений для матрицы коэффициентов этой системы дифференциальных уравнений.
При решении этих задач приходится иметь депо с плохо обуслов. ленной метрицей, разброс собственных значений которой составляет 19— 15 порядков. Заметим, что величина этого разброса мало зависит от способа дискретизации, так как определяется физикой процесса, т.е. константа скорости (минимальное по модулю собственное значение) отличается от других собственных значений обычно на несколько порядков. Рассмотрим процедуру численного решения нестационэрной задачи.
Пос. ле дискретизации кинетического уравнения задача сводится к решению линейной системы обыкновенных дифференциальных уравнений. Таким образом, задача сводится к вычислению матричной экспоненты от плохо обус- 197 ловленной матрицы большой размерности. Вследствие плохой обусловлен- ности матрицы А вычисление ехр (Ас) в виде ряда Тейлора невозможно.
Зта системы имеет большую размерность, поэтому использование функцио- нального уравнения Пойа ехр(А(г! +г!)) = ехр(Аг,) ехр(Аг!) (8.33) оказывается очень трудоемким. Для вычисления матричной экспоненты используется разложение в ряд производящей функции дпя полиномов Ла- герра [95[ иг — ! ,. l Аг! ! ехр Ат~ (Š— Аг!) "' 2, (.,".'(т) ( — 1)!~ ю:е Š— Аг! / (8.34) ИI (Е, Е ') = К ~ (Е) )г (Е', Е) К (Е ) . (8.35) Симметричную часть ядра, равновесную функцию распределения и констан- 198 где (.,'." (т) — полином Лагерра; л! = 20; Г! = гlгл. Как следует из принципа детального равновесия, матрица А обладает чисто действительным отрицательным спектром и вычисление по этой формуле можно производить для любого времени г с большой точностью.
Как видно иэ формулы (8.34), алгоритм вычисления решения уравнения (8.29) сводится к последовательности операций перемножения матриц Аг! и (Š— Ат!) ' на некоторые вектора и сложения получившихся векторов. В случае, когда ядро интегродифференциального уравнения отличается от ядра сильных столкновений и задача занимает промежуточное положение между диффузионной моделью и моделью сильных столкновений, матрица А имеет, как правило, ленточную структуру. Матрица (Š— Аг,) ' не обладает ленточной структурой, поэтому использовалось факторизованное представление матрицы (Š— Ат, ) [137[, т.е.
разложение ее на произведение правых и левых ленточных треугольных матриц и вычисление обратных в виде соответствующего произведения обратных ленточных треугольных матриц. Надо заметить, что ширина правых и левых ленточных треугольных матриц такая же, как и соответствующие ширины исходной матрица А (то же относится и к ширине обратных ленточных матрица) . Однако при перемножении этих треугольных матриц получается полная матрица, поэтому используется процедура последовательного перемножения полученных обратных ленточных треугольных матриц на вектор и компактное хранение этих матриц в памяти ЭВМ. Таким образом, описанная процедура не требует дополнительного объема памяти для вычисления матричной экспоненты, что позволяет выбирать достаточно малый шаг при разбиении энергетического интервала, т.е.
при дискретизации задачи. С использованием описанной методики в работе [139[ был проведен численный эксперимент. В качестве модели бралась условная пятиатомная молекула (число внутренних степеней свободы з . 9) . Величина энергии диссоциации Ее составляла 100 ккалlмоль. Интервал энергии, на котором производилась дискретизация, выбирался равным 2Ез и делился на 200 одинаковых отрезков разбиения г! = 2ЕеДУ ()У= 200) .
Ядро уравнения (8.13) строилось следующим образом: записывалась симметризуемая часть ядра ту спонтанного распада задавали следующим образом: а (Š— Е' + Ы, Е' — !)! < Е < Е', И!(Е,Е) = о(Е' — Е+Ы, Е'(Е~Е'+Ь, 0 , Е<.Е' — Ь или Е> Е'+Ь, К(Е) ехр ( — ), (э+1)! 8 Т 0 , Е< Ео. А(Е) = Е-Ео Т( — ), Е> Ео, ')=10 . Ео (8.36) (8ЗУ) (8.38) ,/ А (Е', Е) г)Е' о примерно равен единице. Интегралы в уравнении (8.28) аппроксимировапись с помощью квадратурных формул трапеций де+ ! Ь / А (Ер Е) Р (Е) г)Е = — (Ар. оТо +Ар. о ! ! ро+ ! ) ° лв 2 (8.39) Тогда получаем систему из )У линейных дифференциальных уравнений !т !у 1р= шЬ ь А'р !Р! — !>Ьгр Х Ар ! — Ар Фр. (8.40) !'= ! )= ! Ненулевые элементы матрицы А этой системы имеют вид !у арр+!ыЬАрр+!ар !ршЬАр!рдррсоЬХАр)Ар )= ! Эта модель рассчитывалась в широком диапазоне значений ы при 1000 К.
Нестационарная функцив 1 (Е, т) — решение этой системы — использовалась для получения зависимости концентрации и от времени. Константа скорости мономолекулярной реакции вычислвлась по кинетической кривой методом наименьших квадратов. затем с помощью функции г (е, г) определялсв так называемый квазистационарный лоправочный фактор ~„(Е) для равновесной функции распределения: Го (Е) = Т(Е, т) ехр(Аг) (А — константа ско. рости мономолекулярной реакции).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
