Математическое моделирование высокотемпературных нелинейных процессов теплопереноса

Восточно-Европейский журнал передовым технологий ЕЗМ 1729-3774

I--Запропоновано наближен аналтико-числовi методи математичного моделю-вання високотемпературних нелтшних процеыв теплоперенесення. Розглянуто чотири моделi теплового режиму гiр-ничих масивiв при пожежах у шахтних виробках - першi та третi крайовi задачi нелтшног теплопровiдностi при сталих та змтних граничних умовах

Ключовi слова: математичне моделю-вання, високотемпературн процеси, гiрничi масиви, аналiтико-числовi методи □-□

Предложены приближенные аналити-ко-численные методы математического моделирования высокотемпературных нелинейных процессов теплопереноса. Рассмотрены четыре модели теплового режима горных массивов при пожарах в шахтных выработках - первые и третьи краевые задачи нелинейной теплопроводности при постоянных и переменных граничных условиях

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ВЫСОКОТЕМПЕРАТУРНЫХ НЕЛИНЕЙНЫХ ПРОЦЕССОВ ТЕПЛОПЕРЕНОСА

А. П. Слесарен ко

Доктор физико-математических наук, профессор, лауреат Государственной премии Украины, ведущий научный сотрудник Институт проблем машиностроения им. Подгорного НАН Украины ул. Пожарского, 2/10, г. Харьков, Украина, 61046 Контактный тел.: (057) 265-51-89, 096-386-30-22

И.Р. Венгеров Кандидат физико-математических наук, старший научный сотрудник Институт физики горных процессов НАН Украины ул. Р. Люксембург, 72, г. Донецк, Украина, 83114 Контактный тел.: (062) 337-80-46, 050-998-34-67

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

Разработка методов профилактики, тушения и локализации последних невозможна без математического моделирования тепловых режимов горных массивов и выработок [1]. Соответствующие математические модели формулируются как краевые задачи нелинейной теплопроводности, решаемые обычно численными методами, требующими анализа сходимости вычислительных алгоритмов в каждом конкретном случае и зачастую малопригодными для проведения оперативных инженерных расчетов.

Более перспективны, на наш взгляд, приближенные аналитико-численные методы решения нелинейных краевых задач, базирующиеся на структурах решений (региональных и глобальных), что было показано в работах А.П. Слесаренко [2].

В настоящей работе такие приближенные анали-тико-численные методы используются для решения ряда задач высокотемпературного нагрева горного массива (пожарными газами с температурой порядка 103 К).

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

Анализ литературных источников показал [2], что нелинейные краевые задачи теплопроводности твердых тел решаются, наряду с численными методами, аналитическими: использованием автомодельных переменных, подстановок, полных и частичных линеаризаций и т.д. [2-5]. Эти методы, как правило, применяются к задачам с краевыми условиями простейшего вида и с ограниченным набором функций температуры - теплофизических параметров систем.

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

В настоящей работе мы ограничиваемся одномерными моделями теплопереноса во внешних областях □+1& = {х е(0, ~)} .

Целью работы является демонстрация предлагаемых методов аппроксимации и стратификации на четырех, практически важных, моделях подземных пожаров. Используем метод экспоненциальной аппроксимации температурных полей в ограниченных оценками локализации полей конечных областях [3] и методы стратификации двух видов: а) стратификации по координате х (разбиение, специальным образом, конечной области на конечное число слоев); б) хроностратификации по переменной t (разбиение заданного конечного интервала времени процесса ^ на ряд хронослоев).

T (t) = Tr (0) +[Tr (ts) - Tr (0)](t/ts)"T, Tr (0) = T<-),T(ts) = T«.

Формулой (1) описывается и случай «термического удара»: при t>0 и nT = 0, Tr(t) = T+& = const. Он является мажорантным, представляя «правую» (верхнюю) границу оценочной вилки для величины зоны локализации поля.

Ответ на второй вопрос следует из того факта, что максимальная температура в горном массиве не может превышать TJ(+) при отсутствии в нём тепловых источников. При наличии их, что характеризуется заданием функции их плотности F(n,t)< F при nsQ+m) и t e(0,ts) , будем иметь:

T(m)(n,t) <^ts = TF+), Cmln = minC (T), T e[T,T]. (2) Здесь:

Задача №1. Предпосылки моделирования. Формулировка последних заключается в ответах на следующие вопросы: 1) какой процесс моделируется - нагревание или охлаждение области О"&; 2) каковы

минимальная (т) и максимальная (т) температуры в области и на её границах в моделируемый период t е[0,^],^ <~; 3) как с изменением температуры изменяются теплофизпараметры; 4) как оценивается зона локализации температурного поля; 5) как находятся аппроксимации температурного поля.

Первый вопрос актуален потому, что в отличие от процессов линейной теплопроводности, при нелинейной теплопроводности нет эквивалентности задач на нагревание и охлаждение из-за явлений термического гистерезиса горных пород. Вид функций С(Т) и ^(Т) при нагревании до высоких температур и остывании затем до умеренных, различен [2]. В частности, при моделировании подземных пожаров, необходимо рассмотрение двух подзадач: первая должна обеспечить прогноз температур горного массива при возникновении, развитии и стабилизации пожара, характеризуемого динамикой температуры нагревающей массив среды (пожарных газов и др.) - е(0,^). Вторая

подзадача решает ту же задачу в период затухания пожара, его тушения и последующей вентиляции выработки ^ е(^,»)). Мы рассматриваем только первую подзадачу, полагая, что горный массив нагревается. Функция Тг(^ при этом должна быть представлена в нормализованном (бимонотонном) виде:

T = minT(m)(n,t) = minф(п), Пе[П-,то),

T = max {TIi+),TF+)}.

Следовательно, в любой точке области О+т) и для любого момента времени t е(0,^), имеет место вилка

T < T(m)(n,t) < T, AT = max AT(m) = T - T.

Здесь ДТ- температурный диапазон модели (ширина вилки).

Для ответа на третий вопрос достаточно нормализовать (представить - бимонотонными аппроксимациями) первичные зависимости с (Т) и ^(Т):

Г й n Т - Т |

c (Т )-c(T)_ 1 AT J

c (T) = c(-),c (T) = c(+) = Kcc(-),

(T) = c(-),c (T)

£ (t)=£(t)+ £(t)-£(t)

&T - T 4 n£

£(T)=£(-), £(t) = £(+) = K££(-).

^Масштабные множители Кс и К, определяемые ДТ, могут принимать значения меньше или больше единицы. Для коэффициента теплопроводности а=а(Т) получаем:

