Основы численных методов для газодинамики с теплопроводностью (стр. 7 )

После того, как найдено экстраполированное давление, дальнейшие вычисления можно вести, как в схеме КРЕСТ или квазиКРЕСТ. Эктраполированное давление используется при счете скорости, в конце шага давление пересчитывается с помощью уравнения состояния. Двухэтапная схема абсолютно устойчива и сохраняет основные достоинства схемы КРЕСТ (погрешность аппроксимации, естественно, несколько изменяется).

2.2.3. Сходимость.

При проведении расчетов с помощью разностной схемы возникают, как мы видели, погрешности аппроксимации и погрешности округлений. При исследовании устойчивости нас интересовал в основном рост последних, так как именно он определял, можно ли вообще считать по данной схеме. Теперь мы, напротив, будем пренебрегать округлениями, считая схему устойчивой, и сосредоточимся на погрешности аппроксимации. Напомним, что эта погрешность определяет, насколько разностное уравнение отличается от дифференциального уравнения. Нас реально интересует другой вопрос: насколько решение, найденное с помощью разностной схемы, отличается от решения дифференциального уравнения. Исследование можно разбить на несколько этапов:

1. Оценка ошибки после одного шага по времени,

2. Суммирование ошибок при счете до интересующего нас момента времени t,

3. Оценка того, как шаг сетки влияет на результирующую погрешность.

Рассмотрение будем вести на примере явной разностной схемы для обыкновенного дифференциального уравнения:

При исследовании аппроксимации мы выяснили, что при подстановке в это уравнение решения дифференциального уравнения возникает невязка

Соответственно, решения разностного и дифференциального уравнений выразятся следующим образом (последнее с точностью до малых более высокого порядка):

Если считать, что на нижнем слое решения совпадают, то разница на верхнем слое (погрешность одного шага) составит

Если считать, что схема устойчива (в рассматриваемом случае достаточно предположить, что τ < 2 / c), то эта погрешность не возрастая по модулю будет передаваться на следующие шаги по времени вплоть до шага N = t / τ. Учитывая, что такие погрешности возникают на каждом шаге по времени и что в силу линейности уравнения они ведут себя независимо и могут быть, таким образом, просуммированы на последний момент времени, получим оценку итоговой ошибки:

Отсюда видно, что погрешность, во-первых, мала, если мал шаг τ, и, во-вторых, она уменьшается при уменьшении этого шага, несмотря на то, что число шагов и число допущенных отдельных ошибок при этом растет. Последнее свойство и означает сходимость.

Мы воспроизвели на этом простом примере схему доказательства знаменитой теоремы Лакса:

Если разностная схема обладает аппроксимацией и устойчивостью, то ее решение сходится к искомому решению дифференциального уравнения при уменьшении шага сетки.

Разумеется, для уравнений с частными прозводными необходимо одновременное, в ряде случаев согласованное, уменьшение шагов сетки по пространству и времени. Само доказательство также усложняется, поскольку необходимо оперировать с векторами и их нормами. Но основные моменты рассуждений - оценка погрешности решения в отдельной точке, ограничение ее роста в предположении устойчивости и суммирование итоговой погрешности - остаются.

Теорема Лакса дает достаточные условия сходимости. Не всегда акцентируют внимание на том, что они не являются необходимыми. Между тем, существуют примеры, когда приходится пользоваться схемами, не обладающими аппроксимацией. Один из них - теплопроводность в комплексе СИГМА.

При совместном счете газодинамики и теплопроводности выбор разностной схемы зачастую определяется соображениями физического характера и вопросами интерпретации. Так, скорость удобно задавать в узлах сетки, тогда естественным образом рассчитываются координаты узлов и их движение. Плотность в этом случае разумно связывать с отрезком в одномерном случае или ячейкой сетки в двумерном. Тем самым для лагранжевой сетки естественным образом выполняется закон сохранения массы. С ячейкой целесообразно таким образом связывать и все остальные термодинамические величины, так как они должны выражаться друг через друга при счете уравнений состояния. Тем самым температуру необходимо связывать с ячейкой или, для определенности, с ее центром. В задачах ядерной физики, где потоки тепла очень велики, абсолютно необходимо обеспечить консервативный счет потоков, для чего их нужно задавать на границах ячеек. В этом случае, сколько тепла вытечет из одной ячейки, столько добавится в соседнюю. Тем самым разностная схема для теплопроводности определяется практически однозначно. Посмотрим, к чему это приводит на неравномерной сетке. Заметим, что даже равномерная сетка при счете газодинамики может стать неравномерной.

