(1-f Y)Ar?+-YAr,%

= а(Л1,® Lxx + ® Lyy) [(1 - р)П k + Ttkl (8.43)

т. е. соотношение, эквивалентное (8.26). Однако здесь массовые операторы (8.30) по координатным направлениям будут представлены в обобщенном виде

Мх = м1 = {б, 1 - 26, б). (8.44)

Ясно, что конечно-разностная схема (п. 8.2.2 и 8.2.3) соответствует значению 6 = 0, тогда как конечно-элементная схема (п. 8.3.1) соответствует значению б = 1/6.

После приближенной факторизации уравнения (8.43) получается следующий двухэтапный алгоритм, заменяющий (8.39) (8.40):

8.3.2, TWDIF: обобщенная схема для реализации расчетов с конечными разностями или конечными элементами

Различные схемы расщепления, рассмотренные в п. 8.2.2, 8.2.3 и 8.3.1, могут быть объединены в одну обобщенную трехслойную схему, охватывающую как конечно-разностную, так и конечно-элементную формулировки. Мы покажем это здесь для области, изображенной на рис. 8.1.

Распределение температуры в области Ojc 1, Оу 1 определяется уравнением (8.1) со специальными граничными условиями

Г (О, у, t) = 20 + 80y, f(l, у, 0=20 + 80[y-e- sin(0.5nf/)],

Г(:с, 0,0 = 20, f{x, 1, 0=20-f 80[(/-в-о-5 8т(0.5яд:)].

Для данного примера существует следующее аналитическое решение:

f = 20 -f 80 [у - е-- sin (0.5ях) sin (0.5я(/)]. (8.42)

Обобщенная схема может быть введена за счет трехслойной конечно-разностной дискретизации уравнения (8.29) при а = ах = oLy, что дает

Первый этап где

ms = [аЩМу®Ux + Мх0Lyy)Tl k + yMx<2>My 11 *]/(! + v).

(8.46)

Второй этап

\Му - ((Г1) lyy] А-П == АГ;. (8.47)

Для уравнений, формируемых в точках вблизи границ, нужно задать граничные значения величин АГ* в (8.45) и АТ - в (8.47). Надлежащий вид этих граничных значений проще всего определить путем свертывания (8.45) и (8.47) в одну комбинированную схему, получаемую за счет исключения АР. Тогда становится ясно, что требуемые для решения (8.47) граничные значения А7* могут быть получены непосредственно из формул (8.41), тогда как граничные значения для ДР должны вычисляться также из (8.41) через посредство обращенного соотношения (8.47); например, в точке j=NX (х=\) имеем

ДГкх. . = [Af, ~ [-f) М А7йх\. (8.48)

Если граничные условия (8.41) не зависят от времени, то для АР и Ar+i получаются нулевые граничные значения. Известные граничные значения зависимых переменных АР и АР+ переносятся в правые части уравнений (8.45) и (8.47) соответственно до решения уравнений.

Реализация вышеописанного алгоритма проводится с помощью программы TWDIF (рис. 8.4). На первом этапе соотношение (8.45) приводит к трехдиагональной системе уравнений, связанных с каждой из сеточных линий, идущих в направлении X. Эти уравнения решаются с помощью подпрограмм BANFAC и BANSOL (рис. 6.18 и 6.19). Вычисление по формуле (8.46) реализуется в подпрограмме REDIF (рис. 8.5). На втором этапе соотношение (8.47) служит средством получения трехдиагональной системы уравнений, связанных с каждой из сеточных линий, идущих в направлении у\ система решается с использованием подпрограмм BANFAC и BANSOL. Основные параметры, используемые в программе TWDIF, описаны в табл. 8.1. Типовая выдача результатов для случая конечно-разностной схемы Кранка - Николсона на грубой сетке показана на рис. 8.6.

2 С TWDIF APPLIES APPROXIMATION FACTORISATION TO SOLVE

3 с THE UNSTEADY HEAT CONDUCTION PROBLEM FOR T(X,Y) .

4 с REDIF EVALUATES THE RIGHT-HAND SIDE.

6 REAL*8 SUNT,RMST,DSQRT,AN

7 DIMENSION T(41,41),DT(41,41),R(41,41),EMX(3),EMY(3),

8 1B(5,41),RRT{41),DDT(41),SIX(41),SIY(41),Y(41),ERR(41,41>

9 COMMON DX,DY,EMX,EMY,NX,NY,R,T,DT

10 OPEN(1, FILE= TWDIF.DAT)

11 OPEN(6, FILE= TWDIF.OUT)

12 READ(1,1)IBC,NX,NY,ME,GAM,BET

13 READ(1,2)ALF,DTIM,DTMCH,TMAX

14 1 FORMAT(4I5,2F5.2)

15 2 FORMAT{4E10.3)

16 С

17 PI = 3.141592654

18 PIH = 0.5*PI

19 NXS NX + 1

20 NYS NY + 1

21 NXP = NX - 1

22 NYP = NY - 1

23 IFdBC .EQ. DNXN = NXP

24 IFdBC .EQ. DNYN NYP

25 IFdBC .EQ. 2)NXN = NX

26 IFdBC .EQ. 2)NYN = NY

27 NXPP = NXN - 1

28 NYPP NYN - 1

29 AN = NXPP*NYPP

30 ANX = NXP

31 - DX = l./ANX

32 ANY = NYP

33 DY = l./ANY

34 SX = ALF*DTIM/DX/DX

35 SY = ALF*DTIM/DY/DY

36 EMXd) = 0.

37 IF(ME .EQ. 2)EMXd) = 1./6.

38 IF(ME .EQ. 3)EMX(1) = 1./12.

39 EMX(2) = 1. - 2.*EMX(1)

40 EMX(3) = EMX(l)

41 DO 3 J = 1,3

42 3 EMY{J) = EMX(J)

43 WRITE(6,4)NX,NY,ME,GAM,BET,SX,SY

44 4 FORMATC UNSTEADY HEAT CONDUCTION WITH NX,NY =\2I3,/,

45 ! ME =M3, GAM =\F5.2. BETA =,F5.2, SX,SY =\2F6.3)

46 IF(ME .EQ. 1)WRITE(6,5)EMX

47 5 FORMATC APPROX. FACT., 3PT FDM, EMX=,3E10.3)

48 IF{ME .EQ. 2)WRITE{6,6)EMX

49 6 FORMATC APPROX. FACT., LINEAR FEM, EMX=,3E10.3)

50 IFdBC .EQ. 1)WRITE(6,7)

51 IFdBC .EQ. 2)WRITE(6,8)

52 7 FORMATC DIRICHLET B.C.)

53 8 FORMATC DIRICHLET AND NEUMANN B.C.)

54 VRITE(6,9)DTIM,DTMCH,TMAX

55 9 FORMATC DTIM =,E10.3, PRINT INT. =,E10.3, TMAX =,E10.3, )

56 С

57 С GENERATE INITIAL SOLUTION

58 С

59 DO 11 J = 1,NX

60 AJ * J - 1

61 X = AJ*DX

62 SIX(J) = SIN(PIH*X)