CC BY
66
Труды Карельского научного центра РАН № 4. 2014. С. 48-61
УДК 519.6:539.2
ОПРЕДЕЛЕНИЕ ПАРАМЕТРОВ ТЕРМОДЕСОРБЦИИ ВОДОРОДА
Ю. В. Заика
Институт прикладных математических исследований Карельского научного центра РАН
Рассматривается дегазация методом термодесорбции (ТДС) образца конструкционного материала, насыщенного водородом. Представлены математические модели ТДС-эксперимента и метод их параметрической идентификации.
Ключевые слова: водородопроницаемость, краевые задачи ТДС-дегазации, параметрическая идентификация.
Yu. V. Zaika. ESTIMATION OF HYDROGEN THERMODESORPTION PARAMETERS
The degassing of a hydrogen presaturated structural material sample by the thermal desorption method is considered. The mathematical models of a TDS-experiment and parametric identification method are presented.
Key words: hydrogen permeability, boundary-value problems of TDS-degassing, parametric identification.
Введение
Интерес к взаимодействию водорода с различными материалами носит многоплановый характер [1-3, 5, 7, 8, 10-12, 18]. Достаточно упомянуть задачи энергетики, защиты металлов от водородной коррозии, проектирования химических реакторов, ракетостроения. В частности, поскольку в термоядерных реакторах (пока в отдаленной перспективе) предполагается использование радиоактивного изотопа водорода — трития, возникает проблема возможных диффузионных утечек трития и его накопления в конструкционных материалах. Основное внимание уделялось проблеме водородной хрупкости металлов, но постепенно акцент смещается в сторону использования полезных свойств водорода. Гидри-
ды позволяют удерживать большое количество этого экологически чистого энергоносителя. С этим связаны перспективы водородных аккумуляторов и двигателей с относительно высоким уровнем безопасности: без высоких давлений и низких температур. На обратимом легировании металлов водородом основаны пластифицирование и термоводородная обработка титановых сплавов. Некоторые частные задачи исследованы в [9, 13-17]. Энтузиасты говорят не только о водородной энергетике, но и о водородной экономике.
Для повышения эффективности экспериментальных исследований, решения прикладных задач и обобщений необходимы математические модели взаимодействия изотопов водорода с конструкционными материалами и методы их параметрической идентификации.
0
Экспериментальный опыт показывает, что лимитирующими являются не только диффузионные процессы внутри металла, но и физикохимические явления на поверхности [1, 2]. Параметры переноса зависят и от технологических особенностей получения партии материала, поэтому вряд ли следует ориентироваться на получение «табличных данных», нужны эффективные алгоритмы обработки экспериментальных кривых. Адсорбция, растворение, диффузия являются предметом теоретических и экспериментальных исследований. Но каждый дополнительный коэффициент приводит к скачку уровня сложности обратной задачи параметрической идентификации. В статье остановимся на методе термодесорбции, учитывая лишь лимитирующие факторы рассматриваемого далее эксперимента.
Диффузия с обратимым захватом и (де) сорбция
Рассмотрим перенос водорода сквозь образец тестируемого металла (пластину толщиной £). Физико-химическую терминологию в дальнейшем используем в минимально необходимом объеме. Кратко говорим о металлической мембране, хотя это может быть многокомпонентный сплав, интерметаллид. Считаем, что нагрев относительно медленный, практически равномерный, так что градиент температуры пренебрежимо мал и диффузионный поток можно считать пропорциональным градиенту концентрации. Концентрация растворенного водорода (в атомарном состоянии) мала и часть его взаимодействует с ловушками. Здесь под ловушками понимаем микродефекты кристаллической структуры, которые могут удерживать водород. В качестве модели диффузии с обратимым захватом внутри мембраны примем систему дифференциальных уравнений
С*(*,х) = Б(Т)схх - а!(Т)с + а2{т)г, (1)
гг(г,х) = <ц(Т)с(Ь,х) - а2{Т)г{г,х), (2)
где 4 — время, (4, х) е = (0,4*) х (0,^); с(£, х) — концентрация диффундирующего водорода (атомарного); г(Ь,х) — концентрация захваченного диффузанта; И — коэффициент диффузии; <ц, <22 — коэффициенты поглощения и выделения атомов Н ловушками. Ограничение емкости ловушек осуществляется заменой произведения а\с на [1 — г/гтах\а\с.
Для определенности полагаем, что величины .О, щ зависят от температуры Т{£) по закону Аррениуса с предэкспоненциальными множителями Ио, аог и энергиями активации Еи, Е\ (Д — универсальная газовая постоянная):
.О = £)0ехр{-£'о/[ДГ(4)]}, щ = ...
Графиками таких зависимостей /(Г) являются б'-образные кривые насыщения.
Известны более детализированные модели переноса. В частности, можно учесть несколько каналов диффузии (транскристаллическую, по границам зерен, вдоль дефектов) с взаимообменом между ними [1]. В контексте прикладных задач, выделяя только лимитирующие факторы, находят компромисс между полнотой модели и реальными возможностями ее параметрической идентификации используемым экспериментальным методом.
Определенные трудности численного моделирования связаны с нелинейными граничными условиями. Пусть поверхность мембраны контактирует с газообразным водородом. С учетом сорбционных процессов краевые условия моделируются следующим образом [1, с. 177-206 (Габис, Компаниец, Курдюмов)]:
с(0, х)=с(х), г(0, х) = г(х), ж € [0, (3)
со = д(т)яо(0> = д(т)яе(*), * е [о, и], (4)
й)(0 = ^(Т)ро(г) - Ь(Т)^(4) + О(т)сх\о, (5)
®(*) = 1Л8(т)ре{ь) - ъ(т)<£(г) - о(т)сх\е, (б)
Ь{Т) = Ъоехр{-Еь/[ЯТ\}, я(Т) = ...
Здесь: со(4) = с(4,0), ц{£) = с{Ь,£) — граничные объемные концентрации диффундирующего атомарного водорода; <7о(£)> <%(£) — концентрации на поверхностях мембраны (х =
0,£); д(Т) — параметр локального равновесия между концентрациями на поверхности и в приповерхностном объеме; // — кинетический коэффициент; й(Т) — коэффициент, отражающий тот факт, что только малая часть «налетающего» водорода окажется в форме атомов на поверхности (будем называть з коэффициентом прилипания, имея, однако, в виду, что он характеризует итог общего процесса физадсорбции-диссоциации-хемосорбции молекулярного газа в атомы на поверхности); ро(4), рс{£) — давление газа (Н2) с соответствующих сторон мембраны; Ь(Т) — коэффициент десорбции.
0
Соотношения (4) означают, что объемные приповерхностные концентрации со(4), ц(Ь) пропорционально «отслеживают» текущие концентрации <7о(£)> Яе{1) на поверхности, скорость растворения относительно велика. В балансе потоков (5), (6) учитывается возможность накопления Н на поверхности.
В модели фигурирует как молекулярный, так и атомарный водород. Для единообразия подсчет будем вести в атомах: [с] = 1/ст3, И = 1/ст2, [Бсх] = [.7] = 1/ст2в (.7 = Ъц2). В терминах кинетической теории газов величина /лр определяет число частиц (в данном случае молекул Н2), соударяющихся с единичной площадкой поверхности в единицу времени. Но за счет безразмерного множителя й удобно в дальнейшем воспринимать цзр как плотность потока атомов, оседающих на поверхности. Это интегральный показатель, без разделения процесса на различные стадии.
Показатели Е в экспонентах называем энергиями активации, хотя они могут быть линейными комбинациями энергий активаций и теплот более элементарных стадий процессов и иметь разный знак. Детально обсуждать границы применимости модели не будем, у нее своя ниша в широком спектре исследований.
Замечание 1. Давления и концентрации считаем относительно малыми. При необходимости можно ввести в рассмотрение зависимость коэффициента диффузии от концентрации. Наиболее распространенный вариант [3]: в уравнении (1) слагаемое _Осжж заменяют на
дх(Одхф,х)), .О = -0(Т,с) = -0*(Т)(1 + ас).
Здесь и в дальнейшем, когда это удобно, в обозначениях частных производных помимо индекса будем использовать символ д. Учет зависимости коэффициента диффузии от концентрации в форме так называемого термодинамического множителя (1 + ас) является одним из возможных вариантов наряду, например, с экспоненциальным представлением Б = Б(Т, с) = -0*(Т) ехр{ас}. Обычно значение а относительно мало, применима техника метода малого параметра. В прямых задачах нетрудно дополнительно учесть неравномерность прогрева пластины, эффекты термодиффузии и переноса тепла частицами.
Замечание 2. Модель (4) быстрого растворения (локального равновесия) на поверхности
получается из более общих соотношений баланса потоков
к-(Т)д0Ц) - к+(Т)с0(1) = -0(Т)сх(1, 0),
к~(Т)деЦ) - к+(Т)се(1) = 0(Т)сх(1,£).
Коэффициенты к~, к+ характеризуют интенсивность процессов растворения в объеме и выхода на поверхность. Если эти процессы в рассматриваемом диапазоне температур существенно быстрее диффузии (Осх ~ 0), то получаем соотношение (4) с д = к~/к+. Если поверхность изотропна (в смысле Е^~ ~ Е^+), то д(Т) слабо зависит от температуры. Формально можно записать аррениусовскую зависимость д(Т) = до ехр{—Ед/[ЯТ}], Ед = Еь~ — Еь+, но «энергия активации» Ед не обязательно положительна. Дополнительно можно учитывать степени заполнения поверхности и насыщенности объема:
к~{Т) [1 - Оо^фс^] *),*(*)-
- к+(Т)[ 1 - тО{Т)сх\х=01
Наличие «порогового» множителя (1 — с/стах) приводит к следующему. Если концентрация с в приповерхностном объеме близка к максимально возможной, то растворение практически прекращается. Величина в(Ь) = Я{^)/Ятах означает степень заполнения поверхности. В уравнениях (5), (6) можно моделировать плотность потока адсорбции атомов Н (диссоциативной хемосорбции водорода на поверхности) выражением цз(Т)р(£)( 1 — в{Ь))2. В диапазоне малых концентраций в -С 1, что согласовано с квадратичной десорбцией, линейностью уравнения диффузии, И ф -О(с). Зависимостью кинетической константы от температуры (/х ~ 1 / \/Т) обычно пренебрегают на фоне экспоненты в з(Т). Формально можно объединить «всю зависимость» от Т, не считая з(Т) аррениусовским параметром.
Замечание 3. Не для всех материалов целесообразно придавать особую роль поверхности. Часто достаточно уравнения (4)—(6) заменить граничными условиями III рода:
МГ)лм(0 - КГ)4№ = ТО(Т)сх|0>, . (7)
Формально (7) получается из (5), (6) при малой (<7 ~ 0) скорости накопления на поверхности. Коэффициент десорбции обозначается
одной буквой, хотя в (7) и (5), (6) это разные величины (включая размерность). С учетом «усреднения» многостадийного процесса коэффициент объемной десорбции b в условии (7) (ОД-модель) является эффективным коэффициентом рекомбинации [5].
Модель нацелена на задачи, в которых необходимо учесть поверхностные процессы: защитные покрытия, или даже стенки, когда для больших промежутков времени желательно писать меньше нулей в граничных условиях. В качестве предостережения отметим, что при чтении экспериментальных работ целесообразно уточнять, как надежно «измеренный» коэффициент на самом деле вычислялся. Множество коэффициентов с одним и тем же названием обычно содержит более одного элемента.
Экспериментальный метод
Термодесорбционная спектрометрия (ТДС) [3]. В камеру с лентой из исследуемого металла или сплава подается водород в газовой фазе при сравнительно большом давлении. Лента нагревается электрическим током с целью увеличения скорости сорбции. После того как образец поглотит достаточное количество водорода (до состояния равновесного насыщения), он быстро охлаждается (отключается ток нагрева). При этом резко падают скорости физико-химических процессов, и значительное количество водорода остается в образце. В режиме последующего вакуумирова-ния камеры лента снова нагревается. Закон нагрева T(t) может варьироваться в широких пределах. С помощью масс-спектрометра измеряется давление молекулярного водорода в вакуумной камере, обусловленное десорбцион-ным потоком с поверхности (плотность потока десорбции обозначим через J(t) = b(t)q2(t)):
p(t) =9i [ J(r) exp{ (r - i)0Q 1} dr. (8) Jo
Здесь и далее для упрощения обозначений примем сокращенную запись (знак тождества трактуется как равенство по определению):
b(t) = b(T(t)), g(t) = g(T(t)), D(t) = D(T(t))...
Для ТДС выполнены условия симметрии:
P{t) =Ро=Ре, q{t) = qo = qe, c0 = Q,
D{t)cx|0 = -D(t)cx\e, c(x) = c(£-x)...
Константа Q\ зависит от площади поверхности ленты S (в\ = S 62), в о и 02 определяются конкретными характеристиками экспериментальной установки, в частности, объемом камеры V и скоростью откачки вакуумной системы v (во = V/v). Выбор модели измерений (8) обусловлен опытом: впрыск порции водорода в камеру (5-импульс) приводит к резкому скачку давления с последующим экспоненциальным затуханием. Уравнение (8) является классическим в теории измерений. Специфику задачи отражает функция J(t). В дифференциальной форме J(t) = (p(t)/6o +p(t))/6i.
Если дождаться равновесного предварительного насыщения, то в начальный момент времени (i = 0), определяемый повторным нагревом, в (3) с(0,ж) = с = const, z(0,x) = z = const, x G [0,4 Значения с, z априори неизвестны. Но в силу уравнений диффузии с обратимым захватом (1), (2) они при t = 0 (Т = Т(0)) связаны соотношением сцс = a,2Z. Время i* окончания эксперимента определим условием p(t) ~ 0, i ^ i*, c(t*,x) = z(t*,x) = 0.
Замечание 4. При медленном нагреве пренебрегают производной давления p(t): J(t) = 9p(t) = p(t) / ($o#i) ■ Коэффициент 0=1/(0o#i) определяется по результатам калибровки прибора. Если точность оценивания в недостаточна, приходится оперировать измерениями лишь в относительных единицах (приближенно известно ТОЛЬКО отношение J{t)/«Лпах) ■
Параметрическая идентификация
(ОБЪЕМНАЯ ДЕСОРБЦИЯ, dj « 0)
Остановимся на модели с эффективным коэффициентом рекомбинации согласно граничному условию (7). Кроме того, считаем обратимый захват и ресорбцию второстепенными факторами (щ ~ 0, fisp ~ 0).
Замечание 5. В ОД-модели граничное условие Dcx(t,0) = b(%(t) в пределе (t —> +0) формально противоречит равномерному распределению с(0, х) = <р(х) = с. На самом деле в процессе вакуумирования перед ТДС-дегазацией функция <р(х) успевает немного «прогнуться» на краях под действием уже активированной десорбции. В алгоритмах идентификации (р(х) будет использоваться только интегрально (в соответствии с обобщенным решением), так что искажением можно пренебречь.
= ± нели-
Уравнение для оценки параметров
Граничные условия _Осж|о нейны. Но по постановке обратной задачи функция J = Ьсд ^ известна как функция времени. Поэтому формально (без учета погрешностей эксперимента и обработки измерений) можем рассмотреть линейную краевую задачу II рода, заменив Ьсц ^ на плотность десорбции J(t). Сделаем замену Ь' = 0(з)с1з и
после преобразований оставим обозначение 4:
дс
дЬ
д2с дх2’
дс
дх
п= 1
П7ТХ ПТТу
х сое —— сое ——,
где С — функция Грина (функция источника). Явное представление со = со (4):
со (0 = с
Jo
г) с2т,
ад = |[1+2Е ехр{-5^8}]-
п=2тп -1
В исходном времени аргумент 4 — т заменится интегралом £о(з)с1з, а вместо 3(т) будет J(т). Выражение со (4) следует подставить в равенство .7(4) = Ь(4)с^(4) (в = со\/б). Функция ,7(4) известна, поэтому получаем семейство уравнений для оценки параметров:
Ф(4; П0, Ев, Ьо, ЕЪ) = 0, 4 е [41,4г] С (0,4*).
Та же зависимость с учетом сх(Ь,£о) = 0 по_ лучается для отрезка [0, £о] (£о = £/2 — центр симметрии). Имеется в виду «эквивалентная»
краевая задача с* = ОсХХ) х е (0, £о), <р(х) = с, Осх(Ь, 0) = Ьсц, сж(*, ^0) = 0.
Формально можно считать, что обратная задача параметрической идентификации существенно упрощена: осталось подобрать константы в явной формуле. Но коэффициент диффузии И входит в выражение под символами интеграла и рада, что делает этот этап нетривиальным. Если иметь в виду итерационную процедуру оценивания, то, ограничиваясь частичной суммой рада, можно вычислять интегралы типа свертки Лп = Jn.it)'-
Щ = В (4).7(4), .7(4) = Ь(4)с0(4) = Ь(4)с^(4).
Можно было бы сразу перейти к безразмерному (характеристическому) показателю времени диффузии Ь' /£2 и нормировкам х/£, с/с, но в дальнейших промежуточных выкладках это не является принципиальным. По постановке ТДС-эксперимента в режиме «быстрое охлаждение — медленный нагрев» имеем .7(0) ~ 0 ИЛИ ПО крайней мере .7(0) «С .7тах.
Решение краевой задачи с = с(4, х):
с = с — [ 3(т){С(х,1',$,т) + С(х,1',£,т)} (1т,
Jo
1 2 ГП27Г2 1
С(М;у,г) = - + ^2^ехр|-^-(г-4)|х
Jn —
У Л(т) ехр|-^^- У йт = 0(п 2)
как решения линеиных начальных задач п2тг2
Jv.it) = —-р-о(г)лп(г) + J{t), ио) = о
(п = 2т, т € М). При численной реализации уравнения следует нормировать на .7тах.
Замечание 6. Образно говоря, рад для функции К (.б) (0 < я < 1) сходится плохо, спасает (почленное) интегрирование. Чтобы при аппроксимации не оперировать частичной суммой с большим числом слагаемых, для малых
з целесообразно использовать другое представление рада [3, с. 179]:
оо 1 00 2
1 + 2 £ ехр {-»V*} = -= £ ехр {- ^ },
п=1 —оо
в = &£т<*- (9) Возникает задача «склейки» двух частичных сумм с небольшим числом слагаемых на рассматриваемом отрезке времени ТДС-пика.
Если равновесная концентрация с известна, а поток регистрируется в относительных единицах (.7(4) /.7тах), то уровень плотности термодесорбции .7тах находится после опыта из материального баланса:
/ 7(т)Ж- => 2 / 7(т)7ш^ссгг = Jo J о
с£ = 2
Определение коэффициента десорбции
Тождество Ф(Р,Во,Еп,Ъо,Еъ) = 0 позволяет составлять любое количество уравнений /(До, Ев, Ьо, Еь) = 0. Оценка значе-
ний Do, Ejj, bo, Ej> по одной экспериментальной кривой возможна, но требует решения системы уравнений численными методами. Поставим в некотором смысле противоположную цель: ценою дополнительной экспериментальной информации свести вычисления к минимуму. Это повысит и точность оценивания.
Будем исходить из того, что техника определения равновесных концентраций в диапазоне с = с(р,Т) ос у/p (закон Сивертса) отработана. Коэффициент диффузии при фиксированной температуре достаточно надежно определяется методом проницаемости, когда за счет большого перепада давлений на входной и выходной сторонах мембраны при относительно высоких температурах и не слишком малых I достигается диффузионный режим проницаемости. Подробный анализ DLR (diffusion limited regime) представлен в [5]. По времени запаздывания to (точка пересечения с осью t асимптоты графика количества Н, проникшего сквозь мембрану) находят D = ^2/ (6io) [3]. При этом не обязательно достижение входной концентрацией уровня со ~ с — главное, чтобы переходные процессы на входе относительно быстро приводили к установлению Co(i) ~ Со = const, t > £ (е -С i*). Неявно предполагается, что этот начальный всплеск co(i) существенно не изменит начальное распределение <р(х) = 0 (хотя в модели входной поток —Dcx(t,0) формально не ограничен) . Варьирование температуры образца позволяет оценить значения Do, Ец.
Гораздо сложнее определять параметры поверхностных процессов, поскольку SLR (surface limited regime) характерен для невысоких температур и давлений. При этом градиент концентрации мал, модель упрощается (в классе обыкновенных дифференциальных уравнений), но резко падает точность измерений. В ТДС-эксперименте (обычно То = Т(0) — комнатная температура) по мере нагрева SLR плавно переходит в DLR. Самое интересное (окрестность пика потока) происходит на этапе активного «соизмеримого» взаимодействия диффузии и десорбции. По этим причинам рассматриваем распределенную модель. Чтобы повысить точность оценивания и упростить математическое обеспечение, считаем, что помимо ТДС-кривой J(t) известны равновесная концентрация с (при заданных условиях насыщения р,Т) и температурная за-
висимость О = -О(Т) (в аррениусовском случае Ио, Ев) - Не обязательно это потребует дополнительных экспериментов. Параметры поверхностных процессов существенно зависят от трудноконтролируемых внешних условий: окислы, примеси, шероховатость. В этом одна из причин разброса оценок. Объемные же параметры с, И «старого» конструкционного материала известны и отражены в справочной литературе. Изменение свойств поверхности может быть и целенаправленным: например, напыление защитного (микро) нанослоя.
По заданному закону нагрева Т = Т(Ь) определим функции
а(г) = ехр 1^- У Я(Т(я)) <1з},
■7п(£) = оГп (Ь) [ J(т)an (т) йт, п = 2т.
Jo
Суммируя достаточное количество Jn, находим концентрацию
со(Ь) &с-2Г13 - 4^_1{,/2(4) + ... + ,72к(Щ,
8{Ь) = /д ,/(т) <1т. Подставляя J = \р + р/6о]/6\ и интегрируя слагаемые с р по частям, приходим к выражению функций через давление р(Ь) без производной. Фиксируем отрезок [£1,^2] С (0,4*), соответствующий пику ТДС-спектра J(T(t)). На начальном отрезке времени [0, измерения малоинформативны и имеется рассогласование краевых условий модели при 4 —> +0 (если мыслить в терминах классических решений краевых задач). Остаточную дегазацию (4 > £2) также не рассматриваем по причине .7 < «Тщах-
Алгоритм идентификации модели сводится к подбору параметров Ьо, Еъ из условия
J(t) = Ъоехр{-Еь/[КГ]}со(г), Ь е [ьм].
Ориентируясь для определенности на значения Еъ в несколько десятков кДж, представим экспоненту в нормированной форме ии{£),
и({) = ехр{—104/[ДГ(4)]} (Еъ = Юг/кДж).
После логарифмирования приходим к линейному по 1пЬо и V соотношению:
Л(4) = 1п {J(t)c^2(t)} = 1пЬо + г/1пи(Ь) =
= \пЬ0- 1/104[ЯГ(*)]-1, Ь € [ЬМ-
На плоскости \Т~1, Л} имеем отрезок прямой, по пересечениям которой с осями координат
находим значения 1пЬо, V (Ьо, Еь) ■ Аппроксимацию данных прямой можно строить и в координатах {Т, АТ}.
Сопряженные уравнения
Рассматриваем эксперимент, когда можно считаться) та с, J(0) и 0 (.7(0)<^С.7тах): насыщение, охлаждение, вакуумирование, медленный нагрев. Обычно Т{£) = Т0 + уЬ (Т(Ьт) = Ттах => у = 0, Ь ^ 4т). В интегральном измерении несущественно рассогласование при
4 —> +0 равномерного насыщения с условием
Всх(Ь,0) = Ъ(Ь)с1(£) (.7 = Ьсц).
Если не ссылаться на теорию обобщенных решений краевых задач, то можно считать, что к началу отсчета времени |(с—со(0))/с | «С 1 и на отрезках [0,е], [£ — £,£], 0 < е <С £, реализуются градиенты концентрации, соответствующие (в классическом смысле) граничным условиям. Запись .7(0) II^ ^1.1'О I ьттач
температура низка (например, комнатная) и десорбция еще слабо активирована (Ь(0) ~ 0).
Соотношение, получаемое преобразованием по частям двойного интеграла от произведения 1р\(% — Осхх], после подстановки решения сопряженного уравнения = —0(Ь)'фхх (см. [4]) имеет следующий вид:
f ip(0,x) dx +
Jo
Ги дс Ги
J [ip(t,Q)+4p(t,£)\D— dt+j cqD
0 = —с
rt
(<9ф
дх
dt.
Подразумеваем ~ 0, т.е. «пробная»
функция г/} если и растет, то значительно медленнее, чем убывает концентрация на этапе остаточной дегазации. Варианты чр = 1, чр = х дают общий материальный баланс. Для функции ip(t,x) = f3(t) ехрсгж с учетом J = Ьсq'.
Х[1 + exp а£] + Y a [exp а£ — 1] = — [exp а£ — 1],
а
c£ = 2j J(r)dr, f3(t) = exp | — a2 JDdr j,
X
= f'pjdt, Y= f'pDViri Jo Jo
-1Jdt.
Параметр а варьируем в пределах сравнимости по порядку а£ ~ 1. Это касается как функции ехр{<7ж}, так и /3(^)- В случае а = \/£
в аргументе экспоненты /3(t) — безразмерное характеристическое «время диффузии». Замена а на —а не меняет функцию /3(t), значения X, Y и уравнение с точностью до множителя. При предельном переходе а —> 0 имеем чр = 1.
Замечание 7. Аналогично предыдущему можно получить линейные по переменным X, Y уравнения для
ф = /3(t) sin ах, ф = P(t) cos ах,
P(t) = exp {a2J D(t) dr j.
Здесь уже функция /3(t) растет, поэтому следует проанализировать условие равенства нулю интеграла от ^(i*, x)c(t*, х) по х G [0,4 В силу симметричности c(t,x) относительно х = £q = £/2 этого можно добиться, например, выбором а£ = 2тт (для sin) и а£ = 7Г (для cos). Но при этом уравнения вырождаются. По сходным мотивам не рассматриваем краевую задачу на отрезке [0, с условием симметрии cx(t,£o) = 0 и ограничение фх{Ь,£$) = 0. Следует воспользоваться приближением c(t*,x), например, параболическим по х выражением c(i*, х) = c(i*, х; D, Ь) (см. далее).
Подытожим вышеизложенное схемой оценки коэффициентов модели D(T), b(T), s(T) (_Do, Ed, bo, Еъ, so, E3). Как минимум нужно два эксперимента с различными температурами насыщения Т, чтобы при известном коэффициенте десорбции Ь(Т) из соотношений
с£ =
= 2 [ J(т) dт, цвр = Ьс2 (Т = Т)
Jo
определялись параметры йо и Е3. Желательно варьировать давление р, чтобы значение с существенно изменилось, и скорость нагрева V. Для оценки величин Т^о, Ец, Ьо, Еь берем по два уравнения (а = а\,а^) вида
f = °^X+°^Y - 1 = 0, с с
к =
ехр а£ -Ь 1 exp а£ — 1
(ф = Р(Ь) ехр ах). Заметим, что функция /3(4) зависит только от Т^о, Ец, а. Параметрами эксперимента являются р, Т, V, То, Ттах, а значения а задаются численно. Во избежание плохой обусловленности системы можно добавить семейство уравнений /£ = 0. Производные выражений X, У по переменной а вычисляются под знаком интеграла. Из-за большого
разброса порядков величин перейдем к комплексам переменных и параметров, например:
* -—- -—- л J / ^ ^ ^
/ = аЫХ + Г- 1 = 0, J=^,b=—^,
с£ £
X = - f 'j(t)j3(t) dt,Y=-f 'pdVb-'jdt, t* Jo t* 7o
.D = Da2t*, al ~ 1. Далее можно следовать методу наименьших квадратов (МНК):
-F(A), До, Ьо, -Еь) = ^ccj/j2 ->• min, а* > 0.
Замечание 8. От предэкспоненты Do в функции /3(t) можно избавиться, выбирая а2 = a2/Do и задавая значения а. При необходимости вместо равномерного распределения <р(х) = с можно использовать уточнение
<р(х) = с — А[х — £о]2к. Интегрирование функции —с(0, х)ф(0, х) приведет к появлению слагаемого, линейного по А. Значение А > 0 определяется материальным балансом дегазации (равновесную концентрацию с считаем известной). Согласование краевых условий D(p\0) = Ь(р2(0) (t = 0) дает дополнительное уравнение связи искомых величин Do, Ed, bo, Ej,.
Уточним алгоритм для варианта Т = Т (D,b = const: после насыщения сбрасываем давление с помощью дополнительной вакуу-мируемой емкости и измеряем поток). В качестве пробной функции фиксируем 1p(t,x) = (3(t) expax, (3(t) = exp{—a2Dt}. Чтобы избежать зависимости подынтегральных функций в X, Y от оцениваемых параметров, удобно положить а2 = a2/Dt*: /3(t) = ехр{—a2t/t*}. Тогда X становится известным числом, а интеграл Y известен с точностью до множителя:
у __ Т I У V _____ -ул — JmaxL*yv, 1 — гг 1 ,
у/Ь
/3 dr.
Уравнение / = 0 для оцениваемых параметров принимает вид
f(x 1, Х2) = аХА\ус{х\)у[х\ + a2Yx 2 — 1 = 0,
Imax л Jmaxt*
(10)
Выбираем два а и, исключая линейно входящую переменную Х2, решаем численно уравнение относительно у = у/х{, Х\ е [х^,х^].
Внесем поправку, возникающую при замене начальной равномерной концентрации <р(х) = с(0,ж) = с на <р(х) = с — А(х — £о)2к, А > 0. Ограничимся к = 1. При вычислении интеграла от —с(0, х)ф(0, х) появится дополнительное слагаемое Ь:
/•е
/ (cф)\t0*dx = — / (р(х),ф(0,х)<1х =
Jo J о
г£ г£
= —с / ^(0, х) dx+L, Ь = А (х—£оdx = Jo J о
А(£2а2 + 8) г п лЛ А£. п лЛ
= „----- [ехр а£ — 1]---^ [ехр а£ + 1].
4(7 а
В левой части уравнения (10) (в функции /) добавится слагаемое
М = „ \£2а2 — 4х£а + 8] =
Аса2 1 1
= [х!О2 - Ах(х{)у/х{ + 8].
В последнем равенстве сделаны указанные выше подстановки _Осж|о = ИА£ = .7(0),
£2 А .7(0)**
.О**’ Х1 .04* ’ с£ ‘
Если точность измерения уровня ,7(0) недостаточна, то преобразуем А2:
.7(0)4*
DAIU А£2 л л , . = — ^А2 = А2(х 1). £с х\с
x(xi) = [expjcr-v/^T} + ^[ехр^сгу^жГ} - 1]
1. Пусть равновесная концентрация начального насыщения с известна. Константа А > 0 определяется материальным балансом:
г£ о /■*»
/ ^p(x)dx= / J(т)dт = S^f,
J О ^0
с£о — 3-1А^§ = й'*. Согласование краевых условий дает Ь = Ь(-О): -О^'(О) = Ь<^2(0),
у/ОМ = \/Ь[с — А£о] => Х2 = х2{х{).
В уравнении (10) / = 0 с добавкой М в левой части и заменой х2(х{) остается одна переменная х\. Начальное приближение х® находим из ОА£ = .7(0). Коэффициент й(Т) определяется начальным насыщением: цвр = Ьс2.
2. Если равновесная концентрация с неизвестна, то материальный баланс дает выражение с = с(А), а согласование краевых условий приводит к функциональной зависимости
Х2 = x2(xi,А). Следовательно, нужны два уравнения при различных значениях а. Помимо исходного семейства уравнений / = 0 целесообразно использовать и соотношения f'a = 0.
Параболическое приближение
Напомним стадии ТДС-эксперимента: насыщение, охлаждение, вакуумирование, нагрев (с(0,х) = <р(х) = с). Сходимость алгоритмов оценивания в нелинейных задачах обычно локальная. Поэтому актуален выбор предварительной аппроксимации. «Куполообразный» характер распределения концентрации c(t, х) известен. Поэтому в качественном плане за первое приближение удобно взять
c(t, х) та c(t, х) = B(t) — A(t)(x — £0)2,
£ = 2£о, В(0) = с. Симметрия выполнена, функция B(t) > 0 аппроксимирует «срединную» концентрацию c(t,£o), A(t) ^ 0, t > 0. Напомним, что в модели с объемной десорбцией и <р(х) = с краевые условия не согласованы при t —> +0. Решение краевой задачи понимается в обобщенном смысле. Значение Л(0) определяем граничным условием, а <р(х) = с означает, что в интегральные соотношения можно подставлять равномерное начальное распределение, пренебрегая «малым прогибом» на краях. Равновесная концентрация с оценивается по итогам дегазации интегрированием плотности J(t).
Из материального баланса и граничного условия имеем:
Ге £3
c£ — 2S(t)= c(t,x) dx = B(t)£ — A(t) — ,
Jo 12
S(t) = f J(t) dr, D(t)cx(t,0) = b(t)(%(t) => Jo
£D(t)A(t) = b(t) [B(t) - A(t)£l]2.
Выразим отсюда функции A(t), B(t) через коэффициенты модели B,b и известные величины с, S(t). В итоге находим аналитическое выражение c(t,x) = c(t,x; D,b). Относительно функции времени \[А получаем квадратное уравнение с корнями разных знаков, выбираем положительный. Далее подставляем выражения А = А(В, Ь), В = B(D, Ъ) в правую часть равенства-определения co(t) = c(t, 0) и оцениваем значения параметров Dq = Dq, Ев = Ed, bo = bo, Еъ = Efj из соотношения
J1/2(t) = b1/2(t)co{t), t G [ti,t2] С (0,i*).
Замечание 9. В качестве J(t) использована экспериментальная плотность десорбции. Если же нацелиться на «независимость» параболической модели, то следует продифференцировать по t соотношение материального баланса с заменой J = Ьсд. Подставив зависимость В = В (А) из граничного условия, получаем дифференциальное уравнение А = f(A). Начальное условие Л(0) находим из соотношения с = В(0) = В(А(0)).
Первое приближение окажется полезным и в ситуации, когда затруднительно варьировать значения Т, р, v. Приведенная выше система уравнений /j = 0 (тр = /3(t) exp ах, а = 01 2) может оказаться плохо обусловленной. Уточним сопряженное уравнение для вариантов ijj = [3(t) sin ах (cos ах):
1—cos а£
Ф + Xsina£ + Ya\casa£ — 11--------—л— = 0,
ас 1
Ф + Х[1 + cos а£] — Ya sin а£ — — sin ai = 0,
а
Ф = J(ipc)\t,dx, P(t) = ехр|сг2у£)(r)dr|.
Вместо неизвестной концентрации c(t*,x) подставляем приближенное выражение c(t*,x] D,b). При этом можно оставить явную формулу по искомым коэффициентам D, Ь, или ограничиться подстановкой оценочных значений D, Ь, полученных решением обратной задачи по квадратичному приближению. Попутно изложен вариант схемы простых итераций. При необходимости, когда регистрируется только относительный масштаб J/Jinaxj следует провести нормировку уравнений на Jmax и соответствующие преобразования. В общем случае при р/в\ «С р/(во, в\) во всех формулах делаем подстановку J = вр и воспринимаем в дополнительной переменной.
Если ограничиться отрезком [0, £q\
(cx(t,£q) = 0), то аналогично получаем
г£о rtо
0 = / (ipc)\udx — с ip\t=0dx+
Концентрацию c(t,£0) аппроксимируем выражением c(t,£o',D,b) или функцией времени c(t, £0; D, b). Для ip = х, ip = f3(t) ехрах первым слагаемым пренебрегаем. В случае роста /3(t) следует вместо неизвестной концентрации
c(t*,x) брать приближение c(t*,x) с переменными D, Ь или фиксированными оценочными значениями D,b. Поскольку такие уравнения / = 0 более грубые, то их целесообразно использовать в качестве стабилизаторов, добавляя к основным слагаемым в целевой функции МНК величину af2 с достаточно малым параметром регуляризации а > 0.
Аналогично используется любая гладкая симметричная относительно х = £q вогнутая функция c(t,x]A,B,...). Помимо полиномов целесообразно использование фундаментального решения уравнения диффузии.
Функции Грина
И УРАВНЕНИЯ РИККАТИ
Перейдем к ТДС-модели с динамическими граничными условиями:
ct(t, х) = D(t)cxx , с(0, х) = ip{x) = ip{£ - х),
со = се = gq(t), q(t) = -b(t)q2(t) + D(t)cx(t, 0).
Выполняются равенства c(t,x) = c(t,£ — х), Cx(t,0) = —cx(t,£). Параметр локального равновесия поверхности и объема д считаем постоянным. Он получается из баланса
k+(T)co4t) - k~(T)q0je(t) = ±D(T)cx\x=oe
в случае изотропности поверхности и относительной малости правой части: д = к~/к+, Ек- та Ек+, с\х=0,е = cote{t) = gqo,e{t). Вследствие симметрии <7о = qe = q, J = bq2 ([Ь] = cm2/s). Считаем <p(x) = с = gq. Равномерное начальное распределение растворенного водорода уже не противоречит граничному условию q = —bq2 + Dcx(t, 0) при t —> +0. В общем случае задача может оказаться жесткой, когда выполнено |<?(0)| = b(0)q2 1.
Изложим несколько подходов к решению обратной задачи параметрической идентификации. Для этого сначала в компактной форме подытожим представленный выше материал для модели с объемной десорбцией, в которой J(t) = b(t)<%f(t) = ±D(t)cx|o,f (И = cm4/s).
Считая потенциально функцию J(t) известной, переходим к линейной краевой задаче. Здесь два аналогичных варианта: фиксируем для уравнения диффузии граничные условия I рода co,e(t) = I{t)/\/Щ), I = л/J, или II рода ^D(t)cx\ote = J{t). Для определенности остановимся на втором варианте. Линейность краевой задачи позволяет найти явное представление решения c(t, х) (а значит, и
интересующих нас концентраций co^(i)) через соответствующую функцию Грина (функцию источников). Следует иметь в виду, что «явность» формулы условна при наличии символа бесконечной суммы. Полученное таким способом формальное выражение co^(i; D, b, J(-)) подставляем в уравнение измерений
I(t) = у/Щсо^г), t е [h,t2] С (0, **).
Здесь уже функция J(t) является конкретной входной информацией. Явное уравнение Fit] D0, Ер, bo, Еь) = 0 (тождество no t) используем для оценки коэффициентов диффузии и десорбции. Заменяем рад частичной суммой, применяем методы моментов (для получения уравнений f(Do, Еа,Ьо, Еь) = 0) и наименьших квадратов ->■ min.
Заметим, что можно расширить возможности для решения обратной задачи, добавив несколько иной вариант линеаризации исходной нелинейной модели: следуя схеме «предиктор-корректор», достаточно представить квадрат граничных концентраций cj^(i)
в форме cote{t)I{t)/y/b(£).
Метод сопряженных уравнений дает дополнительное («независимое») семейство уравнений для оцениваемых параметров. При этом входные данные являются аргументом помехоустойчивых интегральных операторов. Для достаточно тонких мембран, когда основное внимание уделяем моделированию поверхностных процессов и возможно более грубое описание распределения концентрации в объеме, работоспособным оказывается параболическое приближение c(t,x) та c(t,x). При этом переходим от распределенной модели к обыкновенным дифференциальным уравнениям, что существенно упрощает численное решение обратной задачи с точностью, соизмеримой с точностью ТДС-эксперимента. Эффективной представляется комбинация методов оценивания, чтобы избежать плохой обусловленности системы однотипных уравнений.
Адаптируем изложенное для динамических граничных условий. Условие начального насыщения и общий баланс имеют вид
»s(T)p = b(T)[cg-1}2,
c£ + 2q = c[£ + 2g~1] = 2
T*J(r)
Jo
dr.
Для оценки йо, Е3 необходимо как минимум два эксперимента при известных значениях
b, д. Основная проблема — определение коэффициентов D0, ED, b0, Еъ, д.
Применение функции Грина
Использование функции Грина для линейной краевой задачи II рода, как это делалось для модели с объемной десорбцией, остается без принципиальных изменений. Следует только вместо J использовать функцию G:
cx{t, 0) = G(t) = q(t) + J(i).D_1(i), t = t'.
В исходном времени производная q(t) выражается через b(t), J(t), J(t) в силу определения J = bq2. Но дифференцирование экспериментальных данных некорректно с вычислительной точки зрения. Для граничной объемной концентрации co(t) = cg(t) справедливо представление co(t) = g\/b~lJ =
= с
/П7Г\ 2 Мп = 1 ----- 1
dx +
= / (“0с)I*, dx-с 1р{0,х)
J о J о
+ [ (q + J)[ip(t,0) + if>(t,£)]dt Jo
fU ,
+ coD— dt.
Jo dx
+
(И)
Выберем тр^,х) = /3(4) ехрсгж,
/3(£) = ехр{—и2 [ D(т)dт}.
Jo
Тогда первым слагаемым можно пренебречь в силу тр —> 0, с —> 0 с ростом Ь. Слагаемое, содержащее <7, следует проинтегрировать по частям с учетом /3 = —_Осг2/3. Начальная концентрация с = до (д(0) = д) находится с точностью до неизвестного пока параметра г = д£ из материального баланса
с£ + 2<? = с£[г + 2]г-1 = 2 / J(т)dт = Б*.
Jo
После элементарных преобразований соотношение (11) примет вид
/ = а£нХ + \и£к + z\a2Y — S'*
z + a£ z 2
= 0,
f {й(т) + J(T)}K(~f(t,T)) dr,
Jo
ТГ(*> T) = j D(S)dS’
l Г 00
ЯХя) = 1+ 2^exp{-/ins}
^ n= 1
^ = 2^o- Ограничимся частичной суммой ряда и слагаемые с д(т) преобразуем по частям (именно в такой последовательности, иначе появятся расходящиеся рады). В правой части после преобразований заменяем q(t) на выражение у/J(t)/b(t), q(0) = с/д. Интегралы можно подсчитывать как решения обыкновенных дифференциальных уравнений. Далее для оценки величин д, Do, Ed, bo, Ej, применяем методы моментов и наименьших квадратов. Если температура постоянна (Т = Т), то D,b = const, 7(t, т) = (t — t)D.
Сопряженные уравнения
Приведем соотношения метода сопряженных уравнений для варианта <р(х) = с. На решениях уравнения ipt = —Dipxx получаем
Х= f'/3Jdt, Y= f *f$Dy/b~xJ dt, (12)
Jo Jo
к = [expert + 1][expert — 1]_1. Целесообразно уравнение / = 0 нормировать на Jmaxi*:
_ — z гтР —
/ = а£Х + [аЫ + z] Y + S'* = 0,
z Н- 2
Y = -f
t* Jo
max
Л
а D(t) I(t)
^ in ах ^max
/3(t; D)dt.
y/W)Jn
В случае ТДС-эксперимента с температурой T(t) = Т = const можно сделать интегралы не зависящими от оцениваемых параметров D(T), Ь(Т) подстановкой а2 = а2/Dt*. В общем случае удобно фиксировать характерное значение энергии Е* и представлять экспоненту ехр{—E/RT} в форме степени uv(t), где u(t) = ехр{—Е*/КГ}, I/ = Е/Е*, T = T(t).
Как минимум для оценивания пяти параметров Do, Ed, bo, Еъ, д по уравнениям / = 0 при обработке данных одного эксперимента следует взять три значения а, а другого эксперимента — два. Вместо варьирования а лучше дополнительно использовать соотношения ?а = 0. Коэффициент д исключается, если известна равновесная растворимость с. К тому же есть еще параметры so и Bs, для определения которых два эксперимента — минимум (без учета ошибок измерений). Поэтому целесообразно сначала экспериментальные данные обработать по модели с объемной десорбцией и зафиксировать полученные оценки Do,
Ер. Теперь уже в трех экспериментах можно в уравнениях (12) фиксировать только по одному значению а. Имеем также в виду, что при отсутствии информации об абсолютных значениях потока Л(£) и подстановке Л(£) = вр(£) появится дополнительная переменная в.
Интегро-дифференциальное уравнение Риккати
Рассмотрим модель ТДС-дегазации с динамическими граничными условиями:
ct = D(t)cxx , с(0, х) = (р(х), с0 = се = gq(t), q(t) = R(t) - b(t)q2(t) + D(t)cx|0, R(t) = usp,
p(t) = 0i f J(r) exp {(r — i)^1} dr, J = bq2.
Jo
Функцию R = R(t) (ресорбцию) отличаем от универсальной газовой постоянной R по контексту. Начальные данные симметричны (<р(х) = ip(x — £)), д = const. Для определенности (р(х) = с = const. При необходимости можно внести поправку для варианта <р(х) = с — А[х — Ц2, А = Ао > 0, i = 2£q. Константа А при известной равновесной концентрации насыщения с определяется по итогам дегазации. Для обратной задачи параметрической идентификации входными данными являются функции p(t), J(t).
Сделаем замену времени t' = J^Dds (оставив прежнее обозначение t):
ct(t, х) = схх, с(0, х) = с, с0 = се = gq, (13) cx(t, 0) = -cx(t,£) = q(t) + [J - R]D~1. (14)
Считаем q(t) функциональным параметром, a уравнение (14) — дополнительным соотношением к линейной краевой задаче (13). Сделаем замену, приводящую граничные условия к однородным:
c(t, х) = c(t, х) - gq(t), ct(t, х) = схх + fit), f(t) = ~gq(t), с( 0, х) = ф(х) = 0, со,е = 0.
Запишем решение с помощью функции Грина (функции источников):
сЦ,х)= [ £1(3;, £, *)£(£)$; +
./о
+ [ [ Ст1(ж,£,*-т)/(т)<^<*г,
Jо J0
_ ч 2 Г п2 7г2 ч
^(^е,о = 1Еп=1ехр{-^)х
. П7ТХ . П7Г£
Х81П —81П—.
В динамические граничные условия входит производная сх(1,0):
*(*.0) = -^Г 9(гЕ/ехр{^“(г_*)}4гг>
X)' =Х)п=1,з,5...- ПРИ т = 1 РВД расходится, так что подразумевается почленное интегрирование. В исходном времени получаем
Са;(4,0) = сх^,0) = сх^,£) = сх(1,£) =
Окончательно динамическое граничное условие запишется в форме
q{t) = R(t) - b(t)q2(t)
4gD(t)
x
£
2^2 rt
(15)
xE| ^(r)exp {_ J ds}dr-
Аналогично рассматривается вариант квадратичного начального распределения <р(х)
(?(0) = 5“ V(0), Ф(х) = ~А[Х ~ 4]2 + А#о) ■ Полученное уравнение (15) с квадратичной нелинейностью (которая входит и в функцию R(t) при «расшифровке» p(t)) будем классифицировать как интегро-дифференциальное уравнение Риккати нейтрального типа. Уравнение эквивалентно исходной краевой задаче в том смысле, что решение q(t) однозначно определяет решение c(t,x). Аналогия с функционально-дифференциальными уравнениями нейтрального типа [6] связана с тем, что избавиться от производной q в правой части интегрированием по частям нельзя, так как появится расходящийся рад. Для аппроксимации можно ограничиться конечной суммой и проинтегрировать по частям. Нас интересует
отрезок [ii,i2] С (0,i*), соответствующий выраженному пику термодесорбции (измерения при t та 0,i* малоинформативны). Чтобы не оперировать большим числом слагаемых, целесообразна «склейка» частичных сумм двух разложений (9). При этом от интегралов можно избавиться введением новых переменных, переходя к системе уравнений. Для мощной вакуумной системы ресорбцией пренебрегают (R(t) = 0), что упрощает уравнение.
Параболическое приближение. Для тонких мембран в качестве приближения распределения объемной концентрации примем
с та c(t, х) = B(t) — A(t)(x — £0f,
B(0) = с; <р(х) = с => -А(О) = 0; ср(х) = с — А0(х — £0)2 => Л(0) = Aq.
2£о = I. Связь поверхностной и объемной концентраций в форме условия быстрого растворения gq(t) = c(t, 0) и материальный баланс
ftО ftО ft
<70+ / <р(х) dx = q(t)+ / c(t,x)dx+ / J(r)dr
Jo Jo Jo
дают представление функций A(t), B(t) через функцию q(t). Подставляя далее вместо cx(t, 0) выражение cx(t, 0) в исходное динамическое граничное условие, получаем искомую аппроксимацию — обыкновенное дифференциальное уравнение типа Риккати для поверхностной концентрации q(t).
Для обратной задачи оценки параметров модели ТДС-дегазации имеем скалярное уравнение для поверхностной концентрации q(t) и уравнение измерений I(t) = \fb(t)q(t) («фазовое состояние» известно с точностью до неизвестных априори параметров десорбции Ьо,Еъ). При использовании метода моментов интегрируем произведение M(t)q(t) по частям с последующей подстановкой выражения q(t) = 1{Ь)/у/Щ) (или q = у/вр/y/b).
Оценка поверхностных параметров
Пусть в ТДС-эксперименте известна равновесная растворимость (с = k(T)y/p, Т = Т — температура насыщения); температурная зависимость D(T) (Do, Ed)] начальное насыщение равномерное (<р(х) = с); нагрев относительно медленный и практически равномерный (T(t) = To + vt). Поверхность конструкционного материала модифицирована (примеси,
окислы, напыление, шероховатость). Требуется оценить коэффициенты прилипания s, десорбции Ъ (J = bq2) и локального равновесия д. С точки зрения обратной задачи нас интересуют численные зависимости s(T), Ъ(Т), д(Т). Их физическая интерпретация зависит от выделенных стадий взаимодействия водорода с твердым телом. Считаем s,b,g обобщенными интегральными показателями переноса — эффективными коэффициентами абсорбции, рекомбинации и растворения. В изотропном приближении д = const, Т G [Т~,Т+].
В рамках модели динамического взаимодействия поверхности и объема
fisp = bq2 = bg~2c2 => с = gb^^y/JIs у/p
(Т = Т). Зная зависимости к(Т), Ъ(Т) и значение д, находим s(T). Коэффициент д оценивается материальным балансом дегазации
I+2q = ci[\ + 2(gi)-1] = 2 /*J(r)dr,
Jo
что может потребовать высокой точности определения J(t), или непосредственно по известным данным с, q (например, q соответствует монослою). Остается оценить коэффициент десорбции Ь(Т). Ресорбцией пренебрегаем.
Ограничимся в уравнении (15) первым слагаемым суммы (во втором при п = 3 уже на порядок меняется аргумент экспоненты):
Г 7Г2 /"* 1
a(t) ее expj-^ j^D(s)dsj.
Домножая уравнение на a(t), получаем линейное уравнение
7 = 7 - 7(0 = Я(т)а(т) dr.
Интегрируя с учетом 7(0) = 0, последовательно находим функции времени
7(*)> 7(*)> Я{р) = a-1(i)i(t),
q(t) = - J Q!_1(r)7(r)dr = f(t) (q(u) та 0).
Подставляя выражение q(t) = f(t) в соотношение J(t) = b(t)q2(t), имеем
A(t) = \n{j(t)f~2(t)} =
= 1пЬ0 - I/104[JRT(i)]-1(l04I/ = Eb),t€ [ti,t2].
В координатах {Т-1, Л} ({Г, ТА}) это отрезок прямой, по пересечениям которой с осями координат вычисляются оценки значений Ьо, Е^.
Укажем следующий шаг итерационной процедуры. Фиксируя функцию 71(4) = 7(^)1 учтем второе слагаемое в (п2 = 9):
£a(t)
4gD(t) а-п*,
(t) fq{T)an\r)dT-J{t).
Jo
Обозначая интеграл через 72 (i), приходим к линейному уравнению 72 = —4gD£~lryi — ... и соответствующему уточнению q(t) = /2(4)-
Литература
1. Взаимодействие водорода с металлами / Ред. А. П. Захаров. М.: Наука, 1987.
2. Водород в металлах/ Ред. Г. Алефельд и В. Фелькль. М.: Мир, 1981.
3. Кунин Л. Л., Головин А. И., Суровой Ю. И., Хохрин В. М. Проблемы дегазации металлов. М.: Наука, 1972.
4. Марчук Г. И. Сопряженные уравнения и анализ сложных систем. М.: Наука, 1992.
5. Писарев А. А., Цветков И. В., Марен-ков Е. Д., Ярко С. С. Проницаемость водорода через металлы. М.: МИФИ, 2008.
6. Хейл Дж. Теория функциональнодифференциальных уравнений. М.: Мир, 1984.
7. Indeitsev D. A., Semenov В. N. About a model of structure-phase transfomations under hydrogen influence // Acta Mechanica. 2008. Vol. 195. P. 295-304.
8. Gabis I. E. The method of concentration pulses for studying hydrogen transport in solids // Technical Physics. 1999. Vol. 44, N. 1. P. 90-94.
9. Rodchenkova N. I., Zaika Yu. V. Numerical modelling of hydrogen desorption from cylindrical surface // International Journal of Hydrogen Energy. 2011. Vol. 36, N. 1. P. 1239-1247.
10. Sakintuna B., Lamari-Darkrim F., Hirscher M., Dogan B. Metal hydride materials for solid hydrogen storage: a review// International Journal of Hydrogen Energy. 2007. Vol. 32, N 9. P. 1121-1140.
11. Varin R. A., Czujko T., Wronski Z. S. Nanomaterials for solid state hydrogen storage. Springer, New York, 2009.
12. Vlasov N. M., Fedik I. I. Hydrogen segregation in the area of threefold junctions of grain boundaries // International Journal of Hydrogen Energy. 2002. Vol. 27. P. 921-926.
13. Zaika Yu. V., Bormotova E. P. Parametric identification of hydrogen permeability model by delay times and conjugate equations / / International Journal of Hydrogen Energy. 2011. Vol. 36, N. 1. P. 1295-1305.
14. Zaika Yu. V. Identification of a hydrogen transfer model with dynamical boundary conditions // International Journal of Mathematics and Mathematical Sciences. 2004. N 4. P. 195-216.
15. Zaika Yu. V., Rodchenkova N. I. Boundary-value problem with moving bounds and dynamic boundary conditions: diffusion peak of TDS-spectrum of dehydriding // Applied Mathematical Modelling. Elsevier. 2009. Vol. 33, N 10. P. 3776-3791.
16. Zaika Yu. V., Rodchenkova N. I. Hydrogen-solid boundary-value problems with dynamical conditions on surface // Mathematical Modeling. Nova Sci. Publishers. 2012. P. 269-302.
17. Zaika Yu. V., Rodchenkova N. I. Hydrogen-solid boundary-value problems with free phase transition interface // Advances in Mathematics Research. Vol. 15. Nova Sci. Publishers. 2012. P. 128-180.
18. Zajec B. Hydrogen permeation barrier — Recognition of defective barrier film from transient permeation rate / / International Journal of Hydrogen Energy. 2011. Vol. 36. P. 7353-7361.
СВЕДЕНИЯ ОБ АВТОРЕ:
Заика Юрий Васильевич Zaika, Yury
зав. лаб. моделирования природно- Institute of Applied Mathematical Research,
технических систем, д. ф.-м. н. Karelian Research Centre, Russian Academy of Sciences
Институт прикладных математических 11 Pushkinskaya St., 185910 Petrozavodsk, Russia
исследований Карельского научного центра РАН e-mail: zaika@krc.karelia.ru
ул. Пушкинская, 11, Петрозаводск, tel.: (8142) 766312
Республика Карелия, Россия, 185910
эл. почта: zaika@krc.karelia.ru
тел.: (8142) 766312
