Глава 5. Особенности расчета химически реагирующих потоков (1013339), страница 2
Текст из файла (страница 2)
По другому определению свободный радикал — видмолекулы или атома, способный к независимому существованию (то естьобладающий относительной стабильностью) и имеющий один или дванеспаренных электрона.Как было ужереакцийτ chemуказано, характерные времена протекания химическихдля цепных реакций и для реакций рекомбинации-диссоциации очень сильно отличаются друг от друга. Т.е. справедливо:Dachain > Darecomb(6.22)Dachain - число Дамкелера для цепных реакций, Darecomb- числоДамкелера для реакций рекомбинации-диссоциации.Учитывать необходимо и те и другие реакции. Это порождает серьезнуюматематическую проблему, возникающую при решении системы уравненийпереноса (6.9), т.к. эта система является жесткой.В следующем пункте мы рассмотрим понятие жесткости системуравнений.6.3Жёсткие системыРассмотрим пример (пример 1). Пусть имеется система уравнений dy1 dt = λ1 y1 dy 2 = λ2 y 2 dt(6.23)с начальными условиямиy1 (0) = y 2 (0) = 1(6.24)Зададим значения числовых коэффициентов равными: λ1 = −1 и λ2 = −10 6Система легко решается аналитически, и ее решение имеет вид:y1 (t ) = e −t ;y 2 (t ) = e −106t(6.25)Нам интересно проверить решение системы численными методами,поэтому используем для решения системы (6.23) явный метод Эйлера.Из теории численных методом известно, что для сходимости решенияконечно-разностных уравнений к точному решению необходима такназываемая устойчивость разностной схемы.Метод Эйлера применительно к системе (6.23) дает следующеечисленное решение:yin +1 − yinτyin +1=yin= λyin ;+ τλi yin(6.26)=yin(1 + τλi )гдеi = (1, 2) , τ - шаг конечно-разностной сетки,yin - значение yi на n - ом шаге по t ,yin +1 - значение yi на (n + 1) - ом шаге по t .Конечно-разностная сетка задается соотношением:t n = nτ(6.27)y1i = (1 + τλi ), yi2 = (1 + τλi )2 ,....., yin = (1 + τλi )n(6.28)Таким образом,Условием устойчивости разностной схемы (6.26) является:λ1τ + 1 < 1;λ2τ + 1 < 1(6.29)т.е.
τ < 2 ⋅10 −6Предположим, что нам необходимо определить решение системы при6t = 3 . Точное решение (25) равно: y1 (3) = 0,049787 , y 2 (3) = e −3⋅10 . Значениеy 2 выходит за рамки разрядности любой вычислительной машины ифактически равно нулю.Решение конечно-разностного уравнения (6.26) с учетом (6.27) и (6.28)имеет вид:yi (t ) = (1 + λiτ )(t / τ )(6.30)(При шаге τ = 5 ⋅10 −7 и t = 3 получаем: y1 (3) = 1 − 10 −6)(3 / 10 −6)= 0.04978 ,y 2 (3) = 0Для более раннего момента времени t = 10 −4 при том же шаге τ = 5 ⋅10 −7( )y (10 ) = 3.72 ⋅10численное решение - y (10 ) = 0,9999 , y (10 ) = 6.22 ⋅10получаем:точноерешение-y1 10 −4 = 0,9999 ,−41−4−4−442;−612Мы видим, что для y1 выбранный шаг дает хорошее совпадение дляобоих моментов времени; для y 2 относительное совпадение плохое, носамом деле это не имеет большего значения, т.к.
в обоих случаях решениефактически равно нулю.Для получения более точного численного решения второго уравнениянеобходимо еще больше уменьшить шаг интегрирования. Например, еслиτ = 1 ⋅10 −9 , то численное решение при t = 10 −4 гораздо ближе к точному( )решению: y 2 10 −4 = 3.54 ⋅10 −44 .С другой стороны для получения хорошей точности при решениипервого уравнения совсем необязательно задавать такой маленький шаг какτ = 5 ⋅10 −7 . При шаге в 20000 раз больше, т.е.
при τ = 1 ⋅10 −2численное решение:y1 (3) = 0,049 ,получаем( )y1 10 −4 = 0,9999 - очень близкое кточному решению.Систематизируем полученные данные в Таблица 6.2 и Таблица 6.3Таблица 6.2 Сравнение точного и приближенного решения при различныхшагах по времениВремя t = 3Шаг τy1y1y2y2точноечисленноеточноечисленноерешениерешениерешениерешение5 ⋅10 −70,04978700,049041 ⋅10 −20,049787(погрешность0решениеневозможно1.5%)1 ⋅10 −90,0497870Таблица 6.3 Сравнение точного и приближенного решения при различныхшагах по времениШаг τy1точноерешение5 ⋅10 −71 ⋅10 −21 ⋅10 −90,9999Время t = 1 ⋅10 −4y1y2численноеточноерешениерешение0,99990,99990,99993.72 ⋅10 −44y2численноерешение6.22 ⋅10 −61решениеневозможно3.54 ⋅10 −44Основной вывод, который можно сделать, исходя из полученныхрезультатов, состоит в том, что шаги τ оптимальные для первого и второгоуравнения отличаются друг от друга на несколько порядков.
Это приводит кбольшим неприятностям. Решать первоеуравнение с шагом 1 ⋅10 −9непозволительная роскошь с точки зрения расходования вычислительныхресурсов, а использовать шаг 1 ⋅10 −2 для второго уравнения невозможно сточки зрения устойчивости.Вообще говоря, трудности численного решения подобных систем,получивших название жестких, связаны с выбором шага интегрирования.Дело в том, что характерныевремена исследуемых процессов могутразличаться в 10 6 раз. Следовательно, если при численном решении системынеобходимо выбирать шаг по самому быстрому процессу. В данном случаезатраты машинного времени для исследования самых медленных процессовбудут неоправданно велики.По этой причине имеются следующие альтернативы в выборе подхода кчисленному решению рассматриваемых задач.1. Численно решать систему с шагом, выбранным по условию (6.29), т.е.с учетом характерных времен всех процессов, описываемыхданнойсистемой.2.
Решать систему ОДУ с различными шагами, соответствующимифизическимпроцессамссущественноразличнымихарактернымивременами. В этом случае необходимо задавать условия перехода к другомушагу интегрирования.3.«Пренебречь»рассматриватьлишьбыстропротекающимимедленные,проводяпроцессамииинтегрированиечисленносшагом,превышающим характерные времена быстрых процессов. В этом случаепридется конструировать численные методы, позволяющие проводитьрасчеты с шагом большим, чем выбираемым по условию (6.29).Недостатки альтернативы 1 мы уже показали.Второй подход приемлем для только что рассмотренной задачи, когдаправая часть первого уравнения зависит только от первой функции, а праваячасть второго – только от второй. В общем случае в правой части каждогоуравнения могут быть все искомые функции.Системууравнений(6.23)иначальныхусловий(6.24)удобнопредставить в матричном виде:dy= Ay; y (0 ) = y0dt(6.31) y1 y где y = 2 - вектор искомых функции, A - матрица коэффициентов, y Ny0 - начальное значение вектора y , N - число искомых функций.В рассмотренном примере матрица коэффициентов им диагональный вид0 −1A = 6; 0 − 10 1 y0 = 1 (6.32)и каждое уравнение системы можно решать отдельно.В общем случае необходимо решать всю систему совместно.
Поэтомунампредставляетсянаиболееоптимальнымиспользоватьподход,основанный на альтернативе 3.Но сначала давайте подробнее исследуем понятие жесткости системобыкновенных дифференциальных уравнений.Описаннаявышеситуациявозникаетиз-засобственных значений матрицы системы (6.31):большогоразбросаλ2= 10 6 . Компонента сλ1бóльшим (по модулю) собственным значением вынуждает выбирать мелкийшаг и, одновременно, быстро перестает влиять на решение.
Классдифференциальных уравнений с таким поведением выделяется в теориичисленных методов понятием жестких уравнений.Точнее, система линейных автономных дифференциальных уравнений(6.31) называется жесткой, если, во-первых, все собственные значения λiматрицы A имеют отрицательную вещественную часть (т. е. система (6.31)экспоненциально устойчива).Re(λi ) < 0, (i = 1,2...N )(6.33)max Re(λi )S = 1<i < N>1min Re(λi )(6.34)и, во-вторых,1< i < NЧисло S при этом называют коэффициентом жесткости системы (6.31).Значок >> ("значительно превосходит") на практике обычно означает, что S >100, хотя в химической кинетике часто встречаются задачи с коэффициентомжесткости ≈ 10 6 и более.Более подробно с понятиями жёсткости и устойчивости системдифференциальных уравнений можно ознакомиться в [6, 23].Попробуем применить к системе (6.31) неявный метод Эйлера.y n +1 − y n= Ay n +1(6.35)y n +1 = (E − Aτ )-1 y n(6.36)τт.е.где E - единичная матрицаОтсюда получаем:()2()ny1 = (E − Aτ )-1 y0 , y 2 = (E − Aτ )-1 y0 , ..., y n = (E − Aτ )-1 y0Текст выводаПример 2.
Для системы (6.23), рассмотренной в предыдущем примере,0 1 + τматрица E − Aτ = 6 , а обратная ей: 0 1 + 10 τ (E − Aτ )−1 1= 1+τ 016 1 + 10 τ 0(6.38)Тогда решение на n-ом шаге интегрирования получается равным: 1n (1 + τ )ny =1 1 + 106τ(n (6.39))Если взять τ = 1 ⋅10 −2 (это максимальное значение шага, использованногов предыдущем примере), то получим результаты, представленные в Таблица6.4Таблица 6.4 Сравнение численного решения неявным методом Эйлера саналитическимШаг τy1точноерешение1 ⋅10 −20,049787Время t = 3y1численноерешение0,05053(погрешность1.5%)y2точноерешениеy2численноерешение00Основным преимуществом, полученным в результате использованиянеявного метода Эйлера, является возможность совместно решать всюсистему с шагом интегрирования намного превосходящим шаг, полученнымиз критерия устойчивости (6.29).
При этом появляется возможностьориентироваться только на медленные процессы, проводя интегрирование сшагом, превышающим характерные времена быстрых процессов.Рассмотрим еще один пример, менее жесткий, в котором вторая функцияне столь быстро уходит в нуль.Пример 3. Пусть матрица коэффициентов и начальный вектор равны − 1 99 ;A = 0 − 100 Собственныеотрицательнуюзначения1y0 = 1матрицывещественнуюравнычасть,т.е.(6.40)(-1,-100).решениеОниимеютэкспоненциальноустойчиво. Коэффициент жесткости системы равен 100.Точное решение системы: 2e −t + 100e −100t y= e −100t(6.41)Результаты расчета для 2-х моментов времени представлены в таблицах5,6Таблица 6.5 Сравнение численного решения неявным методом Эйлера саналитическимВремя t = 0.1Шаг τy1y1y2y2точноерешениечисленноерешениеточноерешениечисленноерешение1.810.011.8142150 .1(погрешность0.25%)1.73(погрешность5%)9.76 ⋅10 −44.54 ⋅10 −59.1 ⋅10 −2Таблица 6.6 Сравнение численного решения неявным методом Эйлера саналитическимВремя t = 1Шаг τy1y1y2y2точноерешениечисленноерешениеточноерешениечисленноерешение0.739(погрешность0.5%)0.771(погрешность5%)0.010.73575890 .17.9 ⋅10 −313.78 ⋅10 −443.86 ⋅10 −11Таким образом, мы получили неплохую точность для медленныхпроцессов.Алгоритм РозенброкаРассмотренныечисленныесхемыимеютпервыйпорядокаппроксимации.Приведем еще один метод, основанный на популярном алгоритмеРозенброка, реализованного в ряде математических пакетов.
Он основан наприведении системы (6.31)dy= Aydtк разностной схеме(E − ατA − βτ2)A2 ⋅y n +1 − y nτ(= A ⋅ y n + γ ⋅τ ⋅ A ⋅ y n)(6.42)Числовые коэффициенты подбираются таким образом, чтобы обеспечитьмаксимально возможный порядок точности. Они равны:α = 1.077, β = −0.372, γ = −0.577(6.43)Применим этот метод к примеру 3 и получим результаты для 3-моментоввремени:Таблица 6.7 Сравнение численного решения по алгоритму Розенброка саналитическимШаг τy1точноерешение0.010 .11.814215Время t = 0.1y1численноерешение1.80962.192y2точноерешениеy2численноерешение4.54 ⋅10 −53.3 ⋅10 −5-0.382Таблица 6.8 Сравнение численного решения по алгоритму Розенброка саналитическимШаг τ0.010 .1y1точноерешение0.7357589Время t = 1y1численноерешение0.735760.7357y2точноерешение3.78 ⋅10 −44y2численноерешение1.4 ⋅10 −456.7 ⋅10 −5Таблица 6.9 Сравнение численного решения по алгоритму Розенброка саналитическимШаг τ0.010 .1y1точноерешение9.9574 ⋅10 -2Время t = 3y1численноерешение9.9574 ⋅10 -29.9561⋅10 -2y2точноерешение0y2численноерешение03 ⋅10 −13При шаге τ = 0.01 мы получили хорошие результаты как для медленноменяющейся функции y1 , так и для быстроменяющейся функции y 2 .Для большего шага τ = 0.1 получено удовлетворительно согласованиедля y1 , а функция y 2 попеременно «выскакивает» в отрицательную область,что, конечно, вряд ли годится для серьезных задач.6.4Решениежёсткихсистемприменительнокзадачам химической кинетикиВ предыдущем параграфе было введено понятие жесткости системлинейных обыкновенных дифференциальных уравнений (ОДУ) и многовнимания было уделено способам их решения.Спрашивается, какое это имеет отношение к расчету химическиреагирующих течений? Самое прямое.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
