Р{1, г) =

iV-OO

(11.12)

где S (Z, m)/N - число персонажей S {I, т), достигших цели, деленное на общ,ее число персонажей N. Очевидно, проследить за бесконечным числом персонажей было бы неэкономично, а ве-фоятность можно рассчитать, имея дело с конечным числом персонажей N, порядок которого может быть от 10 до 10. При этом нужно оценить ошибку [i при аппроксимации бесконечности этим относительно малым числом.

При N > 20, используя центральную предельную теорему и соотношения, описывающ,ие нормальные распределения вероятностей, можно показать [3, 7], что приведенное ниже соотношение справедливо в тех случаях, когда рассматриваемые S персонажей выходят из некоторого исходного пункта и либо достигают конечного пункта с вероятностью Р, либо не достигают его с вероятностью 1 - р. Вероятность того, что среднее значение S (1, m)lN при конечном числе N отличается от {S {I, m)/iV]iv-co менее чем на некоторую величину \i, определяется выражением

Р

-(4-)

iV-CO

T1/V2

J e-i*dn* = erf-

n

(11.13)

где

(11.14)

Значения функции ошибок (erf) приводятся в стандартных таблицах [14, 15].

Во многих задачах такой анализ ошибок не может быть использован, поскольку персонажи начинают свой путь не из одного пункта. Например, поток излучения в точке, расположенной на границе некоторой замкнутой полосы, может зависеть от энергии излучения, испускаемого многими источниками. В таких случаях самый прямой путь оценки ошибки результата (например, ошибки в значении величины потока излучения в точке) состоит в расчете статистического среднего для группы из / средних значений. Затем используется центральная предельная теорема, согласно которой статистические пульсации средних имеют нормальное (гауссово) распределение относительно статистического среднего для групп средних значений (группового среднего). Эта величина называется вариацией. Например, если делается 200 проб, то групповое среднее Р рассчитывается на основе 200 проб и так же рассчитываются 20 средних значений Pj, Pj- . . . Pi для каждой из 10 проб. Тогда вариация среднего решения Р будет равна

2 (Pi-p)

(11.15)

Такая вариация является просто среднеквадратичным отклонением группового среднего значения Р от истинного; последнее можно было бы получить с помощью бесконечного числа проб. Из свойств нормального распределения ошибок, которому в общем случае подчиняются результаты расчетов, выполненных методом Монте-Карло, следует, согласно большинству учебников по статистике, что вероятность нахождения группового среднего* Р в пределах ± от истинного среднего равна /68%, в пределах ± 2 7 - 95 % и в пределах ± З7 - 99,7%.

Другой мерой статистических пульсаций среднего является 7 - стандартное отклонение. Так как у находится извлечением квадратного корня из (11.15), то, очевидно, чтобы уменьшить у вдвое, необходимо учетверить число проб, используемых при вычислении результатов (тем самым учетверяется число / при неизменном объеме Pi). Это, возможно, означает учетверение требуемого машинного времени, если не удается уменьшить член в квадратных скобках вследствие уменьшения вариации (рассеяния)

11.3.4. Оценка ошибки

Поскольку решения, получаемые методом Монте-Карло, являются средними по результатам ряда отдельных испытаний, то они будут в абш;ем случае содержать пульсации относительно средней величины. Как в любом процессе такого рода, большая точность нахождения среднего может быть достигнута путем увеличения числа величин, используемых при его определении. Хотя и не удается обеспечить 100%-ный доверительный уровень для получаемой величины, можно подойти к нему сколь угодно близко если это позволяет бюджет машинного времени. Обычно формулируются некоторые конкретные правила экономической целесообразности и оценки желаемой точности для данной задачи, которых затем придерживаются при ее решении.

Чтобы установить точность решений, применяют один из сле-дуюн],их тестов. Например, пусть нужно узнать вероятность достижения некоторого пункта в пределах города при случайном блуждании джентльмена, о котором шла речь в разд. 11.3.1. Чтобы точно определить его успех, мы должны были бы проследить за бесконечным числом условных персонажей и определить вероятность Р {I, т) достижения ими цели (Z, т) в следующ,ем виде:

