3.7. Разностная форма уравнений Навье-Стокса и способы решения системы (Математическое моделирование задач газодинамики и тепло-массообмена. Молчанов А.М. 2013)

7. Разностная форма уравнений Навье-Стокса и способы решениясистемы7.1.Разностное уравнениеРазностная форма системы уравнений Навье-Стокса получается изосновного уравнения (4.8)δ U in, +j ,1k∆t+Fi +1/ 2, j ,k − Fi −1/ 2, j ,k∆x+Gi , j +1/ 2,k − Gi , j −1/ 2,k∆y+H i , j ,k +1/2 − H i , j ,k −1/ 2∆z=0(7.1)путем подстановки полученных в предыдущих параграфах представленийпотоков на гранях контрольного объема.Конвективные потоки (формулы (5.70) и аналогичные ей формулы):( 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  ,nnnnn+ α ( A+ )i −1/ 2, j , k δ U in−+1,1 j ,k + ( A− )i −1/2, j ,k δ U in, +j ,1k  ,n( FC )i −1/ 2, j ,k = ( FC )i −1/2, j ,k( GC )i , j +1/2,k = ( GC )i, j +1/2,k + α ( B+ )i, j +1/ 2,k δ U in, +j ,1k + ( B− )i, j +1/ 2,k δ U in, +j +11,k  ,nnnnn+ α ( B+ )i , j −1/ 2,k δ U in, +j −11,k + ( B− )i , j −1/ 2,k δ U in, +j ,1k  ,n( GC )i , j −1/ 2,k = ( GC )i , j −1/ 2,k( H C )i, j ,k +1/ 2 = ( H C )i, j ,k +1/ 2 + α ( C+ )i, j ,k +1/ 2 δ U in, +j ,1k + ( C− )i, j ,k +1/ 2 δ U in, +j ,1k +1  ,n( H C )i, j ,k −1/2(7.2)n(7.3)nnn= ( H C )i , j ,k −1/2 + α ( C + )i , j ,k −1/ 2 δ U in, +j ,1k −1 + ( C − )i , j ,k −1/ 2 δ U in, +j ,1k n(7.4)Вязкие потоки (формулы (6.17), (6.18) и т.п.):n( FV )i +1/2, j ,k = ( FV )i +1/2, j ,k + α Lni +1/2, j ,kn( FV )i −1/2, j ,k = ( FV )i −1/ 2, j ,k + α Lni −1/ 2, j , k1δ U in++1,1 j , k − δ U in, +j ,1k ) ,(∆x1(δ U in, +j ,1k − δ U in−+1,1j ,k ) ,∆x(7.5)1δ U in, +j +11,k − δ U in, +j ,1k ) ,(∆y1δ U in, +j ,1k − δ U in, +j −11,k ) ,(∆y(7.6)1δ U in, +j ,1k +1 − δ U in, +j ,1k ) ,(∆z1δ U in, +j ,1k − δ U in, +j ,1k −1 )(∆z(7.7)n( GV )i, j +1/2,k = ( GV )i, j +1/ 2,k + α N in, j +1/ 2,kn( GV )i, j −1/ 2,k = ( GV )i , j −1/2,k + α Nni , j −1/2, kn( HV )i , j ,k +1/ 2 = ( HV )i , j ,k +1/2 + α K in, j ,k +1/2n( HV )i , j ,k −1/2 = ( HV )i , j ,k −1/ 2 + α Kni , j , k −1/ 2В результате подстановки получаем основное разностное уравнение:Ai , j ,k δ U in, +j ,1k + Di , j ,k δ U in++1,1 j ,k + Ei , j ,k δ U in−+1,1 j ,k + Bi , j ,k δ U in, +j +11,k + Ci , j ,k δ U in, +j −11, k+ Fi , j ,k δ U in, +j ,1k +1 + G i , j ,k δ U in, +j ,1k −1 = ∆U in, j , k(7.8)где∆t ∆t nnnnFC )i +1/ 2, j , k − ( FC )i −1/ 2, j ,k  −GC )i , j +1/ 2, k − ( GC )i , j −1/ 2,k (( ∆y ∆x ∆t∆t nnnn− ( H C )i , j ,k +1/ 2 − ( H C )i , j ,k −1/ 2  −( FV )i+1/ 2, j ,k − ( FV )i−1/ 2, j ,k ∆z∆x∆t ∆t nnnn−GV )i , j +1/2,k − ( GV )i , j −1/ 2, k  −HV )i , j ,k +1/ 2 − ( HV )i , j , k −1/2  ,(( ∆z ∆y ∆U in, j ,k = −Ai , j ,k = I +α∆t (∆x α∆tnnA+ )i +1/2, j ,k − ( A− )i −1/ 2, j , k  − 2 ( Lni +1/ 2, j ,k + Lni −1/2, j ,k ) ∆xα∆tnn( B+ )i , j +1/2,k − ( B− )i , j −1/ 2,k  − 2 ( N in, j +1/ 2,k + N in, j −1/ 2,k )∆y∆yα∆t α∆tnn+C + )i , j ,k +1/ 2 − ( C − )i , j ,k −1/ 2  − 2 ( K in, j , k +1/ 2 + K in, j ,k −1/2 ) ,( ∆z∆z +(7.9)α∆t Di , j , k = +Ei , j ,k = −Bi , j , k = +Ci , j ,k = −Fi , j ,k =G i , j ,kα∆t∆xα∆t∆xα∆t∆yα∆t∆yα∆tn( A− )i +1/ 2, j ,k +n( A+ )i −1/2, j ,k +n( B− )i, j +1/ 2,k +n( B+ )i, j −1/ 2,k +n( C − )i, j ,k +1/ 2 +α∆t∆x 2α∆t∆x 2α∆t∆y2α∆t∆y 2α∆t2(7.10)Lni +1/ 2, j ,k ,(7.11)Lni −1/ 2, j ,k ,(7.12)N in, j +1/ 2,k ,(7.13)N in, j −1/ 2,k ,(7.14)K in, j ,k +1/2 ,∆z∆zα∆tα∆tn=−( C + )i, j ,k −1/2 + 2 K in, j ,k −1/2∆z∆z(7.15)(7.16)Здесь:∆U in, j ,k - приращения вектора U при переходе от n - го шага по времени к( n + 1) -му шагу, вычисленное явным методом;A i , j , k , Di , j , k , Ei , j , k , B i , j , k , Ci , j , k , Fi , j , k G i , j , k- блочные матрицы коэффициентов,каждый из них является матрицей размера ( 5 × 5 ) ;Матрицы A+ , A− , B+ , B− , C + , C − вычисляются по формулам (5.21),(5.75) и имподобной.Матрицы L, N , K вычисляются по формулам (6.15) при использованиипредположения о «тонком слое», либо по формуле (6.25) при использованиипредположения о «малости дивергенции скорости».Разностное уравнение (7.8) имеет точно такую же форму, какрассмотренные в главе 2 уравнения.

Оно содержит 7 неизвестных векторовδU .Вся система линейных алгебраических уравнений описываетсясемидиагональной матрицей блочных коэффициентов.Для его решения можно использовать методы, описанные в главе 2:итерационный метод Гаусса-Зейделя с переменой направлений, методприближенной факторизации, модифицированный метод приближеннойфакторизации Маккормака и т.п.Только надо иметь в виду, что коэффициенты в системе линейныхуравнений теперь являются не числами, а матрицами размера ( 5 × 5 ) .Поэтому при решении системы уравнений с трехдиагональной матрицейследует использовать не скалярную, а векторную прогонку.7.2.Векторная прогонкаМетод векторной прогонки аналогичен методу скалярной прогонки.Единственное существенное отличие состоит в том, что вместо деления начисло используется умножение на обращенную матрицу.Пустьврезультатефакторизацииуравнения(7.8)возникланеобходимость решения системы уравнений с блочной трехдиагональнойматрицейam Φm = bm Φm +1 + cm Φm−1 + d mгде m = 1, 2,..., N − 1, N .(7.17)Здесь неизвестная величина Φ и коэффициент d являются векторами,размер которых равен 5.

Коэффициенты am , bm , cm являются матрицами ( 5 × 5 ) .Для m = 1, m = N коэффициенты(7.18)c1 = 0, bN = 0где 0 - матрица, состоящая из нулевых элементов.Остальные коэффициенты в этих точках определяются из граничныхусловийЗаписанные условия (7.17), (7.18) означают, чтоΦ1известна взависимости от Φ2 . Уравнение для m = 2 представляет собой соотношениемежду Φ1 , Φ2 , Φ3 . Но поскольку Φ1 может быть выражена через Φ2 , этосоотношение приводится к соотношению между Φ2 и Φ3 . Другими словами,Φ2 можно выразить через Φ3 .

Процесс подстановки можно продолжать дотех пор, пока значение ΦN не будет выражено через ΦN +1 . Но поскольку ΦN +1не существует, мы в действительности на данном этапе получим численноезначение ΦN . Это позволит нам начать процесс обратной подстановки, вкотором ΦN −1 получится из ΦN , ΦN − 2 — из ΦN −1 , ..., Φ2 — из Φ3 и Φ1 — из Φ2 .Предположим, что при прямой подстановке имеем зависимостьΦm = Pm Φm +1 + Qm(7.19)Φm −1 = Pm −1Φm + Qm −1(7.20)после того, как полученоПодставляя (7.20) в (7.17), получаем следующее соотношение:−1Φm = ( am − cm Pm−1 ) bm Φm+1 + ( am − cm Pm−1 )( cmQm−1 + d m )Сравниваем это выражение с (7.19) и получаем(7.21)формулы дляпрогоночных коэффициентов Pm , Qm :−1Pm = ( am − cm Pm −1 ) bm ,Qm = ( am − cm Pm −1 )−1( cmQm−1 + d m )(7.22)Эти рекуррентные соотношения определяют Pm , Qm через Pm−1 , Qm−1 .Заметим, что в начале рекуррентного процесса уравнение (7.17) для m = 1 поформе почти совпадает с (7.19).

Таким образом, P1 , Q1 определяются вследующем виде:P1 = a1−1b1 , Q1 = a1−1d1(7.23)На другом конце последовательности Pm , Qm имеем bN = 0 . Это дает PN = 0 ,и из (7.19) получаем(7.24)QN = ΦNС этого момента осуществляется обратная подстановка с помощьюуравнения (7.19).Краткое описание алгоритма1. Рассчитываем P1 , Q1 из уравнения (7.23).2. Используя рекуррентные соотношения (7.22), получаем Pm , Qm дляm = 2,3,..., N3. Полагаем QN = ΦN4. Используя уравнение (7.19) дляΦN −1 , ΦN − 2 ,..., Φ3 , Φ2 , Φ1 .m = N − 1, N − 2,...,3,2,1 ,получаем.