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

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

Здесь k - некоторый коэффициент. В книге Зельдовича-Райзера [[13]] показано, что это предположение оправдано для случая, когда пробег излучения мал по сравнению с размерами области, в которой ведется расчет. При этом

С учетом указанных упрощений получаются уравнения спектральной неравновесной диффузии:

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

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

В таком виде уравнение теплопроводности было использовано во всех основных программах, созданных в 50-60х годах и применявшихся для расчета современных для того времени изделий. Ситуация, как видим, опять напоминает анекдот о пьяном матросе. Однако, в данном случае проблема гораздо серьезнее. При счете детонации приближенная модель макрокинетики при надлежащем подборе параметров достаточно точно описывает штатные режимы. Уточненная двухстадийная модель позволяет рассчитывать процессы возбуждения и затухания детонации. В целом это обеспечивает все основные потребности разработчиков типовых изделий. Иначе обстоит дело с теплопроводностью. Распространение тепла в узком зазоре, ширина которого сравнима с пробегом излучения, является одной из наиболее важных задач и именно она плохо описывается диффузионным приближением. Очень серьезное внимание на это обращал , в свое время он пытался использовать в качестве некоторой поправки так называемый "геометрический пробег". Известно также, что диффузионное приближение плохо передает форму фронта тепловой волны и поток энергии на границе с вакуумом. Усреднение пробегов по частоте заметно влияет на динамику процессов. Заметим, что все эти недостатки типичны также для уравнения переноса нейтронов. Но для нейтронов отказ от диффузионного приближения произошел достаточно быстро, уже в конце 50х годов существовал целый ряд различных методик счета многогруппового (по спектру) уравнения переноса. Для переноса излучения процесс внедрения заметно затянулся, так как задача оказалась более сложной.

Формально диффузионные уравнения для нейтронов и для фотонов выглядят одинаково и относятся к так называемым параболическим уравнениям. Это означает, в частности, что у них нет характеристик и скорость сигнала является бесконечной. Я помню, как при анализе одного из расчетов объяснял мне это таким образом: "Говорят, что если бросить щепотку соли в Одессе, то в Рио-де-Жанейро вода в тот же момент станет чуточку солонее". Тем не менее, есть внутренняя динамика процессов, связанная с величиной коэффициентов уравнений и приводящая к существенному отличию в расчетах. Забегая вперед, скажу, что разностные схемы, используемые для счета диффузии тепла, предъявляют существенно более высокие требования к памяти вычислительных машин и, отчасти, к их быстродействию. Аналогичные различия сохраняются и, более того, многократно возрастают при переходе к уравнению переноса. Это собственно и определило невозможность расчета уравнений переноса излучения на машинах М-20 и БЭСМ-6 с их памятью в несколько десятков килобайт.

Тем не менее, некоторые шаги в этом направлении делались. В частности, мы с в начале 70х годов рассмотрели интегральный метод, обобщающий метод и , реализованный в И-0. Метод оказался достаточно сложным и программно был реализован тогда только для плоского одномерного случая, результаты были изложены в одном из отчетов ВНИИЭФ и в сокращенном виде в моей книге [[14]]. Некоорых результатов удалось добиться также в методе расчета спектральной неравновесной диффузии, результаты того периода также изложены в одном из отчетов и, также сокращенно, в книге [[15]]. В дальнейшем, насколько мне известно, со своими сотрудниками успешно развивает работы по спектральной неравновесной диффузии.

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

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

" Сидел, думал ".

Через некоторое время матрос попросился обратно на корабль. Возможно, что этому матросу, осознанию им своего места, мы обязаны тем, что имеем сегодня труды , , … Мне, к сожалению, с "матросом" повезло меньше и занятия переносом излучения пришлось оставить.

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

Решение будем искать в виде волны, двигающейся с постоянной скоростью V, для чего введем автомодельную переменную x=x-Vt. Тогда уравнение примет вид

Первый интеграл находится легко, поскольку нас интересует пример, возьмем для простоты константу интегрирования равную нулю:

Получили уравнение с разделяющимися переменными, при его интегрировании снова возьмем константу равной нулю. Кроме того, выразим T через x и подставим x=x-Vt. Полученное решение имеет вид:

Решение определено в области x£Vt, в точке x=Vt температура обращается в ноль, эту точку будем считать фронтом волны, идущей по нулевому фону с конечной скоростью. Если задача рассматривается в ограниченной области, то на границе левее фронта нужно задать температуру, растущую по определенному закону. При другом граничном условии скорость волны будет конечной, но не постоянной. При k>1, что реально соответствует физике процесса, производная температуры обращается в точке фронта в бесконечность. Это вызывает определенные проблемы при построении разностных схем, однако, такая особенность связана не со свойствами тепловой волны, а с особенностями диффузионного уравнения. При счете уравнения переноса на фронте волны наблюдается острый носик.

