Обратный характеристический метод на нерегулярных разностных сетках

При численном интегрировании многомерной системы уравнений гиперболического типа (1.2.2) в сложных областях (например, многосвязных и т. д.) далеко не всегда удается преобразованием независимых переменных (1.2.7) Ьривести область интегрирования к простому и удобному для введения регулярной разностной сетки виду (типа (1.2.6) и т. п.). Кроме того, для ряда подвижных систем координат (1.2.10), (1.2.12) в процессе счета возникают большие деформации разностной сетки, а часто происходит и полная потеря ее регулярности (нулевые объемы, ''выворачивание" лагранжевых сеток в задачах механики сплошных сред и т. д.). Поэтому представляет несомненный интерес построение численных методов для случая, когда на слое начальных данных t = tⁿ и рассчитываемом слое t = = /"⁺¹ = tⁿ + т расположение сеточных узлов может быть произвольным (или управляемым с помощью независимых от способа аппроксимации решаемой системы уравнений алгоритмов, и в частности связанных с поведением искомого решения). Примерами подобных численных методов являются метод свободных точек В. Ф. Дьяченко [69] и др.

Ниже рассматривается один из таких подходов, основанный на аппроксимации условий совместности вдоль линий пересечения координатных и характеристических поверхностей (разд. 3 данной главы) и использовании линейной интерполяции значений искомой вектор-функции и в точках пересечения этих линий с плоскостью t = tⁿ по ее значениям в трех ближайших к каждой из этих точек пересечения сеточных узлах, составляющих охватывающий каждую точку пересечения треугольник. Для многомерного уравнения переноса можно показать (см. гл. IV), что такая разностная схема является схемой с положительной, по Фридрихсу, аппроксимацией (монотонной или мажорантной — по другой терминологии), а из монотонных схем (как известно, все они имеют первый порядок аппроксимации [57, 60]) эта схема является наиболее точной. Соответствующие доказательства легко обобщаются также на случай системы (1.2.2), если матрицы Aj попарно коммутируют (AfAj = AjA_t). Как известно, в этом случае система (1.2.2) распадается на / не связанных друг с другом многомерных уравнений переноса.

1. Ограничиваясь в (1.2.2) для определенности случаем двух пространственных координат *2 (/=2),/= 0 и используя построения, аналогичные приведенным в разд. 3, запишем для произвольной внутренней точки Н на слое t = fⁿ⁺¹ с координатами Хщ, х_2И следующую схему расщепления, несколько отличающуюся от схемы первого порядка точности (3.4) :

В (1), как и в разд. 3, u", i = 1,…, / и S**^l₉ i = /+ 1,…, 27- значения вектор-функции и в точках пересечения i-й ''характеристики" с плоскостью t = const, й⁺¹ $ ^{s u}* f, i = !>••*" 2/) — значения параметров в точке плоскости t = const с координатами х_х я, х₂н (^в которую проектируется точка И). Индексы 1 и 2 в соотношениях (1) для симметризации схемы расщепления целесообразно циклически менять местами [9]. Координаты точек пересечения i-й ''характеристики" с плоскостью t = const (рис. 3.8) имеют вид.

Пусть в фиксированной в пространстве системе координат Ху, х₂ (для упрощения выкладок декартовой) на слое t = tⁿ заданы координаты х{, некоторой совокупности определенным образом занумерованных (J = 1, 2,…) узловых точек, а также некоторая дополнительная информация об этих точках (характеризующая, например, их соседство, принадлежность к границе области интегрирования и т. д.). Выберем одну из точек пересечения ''характеристик" с плоскостью t = tⁿ (2) и три узловых точки, таких, что треугольник BiB₂B₂ (рис. 3.8) с вершинами в этих узлах (обозначим их координаты через = I х{, х{ ] ,/ = 1,2,3) содержит внутри себя выбранную точку из (2) с координатами г, — = •.

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

Выбор нужных треугольников для'''больших" временных шагов т является наиболее сложной и громоздкой частью данного и других подобных алгоритмов, так как в общем случае связан с перебором большого числа возможных вариантов. Существенное уменьшение объема вычислений может быть достигнуто, когда т сопоставимо с шагом интегрирования, выбираемым из условий, аналогичных условию Куранта—фридрихса—Леви.

Рис. 3.8.

Линейная комбинация г = { x_tiх₂ — Ti*^{1 +} Чг*^{2 +} 7з^{г3 П}Р^И выполнении условий.

описывает все точки, лежащие внутри и на границе треугольника BiB₂B$. Тогда для проверки условия принадлежности точки г* = уцт¹ + y_2i *^{2 + +} 7з/^г3 треугольнику BiB₂B₃ необходимо решить линейную относительно ЧцуЧцу 7з/ систему.

где и¹, и², и³ — значения вектор-функции и в узловых точках г¹, г², г³.

2. Соотношения (1) — (6) являются расчетными формулами для случая, когда точка Н лежит внутри области интегрирования и ни одна из точек (2) на слое t — tⁿ не выходит за пределы области интегрирования. Если же в какой-либо из рассчитанных точек Я реализуется ситуация, изображенная на рис. 3.9, при которой одна или несколько ''характеристик'* пересекаются с границей области интегрирования В_{В₂В^ в некоторой точке В (ас плоскостью t — tⁿ в точке В, лежащей за пределами области интегрирования), то для определения параметров в такой точке Н необходимо привлечь одну (как показано на рис. 3.9, точка В₃) или две узловых точки на слое t = t^n+l = tⁿ + т, принадлежащих границе области интегрирования. В остальном расчет такой внутренней точки Н не отличается от рассмотренного в п. 1. В связи с подобной возможностью расчет нового временного слоя должен начинаться с расчета всех граничных узловых точек, который будет описан в п. 3.

Для ситуации, аналогичной изображенной на рис. 3.9, предполагается, что три нужные (удовлетворяющие условиям (3), (4)) узловые точки / = 1,2 (лежащие на слое t = tⁿ) и / = 3 (на слое t = f^n+l) найдены, для определения уи, у₂/" 7з/ ^в интерполяционной формуле (6) имеем вместо (5) линейную систему.

где *1/, х₂1 по-прежнему определяются из соотношений (2). Аналогичным образом может быть рассмотрен случай, когда для интерполяции привлекаются две узловые точки на слое t = tⁿ*¹, принадлежащие границе области интегрирования.

3. Для расчета точки Я, принадлежащей границе области интегрирования (рис. 3.10), сделаем локальный переход от переменных t, x_it х₂ к новым независимым переменным Г, ?*,, например составляющим систему координат с базисом.

где точка В₂ (J = 2) с координатами | ^, х, х I является ближайшей к Я точкой на слое t = Г^л, принадлежащей границе области интегрирования, а точки В1 U = 1) с координатами ⁿ, Xi, x₂ J и В_г (J = 2) с координатами {/", *|, х₂ I — соседние с точкой В₂ и также расположены на границе. Тогда.

Рис. 3.10.

и в переменных 7, ?*, |² исходная система (1.2.2) переходит в (1.2.8), в которой ?{,, 7, определяются дифференцированием соотношений.

(7). Если, как это изображено на рис. 3.10, первые р собственных значений матрицы А. 1 неположительны, то расчетными формулами будет разностная аппроксимация условий совместимости (1.5.2), например,.

и I — р граничный условий (1.5.4). Интерполяционными формулами для определения в (8) u?, и/ = 1,…, р, / + 1,…, 2/ будут по-прежнему соотношения (5), (6) с учетом того, что теперь х_и = Ху (? *,? ²), х₂/ ⁼ = х₂ ({/, 5/) в соответствии с (7), а.

(т.е. по направлению ?² используется одномерная интерполяция, рис. 3.10). Здесь не рассматривается алгоритм получения координат узловых точек Xih, х₂н ^на слое t = Г^л+1 через их значения (х{₉ х4, / = 1,2,…) на слое t = tⁿ (и через значения v! для случая, когда положение границы области интегрирования зависит от искомого решения). Как уже отмечалось, это является независимой задачей и решаться она может по-разному для каждого конкретного случая.

Описанная выше явная разностная схема устойчива и монотонна при любых значениях шага интегрирования т (естественно, при А^А₂ — А₂А иау = const, к — 1,2).