3.5. Расщепление потоков (1013317), страница 2
Текст из файла (страница 2)
Это позволяетрассматривать в дальнейшем самые разные формы этой зависимости, т.е.различные формы уравнения состояния.Для идеального газаp = ( γ − 1) ρ e = ( γ − 1)VV1 5,поэтому производные давления равныp,1 = ( γ − 1) e,p,2 = 0,p,3 = 0,p,4 = 0,p,5 = ( γ − 1) ρ(5.46)В этом случае,u u2 + γ −1 e( )∂FC vu=∂V wu u E + ( γ − 1) e) (ρ2ρuρvρw( ρ E + ρu20 ( γ − 1) ρ ρu00 00 ρuρ uv ρ uwγρ u 00+ p)00(5.47)Перемножаем полученные матрицы010002αβ − uu ( 2 − γ ) −β v −β w β ∂FC ∂FC ∂V−uvvu00==∂U∂V ∂U w0u0 −uw u (αβ − H ) H − β u 2 − β uv − β uw γ u (5.48)где введены обозначения:α=1 2u + v 2 + w2 ) , β = γ − 1(2(5.49)Здесь учитывается, чтоH =E+pρp,2) Дляρ= βeтого,чтобынайтисобственныезначенияматрицы∂FC,∂Uнеобходимо решить относительно λ следующее уравнение 5-ойстепени ∂Fdet C − λ I = 0 ∂U(5.50)Избавляя читателя от довольно длинной, но тривиальной, операциинахожденияопределителяматрицы( 5 × 5) , приведем лишьразмераокончательный вид этого уравнения:( u − λ ) ( u − λ )3где a = γpρ2− a2 = 0 ,(5.51)- скорость звука.Таким образом, собственные значения матрицы∂FCравны∂Uλ1 = u , λ2 = u + a, λ3 = u , λ4 = u , λ5 = u − a(5.52)Матрица Λ A , входящая в представление (5.9), имеет видu0ΛA = 00000 0(u + a )0 00u000 u00 00 0 0 ( u − a ) 0(5.53)Полученные значения собственных значений матрицы A имеют глубокийфизический смысл.
Фактически они характеризуют распространениевоздействия в потоке сжимаемой жидкости в направлении оси x. Любаяточка потока может распространять свое воздействие на окружающеепространство следующим образом:• со скоростью u за счет перемещения самого потока -вниз попотоку;• со скоростью u+a за счет перемещения самого потока ираспространения звуковых колебаний - вниз по потоку;• со скоростью u-a за счет распространения звуковых колебаний; еслискорость u превышает скорость звука a, это воздействие сноситсявниз и воздействие вверх по потоку отсутствует3) Определение собственных векторовДля этого необходимо решить систему( A − λI ) x = 0 ,(5.54)которая, как следует из (5.50), имеет бесконечное множество ненулевыхрешений.Для λ = u получаем:−u1000 2u ( 2 − γ ) −β v −β w β αβ − u−uvv000 x=0−uww000 u (αβ − H ) H − β u 2 − β uv − β uw β u (5.55)Нетрудно заметить, что только 2 строки этой системы являются линейнонезависимыми, поэтому получается система−u12 u (αβ − H ) H − β u0 x=0− β uv − β uw β u 00(5.56)Таким образом, для этого собственного значения существует тринезависимых решения.
Задаем произвольным образом 3 компоненты вектораx - x1 , x3 , x4 , а x2 , x5 находим из решения системы (5.56):x1 = 1; x3 = 0; x4 = 0; x2 = u; x5 = u 2 − αx1 = 0;x3 = 1; x4 = 0; x2 = 0; x5 = vx1 = 0;x3 = 0; x4 = 1; x2 = 0; x5 = wПолучаетсятрисобственныхвектора,соответствующиетремсобственным значениям λ = u :1 00 u 00 , r2 = 1 , r3 = 0 r1 = 0 0 01 v u2 − α w(5.57)Однако для нашей задачи удобнее будет использовать не эти векторы, анекие их линейные комбинации:(5.58)r1 + vr2 + wr3 , ρ r2 , ρ r3В результате получаем следующие собственные векторы1 u v , wα 0 0 ρ , 0 ρv 0 0 0 ρ ρw(5.59)Эти векторы и будут столбцами в матрице ( S A )−1 .Для λ = u + a получаем:1 −u − a2−β u − a αβ − u−uvv−uww u (αβ − H ) H − β u 2−β v −β wβ −a00 x=00−a0 − β uv − β uw β u − a 000(5.60)В этом случае имеем 4 линейно независимых строки и, соответственно,только одно решение и только один собственный вектор.Задаем x1 = 1 , а остальные компоненты x находим, решая систему (5.60):x2 = u + a; x3 = v; x4 = w; x5 = α + au +a2β(5.61)Получаем собственный вектор1u + avwa2 α + au + β (5.62)Для формирования матрицы ( S A )−1 более удобно использовать не этотвектор, а умноженный на1:2a 2 1 2 2au+a 2a 2 v 2a 2 w 2 2a α + au 1+22β 2a(5.63)Аналогичным методом получаем собственный вектор, соответствующийλ =u−a 1 2 2au−a 2a 2 v 2a 2 w 2 2a α − au 1+22β 2a(5.64)4) Построение матрицы ( S A )−1 осуществляется из собственных векторов(5.59), (5.63), (5.64).При этом должно учитываться расположение собственных значений λm вдиагональной матрице Λ A .Получается:112a 2u +cu2a 2v−1S=( A ) v2a 2ww2a 2α + uc 1+α2a 22β12a 2u −c002a 2vρ02a 2w0ρ2a 2α − ua 1+ρv ρ w2a 22β00(5.65)5) Обращение матрицы ( S A )−1 можно осуществить методом ГауссаЖордана или с помощью детерминанта и транспонированной матрицыалгебраических дополнений.Эти способы в данном случае примерно одинаковы по трудоемкости, нодостаточно тривиальны.
Опуская эти выкладки, запишем результатαββuβv βwβ − 2222 1 − a2aaaa αβ − ua a − β u − β v − β w β v1−000SA =ρρw1 −000 ρρ αβ + ua − a − β u − β v − β w β (5.66)Читателю предоставляется возможность самому убедиться, что этаматрица является обратной по отношению к ( S A )−1 .Упражнение.Доказать, что матрицы ( S A )−1 и S A , заданные формулами (5.66) и (5.65),являются взаимно обратными.5.4.Разностный аналог рассматриваемой задачиС учетом всех полученных результатов запишем конечно-объемноепредставление задачи, предложенной в начале данного параграфа.Напомним, что мы рассматривали задачу одномерноготеченияневязкого газа, описываемую уравнением (5.1):∂U ∂FC+=0∂t∂xРазностное представление этого уравнение в контрольном объемеописывается формулой (5.23)( FC )i +1/ 2, j ,k − ( FC )i −1/ 2, j ,k ∂U +=0∆x ∂t i , j ,kИ, наконец, конвективные потоки FC на гранях контрольного объемазадаются формулами (5.25):( FC )i +1/ 2, j ,k= A+U L + A−U R(5.67)Используем введенное в параграфе 4 обозначение приращение вектораU при переходе от n- го шага по времени к (n+1)-му шагу через δ U n+1 инеявное представление потоков в виде формулы (4.13)n ∂F FC = FC + α C δ U n +1 = FC n + α Anδ U n +1 , ∂U n(5.68)Для грани ( i + 1 / 2, j, k ) это представление с учетом расщепления потока поформуле (5.67) имеет видnn( FC )i +1/2, j ,k = ( A+U L )i+1/2, j ,k + ( A−U R )i +1/2, j ,knn +1nn +1+α ( A+ )i +1/2, j ,k (δ U L )i +1/2, j ,k + ( A− )i +1/2, j ,k ( δ U R )i +1/2, j ,k (5.69)Для конвективных потоков слева U L и справа U R используются формулы(5.26) - (5.28).В работе [3] показано, что можно использовать для U L и U R в явной частипредставления (5.69) аппроксимации 2-го или 3-го порядка (формулы (5.27)или (5.28) ), а для неявной части – аппроксимацию 1-го порядка, как болеепростую.
При решении задач методом установления такой подходгарантирует аппроксимацию 2-го или 3-го порядка соответственно попространственной координате.В этом случае получаем:( FC )i +1/ 2, j ,k = ( FC )i +1/ 2, j ,k + α ( A+ )i+1/2, j ,k δ U in, +j ,1k + ( A− )i +1/ 2, j ,k δ U in++1,1j ,k ,nnn( FC )i −1/ 2, j ,k = ( FC )i −1/2, j ,k + α ( A+ )i −1/ 2, j ,k δ U in−+1,1j ,k + ( A− )i−1/2, j ,k δ U in, +j ,1k ,nnnn(5.70)nгде для явных представлений ( FC )i +1/2, j ,k , ( FC )i −1/2, j ,k используются формулы(5.27) или (5.28).Подставляем (5.70) в исходную формулу (5.23) и получаем системуалгебраических уравнений для определения приращения δ U in, +j ,1k :A i , jδ U in, +j ,1k + Di , jδ U in++1,1 j ,k + Ei , jδ U in−+1,1 j ,k = ∆U in, +j ,1k,(5.71)где∆U in, +j ,1k = −∆t nn( FC )i +1/ 2, j ,k − ( FC )i−1/2, j ,k ∆x(5.72)приращение вектора U при переходе от n - го шага по времени к ( n + 1) -мушагу, представленное в явной форме,Ai, j = I +Di , jα ∆tn( A+ )i +1/ 2, j ,k −α ∆tn( A− )i −1/ 2, j ,k ;∆x∆xα ∆tα∆tnn=+( A− )i +1/ 2, j ,k ; Ei , j = − ( A+ )i −1/ 2, j ,k∆x∆x(5.73)Уравнение, вернее, система (5.71) по форме совпадает с системой (5.8) изГлавы 2, и для ее решения можно использовать эффективный методпрогонки.
Единственное, но существенное, отличие состоит в том, что вданном случае рассматривается не скалярные неизвестные величины, как вглаве 2, а векторные, и коэффициенты A i , j , Di , j , Ei , j являются не числами, аматрицами.Поэтому для решения системы (5.71)используется нескалярная, авекторная прогонка. Она будет описана в одном из следующих разделов.5.5.Двумерные и трехмерные задачи.Все, что говорилось здесь о расщеплении конвективного потока FC , легкораспространяется на конвективные потоки GC и HC , направленные по осямy, z соответственно.
Объясняется это тем, что основное уравнение (1.1)абсолютно симметрично относительно координат ( x, y, z ) .Например, конвективный поток через верхнюю грань контрольногообъема представляется формулой, аналогичной (5.25):( GC )i , j +1/ 2,k= B+U L + B−U R(5.74)В этой формуле:U L - некоторое значение вектора U ниже рассматриваемой грани, а U R -значение вектора Uвыше этой грани. Матрицы B+ , B− определяютсясоотношениями−1−1B+ = ( S B ) Λ B + S B ; B− = ( S B ) Λ B − S BМатрицы ( S B )−1 , S B определяются преобразованием подобия(5.75)B≡∂GC−1= ( SB ) Λ B SB∂U(5.76)Диагональная матрица Λ B состоит из собственных значений матрицыЯкоби∂GC:∂Uv0ΛB = 00000 0(v + a )00 0v 0000 v0 0Диагональные матрицы Λ B + , Λ B −0 0 0 ( v − a ) 0(5.77)содержат только положительные иотрицательные элементы матрицы Λ B соответственно.В результате получаются представления потоков GC на гранях, подобные(5.70).Так же поступаем с конвективными потоками HC на передней и заднейгранях контрольного объема (рис.1).Полученные выражения подставляются в основное разностное уравнение(4.8), и в результате получается система линейных алгебраических уравненийотносительно неизвестных величин – значений вектора U во всех узлахсетки.УпражнениеОпределить матрицы ( S B )−1 , S B , входящие в представление (5.76).УпражнениеОпределить матрицы ( SC )−1 , SC , входящие в представление матрицыЯкобиC≡∂H C−1= ( SC ) Λ C SC ,∂Uа также диагональную матрицу собственных значений Λ C ..