Книги по разным темам Pages:     | 1 | 2 | Журнал технической физики, 1997, том 67, № 9 01;07 Новый метод решения задач переноса излучения в излучающих, поглощающих и рассеивающих средах й В.С. Юферев, М.Г. Васильев, Л.Б. Проэкт Физико-технический институт им. А.Ф. Иоффе РАН, 194021 Санкт-Петербург, Россия (Поступило в Редакцию 22 мая 1996 г.) В основе предлагаемого метода лежит новый способ аппроксимации угловой зависимости интенсивности излучения. Вся область изменения телесных углов разбивается на N ячеек, симметричных относительно центра сферы. В каждой из ячеек излучение задается в виде P1 приближения, а для определения множества локальных нулевых и первых моментов получается система дифференциальных уравнений. В некоторых специальных случаях предлагаемый подход может рассматриваться как обобщение метода дискретных ординат, позволяющее естественным образом решить проблему выбора весов в квадр атурных формулах.

Эффективность метода демонстрируется на двух тестовых одномерных примерах. Показано, что в этом случае достаточно высокая точность решения задачи достигается уже при N = 2.

Введение условий в работах [7,8] была предложена его модификация, основанная на использовании самостоятельных разложений [9,10] в каждой из полусфер в одномерном Задачи радиационного теплообмена по-прежнему остаслучае и в каждой из четверти сферы Ч в двумерном.

ются одними из наиболее трудоемких с вычислительной точки зрения. Для их решения было предложено боль- Альтернативный подход, который можно назвать лошое количество методов (см., например, обзоры [1Ц3]), кальным, используется в методе дискретных ординат, кооднако ни один из них не может считаться достаточно гда все поле излучения разбивается на дискретное число универсальным и пригодным на все случаи жизни. Более потоков, с каждым из которых связано фиксированное того, несмотря на поразительное развитие вычислитель- направление в пространстве и соответствующее значеной техники, до сих пор, например, отсутствует решение ние весового коэффициента в квадратурных формулах, задачи радиационного теплообмена в такой достаточ- используемых для вычисления интегралов излучения. В но простой геометрии, как полупрозрачный круговой результате уравнение переноса излучения заменяется цилиндр конечной длины с зеркально отражающими конечной системой дифференциальных уравнений, опипрозрачными границами и коэффициентом преломления, сывающих пространственные изменения интенсивности большим единицы. излучения в этих направлениях. SN-метод был впервые предложен Чандрасекхаром [11] как обобщение метода Основные методы решения задач переноса излучения ШустераЦШварцшильда и в дальнейшем нашел широкое можно условно разделить на две группы в зависимости применение для решения различных задач радиационот того, какая форма уравнения переноса излучения ного теплообмена (см., например, [12Ц16]) вплоть до используется для получения численного решения Ч дифиспользования в коммерческих кодах [17].

ференциальная или интегральная. К первой группе относятся методы сферических гармоник, (PN) моментов, Известно, что точность метода дискретных ординат собственных функций Кейса, дискретных ординат (SN). определяется конструкцией квадратурной схемы, т. е.

Ко второй Ч зонный метод и метод конечных элементов выбором весовых коэффициентов. Какие-либо строгие и их разновидности. Отдельную группу образуют метод математические принципы, позволяющие находить знаМонте-Карло и Фray tracingФ методы. чения весов, в настоящее время отсутствуют, хотя недавно в этом направлении и были получены достаточно При решении уравнения переноса излучения в диффеинтересные результаты [18]. Указанное обстоятельство ренциальной форме основная проблема состоит в выборе существенно снижает эффективность данного метода. В способа аппроксимации зависимости интенсивности изнастоящей работе предлагается новый подход к решению лучения от направления. Обычно используется два подзадач переноса излучения, который позволяет естественхода. В первом угловая зависимость интенсивности изным образом решить эту проблему.

учения аппроксимируется набором функций, заданных во всем диапазоне телесных углов. Указанный подход В основе подхода лежит новый способ аппроксиможет быть назван методом глобальной аппроксимации. мации угловой зависимости интенсивности излучения, Практическое его применение ограничено случаями, объединяющий в некотором смысле метод сферических когда оказывается достаточным использование первых гармоник и метод дискретных ординат. Действительно, членов разложения, например P1- или P3-приближения SN-метод можно рассматривать как метод коллокации в методе сферических гармоник [4Ц6]. Для улучшения на сфере, поскольку уравнение переноса удовлетворяточности PN-метода при наличии разрывных граничных ется на конечном множестве фиксированных значений 1 2 В.С. Юферев, М.Г. Васильев, Л.Б. Проэкт углов, определяющих направления распространения из- С другой стороны, в отличие от P1-приближения коэфлучения. В противоположность этому в предлагаемом фициенты Am не совпадают с компонентами плотности i методе уравнение переноса излучения удовлетворяется потока излучения qm, а связаны с ними формулой (2а).