При счете теплопроводности совместно с газодинамикой возникают две скорости: скорость тепловой волны V и скорость звука C. В зависимости от их соотношения наблюдаются два различных типа течения. При V>C тепловая волна распространяется по неподвижному холодному веществу, движение газа, обычно достаточно плавное, происходит за фронтом волны. При V<C в районе фронта возникает ударная волна, которая уходит вперед, тепловая волна в этом случае распространяется по сжатому и движущемуся веществу. Подобный эффект используется, например, при сжатии лазерных мишеней.

Глава 2. Разностные схемы.

В этой главе мы рассмотрим некоторые вопросы численной реализации моделей механики сплошной среды. Основные понятия теории разностных схем в основном излагаются на достаточно простом "разговорном" уровне. Для серьезного изучения можно рекомендовать книги и [[16]], Р. Рихтмайера и К. Мортона [[17]], [[18]], [[19]] или несколько более легкие и [[20]] и [[21]].

2.1. Дискретизация.

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

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

Оба подхода начали применяться в самом начале становления машинного счета. В задачах обтекания, где рассчитывалось движение однородного газа, естественно было применить эйлеровы сетки. В задачах ядерной физики и физики взрыва, напротив, на первых порах использовались лагранжевы сетки. В таких задачах, как правило, имеется много контактных границ, которые можно выделить линиями лагранжевой сетки и в дальнейшем это выделение сохранять в процессе счета, что существенно влияет на точность счета всей задачи. Кроме того, лагранжеву сетку легко строить неравномерной, выбирая более мелкие ячейки в тонких слоях, а также используя в двумерных задачах сильно вытянутые ячейки там, где это позволяет вид геометрии и характер предполагаемого течения. В целом это позволяет экономить число ячеек, что определяет экономию машинной памяти и времени счета. На первых машинах эти соображения играли определяющую роль. Особенно серьезным на наших машинах было ограничение памяти (на СТРЕЛЕ и М-20, соответственно, 10 и 20 Кбайт). В результате первые двумерные газодинамические программы , , использовали именно лагранжевы сетки. Опыт работы с ними показал и недостатки лагранжевых сеток - в сдвиговых течениях сетки "портились" и счет становился невозможным. В одномерных задачах этот эффект отсутствует, в силу чего лагранжевы сетки в одномерных комплексах сохранили свое значение в полной мере. В двумерных комплексах для борьбы с искажением сеток применялись разные подходы - переинтерполяция, использование подвижных сеток, нерегулярных сеток, метод частиц. Об этом мы будем говорить позже.

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

2.1.1. Сеточные функции.

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

Рассмотрим вначале вопрос о форме представления, для простоты изложения будем считать, что мы имеем дело с функцией одного аргумента U(x). Покроем область изменения аргумента некоторой сеткой с узлами x0,…, xi,…, xm. В узлах зададим некоторые численные значения u0,…, ui,…, um. Совокупность этих значений называется сеточной функцией. Часто значения сеточной функции удобнее определять не в узлах сетки, а в центрах ячеек. Обычно их обозначают в этом случае с использованием дробного индекса: u1/2,…, ui+1/2,…, um-1/2. Заметим также, что сеточная функция может определяться в каждом узле или ячейке несколькими параметрами, так поступают, например, при использовании сплайнов, в методе конечных элементов. Мы на этих случаях подробно останавливаться не будем. Для краткости сеточную функцию будем обозначать{u}.

2.1.2. Аппроксимация функций.

Займемся вопросом о точности представления непрерывной функции U(x) сеточной функцией {u}. Сложность их сравнения заключается, прежде всего, в том, что, говоря языком функционального анализа, эти функции принадлежат разным функциональным пространствам. Существует два способа приведения их к общему виду.

Первый заключается в том, что сеточная функция дополняется до некоторой функции непрерывного аргумента u(x). Например, строится ломанная линия, соединяющая значения сеточной функции в узлах отрезками прямых. Это соответствует линейной интерполяции функции u(x) между узлами. Для случая, когда значения сеточной функции определены в центрах ячеек, зачастую используется еще более простая форма дополнения ее кусочно постоянной функцией. Иногда, впрочем, используются и более сложные восполнения. Так или иначе, после того, как функция u(x) определена, можно сравнивать ее функцией U(x), используя те или иные функциональные нормы. Чаще всего используются норма непрерывных функций, коротко - непрерывная норма, норма C, а также среднеквадратическая или, как ее еще называют, энергетическая норма, норма L2:

Последней нормой мы будем пользоваться чаще и для краткости индекс L в дальнейшем будем опускать.

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

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

Другой способ сравнения рассматриваемых функций использует преобразование функции U(x) в сеточную функцию, или, как говорят, проектирование ее на сетку, после чего можно использовать нормы сеточного пространства:

