Бураго Н.Г. Вычислительная механика (1185926), страница 27
Текст из файла (страница 27)
Складывая записанные выше схемы первого и второгопорядков точности с весом 0 ≤ γ ≤ 1 , получаем простейшуюгибридную схему квазивторого порядка точности:Yin +1 = Yin [(1 − C) γ + (1 − C2 )(1 − γ )] ++C / 2(C − 1)Yin+1 (1 − γ ) + [(1 − γ )C / 2(C + 1) + γC]Yin−1159Глава 16. Рaсчeт сжимаемых теченийкоторая при γ = 0 имеет второй порядок точности и немонотонна, апри γ = 1 имеет первый порядок точности и монотонна.
В областимало меняюшегося решения применяют схему второго порядка( γ = 0 ), а в зоне больших градиентов переходят на схему первогопорядка ( γ = 1 ). Поскольку окрестности ударных волн занимаютмалую часть области решения, то записанная гибридная схема имеет“почти второй порядок точности”.Значения параметра гибридности зависят не только от нормыградинентов решения, но и от типа течения, то есть принимаютсяразличными в зонах разрежения/сжатия и в зонах дозвукового/сверхзвукового течения.Важный шаг в развитии гибридных схем был сделан вметоде коррекции потоков (Boris, Book, 1973), которые предложилизадавать коэффициент гибридности так, чтобы в численном решениине возникали новые экстремумы и чтобы уже существующиеминимумы не уменьшались, а максимумы не возрастали.
То естьповышение точности решения путем гибридизации надопроизводить так, чтобы не нарушить монотонность схемы.Гибридные схемы часто имеют вид двухшаговых схемпредиктор-корректор. При этом на предварительном шаге –предикторе (predict – предсказывать) используется заведомомонотонная схема, а на втором шаге – корректоре (correct –исправлять)используетсяантидиффузионнаясхемасотрицательным коэффициентом диффузии, устраняющая излишнююдиффузию, введенную предиктором. Для сохранения монотонностина корректоре антидиффузия ограничивается коэффициентами,называемыми лимитерами (ограничителями).В дальнейшем идея гибридной схемы, обладающейсвойством монотонности, была воплощена в семействах схем TVD,реализующих свойство “неувеличения полной вариации численногорешения” (“Total Variation Diminition”, Harten, 1983).
Далее былипостроены гибридные монотонные схемы высоких порядковточности семейства ENO. Подробное описание схем TVD и ENOдано в книгах Куликовского, Погорелова, Семенова (2001), Петрова,Лобанова (2007).16.6. Схeмы экспоненциальной подгонкиСамарский (1964) усовершенствовал схeму с нaпрaвлeннымирaзнoстями пeрвoгo пoрядкa тoчнoсти, уменьшив физичeскуювязкoсть:160Глава 16. Рaсчeт сжимаемых теченийAin +1− AiA i − A i −11A ni +1 − 2A ni + A ni −1=α2h1 + α* / αh2nτnn+Unгдe α * = 0.5 | U | h .
Вмeстo физичeскoй вязкoсти αиспoльзуeтся уменьшеннaя («пoдогнанная») вязкoстьα corr = αв рaсчeтe11 + α* / αВ рeзультaтe схeмa сoхрaняeт свoйствo мoнoтoннoсти и дaeт болееточное решение в пoгрaничных слoях. Идeя Сaмaрскoгo рaзвитa всeмeйстве схeм экспoнeнциaльнoй пoдгoнки [см. книгу Дулан,Миллер,Шилдерс,1980].Распространенныйварианткoррeктирoвaнной вязкoсти мeтoдa экспoнeнциaльнoй пoдгoнкидается выражениемα corr 1 α* 2 1 α* 4 α * α* = α cth ≅ α 1 + − + ...α α 3 α 45 α Направленные разности можно заменить центральными.
Надотолькоприэтомдобавитьискусственнуювязкость,компенсирующую дестабилизирующие члены, возникающие впервом дифференциальном приближении центрально-разностнойсхемы. Тогда схема экспоненциальной подгонки примет вид:Ain +1− Aiτnn+UnnAi +1 − Ai −1 1| U | h Ain+1 − 2 Ain + Ain−1= α+2h2 h2 1 + α* / αДальнейшее распространение схемы на нерегулярные сеткии конечно-элементную аппроксимацию решения по пространственным переменным основывается на следующем вариационномуравнении Галеркина∫ (( An +1− A ) / τ n + U ⋅∇An )δ AdV =nV= − ∫ αɶ∇An ⋅∇δ AdV + ∫ αɶ n ⋅ ∇Aδ AdSVSгдеαɶ = α11 + α* / α+|U | h2161Глава 16.
Рaсчeт сжимаемых теченийПриведенная запись справедлива и для пространственногослучая, скорость течения U при этом является вектором, а под hподразумевается характерный размер конечного элемента (ячейки).16.7. Схемы уравновешенной вязкостиДля расчета течений с сильными ударными волнами ипогранслоями, характеризующимися наличием узких зон сбольшими градиентами искомых функций весьма эффективными,робастными и простыми в реализации являются схемыуравновешенной вязкости или SUPG-схемы (Streamline UpwindPetrov-Galerkin method) [Hughes, Brooks, 1979, 1982; Tezduyar,Hughes,1982].Рассмотримприменениеэтихсхем кгиперболической системе уравнений3∂ t U + ∑ ∂ xi Fi = fi =1гдеU = ( ρ , ρ ux , ρ u y , ρ uz , ρ E )-векторконсервативныхгазодинамических зависимых переменных (плотность, компонентыимпульса и полная энергия) размерности 5, Fi - векторы проекцийпотоков на оси пространственных координат xi размерности 5, f вектор свободных членов размерности 5.
Неконсервативная формазаписи этой системы уравнений имеет вид:3∂ t U + ∑ A i ∂ xi U = fi =1гдеx = ( x, y, z ) = ( x1 , x2 , x3 )- независимые пространственныепеременные, A i = ∂Fi / ∂U - матрицы коэффициентов размерности{5x5}. Граничные и начальные условия имеют видU = U Γ (x, t ) на ГU × [0, T ]3∑i =1ni Fi = Fn* (x, t ) на Г F × [0, T ]U t =0 = U 0 ( x) в Ω162Глава 16. Рaсчeт сжимаемых теченийгде Γ = ΓU ∪ Γ F - граница пространственной области решения Ω ,части границы ΓU и Γ F , на которых заданы так называемыеглавные и естественные граничные условия, они не пересекаютсяΓU ∩ Γ F = ∅ , векторы Fn* и U Γ размерности 5, а также начальныезначения искомых функций U 0 ( x) являются заданными.В SUPG схеме решение ищется методом Петрова-Галеркинаи, соответственно, используется вариационная (конечно-элементная)формулировка задачи, то есть, иными словами, применяется вариантформулировки в смысле слабого или обобщенного решения.
Пустьвведена дискретизация области Ω конечными элементами Ω e(ячейками простой формы), тогда основное вариационное уравнениеимеет вид3i =1hhh∫ δ U ⋅ ∂ t U + ∑ Ai ∂ xi U d Ω +ΩNe+∑j =1Ne+∑j =1∫Ω( j )3∫∑Ω( j)3i =1τ ( j ) A i( j ) ∂ xiδ U ( j ) ⋅ ∂ t U ( j ) + ∑ A i( j ) ∂ xi U ( j ) d Ω +i =1D( j ) ∂ xi U ( j ) ⋅ ∂ xiδ U ( j ) d Ω = ∫ f ⋅ δ Ud Ω +Ω∫ F ⋅ δ Ud Γ*nΓFгдеτ ( j ) = max α h( j ) / || A( j ) ||1≤i ≤3ν d( j )= M ( j)∑ l =1i1/ 23 3( j)( j)A∂U⋅A (k j ) ∂ xk U ( j ) ∑ ∑ i xik =1 i =133( j)( j)( j)( j) ∑ ∂ xiφl ∂ xi U ⋅ ∑ ∂ xkφl ∂ xk U i =1k =1D( j ) = ν d( j ) Iздесь j – номер конечного элемента, I - единичная матрицаразмерности 5x5, φl( j ) - функции формы конечного элемента, M ( j ) число узлов в конечном элементе.
Коэффициент ν d( j ) имеет смысл иразмерность искусственной кинематической вязкости ("скоростьпомножить на длину"). Таким образом, в нелинейной схеме SUPGкоэффициентискусственнойвязкостиν d( j )определяется163Глава 16. Рaсчeт сжимаемых теченийавтоматически в зависимости от поведения решения по приведеннойвыше формуле, которая выражает "равенство" норм конвективных ивязких членов.
Отсюда следует и название схем такого типа "схемыуравновешенной вязкости".Отметим, что конечно-элементная реализация схемуравновешенной вязкости в случае равномерных прямоугольныхсеток соответствует центрально-разностной аппроксимации.Реализация схем возможна как с явной, так и с неявной повремени аппроксимацией правой части. При расчете вязких теченийможносовместитьсхемыуравновешеннойвязкостииэкспоненциальнойподгонки,аименно,вформулеэкспоненциальной подгонки корректированной физической вязкостидля центрально-разностных схем полагается при этом α * = ν d( j ) .16.8.
Неявные схемыЕсли недифференциальные по времени члены при записиразностных схем относятся к временному слою, отличному отстарого, то есть к новому или промежуточному слоям, на которыхзначения искомых функций еще неизвестны, то разностная,конечно-элементная, конечно-объемная (и так далее) схемыназываются неявными. В таких схемах искомые значения решенияна новом временном слое приходится находить из решения системуравнений с недиагональными матрицами.
Удовольствие этодостаточно дорогое, поскольку возрастает объем вычислительнойработы, а также затраты времени на программирование и отладкуболее сложного (по сравнению с явными схемами) алгоритма. Ноэти затраты нередко компенсируются повышенной устойчивостьюнеявных схем, позволяющей проводить устойчивый расчет сувеличенным шагом по времени.Несмотря на формальную безусловную устойчивостьполностью неявных разностных схем, шаг по времени при ихиспользовании все же приходится ограничивать условиямиточности, которые обычно сводятся к требованию малогоотносительного изменения нормы решения на шаге по времени.
Втех случаях, когда определяемый из условий точности шаг повремени становится сравнимым по величине с временным шагом,обеспечивающим устойчивый расчет по явным схемам,дляэкономии вычислительных затрат имеет смысл предусмотретьпереход к расчету по явной схеме. Этого можно не делать, если приреализации неявных схем воспользоваться эффективнымиитерационными методами. Для курантовских временных шаговчисло итераций для получения решения резко падает (до 1-3) инеявные схемы становятся столь же экономичными как и явныесхемы.164Глава 16.
Рaсчeт сжимаемых теченийТермин "эффективные итерационные методы" означаетметоды, приводящие к точному решению линеаризованныхуравнений неявной схемы за конечное число итераций. Примеромможет служить семейство методов сопряженных градиентов. Кнелинейным алгебраическим задачам метод сопряженныхградиентов применяется как внутренний итерационный процесс впаре с методом квазилинеаризации Ньютона. Эти методы подробнообсуждались в первой части данного курса.Итерационные методы имеют важную особенность, котораясостоит в том, что при наличии хорошего начального приближенияк искомому решению количество необходимых для отысканиярешения итераций снижается. В нестационарных задачах механикисплошной среды решение на старом временном слое являютсявесьма неплохим начальным приближением к решению на новомвременном слое, особенно, если шаг по времени достаточно мал.Поэтому при использовании итерационных неявных схем длябыстрых процессов эти неявные схемы по затратам вычислительнойработы нередко являются асимптотически эквивалентными явнымсхемам.Имеется большое разнообразие неявных схем.