a (Т^У^Ь a(-)

a(+) = a (T) = K£ a(-) =

^ K£ a(-) = £(+) K a = c(+),

а (Т) = а(-) + (а(+) - а(-)) 6"- ,п. = " (п , ,Пс,Кс,Кх).

Для ответа на четвертый вопрос заметим следующее. Как показано в [3], при линейной теплопроводности в отсутствие источников тепла, для 8+ (^ - зон локализации полей в областях О+т) можно использо+

вать выражение: S+(t) = 4л/°Е (a = const). В нелинейных средах, где a = a(Т), как показано в [6], приведенная формула справедлива при учете в ней переменности a . Поэтому возможен «вилочный» подход:

s(t)<s+(t)<s(t), s(t)=n-+Wat, s(t)=n-+Wit, (9)

a = a(-), a = a(+), n-={0, m = i4 t e(0,ts). (10) [r0, m = 2,3, 4 &

На основе (9) осуществляем переход от сингулярной области к конечной: Q+m) ^Q(m) ={ne(n-, S(t)} .

Область qh может быть далее преобразована, процедурой стратификации, в слоистую систему (двух-, трех- или многослойную).

Пятый вопрос корректен только для случая однородного начального условия (ф(п) = Tn = const) при отсутствии источников тепла в Q+m). Для этого случая записываем вилку, граничные значения которой - ранее уже использованные экспоненциальные аппроксимации [3]:

T(m)(n,t)< T(m)(n,t)< T(m)(n,t), neQ(m), t e(0,ts), (11)

T(m)(n,t) = Tn +(Tr(t)-Tn )exp T(m)(n,t) = Tn +(Tr(t)-Tn )exp

_n-n_ S(t)-n_n-n_ S(t)-n, (12)

Из последних формул для полуширины вилки R (n,t) получаем:

R (n,t) = 0,5 [T(m)(n,t> - T(m)(n,t>]

= AT (t){ exp

_n-n-8(t)-n-exp

_n-n-s(t)-n(14)

ДТ^) = тг(t)- Тп <ДТ,

ДT(ts ) = Tг(ts)-Т" = ТГ+)- Т" = ДТ.

Задача №2. Первая краевая задача с постоянной граничной температурой. Для т=1. Р(1) = 3,551 [3] и формулы (12,13) принимают вид:

Т (x,t)= Т„ + (Т]Г+) - Т„) exp

x er0,^Vat J, t > 0,

T (x,t)= T„ + (t<+)- T„) exp

x e[0,Wat I, t > 0.

Положим в (14): т=1, х/-\\^ = г|, а/а = Ка, тогда: Я (x,t) = Я(п,Ка) = 0,5 [ехр (-3,551п) - ехр (-3,551К0,5п)]. (17)

Численные расчеты по (17) при п £ (0,1],Ка е [1,2; 1,5],

Дп = 0,01, ДКа = 0,1 дают значения Я(п,Ка), приведенные в табл. 1

Таблица 1

Относительная полуширина аппроксимационной вилки

Из таблицы следует:

Ка полуширина вилки возрастает (и, как показали расчеты, при Ка > 1,5 достигает значений > 5%). При Ка < 1,5 решение задачи №2 исчерпывается использованием одного из двух практически эквивалентных выражений

T(1)(x,t) = Tn + 0,5(Т^)-Tn) х х [exp (-3,551К°/п) + exp (-3,551n)],

T(1)(x,t) = T„ +(T<+)-T„ )exp

s(t), s(t)=4Vat.

В случае, если условие задачи №2 таковы, что Ка = а/а > 1,5 применим стратификацию области О(1) = {х е(0,8(t))} на две подобласти (слоя) переменной ширины:

оГи^М}, О81) = {х е(0,81 М)},

о812) = {х е(81 (t), 82 М)}.

Здесь параметр 81 (t)- расстояние, пройденное после начала процесса изотермой Т = Тг0 < Т^+) за время t , а параметр 82 (t)- ширина зоны локализации температурного поля к этому моменту времени, т.е. Т(81 (t),t) = Тг0, Т(82^)^) = Тп . Оба эти параметра могут быть представлены в виде: 8! (^= Zi^/at(i = 1,2), Z2 = 4, а Z1 - требует определения.

Поскольку а = а(Т), воспользуемся оценками (т.е., построим вилку):

§1 ^<§1 (^ (t), §2 М<§2 (^ (t), 81 (t) = §1 (t) = Z^^/ft, 82 (t) = 82 (t) = 4Л/1~, а = а (Т0 ),а = а (Т?),а = а (т(+)), К = Ка1 -Ка2,Ка1 = аД Ка2 = а/а.

Поскольку Ка1 < 1,5, Ка2 < 1,5, то Ка < 2,25.. Если Ка > 2,25 , то двухслойную стратификацию надо заменить на трёхслойную, которая будет справедлива при Ка = Ка1 ■ Ка2 ■ Ка3 < (1.5)3 « 3.38 . Если же Ка > 3,38 необходима четырёхслойная стратификация и т.д.

Значения а = а(Тп) и а = а(Т^) находится из (10), а значение а = а(Тг0) определяется (а по нему и значение «промежуточной» температуры Тг0 ) из условия:

К,=а=аШ=1,5.

а2 а" а(Т)

Для оценки параметра 71 воспользуемся решением линейной задачи о термическом ударе [3], принимающем в задаче №2 вид:

Сопряжение температурных полей на общей границе слоёв Щ и Ц(12 будет автоматическим, но условие равенства производных, вытекающее из физики процесса

ЭТ2(1)р

_ 2,ср

-51,срМ

, t > 0,

может оказаться не соблюдающимся. Чтобы обеспечить выполнение (28), воспользуемся для поля в Щ полиномиальной (квадратичной) аппроксимацией:

Т« (x,t) = АИ + B(t)х + C(t)х2, х е [0,8^ , t > 0. (29)

Коэффициенты А(t),B(t),C(t) определяем из (28) и условий

Т(1)р(х^)| = Т<+), Т(1)р (х^)| = Тг0.

В итоге получаем:

Т (x,t)= Т„ + (Т+)- Т„) erfc I

Л 27!:)&

Если в (22) положить Т(х^) = Тг0,х = то получим:

Т -Т ег 0 = Т+0-^=erfc1 г 0 т(+) — Т I 2

Из (23), по известным Тп,ТГГ+),Тг0, параметр 71 легко находится, так что далее считаем его известным. Теперь параметры §1 (t) и §2 (t) в (20) определим как среднеарифметические граничных значений вилок:

А= Тг(+), =

(2Т1Г+) - ТГ0 - Тп )-У(К)(ТГ0 - Тп)

