1611257148-3fcbe9bca21f0b2c84e721e9e7e203b2 (Сиковский ДФ "методы вычислительной теплопередачи"), страница 6

В этом случаеmin u 0′ ( x )x −ξ= u 0 (ξ ) = u 0 ( x − ut ) ,tчтосовпадаетсрешениемуравненияРимана,полученного выше с помощью метода характеристик.При t > t b возникает ударная волна, в окрестности которой существенна вязкость. Ширинуобласти ударной волны δ можно определить из оценок величин членов уравнения:δ ~ ν max uПри ν → 0 толщина ударной волны стремится к нулю, и ударную волну можнорассматривать, как контактный разрыв.Контактный разрыв движется со скоростью, которую можно вычислить, проинтегрировавуравнение Бюргерса в малой окрестности разрыва:xs + 0xs + 0xs + 0xs − 0xs − 0xs − 0∫∂u∂dx =∂t∂t∫dx1udx − (u + − u − ) s = −2dt∫u 2 − u −2∂u 2dx + O(ν ) = − ++ O(ν )2∂xdx s u + + u −=dt2Разностные схемы для нелинейного уравнения переносаРассмотрим нелинейное уравнение переноса скалярной величины u∂u ∂f (u )+= 0,∂t∂x(*)где f (u ) - поток величины u, решение которого аппроксимируется на равномерной сетке.Введём среднее по контрольному объёмуxi +1 21ui =u( x, t )dx∆x xi∫−1 2Которое может быть связано с узловыми значениями ui с помощью какой-либо формулычисленного интегрирования45u i = ∑ Cij u jВ наиболее простом случае, соответствующем формуле прямоугольников,Cij = δ ijИз уравнения (*) следует∂ u i hi +1 2 − hi −1 2+= 0 , где hi +1 2 - т.н.

функции численного потока (numerical flux functions)∂t∆xОсновная цель метода конечных объёмов – правильная аппроксимация численныхпотоков. Такие аппроксимации могут быть записаны в общем виде, как(hi +1 2 = h u i − m +1 ,..., u i + m)(**)Аппроксимации (**) должны удовлетворять условиям согласованностиh(u,..., u ) = f (u )∂h∂f (u )== c (u ) , ∀u∂uj = − m +1 ∂ u jm∑uФункции численного потока можно определить, как 1) взвешенное среднее от узловых() ( )потоков f u i − m +1 ... f u i + m , 2) как поток в некоторой промежуточной точке u* , 3) какгибрид первых двух случаев.Например, для m=1 в (**) функции численного потока можно определить, как1) hi(+f1) 2 =f (ui ) + f (ui +1 )22) hi(+u1) 2 = f (u* ) , где u * = I (ui , ui +1 ) - промежуточная величина между ui , ui +1 , а I –интерполирующий оператор(3) hi +1 2 = H hi(+f1) 2 , hi(+u1) 2)46Схемы первого порядкаСхема Куранта-Изаксона-Риса (противопоточная)uin +1 − uin hi +1 2 − hi −1 2+=0,∆t∆xnnгдеhi +1 2 =f (ui ) + f (ui +1 ) 1− d i +1 2 ,22(*)где d i +1 2 - т.н.

диссипативный потокd i +1 2 = ci +1 2 (ui +1 − ui ) ,Величина ci +1 2 =ci +1 2∂fможет быть определена по следующим формулам∂u i +1 2⎧ f i +1 − f i⎪u − u ,⎪⎪ i +1 i= ⎨or⎪ ∂f⎪, u* ∈ [ui , ui +1 ]⎪⎩ ∂u u *или гибрид этих двух случаевСхема Годуноваhi +1 2⎧⎪( min ) f (u ), ui < ui +1, или= ⎨ ui ,ui+1f (u ), ui > ui +1⎪⎩ (maxu i , u i +1 )hi +1 2 =1⎡sgn(ui +1 − ui ) ⎡min f (u ) + max f (u )⎤ −max f (u ) − min f (u )⎤⎢⎣(ui ,ui +1 )⎥⎦⎥⎦(ui , ui +1 )(ui , ui +1 )2 ⎢⎣(ui ,ui +1 )2Если f (u ) - монотонная функция на интервале [ui , ui +1 ] , то схема Годунова аналогичнасхеме Куранта-Изаксона-Риса (*) с диссипативным потоком47⎛ ∂f ⎞d i +1 2 = sgn⎜ ⎟( f i +1 − f i )⎝ ∂u ⎠Схемы высокого порядкаСхема Лакса-Вендроффаnnuin +1 − uin hi +1 2 − hi −1 2+=0,∆t∆xd i +1 2 =hi +1 2 =f (ui ) + f (ui +1 ) 1− d i +1 2 ,22∆t 2ci +1 2 (ui +1 − ui ) ,∆xМетод Мак-КормакаМодификация метода Лакса-Вендроффа на основе схемы предиктор-корректор.Более прост в реализации, так как не требует вычисления матрицы Au~in +1 − uin f i +n1 − f i n+=0∆t∆x~~uˆin +1 − uin f i n +1 − f i −n1+1+=0∆t∆xun +1iu~in+1 + uˆin+1=2Схема Мак-Кормака хорошо описывает разрывы.

