Содержание к диссертации
Введение
Глава 1. Моделирование процесса терморасширепия графита в одномерном случае 13
1.1. Физическая модель процесса терморасширения графита 13
1.2. Математическая модель 14
1.3. Методы решения задачи Стефана 16
1.4. Численное решение одномерной задачи терморасширения графита методом сквозного счета 21
1.5. Модель терморасширения графита в ограниченном объеме 31
Глава 2. Моделирование процесса терморасширения графита в двумерной области 42
2.1. Математическая модель 42
2.2. Метод сквозного счета численного решения двумерной задачи терморасширения графита 49
2.3. Оценка влияния конвективных членов 58
Глава 3. Совместная задача определения поля температур и скоростей двухкомпонентпой среды 62
3.1. Постановка задачи определения скоростей движения графита при терморасширении 62
3.2. Условие согласованности 64
3.3. Задача определения потенциала скоростей. Конечно-разностная аппроксимация 67
3.4. Численное решение задачи определения потенциала скоростей 71
3.5. Решение задачи с подвижной границей с одновременным определением поля скоростей среды 76
3.6. Постановка и решение задачи в осесимметричном случае 81
Заключение 88
Литература
- Методы решения задачи Стефана
- Численное решение одномерной задачи терморасширения графита методом сквозного счета
- Метод сквозного счета численного решения двумерной задачи терморасширения графита
- Задача определения потенциала скоростей. Конечно-разностная аппроксимация
Введение к работе
Изделия из терморасширенного графита (ТРГ) за последнее десятилетие приобрели большую популярность в различных отраслях промышленности, в первую очередь, как уплотнительный материал для оборудования, работающего в условиях высоких температур и агрессивных коррозийных сред. Терморасширенный графит (другие названия: пенографит, терморасщепленный графит) в течение долгого времени считается одним из наиболее устойчивых к разрушению материалов для практически любых задач герметизации систем, работающих с жидкостями. Прокладка и набивка из ТРГ признана пожаробезопасной в условиях применения высоколетучих жидкостей и исключительно высоких температур.
Производством окисленного графита (ОГ) и ТРГ, а также материалов на их основе занимается ряд фирм и организаций. Ведущую роль в нашей стране в данной отрасли занимает ЗАО «Унихимтек» (г. Москва), образовавшееся на базе лаборатории химии и технологии углеродных материалов МГУ им. М.В. Ломоносова [27].
Общая схема получения ТРГ и его последующей переработки представлена на рисунке 1 [48]. В настоящее время подавляющую долю ТРГ перерабатывают в гибкую графитовую фольгу и прессованные изделия [25]. Фольгу получают прокаткой. Сцепление между частицами ТРГ и гибкость фольги обеспечивает разветвленная пенообразная структура ТРГ.
В меньшей степени используют порошки ТРГ. Работа в этом направлении в настоящее время носит поисковый характер. Известно [3,51] применение ТРГ в качестве электропроводных и структурирующих добавок в электроды химических источников тока, литированием и фторированием ТРГ получают активные реагенты химических источников тока, модифицированием металлами регулируют тепло- и электропроводность [34]. Порошки ТРГ являются эффективными адсорбентами нефтепродуктов и других органических
соединений. Также известна возможность применения углеродных материалов, в том числе из ТРГ, в качестве катализаторов и носителей катализаторов.
Химическая или термическая очистка
Химическое или электрохимическое внедрение
Природный графит
термообработка в замкнутом объеме
Углеродные
изделия с
регулируемой
плотностью
Интсркалированные соединения графита
Огнезащитные материалы
термо эб] іаботка
\/
ТРГ
Уплотнительные изделия
Графитовая фольга
Порошки ТРГ и
композиты на
его основе
Сальники
Прокладки
Адсорбенты
Катализаторы и
носители катализаторов
Рис. 1. Схема получения ТРГ и его последующей переработки
Одним из новых перспективных способов получения изделий из ТРГ является терморасширеиие иитеркалированного (окисленного) графита в газопроницаемой пресс-форме с заданными размерами - так называемое химическое прессование. Технология процесса разрабатывается в Энгельсском
подразделении ЗАО «Унихимтек» под руководством профессора А.И. Финаенова [1]. Данный метод позволяет получать изделия с заданной плотностью заранее определенной геометрии [56]. Некоторые экспериментальные образцы, полученные данным методом, можно увидеть на фотографиях 1 и 2.
Фотография 1. ТРГ, полученный вспениванием в свободном объеме (слева) и изделие полученное химическим прессованием (справа)
На фотографии 1 представлен ТРГ, полученный вспениванием в свободном объеме (слева) и изделие, являющееся результатом химического прессования ТРГ в прямоугольной пресс-форме (справа). На фотографии 2 представлено изделие из ТРГ (слева), полученное вспениванием в технологической форме в виде полого цилиндра.
Изделия, получаемые в результате вспенивания графита в ограниченном объеме, представляют собой легкий теплостойкий материал, который можно использовать в ракетостроении, а также в других областях промышленности и коммунального хозяйства для нужд тепло- и шумоизоляции.
Фотография 2. Некоторые изделия из ТРГ
Важным при разработке технологии изготовления изделий с заданными свойствами и формой является моделирование процесса и определение его основных характеристик. В частности, необходимо знать полное время процесса, так как при его превышении начинается выгорание терморасширенного графита, потеря массы и ухудшение структуры. При прерывании же процесса до его завершения изделие будет некачественным из-за значительной неоднородности.
Механизм вспенивания, а также теплофизические свойства терморасширенного графита изучены достаточно мало. Одна одномерная модель вспенивания в свободном объеме огнезащитного покрытия рассмотрена в [21]. В ней не рассматривается движение вещества и постоянны теплофизические параметры. Роль математического моделирования процесса терморасширения графита достаточно весома. Актуальным является создание математической модели, когда процесс вспенивания представляется как фазовый переход с учетом движения графита.
Задачи теплопереноса с подвижными границами, вызванными изменением агрегатного состояния вещества, получили название задач типа Стефана [35,39]. Этот класс задач относится к одним из наиболее сложных задач математической физики. Классический вариант задачи Стефана, сформулированный для фазовых переходов типа плавление-кристаллизация, сводится к уравнению теплопроводности в области с заранее неизвестной границей, разделяющей твердую и жидкую фазы и имеющей температуру, равную температуре фазового превращения. Нелинейность задачи обусловлена наличием подвижной границы раздела фаз.
Существует небольшое число аналитических решений подобных задач, ограничивающихся простейшими случаями [30]. Наибольшее развитие получили численные методы решения, значительный вклад в разработку которых внесли в разное время А.А. Самарский [42,43], Б.М. Будак [12-15], П.Н. Вабищевич [17], J. Crank [55].
К особенностям математической постановки задачи терморасширения ОГ относится отличие на несколько порядков теплофизических параметров фаз графита и необходимость учета движения вещества. Для эффективного применения математической модели и оценки влияния конвективных членов в этом случае актуальным является разработка модификаций существующих численных методов.
Целью работы является построение математической модели процесса химического прессования окисленного графита, разработка методов расчета характеристик процесса. Для достижения поставленной цели в работе решаются следующие задачи: 1) построение математической модели процесса терморасширения окисленного графита; 2) разработка методов расчета характеристик процесса при существенном различии теплофизических параметров агрегатных состояний графита; 3) разработка модификации численных методов решения задачи типа Стефана в многокомпонентной среде.
Решение данных задач основано на идеях и подходах кинетической теории и механики сплошной среды для описания движущихся сред при учете фазовых переходов и тепловых воздействий. При построении модели процесса с учетом конвективного тепломассопереноса используются классические уравнения гидромеханики вязкой жидкости. Для решения возникающих задач используется численный конечно-разностный метод сквозного счета со сглаживанием теплофизических параметров в окрестности подвижной границы раздела фаз.
Научная новизна работы заключается в следующем: 1) построена математическая модель процесса химического прессования с учетом движения графита; 2) показана возможность замены в построенной модели стенок пресс-формы некоторым граничным условием третьего рода; 3) разработан алгоритм совместного решения задач определения полей температур и скоростей движения частиц среды; 4) осуществлена оценка влияния конвективных членов на процесс теплопередачи; 5) создана программа численного определения полей температур и скоростей.
Методы решения задачи Стефана
Задача Стефана представляет собой краевую задачу для параболических уравнений в областях со свободными границами, на которых заданы условия
Стефана материального или энергетического баланса. Граница раздела фаз движется внутри области по мере отвода или подвода теплоты к поверхности рассматриваемого объекта. При этом предполагается, что выделение или поглощение теплоты фазового перехода происходит в бесконечно тонкой области материала - на движущемся фронте (классическая постановка задачи). В математической постановке эта задача состоит в определении температурного распределения u{x,t) и фронтов фазовых переходов (одного или нескольких) в предположении, что все коэффициенты уравнений (1.1), (1.2) и условия Стефана (1.6), (1.7), все начальные (1.5) и граничные функции (1.3), (1.4) заданы [2,35,39].
Методам решения проблемы Стефана в классической постановке в литературе уделено значительное внимание. Для отдельных простых случаев можно построить точные аналитические решения, например, для задачи с плоским межфазным фронтом, постоянными граничными условиями первого рода и при равенстве начальной температуры температуре фазового перехода [30]. Несколько аналитических решений задачи теплопроводности с фазовым превращением в полубесконечной области с начальной температурой, равной температуре фазового перехода построены в [22,50,52,56,58].
В работах [24,49] предложено обобщение метода интегральных преобразований на нецилиндрическую область х є [О, (/)], / 0 (однофазная задача Стефана) для первой краевой задачи.
В серии работ [31,33,44] рассмотрен другой подход для области с произвольно движущейся внешней границей. Его практическое исследование предполагает вычисление производных любого порядка от выражений специального вида в общем члене ряда. Метод дает возможность получить аналитическое решение задачи при любом виде граничных условий. Его особенная эффективность заключается при решении обратных задач Стефана. Одной из проблем метода дифференциальных рядов является установление скорости сходимости рядов при определении температурной функции.
Метод функций Грина решения краевых задач для уравнений параболического типа также нашел свое применение для решения задач Стефана. В работе [19] предложен численно-аналитический метод, позволяющий находить приближенные решения одномерных нестационарных задач Стефана с учетом зависимости коэффициентов уравнения теплопроводности от температуры. Сущность данного подхода состоит в последовательном применении преобразования Кирхгофа искомой температуры, метода прямых [28] и метода функций Грина. Численное решение получаемой системы интегральных уравнений типа Гаммерштеина ищется проекционно-сеточным методом [33]. В работе [18] рассматриваются особенности применения этого подхода к решению двумерной двухфазной нестационарной задачи Стефана для любых граничных условий и областей, допускающих построение функций Грина первой или третьей краевой задачи для оператора Лапласа.
В работе [26] применительно к задаче промерзания влажного грунта найдено аналитическое решение двумерной нестационарной задачи Стефана с радиальной симметрией. Метод опирается на известное решение двумерной нестационарной задачи Дирихле для ортотропного полупространства с разрывными условиями первого рода.
Аналитические решения задачи плоской, цилиндрической и сферической формы роста кристалла из чистого расплава указаны в [29].
Из-за нелинейности задачи, обусловленной наличием подвижной границы раздела фаз, поиск аналитических решений задач Стефана в более общих случаях (многомерных, многофазных, с любыми граничными условиями) крайне сложен или даже невозможен. Наибольшее развитие получили численные методы, преимущественно в приложении к проблемам мерзлотоведения, геологии, металлургии. Существующие методы численного решения данной задачи можно условно разделить на два класса: 1) с явным выделением подвижной границы и 2) методы сквозного счета. К первой группе методов относятся методы, в которых положение свободной границы отслеживается на каждом временном слое. С этой целью используются численные методы, в которых свободная граница определяется положением соответствующих узлов. Это достигается за счет использования новых динамических независимых переменных или же согласованных динамических сеток в исходных переменных.
Численное решение одномерной задачи терморасширения графита методом сквозного счета
Для решения поставленной в п. 1.2 одномерной задачи терморасширения графита методом сглаживания используется обобщенная формулировка в виде уравнения (1.8) с исходными краевыми условиями (1.3)-(1.5). Особенность задачи в такой постановке заключается в наличии в левой части уравнения теплопроводности дельта-функции б(и-и,). Дискретизация задачи осуществляется с помощью конечных разностей. Для исключения особенности из условия задачи дельта-функция на некотором интервале температур н, -ДрН, + А2] заменяется приближенно дельтаобразной, или размазанной, дельта-функцией 5 0. Вне этого температурного интервала функция 8 равна нулю. Это размазывание, или, сглаживание, эквивалентно замене на интервале [«„-Л,,!/. + Д2] разрывной функции г(г/-г;ф), такой что dr /dE, = 8( ), непрерывной функцией rj(w-M.), такой что rj ( )= 8(Д). После такой замены вводится сглаженная удельная теплоемкость с(и) из условий: 1) c(ii)=cx(u) при и и, А{, с(и) = с2(и) при и и, + Д2 (то есть с{и)= с\и) вне интервала [и, —\,и, +Д2] ); 2) C(K) = C(W) + .8(W-W,) при ие[и, -А,,и» + Д2], где 8 определяется из условия сохранения баланса тепла на интервале и должна удовлетворять условию И.+d, J 5(u-u,)du = l. (1.10) и.-Д.
В построенной одномерной модели процесса терморасширения ОГ теплофизические параметры ОГ и ТРГ постоянны. В результате из (1.10) для с в интервале сглаживания получаем значение с=(Ь + сД + с2Д2)/(Ді+Д2)» (1.11) что соответствует линейной аппроксимации S = 1/(Д, + Д2).
В итоге в значении коэффициента теплоемкости внутри зоны размазывания учитывается скрытая теплота фазового перехода. Возможны другие способы задания сглаженной теплоемкости на интервале [и, -Д,«, + Д2], например, когда используется параболическая аппроксимация дельта-функции [17]: Г2(и-и.))Л ЧД1+А2У c(u) = c(u)+\b(u-u,), 3 1 5 = 2(Д,+Д2) 0, и(ф. -Д,,и +Д2] ,ие[и, -Д,,и, +Д2], (1.12) для которого условие (1.10) выполнено, также как и при линейной аппроксимации дельта-функции.
После такой замены никаких особых элементов в условии задачи уже нет и можно строить консервативную конечно-разностную схему, например, при помощи интегроинтерполяционпого метода [46]. Использовалась чисто неявная схема, которая является абсолютно устойчивой [40]. Дискретный аналог записанных в обобщенном виде (1.8) уравнений теплопроводности, полученный для неравномерной пространственной сетки имеет вид: n+1
В данной модели вспенивания ОГ коэффициенты каждой фазы с{,с2 не зависят от температуры, однако сглаженная теплоемкость всегда зависит от температуры. Поэтому для определения сеточной функции у на новом временном слое п +1 нужно проводить уточнение получаемого решения при помощи какого-либо итерационного процесса. Применялся простой метод итерационного уточнения коэффициентов, когда при вычислении (s +1) -ой итерации искомой функции у коэффициенты сяк вычисляются по значениям у на предыдущей итерации. В качестве начального приближения на новом временном слое естественным является выбор значения температуры, полученного на предыдущем слое п. Величину коэффициента теплоемкости для каждого временного слоя tntl в узлах сетки, в которых у" 1 є [и, -Avu, +А2], определяем по формуле (1.11) или (1.12). Иначе разностный аналог (1.13) уравнения теплопроводности на каждой итерации можно записать в более удобном виде
Данный дискретный аналог справедлив для всех узлов пространственной сетки, как неподвижных, так и подвижных, значения сеточной функции температуры в которых лежит в интервале сглаживания (однородная схема) [40]. Аппроксимированные условия на границах расчетных областей записываются в каноническом виде [10] Уі=ХіУ2+ У = X2 -i+v2, (1.16) где значения коэффициентов х, и v,, %2 и v2 определяются из соответствующих роду граничных условий формул. Для условия теплоизоляции на правом конце х = I при аппроксимации первой производной по х с помощью правой разностной производной получаем %2 = 1, v2 = 0. Аналогично для смешанного граничного условия (1.4) значения коэффициентов следующие:
В общем случае при нестационарных граничных режимах коэффициенты % и v являются функциями времени. Порядок точности описанной схемы зависит от способа аппроксимации производных в граничных условиях. В описанном выше варианте точность схемы имеет первый порядок по пространству и времени.
Метод сквозного счета численного решения двумерной задачи терморасширения графита
Для применения численного метода сквозного счета совместно с экономичными конечно-разностными схемами отметим, что задача Стефана допускает обобщенную формулировку в виде одного нелинейного уравнения, из которого следуют условия на подвижной границе (2.5), (2.6):
Уравнение (2.14) примечательно тем, что сама неизвестная граница фазового перехода явно не присутствует. Учет теплоты фазового перехода эквивалентен заданию на некотором интервале температур [и, -\,и, +Д2] эффективной теплоемкости: с=с + ХЬ(и-и.). (2.15) Аналогично одномерному случаю c(u) = cj (/ = 1,2) вне интервала [и, -Д[,«. + Д2], а внутри интервала величина с определяется по формуле (2.15), в которой функция б удовлетворяет условию (1.10).
Квазилинейное уравнение теплопроводности (2.14) лежит в основе эффективных вычислительных процедур приближенного решения задач типа Стефана. Опираясь на аппарат конечно-разностных схем, построим алгоритм сквозного счета.
В качестве исходной области выбираем прямоугольник D = {Oj ,, а2 у Ьг). В нем вводим произвольную неравномерную сетку h=\xi,yl),i = \,N\,j = \,N2,xx=a[,xNi = b[,yl=a2,y =b2). На временном отрезке определяем сетку cot ={tn,n = l,2,...} с шагом тя+1 =tn+] n. Для аппроксимации краевой задачи для двумерного уравнения (2.14) используем экономичную конечно-разностную локально-одномерную схему расщепления [41]: Значения коэффициентов \xn,\ik2,vt] vk2(k = \,2) определяются из соответствующих роду граничных условий формул (2.4). Причем они могут быть функциями координат и времени.
Поскольку уравнение (2.14) нелинейное, то для определения сеточной функции w на новом временном слое и + 1 нужно проводить уточнение получаемого решения при помощи какого-либо итерационного процесса. Используется простой метод итерационного уточнения коэффициентов, когда при вычислении (s + 1)-ой итерации искомой функции w коэффициенты Г, р и к вычисляются по значениям w на предыдущей итерации
При заданных коэффициентах уравнения (2.14) и известном поле скоростей итерационный алгоритм нахождения решения на новом временном слое w"+i следующий: 1) Находим промежуточное значение w. Для этого последовательно по строкам l j N2 на каждой строке при фиксированном значении индекса прогонкой по / = 2,Л -1 решаем одномерные сеточные уравнения (2.21) с граничными условиями (2.23). При этом матрицы теплофизических коэффициентов вычисляются по полученным на предыдущей итерации S значениям Wij. В качестве начального приближения s = 0 на новом временном слое используем значения, полученные на предыдущем слое w". Для узлов s-l (#,., .)є to, в которых и„-Д, w! +l н, + Д2 вычисляем эффективную теплоемкость по формуле (1.11). 2) По столбцам при каждом фиксированном \ i N1 прогонкой по j = 2,N2 -1 находим на этой итерации сеточную функцию w;"+1 из решения (2.22), (2.24). При этом для определения коэффициентов уравнений (2.22) используем значения теплофизических параметров, вычисленных в п.1). 3) Если выполняется условие max J s-] s, то получено решение s на этом слое wntt = w"+I, иначе J = s +1 и переходим к шагу 1) для вычисления температуры на новой итерации, s - заданная точность выхода из итерационного цикла.
Задача определения потенциала скоростей. Конечно-разностная аппроксимация
На сетке ША = {(xi,yi),i = l,Nl,j = \,N2,xl =anxNi = ,7, = altyNi =62 выберем разностную аппроксимацию производных в уравнении (3.8) такую, чтобы при суммировании (3.8) по узлам соА оставались только нормальные производные на границе, то есть, чтобы был справедлив разностный аналог формулы (3.11).
Для функции одной переменной Y(x) запишем разностный аналог производной произведения (g(x)Y (x)) х в двух вариантах: а) с помощью правой разностной производной, используя левую разностную производную для записи Y : б) с помощью левой разностной производной, используя правую разностную производную для записи Y :
Для простоты ограничимся случаем постоянного шага hi - h = const. Тогда для аппроксимации ( (х)У (л:)) х выберем полусумму величин, даваемых формулами а) и б): (g(x)r(x)) =-L[(g( +gM)YM -(gl+1 +2g, +gM)y +(g( + .,) 1 = -.(3.15) =i, /ft 2ft
Разложением в ряд Тэйлора можно показать, что правая часть формулы (3.15) может быть представлена в виде Y ig + Y"gi +0[h2). То есть, имеем аппроксимацию со вторым порядком.
Для выражения в правой части справедливо равенство ІУ, =( ,, +gyr.K -Vi)-( 2+ iXiWi)- (зле) i=2
Если производные на концах интервала по х равны нулю, то при разностной аппроксимации с первом порядком YN YN_VY2-YX и сумма в (3.16) обращается в нуль. Можно получить аналогичный результат, если использовать аппроксимацию граничных условий второго рода со вторым порядком.
В итоге для уравнения (3.8) на постоянной сетке получим аппроксимацию: [(PMJ +PuW/ (P/ u +2Pu +P -ukj +(Pu +PM,,WJ+ 2/J2 Jj (3.17) + Г[(РІ,У+І +Pu Wi "(Pu+i + +P/j-iKy +G u +Pu-i)Pu-i]BI Ля где ф,.у=ф(х(., -Д і = 2,...,Nl -\J = 2,...,N2 -1 . Просуммировав равенства (3.17) по / от 2 до JV, -1 и по _/ от 2 до JV2 -1: 2 J=2 X[(P ,J + Р нДф ,у -ф ,-и)-(р2 +РіЛч и "Фи)] + +777 Ewv, +Р 1-ІАФ 1 -Ф -IJ-IPW + Р,-ДФІ,2 — РЛІ Л= Z ХЛ 2« , (=2 (=2 у=2
При выполнении условия (3.9) и аппроксимации — на Г с помощью дп разностных производных с первым порядком получим, что должно выполняться условие, являющееся разностным аналогом (3.11): Рассмотрим случай, когда области D{{t) и Dt(t + At) могут пересекаться. Это возможно, если при терморасширения рассматривать движение частиц ОГ и ТРГ, которое присутствует в реальном процессе. Получим разностный аналог сохранения массы в целом. Как обычно межфазная граница о(/) делит исходную область на подобласти ),(/) и D2{}), занятые ОГ и ТРГ соответственно. Пусть A(t) = D2(t)KJv(t),B(t) = D\A, D+ = B{t)\B(t + At), D =B(t + At)\B(t), Dc =B{t)c\B{t + At) (см, рис. 13). Обозначим nl(t\n2(t),n+ n nc число внутренних узлов, лежащих в области A(t\B{t),D+,D,Dc. После аппроксимации уравнения Пуассона для определения потенциала скоростей ф в любой момент времени t t, имеем сеточное уравнение (3.17) и разностный аналог граничного условия (3.9):
При этом при задании правой части (3.17) согласно определению разностей плотностей по формуле (3.20), условие сохранения массы (3.18) выполнено.
Система алгебраических уравнений (3.17), (3.22) является вырожденной, поскольку задача Неймана имеет единственное решение с точностью до аддитивной постоянной. Опишем алгоритм определения частного решения данной задачи в прямоугольной области.
