1D12.5) DO 25 К = 1,NY

DO 27 К = 1,NY

МЗ/ ITERATIONS THE RNS RESIDUAL IS

IF(IRD .EQ. 0)GOTO 30 DO 29 К = 1,NY VRITE(7,6)(U(K,J),J=1,NX) VRITE(7,6)(V(K,J),J=1,NX)

STOP END

Рис. 6.10 (окончание).

Выражения (6.13) вводятся в вектор R в следующем порядке:

-2,2 R2,2y Л%,3> Л2,3> R2,NY-U 2, ОТ-Ь 2> >

а порядок введения скоростей в вектор q таков:

2,2> 2.2> %,3> 2,3> > %,.УУ-Ь 2,ОТ~Ь 3, 2> 3,2>

В качестве начального решения относительно й и v было использовано точное решение уравнений (6.12) со значениями параметров ai = аг = 110.13, аз = а4 = О, as =1.0, А. = 5, Re = 10, Утах = я/30 и Хо= 1.0. Однако с помощью программы NEWTBU, примененной на сетке 5X5, т. е. с 18 узловыми неизвестными значениями и и не удается получить сходящееся решение. Для получения такого решения необходимо

DO 22 J = 2,NXP

DO 21 К = 2,NYP

ИВ = 2*(К-2) + 2*(NY-2)*(J-2)

DU = - RD(KB+1)

DV = - RD{KB+2)

U(K,J) = U(K,J) + DU*OM

V(K,J) = V(K,J) + DV*OK

CONTINUE

CONTINUE

GOTO 16

GENERATE OUTPUT

подвергнуть нижней релаксации поправку к q* Это означает, что применяется следующая формула:

q( +i) = q( )4-coAq(+i)

(6.19)

при со = 0.15. Однако и в этом случае скорость сходимости очень мала (см. данные на рис. 6.14). Такая ситуация не может считаться необычной при применении метода Ньютона к дискре-тизированным уравнениям, описывающим задачи гидроаэродинамики. А именно поправки, определяемые по методу Ньютона, могут оказаться достаточно большими для того, чтобы текущее решение вышло за пределы радиуса сходимости (см. п. 6.1.1), если только эти поправки вводятся без изменений.

Вид формулы (6.19) подсказывает соответствующую стратегию. Она состоит в том, что величина Aq , получаемая в результате решения уравнения (6.17), определяет лишь на-

правление поиска. Для трех значений со строятся новые реше-.ния

q( +i) = q( ) Aq( +i)

л вычисляются соответствующие им невязки Rm \ Для каждого значения озт рассчитывается среднеквадратичная невяз-

2 SUBROUTINE EXBUR(UE,VE,A,AL)

4 С CALCULATES THE EXACT SOLUTION OF THE TVO-DIMENSIONAL BURGERS

5 С EQUATIONS USING THE COLE-HOPF TRANSFORMATION

7 DIMENSION A{5),UE(21,21),VE{21,21),PH(21,21) 3 COMMON DX,DY,RE,NX,NY

9 PI = 3.1415927

10 XZ = 1.0

11 YMAX = PI/6./AL

12 ANY = NY - 1

13 DY = YMAX/ANY

14 ANX = NX - 1

15 DX = 2./ANX

16 DO 2 J 1,NX

17 AJ J - 1

18 X = -1. + AJ*DX

19 XD = X - XZ

20 DEX = EXP(AL*XD) + EXP(-AL*XD).

21 DDX = EXP(AL*XD) - EXP(-AL*XD)

22 DO 1 К = 1,NY .23 AK = К - 1

24 Y = AK*DY

25 SY = SIN(AL*Y)

26 CY = COS(AL*Y)

27 PH{K,J) = Ad) + A(2)*X + A(3)*Y + A(4)*X*Y + A(5)*DEX*CY

28 PHX = A{2) + A(4)*Y + A{5)*AL*DDX*CY

29 PHY = A(3) + A(4)*X - A(5)*AL*DEX*SY

30 UE(K,J) = - 2.*PHX/PH(K,J)/RE

31 VE(K,J) = - 2.*PHY/PH(K,J)/RE

32 1 CONTINUE

33 2 CONTINUE

34 RETURN

35 END

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

ка Rmmsy причем предполагается наличие квадратичной зависимости Rmmvs ®- Затем выбирается такое значение со, а следовательно, и Aq(*+>, которое ведет к минимальному Rrms яля заданного направления поиска Aq(*+>. Эта стратегия, однако, может оказаться слишком доростоящей в вычислительном смысле, если невязки Ru и Rv имеют сложные выражения.