Лучше всего разрывы описываются,когда на шаге предиктор используются разности по потоку.Общий принцип применения аппроксимаций высокого порядка: для функциичисленного потока берётся линейная комбинация аппроксимации низкого порядка(обычно противопоточная) и аппроксимации высокого порядка. Аппроксимациючисленного потока высокого порядка можно представить в виде()hi +1 2 = hil+1 2 + hih+1 2 − hil+1 2 = hil+1 2 + hia+1 2где hia+1 2 - так называемый антидиффузионный поток, т.е.

снижающий численнуюдиффузию противопоточной схемы48Прямое применение аппроксимаций численного потока высокого порядка может привестик появлению нефизичных осцилляций в численном решении, что является следствиемнарушения монотонности схемы.Условие монотонности для решения исходной задачи заключается в том, что для двухзаданных начальных данных таких, что одно везде больше другого: u ( x,0 ) ≥ w( x,0 ) , далееэто свойство сохраняется: u ( x, t ) ≥ w( x, t ) .Для конечно-разностной схемыuin +1 = H (uin− m ,..., uin ,..., uin+ m )условие монотонности имеет вид∂H≥ 0 , j ∈ [− m, m]∂uin+ jДля линейных схемuin +1 = ∑ aij u njjиз этого условия вытекает неотрицательность коэффициентов aij .Согласно теореме С.К.Годунова (1959) не существует линейных монотонных схем спорядком аппроксимации выше первого. Для достижения более высокой точности можноослаблять свойство монотонности в пользу более слабого свойства сохранениямонотонности начального распределения.Это достигается в классе схем, обладающих так называемым свойством уменьшенияполной вариации (Total Variation Diminishing property), или TVD-схем.Этим свойством обладает исходное непрерывное уравнение, для которого справедливоTV ( u ( x, t2 ) ) ≤ TV ( u ( x, t1 ) ) , ∀ t2 ≥ t1где полная вариация определена, как49TV (u ( x, t )) =∞∫−∞∂u ( x , t )dx∂xДоказать выполнение TVD-свойства для решения уравнения Бюргерса∂ 2u∂u∂u=ν 2 ,+u∂x∂x∂tпредлагается в качестве задачи.Для конечно-разностных схем TVD-свойство записывается в видеTV (u n +1 ) ≥ TV (u n ) ,TV (u ) = ∑ ui +1 − ui = ∑ ViiОднимизiраспространённыхспособовпостроенияTVD-схемявляетсяметодограничителей потока (Flux limiter method) (Harten, 1984):(hi +1 2 = hil+1 2 + ψ i +1 2 hih+1 2 − hil+1 2)(*)в котором вводится функция-ограничитель потока ψ i +1 2 , меняющаяся в пределах от 0 до1, близкая к единице в областях с гладким изменением решения и близкая к нулю вобластях с сильными градиентами.

Таким образом, в областях с резкими изменениямивключаетсядиссипативнаяпротивопоточнаясхема,сглаживающаянефизичныеосцилляции.Обычно ограничитель потока выбирается в виде ограниченной функции от отношенияприращений решения в соседних узлахψ i +1 2 = ψ (ri ) , ri =ui − ui −1ui +1 − uiВеличина ri является индикатором различных типов поведения решения:1) ri ≤ 0 - область осцилляций2) 0 ≤ ri < 1 - область экспоненциального роста, характерного при приближении к разрыву3) ri ~ 1 и ri > 1 - области плавных изменений50Исходя из этого ограничитель потока должен иметь следующее асимптотическоеповедение:ψ (r ) = O (r ) , при r < 1ψ (r ) = 0 , при r < 0ψ (r ) = O (1) , при r ≥ 1Дальнейшие ограничения на вид функции ψ (r ) накладывают условие выполнения TVDсвойства. Так, для линейного уравнения переноса и противопоточной схемы для hil+1 2 исхемы Лакса-Вендроффа для hih+1 2 для выполнения TVD-свойства ограничитель потокадолжен удовлетворять неравенству (Davis, 1984; Roe, 1984; Sweby, 1984)0 ≤ ψ (r ) ≤ 2 min (r ,1)при условии справедливости условия Куранта ci +1 2 ∆t ∆x < 1 .Наиболее часто встречающиеся в литературе популярные ограничители имеют вид:Minmodψ (r ) = max[0, min (r ,1)]Superbeeψ (r ) = max[0, min (2r ,1), min (r ,2)]VanLeerψ (r ) = (r + r ) (1 + r )VanAlbadaψ (r ) = (r 2 + r ) (1 + r 2 )Влияние вязкости на устойчивость конечно-разностных аппроксимаций уравненияБюргерсаМетод ВВЦПuin +1 − uin f i n+1 − f i n−1uin+1 − 2uin + uin−1=ν+∆x 22 ∆x∆tАнализ устойчивости для линеаризованного уравненияuin+1 − uinun − unu n − 2uin + uin−1+ c i +1 i −1 = ν i +12∆x∆t∆x 251G = 1 − 2r (1 − cos β ) + iC sin βGФункцияописываетэллипс,которыйдолженнаходиться внутри круга:1) левый край не должен выходить за круг: r ≤122) радиус кривизны правого края должен бытьбольше 1Кривизна =r≤C2∂2x2,x=ReG≈1−rβ,y≈−Cβ⇒≤1∂y 22r1 C2,≤ 1 отсюда получаются довольно узкие границы устойчивости при C ≤ 1 :2 2rC21≤r≤ ,22Но даже при выполнении этих условий может наблюдаться потеря монотонностиC⎞⎛ C⎞⎛схемы.

