ветственно по нормали и касательной к /-й панели. Если rik - единичная нормаль к k-н панели, Akj в (14.7) вычисляется следующим образом:

Akftik-Vkj. (14.14)

После решения системы (14.6) вычисление (14.9) и (14.10) в контрольных точках заменяется на (14.11). Таким образом сравнительно трудоемкая операция вычисления \kj проводится лишь один раз. Другой трудоемкой частью рассматриваемой процедуры является решение системы (14.6).

Можно заметить, что решение для Ф, определяемое (14.2),. точно удовлетворяет уравнению (14.1) при N-oo. Плотности источников О/ выбираются так, чтобы выполнялось граничное условие обращения в нуль нормальной составляющей скорости (14.4) на поверхности тела. Для точек, удаленных от тела, логарифмические функции в (14.2) обеспечивают возврат решения к скорости потока Uoo, параллельной оси х, т. е. Ф = UooXk

Диагональные элементы в системе (14.6) велики, но система не является строго диагонально преобладающей. Однако, если число элементов превышает 1000, итерационные методы, рассмотренные в § 6.3, становятся весьма эффективны и рекомендуются для использования [Hess, 1975]. При малом числе элементов более эффективными оказываются прямые методы (§ 6.2). Поскольку матрица А полная, время решения системы (14.6) будет расти как O(N). Поэтому данная часть расчета будет занимать большую часть времени при больших Л.

Однако при заданном времени расчета панельный метод дает значительно более точное решение, чем конечно-разностные методы или метод конечных элементов, реализованные на соответствующей сетке, окружающей тело. Это означает, что при одинаковом времени выполнения в панельном методе имеется гораздо меньше неизвестных (плотности источников), чем неизвестных узловых значений в конечно-разностных методах.

14.1,2, PANEL: численная реализация

В данном разделе панельный метод будет применен для расчета течения у эллипса. Общая структура программы PANEL изображена на рис. 14.6, распечатки подпрограмм - на рис. 14.7-14.11.

Поверхность эллипса задается уравнением

x + {ylbf=\ или Jt = cos0, = fesin0, (14.15>

где b - длина малой полуоси. В подпрограмме BODY (рис. 14.8) формула (14.15) для определения координат тела используется с постоянным шагом по 6.

Чтобы иметь возможность производить по программе PANEL расчет дозвуковых невязких течений, по координате у проводится преобразование Прандтля - Глауэрта (п. 14.1.6), которое имеет вид yinc = y(l Течение несжимаемой жидко-

сти около поверхности (х, утс) является эквивалентным обтека-

PANEL Ввод данных

BODY

Координаты :х, к Контрольные точки: с Параметры: dSj. cos о:., sin с.

MATELM

Элементы ifjifkji

Рис. 14.6. Структура программы PANEL.

нию исходной поверхности потоком с числом Маха Моо вдали от тела.

Элементы матрицы Akj, входящей в систему (14.6), определяются в подпрограмме MATELM (рис. 14.9) по формулам (14.12) - (14.14). Матрица Akj хранится в массиве FNy/. Приращения касательной составляющей скорости, соответствующие источникам единичной плотности, хранятся в РТ/г, для дальнейшего использования в подпрограмме SURVL (рис. 14.10). Значения Rk, соответствующие (14.8), вычисляются в подпрограмме MATELM при несколько более общем предположении, что набегающий поток имеет две компоненты скорости (f/oo, Voo),

Для решения системы (14.6) используются подпрограммы FACT (рис. 6.15) и SOLVE (рис. 6.16), в результате чего определяются плотности источников а. Скорость и давление вычисляются в подпрограмме SURVL (рис. 14.10) по формулам, эквивалентным (14.11), но с использованием значений FT и FN.

Компоненты скорости и давления в точках (хр.ур), внешних ло отношению к телу, вычисляются в подпрограмме POINT <рис. 14.11) по формулам, эквивалентным (14.11) -(14.13). Ко-

2 С PANEL CALCULATES VELOCITY AND PRESSURES ABOUT

3 С AN ARBITRARY CLOSED BODY USING THE PANEL METHOD.

5 DIMENSION X(50),Y(50)ДС(50),YC(50),DS(50),FN(50,50)

6 1 ,ГТ(50,50) ,RHS(50) ,SDE(50) ,CI(50),SI(50) ,AA(50,50) ,Ш1(50)

7 COMMON X,Y,XC,YC,DS,FN,FT,RHS,PI,CPI,CI,SI

8 1 ,UINF,VINF,SDE

10 OPEN(1,FILE= PANEL.DAT )

11 OPEN(6, FILE= PANEL.OUT)

12 READ(1,1)N,IPR,UINF,VINF,FMN,B

13 1 FORHAT(2I5,4E10.3)

14 с

15 WRITE(6,2)N,B

16 WRITE(6,3)UINF.VINF, FMN

17 2 FORMAT(IX,PANEL METHOD WITH M2, ELEMENTS,SX,

18 1 ELLIPSE MINOR SEMI-AXIS =,F6.3,/)

19 3 FORMAr(IX,ONSET VELOCITY COMPONENTS = ЧгГб.З,

20 1 2X,FREESTREAM ИАСН NUMBER * ,F6.3, )

21 С

22 M=N+1

23 PI = 3.14159265

24 CPI=2.0/PI

25 С

26 С- CALCULATE COORDINATES OF BODY AND CONTROL POINTS.

27 С

28 CALL BODY(N,M,IPR,FMN,B) 2 С

30 С CONSTRUCT THE MATRIX EQUATION. Jl С

32 CALL MATELH(N.M,IPR)

J3 С

34 С TRANSFER FN INTO AA

35 С

36 DO 5 К = 1,N

37 DO 4 J = 1,N

38 4 AA(K,J) = FN(K.J)

39 5 SDE(K) = RHS(K)

40 С

41 С FACTORISE AA INTO L:U

42 с

43 CALL FACT(N,AA,IKS1)

44 С

45 С SOLVE FOR THE SOURCE DENSITIES, SDE(K)

46 С

47 CALL S0LVE(N,AA,IKS1,SD£}

48 С

49 С CALCULATE VELOCITY AND PRESSURE AT ТИБ BODY SURFACE

50 С

51 CALL SURVL(N,B,FMN)

52 С

53 С CALCULATE FLOW AT GIVEN POINTS.

54 С

55 CALL POINT(N,FMN) 5S STOP

57 END