^(К) К2 + 2,51К1 г(t) (ТГ+)- Тп )-^(К)(Тг0- Тп) ^(К)" К2 - К1 , 2 (t)

т« М=Т«(2Т«+)-ТГ0 - Тп )-^(К)(ТГ0 - Тп)

(Т+)-Тп )-У(К)(ТГ0 - Тп)&

/ \\2 х

, хе^К^)].

§1,ср(t) = 0,5(§1 М + 8(t)) = 7(^+ ТЕ) = ^б^, (24) Т.о., функции Т« (x,t) (32) и Т« (x,t) (27) пред/- = \\ 1 + к0,5 - §2и = §2,ср(t) = 0,5(§2(t) + s2(t))= ^^к12 8(t) = k2s(t), (25)

Температурное поле в слое Ц(1) можно представить в виде экспоненциальной аппроксимации для Т(1)(х,^ (15), осуществив в ней замены

Т« ^ТГ0, х^х-81 М, «(^(К-К1 )8(^:

ставляют приближенное решение задачи №2 - о термическом ударе в области Ц(1) с теплофизпара-метрами, зависящими от температуры. Это решение, при t = 0, когда 81ср(t) = K18(t) = 0 и Т-р(x,t) не существует, за счет Т^ (x,t) , как это следует из (27), обращается в Тп , т.е. удовлетворяет начальному условию. При t , когда напротив, решение Т^р (х^) распространяется на область х е(0,, из (32) следует, что НтТ1(1с)р (x,t) = тГ+) , как это вытекает из постановки задачи. Общий вид решения:

Т21 = Тп +(ТГ0 - Тп)ехр

х е[К18,К28].

х - K18(t)

(К2 - К

В слое Ц81 можно воспользоваться аналогичной

аппроксимацией, пот Т1(1с)р (x,t) граничных условий:

Т« (0^) = Т(+),Т« (81,ср(t),t)= ТГ0.

T(1)(x,t) =

[т« (х4 х е(0,81,ср (t)), t > 0,

[т^ (x,t), х е(81,ср (t), 82,ср (t)), t > 0.

выполнения для поля

Задача №3. Первая краевая задача с переменной граничной температурой. В этом случае функция Тг (^ бимонотонновозрастаетот Тг (0) = Тп до Тг (^) = Т]Г+) согласно (1), а решение соответствующей линейной задачи не позволяет найти, в достаточной простой форме, закон движения изотермы Т (х^) = Тг0. В отсутствие этого закона построить аппроксимации полей Т1 (х^)< Тг0 и Т2 (x,t)> Тг0 нельзя. Это обстоятельство

требует отказа от пространственной стратификации, использованной в задаче №2, и осуществления хроно-стратификации - разбиения всего интервала времени задачи tе(0,^] на два: tе(0,tг) и tе^г,ts]. Решение будем искать в виде:

T(2)(x,t) =

T|(2)(x,t), x e(0,S<2P (t)), t e[0,t ], T22)(x,t), x e (0, §21 (t),), t e[t,ts ].

Пусть, как и в предыдущем случае, температуре Тг0 = Тг (tг) будетсоответствоватьзначение а = а (Тг0) такое, что а/а = Ка1 = 1,5. Тогда по аналогии с задачей №2, можем записать:

T<2)(x,t)= Т„ +(Tr (t)-Т„ )exp

x e(0,S<2p (t)),t e(0,tr). Здесь:

S|2p (t) = 0,5(S12) (t) + (12) (t)) =

(12)p (t)

=0,5(^vat+4^ )=i ^ §(2) (t)=s(t)=4Vat.

s(2)(t),

T22)(x,t)= T„ +(Tr(t)-T„ )exp x£(0,82^ (t)),t e(tr,ts). Здесь:

= 0,5(4^+WE)= ^l05 K!°15s(t).

Т.о., решение задачи №3 имеет вид (34), где Т1(2) (х^) и Т2(2) (x,t) определены, соответственно, (35) и (37).

Параметр 8^) = 4Vat известен,

