Марчук Г.И., Лебедев В.И. Численные методы в теории переноса нейтронов

Подождите немного. Документ загружается.

Учитывая (9.8.5)

и то, что J

является одновременно четной функ-

цией относительно г|)

= 0 и

-ф

= л/2,

приблизим

формулой

(Л

Ф)

« ^о W + а

(г)

cos

2 ф,

(9.8.6)

где коэффициенты

выберем

так,

чтобы (9.8.6) превращалось

в тождество при

i|) = 0,

л/2. После такого выбора, подставляя (9.8.6)

в (9.8.5), получаем квадратуру

Фо(П) =7(r

0)/4

/(0,

г,)/2

(9.8.7)

при

г« < у

. Для t/i < г < г/

выведем более точную квадратуру.

Для этого функцию

представим

виде

h (*\

*Ф)

« я

(/-) +

«1

cos 2 ф

+ а

(г)

cos 4 ф. (9.8.8)

В формуле (9.8.8) коэффициенты щ

t = 0, 1,2,

выберем так, чтобы

(9.8.8) превращалось

тождество

узлах построенной сетки,

т. е.

при -ф

= 0,

*ф

= л/2,

-ф

arctg {yjx^j,

где х

УV

—

у\.

Тогда, под-

ставив (9.8.8)

(9.8.5), получим квадратуру

+ {Z-ylM)J(0

r)].

(9.8.9)

Недостатком квадратуры (9.8.9) является

то, что ее

веса могут

менять знак. Ниже приведем формулу

положительными весами.

Теперь получим квадратуру для

< г ^ R. Для

этого обозначим

*' = *i = Vr

—

yl; у' = уъ х

= Уг

—

у\\

у"

= t/

, где k =

{max

j:yj<r}.

Сначала преобразуем функцию ф

(г):

_ г

Фо(г)

= —

Г (С08«ф + 81п«ф)У^

-1-

7 (}/7

у) X

J яг

X Vr

=tJ

dy+§7(х, Vr*^*) Vr^^dx 1

или

Фо

(г) = (Л + /, + /

)/(2 яг

где

j 7 (1/^=7,

(/)

VF=y*~dy+2

(*,

У7*=]?) ]/7^=^dx;

A-

f j{Vr

-y\y)—^='

==2r

(*,

1/T

—^=r

•

220

Пусть полуокружность г

= х

+ у

, у > 0 пересекает харак-

теристики y =

yj9

]=l,

2, ..., k, в точках Pj (x

yj), Р {—x

yj),

где Xj ^ 0. Эти точки бозьмем за узлы квадратур для /

, /

Оба интеграла в 1

заменим формулой трапеций по точкам Pj.

Интегралы /

, /

заменим формулой Симпсона, используя при этом

четность по у функции J и значения У, J соответственно в точках

(Л 0), (*i, f/i) и (0, г), (±*

В результате получим

(r) = S(/)/(2nr

); (9.8.10)

5(^

= -|-г

7(г,0) + 2 B

J{x„y

) + ±rx

J{0

(9.8.11)

где

Bj = Xj(yj+i—yj-i) + yj(Xj-i—Xj+i)>

/=2,3,...,

—

Недостаток квадратуры (9.8.10) в том, что она не является точной

при J = const и не учитывает того, что функция J (х, у) слабо меняется

not]) в окрестности г = 0. Для повышения точности введем при t/

^ )? новые квадратуры вида

<PeM = S(/)/S(l)

(9.8.12)

избавленные от упомянутых недостатков; они являются вполне удов-

летворительными по порядку точности на классе уравнений с разрыв-

ными коэффициентами.

С этой же целью можно модифицировать и квадратуру (9.8.9):

если функция / (х, у) гладкая, она мало меняется в окрестности г = 0;

поэтому можно для Уг<.г^у

взять квадратуру (9.8.10), (9.8.11),

положив в ней к = 1 и

= 2 г [яг/2 — (2/3) (х

+ yj)]. (9.8.13)

Тогда она будет точной при J = const и будет содержать только по-

ложительные веса.

Точность квадратур (9.8.10)—(9.8.12) повышается для гладких