Рассмотрим для простоты плоскую одномерную задачу. Узлы сетки обозначим через xk. Тогда явную схему естественно задать в виде:

Здесь координата с дробным индексом относится к центру ячейки: xk+1/2 = 0.5(xk + xk+1). Обозначим для краткости h = xk+1 - xk, hb = xk+3/2 - xk+1/2, ha = xk+1/2 - xk-1/2. Для исследования аппроксимации разложим решение в окрестности точки (xk+1/2 , tn):

Подставляя эти выражения в разностное уравнение, получим

На равномерной сетке ha = hb = h и мы получаем обычную аппроксимацию. На неравномерной сетке коэффициент при второй производной, вообще говоря, отличен от единицы и аппроксимация нарушается.

Сходимость разностного решения к решению дифференциального, тем не менее, в этом случае можно доказать (с использованием разностного аналога теорем вложения Соболева). Мы не будем останавливаться на доказательстве, так как его сложность превосходит тот уровень, который был намечен для данного обзора.

2.2.4. Первое дифференциальное приближение.

Представим, что существует гладкая функция, значения которой в узлах сетки совпадают с решением разностного уравнения. Нам не важно, построена ли эта функция каким-нибудь восполнением разностной функции, важен сам факт ее существования. Предположим дополнительно, что эта функция может быть разложена в ряд Тейлора по крайней мере в некоторой части области решения разностного уравнения. Подставим эту функцию в разностное уравнение и проделаем те же выкладки, которые мы делали раньше при исследовании аппроксимации. Будем считать на первых порах, что мы сохраняем при этом все члены ряда Тейлора.

Рассмотрим к примеру явную разностную схему для линейного обыкновенного уравнения

Подставим в нее функцию U(t).

Раскладывая эту функцию в ряд, получим

Смысл такой подстановки иной, несмотря на то, что выкладки формально совпадают. Раньше мы подставляли в разностное уравнение решение дифференциального уравнения и вычисляли невязку, то есть смотрели, насколько хорошо (или плохо) решение одного уравнения удовлетворяет другому. Теперь мы подставляем функцию, которая заведомо удовлетворяет разностному уравнению, так что невязка равна нулю, но делая выкладки, мы фактически преобразуем уравнение. В подобных случаях уравнения (имеющие одинаковые решения) называются равносильными. Особенность данного преобразования в том, что одно из равносильных уравнений является разностным, а второе - дифференциальным. Тем не менее, равносильность уравнений наводит на мысль, что свойства разностного уравнения переносятся на дифференциальное, которое исследовать иногда может быть проще. Еще проще исследовать уравнение, в котором оставлены только младшие по каждой переменной члены ряда Тейлора, такое уравнение часто бывает близким по форме к известным уравнениям математической физики, так что его свойства известны. Это уравнение называют первым дифференциальным приближением разностного уравнения. Первое дифференциальное приближение уже не равносильно разностному уравнению, но, тем не менеее, обычно сохраняет некоторые важные его свойства. Для явной схемы первое дифференциальное приближение имеет вид

В некоторых случаях можно пойти на дальнейшие упрощения. Считая f = const, продифференцируем это уравнение почленно и пренебрежем третьей производной. Это даст возможность выразить вторую производную через первую, подставив ее в уравнение, получим еще одну форму первого дифференциального приближения.

Приведем для сравнения графики решения исходного уравнения V(t), первого дифференциального приближения U(t) и разностного уравнения u(k).

Рис. Точное решение, первое дифференциальное приближение, разностное.

Видно, что разностное решение действительно ближе к решению уравнения первого дифференциального приближения.

Подробно свойства дифференциального приближения, в том числе первого, а также различные формы этого приближения, изучены и обоснованы в трудах (см., например, [[29]]) и других авторов, достаточно подробный список литературы содержится там же. отмечает, кстати, что идея использования дифференциального приближения для анализа разностных схем была высказана в 50-х годах . Мне в связи с этим вспоминается, какое сильное впечатление своим умением неординарно взглянуть на вещи произвел на меня в свое время Анатолий Иванович Жуков, в частности, доклад, сделанный им на одной из межобъектовских конференций. В докладе проводился изящный анализ параметрического семейства разностных схем, включающего схемы от КРЕСТа до схемы Годунова. Он показал, что все отличие схем сводится в конечном случае к той или иной (схемной или искуственной) вязкости.

Мы не будем заниматься вопросами обоснования свойств дифференциального приближения, отсылая интересующихся к названной книге, приведем лишь несколько характерных примеров.

Рассмотрим первое дифференциальное приближение для разностного уравнения, соответствующего обыкновенному дифференциальному уравнению. Для этого приравняем нулю оператор, полученный нами раньше при исследовании аппроксимации:

