Введение в вычислительную физику - Федоренко Р.П.
ISBN 5-7417-0002-0
Скачать (прямая ссылка):
^OzO + ^OzI = ^fO
•^mZm-l + ^mzOi + ^oizO. + ! = #«. -s^MzM-1 + ^MzM = ^m-
Здесь использованы обозначения
(13)
I : J_ A0 Ail2
m m 38 = m m
О ^m+1/2 \ ’ m »-1/2 т+1/2 \ рО & m + 1/2 ^
z>l /2 z>l т + 1/2 т + 1/2
Sr..=
( А_
В,
m +1/2
Формулы вычисления элементов этих матриц через значения функций I/, E и их производных очевидны, но громоздки. Нет необходимости их воспроизводить.
Система уравнений (13) имеет «трехдиагональную» форму и решается несложным обобщением алгоритма прогонки. Вывод формул алгоритма отличается от вывода, изложенного в § 10, только тем, что теперь мы работаем с матрицами (некоммутативная алгебра) и надо аккуратно следить за порядком множителей. Решение ищется в форме
338
ПРИБЛИЖЕННЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ
[Ч. II
zm_, = XtnZm + Ym, где Xm — матрица 2x2, Ym- вектор. Опуская простые выкладки, приведем результат (т=1,2, ..., M — 1):
X1 = V0, Хт+1 = -(<*„ + Г„,
Yl = OfO1^f0- = + (&т YJ-
Теперь последнее уравнение (13) и прогоночное соотношение = XMzM + Ym можно разрешить относительно zM. Эта величина вычисляется и позволяет начать «обратную прогонку» — вычисление справа-налево искомых величин zm.
Мы не будем здесь обсуждать проблем разрешимости всех встречающихся в алгоритме задач (существования обратных матриц) и сходимости итерационного процесса. Укажем лишь, что легко угадать тривиальный результат: при достаточно малом т все обстоит благополучно. Это естественное следствие вырождения в пределе г-*0 всех уравнений в тривиальные. Теоретические оценки того малого т, начиная с которого гарантируется успех вычислений, в практике расчетов не используются — это привело бы к неоправданно заниженному шагу по времени. Однако сам факт зависимости, например, скорости сходимости итераций от т (она тем выше, чем меньше т) используется в режиме обратной связи. Считается, что требуемая точность должна достигаться за три-пять итераций.
Если итераций потребовалось больше, следующий шаг интегрирования уравнений выполняется с уменьшенным шагом т. Если точность достигается за меньшее число итераций, т увеличивается в пределах, определяемых другими критериями выбора шага. Что касается сопоставления скорости сходимости методов раздельной и векторной прогонок, то преимущество имеет последняя. Это естественно: оба метода являются комбинацией метода линеаризации (Ньютона) и метода простой итерации (Пикара).
Общая картина в таких алгоритмах такова, что метод оказывается тем быстрее сходящимся, чем больше в нем доля метода Ньютона. Однако итерация метода векторной прогонки требует больших вычислений. Вообще, следует отметить, что основное время работы ЭВМ связано не с прогонкой, а с вычислением коэффициентов системы (13). Одним из ресурсов экономии вычислительной работы является алгоритм с однократным вычислением этих коэффициентов; при этом в процессе итераций пересчитываются только Dm в (13).
Поясним это на промере решения нелинейного уравнения f(x) =0. Упрощенный вариант метода Ньютона с однократным вычислением / имеет форму
Xі + 1 = Xі - f Х(х°) /(Xі).
§22]
РЕАЛИЗАЦИЯ РАЗНОСТНОЙ СХЕМЫ
339
При достаточно хорошем начальном приближении х° он сходится «линейно», т.е. IIZ(X1)II Q1. где ||Е — fx(x) /J1(^0)II (х — Pe" шение системы). Эту оценку предоставим вывести читателю.
И наконец, подчеркнем, что приведенные выше формы организации решения уравнений на верхнем слое образуют некоторую общую схему, в рамках которой возможны различные варианты. Они появляются при различных способах отнесения тех или иных неизвестных к i-й итерации (по таким неизвестным проводится линеаризация) или к (і — 1)-й, а в иных случаях и к п-му слою. Выбор того или иного варианта диктуется особенностями решаемой задачи и здесь, естественно, не конкретизируется.
Расчет ударной волны по недивергентной гибридной схеме.
Сказанное выше может создать у читателя впечатление, что для правильного расчета ударных волн дивергентная форма разностных уравнений является существенным фактором. В общем это верно. He следует только возводить это положение в ранг абсолютного, безусловного требования к используемым в расчетах схемам. Обсудим этот вопрос подробнее, опираясь на результаты вычислительного эксперимента, проведенного автором в 1962 г.
Рассмотрим численное решение задачи о распаде произвольного разрыва в начальных данных. Кроме трудностей расчета ударной волны, мы имеем проблему расчета контактного разрыва, так как задача решается в переменных Эйлера. Итак, уравнения имеют вид
В качестве основных термодинамических величин берутся удельный объем V и величина с = Vpv (которая с точностью до множителя совпадает с адиабатической скоростью звука; рассматривается идеальный газ с 7 = 5/3). Вязкость берется в форме фон-Неймана:
q = (Lh)zux(ux — I их I), где L 4 -н 5.
Выбор «экзотической» переменной с объясняется просто. В зоне размазанной волны (см. § 20) переменная с ведет себя так же, как функции инь, тогда как р и е фактически размазываются на вдвое меньшую длину. Поскольку профиль волны должен быть разрешен четырьмя-пятью счетными точками, при счете в терминах е или р пришлось бы увеличить L раза в два, что приводит к слишком большому размазыванию и и v. При решении задач в лагранжевых координатах этой проблемы нет, так как переменные е и р по х не дифференцируются, а шаг по времени по разным причинам таков,