Г S (I, т) -1

отдельных средних значений. Много времени и изобретательности затрачено на попытки уменьшить рассеяние с помощ,ью методов, имеющих названия: выборка по группам, русская рулетка и выборочные испытания. Эти и другие способы уменьшения вариации рассмотрены в работах [3, 7]. Ощутимая экономия машинного времени при использовании таких методов является щедрой компенсацией затраты времени на их изучение, и читатель, намеревающийся использовать метод Монте-Карло для какой-либо довольно сложной задачи, вынужден применять их.

11.4. ПРИЛОЖЕНИЯ К ПЕРЕНОСУ ТЕПЛОВОГО ИЗЛУЧЕНИЯ

11.4.1. Введение

Как указывалось в гл. 8-10, при составлении баланса энергии в условиях теплообмена излучением получаются интегральные уравнения, в которых неизвестными являются распределения температуры по поверхности или теплового потока. Интегральные уравнения получаются также при рассмотрении теплообмена излучением в излучающих средах типа газов. Эти уравнения достаточно трудно решаются]и являются следствием использования макроскопического подхода при выводе соотношений для тепловых потоков. Принимая вероятностную модель процесса радиационного обмена и метод Монте-Карло, можно использовать иолумакро-скопический ) подход и избежать многих трудностей, присущих процессам осреднения нри записи обычных интегральных уравнений. При таком подходе можно исследовать действие лгалых составляющих интегральной энергии, а не пытаться решить сразу задачу для всей энергии. Рассмотрим микроскопическую модель процесса радиационного обмена, а затем приведем два примера.

11.4.2. Модель процесса радиационного обмена

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

- Поверхность пдощадью Аг, т, = О

Типичный путь ргш

dQe,i = ei{Ti) oT\dAu

где 6i {Ti) - интегральная полусферическая степень черноты, определяемая в этом случае соотношением (3.6а):

) Или иолумикроскопический .

I I ei, 1 Pi. i) *{b i) cos p dco dX 0 Q

(11.19)

фотонов, которая представляет собой пучок, переносящий данное количество энергии w. Мы можем легко представить себе эту модель в виде группы фотонов, связанных вместе. В тех задачах, в которых рассматриваются спектральные величины и длина волны пучка задана, фотоны с этой длиной волны группируются вместе, образуя пучок с энергией w.

При одинаковой энергии всех пучков локальный поток энергии рассчитывается путем подсчета числа пучков, достигающих рассматриваемой единицы площади в единицу времени, и умножения результата на энергию пучка. Пробеги пучков и их истории рассчитываются методом Монте-Карло, как будет показано в одном из примеров,

11.4.3. Пример решения

Рассмотрим, например, довольно простую задачу, приведенную ф^. ц 4 теплообмен излуче-в работе llbj, и исследуем пере- нием между двумя поверхностя-нос энергии между элементарной ми.

площадкой dA, находящейся при

температуре Т^, и бесконечной плоской поверхностью А^, нри температуре = О (фиг. 11.4). Пусть степень черноты элементарной площадки

