CC BY
6
L Научно-технические ведомости СПбГПУ. Физико-математические науки. 15 (3) 2022
St. Petersburg State Polytechnical University Journal. Physics and Mathematics. 2022 Vol. 15, No. 3 -►
Научная статья УДК 519.6
DOI: https://doi.org/10.18721/JPM.15307
ВОССТАНОВЛЕНИЕ МОДЕЛИ ТЕПЛОВОГО ПРОЦЕССА ПО ПРОСТРАНСТВЕННО-ВРЕМЕННЫМ РАСПРЕДЕЛЕНИЯМ
ТЕМПЕРАТУРЫ Н. Ю. Быков 1 2 '
1 Национальный исследовательский университет ИТМО, Санкт-Петербург, Россия 2 Санкт-Петербургский политехнический университет Петра Великого,
Санкт-Петербург, Россия н nbykov2006@yandex.ru
Аннотация. Метод генеративного дизайна модели (ГДМ) применен для восстановления структуры и коэффициентов дифференциального уравнения в частных производных, описывающего процесс нагрева и испарения мишени лазерным излучением. Исходные синтетические данные включают сценарии нагрева, соответствующие поверхностному или объемному поглощению энергии. Показано, что в случае объемного поглощения для корректного восстановления модели требуется применение процедуры препроцессинга, предусматривающей исключение части исходных данных. Предложена модификация метода, позволяющая учитывать зависимость коэффициентов восстанавливаемого уравнения от температуры. Обсуждается влияние различных статистических критериев, применяемых при селекции оптимального подмножества элементов, на точность восстановления структуры уравнения. Эффективность применения метода ГДМ продемонстрирована для широкого диапазона параметров нагрева мишени и разных вариантов задания энергоподвода. Показана возможность генерации модели по зашумленным данным.
Ключевые слова: метод генеративного дизайна, управляемая данными модель, уравнение теплопроводности, лазерный нагрев
Финансирование: Исследование выполнено при финансовой поддержке Российского научного фонда (грант № 21-11-00296, https://rscf.ru/project/21-11-00296/).
Для цитирования: Быков Н. Ю. Восстановление модели теплового процесса по пространственно-временным распределениям температуры // Научно-технические ведомости СПбГПУ. Физико-математические науки. 2022. Т. 15. № 3. С. 83—99. DOI: https://doi.org/10.18721/ JPM.15307
Статья открытого доступа, распространяемая по лицензии CC BY-NC 4.0 (https:// creativecommons.org/licenses/by-nc/4.0/)
© Быков Н. Ю., 2022. Издатель: Санкт-Петербургский политехнический университет Петра Великого.
^Научно-технические ведомости СПбГПУ. Физико-математические науки. 15 (3) 2022 Original article
DOI: https://doi.org/10.18721/JPM.15307
RECONSTRUCTING THE THERMAL PROCESS MODEL USING THE TIME-SPACE DISTRIBUTIONS OF TEMPERATURE
N. Yu. Bykov 1 2 B
1 ITMO University, St. Petersburg, Russia 2 Peter the Great St. Petersburg Polytechnic University, St. Petersburg, Russia
H nbykov2006@yandex.ru
Abstract. The method of generative model design (GMD) has been applied to reconstruct the structure and coefficients of a partial differential equation describing the target's heating and its evaporation by laser radiation. The initial synthetic data includes heating scenarios corresponding to surface energy absorption or to volume one. It was shown that reconstructing the model correctly required the use of a preprocessing technique providing the exclusion of a part of the initial data if the volume absorption took place. A modification of the method that made it possible to take into account the temperature dependence of the coefficients of the reconstructed equation was put forward. The influence of various statistical criteria used in selecting the optimal subset of elements on the accuracy of reconstructing the equation structure was discussed. The efficiency of the GMD was demonstrated for a wide range of target heating parameters and different options for setting the energy input. The possibility of model generating by noisy data was shown.
Keywords: method of generative model design, data-driven model, heat transfer equation, laser heating and evaporation
Funding: The reported study was funded by Russian Science Foundation (Grant No. 2111-00296).
For citation: Bykov N. Yu., Reconstructing the thermal process model using the time-space distributions of temperature, St. Petersburg Polytechnical State University Journal. Physics and Mathematics. 15 (3) (2022) 83-99. DOI: https://doi.org/10.18721/JPM.15307
This is an open access article under the CC BY-NC 4.0 license (https://creativecommons. org/licenses/by-nc/4.0/)
Введение
В настоящее время модели, управляемые данными, широко используются для предсказания параметров процессов, протекающих в обществе, экономике и природе, с целью осуществления последующего контроля над ними [1, 2]. К таким моделям предъявляется ряд требований, среди которых основными являются точность прогнозов и интерпретируемость самой модели. В случае, если объект характеризуется количественными предикторами, то модель представляет собой, как правило, линейную или нелинейную регрессионную зависимость [3]. Альтернативными выступают модели в виде обыкновенных дифференциальных уравнений (ОДУ) или дифференциальных уравнений в частных производных (ДУЧП), восстановленные по имеющимся данным. Предполагается, что модели в виде ДУ, сгенерированные по данным, будут характеризоваться как хорошей интерпретируемостью, так и точностью соответствующих прогнозов [2].
Применение методов восстановления модели процесса в виде ДУЧП по имеющимся данным актуально для исследования процессов тепломассопереноса. Для существенной части тепловых процессов вид соответствующего ДУ хорошо известен и в простейшем случае представляет собой «классическое» уравнение теплопроводности. Однако, в общем случае, уравнение может включать вторую производную от температуры по времени [4, 5], а для движущейся среды — соответствующее конвективное слагаемое. При наличии энерговыделения внутри объекта, в уравнении появляется также дополнительное слагаемое (слагаемые), определяющее мощность внутренних источников теплоты [5].
© Bykov N. Yu., 2022. Published by Peter the Great St. Petersburg Polytechnic University.
Разработка инструмента отбора значимых членов уравнения из большого массива «строительных блоков» сделает возможной индикацию процессов, протекающих во внутреннем объеме объекта и тем самым недоступных для визуализации в ходе эксперимента. Примерами таких процессов могут служить фазовые переходы и химические реакции. Возможность восстановления конвективного слагаемого с соответствующим весом в уравнении энергии позволит сделать выводы о наличии конвективных процессов в исследуемом объекте. При этом можно получить данные о качественном изменении теплового процесса (например, о переходе от режима преобладающей теплопроводности к режиму развитой конвекции при нестационарном нагреве жидкости [6]).
Решение задачи о восстановлении структуры ДУЧП и определения коэффициентов уравнения применительно к тепловым процессам может быть выполнено с использованием предложенного в статье [7] метода генеративного дизайна модели (ГДМ). Метод предполагает реализацию нескольких этапов.
I. Определение наиболее полной возможной структуры восстанавливаемого уравнения.
II. Дискретизация элементов ДУЧП.
III. Вычисление значений элементов векторов дискретизированных шаблонов по имеющимся данным о пространственно-временных распределениях температуры.
IV. Применение статистических методов для определения оптимальной структуры и коэффициентов ДУЧП.
Несмотря на наличие непосредственно «идеи» метода, без ответа остается ряд вопросов.
Во-первых, результативность метода, очевидно, зависит от качества исходных данных. Помимо указанных выше этапов реализации алгоритма ГДМ, необходимо рассмотреть этап препроцессинга имеющихся синтетических или экспериментальных данных. На указанном этапе следует проанализировать качество таких данных и принять решение об их полном или частичном использовании.
Во-вторых, в литературе подробно не рассматривался алгоритм построения дискрети-зированных шаблонов для случая температурозависимых параметров модели.
В-третьих, на этапе выбора оптимальной структуры уравнения возможно использование различных статистических критериев (критерий Мэллоу, информационные критерии и пр.). Необходим анализ адекватности выбора статистического критерия при определении оптимальной структуры модели.
В-четвертых, одним из самых важных можно считать открытый вопрос о возможности применения алгоритма ГДМ к зашумленным данным.
Конечно, в целом необходимо лучшее понимание возможностей метода ГДМ и наработка практических навыков его использования.
Настоящая работа посвящена развитию метода ГДМ, предложенного в статье [7]. Результаты работы направлены на дальнейшее расширение сферы применения метода, в том числе на получение ответов на поставленные выше вопросы.
Генерация исходных синтетических данных
Для более полной демонстрации возможностей метода ГДМ, в исходные данные для восстановления модели теплового процесса, представленные в данной статье, включена информация о пространственно-временной эволюции температуры в материале, учитывающая как наличие либо отсутствие физических процессов, сопровождающих нагрев, так и температурную зависимость теплофизических параметров среды.
Одним из удобных объектов исследования служит процесс нагрева металлической мишени лазерным излучением. Во-первых, данный процесс зависит от параметров излучения и может сопровождаться фазовыми переходами (плавление материала и испарение его поверхности), во-вторых, диапазон изменения температуры материала очень широк и требует учета температурной зависимости параметров среды. В-третьих, данные могут быть получены синтетически, путем численного моделирования [8 — 10].
В настоящей работе в качестве примера рассматривается нагрев ниобиевой мишени импульсами лазерного излучения умеренной интенсивности. Данный нестационарный процесс может считаться одномерным для случая, когда глубина прогрева мишени существенно меньше диаметра лазерного пятна.
Данные о пространственно-временном распределении температуры Т(х^) в мишени для последующего применения метода ГДМ генерируются путем численного решения уравнения теплопроводности [5, 9, 10] вида
/зт
С р
dT дТ -ю ydt дх у
д . дТ +
ТхХ~дХ+qv' (!)
где х, м, — пространственная координата; t, с, — время; р, кг/м3, с, Дж/(кгК), Х(7), Вт/(м-К), — плотность, теплоемкость и коэффициент теплопроводности материала соответственно; qv, Вт/м3, — объемный источник тепла; ю, м/с, — скорость движения испаряющейся поверхности.
Скорость ю зависит от температуры поверхности Т8 (Т = Т (х = 0)) и есть функция только времени: ю = ю(t). Уравнение (1) записано в движущейся системе координат х, с началом координат, соответствующим поверхности мишени. Для температур поверхности, существенно меньших Тъ (Тъ = 5033 К — температура кипения ниобия при нормальном давлении ръ [11]), испарения поверхности не происходит и скорость движения поверхности ю = 0. При температурах, превышающих Тъ или около нее, наблюдается испарение поверхности. Скорость ю становится отличной от нуля. Коэффициент ю Ф 0 может служить индикатором наличия процесса испарения. Плавление ниобия в рассматриваемой постановке не учитывается (температура плавления ниобия Ть = 2750 К [11]).
Граничное условие на поверхности задается в виде
¡д
дх
=qs- lp®> (2)
x=0
где q , Вт/м2, — поток энергии лазерного излучения через поверхность мишени; L, Дж/кг, — скрытая теплота испарения.
На удаленной границе Т (го, t) = Т, начальная температура соответствует Т (x, 0) = Т0 (Т0 = 300 K).
В данной работе рассматриваются два варианта поглощения лазерной энергии:
(i) задание потока энергии лазерного излучения через поверхность мишени в единицу времени (поверхностный поток) согласно выражению
qs = (1 - Rf )Wo, (3)
где W0, Вт/м2, — плотность потока излучения, падающего на поверхность мишени; Rf — коэффициент отражения поверхности материала; при этом объемное поглощение отсутствует (qv = 0);
(ii) считается, что объемное поглощение qV Ф 0, qs = 0.
В последнем случае подвод тепла на расстоянии x от поверхности определяется выражением
qv (x, t) = а J (хt), (4)
1 ( x, t) =10(t )exP(-aaXX (5)
где I(x,t), Вт/м2, — интенсивность излучения в материале на расстоянии x от поверхности; I0 (t), Вт/м2, — плотность потока излучения, проникающего в материал, I0 (t) = (1 - Rf)W0; aa, м—1, — коэффициент поглощения.
Уравнение (5) выражает закон поглощения Бугера — Ламберта — Бэра. Предполагается, что во время импульса мощность излучения не меняется, т. е.
W0 = const, t < t0; W0 = 0, t > t0. (6)
Необходимые значения оптических параметров взяты из книг [12, 13]: aa = 5-107 м-1 (для длины волны лазерного излучения порядка 1 мкм), R = 0,77. Длительность лазерного импульса t0 = 0,1 мке.
Генерируемые данные соответствуют двум температурным диапазонам:
I. температура мишени существенно ниже температуры кипения Тъ (и температуры плавления ТЬ);
II. температура поверхности выше Тъ.
Для I температурного диапазона экспериментальные данные [14] по теплопроводности ниобия в интервале 300 < Т < 2200 К аппроксимируются полиномом третьей степени:
* = *0 + 6,687 • 10-3 Т + 6,652 • 106 Т2 - 2,256 • 10-9 Т3, (7)
где \ = 51,49 Вт/(м^К).
В представленной постановке плотность и теплоемкость ниобия считаются постоянными и составляют р = 8570 кг/м3, с0 = 263 Дж/(кг-К).
Для II диапазона температурные зависимости параметров жидкого ниобия вблизи его точки плавления в литературе существенно разнятся [14]. В этом случае мы считаем коэффициент теплопроводности и теплоемкость постоянными во всем интервале температур. При этом теплофизические параметры соответствуют параметрам жидкого ниобия [11, 14]:
рЬ = 7580 кг/м3, с0Ь = 449,9 Дж/(кг-К), = 65 Вт/(м-К).
Скорость движения поверхности соответствует закону Герца и предположению о том, что 18 % испаренных атомов рассеиваются в результате столкновений в газовой фазе обратно на поверхность [15]. Плотность насыщенного пара следует закону Клапейрона — Клаузиуса [15, 8]:
ю(Х) = 0,82 ^
Г \1/2
т
2пквТ8 у
ехр
Ьт
kв
1 1
\
,т т У1 ъ 1sУ
(8)
где кв, Дж/К, — постоянная Больцмана; т, кг, — масса атома.
Для численного решения уравнения (1) использована схема Кранка — Николсона метода конечных разностей, обеспечивающая второй порядок аппроксимации как по пространственной координате, так и по времени [5].
Теплопроводность в точке между узлами сетки , и ,+ 1 определяли по средней температуре в узлах:
* ,+1/2=*( (т+1+т )/2).
Для каждого последующего временного шага (п + 1) значения коэффициента теплопроводности, в случае учета ее температурной зависимости, подбирались в ходе итерационного процесса.
Схема расчета задачи с объемным поглощением предполагает наличие объема, в котором происходит поглощение, вокруг каждого узла сетки. Для рассматриваемого одномерного случая границы объема для узла, с координатой х , определяются как
х, - Ах / 2 < х < х, + Ах /2.
С учетом закона Бугера — Ламберта — Бэра (см. равенство (5)), плотность мощности теплового источника в объеме выражается как
Чу (Ху, X) = ^Ар(ехр[-аа (х} - Ах /2)] - ехр[-аа (х} + Ах /2)]}. (9)
При Ах ^ 0 выражение (9) соответствует выражениям (4), (5). С учетом ко нечного значения шага Ах, в основном численном алгоритме реализовано выражение (9). Первый узел внутри тела для постановки задачи с объемным поглощением имеет координату х0 = Ах/2 (поверхность мишени соответствует значению координаты х = 0). Помимо объемного поглощения, необходим учет граничного условия в виде (2) при ч = 0. Для реализации данного граничного условия в схему включен фиктивный граничный узел с координатой х = - Ах/2.
Поглощение энергии лазерного излучения в данном узле отсутствует,
Т (х_1, 0 = т (х0, = Т (£).
Полное время наблюдения за процессом прогрева мишени — 2 мкс.
С целью верификации численного алгоритма, предназначенного для получения данных методом ГДМ, в данной работе выполнено сравнение с имеющимися в монографии [16] аналитическими решениями, полученными при отсутствии процесса испарения поверхности (ю = 0) как для случая объемного поглощения (см. выражение (4)), так и для задания потока энергии через поверхность в виде формулы (3).
Данные, полученные путем численного решения уравнения теплопроводности (1), хорошо согласуются с аналитическими решениями (рис. 1 и 2). Для варианта объемного поглощения энергии излучения, использование уравнения (9) и описанного алгоритма существенно повышает точность расчетов, по сравнению с результатами при использовании выражений (4) и (5) для большого пространственного шага (см. рис. 1). Если же задавать наличие поверхностного потока (рис. 2), то решение аналогично случаю объемного поглощения. Решения различаются только в непосредственной приповерхностной области мишени.
г, к
Рис. 1. Распределения температуры Т по глубине мишени X в момент времени t = 10 нс, при задании объемного поглощения лазерного излучения (вариант расчета 6).
На вставке — начальная область X.
Аналитическое решение [16] (кривая 1) сопоставлено с результатами разных расчетных вариантов: с учетом выражения (9), с шагами 2 нм (кривая 2) и 40 нм (кривая 3); с учетом выражений (4), (5) и с шагом 40 нм (кривая 4)
1«0
Рис. 2. Распределения температуры Т по глубине мишени X в момент времени t = 10 нс, при задании потока лазерного излучения через поверхность (вариант расчета 1).
Аналитическое решение [16] (кривая 1) сопоставлено с расчетным результатом (кривая 2)
Варианты, соответствующие сгенерированным данным, приведены в табл. 1, а примеры данных о распределениях температуры — на рис. 3.
Для второго температурного диапазона результаты расчетов соответствуют варианту 11 (см. табл. 1 и рис. 3,Ь). Для интенсивности излучения лазера 1,2 ТВт/м2 и длительности импульса 100 нс температура поверхности к моменту завершения лазерного импульса превышает температуру кипения и наблюдается процесс поверхностного испарения. К концу временного интервала излучения скорость движения фронта испарения составляет примерно 0,3 м/с.
Таблица 1
Расчетные варианты для генерации исходных данных
Вариант Подвод энергии Временной срез, мкс X Лх,нм Nл
1 qv = 0 0,01 Х0 2,0 973
2 1,00 Х0 11333
3 0,01 (7) 998
4 1,00 (7) 11704
5 2,00 (7) 15962
6 qs = 0 0,01 Х0 775
7 1,00 Х0 11333
8 0,01 (7) 801
9 1,00 (7) 11704
10 1,00 (7) 100 228
Вариант 11: qV = 0, временной срез составляет 0,0975 мкс, /п(0 = 1,2 ТВт/м2, с = сге, X = Хт, Лг = 2 нм, N = 3355
Примечание. Для вариантов 1 — 10 10() = 0,3 ТВт/м2, с = с0.
Обозначения: qV — объемное поглощение лазерного излучения; qs — поток энергии через поверхность мишени; I (?) — плотность потока излучения, проникающего в материал; с — теплоемкость, с0 = 263 Дж/(кг-К), с01 = 449,9 Дж/(кг-К); X — коэффициент теплопроводности, Х0 = 51,49 Вт/(м-К), Х01 = 65 Вт/(м-К), (7) — номер формулы полинома; N — число степеней свободы.
а)
Ь)
Рис. 3. Пространственно-временные распределения температуры для вариантов 4 (а) и 11 (Ь) (см. табл. 1)
Метод генеративного дизайна
Первый этап алгоритма метода генеративного дизайна предполагает определение полного числа возможных элементов восстанавливаемого уравнения (в зависимости от типа задачи). Полный шаблон ДУЧП для рассматриваемой задачи включает конвективное слагаемое и записывается в следующем виде:
дТ 1 д ,ЗТ дТ qv „
— +--X— + ю— + — = 0,
дt ср дх дх дх ср
(10)
где коэффициенты X и ю полагаются неизвестными, тогда как коэффициенты р, с и мощность тепловых источников — известными.
Температурная зависимость коэффициента теплопроводности X полагается известной:
х = во + РТ + р2Т2 + взТ3,
(11)
где в — неизвестные коэффициенты.
Второй этап алгоритма ГДМ предполагает дискретизацию элементов уравнения методами конечных разностей (КР) [7] или конечных элементов (КЭ) [17] и вычисление значений дискретизированных шаблонов по имеющимся данным. В настоящей работе использован метод КР.
Разностный шаблон второй производной по пространству для регулярной сетки имеет вид
%У 1Л+1/2Т+1 -(Х+1/2 +Киг)Т1 + -1/2 Т ]
дх
С учетом зависимости (11) выражение (12) принимает вид
д Т
дх
дх
I К1 в,,
,=0
к
',1
0,5К', П + 0,5 К
i ,п+1
(12)
(13)
(14)
К'
Ах2
^ Т'П I Т'П \
Т з+1 + Т 3
(?и - ТП) +
^ Т'П I 'Т'П \
Т з-1 + Т 3
(Т- - Т")
(15)
В выражениях (13) — (15) и далее индекс ' соответствует пространственному срезу, включающему три узла: 3 — 1,3 и 3 + 1; индекс I — временному срезу, включающему два временных слоя: п и п + 1. Коэффициенты К,п+1 определяются формулами (15) с заменой индекса п индексом п + 1.
Дискретизированное уравнение (10) для пространственного среза ' и временного среза I имеет вид р
IX а р = 0,
р=1
(16)
где Р( = 7 для приповерхностных узлов во время действия импульса, при объемном поглощении и Р( = 6 в остальных случаях; а1 = -1, ар = в /(ср) для 2 <р < 5; а6 = ю, а7 = (ср). Параметры а2 - а6 неизвестны и подлежат определению.
Коэффициенты а ' можно найти по имеющимся синтетическим (или экспериментальным) данным:
грп+1 гр
а! = ^ - Т
3
а11 = К'а-2 для 2 < 5 < 5, (18)
грП+1 грп+1 грп грп
а6 = 0,5 ^ ~ 7-1 + 0,5 7+1 ~ 7-1, (19)
6 2Лт 2Лт
а'Т, = 1. (20)
В векторном виде выражение (16) можно представить как
Е а р = Z 0, (21)
р=1
где Z0 — нулевой вектор; вектор V состоит из элементов а 11 , при этом индекс г изменяется от 2 до N — 1 (пространственный индекс у изменяется от 1 до N N — число пространственных узлов)); индекс 1 изменяется от 1 до Ь (Ь — число временных срезов), Ь < п (п — число шагов по времени):
VI = (уР1уР V3...уР ^ -2)х Ь) = : (а 21а 31а41 ...а ^ - 1}1а22 а32 а42 ...а ^ -1)2 ...а 2Ьа ъЬа4Ь ...а ^ - 1)Ь). (22)
^ Р Р Р Р Р Р Р Р Р Р Р Р ^
Выражение (21) можно представить как
V = а 0 V;
а,
' (23)
Р=2
где У = - а^, а1 = -1, а0 = 0.
Оценки коэффициентов а можно сделать с использованием метода наименьших квадратов (МНК).
Наличие исходных данных позволяет вычислить значения компонент векторов V .
Третий этап алгоритма ГДМ предполагает применение к выражению (23) методов статистического обучения. Последние используют алгоритмы выбора оптимального подмножества элементов и статистические критерии для выбора единственной «правильной» комбинации элементов [3, 18]. В качестве таких критериев можно использовать информационный критерий В1С или С (критерий Мэллоу) [3, 19, 20]. Соответствующие критерии вычисляются как
RSS
В1С = п 1п-+ к 1п п, (24)
п
где п = N — 2) Ь — число наблюдений; К88 — сумма квадратов остатков; к = ре + 2, ре — число элементов, входящих в сумму в правой части выражения (23) (максимально возможное число элементов: Р = Р( - 1),
СР - п + 2(Ре + а (25)
где ^ — средний квадрат остатков после регрессии всего множества предикторов.
Результаты применения метода ГДМ
В общем случае точность восстановления модели в виде ДУЧП с коэффициентами при производных, зависящими от решения, определяется количественным параметром — числом степеней свободы, а также качеством данных. В рассматриваемой задаче под числом степеней свободы подразумевается число узлов (точек) с известным значением температуры. Это число зависит от шага сетки и момента времени, которому соответствуют данные. С течением времени тепло распространяется в глубь материала, увеличивая нагретую область и (соответственно) число узлов для случая равномерной сетки. Под качеством
данных подразумевается несколько факторов: температурный диапазон, который охватывается данными; возможная корреляция данных, а также наличие и вид энергоподвода.
Для случая задания потока тепла через поверхность в виде условия (2) и выражения (3), тепловой баланс элементарного объема, связанного с сеточным узлом], определяется только процессом теплопроводности. Ожидается, что для генерации модели можно будет использовать данные во всех узлах рассматриваемой области.
В случае объемного поглощения такой процесс происходит в приповерхностном слое мишени с характерным масштабом порядка 5 ~ 1/аа = 0,2 107 м (20 нм). Для приповерхностных узлов вклад слагаемого qV в уравнении (1) является существенным в отношении теплового процесса. Величину Osi, которая есть отношение энергии, подводимой к элементарному объему, соответствующему узлу ] с границами
х^ -Ах / 2 < х < х^ + Ах /2
при поглощении излучения, к энергии, подводимой/отводимой от элементарного объема за счет теплопроводности, можно определить как
Os L =
qv (х , t )Ах
q-
Áx/2
' qÁx/2
q-
Áx/2
-X
dT_ dx
Xj —Áx/2
, qÁx/2 -
3T_ dx
(26)
xj + Áx/2 '
Граница «влияния» объемного поглощения энергии на тепловой баланс зависит от распределения модуля локального числа Остроградского |Os |. На рис. 4 в качестве примера приведены данные об изменении |OsJ для случая объемного поглощения и момента времени t = 0,1 мкс. Для точек в приповерхностном слое вблизи координаты x = 7-10"8 |OsL| >> 1 и процесс теплопроводности играет второстепенную роль. Последнее обстоятельство влияет на точность восстановления теплофизических параметров и качество генерации модели уравнения.
В настоящей работе для стадии облучения и случая объемного поглощения приповерхностные узлы (200 узлов) исключались из данных, предназначенных для генеративного дизайна, и сумма в выражении (23) включала полное число предикторов P = Pt — 1 = 5 (P = 6). При учете приповерхностных узлов в исходных данных или при рассмотрении только приповерхностных узлов, для которых выражение (23) включало также слагаемое, связанное с объемным энерговыделением (P = 6, P t = 7), погрешность восстановления коэффициентов превышала 50 % и восстановленная модель включала «лишнее» конвек-
|OsJ
Х, ^m
Рис. 4. Распределение модуля локального числа Остроградского по координате X
тивное слагаемое для первого диапазона температур, меньших температуры плавления Т . Для стадии диффузии тепла после завершения импульса, из данных исключался один приграничный узел для учета специфики задания разностных шаблонов.
Варианты с заданием потока энергии через поверхность оказались не чувствительны к наличию в рассматриваемых данных информации о приповерхностных узлах. Для случая задания поверхностного потока энергии вне зависимости от стадии процесса, из данных исключался один узел, непосредственно прилегающий к поверхности мишени.
В процессе нагрева температура мишени изменяется в широких пределах: Т0 = 300 < Т < 2200 К для первого температурного диапазона и Т0 = 300 < Т < 7000 К для второго. Внутренние тесты показали неэффективность использования данных по точкам в глубине мишени с температурой менее 1,003 Т0 = 301 К. Эти узлы также были исключены из рассмотрения.
Данные, сгенерированные для различных временных срезов, коррелируют друг с другом [3]. В настоящей работе модель уравнения восстанавливается по данным, соответствующим одному временному срезу, состоящему из двух временных слоев. Рассмотрены данные, соответствующие временным срезам как на стадии облучения ^ = 10 и 100 нс), так и после завершения лазерного импульса (1 и 2 мкс).
Число степеней свободы, используемых для генерации модели, составляло от 220 до 16 тыс. и зависело от варианта (см. табл. 1).
На первом этапе исследования метод ГДМ применялся к данным, соответствующим вариантам 1, 2, 6 и 7, предполагающим постоянство коэффициента теплопроводности Х0 и температурный диапазон, не превышающий температуру плавления (см. табл. 1). Для варианта 1 результаты применения алгоритма выбора оптимального подмножества элементов приведены в табл. 2. Число в первом столбце означает количество элементов ре, входящих под знак суммы в уравнении (23). Согласно полученным результатам, минимальные значения критериев В1С и С соответствуют уравнению с одним слагаемым с коэффициентом а2 (помимо слагаемого с коэффициентом а1 = -1) для варианта 1. Для всех вариантов статистическая процедура пакета R [18] определяет значение коэффициента а2 равным 2,2845 10-5 (табл. 3), что соответствует значению коэффициента температуропроводности
а2 = Х0/(рс0) = 2,2845 10-5 м2/с
(см. выражение (7)), которое используется при генерации данных.
Следует обратить внимание, что, помимо точного восстановления коэффициента, так же правильно восстановлена и структура уравнения, которая не должна включать конвективного слагаемого с коэффициентом а6. Таким образом, результаты генерации уравнения с постоянными коэффициентами можно признать хорошими.
Аналогично метод ГДМ применяли к восстановлению уравнения по данным, предполагающим температурную зависимость коэффициента теплопроводности параметров в виде (7). В табл. 2 в качестве примера применения процедуры оптимального подмножества элементов представлены результаты для варианта 3 (случай задания поверхностного потока энергии и момента времени t = 10 нс), соответствующего режиму облучения. Для данного варианта процедура правильно воспроизводит шаблон уравнения без конвективного слагаемого, соответствующего коэффициенту а6. Варианты 4 и 5 (см. табл. 1) соответствуют большим временам процесса: t = 1 и 2 мкс соответственно. Температура мишени снижается, и к моменту времени 2 мкс диапазон ее изменения составляет 300 — 550 К. Число степеней свободы возрастает до 16 тыс. Для больших моментов времени качество восстановления структуры модели практически не меняется и конвективное слагаемое не воспроизводится. Для вариантов с заданием объемного поглощения (варианты 8 и 9) структура искомого уравнения также восстанавливается правильно. Изменение вычислительных параметров — увеличение пространственного шага в 50 раз и пропорционального уменьшения числа степеней свободы в варианте 10 (по отношению к варианту 9) не оказывает влияния на качество восстановления структуры.
Результаты применения алгоритма ГДМ для вариантов 1 — 10 представлены в табл. 3. Здесь приведены нормированные коэффициенты, соответствующие выражению (7) (вы-
Таблица 2
Результаты процедуры оптимального отбора подмножества элементов
для вариантов 1 и 3
Число элементов ре а2 аз а4 а5 аб В1С С р
Вариант 1
1 * - - - - - 29196,57 2,59
2 * - - - * - 29190,57 3,71
3 * * - - * - 29183,85 5,54
4 * * * * - - 29178,95 5,56
5 * * * * * - 29173,63 6,00
Вариант 3
1 * - - - - - 3628,85 2,44-1014
2 * - - - * - 1098,50 1,60-10"
3 * * - - * - 12940,78 2,26-1010
4 * * * * - - 29934,64 4,53
5 * * * * * - 29928,26 6,00
Обозначения: а2 - аб — коэффициенты уравнения (23); В1С, С - статистические критерии.
Примечания. 1. В левом столбце приведено число элементов ре, входящих под знак суммы в уравнении (23). 2. Жирным шрифтом выделены выбранные результаты.
делены жирным шрифтом), и коэффициент а6, связанный со скоростью испарения поверхности. В остальных строках приведены значения восстановленных коэффициентов для указанных вариантов. Значения нормированных коэффициентов в зависимости Т) воспроизводятся с хорошей точностью. Погрешность восстановления коэффициентов слабо зависит от временнуго среза, которому соответствуют данные, и вида энергоподвода, при условии исключения части данных о приповерхностном слое для случая объемного поглощения энергии.
Максимальная погрешность в определении коэффициентов полиномиальной зависимости теплопроводности составляет 0,01 %. Погрешность определяется как
е = |ат / а, -1| -100%,
где ат — коэффициент, для которого наблюдается максимальное расхождение с теоретическим значением аг
Для полного коэффициента теплопроводности погрешность оказывается еще меньше и не превышает 0,002 %. Погрешность оказывается одинаковой для вариантов 9 и 10, несмотря на различие в числе степеней свободы в 50 раз.
Исходные данные для варианта 11 предполагают наличие процесса испарения поверхности мишени, поскольку температура поверхности при соответствующем большом значении интенсивности излучения (см. табл. 1) и времени импульса превышает температуру кипения ниобия. Применение метода ГДМ позволяет правильно восстановить структуру модели, которая для данного случая включает конвективное слагаемое (см. табл. 3). Восстановленное значение коэффициента а6 для варианта 11 отлично от 0.
В табл. 2 приведены значения критериев С и В1С для двух вариантов: 1 и 3. Для этих вариантов минимумы критериев соответствуют одному и тому же набору элементов. Однако для нескольких вариантов критерий Ср, в отличие от В1С, предсказывает непра-
Таблица 3
Результаты применения метода ГДМ
Вариант а0 а2 (х105) а3 (х109) а4 (*1012) а5 (х1015) аб
Сравнение теоретических значений с восстановленными
1, 2, 6, 7 0 2,2845 0 0
1 —7,1669'105 2,2845 0 0
2 1,0482-10"7
6 —3,6375-10—5
7 —8,8505-10—8
Сравнение теоретических значений с восстановленными
3—5, 8—10 0 2,2845 2,9668 2,9513 -1,0009 0
3 2,2415-10—5 2,2845 2,9669 2,9513 —1,0009 0
4 1,2132-107 2,9668
5 —3,1960-10 8 2,9514
8 8,7904-10—6 2,9669 2,9511 — 1,0008
9 4,1593-10—8 2,9668 2,9513 — 1,0010
10 —3,8343-10—7 2,9515
Сравнение теоретических значений с восстановленными
11 0 1,9060 0 -
11 1,4656-10 5 1,9060 0 0,3
Примечания. 1. Значения нормированных коэффициентов, соответствующих выражению (7) (теоретические), выделены жирным шрифтом. 2. Для всех вариантов а1 = -1.
вильную структуру модели. Кроме того, использование В1С приводит к более разреженной модели с меньшим количеством элементов.
Таким образом, можно сделать вывод о предпочтительности критерия В1С для рассматриваемого класса задач.
Возможность метода ГДМ работать с зашумленными данными является крайне важной для реальных приложений. Для демонстрации этой возможности ниже рассмотрен дополнительный вариант расчета, аналогичный по исходным параметрам вариантам 1 и 2 из табл. 1 и соответствующий временному срезу X = 0,0975мкс. Как и ранее (варианты 1 и 2), синтетические данные сгенерированы на сетке с шагом Ах = 2 нм. Эффект шума достигался за счет моделирования значений температуры, согласно закону
Т8 (х, х) = Т (х, х) (1 + 05,), (27)
где Т(х,Х) — значение температуры, найденное из численного решения (1); 0 — случайная величина с равномерным законом распределения из диапазона [—1, 1], 5г — относительная погрешность.
При генерации зашумленных данных величина 5г варьировалась от 10-5 до 10-2 (см. табл. 4).
Для заданного уровня шума метод ГДМ не позволяет корректно восстановить как структуру ДУЧП, так и коэффициенты при производных. Поэтому в случае зашумленных данных потребовалось применение дополнительной процедуры регуляризации, в основе которой лежит учет погрешности задания функции при выборе шага сетки [21].
Таблица 4
Результаты применения метода ГДМ к зашумленным данным
Вариант 5 г а2 е, %
1ш 10-5 2,2854-10-5 0,04
2ш 10-4 2,2912-10-5 0,29
3ш 10-3 2,5112-10-5 9,92
4ш 10-2 3,3634-10-5 47,22
Примечание. Для вариантов 1ш - 4ш аб = 0.
Обозначения: 5г - варьируемая относительная погрешность; а2, аб - коэффициенты уравнения (23); е - погрешность, е = |а2 / а( - 1|100% (а^.- теоретическое значение).
Синтетические данные, как отмечалось выше, сгенерированы с малым шагом по пространственной координате, что дает возможность подобрать необходимый больший шаг, кратный Дх. Используя процедуру определения шага (в соответствии с работой [21]) для указанных в табл. 4 значений относительной погрешности входных данных, мы получили значение Ах = 0,1 мкм. Для заданного значения шага метод генеративного дизайна правильно воспроизводит структуру уравнения без конвективного слагаемого (а6 = 0) для всех рассматриваемых значений 5г, а погрешность восстановления коэффициента температуропроводности
е = |а2 / а( - 1|100%
(а{ = 2,2845 10-5 м2/с - теоретическое значение коэффициента) меняется от 0,04 % для относительной погрешности исходных данных 10-5 до 47 % для 5г = 10-2 (см. табл. 4). Значение относительной погрешности 5г = 10-2 соответствует таковому для абсолютной погрешности определения температуры - около 25 К у поверхности мишени и 3 К у удаленной границы.
В случае невозможности варьирования шага по пространственной координате, для синтетических или экспериментальных данных может быть использована альтернативная процедура регуляризации, основанная на интерполяции температуры с использованием сплайн-функций [21].
Выводы
Использование метода ГДМ для восстановления ДУЧП, описывающих тепловые процессы, открывает широкие возможности. ГДМ можно применять как для непосредственного построения математической модели сложного явления по имеющимся данным, так и для индикации сопутствующих процессов, таких, например, как химические реакции или фазовые превращения, а также для уточнения теплофизических параметров материалов. Однако практика применения ГДМ пока довольно небольшая и предполагает необходимость развития самого метода.
В настоящей статье продолжена работа по развитию алгоритма генеративного дизайна модели применительно к восстановлению модели теплового процесса, которая в общем случае может включать конвективное слагаемое и предполагает зависимость параметров материала от температуры.
Представлены результаты восстановления уравнения в частных производных, описывающего процесс прогрева и испарения металлической мишени импульсом лазерного излучения.
Для генерации исходных данных с целью последующего применения метода ГДМ численно решается уравнение теплопроводности, описывающее нестационарный процесс,
для различных способов учета поглощения энергии лазерного излучения. Предложен способ повышения точности расчета для случая объемного поглощения и ограниченных вычислительных ресурсов, подразумевающих большое значение шага по пространственной координате. Правильность работы вычислительного алгоритма верифицирована путем сравнения с имеющимися аналитическими решениями.
Проведенное исследование показало чувствительность метода ГДМ к типу используемых данных. Для случая, когда исходные данные для восстановления модели соответствуют объемному поглощению энергии материалом и стадии облучения мишени, необходима дополнительная оценка доли данных, которые могут использоваться для генерации модели. Исключение данных о температуре в приповерхностных узлах ведет к более качественному восстановлению структуры модели и минимальной погрешности в определении температурозависимых коэффициентов генерируемого ДУЧП. В случае, если данные соответствуют стадии диффузии тепла после завершения лазерного импульса при объемном поглощении или для любой стадии процесса при задании потока тепла через поверхность мишени, отсутствует проблема с исключением приповерхностных узлов и предварительным препроцессингом данных.
Генерация модели в виде ДУЧП с переменными коэффициентами ведет к необходимости учета существенно большего числа предикторов для применения методов статистического обучения при поиске оптимальной структуры и коэффициентов уравнения. Несмотря на данное обстоятельство, результаты исследования показывают, что метод генеративного дизайна модели, примененный к имеющимся данным, приводит к хорошим результатам в отношении восстановления структуры модели. При температурах, существенно меньших температуры кипения материала, структура уравнения не включает конвективного слагаемого. При превышении температурой поверхности мишени температуры кипения, восстанавливается конвективное слагаемое, связанное с процессом испарения поверхности. Для рассмотренного числа степеней свободы (более 200) погрешность восстановления температурозависимого коэффициента теплопроводности оказывается менее 0,002 %.
Показана возможность применения метода ГДМ для восстановления модели по за-шумленным данным. В этом случае для получения корректной структуры уравнения и значений коэффициентов при производных необходимо использование дополнительных процедур регуляризации.
Дальнейшая практика применения метода генеративного дизайна будет связана с апробацией метода на реальных экспериментальных данных.
СПИСОК ЛИТЕРАТУРЫ
1. Hutter F., Kotthoff L., Vanschoren J. (Eds.) Automated machine learning. Methods, systems, challenges. (The Springer Series on Challenges in Machine Learning). Cham, Switzerland: Springer Nature, 2019. 220 p.
2. Maslyaev M., Hvatov A., Kalyuzhnaya A. V. Partial differential equations discovery with EPDE framework: Application for real and synthetic data // Journal of Computer Science. 2021. Vol. 53. July. P. 101345.
3. Гарет Д., Уиттон Д., Хасти Т., Тибширани Р. Введение в статистическое обучение с примерами на языке R. Пер. с англ. С. Э. Мастицкого. Изд. 2-е, испр. М.: ДМК Пресс, 2017. 456 с.
4. Лыков А. В. Теория теплопроводности. М.: Высшая школа, 1967. 600 c.
5. Самарский А. А., Вабищевич П. Н. Вычислительная теплопередача. М.: Едиториал УРСС, 2003. 784 с.
6. Манухин Б. Г., Гусев М. Е., Кучер Д. А., Чивилихин С. А., Андреева О. В. Оптическая диагностика процесса свободной конвекции жидкости // Оптика и спектроскопия. 2015. Т. 119. № 3. С. 418-423.
7. Быков Н. Ю., Хватов А. А., Калюжная А. В., Бухановский А. В. Метод восстановления моделей тепломассопереноса по пространственно-временным распределениям параметров // Письма в Журнал технической физики. 2021. Т. 47. № 24. С. 9-12.
8. Булгаков А. В., Булгакова Н. М., Бураков И. М., Быков Н. Ю. и др. Синтез наноразмерных материалов при воздействии мощных потоков энергии на вещество. Новосибирск: Институт теплофизики СО РАН, 2009. 462 с.
9. Bykov N. Y., Bulgakova N. M., Bulgakov A. V., Loukianov G. A. Pulsed laser ablation of metals in vacuum: DSMC study versus experiment // Applied Physics А. 2004. Vol. 79. No. 4—6. Pp. 1097-1100.
10. Bykov N. Y., Lukyanov G. A. The direct simulation Monte Carlo of cluster formation processes in laser plume // Proceedings of the 25-th International Symposium on Rarefied Gas Dynamics. Novosibirsk: Publishing House of SB RAS, 2007. Pp. 645-650.
11. Ниобий // Химическая энциклопедия. В 5 тт. (1988-1998). Т. 3. Мед-Пол. Под ред. Н. С. Зефирова, И. Л. Кнунянца, Н. Н. Кулова. М.: Советская энциклопедия, 1992. С. 249.
12. Золотарев В. М., Морозов В. Н., Смирнова Е. В. Оптические постоянные природных и технических сред. Ленинград: Химия, 1984. 216 с.
13. Физические величины. Справочник. Под ред. И. С. Григорьева и Е. З. Мелихова. М.: Энергоатомиздат, 1991. 1230 с.
14. Зиновьев В. Е. Теплофизические свойства металлов при высоких температурах. М.: Металлургия, 1989. 384 с.
15. Анисимов С. И., Имас Я. И., Романов Г. С., Ходыко Ю. В. Действие излучения большой мощности на металлы. М.: Наука, 1970. 272 с.
16. Карслоу Г., Егер Д. Теплопроводность твердых тел. Пер. с англ. М.: Наука, 1964. 488 с.
17. Bykov N., Hvatov A., Kalyuzhnaya A., Boukhanovsky A. A method of generative model design based on irregular data in application to heat transfer problems // Journal of Physics: Conference Series. 2021. Vol. 1959. P. 012012.
18. R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria, 2020. Available online at https://www.R-project.org/.
19. Priestley M. B. Spectral analysis and time series (Probability and Mathematical Statistics). Cambridge, UK: Academic Press, 1981. 706 p.
20. Mallows C. L. Some comments on C // Technometrics. 1973. Vol. 15. No. 4. Pp. 661-675.
21. Ватульян А. О. Коэффициентные обратные задачи механики. М.: Физматлит, 2019. 272 с.
REFERENCES
1. Hutter F., Kotthoff L., Vanschoren J. (Eds.), Automated machine learning. Methods, systems, challenges, (The Springer Series on Challenges in Machine Learning), Springer Nature, Cham, Switzerland, 2019.
2. Maslyaev M., Hvatov A., Kalyuzhnaya A. V., Partial differential equations discovery with EPDE framework: Application for real and synthetic data, J. Comp. Sci. 53 (July) (2021) 101345.
3. James G., Witten D., Hastie T., Tibshirani R., An introduction to statistical learning: with applications in R, Springer, New York, 2013.
4. Lykov A. V., Mikhailov Y. A., Theory of heat and mass transfer, Institute for the Promotion of Teaching Science and Technology (IPST), Thailand, 1963.
5. Samarskii A. A., Vabishchevich P. N., Computational heat transfer, in 2 Vols., Willey, Chichester, 1995.
6. Manukhin B. G., Kucher D. A., Chivilikhin S. A., et al., Optical diagnostics of the process of free liquid convection, Optics and Spectroscopy. 119 (3) (2015) 392-397.
7. Bykov N. Yu., Khvatov A. A., Kalyuzhnaya A. V., Bukhanovskiy A. V., A method for reconstructing models of heat and mass transfer from the spatio-temporal distribution of parameters, Technical Physics Letters. 47 (24) (2021) 9-12 (in Russian).
8. Bulgakov A. V., Bulgakova N. M., Burakov I. M., Bykov N. Y., et al., Nanosized material synthesis by action of high-power energy fluxes on matter, Institute of Thermophysics SB RAS, Novosibirsk, 2009 (in Russian).
9. Bykov N. Y., Bulgakova N. M., Bulgakov A. V., Loukianov G. A., Pulsed laser ablation of metals in vacuum: DSMC study versus experiment, Appl. Phys. А. 79 (4-6) (2004) 1097-1100.
10. Bykov N. Y., Lukyanov G. A., The direct simulation Monte Carlo of cluster formation processes in laser plume, Proc. 25-th Int. Symp. Rarefied Gas Dynamics, Publishing House of SB RAS, Novosibirsk (2007) 645-650.
11. Niobiy [Niobium], In Book: Khimicheskaya entsiklopediya, T. 3. Pod red. N. S. Zefirova i dr. [Encyclopedia of chemistry, Vol. 3, Edited by N. S. Zefirov et al.], Soviet Encyclopedia Publishing, Moscow, 1992. P. 249 (in Russian).
12. Zolotarev V. M., Morozov V. N., Smirnova E. V., Opticheskiye postoyannyye prirodnykh i tekhnicheskikh sred [Handbook of optical constants of natural and technical media], Khimiya Publishing, Leningrad, 1984 (in Russian).
13. Grigoriev I.S., Melikhov E. Z. (Eds.), Handbook of physical quantities, CRC Press, Boca Raton, Florida, USA, 1997.
14. Zinoviev V. E., Teplofizicheskiye svoystva metallov pri vysokikh temperaturakh [Thermal properties of metals at high temperatures], Handbook, Metallurgia, Moscow, 1989 (in Russian).
15. Anisimov S. I., Imas Ya. I., Romanov G. S., Khodyko Yu. V., Deystviye izlucheniya bolshoy moshchnosti na metally [Effect of high-power radiation on metals], Nauka Publishing, Moscow, 1970 (in Russian).
16. Carslaw H. S., Jaeger J. C., Conduction of heat in solids, 2-nd ed., Oxford University Press, Oxford, UK, 1959.
17. Bykov N., Hvatov A., Kalyuzhnaya A., Boukhanovsky A., A method of generative model design based on irregular data in application to heat transfer problems, J. Phys. Conf. Ser. 1959 (2021) 012012.
18. R Core Team. R: A language and environment for statistical computing, R Foundation for Statistical Computing, Vienna, Austria, 2020. Available online at https://www.R-project.org/.
19. Priestley M. B., Spectral analysis and time series (Probability and Mathematical Statistics), Academic Press, Cambridge, UK (1981).
20. Mallows C. L., Some comments on C , Technometrics. 15 (4) (1973) 661—675.
21. Vatulyan A. O., Koeffitsiyentnyye obratnyye zadachi mekhaniki [Inverse coefficient problems in mechanics], Fizmatlit, Moscow, 2019 (in Russian).
СВЕДЕНИЯ ОБ АВТОРЕ
БЫКОВ Николай Юрьевич — доктор физико-математических наук, старший научный сотрудник факультета цифровых трансформаций Научно-исследовательского университета ИТМО; ведущий научный сотрудник, профессор кафедры физики Санкт-Петербургского политехнического университета Петра Великого, Санкт-Петербург, Россия. 195251, Россия, г. Санкт-Петербург, Кронверкский пр., 49А nbykov2006@yandex.ru ОВСГО: 0000-0003-0041-9971
THE AUTHOR
BYKOV Nikolay Yu.
ITMO University, St. Petersburg,
Peter the Great St. Petersburg Polytechnic University
49A, Kronverksky Ave., St. Petersburg, 197101, Russia
nbykov2006@yandex.ru
ORCID: 0000-0003-0041-9971
Статья поступила в редакцию 08.04.2022. Одобрена после рецензирования 06.06.2022. Принята 06.06.2022.
Received 08.04.2022. Approved after reviewing 06.06.2022. Accepted 06.06.2022.
© Санкт-Петербургский политехнический университет Петра Великого, 2022
