Бураго Н.Г. Вычислительная механика (1185926), страница 26
Текст из файла (страница 26)
Рaсчeт сжимаемых теченийния;2) описание поведения решения внутри ячеек интерполяцией;3) расчет потоков на границах ячеек с использованием решенийзадачи Римана о распаде произвольного разрыва;4) явное интегрирование по времени уравнений газовой динамики нашаге по времени по какой-либо из схем Рунге-КуттаРассмотрим эти составляющие подробнее.16.4.1.
Консервативная аппроксимация законовсохраненияВ семействе схем Годунова используется следующаяинтегральная форма закона сохранения:∂YdV + ∫ G(Y)dV + ∫ (F(Y) + Y(u − w )) ⋅ ndS = 0∂t V∫*V*S*где V* - ячейка сетки (контрольный объем), Y - консервативнаяпеременная, G(Y) - источник (сток) консервативной переменной,F (Y) - диффузионный поток консервативной переменной,n-вектор внешней единичной нормали к границе ячейки, u - скоростьматериальной сплошной среды, w - скорость произвольныхподвижных координат (скорость сетки). Обычно интегралы погранице ячейки вычисляются на промежуточном временном слое поɶ , найденным из задачи Римана о распаде разрывазначениям Yрешения между смежными ячейками.Простейший вариант метода Годунова на произвольнойдвижущейся сетке получается при аппроксимации искомыхфункций постоянными по пространству значениями в каждой ячейкеи при использовании явной схемы первого порядка точности повремени.
Такая простейшая схема Годунова записывается в виде:(YV)in +1 − (YV)in+ (G(Y)V)in +∆t nm(i)+ ∑ [(F(Y) + Y(u − w )) ⋅ nS]nj +1/ 2 = 0j=1где нижний индекс указывает номер ячейки, верхний индексобозначает номер временного слоя, суммирование проводится пограням, ограничивающим ячейку, m(i) – число граней для ячейки i,суммирование проводится по границам данной ячейки с соседними154Глава 16.
Рaсчeт сжимаемых теченийячейками, дробный верхний индекс обозначает значения величин награнице ячейки в момент времени t n + ∆t n / 2 , полученные изрешения задачи Римана о распаде произвольного разрыва искомыхфункций Y на этой границе между ячейками. Описанный способразностной аппроксимации законов сохранения называется интегроинтерполяционным методом или методом контрольных (конечных)объемов.16.4.2. Расчет значений на границах ячеекРассмотрим задачу Римана (см. Куликовский, Погорелов,Семенов, 2001) для одномерной гиперболической системыуравнений∫ (Udx + F (U )dt ) = 0LU L при x < 0t = 0 : U = U 0 ( x) = U при x > 0 Rгде U и F - векторы плотностей и потоков сохраняемых величин, L– контур области интегрирования контрольного “объема” S вплоскости t-x:S={(t,x): 0 ≤ t ≤ ∆t , −∆x / 2 ≤ x ≤ ∆x / 2 }.Задача Римана имеет два элементарных решения.
Первое отвечаетдвижущемуся сильному разрыву, при этом решение определяетсясоотношениями на скачкахW [U ] + [ F (U )] = 0где W – скорость распространения разрыва, скачки величинобозначены квадратными скобками.Второе элементарное решение является непрерывнодифференцируемым и описывает распространение волн Римана. Этоавтомодельное решение ищется в виде U k = U k (ξ ) , где ξ = x / t иполучить его в общем случае нелинейных задач удается толькоприближенно (см. книгу Куликовского, Погорелова, Семенова2001).
Точное решение имеется для линейной системы уравнений.Оно является комбинацией бегущих волн со скоростями λ pnU p = ∑ Ω Rpk wk ( x − λk t ), p = 1,..., n.k =1155Глава 16. Рaсчeт сжимаемых теченийгдеnU = {U p }np =1 ,wk = ∑ Ω LpkU 0 p ,Ω R = {ΩA = ∂U F (U ) = Ω L ΛΩ R ,p =1nRpk p , k =1 ,}Λ = {λ( p )δ pl }np ,l =1-диагональнаяматрица,составленнаяизсобственных значений λ p матрицы A , Ω L и Ω R - матрицы левых иправых собственных векторов матрицы A .Реализация формул точного аналитического решения задачиРимана, если его удается получить, для вычисления потоков награницах ячеек связано с большими затратами вычислительнойработы.
В зависимости от начальных значений решения по обестороны от произвольного разрыва, общее решение задачи Риманапредставляется набором элементарных решений, описывающихдвижущиеся ударные волны, контактные разрывы и волныразрежения, разделенные областями с постоянными значениямиискомых функций. Для уравнений со сложной физикой, в частностидля уравнений, учитывающих влияние электро-магнитныхэффектов, фазовых превращений, сложной реологии, например,пластичности, аналитическое решение задачи Римана затруднено иво многих случаях не найдено. Поэтому предложено большое числовариантов метода Годунова (десятки вариантов, см. Куликовский,Погорелов, Семенов, 2001), основанных на приближенных решенияхзадачи Римана.Простейшим вариантом является метод Куранта-ИзаксонаРиса (КИР) (Courant, Isaacson, Rees, 1952), который основан наприближенномрешениизадачиРиманадлялокальнолинеаризованной исходной системы уравнений, состояшем издвижущихся разрывов, которые разделяют области с постояннымизначениями искомых функций.Метод Роу (Roe, 1981) основан на точном решении задачиРимана для специальным образом линеаризованной исходнойсистемы уравнений.
Отличие этого решения от используемого вметоде КИР состоит в том, что оно точно воспроизводит решениенелинейной задачи Римана для движущегося сильного разрыва.Метод Ошера (Osher, 1981) основан на приближенномрешении задачи Римана в виде комбинаций волн Римана.Заметим, что различия в модификациях схемы Годуновапроявляются только на нелинейных уравнениях.
Для линейнойсистемы гиперболических уравнений с постоянными коэффициентами, все модификации схемы Годунова одного и того же порядкаточности эквивалентны и совпадают с какой-либо из классическихконечно-разностных схем.156Глава 16. Рaсчeт сжимаемых течений16.4.3. Повышение порядка точностиПостроение схем более высокого порядка точности в методеГодунова достигвется путем сочетания использования кусочнополиномиальной аппроксимации решения внутри ячеек сразличными схемами Рунге-Кутта интегрирования по времени.Отметим, что чаще всего используются кусочно-линейнаяаппроксимация решения по пространству и двухшаговая схемапредиктор-корректор интегрирования по времени, обеспечивающиепочти второй порядок точности по пространству и времени (vanLeer, 1984; Borrel, Montagne, 1985; Родионов, 1987).Первый этап, называемый предиктором (“предсказателем”),реализуется по некоторой обычной конечно-объемной схеме наполушаге по времени:(YV)in +1/ 2 − (YV)in+ (G(Y)V)in +∆t n / 2m(i)+ ∑ [(F(Y) + Y(u − w )) ⋅ nS]nj = 0j=1Второй этап, называемый корректором ("исправителем"),представляет схему Годунова на полном шаге, в которойполученные на предикторе значения используются в качественачальных данных для решения задачи Римана о распадепроизвольного разрыва на гранях и последующего вычисленияграничных потоков:(YV)in +1 − (YV)in+ (G(Y)V)in +∆t nm(i)+ ∑ [(F(Y) + Y(u − w )) ⋅ nS]nj +1/ 2 = 0j=1Для дальнейшего повышения порядка точности методаГодунова используются более точные квадратурные формулы приаппроксимации интегралов и многошаговые процедуры Рунге-Куттаинтегрирования по времени.Явные схемы Годунова устойчивы при обычшых для явныхсхем ограничениях шага по времени.
Это ограничение погиперболичности, или КФЛ-условие, или условие Куранта,Фридрихса, Леви, или просто условие Куранта. Это условиеопределяет шаг по времени как отношение шага по пространству кскорости распространения возмущений.157Глава 16. Рaсчeт сжимаемых течений16.4.4. Расчет вязких теченийИзначальнометодГодуновапредназначендлягиперболических уравнений.
Поэтому при решении уравнений связкостью, которые являются параболическими, применяетсярасщепление по физическим процессам. На первом этапевычислений каждого шага по времени рассчитывается невязкаягиперболическая часть уравнений по схеме Годунова, а на второмэтапе отдельно по явной или неявной схеме производится учетвязких членов уравнений. При расчете вязких членов по явной схемеучитывается диффузионное ограничениеτn ≤h2( D + 1)!αгде D – число независимых пространственных переменных. Влюбом случае надо контролировать рост нормы решения с помощьюусловий точности.16.5.
Гибридные схемыКак правило, схемы высокого порядка точности вокрестности разрывов дают решения с нефизическимиосцилляциями в отличие от двухслойных схем первого порядка,которые могут быть монотонными в окрестности разрывов.Монотонные схемы при отсутствии источниковых членов непроизводят новых минимумов и максимумов решения (осцилляций).Идея почти монотонных гибридных схем повышеннойточности состоит в “скрещивании” немонотонных схем высокогопорядка точности с монотонными схемами низкого порядкаточности. При этом схемы высокого порядка используются вобластях гладкого решения, а монотонные схемы низкого порядкаточности применяются в окрестности разрывов.По определению Годунова двухслойная монотоннаялинейная схема имеет вид:Yin +1 =∑γYj= J(i)jnj, γj ≥ 0где J(i) – множество номеров узлов, образующих шаблон дляаппроксимации уравнения в узле i (множество соседей узла i).
Вмонотонной схеме положительное возмущение решения на старом158Глава 16. Рaсчeт сжимаемых теченийслое δYjn ≥ 0 ( j ∈ J(i) ) вызывает положительное возмущениерешения на новом слое δYin +1 ≥ 0 . Это свойство обеспечиваетсяположительностью коэффициентов γ j ≥ 0 .Для уравнения переноса∂Y / ∂t + U∂Y / ∂x = 0типичным примером монотонной схемы служит схема с разностямипротив потокаYin +1 = (1 − C)Yin + CYin−1 , C = ∆tU / hимеющая первый порядок точности.
Годуновым (1959) былопоказано, что не существует линейных двухслойных схем второгопорядка точности, обеспечивающих монотонное изменениерешения.Типичным примером немонотонной схемы может служитьдвухслойная одношаговая схема Лакса-Вендроффа второго порядкаточностиYin +1 = Yin − C / 2(Yin+1 − Yin−1 ) + C 2 / 2(Yin+1 − 2Yin + Yin−1 )которую можно переписать такYin +1 = Yin (1 − C2 ) + C / 2(C − 1)Yin+1 + C / 2(C + 1)Yin−1Эта схема устойчива при C ≤ 1 , но немонотонна, то есть допускаетобразование нефизических осцилляций решения вблизи разрывов.Для монотонизации, то есть для устранения нефизическихосцилляций, в схемах высокого порядка точности необходимоиспользовать нелинейную вязкость, зависящую от поведениярешения.Гибридные схемы являются нелинейными схемами, которыеобразуются при сложении с весом (параметром гибридности)монотонной схемы первого порядка и немонотонной схемыповышенного порядка точности (например,схемы ЛаксаВендроффа).