т.к. а = а (Тп), 8(2сР (t) = 1,118^) (т.к. Ка1 = 1,5),

б® (t) = 0,612(1 + K025)8(t).

Параметр Ka2 находится из соотношения K°2 = a /a, где a = a(t^+)),° = a(Тг0). Величина T0 следует из условия Ka1 = a /a = 1,5 (a = a (Тг0)). Значение t = tr находится из (1):

Т(+) - Т V 1г Ln/

Задача №4. Третья краевая задача с постоянной температурой греющей среды. Здесь предполагается, что в граничном условии третьего рода на границе х = 0 области Ц(1) в момент времени t = 0 скачком устанавливается и далее сохраняется постоянная температура границей среды

-Т(+): г

£(Т(3))

= а0 [T(3)(x,t)- тГ+) ]x=0, t e(0,ts). (40)

Если Тг0 = Tr(ts) = Т^+), т.е. K° = a/a< 1,5 для всего периода процесса, то (35) будет приближенным решением задачи №3.

Если же Тг(^) = Т^+) > Тг0, Ka2 = a/a < 1,5, то можно находить решения на интервале t е^г0,^]-T2(2)(x,t). Если окажется, что Ka2 > 1,5 , необходимо хронострати-фикацию осуществлять на три временных слоя.

Пусть K°2 < 1.5, т.е. соблюдено требуемое условие достаточности получаемого приближенного решения (т.е. «узости» соответствующей вилки). Тогда решение T^2)(x,t) можно представить в виде:

Здесь: T(3)(x,t)- температурное поле в области Q(1)(T(3)(x,t)< Т^+)), а0 = const - коэффициент конвективного теплообмена, £(т(3))| =£(^(+)(t)), |i(3)(t)- переменная температура на границе области x = 0. Использовать решение задачи №2 нельзя, т.к. нет оценки для движения изотермической поверхности, позволившей осуществить двухслойную стратификацию по координате.

Воспользуемся аналогией с задачей №2, где граничная функция TI^(t) бимонотонно возрастает при t e(0,ts) . В настоящей задаче №4 эту роль играет функция ц(3) (t), которую, в духе метода функций склейки [3] будем считать известной, сводя третью краевую задачу к первой. Осуществляем хроно-стратификацию на два временных слоя: tе(0^г) и t е^г,ts]. Для решений на этих слоях используем, адаптируя к настоящей задаче, выражения (35) и (37):

T<3)(x,t) = Т„ + [^13)(t)-Т„ J exp

x e(0, S(3c)p (t)),t е(0^г),

T<3)(x,t) = T„ + [^23) (t) - T J exp x e(0,8^ (t)),t e(tг ,ts).

(21 (t).

Здесь:

T(3)(x,t) = -(43)

T1(3)(x,t) и(з)(t) = ^^13)(x,t), t e^^^

T23)(x,t/ ^ U = [^3)(x,t), t e^ ,ts ].

Параметры 8(3)р(t)(i = 1,2) формально совпадают с 81 (t), но входящие в пик величины, кроме а = а(Тп), определены иначе: а = а(ц(3) (^)) = а(Тг0), а = а(ц(3) (ts)).

Здесь ц(3)(^) = Тг0 - такое значение, которое обеспечивает выполнение условия а/а = Ка1 = 1,5, а

