Разработана и апробирована методика моделирования акустических сигналов, инициируемых лазерным импульсным излучением в композитах прозрачная матрица - наночастицы металлов с учетом процессов плавления. Методика заключается в расчете функции источников давления в зависимости от времени и координаты и ее свертки с функцией Грина одномерного волнового уравнения. Апробация проведена на практически важных композитах на основе пентаэритрит тетранитрата с наночастицами алюминия радиуса 50 нм. Плавление, характеризующееся повышением удельного объема, приводит к увеличению амплитуды максимума функции источников и появлению области отрицательных значений. Рассчитаны зависимости эффективной константы роста оптоакустического сигнала и его амплитуды от плотности энергии импульса, которые необходимо учитывать при исследовании оптических свойств наночастиц данным методом. Результаты работы важны для развития методов неразрушающего контроля и прогнозирования особенностей функционирования устройств фотоники и оптических детонаторов, содержащих наночастицы.
Method of modeling of the optoacoustic signals in composites based on transparent matrix - metal nanoparticles.pdf Введение Оптоакустическая спектроскопия [1] является перспективным методом исследования, дополняющим классические методы спектрофотометрии. Измерение акустических сигналов при воздействии лазерного импульса на вещество позволяет исследовать оптические характеристики образцов с высокой оптической плотностью [2, 3] и высокой пористостью [4]. Развиваются методы оптоакустической спектроскопии в медицине, основанные на анализе акустических сигналов тканей организмов [5]. Перспективно использование методов оптоакустической спектроскопии при исследовании оптических свойств композитов на основе прозрачной матрицы и светопоглощающих и светорассеивающих наночастиц [6]. В экспериментальных работах [2, 7] были исследованы закономерности оптоакустических сигналов в композитах на основе пентаэритрит тетранитрата (PETN), допированного наночастицами металлов. Лазерное излучение наносекундной длительности эффективно нагревает наночастицы, вызывая появление импульса давления в поглощающем слое вещества. Отличием возникающих эффектов от рассматриваемых в классических моделях оптоакустических процессов [1] является достижение высоких температур локального нагрева [7, 8]. В этом случае нельзя пренебрегать изменением теплофизических характеристик вещества при увеличении температуры и в точках фазовых переходов. Цель работы заключается в разработке методики математического моделирования закономерностей оптоакустического отклика композитов прозрачная матрица - наночастицы металлов в условиях импульсного возбуждения с учетом изменения теплофизических свойств наночастиц и матрицы при плавлении. Апробация методики выполнена с использованием параметров композитов PETN - наночастицы алюминия. Наночастицы алюминия активно используются в устройствах нелинейной оптики [9, 10], перспективны для создания капсюлей оптических детонаторов на основе PETN (штатное бризантное взрывчатое вещество) [7]. Во всех перечисленных приложениях функциональные характеристики определяются оптическими свойствами наночастиц в прозрачной матрице. В настоящей работе будет рассмотрено формирование оптоакустической волны при поглощении импульсного излучения наночастицами металла в прозрачной матрице в одномерном приближении. При импульсном облучении системы кинетическая зависимость акустического сигнала разбивается на три области: 1) «мертвое время», связанное с распространением первичной волны от области возбуждения до детектора; 2) первичная волна, отражающая распределение тепловых источников; 3) вторичные волны, связанные с отражением первичной волны от границ образца и конструктивных элементов экспериментальной ячейки и датчика [1, 2, 7]. В дальнейшем изложении будет рассматриваться только вторая область. Модель расчета оптоакустического сигнала Рассмотрим формирование акустической волны в композитах PETN - наночастицы металла, облученных наносекундным лазерным импульсом. Взаимодействие излучения с веществом включает его поглощение и рассеяние наночастицами, отражение от передней и задней поверхности образца; приводит к нагреванию наночастиц и, в результате кондуктивного теплопереноса, - к повышению температуры слоя вещества матрицы, прилегающего к наночастице. Локальное повышение температуры стимулирует термическое расширение вещества и генерацию продольной акустической волны. Повышение температуры при достаточной плотности энергии импульса приводит к плавлению матрицы и наночастицы. Для большинства веществ плавление сопровождается увеличением объема. Рост объема наночастицы и окружающей матрицы должен вызывать деформацию материала холодной матрицы между наночастицами. Оценим характерное время, требуемое для деформации матрицы и достижения локального механического равновесия вещества в пространстве между наночастицами. Характерное расстояние между наночастицами примерно соответствует половине обратного кубического корня из концентрации (шт./см3). Механическое возмущение в среде распространяется со скоростью звука, поэтому время механической релаксации можно оценить как , (1) где ρPETN и ρAl - плотности матрицы (PETN) и наночастиц (алюминий) соответственно; ωAl - массовая доля наночастиц алюминия в образце; R - радиус наночастиц (50 нм); cs = 2330 м/с - скорость звука в PETN. При ωAl = 0.03 % и R = 50 нм [5-7] типичное время установления локального механического равновесия составляет τ = 350 пс. Длительность импульса на полувысоте современных лазерных стендов, используемых для исследования взрывных характеристик прессованных таблеток PETN, сенсибилизированных наночастицами алюминия, составляет ti = 14 нс [2, 7]. Поскольку ti >> τ, установление локального механического равновесия в пространстве между наночастицами является быстрым процессом, слабо влияющим на кинетику оптоакустического сигнала. В [7] был исследован диапазон массовой доли наночастиц от 0.0025 до 0.5 %. Даже в случае наименьшей ωAl = 0.0025 % величина τ = 0.8 нс. Таким образом, неравенство ti >> τ выполняется во всем актуальном диапазоне концентраций наночастиц. Температура плавления слабо чувствительна к изменению давления в диапазоне около атмосферного, поэтому можно вначале выполнить расчеты распределений температуры и движения фронтов плавления вокруг наночастиц на различных глубинах в образце, затем на их основе рассчитать функции источников давления, после чего получить сам оптоакустический сигнал как свертку функции источников с функцией Грина для акустической волны. В работах [2, 5-7] диаметр пучка лазера на поверхности образца составлял 2.5-3.0 мм, типичный показатель поглощения ~ 100 см-1, что дает возможность использовать одномерное приближение при моделировании акустических эффектов. Для расчета температурных полей и движения фронтов плавления в системе наночастица металла - матрица учитывались процессы нагревания наночастицы импульсом лазера, кондуктивный перенос тепла и плавление материалов матрицы и наночастицы [11, 12]. Плотность энергии, поглощаемая наночастицей, находящейся на глубине z, рассчитывалась из закона Бугера, сечение поглощения излучения считалось равным геометрическому сечению сферической наночастицы. Система дифференциальных уравнений имеет следующий вид: (2) где T - температура; x - пространственная координата, отсчитываемая от центра наночастицы; α и αM - коэффициенты температуропроводности матрицы и металла. Применимость сферической симметрии связана с особенностями взаимодействия света с наночастицами металла в прозрачной матрице. Многократное хаотическое рассеяние света приводит к усреднению освещенности по направлениям, что и позволяет использовать (как и в работах [8, 11]) при расчетах сферическую симметрию. Данное положение было экспериментально доказано в прессованных таблетках PETN - алюминий [12]. На границе наночастица - матрица x = R происходит поглощение энергии излучения, что приводит к граничному условию: , (3) где c и cM - объемные теплоемкости матрицы и металла; J(t) - зависимость мощности излучения лазерного импульса от времени. Зависимость J(t) близка к функции нормального распределения [2, 7]. Принимая за начало отсчета времени положение максимальной интенсивности импульса, для величины J(t) получаем выражение [8] , (4) где H - плотность энергии лазерного импульса, падающего на наночастицу; ki - параметр, определяющий длительность импульса на полувысоте (ti) ( ). На границе рассматриваемой области (включение и слой энергетического материала толщиной 7R) задавалось граничное условие первого рода T = 300 K. Движение фронта плавления матрицы моделировалось в рамках задачи Стефана [13, 14]. Температура на фазовой границе принималась равной температуре плавления, скорость движения границы xL рассчитывалась по выражению [13-15] , (5) где ΔHm - молярная энтальпия плавления вещества матрицы; Vm - молярный объем вещества матрицы. В выражении (5) пренебрегается различием коэффициентов температуропроводности и теплоемкости твердой и жидкой матрицы, поскольку для жидкого состояния литературные данные отсутствуют. Для расчета тепловых потоков в (5) распределение температуры аппроксимировалось выражением по обе стороны от фазовой границы: , (6) где ТL - температура на границе матрица - наночастица, которая в рамках классической задачи Стефана совпадает с температурой плавления; ai - подгоночные параметры. Для аппроксимации использовался взвешенный метод наименьших квадратов, весовые коэффициенты ближайших к фазовой границе трех точек принимались равными единице, далее весовые коэффициенты при четырех следующих значениях уменьшались, составляя 0.95, 0.90, 0.80, 0.70, а затем в геометрической прогрессии со знаменателем exp(-0.5). Расчет градиента с учетом дополнительных точек предлагался в [15], отмечалось, что это повышает точность вычислений. В случае моделирования плавления наночастицы из-за высокого значения температуропроводности, которое приводит к практически однородному распределению температуры, использование разницы тепловых потоков для расчета скорости фронта плавления затруднительно. Поэтому движение фазовой границы в металлической наночастице рассматривалось как квазихимическая реакция перемещения атомов из жидкой фазы в твердую. Для скорости перемещения границы использовалось уравнение, подобное кинетическому уравнению для продукта обратимой реакции типа [16]: , (7) где xLM - положение границы плавления в металле; - постоянная Больцмана; и - энтальпия и температура плавления металла; ћ - постоянная Планка; NL - число Лошмидта; γ обозначает величину , (8) где и - энтропия и энтальпия активации перехода молекулы из жидкой фазы в твердую. При расчетах принималось γ = 0.1 - значение, оцененное в работе [8]. Для расчета температуры на границе плавления в металле температурное поле по обе стороны от фазовой границы аппроксимировалось следующим образом: , (9) где x - пространственная координата, отсчитываемая от центра наночастицы; R - радиус наночастицы; bi - подгоночные параметры. Затем вычислялись экстраполированные величины температуры на границе, которая далее получалась как среднее гармоническое: . (10) Данные подходы к вычислению скорости фазовой границы уточняют методику [8] с повышением устойчивости решения. В частности, значительно сокращается резкое изменение скорости границы при пересечении ею центра расчетной ячейки. В результате в каждый момент времени после начала действия импульса рассчитывались температурные поля в наночастице и окружающем пространстве вместе с положением фазовых границ плавления. Рассмотрим методику расчета кинетической зависимости давления p исходя из приближения локального механического равновесия. Система может содержать материалы матрицы и металла как в твердом, так и жидком состоянии при различных температурах, кроме того, пространство между наночастицами можно считать не нагревающимся, а только испытывающим деформацию со стороны расширяющегося материала. Тогда суммарный объем в расчете на одну наночастицу запишется следующим образом: , (11) где j - индекс, нумерующий фазы (1 - матрица твердая, 2 - матрица жидкая, 3 - металл твердый, 4 - металл жидкий); Vj - область пространства, занимаемая фазой j; T(x) - распределение температуры; wj(T) - температурная зависимость относительного изменения объема для фазы j, возникающая из-за термического расширения; V0 - объем матрицы, в котором увеличение температуры не рассматривалось (за пределами области, рассчитываемой выше). Выражая давление, получим . (12) В качестве функций wj(T) использовались данные [17, 18] для алюминия и коэффициент термического расширения PETN [19] для матрицы, изменение объема матрицы в точке плавления считалось равным δVL = 5-15 %, что соответствует данной величине для ряда органических веществ [17]: (13) где - функция Хэвисайда. Далее рассчитываемое по (12) давление численно дифференцировалось по времени и выполнялось интегральное преобразование вида , (14) где τ - характерное время релаксации давления между наночастицами, вычисляемое по выражению (1). Данное преобразование учитывало запаздывание в выравнивании давления между наночастицами. Функция источников давления S удовлетворяет следующему соотношению, которое можно использовать для проверки результатов расчета: , (15) где величина конечного давления pf оценивалась по уравнению (12) при подстановке в него конечной температуры, которую достигает рассматриваемая область пространства после достижения теплового равновесия: , (16) где R1 - положение границы расчетной области. Это же давление использовалось для построения «идеализированной зависимости» функции источников, которая отвечает среде с непосредственным переходом энергии излучения в давление: . (17) Точность выполнения соотношения (15) составляла около 0.01 %. Для расчета оптоакустического сигнала использовалась функция Грина для одномерного волнового уравнения, в которой уже вычислена свертка по времени [1]: , (18) где zs - точка наблюдения, которая бралась на расстоянии 5/μ от поверхности. Расчет проводился следующим образом. Задавалось значение плотности энергии на поверхности образца. Показатель поглощения определялся концентрацией наночастиц (например, 0.1 %), сечение поглощения которых считалось равным геометрическому значению μ = πR2С. Пространство разбивалось на 85 ячеек равной толщины, составлявшей 0.05/μ, и выполнялся расчет функции источников для каждой ячейки. Для последних нескольких ячеек при использованных параметрах модели фазовых переходов не возникало. Функции источников сохранялись. При вычислении оптоакустического сигнала сохраненные функции источников считывались и подставлялись в (18). Для кинетических зависимостей оптоакустического сигнала рассчитывались амплитуда и эффективная константа скорости нарастания, которую определяли как коэффициент наклона в полулогарифмических координатах в области от 5 до 50 % амплитуды [2]. Результаты и их обсуждение Рассмотрим, каким образом изменяется локальная функция источников в зависимости от плотности энергии импульса и как возникающие особенности связаны с кинетикой нагревания наночастиц в рамках модели (1) - (18) при величине изменения объема матрицы в точке плавления δVL = 15 %. В простейшем случае малой плотности энергии (0.005 Дж/см2), когда отсутствуют фазовые превращения, кинетические зависимости температуры на границе наночастица - матрица и функции источников давления показаны на рис. 1, а. Зависимость температуры от времени не имеет особенностей и подобна кривым, рассчитанным в работе [8]. Соответствующая функция источников отличается от идеализированной зависимости меньшей амплитудой максимума и затянутостью заднего фронта. Данный эффект возникает из-за того, что коэффициент температуропроводности металла больше коэффициента температуропроводности матрицы. Наночастица раньше прогревается и становится основным источником давления. Затем прогревается матрица, что приводит к дополнительному росту давления из-за ее большего коэффициента термического расширения. При повышении плотности энергии до 0.01 Дж/см2 к процессам теплопередачи добавляется плавление матрицы. Плавление приводит к небольшому разрыву производной на стадии охлаждения при температуре плавления матрицы 413 К (рис. 1, б). Функция источников меняется значительно сильнее. Начало плавления приводит к резкому ускорению роста давления, затем на этапе затвердевания появляется область, в которой функция источников отрицательна. Данный эффект связан с уменьшением объема при затвердевании при передаче тепла на нагрев более глубоких слоев матрицы. Рис. 1. Рассчитанные зависимости функции источников в случае нагрева наночастицы (кр. 1), идеализированной системы (кр. 2) и температуры на границе наночастица - матрица (кр. 3) от времени при значениях плотности энергии, поглощаемой наночастицей 0.005 (а), 0.01 (б), 0.05 (в), 0.06 (г), 0.1 (д) и 0.3 Дж/см2 (е) Дальнейшее повышение плотности энергии делает роль плавления более выраженной. При Н = 0.05 Дж/см2 (рис. 1, в) процесс затвердевания матрицы проявляется в виде участка практически постоянной температуры на зависимости T(t), который свидетельствует об увеличении толщины проплавленного слоя по сравнению с результатами предыдущего расчета. Благодаря боль- шой величине изменения объема матрицы при плавлении амплитуда функции источников повышается, достигая 5.07•1014 Па/с, что значительно больше, чем для идеализированного случая (3.31•1014 Па/с). Область отрицательных значений функции источников также расширяется. При плотности энергии 0.06 Дж/см2 происходит неполное плавление наночастицы, которое проявляется как плато в области максимума на кинетической зависимости температуры на границе матрица - наночастица (рис. 1, г). Плавление металла является более энергозатратным, чем матрицы при меньшем изменении объема. Благодаря этому происходит уменьшение скорости роста функции источников. Затем наблюдается резкий спад функции источников S из-за окончания импульса. Развившееся плавление матрицы приводит к росту амплитуды функции источников, и она превышает амплитуду в идеализированном случае на 53.7 %. Плотности энергии 0.1 Дж/см2 (рис. 1, д) достаточно для полного расплавления наночастицы. Во временной области плавления на функции источников появляется область уменьшения с формированием локального минимума из-за неравенств ΔHLM > ΔHL, δVL > δVLM. Амплитуда функции источников выше, чем для идеальной среды на 33.8 % и достигает 8.86∙1015 Па/с. При дальнейшем повышении плотности энергии до 0.3 Дж/см2 (рис. 1, е) область локального минимума на функции источников, который связан с плавлением наночастицы, сильно сокращается. Сохраняются тенденции повышения амплитуды температуры на границе наночастица - матрица и функции источников, как и расширение области отрицательных значений S, связанных с затвердеванием матрицы. При этом относительное увеличение амплитуды по cравнению с идеализированным случаем практически не изменяется, составляя 33.6 %. Примеры рассчитанных оптоакустических сигналов приведены на рис. 2. Значение давления нормировано на величину амплитуды оптоакустического сигнала, возникающего при воздействии на образец с теми же теплофизическими и оптическими характеристиками и 100 %-м мгновенным термическим выходом [1]: , (19) где α(ki) - коэффициент, учитывающий длительность импульса, равный 1 для ki→∞, при длительности импульса на полувысоте 14 нс α = 0.03615. Данный способ представления результатов позволяет сравнивать зависимости, рассчитанные при различной величине плотности энергии. Рис. 2. Рассчитанные кинетические зависимости давления вдали от поверхности (оптоакустические сигналы) при различных значениях плотности энергии на поверхности образца: а - 0.01 (кр. 1), 0.015 (кр. 2), 0.04 (кр. 3) и 0.08 Дж/см2 (кр. 4); б - 0.1 (кр. 1), 0.15 (кр. 2), 0.3 (кр. 3) и 0.5 Дж/см2 (кр. 4) При увеличении плотности энергии оптоакустические сигналы становятся все менее симметричными относительно точки пересечения зависимости оси абсцисс. Зависимость, рассчитанная при плотности энергии 10 мДж/см2, почти симметрична. Плавление затрагивает только нагревание наночастиц в приповерхностной области и сказывается в виде небольших отклонений от экспоненциальной зависимости вблизи максимума. Увеличение плотности энергии с повышением толщин проплавленного слоя матрицы вокруг наночастиц приводит к изменениям оптоакустических сигналов. Хорошо видно увеличение амплитуды сигнала и скорости его нарастания на переднем фронте. Одновременно глубина минимума в отрицательной области возрастает и появляется второй локальный максимум (рис. 2, а). Чем выше плотность энергии, тем дальше он сдвигается вдоль оси времени. При высоких значениях плотности энергии (рис. 2, б) передний фронт сигнала перестает изменяться, глубина минимума относительного давления начинает уменьшаться, амплитуда локального максимума вначале увеличивается, а затем начинает уменьшаться. Наблюдаемые изменения связаны с трансформацией функции источников по мере роста плотности энергии. Особенностью данной функции при нагреве наночастиц является наличие области отрицательных значений, связанных с затвердеванием материала матрицы. По мере роста плотности энергии ширина данной области увеличивается, а минимальное значение становится меньше. Поскольку мы рассматриваем формирование оптоакустического сигнала в образце со свободной передней границей, в решении присутствуют две волны: сжатия и растяжения (см. уравнение (18)). В волне сжатия, которая приходит раньше в точку наблюдения, отрицательная область функции источников приводит к уменьшению давления, в волне растяжения - к росту. Если минимум в волне сжатия, связанный с функцией источников, совпадает с минимумом в волне растяжения, происходит углубление минимума. Дополнительный максимум на сигнале возникает из-за того, что отрицательная функция источников в волне растяжения приводит к росту давления в точке наблюдения. По мере роста плотности энергии минимум на функции источников сдвигается в область больших значений времени, что выражается в соответствующем сдвиге локального максимума. В качестве характерных параметров сигнала использовались амплитуда максимума и эффективная константа роста на переднем фронте, которая определялась при аппроксимации сигнала экспоненциальной функцией. Рассчитанные зависимости данных параметров от плотности энергии приведены на рис. 3. Амплитуда максимума демонстрирует тенденцию к возрастанию со стремлением к стационарному значению в пределе высокой плотности энергии (рис. 3, а). В пределе малой плотности энергии относительная амплитуда меньше единицы из-за меньшего коэффициента термического расширения металла по сравнению с матрицей, что приводит к уменьшению амплитуды функции источников. Локальный минимум в области 0.105 Дж/см2 связан с плавлением металла, которое, ингибируя плавление матрицы, приводит к уменьшению максимума функции источников (рис. 1, д). Рис. 3. Рассчитанные зависимости амплитуды (а) и характерной константы роста (б) оптоакустического сигнала от плотности энергии при значениях относительного изменения объема матрицы в точке плавления 15 (кр. 1), 10 (кр. 2) и 5 % (кр. 3) Зависимость эффективной константы роста сигнала от плотности энергии является немонотонной (рис. 3, б). При малых значениях плотности энергии эффективная константа роста сигнала совпадает со значением для идеализированной системы. Затем в области 0.015 Дж/см2 наблюдается минимум сигнала. Осциллограмма сигнала, соответствующего области минимума, приведена на рис. 2, а. Благодаря плавлению в области, близкой к поверхности, начинается повышение сигнала в области максимума в этом случае. На кинетической зависимости давления возникает область более медленного изменения сигнала. Она появляется на границе между частями сигнала, связанными с областью, где началось плавление и где еще не началось - на волне сжатия. При аппроксимации сигнала данная область попадает в район больших значений, что приводит к уменьшению эффективной константы роста сигнала. По мере роста плотности энергии часть сигнала, связанная с областью образца, в которой плавление не началось, отодвигается в сторону меньших значений времени. В области, где началось плавление, амплитуда функции источников выше, чем там, где плавления не было (дальше от поверхности). При этом амплитуда функции источников, деленная на плотность энергии импульса, зависит от плотности энергии немонотонно: сначала растет, затем наблюдается максимум при Н = 27.5 мДж/см2, после чего уменьшается со стремлением к стационарному значению. Благодаря уменьшению амплитуды в области высоких величин плотности энергии, наблюдается уменьшение эффективной константы роста сигнала. Из рис. 3 следует, что повышение величины изменения объема приводит к росту как эффективной константы, так и амплитуды сигнала. Это естественно, поскольку эффект роста объема при плавлении является причиной описанных выше явлений. В частности, при относительном увеличении объема в точке плавления 15 % максимальное увеличение эффективной константы роста сигнала составляет 34 %, а при δVL = 5 % - рост на 12 %. В случае амплитуды сигнала рост составляет 95 и 35 % в тех же условиях. Выполненный цикл расчетов показывает, что эффективная константа роста сигнала, на основе которой в рамках оптоакустического метода рассчитывается показатель поглощения, может значительно изменяться в результате плавления вещества. Поэтому следует аккуратнее относиться к результатам оптоакустических измерений в системах с наночастицами металлов. В случае композитов PETN-Al возрастание эффективной константы роста сигнала в 5 раз при увеличении плотности энергии было зафиксировано в [2], то есть согласие с данными [2] наблюдается на качественном уровне. Для количественного описания эффекта необходимо учесть зависимость оптических свойств наночастиц от их радиуса, длины волны и температуры, а также стимулируемые нагреванием химические реакции экзотермического разложения материала матрицы, что и планируется сделать в дальнейших исследованиях. Заключение Предложена методика моделирования оптоакустических сигналов, инициируемых лазерным излучением в композитах прозрачная матрица - наночастицы металлов, с учетом процессов плавления матрицы и наночастицы. Показано, что плавление матрицы и наночастицы приводит к изменению вида сигнала. Наибольший эффект вносит плавление матрицы, благодаря которому возникают зависимости эффективной константы роста сигнала и его амплитуды от плотности энергии. Методика не учитывает ряда факторов, перечисленных выше, которые будут включены в модель в следующих работах.
Гусев В.Э., Карабутов А.А. Лазерная оптоакустика. - М.: Наука, 1991. - 304 с.
Адуев Б.П., Нурмухаметов Д.Р., Белокуров Г.М. и др. // Опт. и спектр. - 2018. - Т. 124. - № 3. - С. 404-409.
Schmid Th. // Anal. Bioanal. Chem. - 2006. - V. 384. - No. 5. - P. 1071-1086. DOI: 10.1007/s00216-005-3281-6.
Monchalin J.-P., Bertrand L., and Rousset G. // J. Appl. Phys. - 1984. - V. 56. - P. 190. https://DOI.org/10.1063/1.333751.
Xia J., Yao J., and Wang L.V. // Prog. Electromagn. Res. - 2014. - V. 147. - P. 1-22.
El-Brolossy T.A., Abdallah T., Mohamed M.B., et al. // Eur. Phys. J. Special Topics. - 2008. - V. 153. - No. 1. - P. 361-364. https://DOI.org/10.1140/epjst/e2008-00462-0.
Адуев Б.П., Нурмухаметов Д.Р., Звеков А.А., Каленский А.В. // ФГВ. - 2016. - Т. 52. - № 6. - С. 104-110.
Каленский А.В., Звеков А.А., Никитин А.П., Адуев Б.П. // Теплофизика и аэромеханика. - 2016. - Т. 23. - № 2 (98). - С. 271-279.
Knight M.W. , King N.S., Liu L., et al. // ACS Nano. - 2014. - V. 8. - No. 1. - P. 834-840.
Parashar P.K., Sharma R.P., and Komarala V.K. // J. Appl. Phys. - 2016. - V. 120. - P. 143104.
Буркина Р.С., Морозова Е.Ю., Ципилев В.П. // ФГВ. - 2011. - Т. 47. - № 5. - С. 95-105.
Адуев Б.П., Нурмухаметов Д.Р., Белокуров Г.М. и др. // ЖТФ. - 2014. - Т. 84. - Вып. 9. - С. 126-131.
Gupta S.C. The Classical Stefan Problem (Second Edition). Basic Concepts, Modelling and Analysis with Quasi-Analytical Solutions and Methods. - Elsevier, 2018. - 726 p.
Лыков А.В. Теория теплопроводности. - М.: Высшая школа, 1967. - 422 с.
Groot R.D. // J. Comput. Phys. - 2018. https://DOI.org/10.1016/j.jcp.2018.04.051.
Жвавый С.П. // ЖТФ. - 2000. - Т. 70. - Вып. 8. - С. 58-62.
Таблицы физических величин: справочник / под ред. И.К. Кикоина. - М.: Атомиздат, 1975. - 114 с.
Kirk-Othmer. Encyclopedia of Chemical Technology / ed. by Watcher. - Wiley&Sons, 1998. - V. 2. - P. 93.
Олинджер Б., Кейди Г. // Детонация и взрывчатые вещества: сб. / под ред. А.А. Борисова. - М.: Мир, 1981. - С. 203-219.