Anderson-et-al-1 (1185923), страница 65
Текст из файла (страница 65)
й 6.5. Методы решения уравнения потенциала 359 Разностное выражение (6.147) можно также записать в виде д (Р д ) их т/(РСС.С/2/аф/) =О, (6.148) если рс+1/2 определить следующим образом: РС+1/2 ( С) РС+!/2+ "СР/-1/2 (6.149) (6.150) где (/ и 1/ — контравариантные компоненты скорости, задавае- мые выражениями У= А/фь+ Аафч, 1/= Ааф2+ Аафи, (6.151) где А/=~я+~2, Аа=~хЧх+~аЧн Аа=Ч2+Чу, (6.152) д($, Ч) = д(„„) =зхЧа — $НЧ„ где значения в центрах ячеек рс+1/2 получаются из уравнения энергии (6.134), в которое входит только величина и, которая оценивается как (фа+1 — фс)/с5х. Из уравнений (6.148) и (6.149) видно, что добавление искусственной вязкости эквивалентно тому, что плотность берется с запаздыванием.
Джеймсон [Загпезоп, 1975] добавляет искусственную вязкость явным образом, тогда как в описываемой здесь схеме она включается посредством представления плотности в специальном виде. Если для искусственной вязкости и выбирается представление (6.144), то оба метода дают одинаковые результаты. Если т = О, то схема хорошо работает только в дозвуковых зонах и неустойчива в сверхзвуковых зонах. Однако если и равно положительной константе, то схему можно применять как для дозвуковых, так и для сверхзвуковых течений.
Следует заметить, однако, что в этом случае результирующая схема дает только первый порядок точности и является сильно диссипативной. Недавно были разработаны схемы второго порядка, в которых значения плотности берутся со смещением [5(е(про(1,,1аптезоп, 1981[. Прежде чем применять подход, основанный на введении искусственной сжимаемости или использовании значений плотности со смещением, преобразуем уравнения в вычислительной плоскости ($,11), что позволяет более просто ставить граничные условия. Если выполнить преобразование уравнений к строго дивергентному виду [1/!у(апс(, 1974), то уравнение (6.132) принимает вид 360 Гл. 6.
Численные методы решения уравнений течения а плотность определяется выражением р = ~1 — — ", ' Мт„(иф, + Ȅ— 1)1'л" ". (6.155) Необходимо знать также и метрические козффициенты. Они связаны с производными физических координат следующим образом: 6 =арч т) = таей (6.154) $в = — !кч, пя = Ухе. 6.6Л, Методы дая сверяавуиовых течений где Если плотность выразить из уравнения (6.153), а 0 и )~ в из уравнения (6.151), то изменение плотности запишется в виде Р т~ 8 + 8 ) ф~ р~г д дх а~~ д$ дч (6.156) где йф = ф — фь Лф; = ф~ — ф; ь Выражая плотность из уравне- ния (6.155) и У из (6.151), производную по $ в уравнении по- тенциала можно записать так: РРЯ'т Зе Р,ит, ао,.1 РсУ!Ую ~ д(ЬЮ 1 Р~Аа до~ 1, (6.157) а~~г~ ) дч 1 дч ~ Для демонстрации процедуры решения уравнения полного потенциала для сверхзвуковых течений выберем 5 в качестве маршевого направления, и пусть на слое 1 и всех предшествующих слоях нам известна вся информация.
Наша задача состоит в получении значения ф на слое 1+ 1. Первый член уравнения потенциала аппроксимируем в точке (1+ 1/2, 1), а второй— в (1,1) или в (1+ 1,1) для полностью неявной схемы. Рассмотрим первый член уравнения (6.150): д/дй(рУ/У). Чтобы получить решение в следующих точках по направлению $, разложим р в ряд относительно известных величин 1-го слоя. Так как р явно зависит от компонент скорости, запишем р=р(ф„, фя) и р=р,+Ар, (6.155) 6 6.6. Методы решения уравнении потенциала Если используется конечно-разностная аппроксимация первого порядка против потока д( ) 1 дз .
аз (( )с+! / ( )х !1' (6.158) то главный член ошибки аппроксимации есть (6.159) Это дает положительную искусственную вязкость, когда уз/А, > аз,. (6.160) д (Р у) д ( у (АзФЗ+ Азтч)1 что приводит к неявному алгоритму по маршевой координате. Если это разностное выражение записать через разность потен- циалов Лф, то д (РУ) (Р~+~Аз Ьф) (Рзе,Аз ддаф~ (Р",+|Аз дфз) (6.161) где р,", — плотность на и-й итерации в точке сетки 1+ 1.
На каждом шаге при продвижении по $ новое значение плотности должно определяться путем итераций. В большинстве задач оказывается, что этот процесс состоит из одной-двух итераций. Чтобы корректно ввести искусственную вязкость, шаблон при аппроксимации плотности смещается против потока (см. [Но1з1, 1979; На1ех е1 а1., 1979]). Значение ри, в уравнении (6.161) заменяется выражением +т з+ь у+из ° ° з+! г+1 )+цз +, Г+ з ( )г+1!+на+ (6.162) Последнее условие должно быть выполнено, если решается мар-' шевая задача для уравнения потенциала в случае сверхзвукового течения.
Если условие (6.160) не выполняется, искусственная вязкость становится отрицательной и это ведет к неустойчивости. Отметим, что возникающие в уравнении (6.157) смешанные производные аппроксимируются разностями против потока, зависящими от знаков коэффициентов при них, чтобы обеспечить преобладание диагональных членов в случае неявной схемы. Производная по з) для уравнения потенциала вычисляется в точке (+1: 362 Гл. 6. Численные методы решения уравнений течения где У1+1, 1+па > О У1+1, г+ця < О, ч'1+1,!+1/2=11 ( а ) ~ (6.163) Постановка граничных условий на поверхности тела почти такая же простая. Если система координат предполагается адаптированной к поверхности тела, то в двумерном случае на поверхности тела можно задавать условие типа У =А,ф + А ф„=О.
Для задания, граничных условий на поверхности (1 = 2) тела можно воспользоваться фиктивной точкой (1 = 1), расположенной под поверхностью. Пусть А А ( 1+1)1+1. 3 ( 1+1)1+1, 1 11+12 21 ( + 3 А1 — 1щ А 2 Ьч (6.165) Это выражение используется для исключения из уравнения на границс значения функции в фиктивной точке (/ = 1). В расчете с неявным граничным условием отсутствует ограничение на размер шага по маршевой координате, которое может возникнуть из-за наличия границы. К тому же неявное задание граничных условий ускоряет сходимость итерационного процесса нахождения плотности на каждом шаге по 5. причем 3 = О для Ус+1, 1+122 ) О и 3 = 1 для У1+ь 1+112 ( О. Такое представление плотности (уравнения (6.162) и (6.163)) вводит положительную искусственную вязкость во всех точках, где М1+1 > 1 (скорость потока сверхзвуковая).
Поскольку мы рассматриваем маршевую процедуру решения уравнения потенциала для сверхзвуковых течений, то метод оказывается несостоятельным, когда поток дозвуковой и искусственная вязкость становится отрицательной. Чтобы найти решение на слое 1 + 1, зададим на нем неявные граничные условия.
В любой плоскости симметрии можно воспользоваться простым отражением. Например, в двумерной задаче, когда делают проход по направлению 2), можно ставить это условие на границе (считаем, что симметричные граничные условия мы имеем при 1' = 1' „— 1), задавая 1+1' ~мах 1+1' ~мах (6.164) й 6.5 Методы решения уравнения потенпиала В конечном счете процесс вычислений сводится к решению трехдиагональной системы уравнений для Лф. В этом легко убедиться, рассматривая систему, записанную в виде [1+ яи д + 0 д (Ся)+ р д (Сз д )15ф где Сь Ст Са и Р— коэффициенты при членах с производными, а Я вЂ” известная правая часть уравнения.
В качестве упражнения читателю предлагается найти выражения для этих членов 0,35 0.30 О. 25 0.20 о 0.15 0.10 0.05 0.0 0.2 О.ст 0.6 0.3 1.0 1.2 Расстояния от поверхности Рип 6.23. Решение уравнения потенциала для обтекании клина; — уравнение полного потенциала (скаляр); — — — уравнения Эйлера. (см. задачу 6.15).
Решая трехдиагональную систему, находим Лф. Зная Лф, вычисляем рвьь Это новое значение плотности используется для расчета новой величины Ьф из уравнения (6.166), и этот итерационный процесс повторяется, пока не будет достигнута сходимость. В трехмерном случае следуют той же процедуре с той лишь разницей, что решение для Лф получают, используя приближенную факторизацию, что требует решения трехдиагональных систем по направлениям т1 и ~. Шанкар и Ошер 15йапкаг, Оз)тег, 1982) применилн описанную выше процедуру для решения уравнения полного потенциала.
Вместо простой линеаризации плотности по уравнению (6.155) здесь разлагают в ряд произведение плотности и контра- вариантной компоненты скорости по направлению $ в виде (р(у)„, =(ри), + л(ри). $ 6.6. Методы решения уравнения потенциала Помимо этого использовали разности вверх по потоку, взятые на характеристиках гиперболической системы уравнений. Такое представление производных вверх по потоку обеспечивает наиболее четкое применение искусственной вязкости, нежели искусственной сжимаемости или смещение шаблона вверх по потоку при аппроксимации плотности.