функций, если точки P

Pj расположены равномерно на полуокруж-

ности х

+ у

= г

. Поэтому если характеристики у = yj в сетке

провести по закону

yj = R sin [/я/ (2ЛГ

+ 1)], / = 0, 1, ..., N

то они разделят полуокружность х

+ t/

= -R

; у ^ 0 на равные

дуги. Тогда для г

близких к г = /?, будет получена хорошая точ-

ность квадратур (9.8.12). Для г = # можно воспользоваться квадра-

турой Чебышева повышенной точности

фо(/?)=

1^тг

7(i?

]

8Л4)

(xf + у* = R; (XJ, yj) £ D

), являющейся точной, если функция

J (R cos г|), R sin tf) представима в виде тригонометрического по i|? мно-

221

гочлена степени

+ 1 ПОП. В

формуле (9.8.14) исключена

из

узлов

квадратуры точка

(О, R). Это

выгодно сделать ввиду упрощения рас-

четных формул

для

интегралов, взятых

при г = R по

отрезкам

[О

<Ч> <

я/2],

[я/2 <г|> < я] [см.

(9.8.21), (9.8.24)].

Здесь следует заметить,

что

если коэффициенты уравнения

источ-

ник терпят разрыв

некоторой точке

г — р, то

функция

J (х, у) от

решения

ср (х, у, у)

имеет,

как

правило, разрыв

производной

по у

при

у = р,

распространяющийся слева направо вдоль характеристики

= р при х > 0.

Поэтому

ряде многозонных задач целесообразно

для повышения точности квадратур выбирать такое разбиение

{yj},

которое включало

бы все

точки разрывов коэффициентов уравнения

и источника.

Учитывая,

что

решение

ср

(х, у, у)

является четной

по у

функцией,

интегралы

по у,

входящие

интегральные члены, аппроксимируем

квадратурной формулой наивысшей алгебраической степени точности

для четных функций

[90]:

»>

(*,

У)

= f

<{>dy

- 2

(*.

У*

Уп)

RN>.

(9.8.15)

*-i

Узлы

веса формулы (9.8.15) совпадают

положительными узлами

весами квадратурной формулы Гаусса

[101] для

отрезка

[—1,1]

порядка

2 N

Оценка погрешности формулы (9.8.15) содержится

в ра-

боте

[90]:

если

ср

как

функция

принадлежит классу

C4N3,

(Ю,

1]), то

MiN

2iN

Г ^ Y*-

MiN

' 2-'*

(9.8.16)

' (W.+

l)!

L(2iV,+ l)...(4iV,)J

(4tf,)t

д^(р

—~

, или

|/?дг,|< E

N,-\

(ф, — 1, 1).

Здесь

(ф,

Оу

—

1,1) —

наилучшее приближение

на

отрезке

[—1,1]

функции

ср

многочленом

&-й

степени

по у.

Пусть

для

функции

ф (x

г/, v)

определены интегралы

я/2 1

5±Ф==—Г <W>fq>(±/?cost|>, Rsinty,

y)dy;

(9.8.17)

Я/2

У±Ф

= — f

cos^dt (4(±/?cos^, #sin\|>,

Y)(1—

)

1/2

dy;

(9.8.18)

Я/2

Тф

= Г

(R cos

я|),

R sin

г|),

7) cos

\|?di|).

(9.8Л9)

Легко проверить,

что

Фо (R)

S+ф

S-ф;

(Я)

У+Ф

У^Ф;

l = 1/2; У

1 = —

У^1

= 1/4; Т\ = 1.

(9.8.20)

222

Значения этих интегралов нам потребуются в итерационном ме-

тоде.

Получим квадратурные формулы для интегралов (9.8.17)—

(9.8.19),

узлами которых являлись бы точки у = уи\ tyj =

jn/(2N

+ 1); k = 1, 2, *.., N

; j = 0, 1, ..., N

. От квадратурных формул

потребуем, чтобы они имели достаточно высокую степень точности

случае,

когда функция

ср

(R costp, R

sintf,

у)

является алгебраическим

многочленом по у и тригонометрическим многочленом по <J). Кроме

