Исследование погрешности вихревых методов путем численного эксперимента Текст научной статьи по специальности «Математика»

УДК 532.5

Ю. В. БРАЗАЛУК, Д. В. ЕВДОКИМОВ, Н. В. ПОЛЯКОВ

Днепропетровский национальный университет имени Олеся Гончара

ИССЛЕДОВАНИЕ ПОГРЕШНОСТИ ВИХРЕВЫХ МЕТОДОВ ПУТЕМ ЧИСЛЕННОГО ЭКСПЕРИМЕНТА

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

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

Iu. V. BRAZALUK, D. V. YEVDOKYMOV, N. V. POLIAKOV

Dnepropetrovsk National University named after Oles Honchar

INVESTIGATION OF VORTICAL METHOD ERRORS BY NUMERICAL EXPERIMENTS

Annotation

Vortical methods of computational potential theory attract attention of investigators in computational fluid dynamics last years due to their simplicity and effectiveness. However the traditional vortical methods cannot take into account viscosity of the fluid, what strongly restrict their application field. There were a lot of attempts to incorporate a viscous dissipation into methods of discrete vortices, but they did not provide any sufficient progress. A new approach to modeling of fluid viscosity is proposed in the present paper. It is based on comparison of direct and inverse motion of discrete vortical structure. The obtained difference is considered as a result of diffusive velocity action. A relation of calculated and diffusive shifts of discrete vortices gives relation of so-called computational viscosity and physical viscosity of fluid. Then the proposed approach gives an opportunity to model a physical fluid viscosity by computational viscosity. The proposed approach can be used for improvement of algorithms of vortical methods.

Key words: vortical methods, computational viscosity, physical viscosity, discrete vortex method, direct and inverse motions of vortices, evolution of vortical structures

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

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

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

Постановка проблемы. При теоретическом анализе погрешности конечноразностных алгоритмов решения эволюционных краевых задач [1, 2] были обнаружены две основных формы ошибки расчёта - схемная диффузия и схемная дисперсия. Не останавливаясь подробно на свойствах и природе указанных вычислительных эффектов, отметим только, что в большинстве случаев схемная диффузия преобладает. По аналогии для задач вычислительной гидромеханики вводится понятие схемной вязкости. К сожалению, для вихревых методов - метода дискретных вихрей [3-5] и комбинированного метода граничных элементов и дискретных вихрей [6, 7] аналогичные понятия ввести не удалось, хотя такие попытки предпринимались неоднократно. Причиной указанных затруднений является специфическая форма дискретновихревой аппроксимации: схемная вязкость в конечноразностных алгоритмах вводится в результате анализа локальных аппроксимаций, но дискретные вихри вообще не предполагают локальной аппроксимации. То есть, по-видимому, понятие схемной вязкости не может быть перенесено на вихревые методы. Но с другой стороны, расчётная практика показывает, что влияние погрешности на вихревые структуры качественно аналогично вязким эффектам. Поэтому представляется целесообразным разработка аналогичного подхода, основанного на некоторых эвристических предположениях.

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

Анализ публикаций по теме исследования. Идея вихревых методов, в частности, метода дискретных вихрей, восходит к пионерским работам Гельмгольца и Кельвина (история вопроса подробно изложена в [5, 8]), однако только с появлением достаточно мощной вычислительной техники метод дискретных вихрей был надлежащим образом реализован алгоритмически, поскольку данный метод требует решения достаточно больших систем линейных алгебраических уравнений. Практически применение метода дискретных вихрей началось с работ С. М. Белоцерковского и А. Чорина, современное состояние вопроса достаточно полно отражено в их же монографиях [3-5, 8]. Как отмечалось выше, идея прямого моделирования физической вязкости схемной вязкостью для метода дискретных вихрей была особенно притягательна и неоднократно высказывалась [5], но надлежащего развития не получила. Ситуация практически аналогична и для комбинированного метода граничных элементов и дискретных вихрей, предложенного в работах [6, 7]. Хотя этот подход и точнее традиционного метода дискретных вихрей, следует отметить, что собственно для метода граничных элементов на основе специальной методики численного эксперимента были проведены достаточно обширные исследования точности [9], но они ограничились уравнением Лапласа, то есть, случаем потенциального течения, алгоритмы же метода граничных элементов, предназначенные для расчёта вязких течений, особой популярности не получили из-за низкой эффективности. Следует отметить, что в исследованиях по методу дискретных вихрей иногда за моделирование вязкости ошибочно принимают моделирование граничного условия прилипания по методу Лайтхилла-Чорина [8], но без учёта вязкости в вихревых зонах такой подход, разумеется, недостаточен. Тут необходимо привести объяснение популярности метода дискретных вихрей для расчёта течений при больших числах Рейнольдса: многие задачи обтекания, особенно обтекания неудобообтекаемых тел, автомодельны по числу Рейнольдса, поэтому значение физической вязкости, эффективной вязкости и схемной вязкости на структуру течения особого влияния не оказывают, хотя, конечно, процессы диссипации в рециркуляционных областях за телом для различных значений этих параметров будут существенно различны. Таким образом, представляется важным и нужным восполнить существующий пробел в исследовании погрешности вихревых методов вычислительной теории потенциала.

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

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

