(5.1)
Коррекция алгоритма состоит в следующем: пусть
(5.2)
, где .Тогда треугольник с номером k исключается из расчета.
Фактически это означает, что в дискретном варианте формул для (1.1) и (1.8) исключается некоторая (незначительная) часть области интегрирования, отвечающая таким выродившимся треугольникам. (Заметим, что она покрывается другой парой треугольников). Отметим также, что предлагаемые различными авторами «безавостные» модификации (например, замена обращающегося в нуль якобиана в знаменателе функционала на нечто иное), по-видимому, вполне
работоспособные при расчете стационарных сеток, приводят к неконтролируемым значениям скоростей узлов сетки при расчетах нестационарных сеток.
Описанный выше прием (5.2) исчерпывает меры, связанные с вырождением четырехугольных ячеек в треугольники.
Значительно сложнее обстоит дело с учетом возможного появления невыпуклых ячеек. Их появление обнаруживается в случае, если окажется, что для некоторого («плохого») треугольника
(5.3)
.Гипотетически возможны следующие альтернативные варианты разрешения таких ситуаций.
а) Не обращать внимание (ничего не менять в расчетных формулах). Например, в случае, изображенном на рис.4г, если единственная «плохая» (невыпуклая) ячейка является «угловой» для расчетной области, соответствующий «плохой» треугольник получит номер k=3 и не будет участвовать в расчетных формулах, как уже отмечалось в § 3.
Заметим, что придется позаботиться о том, чтобы программа контроля сетки на выпуклость пропустила «провинившуюся» ячейку, не провоцируя появление новых невыпуклых ячеек.
б) «Плохой» треугольник исключается, т.е. не участвует в описанных расчетных формулах. Такой вариант в принципе допустим. Например, для расчета гармонических сеток в работе [13] рассматривался алгоритм, при котором разрезание невыпуклых ячеек осуществлялось только с помощью одной («хорошей») диагонали, и даже для выпуклых ячеек выбиралась одна диагональ – та, при которой реализуется меньшее значение функционала. Работоспособность такого алгоритма подтверждалась практическими расчетами.
в) В формуле (1.2) заменяем величину g0 на
. Очевидно, что при этом будет сохранена положительная определенность всех матриц Gk. Для функционала с якобианом в соответствии с упоминавшейся теоремой из работы [7], существование невырожденной сетки для такого набора матриц Gk гарантировано. Проблема состоит в том, как одолеть вариационный барьер, за которым оказался соответствующий из узлов сетки, и какая скорость его движения будет при этом выработана.г) Принимаются меры по исключению возможности появления невыпуклой ячейки сетки. В качестве такой меры предлагается, например, использовать следующий прием. В случае (5.3) узел с номером k проектируется на отрезок, соединяющий узлы с номерами
(k-1) и (k+1), и точка с полученными координатами
заносится вместо . Невыпуклая ячейка превращается при этом в треугольник. Для полноты изложения приведем соответствующие расчетные формулы. Полагаем, .
Согласно уравнению прямой, соединяющей узлы
и :(5.4)
, .Согласно уравнению перпендикуляра через точки
и :.
Отсюда получаем формулу для
:(5.5)
.После этого вычисляем
по формулам (5.4).Конечно, следует иметь в виду, что такое волевое превращение невыпуклой ячейки в треугольник может привести к бесконтрольному увеличению скорости соответствующего узла сетки.
Тем не менее, исходя из практического опыта, заметим, что допущение в расчете сеток невыпуклых ячеек может приводить к весьма неблагоприятным последствиям, поскольку они могут приобретать очень «плохую» форму. При этом вряд ли можно ожидать разумной аппроксимации даже первых производных, необходимых для последующего решения основной задачи.
Поэтому, исходя из обсуждения возможных альтернативных вариантов действий в случае появления невыпуклых ячеек, по-видимому, целесообразно остановиться на последнем. А именно, необходимо уже в исходных данных избавиться от невыпуклых ячеек (например, описанным приемом) и в дальнейшем не допускать их появления. Соответствующая технология, опирающаяся на подбор в итерационном процессе значений параметров
, была описана в § 4.Отметим также, что рекомендуемое решение исходит из интересов расчета газодинамических задач на регулярных четырехугольных сетках. Оно может быть другим, например, в случае, если основная задача решается на треугольных сетках или в ситуациях б) и в), когда присутствие невыпуклых четырехугольных ячеек вполне допускается.
§ 6. Проблема контроля скоростей узлов сетки.
В работе [14, стр.17] С.К.Годунов написал: «… расчет выделенных ударных волн, движение которых вызывает перемещение точек двумерной или трехмерной сетки, также не изучен до сих пор ни теоретически, ни экспериментально с исчерпывающей подробностью. Заведомо не исследована обратная связь такого перемещения точек на движение волны, нет ясности в том, какие эффекты эта обратная связь вызывает».
В подтверждение этих слов стоит сказать, что накопленный отрицательный опыт связан в первую очередь с недопустимыми скоростями узлов сетки, появляющимися в расчетах из-за несовершенства алгоритмов конструирования сеток. Безусловно, скорость движения сетки весьма существенно влияет на расчет основной решаемой задачи. Это можно видеть уже на примере элементарной задачи о распаде разрыва, являющейся идейным стержнем методики расчета газодинамических задач, как описано в монографии [2].
В случае подвижных границ алгоритм построения сеток должен непрерывно обеспечивать адекватное изменение положения узлов сетки внутри области. Казалось бы естественным требование, чтобы скорости «внутренних» узлов сетки были заключены в пределах, в которых изменяются скорости «граничных» узлов. (Это аналогично известному «принципу максимума» для решений эллиптических уравнений). Однако, несмотря на кажущуюся его естественность, легко могут быть придуманы примеры, когда такое требование не выполнено, например, при применении для расчета сеток интерполяционных формул.
Сейчас уместно более подробно осветить этот вопрос, как было обещано выше в § 2. В работе [15] представлены результаты сравнения нескольких вариантов построения двумерных разностных сеток посредством интерполяционных формул.
В качестве одного из «оправдавших ожидания» алгоритмов можно рекомендовать использование так называемой трансфинитной интерполяции. В простейшем двумерном варианте интерполяционная формула для функции
на единичном квадрате , , по ее значениям на контуре квадрата, имеет вид:(6.1)
Здесь
- монотонно возрастающие функции своих аргументов на отрезке [0,1], причем, .
Легко проверить, что для произвольных функций
, удовлетворяющих этим условиям, будем иметь равенство интерполянты заданным граничным значениям.Рассматривая последовательности координат узлов сетки, заданных на границах области, как значения в точках
, некоторых непрерывных функций, и применяя формулу (6.1) независимо для каждой из двух пространственных координат х и у, получаем интерполяционные формулы для расчета координат внутренних узлов сетки.