3.3. Искусственная вязкость (Математическое моделирование задач газодинамики и тепло-массообмена. Молчанов А.М. 2013)

3. Искусственная вязкостьДля того, чтобы разобраться с механизмом возникновения искусственнойсхемнойвязкости,рассмотримдвумернуюзадачуконвективногораспространения тепла при малых скоростях газового потока. Переносомтепла за счет теплопроводности пренебрегаем, а поле скоростей считаемзаданным.Уравнение энергии при этих условиях имеет вид∂∂∂( ρT ) + ( ρ uT ) + ( ρ vT ) = 0∂t∂x∂y(3.1)Для получения дискретного аналога этого уравнения используем методконтрольного объема. Для этого расположим узлы сетки, как на рисунке 1.

Ипроинтегрируем уравнение (3.1) по заштрихованному объему. Контрольныйобъем имеет форму параллелепипеда размером ∆x × ∆y × 1 .Рис.1. Контрольный объем (заштрихованная область) для двумерной задачиВ результате интегрирования получаем:∂( ρT ) dV + ∫∫ ( ρuT )i +1/ 2 − ( ρuT )i −1/ 2  ∆y + ∫∫ ( ρ vT ) j +1/ 2 − ( ρ vT ) j −1/ 2  ∆x = 0∂t ∫∫∫Индексы i, j нумеруют ячейку по осямx(3.2)и y соответственно. Величина( ρT ) задается в центре ячейки, а конвективные потоки ( ρ uT ) , ( ρ vT ) – насоответствующих границах ячейки, поэтому они имеют не целочисленныеиндексы.Считаем, что каждый поток постоянен по всей грани. В этом случаеуравнение (3.2) принимает вид( ρuT )i +1/ 2 − ( ρuT )i −1/ 2 ( ρ vT ) j +1/ 2 − ( ρ vT ) j −1/ 2∂+=0( ρ T )i , j +∂t∆x∆y(3.3)В соответствии с требованиями из предыдущего раздела будемпредставлять конвективные члены в неявной форме – на ( n + 1) -ом шаге повремени.

Встает вопрос о представлении потоков на гранях контрольногообъема, т.е. выражении этих потоков через значения в узловых точках. Напервый взгляд кажется, что оптимальным является представление в видесреднего между соседними узлами. Например,1( ρ uT )i + ( ρ uT )i +1  ,21= ( ρ vT ) j −1 + ( ρ vT ) j 2( ρ uT )i +1/ 2 =( ρ vT ) j −1/ 2(3.4)Однако в последующих разделах мы увидим, что это далеко не всегдалучший выбор.Поэтому используем так называемые параметры вверх по потоку. Еслипродольная скорость u ≥ 0 , то поток движется слева направо, и мы полагаем,что( ρ uT )i +1/ 2 = ( ρ uT )i ,( ρ uT )i −1/ 2 = ( ρ uT )i −1Аналогично при поперечной скорости v ≥ 0 поток направлен вверх, и(3.5)( ρ vT ) j +1/ 2 = ( ρ vT ) j ,( ρ vT ) j −1/ 2 = ( ρ vT ) j −1(3.6)В результате уравнение (3.3) представляется в видеΦ in,+j1 − Φ in, j∆tn +1+n +1n +1n +1( uΦ )i , j − ( uΦ )i −1, j ( vΦ )i , j − ( vΦ )i , j −1∆x+∆y=0,(3.7)где Φ = ρTПри решении задачи методом установления при приближении кравновесию Φ in,+j1 ≅ Φin, j , и уравнение (3.7) преобразуется к видуai , j Φ in,+j1 = bi , j Φ in−+1,1 j + ci , j Φ in,+j1−1(3.8)гдеai , j =1 n +1 1 n +11 n +11 n +1ui , j +vi , j , bi , j =ui −1, j , ci , j =vi , j −1∆x∆y∆x∆y(3.9)Рассмотрим случай смешения двух потоков, движущихся параллельно сравной скоростью (см.

рис.2), но с различными температурами в начальномсечении.Если бы теплопроводность λ не равнялась нулю (случай a), то междуэтими потоками образовывался бы слой смешения, в котором функцияΦ = ρT постепенно изменяется от большего значения к меньшему; ширинаэтого слоя в поперечном направлении увеличивается вниз по потоку.В рассматриваемом нами случае теплопроводность отсутствует: λ = 0 ,поэтому слой смешения не образуется, и будет иметь место скачок функцииΦ = ρT в поперечном направлении (случай b).Рис.2. Распределение функции Φ = ρT при наличии (а) и отсутствии (б)теплопроводности: 1 — горячий поток; 2 — холодный потокПостроим сетку так, чтобы ось x была направлена параллельно потоку.