Знак плюс соответствует неявной схеме, минус - явной. Что видно из этого уравнения? Во первых, видно, что его главные члены совпадают с соответствующими членами исходного дифференциального уравнения. Это означает наличие аппроксимации. Младший член пропорционален первой степени шага сетки, что определяет порядок аппроксимации. Решение уравнения легко выписывается в явном виде, мы не будем этого делать, отметим однако, что знак коффициента при производной (1 ± 0.5 cτ) для неявной схемы всегда совпадает со знаком коэффициента при функции (напомним, что коэффициент c мы считаем положительным). Отсюда следует, что решение однородного уравнения является убывающим по величине, то есть решение неоднородного уравнения устойчиво относительно возмущений. Тем самым первое дифференциальное приближение передает свойство безусловной устойчивости неявной разностной схемы. В случае явной схемы при

решение также устойчиво, при нарушении этого неравенства решение однородного уравнения становится растущим, так что малые возмущения в решении неоднородного уравнения неограниченно растут. Условие устойчивости первого дифференциального приближения таким образом совпадают с условиями устойчивости для явной разностной схемы.

Для симметричной схемы первое дифференциальное приближение имеет вид:

Раньше мы видели, что разностная схема безусловно устойчива. О первом дифференциальном приближении этого сказать нельзя. Мы столкнулись, таким образом, со случаем, когда свойство устойчивости не сохраняется при переходе к первому дифференциальному приближению.

Рассмотрим устойчивость применительно к уравнению теплопроводности. Для явной схемы первое дифференциальное приближение имеет вид

Здесь младшие члены содержат производные по разным переменным, что затрудняет анализ их влияния на свойства уравнения. Выразим поэтому вторую производную по времени через остальные члены уравнения. Для этого продифференцируем обе части уравнения по времени и подставим в смешанную производную выражение первой производной из исходного уравнения. Выкладки, конечно, должны быть очень громоздкими, но мы ограничимся только главными членами:

Подставляя полученное выражение в первое дифференциальное приближение получим новую его форму

Легко понять, что если бы мы выполнили выкладки в полном объеме, то в полученном уравнении появились бы члены второго порядка малости. В рамках первого дифференциального приближения мы ими пренебрегаем. Внешне это выглядит так, как будто мы выразили нужную нам производную непосредственно из дифференциального уравнения. На самом деле это не так, функция U здесь другая и мы совершаем эквивалентные (хотя и приближенные) преобразования над уравнениями, содержащими именно эту функцию. Заметим, что не всегда получаемые члены имеют больший порядок.

Для исследования устойчивости нам нужно исследовать поведение возмущений данного уравнения. Поскольку уравнение само является однородным, то возмущения удовлетворяют этому же уравнению. Решение можно искать в виде U = exp(λt + ikx), тогда

λ = - k2 - (1/2 τ - 1/12 h2) k4 .

Для устойчивости возмущений необходимо, чтобы λ было отрицательным при любых значениях k. Не будем останавливаться на анализе полученного выражения, так как очевидно, что результат никак не коррелирует с условиями устойчивости разностной схемы. Уравнение теплопроводности является еще одним примером, когда первое дифференциальное приближение теряет нужные нам свойства разностной схемы.

Перейдем к уравнениям газовой динамики. Рассмотрим явную схему. Для нее первое дифференциальное приближение имеет вид

,

Заменяя производные по времени, подобно тому, как мы это делали для уравнения теплопроводности, получим

Решение для возмущений будем искать в виде U = uּexp(λt + ikx), P = pּexp(λt + ikx), тогда характеристическое уравнение будет иметь вид

Корни этого уравнения имеют положительную вещественную часть

Следовательно, exp(λt) будет расти, что определяет неустойчивость уравнений первого дифференциального приближения. Явная разностная схема, как мы знаем, также неустойчива. В дифференциальном приближении для неявной схемы, как отмечалось выше, члены первого порядка входят с противоположными знаками. Это приводит к изменению знака Re(λ), что означает безусловную устойчивость уравнений дифференциального приближения.

Таким образом мы видим, что анализ устойчивости с помощью первого дифференциального приближения иногда проще, чем исследование нормы оператора шага, но не всегда дает правильные результаты. Подробнее условия применимости рассмотрены в упоминавшейся книге . Там же рассмотрены другие свойства разностных схем, которые можно исследовать с помощью дифференциальных приближений.

2.2.5. Искуственная вязкость.