i в среднем в каждом из элементарных телесных углов Выражая из (2а) Am через qm иподставляяв(1), получим j n разбиения сферы. Угловая зависимость излучения в этих ячейках может быть аппроксимирована различными споAm = pmqm, (3а) i n n собами. В данной работе для этой цели используется Pn=приближение вследствие своей простоты и физической ясности. Таким образом, можно сказать, что предлагаm емый метод занимает относительно SN-метода такое же I = s-1(I0 + amqm), (3б) m n n положение, как метод конечных элементов относительно n=метода конечных разностей или метода коллокаций.

где am = lim pm, матрица {pm} является обратной к in n i=1 j Основная идея предлагаемого подхода была впервые матрице {pm}.

i j сформулирована в работе [19,20] как обобщение метода Представление (1)Ц(3) показывает, что полные момендифференциального приближения и показала прекрасты интенсивности излучения являются суммой локальные результаты при решении некоторых тестовых задач.

ных моментов, т. е.

В настоящей работе дается общее изложение данного метода.

N N m I0 = I0, qi = qm, i m=1 m=Формулировка метода N m Разобьем всю область изменения телесных углов на N I2,nj = I2,nj, i, j, n = 1, 2, 3, (4) областей (ячеек) m таким образом, что m =+ -, m m m=а подобласти +, - являются симметричными отноm m таким образом, все интегральные условия, которые в сительно центра сферы, т. е. для каждого направления методе дискретных ординат используются при конструиm,+ {li, i = 1, 2, 3} +, где lim есть направляющие m ровании квадратурных формул, в данном случае выполкосинусы относительно осей координат xi, существует няются автоматически.

m,симметричное направление {li, i = 1, 2, 3} -, m Разбиение сферы может быть выполнено различm,- m,+ так что li = -li. Представим интенсивность ными способами. В простейшем из них, используя излучения в каждой из ячеек m в виде, аналогичном сферическую систему координат, элементарные ячейки P1-приближению (для простоты изложения рассматриваm =m определяются множествами pq ется серое приближение) (p-1 < p) ( - p < - p-1), m m I(r, ) = s-1 I0 (r) + li Am(r), lim m, (1) m i (q-1

m водится в Приложении. При решении одномерных и некоторых многомерных задач разбиение (5) можно Легко показать, что благодаря свойству симметрии m упростить, положив 0 < 2, т. е. использовать углового разбиения коэффициент I0 оказывается равным только разбиение по углам. В этом простейшем локальному нулевому моменту интенсивности излучения случае в области m, в то время как для первого и второго локальных моментов будем иметь sm = 4(m-1 - m), = cos, m qm = ln I d= pm Am, (2а) n nj j pi j = 0, если i = j, и Am = qm. (6) i i m pii j=Другой способ конструирования {m} состоит в симm m m I2,nj = ln lmI d=pm I0, (2б) метричном разбиении сферы на области одинаковой j nj m формы и равной площади. Однако в этом случае разгде нообразие возможных разбиений оказывается ограниченным. Минимальное из них (N = 3) получается m pm = s-1 ln lmd, j, n = 1, 2, 3.

nj m j в результате проекции куба на сферу. Далее следует m N = 4 Ч проекция октаэдра на сферу, N = 6 Ч Можно видеть, что уравнение (2б) эквивалентно усло- проекция ромбического октаэдра и т. д. В общем случае вию замыкания в P1-методе применительно к ячейке m. разбиение {m} может зависеть от пространственных Журнал технической физики, 1997, том 67, № Новый метод решения задач переноса излучения в излучающих, поглощающих... координат и это является еще одним существенным Система дифференциальных уравнений первого порядка преимуществом данного метода, поскольку позволяет (8), (10) определяет множество локальных нулевых и m использовать конкретные особенности решаемой задачи. первых моментов I0 и qm. В принципе она может быть i Однако в данной работе мы ограничимся рассмотрением сведена к системе уравнений второго порядка относиm ситуации, когда разбиение {m} от пространственных тельно I0 или qm. Последняя система является более i координат не зависит. удобной для численного решения. Поскольку процедура приведения в общем случае является достаточно громоздкой, то мы проделаем эту операцию на нескольких Основные уравнения конкретных примерах.

а) Плоский слой, осевая симметрия. Испольm Чтобы получить уравнения относительно I0 и qm необi зуя разбиение (5), (6) и полагая, что ось 0X3 направлена ходимо вычислить локальные нулевой и первый моменты перпендикулярно границам слоя, будем иметь от исходного уравнения переноса излучения, которое имеет вид mn mn fimn = i3j3 f, hmn = i3l3 f /pn, qm = pm Am.

j il 33 33 I В результате уравнения (8), (10) принимают вид lj + I = xj j=1 dqm N n = (rmn - mn)I0 +(1-)smIB, (11а) dx3 n= (, )I(r, )d + (1 - )IB, (7) m dI0 N где = k +, = /(k + ), k и Ч есть pm = (hmn - mn)qn. (11б) dx3 n=1 коэффициенты поглощения и рассеяния соответственно.

