Основы численных методов для газодинамики с теплопроводностью (стр. 6 )
 | Из за большого объема этот материал размещен на нескольких страницах: 1 2 3 4 5 6 7 8 |
Построенная схема в принципе удобна для расчетов, но обладает одним специфическим свойством, на которое, повидимому, впервые указал . Значения сеточных функций образуют две группы — первая содержит значения u с нечетным нижним индексом и значения p с четным, вторая, наоборот, — с четными индексами у u и нечетными у p. Легко убедиться, что эти группы считаются независимо и не влияют друг на друга. Мы тем самым ищем два самостоятельных решения нашей задачи. Возникают вопросы, нужны ли нам два решения, что мы получаем ценой удвоения вычислительной работы, какое из решений лучше и так далее. Более сложным является вопрос, который мы не хотим подробно обсуждать, но которого не избежать — сколько граничных условий нужно задавать при такой постановке.
Все эти проблемы, а, возможно, и соображения физического порядка привели к выбору схемы, у которой давление задается в полуцелых точках, что по смыслу равносильно выбору более густой сетки и вычислению на ней одного из вышеназванных решений. Схема имеет вид:
Аппроксимация определяется следующими выражениями:
Погрешность аппроксимации имеет первый порядок по времени и второй по пространству. Это связано с тем, что по времени мы используем односторонние разности, а по пространству — центральные. Естественно, возникает желание повысить порядок аппроксимации, сместив одну из функций на полшага по времени. Удобнее это сделать для скорости. В результате получается знаменитая схема «КРЕСТ», сыгравшая огромную роль в проведении одномерных, а затем и двумерных расчетов газодинамики с теплопроводностью, во ВНИИЭФ эта схема связана прежде всего с именем . Схема имеет вид:
Название схемы связано с расположением четырех используемых в ней точек сетки на плоскости x, t. Аппроксимация определяется выражениями:
Порядок аппроксимации действительно увеличился до второго.
Любопытно отметить, что две последние схемы отличаются лишь обозначениями (использование полуцелого индекса вместо целого). Если по ним решать одну и ту же разностную задачу, то результаты совпадут полностью. Как при этом трактовать тот факт, что схемы имеют разный порядок точности? Можно, конечно, ссылаться на то, что при сравнении разностного решения с точным само точное решение следует брать в несколько отличающиеся моменты времени. Поэтому, если решение для одной схемы близко к точному в момент t, то для другой его нужно сравнивать с точным решением на момент t + 0.5 τ, которое может отличаться как раз на величину O(τ). Однако, на деле нам не известно точное решение и мы не занимаемся такими сравнениями. Нас обычно интересуют интегральные величины, такие как скорость разлета осколочной оболочки в задачах метания, максимальная средняя плотность одной из областей в задачах обжатия, время от начала обжатия до момента достижения критической плотности и тому подобное. Все эти результаты в обоих разностных решениях не отличаются. Поэтому к оценкам порядка аппроксимации нужно подходить с некоторой долей скепсиса, понимая, что они могут служить лишь некоторым ориентиром, но не решающим критерием. В дальнейшем мы будем рассматривать в основном первую из этих двух схем, назовем ее «квазиКрестом».
Заметим, что именно эта схема использована в таких комплексах, как СИГМА, ЧАС, ТЕМП. Это связано с тем, что ее удобнее сочетать со счетом других процессов — теплопроводности, кинетики ядерных реакций и т. п.
2.2.2. Устойчивость.
Разностную схему называют неустойчивой, если ошибки растут в ней настолько быстро, что сильно искажают решение. На практике попытка считать по неустойчивой схеме чаще всего приводит после некоторого числа шагов к переполнению. Источник ошибок в данном случае не рассматривается, всегда, в частности, в счете присутствуют погрешности округления. При таком подходе неустойчивость является внутренним свойством разностной схемы и не зависит от того, какое дифференциальное уравнение эта схема аппроксимирует. Устойчивой считается схема, у которой раз допущенная ошибка растет не слишком быстро, или, что чаще стараются обеспечить на практике, не растет. Более точное определение мы дадим позже.
При исследовании устойчивости обычно считают погрешнось малой, что позволяет, пренебрегая малыми второго порядка, линеаризовать схему и тем самым рассматривать линейные однородные преобразования. Дело, таким образом, сводится к исследованию норм или спектра линейных операторов. Подробнее мы все это рассмотрим на примерах.
2.2.2.1. Обыкновенное дифференциальное уравнение
Рассмотрим прежде всего явную разностную схему:
Коэффициент c будем считать положительным. При отрицательном коэффициенте решение дифференциального уравнения неограниченно растет и анализ роста погрешности в сравнении с ростом решения требует учета дополнительных деталей, на которых нам не хотелось бы останавливаться. При положительном c решение дифференциального уравнения ограничено и того же естественно ожидать от разностного решения.
Допустим, что на некотором шаге возникла погрешность (например, округления) Δun. Тогда на следующем шаге мы будем иметь:
Мы получили линейное однородное уравнение, из которого следует
Погрешность, таким образом, при переходе на следующий шаг умножается на множитель q, она будет расти, то есть схема будет неустойчивой, если |q| > 1. Напротив, схема устойчива, если |q| ≤ 1. В этом случае единичная погрешность не растет, или даже убывает. Если на каждом шаге допускаются новые ошибки, то каждая из них (в линеаризованной постановке) будет изменяться независимо, а суммарная ошибка, как мы увидим ниже при анализе сходимости, остается приемлемой. Условие устойчивости явной схемы, таким образом, имеет вид:
Счет по явной схеме можно вести только при достаточно малом шаге по времени. Насколько такой шаг годится для конкретной задачи, зависит от масштаба разрешения, выбранного в этой задаче, к этому вопросу мы вернемся при обсуждении коэффициента потерь.
Несколько примеров поведения погрешности представлено на следующих рисунках. Возьмем для определенности значения c = 1, Δu0 = 0.001, тогда устойчивость будет зависеть от шага τ.
Рис. Явная схема, τ = 0.5.
Рис. Явная схема, τ = 1.8.
Видно, что погрешность убывает в обоих случаях, но во втором случае — немонотонно. Дальнейший рост τ приводит к неустойчивости.
Рис. Явная схема, τ = 2.6.
Погрешность в данном случае растет по модулю на каждом шаге, что вскоре должно привести к переполнению. Для сложных задач исследование устойчивости в аналитическом виде удается провести далеко не всегда. В счете неустойчивость проявляет себя чаще всего подобными осцилляциями.
Для неявной схемы
Очевидно, что он всегда меньше единицы. Подобные схемы называются безусловно устойчивыми. По неявной схеме можно считать с любым шагом — авостов не будет. Какая при этом будет точность решения, зависит от самого решения, погрешности аппроксимации и т. п., но уже не от устойчивости.
Для симметричной схемы
множитель перехода имеет вид:
Очевидно, эта схема также безусловно устойчива. В определенных условиях схема имеет более высокую точность, чем две предыдущие, к этому анализу мы еще вернемся. Следует, впрочем, отметить что при больших τ погрешность убывает достаточно медленно и при этом осциллирует.
Рис. Симметричная схема, τ = 100.
2.2.2.2. Уравнение теплопроводности
Рассмотрение начнем с явной схемы:
Счет идет также шагами по времени, поэтому при исследовании устойчивости необходимо проследить за изменением погрешности при переходе на очередной шаг. Пусть на n-ом шаге погрешность задана во всех пространственных точках и равна ∆ukn. Тогда, как легко проверить, на следующем шаге погрешность будет равна
Если 
≤ 0.5, то все коэффициенты в правой части положительны и их сумма равна единице.
Подобные уравнения обладают интересной особенностью: они удовлетворяют так называемому принципу максимума. Предположим, что все компоненты погрешности на n-ом шаге не превосходят некоторого числа M, тогда вся правая часть не превосходит этого числа. Следовательно, погрешность на n+1-ом шаге также не превзойдет этого числа: ∆ukn+1 1 могут быть отрицательными, что нефизично и приводит, в частности, к нарушению принципа максимума. Это является еще одним примером того, что вопрос о точности не сводится только к порядку аппроксимации.
2.2.2.3. Прогонка.
Условие устойчивости явной схемы имеет вид t/h2≤0.5. Много это или мало? Все зависит от задачи. Мы рассматривали упрощенное уравнение, реальное уравнение теплопроводности включает в своих членах коэффициенты теплоемкости Cv, теплопроводности κ, а также плотность ρ. Вместе с шагами по времени и пространству они образуют безразмерную комбинацию, которую называют числом Куранта, или, сокращенно, курантом:
Условие устойчивости явной схемы при этом имеет вид Cur ≤ 0.5. Число Куранта характеризует физическую задачу, как таковую, независимо от разностной схемы. Под величинами τ и h при этом понимаются масштабы разрешимости задачи. В «обычных» тепловых задачах тепло передается посредством диффузии (в газах) или передачи энергии колебаний атомов от одной части тела к другой (в твердых телах). Оба процесса протекают достаточно медленно, число Куранта невелико и применение явной схемы обычно не вызывает проблем. Иначе обстоит дело в задачах ядерного взрыва. Перенос энергии осуществляется квантами света (точнее, квантами рентгеновского диапазона и γ — лучами). Коэффициент теплопроводности κ при этом достаточно велик и курант зачастую достигает значений 105 и больше. Применение неявной или какой-нибудь другой устойчивой схемы в этих условиях становится неизбежным. Ситуация осложняется еще тем, что в различных частях задачи в один и тот же момент времени плотность, температура и зависящий от них коэффициент теплопроводности могут сильно отличаться, так что курант принимает значения от величин порядка единицы и меньше до сотен тысяч. Разностная схема должна обеспечить точный счет во всем этом диапазоне чисел Куранта. Это накладывает дополнительные требования на разностные схемы. Подобная проблема позднее была изучена применительно к так называемым жестким системам обыкновенных дифференциальных уравнений (см., например, [[23]]). Тем интереснее, что неявная схема изначально обладает, как мы позже увидим, этими свойствами. Это и определило ее успешное применение в задачах ядерной физике еще в 50-х годах.
Реализация неявной разностной схемы требует на каждом шаге по времени решения системы уравнений. Собственно, это решение и называется прогонкой. Мы рассмотрим систему уравнений более общего вида, чем та, которая порождается простейшей схемой. Тем самым мы учтем случай переменных коэффициентов, который типичен для реальных задач. Для простоты обозначений мы не будем указывать верхний индекс у искомых величин, тогда система уравнений может быть записана в следующей форме:
Свободный член здесь ненулевой, в него входят значение u с нижнего слоя по времени, а также члены от правой части уравнения теплопроводности, которые в реальных задачах обычно присутствуют. Кроме того, должны быть заданы граничные условия. Мы будем считать, что на обеих границах задана температура, то есть известны u0 и un+1.Таким образом, система содержит n уравнений с n неизвестными. Уравнения такого рода мы будем называть трехчленными.
Первый этап прогонки состоит в том, чтобы заменить систему трехчленных уравнений двучленными. Последние удобно брать в виде:
Заметим, что первое из трехчленных уравнений фактически является двучленным, поскольку u0 известно. Остается только привести его к нужному виду, то есть определить коэффициенты P и Q:
После этого остальные двучленные уравнения строятся рекуррентно в порядке возрастания индекса (как говорят, слева направо). Пусть двучленное уравнение известно для индекса k-1. Тогда с его помощью можно исключить uk-1 из трехчленного уравнения:
Для того, чтобы привести его к нужному виду, остается выразить коэффициенты:
Собственно, рекуррентное вычисление коэффициентов по этим формулам и составляет первый этап, или, как его называют, прямую прогонку.
Последнее двучленное уравнение соответствует индексу k = n. Поскольку un+1 нам известно, мы можем из этого уравнения найти un. Двигаясь дальше в направлении убывания индекса (справа налево) мы из двучленных уравнений определим все неизвестные. Этот ход называют обратной прогонкой.
Казалось бы, что разработку столь простого алгоритма можно рассматривать, как текущий эпизод в программировании. Тем не менее, в 60-х годах разгорелся нешуточный спор за приоритет в создании прогонки. Мне пришлось быть свидетелем этого. Мы подготовили статью [[24]], содержащую ссылку на один из авторских коллективов [[25]]. Когда я попросил подписать представление статьи в ДАН, он сказал, что, как физик, не может рекомендовать статью, в которой приоритет создания прогонки отдается математикам. На это я ответил, что лично не знаю историю вопроса и предложил сослаться на оба коллектива (включая [[26]]). Это устроило Я. Б. и статья была рекомендована. Только позднее я понял, что вопросы приоритета связывались, конечно, не с самим алгоритмом (являвшимся, кстати, частным случаем метода Гаусса решения систем линейных уравнений), а с понятиями устойчивости, разработки неявной схемы как способа обеспечения устойчивости и другими элементами, заложившими фундаментальные основы теории разностных схем.
Перечитывая этот фрагмент, я заметил, что по привычке допустил краткое упоминание «Я. Б.». Мне вспомнилось, что на инструктаже при поступлении на работу нас предупреждали, что в разговорах нельзя называть фамилии ведущих ученых и рекомендовали пользоваться сокращениями «ЯБ», «АДС». Нужно сказать, что это оказалось удобным и мы охотно пользовались этими сокращениями даже тогда, когда не было посторонних. Разумеется, в этом не было ни тени фамильярности, поскольку авторитет этих людей был и остается необычайно высоким. Впрочем, Николая Александровича Дмитриева почему-то называли по имени-отчеству или фамилии. Возможно, в режиме не знали, каким авторитетом он пользовался.
2.2.2.4. Уравнения газодинамики.
Для уравнений газодинамики подход к исследованию устойчивости основывается на тех же принципах, которые мы уже рассмотрели применительно к уравнению теплопроводности. Необходимо определить множитель шага, характеризующий рост погрешности при переходе от одного временного слоя к следующему. Различие в технике заключается в том, что раньше у нас на каждом слое был один вектор, теперь их будет два — погрешность скорости и погрешность давления. Разумеется, эти векторы можно объединить в вектор двойной размерности, но тогда возникнут трудности с определением собственных векторов. Оставим поэтому два вектора такими, как они есть. Множитель перехода λ будем искать в предположении, что
Здесь a и b — неопределенные коэффициенты, λ — искомый множитель шага. Рассмотрим явную разностную схему:
Напомним, что погрешности удовлетворяют той же системе уравнений. Подставляя выбранные выражения, получим:
Используя формулу Эйлера, эти уравнения можно преобразовать к виду:
Для того, чтобы полученная система двух однородных уравнений относительно неопределенных коэффициентов a и b имела ненулевое решение, необходимо, чтобы ее определитель равнялся нулю:
Отсюда (λ -1) 2 + (τ/h) 2 sin2ω = 0, или λ2 -2 λ + 1 + (τ/h) 2 sin2ω = 0 . Значение ω попрежнему следует рассматривать в диапазоне [0, π]. Среди этих значений найдутся такие, например, π/2, что свободный член уравнения будет больше единицы при любых шагах по времени и пространству. Но, согласно формулам Виета, свободный член приведенного квадратного трехчлена равен произведению корней этого трехчлена. Отсюда следует, что по крайней мере один из корней по модулю больше единицы. Схема, таким образом, безусловно неустойчива. Ею, тем не менее, иногда пользуются, например, в методе частиц [[27]]. Рост погрешности при этом в области интенсивного течения подавляется конвективными членами и проявляется только в застойных зонах.
Анализ других схем для уравнений газодинамики производится аналогично. Необходимо, прежде всего, сообразить, куда в определителе следует поставить λ, обозначающую, что этот член берется с верхнего слоя. Для неявной схемы — это оба члена, соотвествующих пространственным производным. Определитель при этом выглядит так:
Уравнение приобретает вид (λ -1) 2+ λ2 (τ/h) 2 sin2ω = 0. Заметим сразу, что это уравнение не имеет вещественных корней, так как оба слагаемых в левой части неотрицательны и одновременно в ноль не обращаются. Комплексные корни квадратного уравнения с вещественными коэффициентами сопряжены и, следовательно, равны по модулю. Согласно формулам Виета их произведение, то есть квадрат модуля равен отношению свободного члена (здесь 1) к коэффициенту при старшем члене: 1+ (τ/h) 2 sin2ω. Это отношение при всех значениях параметров не превосходит единицы, следовательно, неявная схема газодинамики безусловно устойчива.
Рассмотрим теперь схему КРЕСТ, точнее квазиКРЕСТ:
Она отличается тем, что пространственная производная для давления в ней аппроксимируется на нижнем слое, а для скорости — на верхнем. Кроме того шаг аппроксимации равен h, а не 2h. В результате определитель для этой схемы выглядит следующим образом:
Соответственно, уравнение имеет вид
λ2 -τ/h) 2 sin2 (ω /2 )) λ + 1 = 0 .
Произведение корней этого уравнения (опять по формулам Виета) равно единице. Если корни комплексные, то каждый из них по модулю равен единице. То же относится к случаю кратных вещественных корней. В остальных случаях вещественные корни различны по модулю и, следовательно, один из них больше единицы. Таким образом, условие устойчивости схемы сводится к неположительности дискриминанта:
(1 — 2(τ/h) 2 sin2 (ω /2≤ 0.
Для этого требуется, чтобы
-1 ≤ 1 — 2(τ/h) 2 sin2 (ω /2 ) ≤ 1.
Правое неравенство выполняется всегда, а для выполнения левого при любых значениях ω необходимо, чтобы
Это неравенство представляет собой условие устойчивости схемы квазиКРЕСТ. Поскольку схема КРЕСТ отличается только обозначением индексов, что неважно с точки зрения устойчивости, то условие устойчивости для нее такое же. Отношение, стоящее в левой части неравенства, назывется числом Куранта для газодинамики. В случае реальных (неупрощенных) уравнений число Куранта содержит скорость звука C:
Это число является безразмерным и имеет смысл независимо от разностной схемы. Под τ и h при этом понимаются масштабы разрешимости задачи. В динамических задачах масштабы обычно согласованы так, что Cur ≈ 1, в результате условие устойчивости схемы КРЕСТ вполне приемлемо, что и определяет (вместе с ее простотой и точностью) исключительную значимость этой схемы.
2.2.2.5. Предиктор-корректор.
Разностная схема должна, как мы видели, удовлетворять одновременно нескольким различным требованиям. Иногда эти требования трудно совместить в рамках одной схемы. Типичный пример — счет теплопроводности в задачах ядерного взрыва. Схема должна быть, естественно, устойчивой и, кроме того, очень точно сохранять балланс тепла. В некоторых легких областях поток тепла через область за шаг по времени бывает сравним с количеством тепла в этой области. Ясно, что даже небольшая ошибка в баллансе может привести через несколько шагов к недопустимой погрешности. Выполнить подобные требования иногда легче поочередно. Счет шага разбивается на два полушага, один из них — предиктор — обеспечивает устойчивость, второй — корректор — баллансность. Подобный способ используется, в частности, в комплексе СИГМА. Возможно и другое распределение свойств между полушагами.
Впервые я столкнулся с этой идеей буквально в первые дни работы во ВНИИЭФ. Нужно было решить дифференциальное уравнение u'(t) = f(u) со сложной нелинейной правой частью. Явная схема
не годилась, а неявная была сложна в реализации. Мой начальник Юрий Михайлович Шустов (не могу не вспомнить с удовольствием время работы в его отделе) сразу четко произнес «Эйлер с пересчетом»:
Здесь промежуточное значение u* вычисляется для того, чтобы подставить его в правую часть и тем самым избежать решения нелинейного уравнения.
Рассмотрим еще в качестве примера квазинеявную схему для счета газодинамики. Раньше мы говорили, что масштабы разрешимости по времени и пространству в задачах газодинамики, как правило, согласованы, что позволяет с успехом пользоваться схемой КРЕСТ. Но бывают задачи, где в различных областях пространственный масштаб отличается. Типичной задачей такого рода является взрыв конденсированного ВВ или ядерного заряда в атмосфере. В самом заряде необходимо брать мелкую сетку, в воздухе вдали от заряда сетку можно брать значительно крупнее. На начальном этапе расчет следует вести с мелким шагом по времени, соответсвующим скорости процессов в заряде. На заключительном этапе временной масштаб можно брать более крупным, соответствующим движению фугасной волны в воздухе. Точность расчета процессов в самом заряде и его окрестности в этот момент мало существенна, так что шаг по времени может определяться сеткой в воздухе. Поскольку сетка в заряде остается мелкой, такой шаг требует применения абсолютно устойчивой схемы. Выше мы рассматривали неявную схему, однако, ее реализация требует прогонки двух величин — скорости и давления, что затруднительно, особенно в двумерном случае. Кроме того, возникают проблемы в случае нелинейного уравнения состояния. Более удобна двухэтапная схема, реализованная в одномерном аналоге комлекса СИГМА [[28]], а также в программах РОЗГА и ЧАС.
На первом этапе рассчитывается экстраполированное значение одной из величин (давления или скорости) на верхнем слое по времени. В двумерном случае удобнее брать давление, которое является скалярной величиной. Для простоты мы рассмотрим плоский одномерный случай, но именно с экстраполяцией давления. Первый этап обеспечивает безусловную устойчивость схемы, второй — консервативность и реализацию нелинейных уравнений состояния.
Будем исходить из дифференциальных уравнений:
На первом этапе возьмем линеаризованное уравнение для давления:
Экстраполированную скорость возьмем из разностного уравнения:
Подставим u* в уравнение для давления на место u:
Полученное уравнение можно заменить неявным относительно p разностным уравнением и решать его с помощью прогонки. Скорость при этом берется с нижнего слоя. При желании член, содержащий скорость, можно вообще опустить, на устойчивость это не влияет. В таком случае первый этап можно трактовать, как применение сглаживающей процедуры.