того,

потребуем, чтобы точно выполнялись соотношения (9.8.20).

Учитывая квадратурную формулу (9.8.15), получаем

S±

= S

(±

Vi> У&

(

где

aj = (2N

+ l)"

при / ф 0; 2"

(2N

+ I)"

при / = 0; (9.8.22)

|i,

= xj R'

cos

*, = (R

- y})W. (9.8.23)

Формула (9.8.21) становится точной [101], если функции

J (±R cos ф, R sin

*ф)

представимы при 0 <

< я/2 четными три-

гонометрическими многочленами степени 2^2 + 1; соотношение

(9.8.20) для S

<p, S-Ф выполнено.

Интеграл Тер заменим формулой

Тф

Ту (ф)

2 «i ^

Ф (Л cos

fy,

sin

v), (9.8.24)

/=о

где Pi=l 2

il*i) *

Формула (9.8.24) без множителя р

дает хорошую точность для

интегрирования функций вида ц,ф, ибо в ней взяты чебышевские узлы

и веса. Множитель же p

который не изменяет общего порядка точ-

ности формулы (9.8.24) на классе достаточно гладких функций, нужен

для выполнения соотношения (9.8.20) для 7ф.

Аналогично, требуя для У

выполнения (9.8.20), полагаем

+Ф

2 ЛПФЫ)(1-Т1)

,/2

, (9.8.25)

где