На начальном этапе развития машинной математики конструирование удачной разностной схемы требовало определенной изобретательности, зачастую этому помогали соображения физического характера. Примером может служить известная схема [[30]]. Со временем пришло понимание, что основные свойства схемы - аппроксимация, зачастую устойчивость и другие определяются остаточными членами первого дифференциального приближения. Это привело к идее строить разностную схему, исходя в основном из удобства реализации, а затем добавлять в нее малые члены, улучшающие свойства схемы или компенсирующие ее недостатки. Такой подход, в частности, систематически использовался при разработке комплексов СИГМА и ЧАС. В газодинамических задачах добавочные члены часто берутся в форме вязкости, впервые такой прием был использован Дж. фон Нейманом и Р. Рихтмайером [[31]] для размазывания ударных волн, отсюда и произошел сам термин "искуственная вязкость".

Рассмотрим в качестве примера схему квазиКРЕСТ:

Характеристический определитель для этой схемы имеет вид:

Корни этого уравнения при выполнении условия устойчивости, как мы видели, являются комплексными и по модулю все равны единице. Это, с одной стороны, способствует высокой точности разностного решения. С другой стороны, это приводит к тому, что возникающие в результате погрешностей счета осцилляции не затухают со временем. График решения при этом выглядит неприглядно и, несмотря на высокую точность в среднем, погрешности в отдельных точках могут быть значительными. Введем в первое уравнение искуственную вязкость:

Здесь ν - коэффициент вязкости. Весь член, очевидно, является малым второго порядка, для того, чтобы убедиться в этом, достаточно его умножить и поделить на h2, разделенная разность при этом аппроксимирует вторую производную. С учетом добавленного члена определитель примет вид

Характеристическое уравнение становится следующим:

λ2 -2{1 - 2[(τ /hνh2]sin2 (ω /2 )} λ + 1- 4 ν τ h2sin2 (ω /2 ) = 0 .

Изменение во втором коэффициенте несколько ослабляет условие устойчивости, мы не будем на этом останавливаться. Свободный член теперь (при очевидных ограничениях на ν τ h2) меньше единицы, так что таковыми будут и модули корней. Иными словами, возмущения будут затухать. При этом скорость затухания зависит от ω, она больше для высокочастотных осцилляций и мала для гладких возмущений. За счет подбора величины коэффициента искуственной вязкости можно добиться требуемого качества решения, сделав график достаточно плавным и не слишком сильно исказив гладкую часть решения, определяемую низкими гармониками. В этом отношении искуственная вязкость удобнее схемной вязкости, так как содержит в явном виде параметр, позволяющий регулировать ее величину.

Дополнительная степень свободы появляется, если вязкость считать независимо от разностной схемы на добавочном шаге. Например, применительно к газодинамике, можно ввести пересчет

После этого значение us используется вместо un в обычной схеме, например, квазиКРЕСТ. Такую процедуру можно трактовать, как сглаживание сеточной функции. Технически сглаживание обычно применить легче, чем искусственную вязкость, особенно, если его нужно вводить в готовую программу. Именно так оно было применено, по-видимому, впервые в комплексе СИГМА, когда мы с совместно с и решали одну сложную задачу. При этом после ряда экспериментов мы остановились на использовании четвертых разностей (вместо вторых, типичных для физической вязкости). Четвертые разности сильнее дискриминируют высокочастотные осцилляции и, тем самым слабее изменяют основное решение. Позднее сглаживание было более подробно исследовано .

Отметим, что сглаживание, вообще говоря, дает несколько иной эффект, чем прямое введение вязкости в разностную схему. Это связано с тем, что фактически мы применяем здесь схему расщепления. Подробнее на схемах расщепления мы остановимся позже, здесь только заметим, что из-за малости поправок их влияние сопоставимо.

2.2.6. Коэффициент потерь.

Аппроксимация и, при соответствующих условиях, сходимость разностной схемы означают, что выбирая достаточно мелкий шаг сетки мы можем обеспечить сколь угодно высокую близость разностного решения к решению дифференциального уравнения. Означает ли это точность решения физической задачи? Не всегда. Во-первых, как мы видели, в самих дифференциальных уравнениях используются подчас довольно грубые модели, такие как диффузионное приближение переноса излучения, модели детонации, упругопластичности. Уравнения состояния задаются некоторыми аппроксимационными выражениями, также имеющими ограниченную точность. Константы уравнений состояния и других моделей измеряются, как правило, в эксперименте и их точность тоже невелика. Все эти погрешности остаются и их нельзя уменьшить, изменяя шаг сетки.

Кроме того, при измельчении шага сетки растет число операций, что приводит к накоплению погрешностей округления, в результате, начиная с некоторой величины шага суммарная погрешность начинает расти. Этот факт хорошо известен для алгебраических задач. Применительно к системам разностных уравнений ситуация усугубляется ухудшением обусловленности при увеличении порядка системы. Существуют и другие причины. Уравнения идеальной (невязкой) жидкости идеализируют реальную физическую задачу. На конечной сетке разностное уравнение, имеющее определенную аппроксимационную вязкость, может оказаться ближе к физической задаче, чем к идеализированной. С уменьшеним размера сетки в разностном решении могут появиться такие эффекты, как турбулентность, перемешивание и другие, которые не закладывались в рассматриваемую модель.

Так или иначе, расчет всегда проводится на конечной сетке и желательно оценить погрешность именно для данного шага. В явном виде для большинства схем это невозможно, но существуют некоторые подходы, позволяющие сравнивать схемы между собой и выбирать разумный шаг сетки. Один из таких подходов, связанный с понятием коэффициента потерь [[32]], мы вкратце рассмотрим.

Будем исходить из того, что при решении конкретной задачи мы ориентируемся на некоторое значение масштаба разрешения L и допустимую погрешность ε. О погрешности следует поговорить особо. Обычно применяются, как известно, абсолютная или относительная погрешности. При оценке точности решения физической задачи естественно ориентироватся на абсолютную погрешность, однородную во всей рассматриваемой области. Однако, для того, чтобы охватить задачи, отличаюшиеся масштабами рассматриваемых величин, погрешность следует соотносить с характерным или максимальным значением величины в задаче. Таким образом, мы будем оценивать качество решения с помощью относительной погрешности, но отнесенной не к значению рассчитываемой величины в каждой точке, а к максимальному значению в задаче.

Рассмотрим в качестве примера обыкновенное уравнение

Допустим, что нас интересуют значения функции на сетке с шагом, равным масштабу разрешения L. Функцию f будем для простоты считать постоянной (или почти постоянной на расстояниях масштаба разрешения). Тогда точное (или почти точное) решение на одном интервале может быть представлено следующим соотношением

Возьмем неявную разностную схему

Естественным и экономичным было бы решение этого уравнения на сетке с шагом L. Однако, такой подход не обеспечивает, вообще говоря, нужной точности. Возьмем поэтому более мелкий шаг h = L/p, где p будем считать целым и назовем его коэффициентом потерь. Он показывает, во сколько раз должно возрасти число операций для достижения нужной точности. Коэффициент потерь зависит, таким образом, от масштаба разрешения и требуемой точности. Для рассматриваемой схемы решение может быть представлено в следующем виде

Для обеспечения заданной точности (в предположении, что u и f неотрицательны) достаточно, чтобы

Мы оцениваем погрешность, приобретаемую на одном этапе расчета и не рассматриваем суммирования погрешностей. Это связано с тем, что мы не стремимся к получению строгих (сильно завышенных) оценок, а хотим выделить определяющую причину потери точности. Это позволит нам сравнивать отдельные схемы между собой.

Левая часть неравенства стремится к нулю при c → 0 и c → ∞. Наибольшее значение погрешности достигается при c ≈ 1/L, при этом коэффициент потерь должен удовлетворять условию |1/e - 1/(1+1/p)p|< ε. График левой части, как функции от p имеет следующий вид (синие столбцы)

Погрешность, равная 1% достигается при p = 18.

Для симметричной схемы погрешность оценивается выражением

График этой погрешности при c = 1/L изображен коричневыми столбцами. Видно, что при этом значении c погрешность значительно меньше, точность порядка 1% достигается уже при p = 2. Однако, для симметричной схемы максимум погрешности достигается не в этой точке, а при c → ∞. Более полное сравнение неявной и симметричной схем дает следующая таблица

Коэффициент потерь
	Eps\ cL	0.01	0.1	1	10	100	104	106	108
Неявн. сх.	4%	1	1	4	2	1	1	1	1
1%	1	1	18	3	1	1	1	1
0.1%	1	50	180	8	2	1	1	1
Симм. сх.	4%	1	1	1	3	9	90	900	9000
1%	1	1	2	4	11	107	1070	10700
0.1%	1	1	6	9	13	130	1300	13000

Видно, что у симметричной схемы коэффициент потерь неограниченно растет с ростом произведения cL. Это ограничивает круг задач, для которых можно рекомендовать симметричную схему. Но зато для задач, у которых cL ≤ 1, симметричная схема значительно экономичней, особенно, если требуется высокая точность. Это объясняется более высоким порядком аппроксимации этой схемы, что соответствует более быстрому стремлению погрешности к нулю при уменьшении L.