€1,1 = €{.1(Я, Pi, ГО, (11.16)

а поверхности 2

6i,2 = €l2(>, р2, г^а). (11.17)

Допустим, что степень черноты обеих поверхностей не зависит от азимутального угла 6.

Для элементарной площадки dA интегральная энергия излучения, испускаемого в единицу времени, равна

(11.18)

а Ч.ъ, 1 Ti) - спектральное распределение Планка интенсивности излучения черного тела ири температуре Т^.

Если предположить, что dQe,i - интегральная энергия излучения, испускаемого в единицу времени элементарной площадкой dAl,- состоит из N пучков энергии, то энергия каждого пучка w будет равна

dQe,l

JV = -

(11.20)

Чтобы определить энергию, передаваемую от элементарной площадки dAl поверхности А^, проследим за N пучками энергии, испускаемыми dAi, и определим число пучков 5j, поглощаемых А^- Если энергия излучения, отраженного от А^ к At, а затем обратно к Л2, пренебрежимо мала, то ноток энергии от dAi к А^ будет равен

1 поглощенная 2

€1(1) ndAj

(11.21)

Теперь рассмотрим, как определить направление распространения нучка и его длину волны. Это необходимо сделать таким образом, чтобы направления и длины волн N пучков соответствовали ограничениям, определяемым степенью черноты новерхности и законами, описывающими радиационные процессы. Например, если мы задаем длины волн N пучков, то спектральное распределение энергии испускаемого излученпя, задаваемое процессом Монте-Карло (распределение энергии wNi АХ для дискретных интервалов АА,), должно достаточно точно воспроизводить действительный спектр энергии испускаемого излучения (в виде зависимости nx,ii№,idX от X). Чтобы гарантировать это, применяют методы, изложенные в разд. 11.3.2.

Поток излучения, испускаемого элементарной площадкой dAi в интервале длин волн dX, включающего длину волны А в пределах телесного угла di, осью которого является нанравление р^, равен

dQie. 1 {К Pi) = 2ла, 1 (Х, Pi, Tl) ilb, 1 {X, Tl) cos Pi dAl sin pi dp, dX.

(11.22)

Интегральный ноток излучения, испускаемого dAi, определяется уравнением (11.18). Вероятность Р {X, р^) dp dA, исиускания излучения в интервале длин волн, включающем X, в телесном угле, осью которого является нанравление р^, равна тогда потоку излучения в интервале didX [уравнение (11.22)], деленному на

интегральный ноток пзлучения

die, 1 {X, Pi)

уравнение (11.18)]:

Р {X, pi) dPi dX =

dQe,i

nex, 1 {X, Pi) 1 W COS Pi sin Pi d?,i dk 6iorf

(11.23)

(Для простоты указание на функциональную зависимость от Ti опущено.)

Предположим теперь, что нанравленную спектральную степень черноты можно представить в виде произведения функции двух неременных: угла и длины волны, т. е.

€я,1(А, Pi) = ®i(l)©2(Pi). (11.24)

Это допущение, вероятно, несправедливо для большинства реальных поверхностей, так как в общем случае угловое распределение степени черноты зависит от длины волны, как показано, например, на фиг. 5.1. Из (11.24) следует, что зависимость стенени черноты от каждой неременной может быть найдена интегрированием по другой переменной [см. уравнение (11.9)]. Тогда нормализованная вероятность исиускания излучения в интервале dX будет равна

P{X)dX = dX J Р{Х, pi)dPi = о

2ndX 1 €х i(A., Pi) xb i(X.)smPiCOspidPi

=---- ( 2>

Подставляя результат в (11.4) и учитывая, что Р (Х) dX = О в интервале - оо <;А <iO, получим

я, л/2

2 I 1 Pi)4b 1 (>b*)sinp,cosPidpidX*

R.=---. ( -256)

где звездочкой отмечена неременная интегрирования. Если число пучков N велико и это уравнение решается относительно X для каждой случайной величины R, то для практических приложений время счета оказалось бы слишком большим. Чтобы обойти эту трудность, уравнения тина (11.256) можно заранее численно проинтегрировать для ряда значений А и но полученным результатам построить кривую. Часто хорошим ее приближением оказывается полином

X = A + BRf,-CRl+ ... . (11.26)

Полином (11.26) чаще используется в программе решения задачи, чем уравнение (11.256).

Аналогичным образом получим соотношение для полярного угла Pi

Pi 00

йр,= j J/>(pf, ;i)ddp! =

о о

2я U €х, 1 Pi) хь, 1 W sin Pf cos p* dk dpf

-, (11.27)