Вэтом случае поперечная составляющая скорости равна нулю. Исходя изформул (3.9) и (3.8), получаем, чтоai , j = bi , j , ci , j = 0(3.10)Φ in,+j1 = Φ in−+1,1 j(3.11)иВ результате данные значения будут устанавливаться во всех точкахвдоль каждой горизонтальной линии. Следовательно, скачкообразныйхарактер профиля Φ и, соответственно, температуры, имеющий место вначальном сечении, будет сохранен. Не происходит размытия профиля и,следовательно, в этом случае искусственная вязкость не появляется.Ситуация значительно изменяется, когда та же задача решается на сетке,узловые линии которой составляют 45° с направлением потока.

Для удобстваиспользуем равномерную сетку с ∆х=∆у. Скорости потока в направленияхосей х и y равны. Исходя из формул (3.9) и (3.8), получаем, чтоai , j = bi , j + ci , j , bi , j = ci , jи(3.12)Φ in,+j1 = 0.5 ( Φ in−+1,1 j + Φ in,+j1−1 )(3.13)Рис.3. Случай потока, направленного под углом 45° к сеточным линиям:1 — горячий поток; 2 — холодный потокДля сетки, показанной на рис.

3, разрыв начального профиля Φобеспечивается равенством всех значений Φ = 100 вдоль левой и Φ = 0 вдольнижней границ. Результат решения во внутренних точках записан околокаждой узловой точки. Если схемная искусственная вязкость отсутствует, тосверху от диагонали, направленной из нижнего левого угла, должныполучиться значения 100 и значения 0 вниз от диагонали. Фактическоерешение дает неточный профиль Φ , подобный тому, какой изображен на рис.2, а.Исходя из этого примера, можно сделать следующие замечания по поводуискусственной вязкости [1].1. Искусственная вязкость имеет место, когда поток наклонен поотношению к линиям сетки и существует ненулевой градиент зависимойпеременной в направлении по нормали к потоку.2.

Приближенное выражение для коэффициента искусственной вязкостидвухмерном случае дано в [1]:Γa =ρU ∆x∆y sin ( 2θ )4 ( ∆y sin 3 θ + ∆x cos3 θ )(3.14)где U — модуль вектора скорости; θ — угол наклона (от 0 до 90º)вектора скорости к направлению оси х. Из этого уравнения видно, чтоискусственная вязкость не появляется, если результирующий потокнаправлен вдоль одной из сеточных линий; кроме того, искусственнаявязкость является максимальной, когда направление потока составляет угол45° с линиями сетки.3.

Вклад искусственной вязкости можно уменьшить, используя меньшиешаги ∆х и ∆y и располагая сетку (если это возможно) так, чтобы сеточныелинии более или менее совпадали с направлением потока.4. Поскольку реальная вязкость имеет место во многих задачах, тодостаточно сделать искусственную вязкость малой по сравнению с реальной.5. Основной причиной возникновения искусственной вязкости являетсяпрактика обращения с потоком через каждую грань контрольного объема какс локально-одномерным. Схемы, которые обеспечили бы меньший вкладискусственной вязкости, должны учитывать многомерную природу потока.Для этого также необходимо включать большее число соседних точек вдискретныйаналог.Внашемслучаеэтоточки( i − 1, j − 1) , ( i − 1, j + 1) , ( i + 1, j − 1) , ( i + 1, j + 1) .В заключение отметим, что искусственная вязкость носит не тольконегативный характер.

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