1612725465-542b3179b36a4849700e0b2ecf7da111 (828844), страница 12
Текст из файла (страница 12)
Рассмотрим все возможные случаи.Пусть(4.44)unj > unj+1 .n+1nxjxj+1/2xj+1Рис. 15. Общая граница xj+1/2 соседних элементарных ячеек, на которой вычисляются потоки в схеме С. К. Годунова76Тогда точное решение задачи (4.5), (4.43) дается формулой (4.25), т. е.при t > tn имеем{unjпри x < xj+1/2 + D(t − tn ),u(x, t) =(4.45)nuj+1 при x > xj+1/2 + D(t − tn ),где D – скорость (4.23) движения разрываD=unj + unj+1.2(4.46)Нас интересует решение только на линии x = xj+1/2 .
В соответствиис формулой точного решения (4.45) имеем (рис. 16, а и б){unjпри D > 0,u(xj+1/2 , t) =(4.47)unj+1 при D < 0.Если скорость движения точки разрыва равна нулю, то такой разрыв называется стационарным скачком (рис. 16, в). Он возникает приусловииunj+1 = −unj < 0.(4.48)В этом случае в качестве u(xj+1/2 , t) можно взять любое из чисел unj+1или −unj . Для определенности положим u(xj+1/2 , t) = unj+1 .n+1nn+1xjxj+1/2nxj+1xjаxj+1/2xj+1бn+1nxjxj+1/2xj+1вРис.
16. Характеристики в задаче о распаде произвольного разрываnпри unj > uj+1 и D > 0 (а); D < 0 (б); D = 0 (в)77Пусть теперь выполняется обратное неравенствоunj < unj+1 .(4.49)Тогда решением задачи (4.5), (4.43) будет центрированная волна разрежения (4.26), т. е. непрерывная при t > tn функция nuj , x−xj+1/2u(x, t) =,t − tn un ,j+1x ≤ xj+1/2 + unj (t − tn ),xj+1/2 + unj (t − tn ) < x << xj+1/2 + unj+1 (t − tn ),(4.50)x ≥ xj+1/2 + unj+1 (t − tn ).Эта функция на линии x = xj+1/2 принимает следующие значения(рис.
17): n uj ,0,u(xj+1/2 , t) = nuj+1 ,n+1n0 ≤ unj ≤ unj+1 ,unj < 0 < unj+1 ,unj ≤ unj+1 ≤ 0.(4.51)n+1xjxj+1/2nxj+1xjаxj+1/2xj+1бn+1nxjxj+1/2xj+1вРис. 17. Характеристики в задаче о распаде произвольного разрываnnnnnпри unj < uj+1 и uj ≥ 0 (а); uj+1 ≤ 0 (б); uj < 0 < uj+1 (в)78Полученное решение u(xj+1/2 , t) и берется в качестве величины u∗j+1/2 : nuj ,еслиunj > unj+1иD > 0, un , еслиnnu>uиD ≤ 0,j+1jj+1 nuj ,еслиunj < unj+1и unj ≥ 0,∗(4.52)uj+1/2 =unj+1 , еслиunj < unj+1и unj+1 ≤ 0,0,если unj < 0 < unj+1 , nuj ,еслиunj = unj+1 .Итак, в схеме распада разрыва вначале вычисляется величина u∗j+1/2 ,∗затем потоки fj+1/2= f (u∗j+1/2 ), и, наконец, используется консервативная схема (4.30).Легко проверить, что в случае знакопостоянной сеточной функцииunj схема С.
К. Годунова совпадает с противопоточной схемой (4.42) (см.задачу 4.3).4.6. Обобщим схему предиктор-корректор (3.3), (3.4), построеннуюдля линейного уравнения переноса (3.1), на случай нелинейного уравнения (4.1). Как и во всех рассмотренных ранее консервативных схемах,∗потоки fj+1/2будем вычислять в полуцелых узлах xj+1/2 = xj + h/2:∗nfj+1/2− 12 (fj+1+ fjn )∗τj+1/2+ anj+1/2nfj+1− fjn= 0.h(4.53)∗nВ уравнении (4.53) τj+1/2= 0, 5τ (1 + θj+1/2), τ – шаг по времени,)(nnnnопределяется поθj+1/2 – параметр схемы, fj = f uj , функция aj+1/2формуле (4.35). Уравнение (4.53) есть результат аппроксимации уравнения для потоковft + a(u)fx = 0,(4.54)которое получается после умножения (4.1) на функцию a(u).
Первыйшаг называется «предиктором».Для того, чтобы схема была консервативной, на шаге «корректор»используем уравнение (4.30):∗∗fj+1/2− fj−1/2un+1− unjj+= 0.τhДля параметра θ = θ0,j+1/2 , гдеθ0,j+1/2 =1− 1;Crj+1/2Crj+1/2 = æ|anj+1/2 |;79(4.55)æ=τ= const,hполучаем противопоточную схему (4.40) (см. задачу 4.4). Если θ = O(h),то схема (4.53), (4.55) аппроксимирует уравнение (4.1) со вторым порядком относительно τ и h. В частности, при θ = 0 выписанная схемапереходит в схему Лакса – Вендроффа, которая дает нефизичные осцилляции на разрывных решениях. В общем случае (для произвольногопараметра θ) схема не сохраняет монотонность численного решения. Нооказывается, что этот параметр можно подобрать так, что схема будетобладать указанным свойством, т.
е. схему можно монотонизировать.В работе [29] на основе анализа дифференциального приближения схемы предиктор-корректор (4.53), (4.55) была предложена следующаяn:формула для сеточной функции θj+1/20g̃k − g̃k−sθ0,kθk =g̃kθ0,kпри |g̃k | ≤ |g̃k−s | , g̃k · g̃k−s ≥ 0,при |g̃k | > |g̃k−s | , g̃k · g̃k−s ≥ 0,(4.56)при g̃k · g̃k−s < 0,где s = sgn(anj+1/2 ), для сокращения записи верхний индекс n опущен,а через k обозначен дробный индекс j + 1/2,()g̃j+1/2 = |anj+1/2 | 1 − æ|anj+1/2 | unx,j+1/2 ;(4.57)θ0,j+1/2 =1 − æ|anj+1/2 |æ|anj+1/2 |.(4.58)Формула (4.56) является аналогом выведенной ранее формулы (3.34),поскольку при a = const (4.56) переходит в (3.34).Теорема 4.1.
При выполнении условияæ · max |anj+1/2 | < 1j(4.59)схема (4.53), (4.55) с переменным параметром θj+1/2 , заданным формулой (4.56), сохраняет монотонность численного решения.Д о к а з а т е л ь с т в о. Преобразуем выражение для потока[]nfj+1− fjn1 nnn∗+ fj − τ ak (1 + θk )fk =f=2 j+1h=]1[ n2fj+1 + fjn − τ (ank ) (1 + θk )unx,k =280=]1[ n22fj+1 + fjn − τ (ank ) (1 + θ0,k )unx,k + τ (ank ) (θ0,k − θk )unx,k =2()]1[ n= fj+1+ fjn − h|ank |unx,k + h gj+1 + gj − s(gj+1 − gj ) =2()]1[ n= fj+1+ fjn + h gj+1 + gj − s(ank + hγk )unx,k ,2гдеgj =)(1minmod g̃j−1/2 , g̃j+1/2 ,2(4.60)функция γj+1/2 задана формулой (3.45). При получении формулы для fk∗учтено, что числа ank и ank + hγk одного знака – это следует из оценки (3.48).
Значит, формулу для потока можно переписать такfk∗ =()]1[ nfj+1 + fjn + h gj+1 + gj − |ank + hγk |unx,k ,2поэтому уравнение (4.55) шага «корректор» запишется в видеun+1− unj1 [ nj+f+ fjn +τ2h j+1()+h gj+1 + gj − anj+1/2 + γj+1/2 hunx,j+1/2 −()]n−fjn − fj−1− h gj + gj−1 − anj−1/2 + γj−1/2 hunx,j−1/2 = 0илиun+1− unjj+τ)()1[ ( naj+1/2 + γj+1/2 h unx,j+1/2 + anj−1/2 + γj−1/2 h unx,j−1/2 −+2]−anj+1/2 + γj+1/2 hunx,j+1/2 + anj−1/2 + γj−1/2 hunx,j−1/2 = 0.Таким образом, схема предиктор-корректор (4.53), (4.55), (4.56) приведена к виду (2.22) с коэффициентами (3.49). Справедливость условий (2.23) теоремы 2.3 устанавливается так же, как в пп.
3.5 при доказательстве теоремы 3.1. Следовательно, при условии (4.59) рассмотреннаясхема будет сохранять монотонность численного решения.81ЗАДАЧИ4.1. Покажите, что выполнения условияæ · max unj ≤ 1(4.61)n,jдостаточно для устойчивости недивергентной схемы (4.38).4.2.
Докажите, что при условии (4.59) противопоточная схема (4.40)сохраняет монотонность численного решения.4.3. Докажите, что если unj > 0, то схема С. К. Годунова для уравнения Хопфа совпадает с противопоточной схемой (4.42).4.4. Покажите, что если в схеме предиктор-корректор (4.53), (4.55)выбран параметр θ = θ0,j+1/2 , где величина θ0,j+1/2 определена в формуле (4.58), то получается противопоточная схема (4.40).§ 5.
Схемы на адаптивной сеткедля уравнения переноса5.1. Ранее, при решении стационарных задач, мы уже использовали неравномерные сетки (см. § 2.7, § 3.7) и убедились, что точностьрасчетов может заметно возрасти за счет подходящего сгущения узловв подобластях сосредоточения особенностей решения. В настоящем параграфе мы продемонстрируем применение неравномерных сеток длячисленного решения нестационарных задач. В качестве примера возьмем начально-краевую задачу для линейного однородного уравненияпереноса с постоянным коэффициентом a > 0:ut + aux = 0,x ∈ (0, l],0 < t ≤ T,x ∈ Ω̄ ≡ [0, l],0 ≤ t ≤ T.u(x, 0) = u0 (x),u(0, t) = µ0 (t),(5.1)(5.2)(5.3)Пусть начальная функция является разрывной в точке x0 ∈ (0, l):{u1 ,x < x0 ,u0 (x) =u2 ,x > x0 .Если µ0 (t) ≡ u1 , то точное решение рассматриваемой задачи вычисляется по формуле u(x, t) = u0 (x − at) или{u1 ,x < x0 + at,u(x, t) =(5.4)u2 ,x > x0 + at.82u2u21.01.0110.50.50.00.001020300xа102030xбРис.
18. Графики точного (1) и численного (2) решений в моментвремени t = 10. θ = 0. Сетка равномерная, N = 60 (а); N = 600 (б)Вначале приведем результаты численного решения этой задачи наравномерной неподвижной сетке с узлами xj = jh (j = 0, . . . , N ) и шагом h = l/N .
Для поиска приближенного решения будем использоватьсхему предиктор-корректор (3.3), (3.4). При параметре θ ≡ 0 она совпадает со схемой Лакса – Вендроффа. Результаты расчетов для этогослучая представлены на рис. 18, на котором показано численное решение, полученное по схеме Лакса – Вендроффа, и точное решение (5.4)при следующих значениях входных данных:T = 10, l = 30, x0 = 10, a = 1, u1 = 1, u2 = 0.(5.5)Во всех численных экспериментах шаг по времени задавался по формулеhτ = kзап ,(5.6)aгде kзап – коэффициент запаса, который в расчетах полагался равным 0.5, что обеспечивало с запасом выполнение условий устойчивости (3.17).
Видно, что в численном решении имеются осцилляции, причемони остаются и на мелкой сетке. Кроме того, видно, что на мелкой сеткеположение разрыва передается точнее.На рис. 19 показано численное решение, полученное по схеме предиктор-корректор с параметром θ = θ0 , где постоянная θ0 определенав формуле (3.5). В § 3 было показано, что для такого значения параметра схема предиктор-корректор превращается в противопоточную83схему (3.7). Из рисунка видно, что противопоточная схема не дает осцилляций в численном решении, однако оно, в силу первого порядкааппроксимации, сильно «размазывается» в окрестности разрыва. Приизмельчении сетки «размазывание» численного решения уменьшается,и разрыв передается точнее.uu1.011.010.50.5220.00.001020300x10а2030xбРис.
19. Графики точного (1) и численного (2) решений в моментвремени t = 10. θ = θ0 . Сетка равномерная, N = 60 (а); N = 600 (б)uu11.011.00.50.5220.00.001020300xа102030бРис. 20. Графики точного (1) и численного (2) решений в моментnвремени t = 10. Параметр θj+1/2вычисляется по формуле (3.34). Сеткаравномерная, N = 60 (а); N = 600 (б)84xНа рис. 20 приведены результаты расчетов при использовании переnменного параметра θj+1/2, значения которого определяются по формуле (3.34). Видно, что, в отличие от схемы Лакса – Вендроффа, cхемапредиктор-корректор (3.3), (3.4) с переменным параметром θ сохраняетмонотонность численного решения.