1 SUBROUTINE rCT(NXH,ETA,ET2,DTR,Q,Dr,U)

3 С APPLIES FLUX-CORRECTED TRANSPORT ALGORITHM TO Q(J,K)

5 DIMENSION Q(101,3),DQ(101,3),U(101),DF(101,3),ADF(101,3)

7 С COMPUTE ANTIDIFFUSIVE FLUXES

9 EMU = ETA + 0.25*ET2*((U(1)+U(2))*DTR)**2

10 DO 1 К 1,3

11 1 ADr(l,K) EMU*(Q(2,K)-Q(1,K))

12 DO 3 J * 2,NXM

13 EMU ETA + 0.25*ET2*((U(J)+U(J+1))*DTR)**2

14 DO 2 К 1,3

15 ADF(J,K) = EMU*(Q(J+1,K)-Q(J,K))

16 С

17 С DIFFUSE THE SOLUTION

18 С

19 Q(J,K) Q(J,K) + DF(J,K) - DF(J-1,K)

20 2 DQ(J-1,K) = Q(J,K) - Q(J-1,K)

21 3 CONTINUE

22 DO 4 К 1,3

23 4 DQ(NXM,K) Q(NXH+1,K) - Q(NXM,K)

24 С

25 С LIMIT ANTIDIFFUSIVE FLUXES

26 С

27 DO б J s 2,NXM

28 DO 5 К - 1,3

29 S SIGN{1.0,ADF(J.K))

30 ADF(J,K) ABS(ADF(J,K))

31 DUM = S*DQ(J-1,K)

32 ADF(J.K) AMIN1(ADF(J,K),DUM)

33 DUM = S*DQ(J+1,K)

34 ADF(J,K) AMIN1(ADF(J,K),DUM)

35 ADF(J,K) = AMAX1(ADF{J,K),0.)

36 ADF(J,K) = S*ADF(J,K)

37 С

38 С ANTIDIFFUSE THE SOLUTION

39 С

40 Q{J,K) = Q(J,K) - ADF{J,K) + ADF(J-1,K)

41 5 CONTINUE

42 б CONTINUE

43 RETURN

44 END

Рис. 14.24. Распечатка подпрограммы FCT.

где i +i/2 = 0.5( /4-и, как правило, 110=1/6, tii=1/3, 112 = -1/6. Очевидно, что шаги (1) - (6) являются прямым обобщением описанного в п. 14.2.6 скалярного алгоритма с постоянными коэффициентами.

Диффузионные потоки на временном слое п определяются в программе SHOCK (рис. 14.17, строки 98 и 99). Все остальные шаги проводятся после вычисления q** и оформлены в виде

подпрограммы FCT (рис. 14.24). Различные параметры, используемые в подпрограмме ЕСТ, описаны в табл. 14.4.

Типичные результаты расчета распространения сильного скачка приведены на рис. 14.25. Использовались 100 узлов по пространству и Ах =0.01. Интенсивность скачка pi/p2 = 5.0 и

Таблица 14.4. Параметры, используемые подпрограммой FCT (и программой

SHOCK)

ф ♦ ♦ ♦

О метод коррекции потока

+ Лаке - Вендрофф (искусственная вязкость)

- точное решение

0.55

0.60

0.65

0.70

0.75

0.80

Рис. 14.25. Распространение сильной волны, Pi/pi = 5.0, у = 1.4.

7 = 1.4. При этих параметрах скорость распространения скачка г = 2.104. В начале расчета скачок располагается в точке л: = 0.501; после 100 шагов по времени с шагом А = 0.001 - в точке а? = 0.711. Это изображено на рис. 14.25.

Очевидно, что алгоритм FCT позволяет получить крутой профиль скачка с незначительными осцилляциями. Для сравнения на рис. 14.25 приведено также решение, полученное по схеме Лакса - Вендроффа с искусственной вязкостью v=1.0

В (14.54). Эти решения были получены после 200 шагов по времени, равных 0.0005, что обусловлено условием устойчивости. Хотя искусственная вязкость эффективно подавляет ложные осцилляции, очевидно, что скачок размазывается на большее число узлов. Это можно преодолеть путем перехода к более мелкой сетке, однако время счета при этом может стать неприемлемо большим. Если используется локально мелкая сетка, как правило адаптивная, чтобы избежать жесткого условия на А/, вся схема должна быть неявной.

Сравнение схем с искусственной вязкостью, схемы ЕСТ и схемы Годунова высокого порядка для одно- и двумерных стационарных и нестационарных течений с сильными ударными волнами и контактными разрывами проведено в работе [Woodward, Colella, 1984]. Контактный разрыв - это поверхность, при переходе через которую значение плотности претерпевает разрыв, а скорость и давление непрерывны. При развитии течения в ударной трубе граница между первоначально сжатым газом и газом низкого давления в дальнейшем становится контактным разрывом.

В работе [Woodward, Colella, 1984] обнаружено, что схема Годунова высокого порядка точности дает более точное решение, чем алгоритм ЕСТ, однако схема Годунова более слсж-на для программирования и требует большего объема вычислений. Как и можно было ожидать, алгоритмы ЕСТ оказались более точными, чем схемы с искусственной вязкостью. В работе [Woodward, Colella, 1984] не исследовались появившиеся позднее схемы ограничения потока (см., например,. [Yee, 1986]), которые можно рассматривать как развитие подхода ЕСТ. Можно ожидать, что решения, полученные по этим схемам, будут гораздо экономичнее, чем полученные по схемам Годунова, и, возможно, практически столь же точными.

Различные методы, от методов с искусственной вязкостью до использующих полную формулировку Годунова, требуют для реализации существенно различного времени. Для нестационарных течений, где важно правильно рассчитать и скорость, и интенсивность скачка, оправдано применение более сложных алгоритмов. Для стационарных течений с ударными волнами точное решение может быть получено и при помощи менее сложных алгоритмов.

14.2.8. Неявные схемы для уравнений Эйлера

В предыдущих разделах были рассмотрены явные схемы решения уравнений Эйлера, описывающих течения с сильными скачками или всюду сверхзвуковые течения, для которых в ста-