Это следует из записи схемы в виде uin+1 = ⎜ r − ⎟uin+1 + (1 − 2r )uin + ⎜ r + ⎟uin−12⎠2⎠⎝⎝При∂uin+1Cc∆x≥ r , или Re ∆ =≥ 2 имеет место≤ 0 – немонотонное поведение. При2ν∂uin+1расчётах возникают осцилляции. Поэтому рабочий диапазон схемы ВВЦП сужаетсяC1≤r≤22Для неявной по диффузионным членам схемыuin+1 − uinun − unu n+1 − 2uin+1 + uin−+11+ c i +1 i −1 = ν i +12∆x∆t∆x 2остаётся только одно условие устойчивостиC2≤r2Метод Мак-Кормака для уравнения Бюргерсаu~in +1 − uin Fi n+1 − Fi nuin+1 − 2uin + uin−1=ν+∆x 2∆x∆t~~uˆin +1 − uin Fi n +1 − Fi n−1+1u~in++11 − 2u~in +1 + u~in−+11=ν+∆x 2∆x∆t52uin+1 =Устойчив при ∆t ≤u~in+1 + uˆin+12∆x 2u max ∆x + 2νМетод Кранка-Николсона с линеаризациейuin+1 − uin 1ν+ ∇ x Fi n + ∇ x Fi n+1 = ∇ 2x uin + ∇ 2x uin+122∆t()(Квазилинеаризация по Ньютону: Fiдля F = u 22(u ) + (u )n +1 2in 2i2n +1)⎛ ∂F≈ Fi + ⎜⎜ i⎝ ∂uinn⎞ n+1⎟⎟ ui − uin⎠()≈ uin u in+1uin+1 − uin u in+1uin++11 − uin−1uin−+11 ν ⎛ uin+1 − 2uin + u in−1 uin++11 − 2u in+1 + uin−+11 ⎞⎟⎟+= ⎜⎜+4 ∆x2⎝∆t∆x 2∆x 2⎠УРАВНЕНИЕ ТЕПЛОПРОВОДНОСТИ∂ 2T∂T=a 2∂x∂tВиды граничных условий (n– внутренняя нормаль):1) Условие Дирихле: T (0, t ) = b(t )2) Условие фон Неймана:∂T (0, t )= c (t )∂n3) Условие теплообмена с окружающей средой: λ∂T (0, t )= hw [T (0, t ) − Ta ] , где n –∂nвнутренняя нормаль к границе области.∂2Задача корректно поставлена если λ > 0 , hw > 0 , тогда − a 2 является положительно∂xопределённым эрмитовым оператором, имеющим положительные собственные значения(малые начальные возмущения будут экспоненциально затухать).53Явная схемаTi n+1 − Ti nT n − 2Ti n + Ti −n1= a i +1∆t∆x 2Ti n+1 = (1 − 2r )Ti n + rTi −n1 + rTi +n1 ,(Погрешность аппроксимации: O ∆t , ∆x 2r=a∆t– «сеточное» число Фурье∆x 2)Устойчивость: подставляем T ~ exp(λt + ikx )e λ∆t = 1 − 2r + 2r cos k∆x = 1 − 4r sin 2k∆x1– схема устойчива при 0 ≤ r ≤22Дифференциальное приближение:∂T∂ 2T a∆x 2∂ 4T ⎡ a 3∆t 2 a 2 ∆t∆x 2 a∆x 4 ⎤ ∂ 6T()−a 2 =+−+1 − 6r+ ...∂t∂x∂x 4 ⎢⎣ 31212360 ⎥⎦ ∂x 6при r =(1точность аппроксимации повышается до O ∆t 2 , ∆x 46)Численная дисперсия отсутствует, т.к.