термодиффузии. Для анализа численного решения уравнения (9.56) будем предполагать, что и и а - постоянные величины. Как и уравнение диффузии (гл. 7), уравнение переноса является строго параболическим (§ 2.3) и для него требуются граничные и начальные условия того же типа, например (7.2) -(7.4), как и для уравнения диффузии.

Однако при больших значениях и/а можно ожидать, что в уравмении (9.56) будут преобладать первые два члена. Тогда уравнение переноса будет проявлять свойства, аналогичные уравнению конвекции (9.2). Следует напомнить, что точные решения уравнения конвекции обычно имеют вид бегущей незатухающей волны (без уменьшения амплитуды). Следовательно, при больших значениях uja можно ожидать, что решения уравнения переноса (9.56) будут иметь вид волны со слабым затуханием.

На основе анализа уравнения конвекции с диффузией, которое является предельным для уравнения переноса, можно также ожидать, что в приближенных решениях уравнения (9.56) появятся прос-ьранственные волны, если будут использованы симметричные трехточечные алгебраические представления для конвективного члена при больших uja.

В этом параграфе мы рассмотрим алгебраические схемы (ВВЦП, Дюфорта - Франкела и т. д.), исследованные ранее для уравнений диффузии и конвекции (§ 9.1), для того, чтобы понять, вносит ли какие-либо новые трудности в численную схему одновременное включение диффузии и конвекции в уравнение.

9 А А. Явные схемы

Применение схемы ВВЦП для аппроксимации (9.56) приводит к алгебраическому уравнению

-It- +-ш---1;?-- -

которое приводит к вычислительному алгоритму вида

Tt =--{s + 0.5С) ти + (1 - 2s) Г? + (s - 0.5С) Tj+u (9.58)

где s = аМ/Ах и С = uAt/Ax.

Из разложения в ряд Тейлора в окрестности узла (/, п) следует, что схема (9.58) аппроксимирует уравнение (9.56) с погрешностью порядка О (At, Ах). Если алгоритм (9.58) рассматривать с точностью до членов порядка 0{At, Ах), то он

аппроксимирует уравнение

Ж + Ж- + - Т1 = (9-59)

т. е. в уравнение в частных производных, которое аппроксимируется алгоритмом (9.59), включается дополнительный член, интерпретируемый выше как главный член ошибки аппроксимации. Следуя методу модифицированных уравнений (§ 9.2), вторую производную по времени можно исключить, заменив ее производными по пространству, что дает

дТ , . дТ дТ

+ {O.SaAt-auAfi + 0.25uAfi-a + u- (9.60)

где а = Таким образом, видно, что использование раз-

ностной схемы первого порядка точности по времени вносит в уравнение переноса диссипацию и дисперсию первого порядка. Для больших значений и/а член искусственной диффузии адТ/дх может оказаться сравнимым с физической диффузией, если не ограничить величину М. Если величина дТ/дх не является слишком большой, то маловероятно, что последний член в (9.60) будет существенным. Следовательно, ограничение на шаг по времени At дает

А/< или С2<25 или Rceu[=-J)-- (9.61)

Ясно, что имеется веская причина использовать схему второго порядка точности по времени при решении уравнения переноса.

Напомним, что схема с разностями вперед по времени и центральными по пространству приводит к алгоритму, который условно устойчив для уравнения диффузии и всегда неустойчив для уравнения конвекции. В табл. 9.3 приведен коэффициент усиления G для уравнения переноса, полученный на основе анализа устойчивости схемы по Нейману. Для устойчивых решений G 1, что приводит к следующим ограничениям:

0<С2<25<1. (9.62)

Соотношение (9.62) включает в себя диффузионный предел (п. 7.1.1). Строго говоря, ограничение (9.62) допускает решения, для которых /?сеи = uAx/a = C/s > 2.0, т. е. можно ожидать появления осциллирующих решений. Из (9.61) также следует, что такие решения будут недостаточно точными.

Таблица 9.3. Алгебраические схемы (дискретизация)

Схема

Алгебраическая форма

Ошибка аппроксимации Е * (главные члены)

ВВЦП

л + 1

X (1 + 2с2)

Вверх по потоку

Д< Ах

-[Са.х-и()х X (1 -ЗС + 2с2)

1 дТ дх

Дюфорт - Франкел

(Г+ - Tl-)/2At+uLJ4-

Х[иАхУб-2аЮУи]-

Лаке - Вендрофф

= 0,

где а* = а + O.SwC Aj:

- [Са Ах-и (АхУб) (1 -С2)]Х X-§+lCaV (А/2)--аАА:2/12-мС (Алг/З) X

Х(С2-1)]

Кранк - Николсон

д7- +

и (АхУб) (1 + 0.5С2)

-а(Дл:2/12)(1 + ЗС)