ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2015, том 55, № 8, с. 1356-1362
УДК 519.63
Посвящается светлой памяти А.П. Фаворского
МЕТОД АДАПТИВНОЙ ИСКУССТВЕННОЙ ВЯЗКОСТИ ДЛЯ РЕШЕНИЯ СИСТЕМЫ УРАВНЕНИЙ НАВЬЕ-СТОКСА1)
© 2015 г. И. В. Попов*, **, И. В. Фрязинов*
(* 125047Москва, Миусская пл., 4, ИПМ РАН;
** 115409Москва, Каширское ш., 31, НИЯУ МИФИ) e-mail: piv2964@mail.ru Поступила в редакцию 09.02.2015 г.
Рассматривается численный метод решения двумерных задач вязкого сжимаемого газа для системы уравнений Навье—Стокса. Предлагаемый численный метод основывается на методе адаптивной искусственной вязкости. Система уравнений Навье—Стокса в разностном виде аппроксимируется на неструктурированных сетках, а именно на треугольных или тетраэдральных элементах. Монотонность разностной схемы по критерию Фридрихса достигается введением в разностную схему слагаемых, содержащих адаптивную искусственную вязкость. Адаптивная искусственная вязкость находится из выполнения условий принципа максимума. В качестве численного эксперимента приводится серия расчетов внешнего обтекания цилиндра для различных чисел Рейнольдса. Библ. 5. Фиг. 2.
Ключевые слова: разностная схема, система уравнений Навье—Стокса, адаптивная искусственная вязкость, численный метод.
Б01: 10.7868/80044466915080141
1. ВВЕДЕНИЕ (ПОСТАНОВКА ЗАДАЧИ)
В работе предлагается новый численный метод решения системы уравнений Навье—Стокса для сжимаемого вязкого газа, основанный на ранее разработанном методе адаптивной искусственной вязкости для задач газовой динамики (метод АИВ, см. [1]—[4]).
В уравнениях Эйлера отсутствует физическая вязкость, в то время как в системе уравнений Навье—Стокса диссипативные слагаемые присутствуют. Поэтому при построении разностной схемы необходимо было обеспечить отсутствие искусственной вязкости в области пограничного слоя, где физическая вязкость существенна. Отсутствие искусственной вязкости в области пограничного слоя гарантирует правильность получаемых результатов в этой области и во всей расчетной области в целом. Данный подход к решению задач газовой динамики на основе системы уравнений Навье—Стокса был протестирован на основе численного моделирования отрывного течения несжимаемой вязкой жидкости в квадратной каверне с подвижной границей при различных числах Рейнольдса (см. [5]). В этой работе проведен численный эксперимент по внешнему обтеканию цилиндра для различных чисел Рейнольдса с целью определения диапазона работоспособности предложенного метода для сжимаемого газа.
Рассмотрим систему уравнений Навье—Стокса в декартовой системе координат в операторном виде:
— + Шу (ур) = 0,
д?
др-^ + Шу (ур Vа - цта) + gradа р = 0,
д?
— + Шу (у (Е + р)-ц (у, т)-^гаёГ) = 0,
1) Работа выполнена при поддержке РФФИ (код проекта 15-01-04620).
МЕТОД АДАПТИВНОИ ИСКУССТВЕННОЙ ВЯЗКОСТИ 1357
где р — плотность, р — давление, V = ч2) — скорость, pva — компоненты импульса, а = 1,2,
2 2
Е = рб + р -—--полная энергия, е — внутренняя энергия. Вектор та содержит компоненты
т ар, в = 1, Б, Б — размерность пространства. Также
/
т ар - 2
V
д V а . д^£
^дхр дх а) 3
— 8 ар ¿¡У v
можно трактовать как компоненты тензора вязких напряжений, а = 1,2, р = 1,2, ц — динамическая вязкость, X = ц——--коэффициент теплопроводности, Т = — =-р--темпера-
(У- 1)Рг Сг (у- 1)Скр
тура, Рг — число Прандтля. Эти уравнения решаются в области общего вида I > 0. Система уравнений замыкается уравнением состояния идеального газа р = (у - 1)рб, или давление задается из экспериментальных таблиц, у — показатель адиабаты Пуассона.
На границе области О и в начальный момент времени I = 0 задаются значения функций р, V и Е (или р).
2. АППРОКСИМАЦИЯ СИСТЕМЫ УРАВНЕНИЙ НАВЬЕ-СТОКСА
Исходную область О покроем произвольной треугольной сеткой, назовем треугольники ячейками расчетной сетки 0.к. Для построенных треугольников найдем центры описанных окружностей. Если центр лежит строго внутри ячейки, то назовем его центром ячейки и обозначим точкой х0 (см. фиг. 1). Если же центр описанной окружности лежит на границе или вне треугольника, то центром ячейки назовем ее центр тяжести. К центру ячейки отнесем все искомые значения сеточных функций р0, V 0, Е0, р0. В рассмотрение еще введем потоковые узлы Х1, получаемые как пересечение отрезка, соединяющего центры соседних ячеек и общего ребра или грани ячеек.
-о, -% центры описанных окружностей
- — или центры тяжести
- е
Фиг. 1. Сеточный шаблон для разностной схемы.
Рассмотрим аппроксимацию уравнения неразрывности. Для построения разностной схемы применим метод конечных объемов. Проинтегрируем дифференциальное уравнение неразрывности
др дг
+ ё1у (ур) = 0
по площади ю0, получим
1(1^ + Шу(ур))йъ = 0, (^ ю0 + |Шу(ур)^ю = 0,
®0
®0
\ а1у (ур) ¿ш = ( (п, у рУ1 = X(рУ} 0 +(ру^ г.
®0
После дискретизации по времени получим разностную схему
о+1 о Т^ (р^о, )0 + (pVП )р, .
Ро =Ро--X-о-е'.
ш^ 2
I
Остальные уравнения аппроксимируются аналогично. Запишем окончательный результат: — уравнения для импульса:
{р^о, ) 0 + (р V,, )0 , К )0 + (а ^ , Ц 0 + Й0,
(Р^а)0,+1 = (Р^а)00 X
1 (о, )0, ~ (о, )с
3 Ао,/Б
^ (оь Ха ) +
(Уа)0, - К)0 Ао,
Р0, + Р0 ( \ + —^-COS (о,, Ха)
где а = 1,2;
— уравнение полной энергии
[К р)0 + К р)0
ш Х
Е"°+1 = к -ШX
г.л < 4
+ V 0У 0, + 1 _У_ 2 2 у-1
V
Рк. + Ре
Р 0, Р0
Ц 0 + Ц 0
1 К, ) - (о, )2 1 (( - Щ ) + ( ( - V02 )
2
Ао,/Б
До,
X0 + X0, У 1 1 (р/р)0, - (Р/р)0
2 у- 1РгЯе
До,
где ^ — длина /-й стороны треугольника, До, — длина вектора х 0х 0 на нормальный вектор п, для 1-й стороны треугольника, Яе — число Рейнольдса. Здесь предполагается, что физическая вязкость ц = 1/ Яе, Б — размерность пространства, в данном случае Б = 2.
Предложенная разностная схема имеет малую диссипацию, и если в решаемой задаче есть разрывы решения, например присутствует ударная волна, то требуется введение дополнительной искусственной вязкости.
3. ИСКУССТВЕННАЯ ВЯЗКОСТЬ
Введем искусственную вязкость п в систему уравнений Навье—Стокса аналогично тому, как это было сделано для уравнений Эйлера (см. [1]—[4]), система примет вид
— + Шу (ур) = Шу (п grad р),
дг
+ div (ур Vа - цта) + gradа р = Шу ( grad р Vа),
дг
— + ё1у (у (Е + р)-Ц (у, та) - Xgrad Т) = Шу ( grad Е).
дг
МЕТОД АДАПТИВНОИ ИСКУССТВЕННОЙ ВЯЗКОСТИ 1359
В разностном виде диссипативный оператор с искусственной вязкостью п записывается в виде
п п , #0, - #0 й
(л. )о "¡1; Iп
Ап,
где # = р, pva, Е, а = 1,2.
Найдем искусственную вязкость п, для этого запишем разностное уравнение неразрывности в следующем виде:
п+1 п
Ро = Ро
1 -11 Ы)0 £, -11
ю 2 ю А П{
ю
I
Ап,
п,) о
п
Ро,
=А
Для того чтобы предложенная разностная схема была монотонной по критерию Фридрихса, потребуем, чтобы коэффициенты в разностном уравнении были положительны, т.е.
'(V п )
Ао = 1 -1г,-11 ^Ц, > о И В0 =п--
0 ^¿—1 П ' „ ¿—I А „ ' 0 А м
_ П ^п, ) о,
> о.
2 ю А п, , Ап, 2
г г 1
Разрешая эти неравенства относительно искусственной вязкости, для первого неравенства имеем
а для второго имеем
следовательно,
п, К'.),
п, ^ К' ■ ><1,
К' ■) •
Полученная искусственная вязкость используется и в остальных уравнениях.
в
4. ОБЛАСТИ ВВЕДЕНИЯ ИСКУССТВЕННОЙ ВЯЗКОСТИ
Если всюду вводить предложенную выше искусственную вязкость, то будет сильное размывание численного решения, поэтому будем вводить вязкость только на ударной волне (волне сжатия) и на осцилляциях численного решения. Для определения области, занятой ударной волной (УВ), используется неравенство
^ < о,
ее
. Ур
где I = —с- — единичный вектор по направлению градиента плотности, а ч, — проекция вектора
Ур|
скорости на направление вектора I. Другая область введения вязкости — область осцилляции численного решения (ОСЦ), когда значение функции плотности ро больше или меньше значений плотностей примыкающих к данной ячейке юо, т.е. проверяются два условия:
Ро > шах_(Ро,) или Ро < тт_ (Ро,).
о, *о,г=1,з о^о.м.з
Для уменьшения искусственной вязкости на ударной волне бралась минимальная искусственная вязкость, которая была получена для уравнений Эйлера (см. [4]), так как в области вне пограничного слоя физическая вязкость стремится к нулю и уравнения импульса и полной энергии переходят в пределе в уравнения Эйлера. Второе соображение введения искусственной вязкости по уравнению неразрывности для сжимаемого вязкого газа основывалось на отсутствии в этом уравнении членов, содержащих физическую вязкость.
Таким образом, искусственная вязкость по областям имеет следующий вид:
Б ),|, Х0 е ОСЦ, 0.5к| л/Б (V о) - ^ (Б (V о ),)2 0, х0 й УВ и ОСЦ,
42 + с2|, х0 е УВ,
где к = , Б — размерность пространства, в данном случае Б = 2, с — скорость звука.
5. ЭТАПЫ РЕШЕНИЯ ЗАДАЧИ
Метод адаптивной искусственной вязкости состоит из трех этапов решения системы уравнений Навье—Стокса.
На первом этапе (этап "предиктор") по значениям до = {ро, (р^а)о, а = 1,2, Ео, ро} в момент времени г = го, в отсутствие искусственной вязкости п = 0, по явной разностной схеме находим
~о+1 (~о+1 ---\о+1 1 т^о+1 ~о+1)
предикторное решение д = {р , (р уа) , а = 1,2, Е , р }.
На втором этапе по полученным предикторным значениям до+1 определяются области, в которых следует ввести искусственную вязкость и для найденных областей по значениям дп вычисляется искусственная вязкость п°+1 = ). Области введения искусственной вязкости описаны выше.
На третьем этапе (этап "корректор") по значениям до+1 и искусственной вязкости п°+1 находим все значения функций до+1 = {ро+1, (р)о+1, а = 1,2, Ео+1, ро+1} на момент времени г = го+1.
6. ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ
Тестирование и верификация описанного метода приводится в [5]. В рассматриваемом ниже численном эксперименте исследуется обтекание вязким потоком газа кругового цилиндра при различных значениях числа Рейнольдса. Тем самым исследуется вопрос, как предложенная методика работает при различных режимах течения. Для всех приведенных расчетов была взята неструктурированная треугольная сетка с числом узлов 151840 и числом расчетных элементов (треугольников) 301530. Сетка была адаптирова
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.