а компонента вектора G задается формулой

1 1

gm=- \ \фт\(1хйу. (5.102)

-1 -1

Параметр т соответствует номерам всех внутренних узлов т. е. 2jNX-\, 2kNY-L Интегралы в (5.101) лучше всего вычислять для каждого элемента по отдельности используя элементные координаты (g, т]), как в (5.64). Нулевые вклады в глобальные интегралы в правых частях (5.101) и (5.102) вносят лишь по четыре соседних прямоугольных элемента, окружающих узел т.

Структура выражений (5.101) -и (5.102) аналогична структуре выражений (5.77) и (5.78), если не считать того, что интегрирование теперь ведется по двум измерениям вместо одного. Интерполяционные функции ф1{х, у) оказываются такими, что вклады в интегрирование могут быть определены в форме произведений одномерных интегралов (см. приложение А.2).

В случае однородной сетки уравнение (5.100) приобретает особенно компактную форму

{тУ ® + x<S> Lyy] , = -1, (5.103>

где индексы /, k соответствуют глобальной сетке (рис. 5.11). Операторы Lxx и Lyy в уравнении (5.103)-это направленные разностные операторы

L = (l/Ax2, -2/Ах 1/Ах\

Lyy = my\ -2/А2, 1/АуУ,

а Мх и My - направленные массовые операторы

М, = м1 = {\/6, 2/3, 1/6); (5.105)

наконец, символ ® обозначает тензорное произведение. Операторы имеют следующий смысл:

LxxWf k = -2- .

Myw k = (1/6) k-l + (2/3) w k + (1/6) w

/,fe=-6-t----

My(S>LxxWi

(5.106>

Комплекс Мх ® LyyWj, k имеет аналогичную структуру.

Трехточечная конечно-разностная дискретизация уравнения (5.97), эквивалентная данной, имеет вид

(5.107)

и может быть представлена в форме (5.103), если определить Мх,ы = м1,ы = {0, 1, 0).

Сами по себе массовые и разностные операторы, фигурирующие в уравнении (5.103), идентичны тем операторам, которые фигурировали в уравнении (5.93). Главным отличием первого из указанных уравнений от второго является то, что в результате интегрирования по двум измерениям возникает тензорное произведение одномерных операторов. Указанное отличие и обозначения массовых операторов были использованы в § 8.3. Более подробное обсуждение этих обозначений дается в приложении А.2. Важность введенного обозначения состоит в том, что оно позволяет осуществлять почленную конечно-элементную дискретизацию вполне аналогично тому, как это де--лается в методе конечных разностей.

Систему уравнений, получаемую путем записи уравнения (5.103) для каждого внутреннего узла (/ = 2, ..., / - 1; k = 2, /С-1), очень удобно решать с применением ПВР (см. § 6.3) подобно тому, как это делалось в п. 5.2.3. Чтобы обеспечить согласование с распечаткой программы DUCT (рис. 5.22), соответствующий алгоритм записывается в форме

= 0.75 [1. +(PAR 1/6) X

+ PAR2 {wl + + PARS .)]/PARI.

(5.108)

PARI =iblaf/Ax + l/Ay\ PAR2 = [- (b/af/Ax + 2/АуЦ/З, PAR3 = [2 {bjaflAx - l/АуЦ/З.

€ помощью соотношения (5.108) определяется предварительное решение, а формула (5.109) дает улучшенное решение

tk =Kk + к - к k) (5-109)

где X - параметр релаксации согласно ПВР.

Вышеприведенная формулировка реализуется в программе DUCT (рис. 5.22) с целью определения распределения скорости в канале прямоугольного сечения (рис. 5.21). Программа

3 С DUCT.FOR COMPUTES THE FULLY DEVELOPED VISCOUS LAMINAR FLOW IN A

4 с RECTANGULAR CHANNEL OF ASPECT RATIO, B/A ( BAR), BY FEM (LINEAR

5 С ELEMENTS) AND FDM (3PT CENTERED DIFFERENCES).

6 С SOR IS USED TO OBTAIN ITERATIVE SOLUTION

8 DIMENSION VEX(41,41) VA(41,41)

10 OPEN(1, FILE- DUCT.DAT)

11 OPEN(6,FILE DUCT.OUT)

12 R£AD(l,l)NX,ME,NEM,IPRrITMX,BAR,EPS,OM

13 1 F0RMAT(5I5,3E10.3)

14 С

15 IF(ME .EQ. 1)VRITE(6,2)

16 IF(ME .EQ. 2)VRITE(6,3)

17 2 FORMATC DUCT FLOW BY FEM: LINEAR ELEMENTS)

18 3 FORMATC DUCT FLOW BY FDM: 3PT CEN. DIFF.)

19 WRITE(6,4)NX,NEM,BAR,ITMX,EPS,0M

20 4 FORMATC NX M3/ NEM-M2/ B/A \F5.2/ MAX ITER M4

21 1 CONV TOL ЧЕЮ.З,* ОМ \TS.3J/)

22 С

23 PI 3.1415927

24 NH NX/2 + 1

25 NY NX

26 NXP NX - 1

27 NYP NXP

28 ANX NXP

29 DX 2./ANX

30 DXS DX*DX

31 DY DX

32 DYS DXS

33 BAS BAR*BAR/DXS

34 PARI BAS + i./DYS

35 PAR2 (2./DYS - BAS)/3.

36 PAR3 (2.*BAS - l./DYS)/3.

37 С

38 С SET UP EXACT AND INITIAL SOLUTIONS

39 С

40 DO 7 J 1,NX

41 AJ J - 1

42 X -1. + AJ*DX

43 DO 5 К 1,NY

44 AK К - 1

45 Y -1. + AK*DY

46 С

47 CALL VEL(X,Y,NEM,PI,BAR,VD,QLX}

48 С

49 WEX(J,K) WD

50 WA(J,K) WD

51 5 CONTINUE

52 IFdPR .EQ. 0)GOTO 7

53 WRITE(6,6)(WEX(J,K),K 1,NY)

54 6 FORMATC WE=M0F7.4)

55 7 CONTINUE

56 WRITE(6,8)

57 e FORMATC/)

Рис. 5.22. Распечатка программы DUCT (начало).