УДК 519.63:533.9.07
Д-р техн. наук С. И. Планковский, Е. В. Цегельник, В. О. Гарин
Национальный аэрокосмический университет им. Н. Е. Жуковского «ХАИ»
МОДЕЛИРОВАНИЕ ТЕЧЕНИЙ В КАТОДНЫХ УЗЛАХ ПЛАЗМЕННОГО ОБОРУДОВАНИЯ ПРИ ОПИСАНИИ ДУГИ ОБЪЕМНЫМ ИСТОЧНИКОМ ЭНЕРГИИ
Рассмотрены вопросы численного моделирования состава атмосферы в катодных узлах электродугового плазменного оборудования с использованием современных коммерческих С¥Б пакетов. Для задач проектировочных расчетов катодных узлов электродугового плазменного оборудования обосновано применение моделей с заданием дуги в виде объемного источника тепла. Проведено сравнение эффективности численной реализации при использовании двух вариантов построения расчетной модели: с выделением подобластей в расчетной области и при задании интенсивности источника при помощи метода Я-функций. На примере тестовой задачи показано, что алгоритм с использованием сетки конечных элементов без деления на подобласти обладает большей устойчивостью и эффективностью с точки зрения времени получения решения. Проведено сравнение результатов расчета состава атмосферы в катодном узле с полым катодом с учетом и без учета влияния дуги на характер течения.
Катодныш узел, математическое моделирование, объемный источник, метод Л-функций
Введение
Многие авторы работ, посвященных разработке и исследованию плазменного оборудования, отмечают, что изучение течений холодного газа при продуве им каналов плазмотронов является важ-ньш предварительным этапом в изучения характеристик и параметров плазменных установок [1]. Оно позволяет выявить структуру течения, формирующегося в канале в зависимости от конструкции и свойств его поверхности, конструкции электродов, способа ввода газа в рабочую зону и других параметров.
Традиционными целями моделирования холодной продувки плазменного оборудования являются предварительная оценка характера воздействия потока газа на дуговой разряд. Таким образом, можно составить некоторое представление о свойствах турбулентности в плазменной установке на основе изучения вихревой структуры холодного потока.
В работах [2—4] математические модели течения холодного многокомпонентного газа предложено использовать для определения режимов подачи защитного газа в плазменных генераторах с термокатодами. Такой подход является оправданным для стартовых режимов работы плазменного оборудования. Однако на рабочих режимах наличие интенсивного источника тепла (электрической дуги) существенно влияет как на картину течения, так и на свойства компонент газовой смеси. Поэтому для обоснованного выбора расхода и способа подачи защитного газа на рабочих режимах плазменного оборудования ак-
© С. И. Планковский, Е. В. Цегельник, В. О. Гарин, 2010
туальным остается вопрос разработки математических моделей, которые учитывали бы основные особенности течения газов в плазменном оборудовании — неравномерность температур, нелинейную зависимость транспортных коэффициентов от температуры и др.
В то же время, с учетом специфики проектировочных расчетов при разработке таких математических моделей необходимо тщательно учитывать требования к вычислительным ресурсам.
С учетом этого целью работы является разработка эффективных, с вычислительной точки зрения, моделей для расчета состава атмосферы в катодных узлах плазменного оборудования.
Анализ известнык работ в области моделирования течений в плазменном оборудовании
В настоящее время сформировалось два основных подхода к моделированию течений в электродуговых плазменных генераторах. Первый основан на представлении дуги в виде объемного источника энергии. При таком подходе учитывается зависимость свойств плазмообразующего газа от температуры, а интенсивность объемного источника определяется с учетом джоулева нагрева и излучения дуги [5, 6]. Полученные решения достаточно хорошо описывают течения на установившихся участках (например, в задачах о дугах в длинном цилиндрическом канале). Однако применение такого подхода для изучения характеристик течения вблизи электродов неизбежно приведет к большим погрешностям. Кроме того, в известных работах, использующих
представление дуги в виде объемного источника тепла, как правило, рассматриваются случаи ламинарного течения в канале [5—7] и не учитывается вихревой характер подачи газа, наиболее распространенный в промышленном плазменном оборудовании.
В последнее время все большее распространение получает другой подход к моделированию течений в плазменных генераторах, основанный на применении МГД приближения. Популярность такого подхода во многом вызвана введением возможностей моделирования МГД течений в последние версии коммерческих CFD пакетов (ANSYS CFX, FLUENT и др.).
Результаты, полученные с использованием трехмерных МГД моделей, в ряде случаев показали хорошее качественное и количественное совпадение с экспериментальными данными. Так, например, в работе [8] в ходе изучения течения в электродуговой плазменной горелке без введения каких-либо искусственных приемов смоделирован процесс привязки дуги к стенке цилиндрического анода при полностью осесимметрич-ной геометрии и граничных условиях. В работе [9] с применением данного подхода смоделирован процесс шунтирования дуги. В рамках исследования [10] был проведен комплекс работ по моделированию процессов в плазменных горелках, применяемых для сварки и наплавки с экспериментальной проверкой, полученных в ходе расчетов рекомендаций.
Модели, основанные на МГД приближении, дают более реалистичную картину течения в каналах плазменного оборудования. Однако это их преимущество достигается за счет высоких вычислительных затрат. Для случая проектировочных расчетов оборудования, особенно в задачах оптимизации, это становится серьезным недостатком.
В работе [11] проведено сравнение результатов моделирования течения газов при использовании МГД приближения и при задании дуги в
виде объемного источника тепла. Результаты расчетов показывают, что учет МГД эффектов приводит к большим температурам и скоростям течения на оси канала (рис. 1). Расхождение результатов моделирования увеличивается с ростом силы тока дуги. Однако, как было показано в работах [2, 3], активные газы проникают в при-катодное пространство в основном за счет диффузии вдоль стенок канала, где температуры невелики и оценка параметров течения рассматриваемыми методами дает близкие результаты. Рост температуры на оси канала и соответствующий рост величины коэффициентов диффузии компенсируется существенным увеличением скорости потока газа.
Таким образом, на этапе проектировочных расчетов при определении рекомендуемых режимов подачи защитного газа на основе оценки максимальных уровней парциального давления активных газов вблизи термокатода, применение упрощенных моделей без учета МГД эффектов может быть оправдано. Более сложные модели целесообразно применять на этапе поверочных расчетов, используя данные моделирования на упрощенных моделях в качестве начальных условий. Для реализации такого подхода необходимо усовершенствование традиционных моделей, использующих задание дуги в виде объемного источника тепла, с целью учета вихревого характера течения вблизи стенок канала.
Описание тестовой задачи
Продолжая исследования, результаты которых описаны в работах [2, 3], будем рассматривать течение в катодном узле с полым катодом. Плаз-мообразующий газ (аргон) подается через полый катод с расходом и через коаксиальный зазор с расходом 62 между катодом и стенкой канала (рис. 2). Подача газа через полость катода производится аксиально, через зазор — с возможностью закрутки. Дуга (на рисунке затонирована) горит между полым катодом и плоской поверхностью,
— дуга в виде объемного источника тепла;
— модель в МГД постановке
Рис. 1. Распределение осевой компоненты скорости и температуры на оси струи, выходящей из плазменной горелки [11]
расположенной перпендикулярно оси катодного узла. Стенки канала катодного узла водоохлаж-даемые.
В тестовой задаче не учитывается теплообмен между потоком газа и стенками, стенки полого катода и нагреваемой поверхности считаются адиабатическими, для стенок канала катодного узла задано условие постоянной температуры.
Газ внутри катодного узла моделировался как идеальная двухкомпонентная смесь (воздух + аргон). Учитывалась нелинейная зависимость характеристик компонент от температуры (рис. 3-6). Данные для расчета были получены путем интерполяции табличных значений [12, 13].
г/Д
КЧЧЧЧЧЧЧЧЧЧ1
г, Л 6.0 4,5 3,0 1,5
к\чч\\\чч\м
При задании интенсивности объемного источника тепла, имитирующего электрическую дугу, учитывались две компоненты: тепловыделение за счет джоулева нагрева и излучение, причем плазма считалась оптически прозрачной. Тепловыделение задавалось по методике, предложенной С. В. Дресвиным [4]. В качестве исходных данных задавался ток дуги и ее диаметр.
Напряженность поля, принимаемая по сечению дуги постоянной, по длине дуги согласно [4] определялась по формуле:
Е =
I
(1)
Рис. 2. Схема тестовой задачи
где I — ток дуги; ст - электропроводность плазмы; ^(г) — площадь сечения дуги.
Интенсивность объемного источника энергии задавалась выражением:
Ж =аЕ2 - пгаа, (2)
где ыгш1 — удельная мощность излучения (Вт/м3).
Рис. 3. Зависимости теплоемкости компонент от температуры
аргон
Рис. 4. Зависимости динамической вязкости компонент от температуры
10000 15000 воздух;
20000 25000 аргон
Рис. 5. Зависимости теплопроводности компонент от температуры
Рис. 6. Зависимость коэффициента диффузии компонент от температуры
Электропроводность (рис. 7) и удельная мощность излучения плазмы (рис. 8) нелинейно зависят от температуры. При моделировании электропроводность задавалась согласно данным [12], а значения удельной мощности излучения для аргона и воздуха определялись по данным [14, 15] и задавались при помощи интерполяции.
Описание численной процедуры решения тестовой задачи
Рассматриваемая задача является достаточно сложной из-за существенной нелинейности свойств компонент газовой смеси. Дополнительной сложностью является необходимость создания поверхностей раздела расчетной сетки, вызванная тем, что при традиционном подходе для задания объемных источников энергии в СББ пакетах требуется выделение в расчетной области трехмерных подобластей.
Сшивка сеток по этим подобластям происходит при помощи различных алгоритмов, однако эта операция приводит к появлению дополнительных вычислительных погрешностей, требует повышенных вычислительных затрат и в ряде случаев может приводить к численной неустойчивости при решении задач, особенно в случае резкого изменения параметров течения вблизи поверхности слияния сеток.
Поэтому в настоящей работе была изучена возможность другого способа задания объемных источников, основанного на применении метода ^-функций. Данный метод, впервые предложенный В. Л. Рвачевым, позволяет строить уравнения ю(ж, у, г) для областей О произвольной геометрической формы, обладающие следующими свойствами:
ю( х, у, г) =
ю > 0, (х, у, г) е О ю = 0, (х, у, г) е дО , ю < 0, (х, у, г)
(3)
Для построения такого уравнения ю(ж, у, ¿) достаточно описать О в виде выражения, образованного операциями алгебры логики с некоторыми опорными областями /¡(ж, у, г) > 0 и заменить операции алгебры логики соответствующими ^-функциями [16]. Наиболее распространенной системой ^-функций, которая применяется для этих целей, является следующая:
2 2
хлу=х+ух +у ;
2 2 хуу=х+у+у
После построения указанным способом уравнения области О, в которой действует объемный источник, для задания его интенсивности Ж можно воспользоваться, например, следующей зависимостью:
Ж = Ж * Н(Н (ю) -1/4),
(4)
где Н — функция Хевисайда.
Легко убедиться, что с учетом (2) и (3) выражение (4) будет давать следующие значения для Ж
Ж =
оЕ -ыга^,(х,у,г) еОудО 0, (х, у, ¿) £ О .
где О — граница области О.
При этом нет необходимости в выделении каких-либо подобластей при построении сетки конечных элементов. Очевидно, что такой подход может быть распространен и на случаи любых поверхностных, линейных и точечных источников.
Для проверки эффективности предложенного подхода тестовая задача решалась двумя способами: при традиционном задании объемного источника с выделением подобласти и сшивкой сеток и предлагаемым методом, основанном на применении ^-функций.
х = -х
■ воздух;
аргон
воздух;
аргон
Рис. 7. Зависимости электропроводности компонент от температуры
Рис. 8. Зависимости удельной мощности излучения компонент от температуры
В качестве начальных условий в области, занятой дугой, задавалась постоянная температура, равная 10000 К для обеспечения среде некоторой начальной электропроводности, необходимой для начала итерационного процесса. При задании начальных условий использовались те же подходы, что и при задании интенсивности объемного источника.
Результаты моделирования и их обсуждение
При решении задачи с традиционным способом задания объемного источника итерационный процесс несколько раз прерывался из-за вычислительных ошибок. В результате решение удалось получить, последовательно решая три задачи: для действия объемного источника постоянной по области дуги интенсивностью 108 Вт/м3; для действия объемного источника постоянной по области дуги интенсивностью 109 Вт/м3; для действия объемного источника с интенсивностью, задаваемой формулой (2). При этом в качестве начальных условий для каждой последующей задачи использовалось решение предыдущей. Для получения решения потребовалось более 6000 итераций.
При задании источника при помощи формулы (4) решение на аналогичной сетке и с тем же, что и в первой задаче, размером шага по времени решение было получено менее чем за 2000 итераций (рис. 9). Решения, полученные при различном способе задания объемного источника, количественно и качественно совпадали. Таким образом, для тестовой задачи предложенный способ задания объемного источника на основе применения метода Я-функций оказался более эффективным с точки зрения устойчивости вычислительного процесса и времени получения решения.
В дальнейшем было проведено сравнение результатов по определению состава атмосферы в канале катодного узла, полученных по результатам холодных продувок [3] и с учетом влияния дуги на параметры течения и свойства компонент газовой смеси. При одинаковых геометрических параметрах задачи и расходах инертного газа результаты существенно отличались. Парциальное давление воздуха вблизи катода для задачи с учетом влияния дуги составило более 10 Па (рис. 10), что превышает критический уровень для всех известных эмиссионных материалов. Для случая холодной продувки уровень парциального давления на срезе катода не превышал величины 10-4 Па.
В работе [4] показано, что для случая холодной продувки катодного узла с полым катодом ламинарным потоком инертного газа парциальное давление воздуха в канале обратно пропорционально величине ехр(Уг/П), где V — скорость вдоль оси канала; Б — коэффициент взаимной диффузии аргона и воздуха.
Сравнение распределения указанной величины для случая холодной продувки и течения с горящей дугой показало, что величина комплекса ехр(К2/Д) для рабочего режима существенно ниже, чем для холодного течения (рис. 11). Это вызвано более быстрым ростом величины коэффициента диффузии по сравнению с осевой скоростью в диапазоне рабочих температур (300 12000 К).
Обеспечение эффективной защиты термоэмиссионного полого катода при рабочих режимах возможно за счет увеличения расхода инертного защитного газа или за счет изменения способа его подачи. Для достижения допустимой величины парциального давления воздуха (»10-3 Па
N
Рис. 9. Распределение температур и векторы скорости течения для тестовой задачи
для современных эмиссионных материалов) в тестовой задаче для этого необходимо увеличить
значение ехр (Уг/П) в 10000 раз, что соответствует десятикратному росту осевой скорости и расхода аргона соответственно.
100 :Р,Па
10 1=4
1,0
0,75 0.5 0,25 0 0,25 0,5 0,75 г/К
Рис. 10. Распределение парциального давления воздуха в канале с учетом влияния дуги
: Vz
-CD
1Е4
* * l l i 1ЕЗ ч % \ 1 I f
l ■ 1Е2 1 1
0,75 0,5 0,25 0 0,25 0.5 0,75 r/R
Рис. 11. Распределение величины exp(Vz/D)в канале для сечения z/R = 6,25
Однако увеличение расхода аргона приводит также к снижению температуры газа в канале вблизи катода, а, следовательно, и к снижению коэффициента диффузии аргон-воздух. Поэтому указанного уровня парциального давления для тестовой задачи удалось достичь при увеличении расхода аргона через коаксиальный зазор в 4,67 раза (с 1,5х10-5 до 7х10-5 кг/с). На рис. 12 приведены графики распределения парциального давления воздуха в канале для этого расчетного случая.
На рис. 13 приведено распределение величины exp(Vz/D) в сечении z/R = 6,25 при увеличенном расходе аргона. Значение этого комплекса выросло только у стенки канала (г/R = 0,75...1), однако это обеспечило снижение парциального давления воздуха до необходимого уровня.
Р,Па -и»
1
Ш-1
1E-2
/ 1E-3 f V s \
0,75 0,5 0,25 0 0,25 0,5 0,75 г/R
Рис. 12. Распределение парциального давления воздуха в канале с увеличенным расходом аргона
: V7
;eD
1E4
I 1 * u f\ 1E3 N / % 1 u
1 \ 1E2 у \
0,75 0,5 0,25 0 0,25 0,5 0,75 r/R
Рис. 13. Распределение величины exp( V/D) в канале для сечения z/R = 6,25 с увеличенным расходом аргона
Это подтверждает сделанный в [2—4] вывод о том, что поступление воздуха в прикатодное пространство в катодных узлах рассматриваемой схемы происходит в основном за счет диффузии вдоль стенок канала.
В работе [2] определены допустимые режимы подачи защитного газа в катодных узлах с полым катодом исходя из условия недопущения циркуляционных вихрей у среза катода. Установлено, что такие вихри неизбежно образуются при значении Re^Re^ > 1,7, где Re2 — число Рейнольдса для потока, проходящего через зазор между стенками катода и канала, а Rej — через полый катод.
Для рассматриваемого случая, несмотря на значительное увеличение расхода через коаксиальный зазор из-за крайне неравномерного распределения температур это соотношение составляет 12,434/15,979 » 0,778.
Для оценки влияния способа подачи защитного газа на картину течения и состав атмосферы на рабочих режимах были проведены допол-
нительные численные эксперименты. Рассматривался комбинированный способ подачи инертного газа - аксиальный с постоянным расходом через полый катод и вихревой через зазор между стенками катода и канала катодного узла.
Установлено, что повышение частоты вращения подаваемого в коаксиальный зазор газа вплоть до 200 Гц существенного влияния на состав атмосферы не оказывает. В большей части канала течение остается аксиальным и ламинарным (рис. 14).
Рис. 14. Линии тока для случая комбинированной подачи газа
Выводы
1. На примере тестовой задачи для катодного узла с полым катодом проведена оценка вычислительной эффективности двух способов задания интенсивности объемного источника энергии — традиционного, основанного на выделении трехмерных областей и последующей сшивкой сеток и основанного на методе R-функций. Показано, что в последнем случае вычислительный процесс более устойчив, а время решения сокращается более чем в 3 раза.
2. Показано, что для правильного назначения режимов подачи защитного газа на рабочих режимах необходимо учитывать влияние дуги. Назначение величины расхода по результатам расчетов состава атмосферы при холодной продувке не обеспечивает необходимого уровня защиты. Требуемая величина расхода может быть откорректирована на основе учета значения комплекса
exp (VZID) .
3. На рабочих режимах плазменного оборудования состав атмосферы в прикатодном пространстве определяется процессом диффузии воздуха вдоль стенок канала катодного узла. Обеспечение заданного состава атмосферы вблизи эмиттера возможно за счет увеличения осевой скорости потока защитного газа у стенки канала (r/R = 0,75...1).
Перечень ссылок
1. Слободянюк В. С. Моделирование вихревых и турбулентных явлений в электродуговых устройствах: дис. д-ра техн. наук : 01.04. 14 / Слободянюк Валерий Сергеевич. — Бишкек, 1996. — 376 с.
2. Планковский С. И. Газодинамическое проектирование катодных узлов в интегрированных CAD/CAE системах / С. И. Планковс-кий // Радюелектронт i комп'ютерш систе-ми. — 2008. — № 4. — С. 60—65.
3. Кривцов В. С. Опташзащя газодинамiчних характеристик дугових плазмотронв з порож-нистим катодом / В. С. Кривцов, С. I. План-ковський, 6. В. Цегельник та ш. //Hауковi вкт НТУУ «КП1». — 2006.— № 3. — С. 106—113.
4. Планковский С. И. Диффузия атмосферного воздуха через встречный ламинарный поток инертного газа в плазмотроне / С. И. Планковский, Е. В. Цегельник, Е. К. Островский // Вестн. Нац. техн. ун-та «ХПИ». Тематический выпуск : Проблемы совершенствования электрических машин и аппаратов. — 2007. — № 25 — С. 72—84.
5. Дресвин С. В. Основы математического моделирования плазмотронов. Ч. 1 : Уравнение баланса энергии. Метод контрольного объема.
Расчет температуры плазмы / С. В. Дресвин, Д. В. Иванов : Учеб. пособие. — СПб : Изд-во Политехн. ун-та, 2004. — 227 с.
6. Дресвин С. В. Основы математического моделирования плазмотронов. Часть 2 : Электромагнитные задачи в плазмотронной технике/ С. В. Дресвин, Д. В. Иванов : Учеб. пособие. — СПб : Изд-во Политехн. ун-та, 2006. - 295 с.
7. Нгуен Куок Ши. Исследование индукционных и дуговых плазмотронов: автореф. дис. д-ра техн. наук : 05.09. 10 / Нгуен Куок Ши. — Санкт-Петербург, 2002. — 40 с.
8. He-Ping Li. Three-dimensional modelling of a DC non-transferred arc plasma torch / He-Ping Li, Xi Chen // Journal of Physics D : Applied Physics. — 2001. — Vol. 34. — N. 17. — P. L99— L102.
9. Trelles J.P. Multiscale Finite Element Modeling of Arc Dynamics in a DC Plasma Torch/ Juan Pablo Trelles, Emil Pfender, Joachim Heberlein // Plasma Chemistry and Plasma Processing. — 2006. — Vol. 26. — N. 6. — P. 557—575.
10. Schnick M. Simulation of plasma and shielding gas flows in welding and cutting arcs with ANSYS CFX / Michael Schnick, Uwe Fbssel, Jцrg Zschetzsche // International Scientific Colloquium Modelling for Material Processing. — Riga, June 8-9, 2006. — P. 143—148.
11. Kim Y. Numerical analysis of flow characteristics of an atmospheric plasma torch/You-Jae Kim, J.-G. Han,Youn J. Kim // Center for advanced plasma surface technology Sungkyunkwan university — 2005. — 17 p.
12. Murphy A.B. Transport coefficients of Argon, Nitrogen, Oxygen, Argon-Nitrogen, and Argon-Oxygen Plasmas / A.B. Murphy, C.J. Arundell // Plasma Chemistry and Plasma Processing. — 1994. — Vol. 14. — N. 4. — P. 451—490.
13. Murphy A.B. Transport coefficients of Air, Argon-Air, Nitrogen-Air, and Oxygen-Air Plasmas / A.B. Murphy // Plasma Chemistry and Plasma Processing. — 1995. — Vol. 15. — N. 2. — P. 279—307.
14. Menart J. Net emission coefficients for argon— iron thermal plasmas / J. Menart, S. Malik // Journal of Physics D : Applied Physics. — 2002. — Vol. 35. — P. 867 — 874.
15. Naghizadeh-Kashani Y. Net emission coefficient of air thermal plasmas / Y. Naghizadeh-Kashani, Y. Cressault, A. Gleizes // Journal of Physics D : Applied Physics. — 2002. — V. 35. — P. 2025— 2934.
16. Рвачев В. Л. Геометрические приложения алгебры логики / В. Л. Рвачев. — К. : Техшка, 1967. — 212 с.
Поступила в редакцию 10.11.2009
S. I. Plankovsky Doctor of Engineering Science, Ye. V. Tsegelnyk, V. O. Garin
DESIGN OF FLOWS IN CATHODE KNOTS OF PLASMA EQUIPMENT AT DESCRIPTION OF ARC BY VOLUME POWER SOURCE
Розглянуто питания числового моделювання складу атмосфери в катодных вузлах плаз-мового обладнання з використанням сучасних CFD пакет1в. Для проектувалъних розра-хунк1в обТрунтовано використання моделей з завданням дуги у вигляд1 об'емного джерела енергИ. Проведено пор1вняння обчислювалъног ефективност1 при використання двох вар— ант1в побудови розрахунковог модели: з видленням подобластей в розрахунковш зош та з завданням ттенсивност1 джерела з використанням методу R-функцш. На приклада тес-товог задач1 показано, що алгоритм з використанням R-функцш е быъш сталим та ефективним з точки зору часу розв 'язання задач1. Проведено пор1вняння резулътат1в роз-рахунку складу атмосфери в катодному вузл1 з порожнистим катодом з урахуванням та без урахування впливу дуги на характер течи.
Катодний вузол, математичне моделювання, об'емне джерело, метод R-функцш
There are reviewed questions of numerical modeling of the atmospheric composition in cathode assembly units of electric-arc plasma equipment by means of modern commercial CFD packages. There is substantiated application of arc-setting models in the form of volume heat source for the development of designs of cathode assembly units of electric-arc plasma equipment. There was compared efficiency of two numerical models: the first one —selection of subareas in calculation area, the second one — setting of the source intensity by method of R-functions. As an example of model problem, it is shown that the algorithm of finite element pattern indivisible into subareas has higher stability and efficiency in from the point of view of solution time. There are compared results of calculations of the atmospheric composition in hollow cathode assembly unit taking and not taking into account arc influence upon the flow pattern.
Cathode assembly unit, mathematical modeling, volume source, method of R-functions