Основы вычислительного теплообмена и гидродинамики - Аникеев А.А., Молчанов А.М., Янышев Д.С., страница 11

В этом случаепридется конструировать численные методы, позволяющие проводитьрасчеты с шагом большим, чем выбираемым по условию (6.29).Недостатки альтернативы 1 мы уже показали.Второй подход приемлем для только что рассмотренной задачи, когдаправая часть первого уравнения зависит только от первой функции, а праваячасть второго – только от второй. В общем случае в правой части каждогоуравнения могут быть все искомые функции.Системууравнений(6.23)иначальныхусловий(6.24)удобнопредставить в матричном виде:dy= Ay; y (0 ) = y0dt(6.31) y1  y где y =  2  - вектор искомых функции, A - матрица коэффициентов, y  Ny0 - начальное значение вектора y , N - число искомых функций.В рассмотренном примере матрица коэффициентов им диагональный вид950  −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 >96100, хотя в химической кинетике часто встречаются задачи с коэффициентомжесткости ≈ 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 01+10τ (E − Aτ )−1 1=  1+τ 011 + 10 6 τ 0(6.38)Тогда решение на n-ом шаге интегрирования получается равным: 1n (1 + τ )ny =1 1 + 106τ(n (6.39))Если взять τ = 1 ⋅10 −2 (это максимальное значение шага, использованногов предыдущем примере), то получим результаты, представленные в Таблица6.497Таблица 6.4 Сравнение численного решения неявным методом Эйлера саналитическимШаг τy1точноерешение1 ⋅10 −20,049787Время t = 3y1численноерешение0,05053(погрешность1.5%)y2точноерешениеy2численноерешение00Основным преимуществом, полученным в результате использованиянеявного метода Эйлера, является возможность совместно решать всюсистему с шагом интегрирования намного превосходящим шаг, полученнымиз критерия устойчивости (6.29).

При этом появляется возможностьориентироваться только на медленные процессы, проводя интегрирование сшагом, превышающим характерные времена быстрых процессов.Рассмотрим еще один пример, менее жесткий, в котором вторая функцияне столь быстро уходит в нуль.Пример 3. Пусть матрица коэффициентов и начальный вектор равны − 1 99 ;A =  0 − 100 Собственныеотрицательнуюзначения1y0 =  1матрицывещественнуюравнычасть,т.е.(6.40)(-1,-100).решениеОниимеютэкспоненциальноустойчиво. Коэффициент жесткости системы равен 100.Точное решение системы: 2e −t + 100e −100t y= e −100t98(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.7390.010.73575890 .1(погрешность0.5%)0.771(погрешность5%)7.9 ⋅10 −313.78 ⋅10 −443.86 ⋅10 −11Таким образом, мы получили неплохую точность для медленныхпроцессов.99Алгоритм РозенброкаРассмотренныечисленныесхемыимеютпервыйпорядокаппроксимации.Приведем еще один метод, основанный на популярном алгоритмеРозенброка, реализованного в ряде математических пакетов.

Он основан наприведении системы (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.192100y2точноерешение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 попеременно «выскакивает» в отрицательную область,что, конечно, вряд ли годится для серьезных задач.1016.4Решениежёсткихсистемприменительнокзадачам химической кинетикиВ предыдущем параграфе было введено понятие жесткости системлинейных обыкновенных дифференциальных уравнений (ОДУ) и многовнимания было уделено способам их решения.Спрашивается, какое это имеет отношение к расчету химическиреагирующих течений? Самое прямое.

Дело в том, что система уравненийнеразрывности химических компонентов (6.9), во-первых, также являетсяжесткой, а, во-вторых, для ее решения применимы все те идеи и методы,которые были изложены в п. 6.3.Рассмотрим для начала систему (6.9) без учета конвекции и диффузии:d ( ρYI )= SIdt(6.44)В векторной форме это уравнение можно представить в виде:dY= f (Y )dt Y1 Y2где Y =  YN C - вектор массовых долей компонентов, S1 / ρ  S2 / ρ f = - вектор источников SN / ρ  C102(6.45)К уравнению (6.45) можно применить те же методы, что и к линейнымОДУ, только вместо матрицыпостоянных коэффициентов необходимоиспользовать матрицу Якоби: ∂f1∂Y∂f  1= A=∂Y  ∂f Nc∂Y 1∂f1∂YN∂f Nc∂YN…(6.46)Неявный метод Эйлера применительно к системе (6.45) имеет вид:y n +1 − y nτ(= f n + A y n +1 − y n)(6.47)Откуда(E − τA ) yn +1n− ynτ= fn(6.48)А метод Розенброка:(E − ατA − βτ2)A2 ⋅y n +1 − y nτ(( ))= f y n + γ ⋅τ ⋅ f y n(6.49)Решение (6.49) основано на следующих действиях, выполняемых накаждом n-ом шаге интегрирования:1.Вычисляется матрица производных (6.48) в точке y nСледующая точка y n +1 находится из матричного уравнения(6.49) с коэффициентами (6.43)Приведем классический пример – cистему Робертсона.

Рассмотрим2.систему трех уравнений:103 dy1 dt = −a1 y1 + a 2 y 2 y3 , dy 2= a1 y1 − a2 y 2 y3 − a3 y 2 , dt dy3 dt = a3 y 2Начальные условия: y1 (0) = 1,y 2 (0 ) = 0,(6.50)y3 (0) = 0Система (6.50) представляет модель химического взаимодействия трехвеществ: вещество «1»медленно превращается в «2»:«1» → «2» (соa1 = 0.1 ),скоростьювещество «2» превращается очень быстро в вещество «3»:«2» → «3»( a3 = 103 ).И, наконец, вещество «2» при каталитическом воздействии вещества «3»,превращается в вещество «1» ( a 2 = 10 2 ) : «2» + «3» → «1» + «3»Используем метод Розенброка.Матрица Якоби имеет вид: − a1A =  a1 0Результаты расчета представлены на рисунке104a2 y3− a 2 y 3 − a3a3a2 y 2 − a2 y2 0 Рисунок 6.1 Результат решения системы (6.50)6.5Методрасщеплениядлясистемыуравненийпереноса химических компонентовОсновную систему (6.9) можно представить в векторной форме∂Y= L(Y ) + f (Y )∂t(6.52)где L(Y ) - оператор, учитывающий диффузию и конвекцию:L(Y ) =∂∂x j ΓIeff ∂Y∂x j ∂ ( ρu j Y )−∂x j(6.53)f (Y ) - вектор источников (см.

п. 6.4). В системе (6.52) лучше отдельнорассматривать процессы переноса (конвекция и диффузия) и процессобразования компонента в результате химических реакций. Поэтому для еерешения имеет смысл применить метод расщепления по физическим105процессам. Он состоит в следующем. На каждом шаге по времени вместосистемы уравнений (54) решается последовательность уравнений:∂Y1= L(Y1 )∂t(6.54)∂Y2= f (Y2 )∂t(6.55)при выполнении условий:( ) ( )( ) ( )( ) ( )Y1 t n = Y t n ,nn +1,Y2 t = Y1 t n +1n +1Yt= Y2 t(6.56)В работе [24] показано, что для сходимости такой системы (т.е.аппроксимации и устойчивости) достаточно, чтобы сходилась каждая изсистем (6.54) и (6.55).Схема расщепления позволяет существенно сократить вычислительныересурсы.В каждое уравнение системы (6.54) входит только один химическийкомпонент; таким образом, можно решать каждое уравнение отдельно.