Основная часть. Рассмотрим формулировку задачи о плоском течении несжимаемой жидкости в терминах завихренность - функция тока. Уравнение для функции тока у

2 2 д у д у

-Т" +—Т" = _ю,

дx2 дy2

уравнение переноса завихренности ш для идеальной жидкости запишем в консервативной форме

дш д/чд/ч — + — (иш)+— и<ю) = 0,

дг дх дy

аналогично для вязкой жидкости

дш д ( \ д / ч --1--шш)+--^ш) = уДш,

дг дх ду

(1)

(2)

(3)

где, как обычно г - время; (х,у) - декартовы ортогональные координаты точки; компоненты скорости

u = -

ду

v = —

ду

дx '

V - кинематическая вязкость жидкости.

Отметим, что вводя диффузионную скорость [10]

gradш

ш

уравнение (4) может быть записано в виде, аналогичном (3)

ТТ + ^и + и0)ш) + -д^ + vo)ш) = 0.

дг дх ду

(4)

(5)

(6)

Согласно общему подходу вычислительной теории потенциала [11, 12] заменим уравнение (1) его граничноинтегральным аналогом

Ф0 )у(Xo )=/Ф0 (X,Xo ^ )Ъ (^-{у^ )дф° (X,X0 )Ъ (X )

Г

дп

+

+ jjш(X )Ф0 (X,Xo )dX,

Б

(7)

где Б - область решения; Г - ее граница; X, Х0 - точки источника и наблюдения, соответственно; функция С определяется как

IX х0 е Б,

С(Х0 ) =

-, Х0 еГ, 2 0

0, X й О, X й Г;

Ф0 - фундаментальное решение решение уравнения Лапласа в плоском случае

( \

1

Ф0 (Х,Х0 ) = ^1п 2л

д/(х " х0 )2 +(У " У0 )2

(8)

Рассмотрим последний интеграл в (7), для чего разобьем область решения Б на достаточно малые подобласти Б; ив каждой из этих подобластей введем точку X;, в окрестности которой разложим фундаментальное решение Ф0 (8) в ряд Тейлора по переменным (х,у). Тогда, полагая, что точка Xo не принадлежит подобласти Б;, составляющая последнего интеграла в (7), связанная с подобластью Б;, может быть представлена

| ш^)ф0 (^0 )К *Ф0 (ЗД ^ш^)« = Ф0 (X,Xo К, (9)

где очевидно ш; = |ш(x)dx .

Г

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

С другой стороны, если поставить в (7) аппроксимацию завихренности системой дискретных

вихрей

М

Ш(Х) = ^Ш15(Х,Х1), (10)

1=1

где 5 - дельта-функция Дирака, получим представление последнего интеграла в (7) в виде

М

|ш(х)фо (Х,Хо )ах =Х®1Ф0 (Х1,Х 0 ). (11)

Б 1=1

Сравнивая представления (11) и (9), можно заключить, что метод дискретных вихрей имеет очевидную связь с распределенной завихренностью, которую система дискретных вихрей аппроксимирует. Точность представления (9) зависит от расположения точки Х1 в подобласти Б1 [13].

В отличие от традиционного представления (10), представление (9) подразумевает возможность обратного перехода от дискретных вихрей к распределенной завихренности при помощи некоторой интерполяции.

При использовании представления (10) уравнение (2) вырождается в уравнение движения дискретных вихрей

X1 = и(хЬУ1), у 1 = у(хЬУ1), (12)

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

По теореме Синга [14], если в (10), (12) после некоторого промежутка времени т* в представлении (10) поменять знаки вихрей, то через время т* вихревая структура вернется в исходное

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

Рассмотрим теперь возможность моделирования физической вязкости при помощи схемной вязкости. Для вязкого течения вместо уравнения (2) следует решать уравнение (3) или следующее из него уравнение (6), которое отличается от уравнения (2) наличием диффузной скорости. Восстанавливая непрерывное распределение завихренности по дискретным вихрям, можно определить (с известными оговорками) диффузионную скорость в любой точке потока. Тогда для рассмотренного выше набора жидких частиц можно ввести Ах1 , Ау1, которые определяются как решения следующих дифференциальных уравнений

^ = VDx(Xl,Уl), ^ = VDy(Xl,Уl), (13)

ат х ат у

которые решим с нулевыми начальными условиями. Вектор 5x1,5у1 можно трактовать как решение аналогичной системы, с той лишь разницей, что в определении диффузионной скорости (5) физическая вязкость заменена схемной вязкостью, с теми начальными условиями. Поскольку величина диффузионной скорости VD пропорциональна кинематической вязкости V , а, следовательно, ей пропорциональна и длина вектора (Ах1, Ау1), то из отношения между векторами 5x1,5у1 и Ах1, Ау1 может быть получено отношение схемной вязкости и физической вязкости. Однако указанные вектора не обязательно сонаправлены, потому как аналог схемной вязкости целесообразно рассматривать только величину, определенную по смещению на величину (5x1,5у0 (Ах1, Ау1) вдоль последнего вектора в скалярном произведении. Перпендикулярную же составляющую смещения можно трактовать как аналог схемной дисперсии. Возможна и альтернативная трактовка, рассматривающая схемную вязкость как векторную величину, но она влечет дополнительные трудности.

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

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

где в левой части стоит приближенное распределенное значение завихренности, =Ш; /Би -осредненная завихренность в области Б;, моделируемая дискретным вихрем , К - количество ближайших к точке (х,у) дискретных вихрей, используемых для интерполяции. Определение по формуле (3) производных от завихренности, входящих в градиент в формуле (5), достаточно просто, хотя и несколько громоздко.

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

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

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

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

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

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

(14)

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

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

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

В целом же предложенная методика представляется универсальным средством совершенствования вихревых методов гидромеханики.

Литература

1. Самарский А.А. Теория разностных схем / А.А. Самарский. - М.: Наука, 1989. - 576 с.

2. Шокин Ю.И. Метод дифференциального приближения. Применение к газовой динамике / Ю.И. Шокин, Н.Н. Яненко - Новосибирск: Наука, 1985 - 364 с.

3. Белоцерковский С.М. Отрывное и безотрывное обтекание тонких крыльев идеальной жидкостью /

C.М. Белоцерковский, М.И. Ништ. - М.: Наука, 1978. - 352 с.

4. Белоцерковский С.М. Математическое моделирование плоскопараллельного отрывного обтекания тел /С.М. Белоцерковский, В.Н. Котовский, М.И. Ништ, П.М. Федоров. - М.: Наука, 1988. - 309 с.

5. Белоцерковский С.М. Моделирование турбулентных струй и следов на основе метода дискретных вихрей / С.М.Белоцерковский, А.С. Гиневский. - М.: Издательская фирма «Физико-математическая литература», 1995. - 368 с.

6. Yevdokymov D. V. Boundary element and discrete vórtices method for ideal fluid flow calculations. //

D. Durban and A. R. J. Pearson (Eds.) Non-linear singularities in deformation and flow. Proceeding of IUTAM Symposium held in Haifa, Israel, 17-21 March, 1997. Kluwer Academic Publisher. - P. 217 - 230.

7. Бразалук Ю.В. Применение комбинированного метода граничных элементов и дискретных вихрей для решения некоторых задач гидродинамического взаимодействия в плоских потоках / Ю.В. Бразалук, Д.В. Евдокимов, Н.В. Поляков // Вестник Харьковского национального университета. - 2003. - № 590. Сер. «Математическое моделирование. Информационные технологии. Автоматизированные системы управления». Вып. 1. - С. 55-60.

8. Chorin A. J. Vorticity and Turbulence. /A.J. Chorin - Springer-Verlag, Berlin, New York. 1994. - 312 р.

9. Бразалук Ю.В. Совместное применение метода малого параметра и метода граничных элементов для численного решения эллиптических задач с малыми возмущениями / Ю.В. Бразалук, Д.В. Евдокимов, Н.В. Поляков. - Вестник Харк. нац. ун-та., - 2005. - № 703. Серия «Математическое моделирование. Информационные технологии. Автоматизированные системы управления». Вып. 5. - С. 50-66.

10. Андронов П.Р. Вихревые методы расчета нестационарных гидродинамических нагрузок / П.Р. Андронов, С.В. Гувернюк, Г.Я. Дынникова. - М.: Изд. МГУ, 2006. - 184 с.

11. Бенерджи П. Метод граничных элементов в прикладных науках / П. Бенерджи, Р. Баттерфилд. -М.: Мир, 1984. - 494 с.

12. Бреббия К. Методы граничных элементов / К. Бреббия, Ж. Теллес, Л. Вроубел. - М.: Мир, 1987. -524 с.

13. Бразалук Ю.В. Расчет обтекания сложных гидродинамических конфигураций комбинированным методом граничных элементов и дискретных вихрей / Ю.В. Бразалук, Д.В. Евдокимов, В.Г. Решняк. - Вюник Дншропетровського ушверситету. Серiя «Мехашка». - 2012. - Вип. 16, Т. 1. - С. 50-67.

14. Synge J.L. On the motion of three vortices / J.L. Synge - Can. J. Math., 1949, V.1. - P. 257-270.