Реально функция U(x) обычно бывает неизвестна, поэтому подобные оценки как в функциональных, так и в сеточных нормах, чаще используются в теоретических рассуждениях, чем на практике.

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

» 0.5 |U(x)'| h.

Соответсвенно этому погрешность на равномерной сетке можно оценить величиной

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

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

2.1.3. Аппроксимация уравнений.

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

Линейного обыкновенного уравнения

Уравнения теплопроводности

Уравнений газодинамики

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

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

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

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

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

Явная схема

Неявная схема

Симметричная схема

Явная схема и здесь позволяет шагами по времени легко вычислить все требуемые значения сеточной функции, если заданы начальные (при n=0) и граничные (при k=0 и k=m) значения этой функции. Неявная и симметричная схемы требуют на каждом шаге по времени решения системы линейных уравнений.

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

2.2. Свойства разностных схем.

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

2.2.1. Аппроксимация.

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

Рассмотрим вначале основные понятия на примере обычного уравнения f(x) = 0. Предположим, что мы нашли приближенное значение xp корня x* этого уравнения. Нас интересует погрешность найденного значения Δx = xp - x*. Чтобы найти погрешность, нужно знать точное значение, а оно, естественно, неизвестно. Чтобы оценить качество найденного приближения, можно подставить его в уравнение и посмотреть, удовлетворяет ли оно уравнению. Как правило, оно удовлетворяет уравнению неточно, остается некоторый остаток R, который называют невязкой уравнения:

f(xp) = R.

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

f(xp) = f(x*) + f '(x*) (xp - x*) + …≈ f '(xp) Δx.

получим

Δx ≈ R / f '(xp).

Если производную можно вычислить, то мы получаем таким образом оценку погрешности. В противном случае приходиться ориентироваться непосредственно на величину невязки. Смысл такого подхода заключается в том, что мы трактуем найденное приближенное решение как точное решение возмущенного уравнения f(x) = R. Зачастую это даже удобнее, так как невязку можно сопоставить с другими погрешностями в уравнении, например, с первоначальными погрешностями входящих в него коэффициентов. Если величина R сравнима с этими погрешностями, то точность решения следует считать удовлетворительной. Я намеренно здесь не уточняю ряда понятий, так как хочу пояснить только смысл невязки. Отмечу, что очень продуктивно такой подход использован в работах Дж. Х. Уилкинсона (см., например, [[22]]).

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

2.2.1.1. Обыкновенное дифференциальное уравнение.

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

Будем считать, что U(t) - решение дифференциального уравнения, причем

U = U(tn), U(t + t) » U + U' t + 0.5 U''t2 +…

Подстановка в уравнение дает

где R - искомая невязка. После очевидных преобразований получим

Поскольку U - решение дифференциального уравнения, то первые три члена в левой части в совокупности обращаются в ноль, так что

Невязка пропорциональна шагу сетки t и имеет первый порядок малости относительно этого шага. Соответственно, порядок аппроксимации разностного уравнения считается равным единице. Членами более высокого порядка малости мы пренебрегли. Невязка, напомним, характеризует отличие разностного уравнения от исходного дифференциального.

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

После преобразований получим:

Членом второго порядка малости мы пренебрегли. Внешне невязка имеет другой вид, хотя порядок малости не изменился. Однако, если считать функцию U(t) достаточно гладкой, то дифференцируя почленно обе части исходного уравнения, получим

откуда следует, что невязки в сущности одинаковы. Выкладки, однако, получаются проще в первом случае за счет более удачного выбора опорной точки.

Для неявной схемы выкладки аналогичны, а уравнение с невязкой имеет вид:

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

Знак плюс здесь соответствует неявной схеме, минус - явной.

Для симметричной схемы разложение удобно делать в окрестности точки tn+ 0.5t, при этом нужно удерживать больше членов. Результат представляется в виде:

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

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

или ct < 2.4. На более грубой сетке невязка схем первого порядка меньше. Таким образом, более высокий порядок аппроксимации свидетельствует не столько о величине ошибки, сколько о скорости ее уменьшения при измельчении сетки. На конечной сетке соотношение погрешностей может быть другим и это важно понимать.

2.2.1.2. Уравнение теплопроводности.

В уравнении теплопроводности искомая функция зависит от двух независимых переменных x и t, что необходимо учитывать при разложении ее в ряд в окрестности выбранной опорной точки. Для явной схемы

опорной удобно взять точку U =, тогда

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

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

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

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

2.2.1.3. Уравнения газодинамики.

Для уравнений газодинамики (мы рассматриваем упрощенную форму этих уравнений)

явную разностную схему естественно взять в виде:

Исследование аппроксимации, выполненное подобно тому, как это делалось для уравнения теплопроводности, приводит к следующим выражениям:

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

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

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