CC BY
23Математика к Математическое
моделирование
ХДК 536.2
Ссылка на статью:
// Математика и математическое моделирование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2016. №5. С. 29-45.
Б01:10.7463/шаШш.0516.0847523
Представлена в редакцию: 23.08.2016 © МГТУ им. Н.Э. Баумана
Вариационный подход к анализу модели теплового взрыва в твердом теле
Зарубин В. С.1'*, Кувыркин Г. Н.1, Савельева И. Ю.1 *[email protected]
1МГТУ им. Н.Э. Баумана, Москва, Россия
Предложен и обоснован вариационный подход к количественному анализу одно- и двухпараме-трической математических моделей теплового взрыва в твердом теле, опирающийся на двойственную вариационную формулировку нелинейной задачи стационарной теплопроводности, содержащую два альтернативных функционала, принимающих совпадающие значения в своих стационарных точках, соответствующих истинным распределениям температуры. Этот подход в сочетании с использованием приближенных распределений температуры, допустимых для указанных функционалов, позволил установить предельные сочетания параметров, соответствующие состоянию теплового взрыва в таких телах канонической формы, как длинный круговой цилиндр и пластина, неограниченная в своей плоскости. Близость полученных результатов к известным точным решениям дает основание считать целесообразным применение данного подхода к анализу состояния теплового взрыва в телах произвольной формы.
Ключевые слова: математическая модель теплового взрыва; двойственная вариационная формулировка нелинейной задачи; альтернативные функционалы
Введение
Температурное состояние твердого тела может зависеть как от условий теплообмена с окружающей его поверхность внешней средой, так и от выделения энергии в объеме этого тела, вызванного, например, протеканием процессов в элементах ядерного реактора или экзотермических химических реакций [1], при поглощением энергии проникающего излучения или переходом в теплоту части электрической энергии при прохождении электрического тока (так называемая джоулева теплота) [2, 3]. Если интенсивность объемного энерговыделения возрастает с увеличением температуры, то возможно возникновение предельного установившегося температурного состояния, при котором отвод к поверхности тела выделившейся в его объеме тепловой энергии достигает максимума. При этом малые приращения температуры приводят к увеличению выделения тепловой энергии, которую уже нельзя отвести к поверхности тела путем теплопроводности без дальнейшего возрастания температуры. В итоге установившееся распределение температуры в теле становится невозможным, что и
определяет состояние теплового взрыва [1], получившее такое название в силу того, что соответствующая математическая модель предсказывает в этом случае неограниченное возрастание температуры.
Анализу состояния теплового взрыва посвящено достаточно много работ, связанных с исследованием процессов горения и взрыва в неподвижной среде и проанализированных в монографиях [1,4]. В большинстве известных работ рассматривают математическую модель, описывающую распределение температуры в случае, когда энерговыделение вызвано экзотермическими химическими реакциями, скорость протекания которых возрастает с увеличением температуры. Зависимость скорости химической реакции от температуры обычно описывают экспоненциальным законом Аррениуса, что приводит к необходимости рассматривать существенно нелинейную математическую модель, содержащую дифференциальное уравнение, в которое входит слагаемое, нелинейно возрастающее с ростом температуры. Даже при упрощающих допущениях эта модель позволяет получить точное решение в аналитическом виде лишь в случае одномерных распределений температуры в двух областях канонической формы: в неограниченной в свой плоскости пластине и в неограниченном по длине круговом цилиндре.
Приближенным численным решением дифференциального уравнения, входящего в нелинейную математическую модель теплового взрыва, удается получить количественные оценки сочетания определяющих параметров, при котором наступает предельное состояние в областях не только канонической формы. Возможности исследования состояния теплового взрыва можно расширить в связи с развитием методов математического моделирования [5, 6], в том числе методов анализа моделей, описывающих температурное состояние твердых тел [7, 8].
В данной работе для анализа математической модели теплового взрыва в однородном твердом теле использован вариационный подход, основанный на двойственной вариационной формулировке соответствующей нелинейной задачи стационарной теплопроводности в таком теле [9]. Эта формулировка содержит два альтернативных функционала, достигающих совпадающих значений в своих стационарных точках, соответствующих истинным распределениям температуры. Такое свойство функционалов позволяет не только получить приближенную количественную оценку сочетания параметров, определяющих состояние теплового взрыва, но и установить возможную наибольшую погрешность такой оценки.
1. Варианты математической модели теплового взрыва в твердом теле
Установившееся распределение температуры Т(М) в однородном твердом теле, занимающим область V, при наличии в теле внутренних источников теплоты с объемной мощностью энерговыделения qv (М, Т), зависящей (в общем случае) не только от координат точки М е V, но и от температуры, удовлетворяет дифференциальному уравнению [9, 10]
V ■ (Л(Т) УТ(М)) + qv(М, Т) = 0, М е V, (1)
где V — векторный дифференциальный оператор Гамильтона; Л — коэффициент теплопроводности материала твердого тела, в общем случае зависящий от температуры. Для однозначного решения этого уравнения в точках P Е S поверхности S, ограничивающей область V, должны быть заданы граничные условия. Примем простейший вариант таких условий в виде
T(P) = To = const, P Е S. (2)
Если выделение теплоты в теле вызвано протекающими в его материале экзотермическими химическими реакциями, скорость которых описывает закон Аррениуса [1], полученный с использованием распределения Больцмана [11] реагирующих частиц по уровням энергии, то объемную мощность энерговыделения можно представить в виде [12]
qv (T) = qV exp( - k^T)' (3)
где qV — предэкспоненциальный множитель, имеющий размерность Вт/м3; E — энергия активации реакции; k0 = 1,38 • 10-23 Дж/К — постоянная Больцмана. Для экзотермических реакций qV > 0. Поэтому при фиксированном значении T0 температуры на границе тела для температуры в его внутренних точках справедливо неравенство T(M) > T0, M Е V. По-
dq (T)
скольку E > 0, из формулы (3) следует, что ——— > 0. Это означает, что установившееся распределение температуры в теле возможно лишь при определенном ограничении на величину qv (T), при котором выделившаяся в теле тепловая энергия еще может быть отведена к его границе путем теплопроводности. При нарушении такого ограничения существование установившегося распределения температуры в теле становится невозможным. Именно в таком случае говорят о состоянии теплового взрыва.
Соотношения (1)-(3) входят в математическую модель, которая является основной в стационарной теории теплового взрыва [1]. Эта модель может быть представлена в нескольких вариантах.
Если коэффициент теплопроводности материала тела допустимо считать не зависящим от температуры, т.е. принять Л = const, то для удобства количественного анализа математической модели теплового взрыва уравнение (1) с учетом граничного условия (2) и равенства (3) можно привести к так называемой квазилинейной форме
V 20 + в exp( ^ ) =0, (4)
где V2 = L2V2 — оператор Лапласа, определенный в системе безразмерных координат, нормированных по характерному линейному размеру L области V;
T - To д qVL2 ( 1) koTo
0 = —; в = exP — ; Y
7Т0 Л7Т0 47/ Е
В уравнение (4) входит два безразмерных параметра в и 7, функциональная связь между которыми, устанавливаемая из решения уравнения (4) при нулевом значении безразмерной
температуры 9 на поверхности Б, определяет состояние теплового взрыва. Следует отметить, что уравнение (1) с учетом равенства (3) можно представить соотношением
V Ч + С ехр( - ^)=0, (5)
содержащим лишь один безразмерный параметр С = Ь2к0/(\Е), но поскольку теперь безразмерная температура определена равенством = к0Т/Е, необходимо при анализе математической модели иметь в виду и второй отличный от нуля безразмерный параметр
= к0Т0/Е, который войдет в граничное условие для уравнения (5). Следовательно, и для этого варианта математической модели состояние теплового взрыва будет определять функциональная связь двух безразмерных параметров С и $0.
В теории горения [1], включающей математическую модель теплового взрыва, принимают упрощающее допущение о возможности в уравнении (4) пренебречь в знаменателе показателя экспоненты слагаемым 9по сравнению с единицей, т.е. от двухпараметриче-ской модели перейти к более простой однопараметрической. Это допущение равносильно выполнению условия (Т (М) — Т0^ /Т0 ^ 1, М Е V, и позволяет вместо функциональной связи между двумя безразмерными параметрами состояние теплового взрыва характеризовать значением лишь одного безразмерного параметра в. Именно для такого упрощенного варианта модели получены упомянутые выше решения в аналитической форме для пластины и цилиндра.
В варианте математической модели, содержащей уравнение (5), также применительно к теории горения при условии $0 ^ 1, выполняемом в случае больших значений энергии активации химической реакции, рассматривают предельный случай $0 = 0 [1], что дает возможность и для этого варианта модели характеризовать состояние теплового взрыва значением одного безразмерного параметра С.
Каждому варианту математической модели теплового взрыва, количественный анализ которой связан с необходимостью решения квазилинейного дифференциального уравнения, можно поставить в соответствие вариационную формулировку модели.
2. Вариационная формулировка модели теплового взрыва
Построение вариационных формулировок для рассмотренных выше вариантов модели теплового взрыва можно объединить, если уравнения (4) и (5) представить в единой форме
V 20 + / (0) = 0, (6)
где 0 соответствует безразмерной температуре, а функция /(0) — второму слагаемому в левой части каждого из указанных уравнений.
Умножим уравнение (6) на вариацию 50 и после интегрирования полученного произведения по объему области V запишем
У(V20) 50 ¿V + У f (0) 50 ¿V = 0.
V V
Применим к этому равенству первую формулу Грина в виде [ 10]
У(V20) ¿0 ¿V = У п ■ ('0) ¿0 ¿Б - I '0 ■ '(¿0) ¿V, (7)
V в V
где п — единичный вектор внешней нормали к поверхности Б. Тогда с учетом равенства нулю вариации ¿0 для фиксированный значений 0 на поверхности Б области V получим
О
1
¿/и^0)2 -//(0)^0 ^(М) = 0,
2 ./ /
V 4 ©* 7
где 0* — нижняя грань множества ожидаемых значений 0 в области V (для модели, включающей уравнение (4), 0* = 0, а для модели, содержащей уравнение (5), 0* = $0). Полученное
равенство соответствует условию стационарности функционала
©
3 [0] = / (1(^)2 - // (0) ¿0)
1(У0)2 - I / (0) ¿0) ¿V. (8)
V 4
Этот функционал допустимо рассматривать на множестве непрерышнык в области V и кусочно дифференцируемых в ней функций 0(М), М € V, принимающих на поверхности Б значение 0*.
Для построения функционала, альтернативного по отношению к функционалу (8) и имеющего с ним совпадающие значения в стационарных точках, расширим область определения функционала (8) путем введения векторной функции q, удовлетворяющей дополнительному условию
q(M) + '0(М) = 0, М € V = V и Б, (9)
где 0 — нулевой вектор. После подстановки q = —У0 в функционал (8) запишем
3 [0, ^ = 1 / (т - // (0) ¿0) ¿V, (10)
а затем при помощи векторного множителя Лагранжа Ь(М) (М € V) введем условие (9) в соотношение (10) и получим
3[0, q, Ь] = 3[0, ч] - IЬ ■ (Ч + '0) ¿V. (11)
V
Равенство нулю вариации функционала (11) (с учетом ее преобразования по первой формуле Грина вида (7))
¿3[0, q, Ь] = У - Ь) ■ ¿Ч - (q + '0) ■ ¿Ь + ('V ■ Ь - /(0)) ¿0^ ¿V = 0
V
приводит к условиям стационарности этого функционала, в которые помимо условия (9) входят равенства Ь(М) = q(M) и V ■ Ь(М) = /(0), М € V. Используя первое из этих равенств и соотношения (8), (10) и (11), приходим к функционалу
3'[0, ч] = 3[0] - 11(Ч + '0)2 ¿V (12)
V
с дополнительным условием 0(Р) = 0* и условиями стационарности (9) и
V ■ я(М) = f (0), М Е V. (13)
Из соотношения (12) следует, что 3'[0, я] ^ 3[0], а при выполнении условия стационарности (9) стационарные значения этих функционалов совпадают: 3'[0*, я*] = 3[0*], причем Я*(М) = —У0*(М), М Е V.
Если потребовать выполнения условия стационарности (13) в качестве дополнительного, то после преобразования (с применением первой формулы Грина вида (7)) функционала (12) с учетом функционала (8) и его дополнительного условия 0(Р) = 0* получим
I [0, я] = — I (я2 — |0 ¿0^ ¿V — 0* I я ■ п ¿Б. (14)
V ©* ' Б
Этот функционал допустимо рассматривать на множестве непрерывно дифференцируемых в области V векторных функций я(М), М Е V, удовлетворяющих условию (13) в качестве дополнительного и интегрируемых по поверхности Б. Требования к функциям 0(М) (М Е V) остаются такими же, как и к допустимым для функционала (8).
Альтернативные функционалы (8) и (14) составляют двойственную вариационную формулировку модели теплового взрыва. Для них справедливо неравенство I[0, я] ^ 3[0], переходящее в равенство при выполнении условия стационарности (9). Это неравенство дает возможность при количественном анализе модели оценивать наибольшую возможную погрешность, возникающую в случае приближенного представления искомого распределения температуры в рассматриваемой области V.
3. Пример применения вариационного подхода к анализу модели
Рассмотрим особенности применения вариационного подхода на примере известной задачи о тепловом взрыве в бесконечно длинном круговом цилиндре радиусом г0 [1], для которой в упрощенной постановке существует точное аналитическое решение [12]
9(Г) = ^7 = 2.^, В = 4 — 1 | - 1)2 — 1. (.5)
где г = т/т0 — безразмерная радиальная координата и вг = (цУут^ДАТо)) ехр(—7). Это
решение удовлетворяет граничному условию 9(1) = 0 на поверхности цилиндра, условию
бв(г)
= 0 осевой симметрии и уравнению
г=0
dr
1 d-( -d-) + вг exp(0) = O, (16)
r dr V dr J
которое следует из уравнения (4), если в этом уравнении принять L = r0 ив знаменателе показателя экспоненты пренебречь слагаемым по сравнению с единицей.
Согласно второму равенству (15), параметр В принимает действительные значения при в ^ 2. Значение в* = 2 является критическим для сплошного кругового цилиндра. Ему соответствует значение В * = 1 и предельно возможное установившееся распределение в цилиндре безразмерной температуры в*(г) = 21п^2/(1 + г2)) при наибольшем значении безразмерной температуры в*(0) = 21п2 « 1,3863 на его оси. Значению в* = 2 отвечает критический радиус г0 = ^2ЛТ0 ехр(7)/(75^) цилиндра. Установившееся распределение температуры в цилиндре возможно лишь при выполнении неравенства г0 ^ г*.
Связь между параметром вг и безразмерной температурой в(0) на оси цилиндра, определяемая формулами (15), представлена на рис. 1 графиком зависимости этого параметра от в(0) (сплошная кривая). Видно, что фиксированному значению вг < в* = 2 соответствуют два значения в(0), одно меньше предельно возможной безразмерной температуры в*(0), а другое больше ее. Ясно, что из двух распределений температуры в цилиндре, формально следующих из решения задачи, реализуемо лишь одно, для которого в(0) ^ в*(0).
Рис. 1. Связь параметра @г и безразмерной температуры на оси цилиндра
Функция в(г), определяемая формулами (15), является допустимой для функционала (8), который в случае одномерной задачи примет вид
1 (1 ( йв(г)\2
3 [в]
2 4 йг
+ вг( 1 — ехр(в(г))) 1г йг.
(17)
После подстановки указанной функции получим
3 = _
О >к —
вг В
4В
1+В
+ 41п(1 + В).
Для функционала (14) (записанного также для одномерной задачи при в(1) = 0)
I[в, дг] = — / (М^ — вг + вг(1 — в(г)) ехр(в(г))) г йг
(18)
2
достаточно рассматривать лишь радиальную составляющую qr (г) векторной функции я, удовлетворяющую равенствам
1 d(fqr (г))
ßr exp(0(r)), qr(0) = 0,
(19)
r dr
которые следуют из условия (13) как дополнительного для этого функционала. Отсюда с учетом первого равенства (15) получим
ßr (1 + B)2r
qr (r) = T TTBW
и после подстановки допустимых функций в функционал (18) запишем
I- = -А(1 + 2в - er^ + (ЦВ)>^ - в) ШС1 + B)).
Поскольку функции 0(r) и qr (r), использованные при записи формул для функционалов J- и 1-, являются допустимыми при любых значениях параметра В, определяющего безразмерную температуру на оси цилиндра, можно рассмотреть зависимости этих функционалов от 0(0). На рис. 2 такие зависимости для трех значений ßr, равных 16/9, 1,9 и 2, представлены сплошными кривыми и пунктирными линиями соответственно для функционалов J и .
Рис. 2. Графики зависимостей функционалов 3* и I* от безразмерной температуры на оси цилиндра при различных значениях параметра ¡Зг
При вг = 16/9 оба функционала имеют совпадающие стационарные точки с совпадающими стационарными значениями при В' = 0,5 и В" = 2, которым, согласно формулам (15), соответствуют значения 0'(О) = 21п(1 + В') и 0,811 и 0"(О) = 21п(1 + В") ъ 2,197. Сравнение полученных результатов с предельно возможным значением 0* (0) ъ 1,3863 приводит к выводу, что стационарная точка, в которой оба функционала достигают максимума, соответствует нереализуемой ситуации.
При значении Д. = 1,9 происходит сближение стационарных точек функционалов по сравнению с предыдущим случаем (см. рис. 2), а при критическом значении в* = 2 эти точки объединяются в общую точку перегиба на графиках зависимостей функционалов от параметра В, абсцисса которой равна указанному выше значению 9*(0), определяемому значением В* = 1. Таким образом, результаты анализа зависимостей функционалов от параметра В полностью согласуются с точным решением данной задачи.
Представленная выше вариационная формулировка модели дает возможность построить приближенное решение и сравнить его с точным решением. Искомое распределение безразмерной температуры в поперечном сечении кругового цилиндра аппроксимируем квадратичной функцией 9'(г) = — г2), допустимой для функционала (17), который после подстановки этой функции принимает вид
Л = 2 (—2 + А-( 1 + — — ехр(А)). (20)
Параметр имеющий смысл искомой безразмерной температуры на оси цилиндра, находим из необходимого условия = 0 существования стационарной точки функционала ^.
ои\
Это условие приводит к трансцендентному уравнению относительно -1, представимому в виде
2—3
1 + (-1 — 1)ехр(-1)"
По этой формуле на рис. 1 построен график (штриховая линия), наибольшая ордината вГ ~ 2,0939 которого соответствует абсциссе — = 9'(0) ~ 1,3608. Из сопоставления с точным решением задачи следует, что предельное значение параметра вг оказалось выше и отличается менее, чем на 5 %, а предельное значение безразмерной температуры на оси цилиндра ниже примерно на 2 %.
Для сравнения используем в качестве допустимой тригонометрическую функцию 9" (г) =
пг
= —2 сое — и после ее подстановки в функционал (17) получим
22
(4 + n2)D вг п I' (^ пг\ , J = ^-+ ^ - вг J exp D2 cos r dr. (21)
0
Необходимое условие = 0 существования стационарной точки функционала приводу
дит к равенству
4 + п
4 + п2 / } ( nr \ nr \ 1 вг = D2\J exp( D cos J cos rdrl ,
^ 0 '
которое использовано для построения на рис. 1 графика (штрихпунктирной линией), имеющий наибольшую ординату в" ~ 2,0204 при значении абсциссы D ~ 1,3937. Описание распределения безразмерной температуры функцией 0"(f) приводит к результатам, более близким к точному решению по сравнению с использованием функции 0'(f): наибольшая
безразмерная температура ^ = в"(0) на осп цилиндра превышает точное значение примерно на половину процента, а предельное значение параметра в" больше в * примерно на один процент.
Взаимное расположение при различных значениях параметра вг графиков зависимости различных функционалов от безразмерной температуры на оси цилиндра иллюстрирует рис. 3. При значениях вг = 16/9, вг = 1,9 и в* = 2,0 с целью сравнения приведены зависимости функционалов 3* и I* (как и на рис. 2 сплошными кривыми и пунктирными линиями соответственно), а также зависимости функционалов 3\ (штриховые линии) и 32 (штрихпунктирные линии). Видно, что при всех указанных значениях параметра вг графики функционала 32 в окрестности стационарных точек функционалов 3* и I* расположены ближе к стационарным значениям этих функционалов по сравнению с графиками функционала 3\, что позволяет считать, что тригонометрическая функция 0"(г) лучше аппроксимирует точное решение, определяемое формулами (15), по сравнению с квадратичной функцией 0'(г).
/*, л, 32
0,4 0
-0,4 -0
-1,2 -1,6
-2,0 -2,4 -2,8
--
/ /
/
/ /
/ N
/ /
/
/ /
/ / N Ч
/' ¿/у / \ *
ч * \\
...
ч \ • 1
\
Ч— п 1\
"О. //Л
Рг 9
(Зг= 1,9
#=2,0
0 0,4 0,8 1,2 1,6 2,0 0(0), Д, £>2
Рис. 3. Графики зависимостей функционалов от безразмерной температуры на оси цилиндра при различных значениях параметра вг
Сплошная кривая и штриховая линия, отмеченные квадратами, характеризуют зависимости от безразмерной температуры на оси цилиндра функционалов 3* и 3\ при значении в'г = 2,0939. Следует отметить, что в этом случае функционал 3* в силу неравенства в'г > в * не имеет стационарной точки в отличие от функционала 3\. Сплошная кривая и штрихпунк-тирная линия, отмеченные кружками, соответствуют зависимостям функционалов 3* и 32 при значении в'г = 2,0204. Достаточно близкое расположение графиков этих функционалов связано с близостью значений в* и в" по сравнению со значениями в* и в'г, но и в этом случае функционал 3* не имеет стационарной точки, поскольку в" > в*.
4. Анализ двухпараметрической модели теплового взрыва
Рассмотрим двухиараметрическую модель, включающую дифференциальное уравнение (4), ирименительно к иластине, одна из иоверхностей которой идеально теилоизолирована, а на другой задана темиература Т0. Для такой иластины функционал (8) иримет вид
3 [9]
1 ( ¿9(г)\2
2\ ¿г
т
+ вн ехр
9
1 + 7 9
¿г,
(22)
где г = г/к — безразмерная координата, отсчитываемая ио нормали от теилоизолированной иоверхности иластины толщиной к, а вн = к27/(АТ0)) ехр(—7).
В качестве функций, доиустимых для функционала (22), иримем = С\(1 — г2) и 92(г) = С2 сов(пг/2). Отметим, что эти функции удовлетворяют не только доиолнитель-ному для этого функционала условию 9(1) = 0, но и естественному для него граничному
условию
дВ(г)
йг
¿=0
0 на идеально теилоизолированной иоверхности иластины. В случае
квадратичной функции из условия
0.1 [01] 0С\
вн = 4 4 / (1 — г2)ех^ С1
\ п
0 стационарности функционала (22) следует
1- г2
1 + 7С1(1 — г2)
1
¿г
(23)
д.1 [02]
дС2
0 стационарности этого же функционала ири исиользовании тригономе-
а условие
трической функции ириводит к равенству
2 / 1
в'1 = у СЛ /ехр( С2
\ п
сов(пг/2)
^ , , 1 соМ — ) ¿г 1+ 7С2 сов(пг/2У V 2
1
(24)
Интегралы в формулах (23) и (24) (как и внутренний интеграл в иравой части соотношения (22)) не удается выразить в элементарных функциях и ириходится ирибегать к численному интегрированию. На рис. 4 сплошными кривыми ириведены зависимости вн от С1, а
1,0
0,6
0,4
0,2
0
0,4 0,
1,2 1,6 СЪС
Рис. 4. Связь иараметра вн и безразмерной темиературы на идеально тенлоизолированной иоверхности иластины
1
штриховыми — в'' от С2 при различных значениях 7 (кривые без символов соответствуют значению 7 = 0, т.е. однопараметрической модели теплового взрыва, кривые со светлыми и темными кружками — соответственно 7 = 0,1 и 7 = 0,2). Видно, что с увеличением параметра 7 возрастают ордината и абсцисса точки максимума на каждой из представленных кривых. Это означает, что с ростом этого параметра возрастают как предельное значение в* , соответствующее состоянию теплового взрыва, так и соответствующие значения в'(0) и в'' (0) безразмерной температуры на теплоизолированной поверхности пластины. Однако графическое представление указанных зависимостей не позволяет с достаточной точностью установить связь между параметрами в* и 7, определяющую предельное установившееся температурное состояние пластины, предшествующее тепловому взрыву.
Используя при фиксированном значении 7 необходимое условие экстремума правых
частей равенств (23) и (24) в виде ^^ = 0 и ^^ = 0, получаем
дС} ОС2
= /0 - (Г+^СоГ-Ь)(1 - Д-1 ^2)) * = 0.
, }( С2 сов(п;72) \ ( С2 со8(пг/2) ) (пг) ,
^2(^2.7) = / 1 - 7-^ / , N^ ехр(-2 у ', , сое — Ыг = 0.
А 2' ¿V (1+ 7С2 сов(п;г/2))^ *Л1 + 7С2 со8(пг/2)У \2)
Задачу нахождения удовлетворяющих этим равенствам значений С* и С* при заданном значении 7 целесообразно свести к решению (с контролируемой погрешностью) задачи безусловной минимизации функций (С}) = (Т}(С1.7))2 и Т"(С2) = (Т2(С2,7))2, используя один из известных эффективных численных методов [13, 14]. По найденным значениям СИ и С** затем несложно вычислить по формулам (23) и (24) соответствующие предельные значения (в/)* и (в/')*.
В табл. 1 приведены результаты вычисления указанных параметров с четырьмя верными десятичными знаками. Для сравнения в этой таблице представлены пары значений в*(0) и в/И, полученные с такой же точностью на основе известного аналитического представления [1, 12]
0(0)/ 0(0) \-}/2 ^ = / [2 I ^ (25)
Ф) V 0 ' /
Таблица 1
7 С* (А )* си (в")* 0*(0) в** в°н
0 1,1819 0,8923 1,1957 0,8773 1,1860 0,8785 0,880
0,01 1,2071 0,9017 1,2213 0,8867 1,2122 0,8878 0,880
0,02 1,2336 0,9115 1,2482 0,8964 1,2390 0,8975 0,892
0,025 1,2474 0,9165 1,2622 0,9013 1,2520 0,9024 0,904
0,04 1,2912 0,9321 1,3068 0,9167 1,2974 0,9178 0,915
0,05 1,3227 0,9429 1,3388 0,9274 1,3292 0,9284 0,927
0,10 1,5154 1,0035 1,5347 0,9873 1,5244 0,9882 -
точного решения дифференциального уравнения вида (4) для рассматриваемой пластины. Наконец, в последний столбец таблицы помещены значения ßh [1], полученные численным решением в работе [15]. Следует отметить, что при некоторых значениях y разность в — ßh I существенно превышает погрешность вычислений с использованием точного решения (25).
Заключение
Вариационный подход к количественному анализу одно- и двухпараметрической математических моделей теплового взрыва, опирающийся на двойственную вариационную формулировку нелинейной задачи стационарной теплопроводности в твердом теле, позволил при использовании приближенных распределений температуры в длинном круговом цилиндре и в пластине, неограниченной в своей плоскости, получить результаты, согласующиеся с известными в литературе, что дает основание для использования данного подхода к анализу моделей теплового взрыва в твердых телах произвольной формы.
Работа выполнена по гранту МК-6573.2015.8 программы Президента РФ государственной поддержки молодых кандидатов наук, а также в рамках проекта 1712 по государственному заданию № 2014/104 Минобрнауки РФ и государственного задания по проекту № 1.2640.2014.
Список литературы
1. Франк-Каменецкий Д.А. Диффузия и теплопередача в химической кинетике: 3-е изд., испр. и доп. М.: Наука, 1987. 502 с.
2. Зарубин В.С., Кувыркин Г.Н., Савельева И.Ю. Температурное состояние диска униполярного генератора // Инженерно-физический журнал. 2014. Т. 87, № 4. С. 796-801.
3. Зарубин В.С., Котович А.В., Кувыркин Г.Н. Устойчивость температурного состояния диска униполярного генератора//Известия РАН. Энергетика. 2016. № 1. С. 127-133.
4. Физика взрыва / Под ред. Л.П. Орленко: 3-е изд., перераб. В 2 т. Т. 1. М.: ФИЗМАТЛИТ,
2002. 832 с.
5. Зарубин В.С. Моделирование. М.: Издательский центр «Академия», 2013. 336 с.
6. Зарубин В.С., Кувыркин Г.Н. Особенности математического моделирования технических устройств //Математическое моделирование и численные методы. 2014. № 1(1). С. 5-17. DOI: 10.18698/2309-3684-2014-1-517
7. Самарский А.А., Галактионов В.А., Курдюмов С.П., Михайлов А.П. Режимы с обострением в задачах для квазилинейных параболических уравнений. М.: Наука, 1987. 480 с.
8. Зарубин В.С., Кувыркин Г.Н. Математическое моделирование термомеханических процессов при интенсивном тепловом воздействии // Теплофизика высоких температур.
2003. Т. 41. №2. С. 300-309.
9. Зарубин B.C., Кувыркин Г.Н. Математические модели механики и электродинамики сплошной среды. М.: Изд-во МГТУ им. Н.Э. Баумана, 2008. 512 с.
10. Власова Е.А., Зарубин B.C., Кувыркин Г.Н. Приближенные методы математической физики: 2-е изд., стереотип. М.: Изд-во МГТУ им. Н.Э. Баумана, 2004. 700 с.
11. Глаголев К.В., Морозов А.Н. Физическая термодинамика: 2-е изд., испр. М.: Изд-во МГТУ им. Н.Э. Баумана, 2007. 272 с.
12. Зарубин B.C. Инженерные методы решения задач теплопроводности. М.: Энергоатом-издат, 1983. 328 с.
13. Аттетков А.В., Зарубин B.C., Канатников А.Н. Введение в методы оптимизации. М.: Финансы и статистика, ИНФРА-М, 2008. 272 с.
14. Аттетков А.В., Зарубин B.C., Канатников А.Н. Методы оптимизации. М.: Издательский центр РИОР, 2012. 270 с.
15. Parks J.R. Criticality Criteria for Various Configurations of a Self-Heating Chemical as Functions of Activation Energy and Temperature of Assembly // J. Chem. Phys. 1961. Vol. 34, no. 1. P. 46-50. DOI: 10.1016/0022-247X(81)90213-4
Mathematics i Mathematical Modelling
Electronic journal of the Bauman MSTU
Mathematics and Mathematical Modelling of the Bauman MSTU, 2016, no. 5, pp. 29-45.
DOI: 10.7463/mathm.0516.0847523
Received: 23.08.2016
© Bauman Moscow State Technical University
http://mathmjournal.ru
Variational Approach to the Thermal Explosion Model Analysis in a Solid
Zarubin V. S.1*, Kuvyrkin G. N.1, Savel'eva I. Yu.1 * [email protected]
1 Bauman Moscow State Technical University, Russia
Keywords: mathematical model of thermal explosion, the dual variational formulation of nonlinear problem, alternative functional
A temperature state of the solid body may depend both on the conditions of heat exchange with external environment surrounding its surface and on the energy release within the body volume, caused, for example, by the processes in nuclear reactor elements or exothermic chemical reactions, absorption of penetrating radiation energy or transformation of a part of the electrical power into heat with flowing electric current (so-called Joule heat).
If with growing temperature the intensity of bulk power density increases, a limited steady temperature state can emerge at which heat extracted to the body surface and released within its volume reaches maximum. Thus, small increments of temperature lead to an increase of heat release, which can not be extracted to the body surface by conduction without further temperature increase. As a result, the steady temperature distribution in the body becomes impossible that determines the state of the thermal explosion, so named due to the fact that in this case the appropriate mathematical model predicts an unlimited temperature increase.
A lot of published papers and monographs concerning the study of the combustion and explosion processes in a stationary medium analyse the thermal explosion state. The most famous papers consider a mathematical model to describe a temperature distribution in the case when heat release is because of exothermic chemical reactions the rate of which increases with temperature growth. The dependence of the chemical reaction rate on temperature is usually described by the exponential Arrhenius law, which makes it necessary to consider an essentially nonlinear mathematical model containing differential equation, which includes the term, nonlinearly rising with increasing temperature. Even with simplifying assumptions, this model allows an exact closed form solution only in the case of one-dimensional temperature distributions in the two areas of the canonical form: in the plate, infinite in its plane, and in the circular cylinder unlimited in length.
An approximate numerical solution of the differential equation that is included in a nonlinear mathematical model of the thermal explosion enables us to obtain quantitative estimates of combi-
nation of determining parameters at which the limit state occurs in areas of not only canonical form. A capability to study of the thermal explosion state can be extended in the context of development of mathematical modeling methods, including methods of model analysis to describe the thermal state of solids.
To analyse a mathematical model of the thermal explosion in a homogeneous solid the paper uses a variational approach based on the dual variational formulation of the appropriate nonlinear stationary problem of heat conduction in such a body. This formulation contains two alternative functional reaching the matching values in their stationary points corresponding to the true temperature distribution. This functional feature allows you to not only get an approximate quantitative estimate of the combination of parameters that determine the thermal explosion state, but also to find the greatest possible error in such estimation.
References
1. Frank-Kamenetskiy D.A. Diffuziya i teploperedacha v khimicheskoi kinetike [Diffusion and Heat Transfer in Chemical Kinetics]. Moscow, Nauka Publ., 1987. 502 p. (In Russian).
2. Zarubin V.S., Kuvyrkin G.N., Savelyeva I.Yu. Temperature state of a unipolar generator disk. Inzhenerno-fizicheskii zhurnal, 2014, vol. 87, no. 4. pp. 796-801. (English version of journal: Journal of Engineering Physics and Thermophysics, 2014, vol. 87, no. 4, pp. 820-826. DOI: 10.1007/s10891-014-1077-2).
3. Zarubin V.S., Kotovich A.V., Kuvyrkin G.N. Temperature condition stability of the disc unipolar generator. Izvestiya RAN. Energetika = Proceedings of the Russian Academy of Sciences. Energy, 2016, no. 1, pp. 127-133. (in Russian).
4. Orlenko L.P., ed. Fizika vzryva. T. 1 [Explosion Physics. Vol. 1]. Moscow, Fizmatlit Publ., 2002. 832 p. (in Russian).
5. Zarubin V.S.Modelirovanie [Modeling]. Moscow, "Academia" Publishing Center, 2013. 336 p. (in Russian).
6. Zarubin V.S., Kuvyrkin G.N. Special features of mathematical modeling of technical instruments. Matemeticheskoe modelirovanie i chislennye metody = Mathematical Modeling and Computational Methods, 2014, no. 1(1), pp. 5-17. DOI: 10.18698/2309-3684-2014-1-517 (in Russian).
7. Samarskii A.A., Galaktionov V.A., Kurdyumov S.P., Mikhailov A.P. Rezhimy s obostreniem v zadachakh dlya kvazilineinykh parabolicheskikh uravnenii [Modes with peaking in problems for quasi-linear parabolic equations]. Moscow, Nauka Publ., 1987. 480 p. (in Russian).
8. Zarubin V.S., Kuvyrkin G.N. Mathematical Modeling of Thermomechanical Processes under Intense Thermal Effect. Teplofizika vysokikh temperatur, 2003, vol.41, no. 2, pp. 300-309. (English version of journal: High Temperature, 2003, vol.41, no. 2, pp. 257-265. DOI: 10.1023/A:1023390021091).
9. Zarubin V.S., Kuvyrkin G.N. Matematicheskie modeli mekhaniki i eleknhjdinamiki sploshnoi sredy [Mathematical models of continuous media mechanics and electrodynamics]. Moscow, Bauman MSTU Publ., 2008. 512 p. (in Russian).
10. Vlasova E.A., Zarubin V.S., Kuvyrkin G.N. Priblizhennye metody matematicheskoi fiziki [Approximate methods of mathematical physics]. Moscow, Bauman MSTU Publ., 2004. 700 p. (in Russian).
11. Glagolev K.V., Morozov A.N. Fizicheskaya termodinamika [Physical thermodynamics]. Moscow, Bauman MSTU Publ., 2007. 272 p. (in Russian).
12. Zarubin V.S. Inzhenernye metody resheniyazadach teploprovodnosti [Engineering methods for solving heat conduction problems]. Moscow, EnergoatomisdatPubl., 1983. 328 p. (in Russian).
13. Attetkov A.V., Zarubin V.S., Kanatnikov A.N. Vvedenie v metody optimizatsii [Introduction to Optimization Methods]. Moscow, Finance and Statistics Publ., SIC INFRA-M Publ., 2008. 272 p. (in Russian).
14. Attetkov A.V., Zarubin V.S., Kanatnikov A.N. Metody optimizatsii [Optimization methods]. Moscow, IC RIORPubl., 2012. 270 p. (in Russian).
15. Parks J.R. Criticality Criteria for Various Configurations of a Self-Heating Chemical as Functions of Activation Energy and Temperature of Assembly. J.Chem. Phys., 1961, vol. 34, no. 1, pp. 46-50. DOI: 10.1063/1.1731612