±(%A

{l-yty'*Y

Изложим еще один способ получения квадратур для ф

[162],

для этого преобразуем интеграл для ф

, сделав в нем замену ф =

= arccos \i:

Фо

= ^-Г^(Ф,г(1~

[

)

)(1~|х

1/2

^. (9.8.26)

223

Пронумеруем все узлы сетки D

, лежащие на радиусе

г при

х > О,

в порядке возрастания координаты х.

Пусть этим

узлам соответствуют

значения [х

= 0 <

< ... <

\in

< fi«

= 1, а значению \ij

соответствует точка (xj, yj). Интеграл (9.8.26) заменим суммой ин-

тегралов, взятых по отрезкам

[\ij„

[xj], / = 1, 2, .., n

Предполагая, что J (ф, rY\— (г

) линейно зависит от \i на

1И7-1> Pj)h

интегрируя по отрезку [(A/-i» М» получаем

»*/-

7-i

где

Cj J

yj)

+ BjJ (xj_

yj-г),

= (|i

-|i

)-4(l-|A/-i)

-(l-|*})

,/2

—p;.!(arcsin [iy—arcsin

fx/.j)];

«i= —(i*i—^-^-ЧО—i*/-i)

1/2

—(i—i*9)

l/2

—

—(iy (arcsin \ij—arcsin

(i/-i)].

Объединяя эти квадратурные формулы

одну, получаем квадрату-

ру для ф

(г).

§ 9.9. ПОСТРОЕНИЕ СЕТКИ И РАЗНОСТНАЯ АППРОКСИМАЦИЯ

ПОВЫШЕННОЙ ТОЧНОСТИ КИНЕТИЧЕСКОГО УРАВНЕНИЯ

для плоского слоя

Пусть в слое 0<д;<# задано кинетическое уравнение

~1х~

в котором

= 2

(Фо

З^ФО

+ / (*, р);

[

i-^.+ 2

(9.9.1)

ах

1 1

Фо=у £ф(*'1*')Ф';

Ф1

= 4" J И-'ФС^ ^')Ф% (9.9.2)

а коэффициенты 2, 2

, [г и источник / (#, р) —кусочно-постоянные

функции от*. Для уравнения (9.9.1) поставим краевую задачу (5.2.11).

На [О, Я] построим сетку D

}

= 0 < х

< ... < x

_i <

< *

... <

= Я, t = 0, 1, 2, ..., N, в узлы которой включим все

точки разрыва функций 2, 2

, \х, f.

Пусть в зоне U

_i, x

] коэффициенты и источник уравнения по-

стоянны, тогда внутри этой зоны целесообразно сетку из п

точек

взять сгущающуюся к краям по закону

km ^

т-1~\~

(Хт—

m-l)

sin

~Г}

ТТ\ •

*,..., П

m~r i)

чтобы учесть возможные быстрые изменения решения в окрестности

границ зон,

224

На отрезке [—1 ^

\х

<! 1] построим сетку D

, образованную точ-

ками ±\i

, k = 1, 2, ..., М, в которой 0< \i

< 1 и \i

<: [i

k+1

Пусть D

= D

X D

означает декартово произведение сеток D

Рассмотрим любую из характеристик кинетического уравнения,

косинус угла которой с осью х равен \i

. Вдоль этой характеристики

уравнение (9.9.1) имеет вид

-^ + 2

<3,

(9.9.3)

где dyldl = \i

dtp

(х, \i

)/dx; dx = |г

d£.

Рассмотрим тот отрезок характеристики, который заключен меж-

ду плоскостями х = #!_!, х = #*, /=1,2,... Тогда согласно (9.6.13)

(*ii I

*) =

<*о

(^-i, li

) +

Ро

Qi-i/2,

(9.9.4)

где а

, р

, Q/_i/2 выражаются формулами (9.6.14)—(9.6.16), в ко-

торых

/ ** _

Qi-i/2

+ 3jTfi

J Ф1^)+/(^_

1/2

);

h = xi —х^

; т = 2Л/(2щ); *!+i/a = (*i + ^i-i)/2.

Аналогичные формулы получаются для характеристик с косину-

сом угла наклона к оси х, равным —\i

; в них ф (х

—\i

) выра-

жается через ф (хи —

ц*)-

Перед выводом квадратурных формул для Q заметим, что функции

и фх связаны между собой уравнением баланса

-^ +

So Фо

=/„(*). (9.9.5)

Поэтому

если мы применим квадратуру Эйлера к интегралу от ф

учтем

(9.9.5),

то получим, что при ц

£ С

, а ( U*-i, **])

J ф1^ = у(фи +

Фп-1)

+ -^-2о(ф

ог

-фое.1) + 0(А

+Э), (9.9.6)

i-i

где

p = min (4, n + а). В (9.9.6) мы учли, что функция / (х, \ik) по-

стоянна

внутри U|_

Xi].

Для

получения квадратур от ф

предположим, что функция ф

(х)

между

точками хи *i-\ продолжена так, что она является решением

Рх-уравнений

метода сферических гармоник на lx^

xi\ — уравне-

ния

(9.9.5) и

-LJ2L

lVl

= f

(9.9.7)

3 dx

принимающим

заданные значения на концах отрезка

[**_!,

] :

Фо

(xt) =

;

(*i-j) =

|-1#

225

Если проинтегрируем (9.9.5) по отрезку Uj-i, x

]

а затем вос-

пользуемся уравнением (9.9.7) для исключения функции ф

и учтем,

что /

(х)

постоянна внутри [х

], то получим

Функцию ф

(х)

будем приближенно искать в виде

Фо

(*) =

Л"

(фо!-1

(**

—

Фо|

(x—Xi-г)

Ь (х—Xt-г)

(х

—х)),

(9.9.9)

где постоянную

выберем такой, чтобы для функции

(*), представ-

ленной формулой (9.9.9), было выполнено равенство (9.9.8). Подстав-

ляя (9.9.9) в (9.9.8), находим

* = (! +

) —^"(Л

^"(Фог

+ Фог-!)]-

При этом

J ФоЛс- ^(1+^^^)"

(ФО

-1 + ФО£+-^-/О). (9.9.10)

Подставляя (9.9.6), (9.9.10) в выражение для Q/_i/2, получаем

"**'** (ф

|-Фо^1) + ^^(Ф1е +

Фп-1)

+ /(^--1/2,^);

Qi-i/i =

от

)

1/2

+ т

-Qi^/3.

Краевые условия (5.2.11) аппроксимируются естественным образом:

в них следует считать fx = ±\ik-

Теперь изложим метод получения разностного уравнения, пред-

ложенный В. И. Лебедевым и С. А. Фроловой, который точнее учи-

тывает экспоненциальный рост или убывание решения вдоль характе-

ристик. Для этого воспользуемся формулой (9.6.2), которая для рас-

сматриваемого случая примет вид

ф(*ь

Р*)

= ехр(—2т

)ф

(*!_!,

) +

+ — Г Q(^*)exp((*-**)— )dt

(9.9.11)

i-l

где T

= 2ft/ (2|ii

). Подынтегральный член в (9.9.11) вычислим в

предположении, что функции ф

, ф

входящие в него, связаны между

собой не только уравнением баланса (9.9.7), но и продолжены внутрь

отрезка [х

] по значениям на концах его как точные решения

226

- или Я

-УРавнений метода сферических гармоник. Это значит, что

к уравнению (9.9.5) мы добавляем второе уравнение

i_^L

l9l

= /;

, (9.9.12)

в котором JSj = 2

; 7i = /1 для

приближения и 2

= S

/(1 +

+ 4/2

/52); /

= /

/ (1 + 42

/52) для Р

-приближения. Рассмотрим

два решения системы (9.9.5), (9.9.12): первое — такое, в котором

функция ф

(х) удовлетворяет условиям

Фо

(**-i)

Ф01-1;

Фо

(**) =

ф<н,

(9.9.13)

второе

—

такое, в котором функция ф

(х) удовлетворяет условиям

Ф1

(*t-i)

= Фп-i;

Ф1

(*«) = фц. (9.9.14)

Заметим, что можно было бы применить для продолжения функ-

ций ф

, ф

Pg-приближение. Тогда можно потребовать от построенных

таким способом продолжений одновременного выполнения краевых

условий (9.9.13), (9.9.14). Однако это усложнило бы расчетные фор-

мулы.

Итак, решение задач (9.9.5), (9.9.12), (9.9.13) и (9.9.5), (9.9.12),

(9.9.14) запишется в виде явных формул, которые можно рассматри-

вать как интерполяционные формулы, учитывающие специфику ре-

шений:

где

Фо (х)

(я—Xi-i) +

X {х

—x

) +

;

Ф1 (х)

Сп

(x—Xt^i) + c

i sh к (x—Xi) +fi/2

Cij

P(4>n{x)—q)'>

= —

р(ф

-1—g)\

(9.9.15)

(9.9.16)

A,= K32oV. p=

1/shX/i;

/o/S

при / = 0 и <7=7i/2i при /=1.

Подставляя выражения (9.9.15), (9.9.16) во второе слагаемое фор-

мулы (9.9.11) (обозначим его Q) и проводя точное интегрирование,

получаем, что

= (Л/|х*) {S

[/i(Cio+ 3|xiJt

) + y

+ 3|x|i

) +

+ /3 (/o/S

+ 3|i|i

)]

+ J

+ 3|i

/0),