Проинтегрируем (7) по области m. Тогда, подставляя Дифференцируя (11а) по x3 и исключая I0, получаем в интеграл рассеяния вместо I представление (1) и искомую систему уравнений второго порядка относиучитывая свойства симметрии индикатриссы рассеяния, тельно qm будем иметь d2qm N N qm N = 2 (rmn - mn) j n dx2 n=1 l=1 pn = (rmn - mn)I0 + (1 - )smIB, (8а) 3 xj n=j=dIB (hnl - nl)ql + (1 - )sm. (12) где dxrmn = (, )d d, (8б) 4sn n m Отметим, что уравнение (12) справедливо для любой функции рассеяния.

mn Ч символ Кронекера.

б) Многомерный случай, одномерное разУмножим далее (7) на li и проинтегрируем снова по б и е н и е. При некоторых условиях одномерное разбиеm. В результате с учетом (1), (2б) и свойств симметрии ние, использованное в предыдущем подразделе, может индикатриссы рассеяния получим оказаться достаточно точным и в многомерных задачах.

N m Поскольку в этом случае по-прежнему pi j = 0 при i = j, Ipm +qm = fimnAn, (9а) то вместо уравнения (12) будем иметь i j xj i n=1 j=1 j j j 2qm N N 3 j где = 2 (rmn - mn) xjxi n=1 l=1 k=1 pn ii j=fimn = limln(, )d d. (9б) j 4sn m n j dIB (hnl - nlik)ql + (1 - )sm. (13) ik k dxi Подстановка сюда выражения (3а), связывающего коэффициенты Am с плотностью потока излучения qm, дает j в) Многомерный случай, разбиение общего вида. Рассеяние изотропное. Из (8б), (9б) следуm I0 N 3 ет, что fimn = 0; rmn = sm/(4). В результате, исключая j pm = (hmn - mnil)qn, (10а) i j l qm из (8), (9), будем иметь xj n=1 l=1 il i j=3 N m где 2I0 sm n pm = -2 -mn I0 -2(1-)smIB.

i j xixj n=1 i, j=hmn = fimn pn. (10б) il j jl (14) j=1 Журнал технической физики, 1997, том 67, № 4 В.С. Юферев, М.Г. Васильев, Л.Б. Проэкт Сравнение предлагаемого подхода где Ч коэффициент черноты; s, d Ч коэффициенты зеркального и диффузного отражения: n Чвектор внус методом дискретных ординат тренней нормали к границе области, а направление связано с условием, что угол отражения равен углу Как уже отмечалось, в методе дискретных ординат уравнение (6) заменяется системой уравнений (для про- падения.

Проектируя уравнение (18) на внутреннюю нормаль и стоты изложения рассматривается случай изотропного интегрируя по области m, а более точно, по тем напрарассеяния) влениям, принадлежащим m, которые удовлетворяют 3 N условию n > 0, получим уравнения, связывающие Im m li + Im = nIn + (1 - )IB, m I0 и qm на границе области xi 4 i i=1 n=-N (n )I()d=IB(T) (n )d m = 1, 2,..., N, (15) m,n>0 m,n>где индекс m указывает соответствующее направление в пространстве, n Ч веса квадратурной формулы.

d N + (n )sI( )d+ Обычно направления выбираются таким образом, что m m n=li = -li-m. Тогда, обозначая J0 = I-m+Im, qm = Im-I-m, m,n>после стандартных преобразований уравнений (15), получим |n |I( )d (n )d. (19) m Jn,n <0 m,n>lim + qm = 0, xi i=Вычисление интегралов, входящих в уравнение (19), представляет наиболее трудоемкую часть при подготовqm m lim + Jке численного алгоритма, основанного на применении xi i=данного метода. Наиболее легко указанные вычисления проводить в одномерном случае для плоского слоя.

N n Используя разбиение (5), (6), получаем = 2 nJ0 +(1-)IB. (16) n=1+ s d m m m I0 (1 - s ) +2qm =IBsm - sm Исключая из (16) qm, будем иметь m-1 +m m N 2J0 m qn n n +n-limlm + Jj -I0, xixj 2 i, j=n=N где m-n = -22 nJ0 +(1-)IB. (17) s = m 2 2 sd, n=m-1 - m m m-Сравнение уравнений (14) и (17) показывает, что 1 2s() s = d.

m они имеют одинаковую структуру. Следовательно, предm-1 - m m pлагаемый подход можно рассматривать как обобщение метода дискретных ординат, позволяющее естественным Численные примеры образом решить проблему выбора весов в квадратурных формулах. Более того, в данном методе невозможно Чтобы продемонстрировать эффективность данного возникновение Фray effectФ, описанного в [16] примениметода были рассмотрены две задачи.

тельно к SN-методу.

Pages:     | 1 | 2 |    Книги по разным темам