ц(% ) = ц!3)< Т+>.

Функция Ц(3)^) должна быть определена с помощью граничного условия (40) с учетом (43), т.е. представленного для каждого из двух хронослоёв:

= а0И%)-Т<+)), tе(0,tг), (44)

Х^))^ =а0(ц^)-Т«), tе[tr,ts]. (45)

Подстановка в (44) и (45) решений (41) и (42) даёт:

е,1 м=

*0"1ср

е,1 М = tе(0,tг),

Т(+) -Т 1Г п

е,2 (t) =

, ч ц2%)- Тп г ,

е,2 Т (:)-Tn, t е[tг Л ].

В температурном диапазоне задачи

Т^х^е^Т^]

зависимость Х = Х(Т) в нормализованном виде даётся (7).

В более узких температурных диапазонах ц13)(^е[Тп,Тг0] и ^(^[т^дМ] можно воспользоваться линейными аппроксимациями Х=Х(Т) :

АХ1 = Х - Х, АХ2 = Х - Х, т(+) - Т Т(+) - т

р _ 1г 1п р _ 1г п

я- =-, =Т -Т

Т(+) -Т 1Г 1Г0

АХ1Я1

, е02 =

АХ(Я(

Поскольку, в соответствии с ранее сказанным, (t) = 4,45^ а 823ср(^ = 5,45>/а, выражения (46) и (47) с учетом (50) и (51) приводятся к виду:

е^)+(х,(^0,)е^)-х,(t)=0, i = 1,2,

/ \\ а,

Х (t) = 1,267 (АОЯ ) Х2 (t) = 1,553

(АХ2Я2 )

Отбросив вторые корни уравнений (52) как не имеющие физического смысла, решения последних приводим к виду:

е^ М =

х, (t)+еo

(х,^,)2

Т.о., функции ц,3)(^(, = 1,2) найдены:

цр^) = Тп + (Т+)-Тп )е„ (t).

Поскольку при t ^ 0 ец1 (t) = 0, а при t ^ ец2 (t) ^ 1, найденные функции Ц(3)(^ удовлетворяют условиям начального и конечного состояний

Задача №5. Третья краевая задача с переменной температурой греющей среды. Ход решения повторяет задачу №4, но в итоге формулы (52) несколько видоизменяются, принимая вид:

е^)+(х,)е^)-х,(t)еr(t)=0, , = 1,2, (56)

Х^ф^+р -Х(-))

Ц^ЬТп

х(ц<3)^))=х+(х(+)-х )

ц(3)М-Тп т(+) - т

-■-Г -&-ГЛ

ег (t) = ег(t)е[o,l]

при Тг (0) = Тп и t е(0, .

Последние зависимости приводятся, как легко проверить, к видам:

х(ц<3)м) = аха (еЦ1 (t)+еol), х(|43)(^) = ах2я2 (еЦ2 (t)+е02),

Все другие обозначения соответствуют задаче №4. Решения уравнений (56) соответствуют (54), где под знаком квадратного корня в числителе дроби надо осуществить замену: х, (^ ^ Х, (t)ег (^ . Сравнение формул для в задачах №4 и №5 показывает, что в последнем случае эти функции в задаче №5 меньше таковых в задаче №4 при всех t е(0,, что соответствует физике процесса и подтверждает мажорантный характер решений задач о термическом ударе.

Литература

Abstract

Mathematical modeling of thermal conditions of massifs during underground fires is the theoretical basis for engineering calculations, planning and firefighting.

The article suggests approximate analytical and numerical methods of mathematical modeling of high-temperature heating of a semi-infinite massif taking into account the dependence of its thermalphysic parameters on temperature.

We have examined four mathematical models based on the first and third boundary problems of nonlinear heat transfer with constant and variable temperatures of the heating medium (gas fire).

We have used the method of approximation of unsteady temperature fields, the method of stratification of the decision area, and the method of chronostratification (partitioning of process on the final time intervals).

We have obtained the solutions in forms, which can reduce them "to the number" with a PC; we have estimated the error of the approximate methods, which is less than 5%. The results can be used to improve the reliability, efficiency and quality of engineering calculations in mining industry.