Integral formulation of solutions of boundary-value problems of heat and mass transfer in domains with moving boundaries
- Authors: Shevelev V.V.1,2
-
Affiliations:
- MIREA – Russian technological University, Institute of fine chemical technologies. M.V. Lomonosov
- National research technological University "MISIS"
- Issue: Vol 29, No 1 (2021)
- Pages: 73-91
- Section: Informatics, Computer Science and Control
- URL: https://journals.eco-vector.com/1991-8542/article/view/65245
- DOI: https://doi.org/10.14498//tech.2021.1.6
Cite item
Full Text
Abstract
Using the Fourier transform, integral representations of solutions to boundary value problems of heat conduction and diffusion in a two-phase region with a moving interface are obtained. The proposed approach makes it possible to obtain the equation of motion of the interface without the need to first find the temperature and (or) concentration fields. This makes it possible to study the stability of the interface with respect to disturbances in its shape. The validity of the proposed approach is demonstrated by the example of self-similar growth of a spherical crystal in a supercooled melt and crystallization of the melt on a substrate of the same substance. On the basis of the obtained equation, which determines the rate of self-similar motion of the interface, the features of the kinetics of crystallization of the melt on the substrate are analyzed. The conditions of applicability of the developed approach to the solution of boundary value problems of heat conduction and diffusion in regions separated by a moving boundary are briefly discussed.
Full Text
Введение
Математические модели процессов тепломассопереноса играют важную роль в исследовании и прогнозировании технологических процессов, в которых для получения материалов с требуемыми свойствами используются процессы теплопроводности (термообработка, кристаллизация, плавление) и диффузии (процессы при фазовых превращениях в твердом состоянии в многокомпонентных системах, затвердевание растворов в расплаве, процессы насыщения материала технологически важным компонентном) [1–7].
Это связано с тем, что указанные математические модели позволяют исследовать влияние кинетики процессов тепломассопереноса на формирование структуры получаемого материала в зависимости как от способа внешнего воздействия на материал, так и от его спектра времен релаксации в условиях указанного внешнего воздействия.
Одним из важнейших технологических приемов создания материалов с заданными свойствами является формирование структуры материала в результате протекания в нем фазовых превращений. При этом кинетика фазовых превращений является одной из важнейших компонент в технологии получения материалов различного назначения, в том числе используемых в электротехнике, приборостроении и вычислительной технике. Это связано прежде всего с тем, что она, как правило, играет существенную роль в формировании структуры материала, определяющей его различные эксплуатационные свойства.
Структура материалов, формирующаяся в процессе протекания в них фазовых превращений, зависит от конечного распределения в них центров новой фазы по размерам. Это распределение определяется следующими двумя процессами: зарождением центров новой фазы и процессом их последующего роста. Первый из указанных процессов моделируется начально-краевыми задачами для уравнения типа Фоккера – Планка, второй – краевыми задачами для уравнений теплопроводности или/и диффузии в областях с движущимися границами или моделями фазового поля. Второй из указанных процессов играет не менее важную роль в формировании конечной структуры материала, чем процесс зарождения центров новой фазы, ввиду известной проблемы устойчивости межфазной границы раздела фаз, приводящей к формированию дендритной структуры материала. Имеется огромное число работ, посвященных исследованию устойчивости межфазной границы, разработке аналитических и численных методов решения соответствующих краевых задач [8–18].
Одним из эффективных аналитических методов решения краевых задач теплопроводности и диффузии, моделирующих кинетику фазовых превращений, является метод интегрального представления соответствующих процессу температурных и концентрационных полей. Этот метод предпочтителен в том случае, когда основной задачей является определение скорости движения межфазной границы, так как он позволяет получить уравнение движения межфазной границы без необходимости предварительного нахождения температурных и(или) концентрационных полей.
В настоящее время в многочисленных работах, посвященных интегральным представлениям краевых задач для областей с движущимися границами, моделирующих кинетику фазовых превращений, указанные представления получены в основном для случая, когда процессы тепло-и/или массопереноса, обусловленные фазовым превращением, происходят только в одной из областей двухфазной области.
Однако в общем случае процессы тепло- и/или массопереноса, обусловленные фазовым превращением, происходят, как правило, в обеих областях двухфазной области и их необходимо учитывать в математических моделях кинетики фазовых превращений, особенно в тех случаях, когда характерные времена процессов переноса – одного порядка. Насколько нам известно, для случая двухфазных областей решения в интегральной форме получены только в отдельных частных случаях [19, 20].
В данной работе на основе интегрального преобразования Фурье предложена интегральная формулировка решения связанных краевых задач для уравнения теплопроводности для двухфазной области с движущейся межфазной границей, моделирующих процесс кристаллизации однокомпонентного расплава. Аналогичными краевыми задачами моделируются и процессы фазового превращения в растворах (твердых и жидких). Так что излагаемый ниже подход может быть применен и для интегральной формулировки решения связанных краевых задач для уравнения диффузии.
Формулировка модели
Типичная математическая модель роста центра новой фазы из расплава в результате теплоотвода через переохлажденный расплав, когда начальная температура в обеих фазах постоянна, а исходная фаза имеет неограниченный объем, имеет вид:
(1)
(2)
(3)
(4)
(5)
(6)
(7)
(8)
(9)
Здесь , - температура в точке области , занятой твердой фазой в момент времени t;
- начальная температура твердой фазы;
, - температура в точке области , занятой жидкой фазой (расплавом) в момент времени t;
- начальная температура жидкой фазы;
; есть переохлаждение жидкой фазы у межфазной поверхности;
температура равновесия между жидкой и твердой фазами, зависящая в общем случае от кривизны межфазной границы;
температуропроводности твердой и жидкой фаз соответственно;
теплопроводности твердой и жидкой фаз соответственно;
нормальная скорость движения межфазной границы;
– плотность твердой фазы;
удельная теплота плавления на единицу объема;
поверхность межфазной границы.
Начальное условие (2) сформулировано в предположении, что в начальный момент времени твердая фаза отсутствовала.
Уравнение (3) есть граничное условие, определяющее температуру на межфазной поверхности. В общем случае это кинетическое условие, которое формулируется в зависимости от механизма роста кристалла. В данном случае граничное условие (3) записано в предположении нормального механизма роста кристалла [1, 2], который реализуется при достаточно больших переохлаждениях жидкой фазы. Условие (4) есть условие равенства температур твердой и жидкой фаз на межфазной поверхности.
При заданной скорости роста условий (1)–(6), (8), (9) достаточно для однозначного определения соответствующих температурных полей. Однако если скорость роста не задана, то необходимо задавать дополнительное условие для ее определения. Таким условием является условие теплового баланса на межфазной поверхности, которое в общем виде записывается в виде, аналогичном уравнению (7).
Условие теплового баланса существенно усложняет возможность получения аналитического решения соответствующей краевой задачи, так как необходимо согласованным образом найти и температурное поле, и зависимость , связанные друг с другом. В результате для определения зависимости необходимо в общем случае решить интегро-дифференциальное уравнение, получающееся после подстановки выражений температурных полей в условие массового баланса [1]. Приближенные аналитические методы решения самосогласованных краевых задач для уравнений теплопроводности и диффузии представлены в [1].
Метод решения
Решение краевых задач (1)–(9) ищем с помощью преобразования Фурье, полагая, что функции и вне областей определения и в своих краевых задачах соответственно равны нулю, а начало координат находится в области .Тогда изображения этих функций определяются следующими формулами:
(10)
(11)
Преобразуем далее уравнения краевых задач (1)–(9) по формулам (10), (11) преобразования Фурье с учетом формул Грина и того, что
,
где граница области ;
радиус-вектор точки М, принадлежащей поверхности ;
нормальная скорость движения границы области ;
внешняя нормаль к .
Получим в результате уравнения
(12)
(13)
Здесь внешняя нормаль к относительно области ,
радиус-вектор точки М.
Начальные условия для этих задач Коши имеют вид, в рамках предположения о постоянстве начальных температур, , . Тогда решения уравнений (12), (13) можно представить в виде
,.
Подставляя в эти формулы соответствующие выражения и из формул (12) и (13), и далее и в формулы обращения преобразования Фурье
,
получим следующее интегральное представление решений краевых задач (1)-(9):
(14)
(15)
где ,, – радиус-вектор точки M ;
радиус-вектор точки Р, а интегрирование по выполняется по каждой компоненте вектора по промежутку .
Выполняя далее в (14), (15) стандартное интегрирование по переходом к сферической системе координат в -пространстве, получим
.
С учетом этого результата выражения (14) и (15) преобразуются к виду
(16)
(17)
Формулы (16), (17) являются искомыми интегральными представлениями температурных полей в каждой из фаз двухфазной области. Они особенно удобны для расчета температурных полей, когда заданы температура на межфазной границе и нормальная скорость ее движения. В случае процесса роста центра новой фазы, когда температура на межфазной границе и скорость ее роста не заданы, а подлежат самосогласованному определению из условий равенства температур на межфазной границе (4) и теплового баланса (7), мы приходим к системе двух интегро-дифференциальных уравнений, из которой необходимо найти , и нормальную скорость движения межфазной границы .
Для аналитического решения таких систем обычно используют либо метод разложения по малому параметру, либо метод линеаризации [3], когда отклонение температурных полей в обеих фазах от известных решений является малым возмущением задачи.
Иллюстрация метода
Однофазная область. В случае однофазной области изложенный подход приводит к следующему выражению для температурного поля:
(18)
где – удельная теплоемкость расплава;
плотность расплава.
В частности, в случае роста сферического кристалла в переохлажденном расплаве из формулы (18), полагая и интегрируя по от 0 до , получим:
(19)
Здесь есть угол между и ; , есть значение на межфазной границе, то есть . Проинтегрируем далее по :
Подставив полученные результаты в (19), получим
(20)
Выражение (20) для температурного поля является решением следующей краевой задачи, вытекающей из формул (1)–(9) и рассмотренной в [1]:
(21)
(22)
(23)
(24)
(25)
Краевая задача (21)–(25) является математической моделью процесса роста сферического кристалла из переохлажденного расплава [1].
В том случае, когда процесс затвердевания расплава контролируется теплопроводностью, температура на межфазной границе близка температуре кристаллизации, то есть переохлаждение на межфазной границе практически снимается и , где есть исходное переохлаждение расплава, из формулы (19) получается известное автомодельное решение с законом движения межфазной границы . Для этого необходимо найти исходя из выражения (20) и полагая затем , далее вычислить достаточно сложные интегралы и подставить результат в уравнение (24). В результате придем к следующему уравнению для параметра β, определяющего кинетику роста сферического кристалла из переохлажденного расплава (выкладки мы не приводим ввиду их громоздкости):
. (26)
Здесь .
Значения левой части уравнения (26) заключены между нулем и единицей [1]. Следовательно, уравнение (26) имеет решение только при условии , которое, таким образом, является и условием существования автомодельного решения.
Двухфазная область. В качестве другой иллюстрации метода рассмотрим двухфазную одномерную краевую задачу для теплопроводности, моделирующую процесс кристаллизации на подложке из того же вещества, что и расплав, лежащий в основе различных технологий получения материалов с требуемыми свойствами [1]. Математическая модель указанного процесса имеет следующий вид:
(27)
(28)
(29)
(30)
(31)
(32)
, (33)
, (34)
. (35)
Здесь , температура в твердой фазе в точке x в момент времени t, начальная температура в твердой фазе; , температура в жидкой фазе в точке x в момент времени t, начальная температура в жидкой фазе; ,.
Запишем решения краевых задач исходя из формул (16), (17) и предполагая, что температура на межфазной границе не зависит от координат y и z. Тогда, интегрируя в (16), (17) по y и z, получим
(36)
(37)
Рассмотрим автомодельный случай, когда процесс затвердевания контролируется теплопроводностью, а не кинетическими процессами на межфазной границе. В этом случае , .
Тогда, так как , то
, (38)
. (39)
Подставив (38), (39) в уравнение (33), получим
(40)
Уравнение (40) содержит три неизвестных величины . Для того чтобы найти параметр , определяющий кинетику процесса кристаллизации, необходимо определить величины . Их мы найдем, удовлетворяя уравнения (38), (39) с помощью полученных решений (36), (37). Однако прежде чем найти , преобразуем выражения (36), (37) с учетом того, что .
Это позволяет объединить второе и третье слагаемые в (36), (37). В результате такого объединения формулы (36), (37) можно привести с учетом того, что , к следующему виду:
(41)
(42)
Тогда
(43)
(44)
В выражения (43), (44) входят однотипные интегралы, которые также появляются при нахождении автомодельного решения рассмотренного выше однофазного случая. Интегралы такого типа подробно рассмотрены в [1] и не приводятся в известных таблицах определенных интегралов. Поэтому рассмотрим метод их вычисления. Для этого ввиду их однотипности ограничимся рассмотрением интегралов
и для его вычисления вводим переменную интегрирования , после чего получим
.
Полученный интеграл вычисляется подстановкой .
Тогда, учитывая, что
,
получим после соответствующего преобразования подынтегрального выражения
,
где .
Учтем полученный результат в формулах (43), (44). Получим с учетом (38), (39) соответственно
(45)
(46)
Из уравнений (45), (46) получим следующие выражения для искомых величин :
(47)
(48)
Подставив выражения (47), (48) в (40), получим искомое уравнение для параметра , определяющее кинетику процесса кристаллизации:
(49)
Уравнение (49) позволяет сделать некоторые выводы относительно кинетики кристаллизации расплава на подложке. В частности, из уравнения (49) следует, что при , то есть когда температура кристаллизации (плавления) меньше и начальной температуры подложки , и начальной температуры расплава , кристаллизация невозможна, а происходит, наоборот, процесс плавления подложки.
При условии же, когда , то есть когда температура кристаллизации выше и начальной температуры подложки , и начальной температуры расплава , процесс кристаллизации гарантированно имеет место и происходит с тем большей скоростью, чем больше указанные начальные переохлаждения подложки и расплава . В этом случае теплоотвод от межфазной границы происходит и в подложку, и в переохлажденный расплав.
Рассмотрим случай, когда , то есть когда начальная температура подложки больше температуры кристаллизации , а начальная температура расплава меньше температуры кристаллизации . Тогда при достаточно малых уравнение (49) примет вид
. (50)
Здесь – соответственно удельная теплоемкость и плотность твердой и жидкой фаз.
Из уравнения (50) следует, что расплав кристаллизуется на подложке при условии , то есть при условии
.
Это условие означает, что расплав кристаллизуется на подложке в результате теплоотвода через расплав. При условии
,
наоборот, плавится подложка, а в случае, когда
,
система расплав – подложка находится в состоянии кинетического равновесия, так как .
В противоположном случае, когда , расплав кристаллизуется на подложке в результате теплоотвода через подложку при условии , то есть при условии
.
При условии
,
наоборот, плавится подложка.
Таким образом, рассмотренные выше примеры показывают, что интегральные представления решений краевых задач теплопроводности и диффузии, моделирующих технологические процессы получения материалов с требуемыми свойствами, позволяют предсказывать кинетику указанных процессов и условия ее реализации без предварительного нахождения температурных и(или) концентрационных полей, что представляет определенные преимущества перед численными методами решения указанных задач.
Выводы
Изложенный подход к интегральной формулировке решений краевых задач теплопроводности или диффузии может быть легко обобщен на случай неоднородного по объему распределения начальной температуры в твердой и жидкой фазах, обусловленного разного рода случайными факторами в технологии производства материалов с помощью затвердевания.
Развитый метод нахождения интегральных представлений решений указанных задач может быть также распространен и на случай, когда исходный объем новой фазы не равен нулю.
Кроме того, метод позволяет также получить интегральные представления решений задач, моделирующих кристаллизацию слитков, то есть кристаллизацию в больших объемах [1, 2], что представляет отдельный технологический интерес.
Предложенный метод проиллюстрирован на примере решения краевых задач для уравнения теплопроводности, однако без существенных изменений он легко переносится на аналогичные краевые задачи для уравнения диффузии.
Развитый в работе метод получения интегральных представлений краевых задач для уравнения теплопроводности и диффузии применим прежде всего в тех случаях, когда на подвижной границе раздела областей переноса заданы значения как искомых функций, так и их нормальных производных, что как раз и имеет место в математических моделях кинетики фазовых превращений.
About the authors
Valentin V. Shevelev
MIREA – Russian technological University, Institute of fine chemical technologies. M.V. Lomonosov; National research technological University "MISIS"
Author for correspondence.
Email: valeshevelev@yandex.ru
Professor
Russian Federation, 78, pr-ct Vernadsky, Moscow, 119454; 4, Leninsky pr-ct, Moscow, 119049References
- Lyubov B.Ya. Teoriya kristallizatsii v bol'shikh ob"yemakh. M.: Nauka, 1975. 256 р. (In Russian).
- Lyubov B.Ya. Kineticheskaya teoriya fazovykh prevrashcheniy. M.: Metallurgiya,1969. 264 р. (In Russian).
- Lyubov B.Ya., Shevelev V.V. Analiticheskiy raschet kinetiki diffuzionnogo rastvoreniya sfericheskogo vydeleniya inoy fazy. Fizika metallov i metallovedeniye. 1973. Vol. 35. № 2. Рр. 330–335. (In Russian).
- Temkin D.E., Shevelev V.V. Kinetics of nucleation in two-component systems. Journal of Crystal Growth. 1984. Vol. 66. № 2. Pp. 380–400.
- Ivantsov G.P. O roste sfericheskogo i igloobraznogo kristallov binarnogo splava. DAN SSSR. 1952. Vol. 83. № 4. Рр. 573–577. (In Russian).
- Temkin D.E. Crystallization process. Consultant Bureau, New York, 1966. Р. 15.
- Kristian Dzh. Teoriya prevrashcheniy v metallakh i splavakh. ch. 1. Termodinamika i obshchaya kineticheskaya teoriya. M.: Mir, 1978. 806 р. (In Russian).
- Mullins W.W., Sekerka R.F. Morphological stability of a particle growing by diffusion or heat flow. J. Appl. Phys., 1963. Vol. 34. № 1. Pp. 323–329.
- Mullins W.W., Sekerka R.F. Stability of a planar interface during so-lidification of a dilute binary alloy. J. Appl. Phys, 1964. Vol. 35. № 2. Pp. 444–451.
- Caroli B., Caroli C., Roulet B., Langer J.S. Solvability condition for needle crystals at large undercooling in a nonlocal model of solidification. Phys. Rev. A, 1986. Vol. 33. № 1. Pp. 442–452.
- Benamar M., Bouissou Ph., Pelce P. An exact solution for the shape of a crystal growing in a forced flow. Journal of Crystal Growth, 1988. Vol. 92. № 1–2. Pp. 97–100.
- Saville D.A., Beaghton P.J. Growth of needle-shaped crystals in the presence of convection. Phys. Rev. A, 1988. Vol. 37. № 9. Pp. 3423–3430.
- William R.A., Gill N. Dendritic growth with thermal convection. Journal of Crystal Growth, 1988. Vol. 91. № 4. Pp. 587–598.
- Jian-Jun Xu. Global neutral stability of dendrite growth and selection condition of tip velocity. Journal of Crystal Growth, 1990. Vol. 100. № 3. Pp. 481–490.
- Jian-Jun Xu. Interfacial wave theory of dendrite growth: global mode solution and quantum condition. Canadian Journal of Physics, 1990. Vol. 68. № 11. Pp. 58–66.
- Langer J.S. Dendritic sidebranching in the three-dimensional symmetric model in presence of noise. Phys. Rev. A, 1987. Vol. 36. № 9. Pp. 3350–3358.
- Jian-Jun Xu. Dendritic growth from a melt in an external flow: uniformly valid asymptotic solution for the steady state. Journal of Flu-id Mechanics, 1994. Vol. 263. № 3. Pp. 227–244.
- Lia J., Tang Y.L., Shenc N., Pana W. Effects of solidification kinetics on phase selection of Ni–Al alloys. Journal of Alloys and Com-pounds, 2001. Vol. 329. № 1–2. Pp. 157–161.
- Shevelev V.V., Lokshin Dzh.L. O nakhozhdenii zakona dvizheniya mezhfaznoy granitsy v zadachakh, modeliruyushchikh kinetiku fazovykh prevrashcheniy. Izvestiya vuzov i energeticheskikh ob"yedineniy SNG. Seriya «Energetika». 2004. № 3. Рр. 44–51. (In Russian).
- Shevelev V.V. Analiticheskiy raschet vremeni, limitiruyemogo diffuziyey rastvoreniya sfericheskoy chastitsy vtoroy fazy. Izvestiya AN SSSR. Seriya Metally. 1988. № 6. Рр. 57–60. (In Russian).
