Производной может быть заменен выражением

-Ч =(. -:)Ч . (18.91)

где Lt и Lx - односторонние двухточечные разностные операторы, приведенные после формулы (18.51). С точностью до О (Л/) оператор LAAqj- может быть заменен оператором

LAAq.--L; (л; Л + (/ (18-92)

Л+ = 0.5 (Л +1 ЛI), Л~ = 0.5 (Л -1 Л I).

(18.93)

Очевидно, что в зависимости от знака Л либо Л+, либо Л- будет равно нулю. Подстановка в (18.90) позволяет получить систему уравнений

A<7;+=AiRHS .

(18.94)

С точностью до 0(Af) (18.94) можно заменить системой

1 + Р А1;(л; -f £)][l + PAZ:( Af\ + )] A?; = А/ RHS

(18.95)

которая является LU-разложением, поскольку первый множитель является нижнетреугольным, а второй - верхнетреугольным. Уравнение (18.95) решается в два этапа:

а,;=[а/ rhs + р [Af + ) а,! ,]/ 1 + р (л;.

(18.96)

-К+р в- (М; I+i7) ч:;И+pj (И;

(18.97)

Уравнение (18.96) решается последовательно от левой границы в сторону уменьшающихся значений /. Уравнение (18.97) решается от правой границы в сторону увеличивающихся значений /. Можно заметить, что приближения, введенные для определения неявных членов, не влияют на стационарное решение, RHS = 0. Таким образом, LU-факторизация является работоспособным и экономичным алгоритмом расчета стационарных решений.

Для сжимаемых уравнений Навье -Стокса LU-факторизация может быть проведена следующим образом. Уравнение

(18.6) записывается в виде

Ж-+~у---Тх---W -

Для упрощения изложения предполагается, что цг = О и = 0. Приближенно факторизованное дискретное представление уравнения (18.98) может быть записано в виде

{1 + р ML {А - Р) ] [I + р MLy {В - Q} ] Aq = RHS , (18.99) где

RHS = LF - F0 + L,(G - G), (18.100)

А=, В = , Р = , Q =

Уравнение (18.99) аналогично алгоритму Бима - Уорминга (18.78), (18.79), за исключением лишь того, что в (18.99) а = = 0 и по-другому аппроксимируются вязкие члены, приводящие к смешанным производным. Это упрощает алгоритм и может быть использовано при определении стационарного решения методом установления, поскольку такое упрощение не влияет на стационарное решение.

Уравнение (18.100) можно также рассматривать как эквивалентное уравнение (18.90). Для проведения приближенной LU-факторизации уравнения (18.99) необходимо, используя разложение (18.60), выделить собственные числа матриц А и В. Из (18.61) следует, что значения собственных чисел могут быть или больше или меньше нуля. Следовательно, матрицы собственных чисел Лл и Ав можно разделить на положительную и отрицательную части:

Лл = А1 + ЛЯ, Ab = AJ+Ab, (18.101)

где, как и в (18.93) для скалярного случая. Ал = 0.5 {А + А 1} и Ал=0.5{Ад~Ал}.

В результате разложения (18.101) выражения LA и L/B в уравнении (18.100) можно заменить, как и в скалярном случае, суммой односторонних разностных операторов, связанных с положительными или отрицательными собственными числами. Например,

и аналогично для Ly. Поскольку F = Aq, можно провести и расщепление потока [Steger, Warming, 1981], в результате чего поток представляется в виде Р = F+ + F~. Данный вопрос кратко обсуждался в п. 14.2.5 в связи с расчетом сверхзвуко-

вых течений со скачками. В рассматриваемом алгоритме, следуя работе [Obayashi, Kuwahara, 1986], расщепление сводится лишь для проведения приближенной LU-факторизации неявных членов.

Благодаря приближенно неявному рассмотрению вязких членов можно провести также расщепление LxP и LyQ:

Lp = LxxP = {Lt - Lx ) P/Ajc, (18.102)

и аналогично для LyQ, Следовательно, в уравнении (18.100)

Lx {А - Р} Lx {IaAXIa + P/Ajc) - Lt (Тд Ад + PJAx),

(18.103)

В аналогичном виде представляется Ly{B - Q}. Однако в работе [Obayashi, Kuwahara, 1986] P/Ajc заменяется выражением Тдй1Тл\ где k выбирается из условия устойчивости, т. е.

= v/(RepAA:), v = max(2M Y!/Pr). (18.104)

Такое приближение Р напоминает скалярную форму (18.95). При этом приближении, а также используя приближенную факторизацию, можно провести дальнейшую факторизацию уравнения (18.99), в результате чего получим LU-форму

[i + PA/L;A] [i - PALA ] [i + PA/L;b] X

X [I - PA/LB~] Aq = А/ RHS (18.105)

А = Тл (лг + /1)Тл. (18.106)

Каждый множитель в (18.105) двухдиагональный и может быть разрешен относительно Aq за один проход в положительном или отрицательном направлениях х или у. Таким образом, весь алгоритм решения системы (18.105) может быть представлен в виде

[I + MLxk\ Aq* = At RHS\ [l - pA/LJ A ] Aq** = Aq*,

(18.107)

[I + PA/L;B+] Aq*** = Aq**, [l - pA/LB ] Aq = A***.

В качестве примера при решении второй системы осуществляется проход справа налево в направлении х. В каждой точке