где

-2ч

J^QJichkh — T

shXh— ХЛехр(—r

))/(X

—

т0;

= (ЯЛ—exp (—x

)

(Mi

ch tfi + x

sh Щ/(1

—%l);

—ехр (—т"

))/тд, т^ = SA/(i

(9.9.17)

(9.9.18)

Заменяя в У

экспоненту рациональной дробью (9.6.14) при а = ог

разлагая входящие в J

и У

выражения ехр (—r

), ехр (—ЯЛ), sh АЛ,

227

ch kh в ряды й пЬДсТавлйя эти выражения Ь равенство (9.9.17), в ко-

тором ехр

(—2х

)

заменена рациональной дробью (9.6.14), получаем

разностное уравнение, определяемое формулами (9.9.11), (9.9.17).

Для расчета <р

, Фа применим квадратурные формулы вида

Л=1

(9.9.19)

где Л

, \i

, k = 1, 2, ...,

УИ,

—

веса и узлы некоторой квадратурной

формулы для [0, 1].

§ 9.10. РАЗНОСТНЫЕ СХЕМЫ МЕТОДА ХАРАКТЕРИСТИК

ДЛЯ САМОСОПРЯЖЕННОГО УРАВНЕНИЯ (ПЛОСКИЙ СЛОЙ)

Рассмотрим уравнение (5.2.12). В нем, не уменьшая общности,

можно предположить, что / = const, так как общий случай одномер-

ной задачи приводится к задаче, в которой величина / постоянна. В

самом деле, вводя среднюю оптическую длину по формуле

x^^l^(t)dt

делая

соответствии с этой формулой замену переменной, получаем

кинетическое уравнение с постоянной функцией /. Итак, уравнение

запишем в виде

-'V-Tr +

cSu

+ F(x, |i)

(9.10.1)

ox*

где 0 < cjx) ^ 1, а оператор Su имеет при изотропном рассеянии вид

Su= f и(х, \i)d\i. (9.10.2)

Напомним, что и

(лг,

\i) и F (*,

\л) —

четные по

ц,

функции, и поэто-

му уравнение (9.10.1) достаточно рассмотреть лишь при 0 < \i < 1.

Пусть для простоты выкладок краевая задача для уравнения

(9.10.1) ставится в слое |л:| <: я. Построим на отрезке [—я ^ х ^ я]

сетку D

= {x

= ih, i = 0, ±1, ..., ±Л/}, где h = яЛ/"

, а N —

некоторое целое положительное число.

На отрезке [0 < ц < 1] построим сетку D

= {\i

k = 1,2, ...

..., М}, образованную точками \i

, k = 1, 2, ..., М, где М

—

неко-

торое положительное целое число. Пусть \х

ft+1

;

/ = М~\ а

= й

х D

—

декартово произведение сеток D

и D

. Введем в

рассмотрение два разностных оператора

и =

и(х—A,

ji)—2и(х,

\i)

+ u(x +

\i);

(ou=*(u(x—h,

ii)

u(x

\i))/\2

5u(x,

|x)/6. J

228

Интегралу

поставим

соответствие квадратурную формулу

5«а«

u(x

)

(9.10.4)

в которой

> 0, k = 1, 2, ...,

Ми2^/1=

Далее будет выясне-

но,

что за

узлы

выгодно брать положительные корни некоторых

ультрасферических многочленов степени

2М [119, 210].

Если

за эти

многочлены взять многочлены Лежандра,

то

квадратура (9.10.4)

бу-

дет квадратурой Гаусса, согласованной

выбором операции

Р в КР-

методе

(см. гл. 11).

Предполагая функцию

(х) непрерывной, построим

во

внутренних

узлах сетки

следующие конечно-разностные аналоги

для

уравне-

ния (9.10.1):

—/

/ryS

+ и =

+ F\

(9.10.5)

—Ph-y&u + ии =

со

(cS

u + F)

(9.10.6)

при \x

= \i

k = 1, 2, 3, ..., M. В

случае периодической задачи

эти

уравнения составляем

и для

точки

х = п.

Уравнения (9.10.5), (9.10.6) являются соответственно уравнения-

ми (9.6.23), (9.6.22), переписанными

для

плоского слоя.

На

классе

локальная аппроксимация

по х

уравнения (9.10.1) уравнением

(9.10.5) будет величиной

(ДО

где р

= min (2, п + а), а

уравне-

нием (9.10.6) — величиной 0(№),

где р

= min (4, п +

оь).

Для

урав-

нения (9.10.6) дополнительно потребуем выполнения неравенства

(9.6.18),

которое

данном случае будет иметь

вид

2/3^. (9.10.7)

§ 9.11. НАХОЖДЕНИЕ ОБЩЕГО РЕШЕНИЯ РАЗНОСТНЫХ

УРАВНЕНИЙ для плоского слоя

Сначала найдем общее решение однородных систем (9.10.5),

(9.10.6)

при с =

const.

Для

определенности найдем общее решение

системы (9.10.6), ибо для системы (9.10.5),

также для систем, аппрок-

симирующих несамосопряженное уравнение, соответствующие

вы-

кладки проводятся аналогично. Итак, однородная система (9.10.6),

записанная

для

точек

имеет

вид

—№~*\i4*ii

+ cor/ =

u. (9.11.1)

Решение

ее

ищем

виде

|i) = Ф

(v) ехр

(—ух/1), (9.11.2)

где

у —

неизвестный параметр. Пусть

yh/21,

а ф (Р)

обозначает

функцию

229