CC BY
27
УДК 537.523
Б.Д. Цыдыпов
НЕЛИНЕЙНАЯ ТЕРМИЧЕСКАЯ ЗАДАЧА ДЛЯ СИСТЕМЫ СОПРЯЖЕННЫХ ЭЛЕМЕНТОВ. МЕТОД РЕШЕНИЯ
В работе представлен метод решения нелинейной тепловой задачи для осесимметричного катодного узла генераторов низкотемпературной плазмы. Тепловая задача основана на решении двумерного уравнения нестационарной теплопроводности с нелинейными граничными условиями для системы «вставка - обойма» с учетом основных видов энергообмена с внешней средой.
Ключевые слова: оптимизационная задача, тепловой поток, энергообмен, генератор, низкотемпературная плазма.
B.D. Tsydypov
NONLINEAR THERMAL PROBLEM FOR THE SYSTEM OF CONJUGATE ELEMENTS. METHOD OF SOLUTION
The article deals with the method of solution nonlinear thermal problem for the axis-symmetrical cathode assembly in low-temperature plasma generators. The thermal problem described above reduces to a nonstationary heat conduction equation with nonlinear boundary condition in the two-dimensional unsert -sleeve system with consideration of main types of heat transfer between the electrode assembly and the ambient medium.
Key words: optimization problem, heat flux, energy transfer, generator, low-temperature plasma.
Введение
В предыдущей работе [1] обоснована и приведена математическая постановка задачи о теплофизическом состоянии составных катодных узлов сильноточных плазменных систем (СПС). Она основана на решении нелинейного уравнения нестационарной теплопроводности
dT
ср— = div(^gradT) + qV dt
(1)
с учетом основных видов теплообмена электродного узла с внешней средой, где Т - температура, с -удельная теплоемкость, р - плотность материала, X - коэффициент теплопроводности, - объемная
плотность внутренних источников и стоков, обусловленных различного рода физико-химическими процессами выделения и поглощения энергии.
Температурное поле Тк(г, ¿) в составном катодном узле цилиндрической симметрии (см. рис. 1 [1]) находится совместным решением уравнения (1) в виде:
ск Рк
dT.
dt
1_d_ r dr
З (т\ dTk гЛк(Т )-к-
dr
d + —
dz
+ Л/ (ук (т )
(2)
и уравнения непрерывности тока
1 d r dr
У (Т) 1Г
dr
d + — dz
Ук (Т)dU
dz
= 0
(3)
для двух сопряженных элементов конструкции: вставки (к = 1) и обоймы (к = 2) с нелинейными граничными условиями I - IV родов комбинированного энергообмена. Здесь
]к = (Л2 + Л2)1/2, Л =~ак (Т)дик / дг, ]1 =~ок (Т)дик / dz - плотность тока и ее компоненты по
цилиндрическим координатам г и z, ик - потенциал электрического поля, ск (Т) - удельная электрическая проводимость. Система уравнений нелинейна, тепло- и электрофизические коэффициенты Хк (Т) и ск (Т) зависят от температуры, являющейся функцией координат.
Метод решения
Тепловая задача в сформулированной постановке [1] аналитически не решается. Поэтому используется численный метод - метод конечных разностей [2]. Нестационарное уравнение теплопроводности (2) решается методом установления (стационирования). Для придания алгоритму
решения универсальности следует перейти к безразмерным переменным. После обезразмеривания уравнение (2) запишется следующим образом:
дУ _ 77 I (д 2 У V 77 I ( 1 дУ , д 2 У
2
дх
дт = г
II
О z
vdz / ( д2 х ^
Г дГ дГ2
+ ^ (Т) (4)
2
VdZ у
+ F 11
+ Г0г
(1 дх + д2 х ^ V Г дГ дг2 у
+ С ~]22р2 (Т) (5)
соответственно для элементов I и II электродной структуры. Здесь
Ъ0z = Хт0 / срИ,,, Ъог = Хт0 / срЯ2 , G0 = }02т0 / а0рсТ0 - безразмерные параметры; у = Т^/То, х = Т2/Т0 -относительные температуры в элементах; т = г/т0, Z = z / Ц, Г = г / ^2 - безразмерные координаты;
т0, Ц, - масштабы координат.
В дальнейшем индекс «тильда» в переменных опускаем, используя в формулах безразмерные величины. Область интегрирования I - II, занимаемую катодным узлом, разбиваем пространственной сеткой
_ _ Г ^ = Н,, К> 0, ; = 0,1,2,...М; Щ = Ц
со, = со, ха, \
Н \гк = кН2, Н2 > 0, к = 0,1,2,..М; МН2 = Я2
с шагами Н1 по координате z и по г. По временной координате введем равномерную сетку а ={ г} = }т0, т0 > 0, } = 0,1,2,...} с шагом т0 .
Значение сеточной функции у (или х) в некотором узле сетки гк, г} ) обозначим у^к.
Соответственно у(Zi±^ гк±Р г} ) = у}к^ у(z¡, гк , г}+1) = у/д1.
На 4-х точечном шаблоне расчетной пространственно-временной сетки дифференциальное уравнение (4) аппроксимируется разностным уравнением:
у}+ - у/,'к = Ъz (уЙ1 - 2у}+1 + у^С!) /(Н / Ц )2 + ^ ( уЛ:1 - у}+1) / 2к(Й2 / ^)2 +
ч (6)
+^ (у}-1 - 2 у}+1 + 2 у}+;) / (V *2)2 + ^} (у/) / *( у}),
где ст(у/) - удельная проводимость на }-том временном слое. Погрешность аппроксимации имеет
порядок 0(Н2+т) [2]. Применяя локально-одномерную схему прогонки метода дробных шагов [3],
уравнение (6) разбиваем на два одномерных
у}+1 - у} = 2*0z (у}1 - 2у/+1 + у/+1) / (Н / Ц )2; (7)
у}+1 - у}=2^ ((у/;; - у}-1 ) / 2кн / ^)2+(у}';; - 2у}+1+2у}++;) / (н / ^)2)+
(8)
+2^} (у/)/*( у})
соответственно по координатам z и г. При прогонке разностные уравнения (7) и (8) приводятся к алгебраической системе уравнений типа
А у;-1 - С; у; + ^ ум = - ^ (9)
с условиями А;, В; > 0, С; > А;+В;, обеспечивающими разрешимость системы методом прогонки. Решение задачи ищем в виде
у; =^+1 у;+1 ;Дчl, ; = 0,1,2,...М -1, (10)
где коэффициенты а;+1 и в ¡+1 вычисляются по реккурентным соотношениям
а+1 = В; / (С; - Аа), в+1 = (ъ - А в)/ (С; - Аа). (11)
Значение начальных прогоночных коэффициентов а1 и в определяется при помощи одного из граничных условий области интегрирования. Затем из второго граничного условия находим значения сеточной функции уМ и по формуле (10) вычисляем все остальные значения у; вплоть до у0. При этом переход от временного слоя } на слой }+1) осуществляется последовательным решением однородных уравнений (7) и (8) по соответствующим координатам.
Так как катодный узел состоит из сопряженных элементов разной геометрии, необходимо разбить его на несколько простых областей и для каждой решать свой прогоночный цикл, «сшивая» соседние
области удовлетворением единым граничным условиям. Для прогонки по координате z выделяются области ОАЪК и ВСЕЪ, а по координате г - области ОАВО, ОСБЬ и ЬВЕК (рис. 1 [1]). Следует отметить особенность прогонки по областям ОАЪК и ОСВЬ. Здесь используется метод встречных прогонок [2]: из граничных условий вычисляются начальные коэффициенты а1 и Д1 для первого элемента, а1 и Д1 для второго элемента, а искомая функция на их границе определяется из условий сопряжения (14-16)
[1].
Рассмотрим решение уравнений (7) и (8) по координатам z и г.
Расчет по координате z
Разностное уравнение (7) преобразуем к виду:
Foz(VЦ)-2у—1 -(2 + 2^(VЦ)-2)у/+1 + Гz(Й1 /Ц)-2у/+ =-2у} .
Обозначая А = В1 = Г0z(Н1/ Ц)-2 = Г0N2, С; = 2 + 2Г0^2,Г = 2у?, получим систему разностных уравнений вида (9), которая решается для каждой области интегрирования. Поскольку вычислительная процедура прогоночных коэффициентов и граничных значений функций достаточно громоздка, ниже приведем их последовательный расчет лишь для сложной области ОАЪК. Расчеты по другим областям проводятся аналогично данной области.
Область ОА¥К
1. Левое граничное условие (охлаждаемая поверхность элемента II). Сравнивая у0 +1 = а1 у/+1 + Д с условием (6) [1], находим = 0, = 1.
2. Правое граничное условие (активная поверхность элемента I). После обезразмеривания условие (4) [1] запишется следующим образом:
у}+1 - у-1 = МЛ-1^ ;Яl-1£laLsTo2(y0)3у0+1, 0< г < г0 (12)
2Н1 / Ц Ца^41 + Я;1£1аЬХ(у0)3у0'+1, г0 < г<^1 .
Здесь удовлетворение граничных условий со вторым порядком точности приводит к необходимости постулировать выполнение уравнений и в граничных точках (расширенной зоне). Вводится неизвестное значение искомой функции у -1 в фиктивной точке. Исключая значение у -1 из уравнения (12) при помощи
Л у-1 - С0 у0 +1 + В0 у1}+1 = -Г0 , у0 +1 = а1 у1+1 + Д ,
определяем начальные прогоночные коэффициенты:
а = |(1 + (V Ц )2/Г£ +Л-1£1^Н1Т03( у0)3)-1, 0 < г < г ;
1 1(1+(*1 / ц )2 / г^+л^мт,3 (у0 )3+к1ан )-1, г < г < ^1;
в =
• - 01
у0+1(У Ц )2/ ^ + дА%%
0 < г < гп
1 + (V Ц )2/^ +Лі-1ЄіаИ1Го3( у0)3
, Го < г < ^
1 + (к /Ц)2 /+^1^1^й1Го3(уО)3 + КтН
3. Условие на границе сопряжения в разностной форме запишется:
у
с другой стороны
Ул+=у1:', кк(2к1)-1(у;+1+уц+\)=к4ад-чх;+1 -4+-1), (13)
уМ}+-1 = + ДЛ'1 , уМ}+-1 - С№1 уЛ'+1 + В№1 уМ2 +1 = -Г№2 ; (14)
+1 = у а11 + Дп Ап +1 — С11 +1 + V* = — Г11 (Л
хЩ -1 = уМ2аМ2 + ДМг, АМ12 хМ2 -1 С М2 хМ2 + ВМ2 хМ2 +1 = Г М2 . (15)
Решая совместно уравнения (13-15), получим граничные значения искомых функций:
\2/
7+1 = х+1 = в + Кв + (к / Ц )2(К / ^ +у О у^ Ущ Хщ К + К - а - К<2 + (к / ц)2(к / <+к / ^) ■
Расчет по координате г
После несложных преобразований уравнение (8) приводится к виду:
^(1 -1/2к)у£ -2(^ + М-2)уГ1 + ^(1 + 1/2к)< =
= Со М -2 ]2( ук )а\ у]к) + 2М -2 у1
или
Аук-1 - Ску1+1 + ад*!1 —-^',
Ак = ^ (1 - 1/2к), Ск = 2(^ + М-2), Бк — ^(1 +1/2 к), (16)
^ = Со М -2 г 2( ук )о--'( ук) + 2М -2 у/.
Аналогично ранее записанному соотношению по z решение уравнения (16) ищем в форме
y k+1 =ак+1 Уш +Ä+1. где
ак+1 =(1 + 1/2k)/(2 + 2F0;1M-2 -(1 - 1/2k)ak),
Ä+i = (For (1 -1/ 2k)ßk + Go jc2M-2a + 2M 2yj+1) / (2^ + 2M 2 - (1 -1/ 2*)^ )•
Область OABG
1. На оси симметрии (k = 0) радиальный тепловой поток равен нулю. Следовательно, в уравнении
1 дТ
(2) слагаемое------при r^0 представляет собой неопределенность типа 0/0. Раскрытие ее по правилу
r дг
Лопиталя
1 дТ д 2Т
lim-------= —-
r^0 r дг дг
дает следующее разностное выражение:
У0+1 - У0 = R2 V2 F^r (У-+1 - 2 y0+1 + y1+1)+G0 jc 2( У0 )о-'( у0 )•
Полагая y_1 = y1, имеем:
y0+1 = (1 - M 2/2F0r )-1 у/+1 = (yj + G0 jc 2( У0 )/2g( y0 ))/(1 + 2F0.M2).
Сравнивая с у0 +1 = a1 yj+1 + ß1, находим прогоночные коэффициенты на оси симметрии:
а = (1 -M2 / 2F0r)-1, ß1 = (y0 + G0jc2(У0)/ 2g(y0)) / (1 + 2FürM2) •
2. Определение граничной температуры yM+1 на цилиндрической поверхности АВ сводится к решению следующей системы уравнений:
(yM, +1 - yM+ч) / (2h2 / R2) = -^1-1а^2yM+1 - y]Ml )3yM!1
A1 yj+1 - ^1 yj+1 + ft1 y* = -Fj
AM1 yM1 -1 ^M1 yM1 + nM1 yM1 +1 1 M1'
После преобразований окончательно получим:
j+1 _ F0:rß^1 + G0 jc2 (yM 1) / (2M 2^( y^ 1)) + M 2 yM 1
/М‘ ^ + М-2 - ^оХ + (0,5^1-1^)1г)(ак1 + £1ак2То3(у, )3)
Расчеты по областям СОЬС и ВЕКЬ проводятся аналогично расчетам по области ОАГК, показанным при прогонке по координате г.
Рассмотренный выше метод применяется также для решения уравнения (3) с соответствующими граничными условиями [1]. Для этого в алгоритме значения е1,2) а, ]с и То необходимо приравнять к нулю, а Х(Т) и Т(г,г) заменить соответственно на о(Т) и и(г,г). Таким образом, задача о тепловом состоянии катодного узла решается в следующем порядке. Сначала из уравнения неразрывности тока находится поле электрического потенциала и(г,г) и определяется распределение плотности тока по формулам:
Гк — Г + Л2)1/2, Л —■О(Т)дик /дг, ]г — -ок(Т)дик /дг.
Затем, используя полученное распределение тока в объеме электрода, решается уравнение теплопроводности и вычисляется поле температур Т12(г, г).
Заключение
Составленный численный алгоритм позволяет рассчитать стационарные температурные поля Tia(r, z) во всей электродной системе «вставка - обойма» цилиндрической симметрии. В постановке корректно учтены такие факторы, как двумерность задачи, нелинейность граничных условий, зависимость тепло- и электрофизических свойств материалов сопряженных элементов от температуры, неравномерность распределения тока в объеме конструкции. Особенностью алгоритма является определение стационарного температурного поля решением нестационарного уравнения теплопроводности методом установления. Эволюционный характер задачи позволяет при решении уравнения (2) по экономичной локально-одномерной схеме прогонки произведение ckpk принять постоянным, хотя оно в общем виде является функцией температуры, координат и т.д. При изучении сугубо нестационарных и переходных режимов системы такое упрощение исключается. В последующей работе будут представлены результаты математического моделирования по разработанному алгоритму.
Литература
1. Цыдыпов Б.Д., Баргуев С.Г. Постановка нелинейной термической задачи для сопряженных элементов // Вестник БГУ. 2010. - Вып. 9. - С. 189-193.
2. Самарский А. А. Теория разностных схем. М.: Наука, 1977. - 656 с.
3. Тихонов А.Н., Самарский А. А. Уравнения математической физики. М.: Наука, 1977. - 736 с.
Цыдыпов Балдан Дашиевич, кандидат технических наук, старший научный сотрудник Отдела физических проблем Бурятского научного центра СО РАН.
Tsydypov Baldan Dashievich, candidate of technical science, senior researcher of Physics Problems department of Buryat Scientific Centre SB RAS.
