1612725465-542b3179b36a4849700e0b2ecf7da111 (828844), страница 22
Текст из файла (страница 22)
е.(8.53) на гидравлическом прыжке. Тогда во всех узлах un+1jсхема (8.100), (8.101), (8.106) сохраняет гидравлический прыжок.Д о к а з а т е л ь с т в о. Пусть j = j0 . Согласно формулам (8.94),(8.95), для собственного значения λn1,j0 +1/2 имеем выражение√()2u+uH1 + H2u1 + u212nλ1,j0 +1/2 =−− u1 u2 +.(8.117)222Учитывая соотношения (8.52), (8.53), получаем, что√H12 + 8H1 u21 + H12H122 H12= 2u1 H1=u1 u2 = u1= u1 √ 22H28H1 u21H1 + 8H1 u1 − H1=H1 + H2,2поэтомуλn1,j0 +1/2 = 0.(8.118)Тогда из формулы (8.112) следует, что[ ()]n0()n1ΛLux j0 +1/2 = 2cλ2 (Hu)x.j0 +1/2Кроме того, в силу равенства (8.52), имеющего место на гидравлическомпрыжке, получаем(Hu)nx,j0 +1/2 =H2 u2 − H1 u1= 0,hnпоэтому (ΛLux )j0 +1/2 = 0.151Очевидно, что unx,j+1/2 = 0 для всех j ̸= j0 .
Это равенство являетсяследствием условий (8.116). Таким образом, при всех j будет выполняться равенство()nΛLux j+1/2 = 0.(8.119)Тогда из уравнения шага предиктор (8.100) следует, что при всех j вектор потоков вычисляется по одной и той же формулеf ∗j+1/2 =f nj+1 + f nj,2вследствие чего уравнение шага корректор (8.101) запишется как( n)nun+1− unjf nj + f nj−11 f j+1 + f jj+−= 0,τh22или)− unjun+11( nj+f x,j+1/2 + f nx,j−1/2 = 0.τ2В силу равенств (8.92), (8.99), (8.119) имеемnf nx,j+1/2 = Anj+1/2 unx,j+1/2 = (RΛL)j+1/2 unx,j+1/2 = 0,поэтому un+1= unj , т.
е. гидравлический прыжок действительно остаjется на месте и не «размазывается» схемой предиктор-корректор.Замечание. При доказательстве теорем 8.1 и 8.2 мы не использовали конкретный вид функций θk (k = 1, 2). Это означает, что теоремыверны для произвольных схемных параметров θk .8.8. Схема Лакса – Вендроффа является частным случаем описанной выше схемы предиктор-корректор и получается из последнейпри θk ≡ 0 (k = 1, 2). Таким образом, в схеме Лакса – Вендроффа«предикторные» величины будут вычисляться с помощью разностныхуравнений (8.84), (8.106), в которых следует положить D = E, θk = 0:()f ∗j+1/2 − 21 (f nj+1 + f nj )+ Anj+1/2 f nx,j+1/2 − Gnj+1/2 = 0,τ /2)( nnnnHj∗ − 12 Hj−1+ Hj+1Hj+1− Hj−1+= 0.τ /22h(8.120)(8.121)Эти величины используются затем в уравнении второго шага (8.101).152Согласно замечанию, приведенному в предыдущем пункте, схемаЛакса – Вендроффа сохраняет гидравлический прыжок при h(x) =const и сохраняет состояние покоя жидкости при любом профиле дна.Однако, как и в случае системы линейных уравнений мелкой воды, онане сохраняет монотонность численного решения.
На рис. 33 изображены графики (штриховые линии) точного решение задачи о прорывеплотины (см. формулы (8.63), (8.64)). Все входные параметры взятытакими же, как в линейной задаче о прорыве плотины (см. пп. 7.4).На рис. 33, а показан профиль свободной границы (сплошная линия),полученный с помощью схемы Лакса – Вендроффа.
Видно, что в численном решении возникают осцилляции за фронтом бора, а также перед волной понижения. Напротив, схема предиктор-корректор (8.100),(8.101), (8.106) с переменными параметрами θk дает решение без «паразитических» осцилляций (см. сплошную линию на рис. 33, б).8.9. Противопоточная схема. Естественным обобщением противопоточной схемы (7.23), аппроксимирующей линейные уравнения мелкой воды, могла бы стать в нелинейном случае следующая схема:()n()n− unjun+1j+ A+ ux j−1/2 + A− ux j+1/2 = Gnj .τηη0.00.0-0.2-0.2-0.4-0.4-0.605-0.6010xа(8.122)5бРис.
33. Графики свободной границы при θk = 0 (а) и при заданиифункций θk по формуле (7.39) (б). x0 = 5, H1 = 1, H2 = 0, 5, t = 3,N = 10015310xЗдесь, в отличие от схемы (7.23), матрицы A± уже не являются постоянными, однако определяются они аналогично прежним формулам (7.21):( ± )n()nA j+1/2 = RΛ± L j+1/2 ,(8.123)при этом диагональные матрицы Λ± являются аналогами матриц (7.20) ( )nλ±0( ± )n1 j+1/2,(8.124)Λ j+1/2 = ( ± )nλ2 j+1/20где(λ±k)nλnk,j+1/2 ± |λnk,j+1/2 |, k = 1, 2,(8.125)2λnk,j+1/2 – собственные числа (8.94) матрицы Anj+1/2 .Осталось указать, как выбирается правая часть уравнения (8.122).Возьмем ее, например, в следующем виде:()0h(xj+1 ) − h(xj−1 )Gnj =, hx,j =.(8.126)n2hHj hx,jj+1/2=Но оказывается, что схема (8.122), (8.126) не будет сохранять состояниепокоя жидкости (см.
задачу 8.1). Следовательно, вектор Gnj необходимозадать каким-то другим образом.Получим теперь противопоточную схему способом, который мы уженеоднократно применяли ранее. В схеме предиктор-корректор параметры θk положим равными величинам θ0k из формулы (8.89). Тогда получаем, что()σ1 01DΛ = Σ, Σ =,æ0 σ2()|λ1 |0))h( +11( +2Λ − Λ− , RDΛ2 L =A − A− ,DΛ ==ææτ0|λ2 |где для краткости индексы n и j+1/2 опущены, σk = sgn(λk ), а матрицыA± определены по формулам (8.123).С учетом приведенных выражений, из «предикторного» уравнения(8.100) получаем]()n()1[nf ∗j+1/2 = f nj+1 + f nj − A+ − A− j+1/2 unj+1 − unj + h (RΣLG)j+1/2 .2154Тогда шаг корректор (8.101) примет вид()n()un+1− unj1 [ nj+f j+1 + f nj − A+ − A− j+1/2 unj+1 − unj +τ2h() ()n()n+h (RΣLG)j+1/2 − f nj + f nj−1 + A+ − A− j−1/2 unj − unj−1 −]n−h (RΣLG)j−1/2 = G∗j .Используя равенство (8.92) и представление A = A− + A+ , получаемокончательный вид противопоточной схемы()n()n− unjun+1j+ A+ ux j−1/2 + A− ux j+1/2 =τ(8.127)]1[nn= G∗j −(RΣLG)j+1/2 − (RΣLG)j−1/2 ,2nгде векторы Gj±1/2 и G∗j вычисляются по формулам (8.86) и (8.102)соответственно.Видим, что схемы (8.122) и (8.127) отличаются только правыми частями.
Преимуществом схемы (8.127) является то, что она сохраняетсостояние покоя жидкости. Это свойство следует из замечания, приведенного в конце пп. 8.7. В справедливости этого свойства можно убедиться и непосредственной проверкой (см. задачу 8.2).Для ровного дна схемы (8.122) и (8.127) совпадают. Численное решение задачи о прорыве плотины, полученное с помощью противопоточной схемы, не имеет осцилляций, но является более сглаженным, чемпри использовании схемы предиктор-корректор.ЗАДАЧИ8.1. Покажите, что для противопоточной схемы (8.122), (8.126) выполнение условия (8.109) теоремы 8.1 не приводит к выполнению тождеств (8.110).8.2.
Покажите, что для противопоточной схемы (8.127) выполнениеусловия (8.109) теоремы 8.1 приводит к выполнению тождеств (8.110),т. е. схема (8.127) сохраняет состояние покоя жидкости.155§ 9. Разностные схемы для задачгазовой динамики9.1. Уравнения газовой динамики. Рассмотрим течение идеального (нетеплопроводного и невязкого) газа в предположении, что в некоторой системе координат Oxyz движение газа происходит только вдольоси Ox и все параметры газа не зависят от других пространственныхкоординат y, z. Cистема уравнений, описывающая такое «одномерное»течение газа имеет следующий вид [18]∂u ∂f+= 0,∂t∂xx ∈ Ω.(9.1)Здесь t – время, Ω = (0, l) – область решения, u – вектор решения,f – вектор потоков,ρρu,ρu2 + pu = ρu , f (u) = (9.2)ρEρu(E + p/ρ)ρ – плотность газа, u – скорость, p – давление, связанное с температуройи плотностью газа уравнением Клапейронаp = ρRT,(9.3)R – газовая постоянная, R = cp −cv , cp – удельная теплоемкость газа припостоянном давлении, cv – удельная теплоемкость газа при постоянномобъеме, E – удельная полная энергия газа, равная сумме внутреннейe и кинетической u2 /2 энергий, при этом e = cv T .
Используя уравнение Клапейрона (9.3), мы можем выразить внутреннюю энергию черездавление и плотностьcv pcv ppe = cv T ===,(9.4)ρR(cp − cv )ρ(γ − 1)ρгде γ = cp /cv – показатель адиабаты, γ > 1. Тогда полная энергия выражается через скорость, давление и плотность по следующей формулеE=pu2+ ,(γ − 1)ρ2(9.5)E=c2u2+ ,γ(γ − 1)2(9.6)поэтому156где c – скорость звука,√c=γp.ρ(9.7)Недивергентная форма уравнений (9.1) имеет такой вид:∂u∂u+A= 0.∂t∂x(9.8)Здесь A = ∂f /∂u – матрица Якоби,01γ−3 2A=u(3 − γ)u222γ−2 3c3 − 2γ 2uc+u+u−γ−12γ−120γ − 1 .γu(9.9)Легко проверить, что она имеет три действительных собственных числаλ1 = u − c,λ2 = u,λ3 = u + c,(9.10)которые при условии c > 0 являются различными, поэтому при c > 0система уравнений (9.8) является гиперболической.Отметим, что уравнения (9.1) могут быть записаны и в других недивергентных формах, например, в виде системы уравнений∂v∂v+ Ae= 0,∂t∂xгдеρv = u ,pu ρAe = 0 u0 γp01/ρ ,uкоторая в покомпонентном виде выглядит следующим образом:ρt + uρx + ρux = 0,(9.11)ut + uux + px /ρ = 0,(9.12)pt + γpux + upx = 0.(9.13)Как и в линейном случае, количество краевых условий на концахотрезка [0, l] зависит от количества входящих в область Ω характеристик.
Например, если при x = 0 скорость положительна, но течение157дозвуковое, т. е. 0 < u < c, то λ1 < 0, λ2 > 0, λ3 > 0, и, следовательно,две характеристики входят в область Ω через ее левую границу x = 0,поэтому на этой границе надо задавать два краевых условия (например, скорость и давление). Если, например, на правой границе x = lтечение сверхзвуковое (u > c), то на этой границе все собственные значения положительны, поэтому все характеристики выходят из области,и краевых условий при x = l задавать в этом случае не надо.9.2. Характеристическая форма уравнений. Обозначим левыесобственные векторы матрицы A, соответствующие ее собственным значениям (9.10), через lk (k = 1, 2, 3), правые – r k .
Составим из этих векторов матрицы L и R, аналогичные (8.11):γ−1 2uc +u−c − (γ − 1)uγ−12γ − 1 u2uγ−1 ,(9.14)L= 1−(γ − 1) 2− 2 22cccγ−1 2−uc +uc − (γ − 1)uγ−121112c22c2u−cu+cuR= . (9.15)222c2c22 u2u1uuu1−+++224c2c 2(γ − 1)24c2c 2(γ − 1)Легко проверить, что для матриц L, R и A выполняются равенствавида (7.11), (8.12), (8.13):RL = LR = E,LAR = Λ,A = RΛL,(9.16)где Λ – диагональная матрица с элементами λk на диагонали.Умножим систему уравнений (9.8) слева на матрицу L и учтем равенства (9.16). В результате получим систему уравнений в характеристической форме∂u∂uL+ ΛL= 0,(9.17)∂t∂xкоторая в покомпонентном виде выглядит следующим образом:−ρc ut + pt + λ1 (−ρc ux + px ) = 0,()c2 ρt − pt + λ2 c2 ρx − px = 0,(9.18)ρc ut + pt + λ3 (ρc ux + px ) = 0.(9.20)158(9.19)В случае модели мелкой воды уравнение (9.17) записывалось покомпонентно в виде уравнений (8.15) для инвариантов Римана.