CC BY
21
Труды Карельского научного центра РАН № 8. 2017. С. 48-59 DOI: 10.17076/mat666
УДК 519.6:539.2
АЛГОРИТМ ПАРАМЕТРИЧЕСКОЙ ИДЕНТИФИКАЦИИ МОДЕЛИ БЫСТРОЙ ВОДОРОДОПРОНИЦАЕМОСТИ
Ю. В. Заика, Н. И. Родченкова
Институт прикладных математических исследований Карельского научного центра РАН, Петрозаводск
Производство высокочистого водорода необходимо для экологически чистой энергетики и различных химико-технологических процессов. Методом измерения удельной водородопроницаемости исследуются различные сплавы, перспективные для использования в газоразделительных установках. Требуется оценить параметры диффузии и сорбции с тем, чтобы иметь возможность численно моделировать различные сценарии и условия эксплуатации материала (включая экстремальные), выделять лимитирующие факторы. В статье представлены нелинейная модель быстрой водородопроницаемости в соответствии со спецификой эксперимента и алгоритм ее параметрической идентификации.
Ключевые слова: водородопроницаемость; нелинейные краевые задачи; разностные схемы; численное моделирование.
Yu.V. Zaika, N.I. Rodchenkova. PARAMETRIC IDENTIFICATION ALGORITHM FOR A FAST HYDROGEN PERMEABILITY MODEL
High-purity hydrogen is required for clean energy and a variety of chemical technology processes. Various alloys, which may be well-suited for use in gasseparation plants, were investigated by measuring specific hydrogen permeability. One had to estimate the parameters of diffusion and sorption to numerically model the different scenarios and operational conditions of the material usage (including extreme ones), and identify the limiting factors. This paper presents a nonlinear model of fast hydrogen permeability in accordance with the specifics of the experiment and the parametric identification algorithm. This work was supported by the Russian Foundation for Basic Research (Project No. 15-01-00744).
Keywords: hydrogen permeability; nonlinear boundary-value problems; difference schemes; numerical simulation.
Введение
Исследования взаимодействия изотопов водорода с конструкционными материалами инициированы прежде всего задачами энергетики, защиты металлов от водородной коррозии, проектирования химических реакторов [2, 3, 6, 8, 10, 13, 14, 17, 20, 21]. Методом измерения удельной водородопроницаемости исследуются различные сплавы, перспективные для использования в газоразделительных установках. Требуется оценить парамет-
ры диффузии и сорбции с тем, чтобы иметь возможность численно моделировать различные сценарии и условия эксплуатации материала (включая экстремальные), выделять лимитирующие факторы. Некоторые частные задачи водородного материаловедения, созвучные теме данной работы, представлены и исследованы в [12, 15, 23, 27, 28].
Лимитирующими являются не только диффузионные процессы, но и физико-химические явления на поверхности [2, 3]. Параметры пе-
®
реноса зависят и от технологических особенностей получения материала, поэтому вряд ли следует ориентироваться на «табличные значения», нужны эффективные алгоритмы решения обратных задач параметрической идентификации адекватных математических моделей по экспериментальным данным. В статье остановимся на методе проницаемости, учитывая лишь основные факторы для прикладной задачи мембранной фильтрации и информационные возможности рассматриваемого эксперимента. Основой для проведенных математических исследований послужили работа [16] и данные по водородопроницаемости некоторых перспективных сплавов [11].
МОДЕЛЬ ВОДОРОДОПРОНИЦАЕМОСТИ
Распределенная модель переноса
Вначале кратко опишем эксперимент. Образец конструкционного материала, нагретого до фиксированной температуры, является перегородкой вакуумной камеры. Предварительно проведена дегазация. В начальный момент времени на входной стороне скачкообразно создается давление молекулярного водорода. Измеряется падающее давление на входе и растущее давление в выходной емкости.
Рассмотрим перенос водорода сквозь образец материала (пластину толщиной £ и площадью 5). Температура Т постоянна в течение одного эксперимента. Концентрация растворенного водорода (в атомарном состоянии) относительно мала и диффузионный поток можно считать пропорциональным градиенту концентрации. Часть атомов Н взаимодействует с ловушками (микродефекты различной природы, включая микрополости), которые могут удерживать водород. Ограничимся представлением об «ограниченном стоке» без дополнительной детализации. В качестве модели диффузии с ограниченным захватом в объеме примем нелинейную систему уравнений
дс д 2 с
« = D(T > â?- f (T
I = f s «<T >
1 -
z(t, x)
(1)
c - flout(T)z, (2)
поправки и не требует более детального моделирования. Для конструкционных материалов (рассматривается металлический сплав) в рабочем диапазоне температур T е [500, 800] K, как правило, aout ^ а-ш ив процессе насыщения (проницаемости) нет необходимости усложнять модель (далее полагаем aout = 0). Величины D, a зависят от температуры T образца по закону Аррениуса с предэкспонен-циальными множителями Do, ao и энергиями активации Ed , Ea (R — универсальная газовая постоянная): D = D0 exp{ — ED/[RT(t)]}, a = a0exp{ —Ea/[RT(t)]}. Начальные данные: в силу предварительной дегазации
c(0, x) = 0, z(0, x) = 0, x е [0,4 (3)
Из баланса потоков получаем следующие нелинейные граничные условия:
dQu dt
= [2s(T )»(T )po(t) —
c
— b(T )c2(t)] S = —SD-
x=0
dQout dt
= [2s(T)MTMt) —
дс
— b(T )c?(t)]S = SD-
x=£
(4)
(5)
где c(í, х) — концентрация диффундирующего водорода (атомарного); г(Ь,х) — концентрация захваченного диффузанта; Б — коэффициент диффузии; а = а-ш и а0^ — коэффициенты поглощения и высвобождения атомов Н ловушками. Знак тождества часто используем в смысле равенства по определению. Величину zmax считаем малой, захват носит характер
Здесь Qin(t), — количество атомов во-
дорода во входной емкости объема У|п и выходной емкости объема Р0иъ £о(^ = c(t, 0), сц(£) = с(^£). Газообразный водород в рассматриваемом «рабочем» диапазоне температур находится в молекулярной форме, но для единообразия, поскольку сквозь металлическую мембрану диффундирует атомарный водород, подсчет будем вести в атомах. Согласно кинетической теории газов плотность ,1Р падающего на поверхность потока частиц связана с давлением р по формуле Герца-Кнудсена: ,1р = рД/2пткТ (к — постоянная Больцма-на, т — масса молекулы водорода). В контексте рассматриваемой методики эксперимента удобно в качестве единиц измерений выбрать [£] = см, [р] = торр. Тогда численно получаем зависимость = до, ^(Т) ^ 2.474 ■ 1022/уТ ([^] = 1н2/(торрсм2с), [Т] = К. На поверхности происходят процессы физической адсорбции, хемосорбции, диссоциации молекул на атомы, растворения. Лишь малая часть «налетающих» атомов Н окажется в абсорбированном состоянии в объеме. Это отражается множителем 2в. Можно вместо 2в написать в (как параметр модели), но далее удобно в записи 2в безразмерный вероятностный множитель в интерпретировать как долю абсорбируемых атомов Н. Итак, 2вдо — результирующий
z
max
поток атомов в объем сквозь поверхность без разделения на более элементарные стадии. По контексту слово «плотность» опускаем.
Далее, J0/ = bc0 е — это плотности потоков десорбции из образца (отклонение от квадра-тичности существенно лишь при экстремальных температурах), b — коэффициент десорбции. Для s и b также предполагаем арре-ниусовскую зависимость от температуры. По крайней мере формально: в экспоненте «энергия активации» Es может оказаться и отрицательной величиной как линейная комбинация энергий активаций и теплот поверхностных процессов на пути «из газа в раствор». Если с обеих сторон мембраны поддерживать постоянное давление насыщения ps = const (s — saturation) молекулярного водорода при постоянной температуре T, то в итоге установится равновесная концентрация с растворенного атомарного диффузионно подвижного водорода. Из модели (4), (5), приравнивая производные к нулю, получаем с а ^pS:
с = Г/pS, Г = y/2s^/b.
Уточним экспериментальные условия. Объемы Vin,out — несколько литров, толщина мембраны £ меньше мм, площадь S — около см2, давление напуска po(0) — несколько десятков торр. Диапазон [pmin,pmax] невелик, ограничимся zmax = ас, а ^ 0.1. Это не приведет к нарушению закона Сивертса (с + Zmax а y/fiS) , причем С + Zmax ~ С = Ц/pS в пределах экспериментальной точности.
Остается определить величины Qin, Qout. В масштабе времени переноса сквозь мембрану газ находится в термодинамическом квазиравновесии, поэтому воспользуемся формулой N = pV/(kT). Здесь N — количество частиц газа, занимающего объем V при температуре T и давлении p (в системе СИ [p] = Па, [V] = м3, [к] = Дж/К). С учетом соотношений торр = 133.322 Па, Па = Дж/м3 (формально), получаем для соответствующих давлений и объемов в граничных условиях (4), (5) Q = 2N = apV/T, а и 1.9311019. Здесьp, V,T означают численные значения в выбранных единицах (торр, см3, К).
Замечание. Конструктивно мембрана находится в трубке (которая нагрета до заданной температуры) между входной и выходной емкостями. Диаметр трубки достаточно велик, чтобы критерием термодинамического квазиравновесия газа в трубке и емкостях считать равенство давлений. В кинетической константе ц.(T) следует брать температуру мембраны. В объемах Vin,out (массивные оболочки которых находятся при комнатной температуре) газ может нагреваться. Целесообраз-
но в предварительном эксперименте напустить в емкости «комнатный» газ при практически непроницаемой металлической мембране. Затем нагреть трубку и зафиксировать повышения давлений. В рамках приближения (уравнения состояния) идеального газа это позволяет оценивать приращения температуры газа в емкостях. Соответствующие температуры газа и нужно использовать в приведенной для Q формуле (и последующих, за исключением значения ц). Необходимость такого уточнения определяется особенностями конкретной экспериментальной установки. Соответствующая корректировка значений T в дальнейших выкладках не должна вызывать затруднений. К тому же эта процедура относительно слабо влияет на итоговый расчет модельных давлений с учетом погрешностей измерений и при относительно больших объемах V.
Модель быстрой проницаемости
Из физических соображений ясно, что, когда мембрана тонкая и материал обладает высокой водородопроницаемостью, быстро установится квазистационарный режим: ловушки насыщены и распределение диффузанта практически линейное по толщине. В этом смысле расчеты в рамках приведенной выше «общей» модели подтверждают ее адекватность. При отсутствии возможности измерять приповерхностные концентрации на практике часто используется приближение Ричардсона для проникающего потока:
J(t) = -DOxc = Di-1 [ce(t) - c£(t)] »
» Jr(t) = Dri-1 [VMtt) — /¡(i)] .
Поставим задачу моделирования концентраций ce/ по давлениям pe/ (что представляет и самостоятельный интерес) без квазиравновесного упрощения c(t) = Г/p(t).
Спустя малое в масштабе эксперимента время te наступает квазистационарный режим (ловушки насыщены и dx c = — [ce(t) — c^(t)]/i), что существенно упрощает исходную модель (1)-(5) (с учетом формулы Q = apV/T):
Ре/(t) = Тве/Mt) — ce(t)], (6)
t ^ te > 0, ве/ = SD [aVin,outi]-1T,
2s^pe/(t) — bc2,e(t) = ±Di-1[ce(t) — c*(t)]. (7)
Поскольку в силу баланса «вход-выход» выполняются равенства 'p>e(t) = — VinV—t pe(t) ^ pe(t) = pe(te) + VinV—1 [pe(te) — pe(t)], то достаточно из граничных условий (7) выразить концентрации ce/(t) = ce/(pe(t)) и подставить
©
в первое уравнение (6) (выбор знаков соответствует порядку индексов 0,£).
Для численного моделирования удобны безразмерные нормированные переменные
Хо^ ф = 1 + 2£со;^)Ш-1, оо,* (^ = 4£2Г2ро,^)Ь2Б—2 - 1.
(8
При этом система уравнений (7) компактно записывается в симметричном виде ао + 2Х* = X), а* + 2Х0 = Х|. Для переменной Х = Х* получаем неполное уравнение четвертой степени [X2 — а£]2 = 4[2Х + а0], которое решается в радикалах (нас интересует положительный корень из физических соображений). Правда, явное выражение несколько громоздко и все равно придется численно интегрировать первое уравнение (6) вида р0 = /(р0). Поэтому нацелимся на вывод дифференциальных уравнений для Х0,£, поскольку информация о динамике граничных концентраций со,* представляет и самостоятельный интерес.
Продифференцируем по времени уравнения (7) и подставим производные давлений из (6). В переменных Хо/ получим систему
х _V V—1
Х0СО = —2вМ0 [Х0 — Х*] ■ ^ 0"', (9)
Х* (^ = 2вМ* [Х0 — Х* М0 ^ ^Т, М* ^
Х0 Х£ — 1
Х0 — V—1Ко"' ХоХ* — 1
, (10)
а^п
аУО
t ^ tо.
О"'
Замена переменных (8) определяет концентрации с0;£^), по которым из уравнений (7) вычисляются модельные давления ро,^).
Сформулируем поэтапно алгоритм численного моделирования давлений р0,£^) ^ ^ t0) при текущих значениях коэффициентов Б, Ь, в (авторы пользовались свободно распространяемым пакетом ВеПаЪ). Ориентируемся на «обычные» условия эксперимента [11, 16, 18, 19, 30], включая значения р, Т, £, V, Б.
1. Фиксируем t = пропускаем быстрые переходные процессы (десятки секунд в масштабе часов эксперимента). Для переменной Х = Х* выбираем корень полинома четвертой степени [Х2 — а^0)]2 — 4[2Х + а0(^)]. Из физических соображений с^о) > с^0) и, следовательно, Х > 1 + л/Г+а = 1 + 2£Г^р^о)Ш—^
2. Из системы уравнений а0 + 2Х* = Хо2, а* + 2Х0 = Х| (t = t0) находим недостающее значение Хо^о). Достаточно и одного уравнения, но учитываем усредняющие процедуры, включая определение значений р0 *(t0).
3. С полученными начальными данными численно интегрируем систему ОДУ (9), (10). Замена переменных (8) определяет концентрации со/^), по которым из уравнений (7) вычисляются модельные давления ро/^), t ^ ^.
Вычислительные эксперименты показывают, что модельные кривые уже практически неотличимы (при t ^ tо) от тех, которые генерируются исходной принятой моделью — нелинейной распределенной краевой задачей.
Отметим принципиальное отличие от квазиравновесной модели (приближения Ричардсона), где единственным параметром для аппроксимации экспериментальных давлений является комплекс Ф = БГ. При реализации описанного выше алгоритма существенными являются все варьируемые параметры исходной модели, влияющие на проницаемость: Б, Ь, в. Тем самым модель быстрой водородопро-ницаемости не теряет информативность о рассматриваемых параметрах переноса.
Моделирование Н-ПРОНИЦАЕМОСТИ
Комплексирование экспериментов
В экспериментальной практике широко применяются различные модификации методов прорыва и термодесорбционной спектрометрии (ТДС). Результаты измерений зависят не только от особенностей установок, но и от тщательности подготовки образцов новых перспективных материалов, подлежащих тестированию на водородопроницаемость. Нередко последовательная смена методик приводит, например, к появлению дополнительных примесей на поверхности образца, что существенно влияет на воспроизводимость результатов. Эти данные являются входными при решении обратных задач параметрической идентификации, чувствительных к уровню различного рода погрешностей. В этой связи целесообразно комплексирование экспериментов с целью повышения точности и информативности измерений. Предлагается следующий вариант «каскадного» эксперимента.
Мембрана, нагретая до фиксированной температуры, является перегородкой вакуумной камеры. Предварительно проведена дегазация образца. С входной стороны скачкообразно создается достаточно высокое давление газообразного водорода. С выходной стороны производится непрерывная вакуумная откачка и с помощью масс-спектрометра определяется проникающий поток. Это эксперимент прорыва. Его достоинством является надежное определение коэффициента диффузии методом Дайнеса-Бэррера (по так называемому
времени запаздывания). Это позволяет «отделить» в модели объемные процессы от поверхностных, параметры которых значительно труднее оценить. По достижении стационарного уровня проникающего потока повышаем входное давление и дожидаемся установления нового стационарного значения. По трем (как минимум) скачкам на входе фиксируем стационары на выходе, что позволяет оценить «степень прямолинейности» изотермы. Затем прекращаем откачку, и эксперимент переходит в стадию «сообщающихся сосудов», как описано выше. После выравнивания давлений (образец практически равномерно насыщен водородом) можно отключить нагрев, затем произвести откачку с обеих сторон мембраны и начать медленный нагрев образца (ТДС-эксперимент). При этом не происходит разгерметизации диффузионной ячейки и поверхность образца не загрязняется дополнительными примесями. Детали первых двух стадий далее будут уточняться по мере изложения метода решения обратной задачи параметрической идентификации.
Численное моделирование прорыва
Ориентируемся для определенности на достаточно подробные описания условий экспериментов и данные по сплавам с высокой водо-родопроницаемостью на основе металлов 5-й группы, в частности ванадия, представленные в публикациях [1, 5, 11, 16, 18, 19, 22, 30]. Фиксируем Т = 673 К, £ = 0.05 см, р(Т) = 2.474 х 1022/уТ 1н2/(торрсм2с), Б = 2 х 10-5 см2/с, Г = 2 х 1020 1н/(см3 луГорр), Ф = БГ = 4 х 1015 1н/(см с ^торр). Задаем значение в = 1.2 х 10-4 и вычисляем соответствующий коэффициент десорбции Ь = 2ц,в/Г2 = 5.7 х 10-24 см4/с. Задаем тройку входных давлений р)12 3 = {30, 50, 70} торр, которые скачкообразно устанавливаются и поддерживаются до достижения стационарных потоков на выходе. Предварительно мембрана дегазирована, с выходной стороны ведется постоянная вакуумная откачка. Температуру газа во входной и выходной емкостях (достаточно большого объема) считаем равной 300 К. Эта незначительная поправка к комнатной температуре обусловлена нагревом диффузионной ячейки с образцом (уточняется характеристиками конкретной экспериментальной установки).
Условия эксперимента таковы, что на выходе мембраны концентрация практически нулевая, а на входе относительно быстро устанавливается стационарная концентрация (но она ниже равновесной): с < с. В рамках модели
определяем с и ci по формулам
С = л/2^всЬ"Т,
2дврс - Ьс2 = 3 = Бс£-1 ^
^ с = - а + л/а2 + Г2р, а = Б[2Ь£]-1.
Для фиксированных выше модельных значений параметров
С1 = 1.06 х 1021 < С1 = 1.09 х 1021,
с2 = 1.38 х 1021 < с2 = 1.41 х 1021,
сз = 1.64 х 1021 < сз = 1.67 х 1021.
Далее приведем решения стандартных краевых задач с граничными условиями I рода, соответствующих скачкам входного давления.
I этап. Краевая задача метода прорыва:
вг = Бсхх, 0 <х <£, г> 0,
с(г, 0) = с1, с(Ь,£) = 0, í ^ 0,
с(0, х) = 0, х е [0,£].
Решение в форме ряда Фурье:
с(Ь,х) = С1 (£ - х)£-1-
Е^ 2с1 ( п2п2В1 л ппх
— ехЫ--— ^ вт ——.
пп I £2 ) £
п=1
Проникающий поток:
31 (*) = -Бсх\1 = Бс1£-1х
1 + 2 ^(-1)n exp
n=1
n2n2Dt £2
С вычислительной точки зрения целесообразно перейти к безразмерному времени t' = Dt/£2, ориентированному на характерное время диффузии £2/D. При малых t ^ 0 возникает особенность при непосредственном использовании частичной суммы представления
f (t') = 1 + 2 ^(-1)n exp {-n2n2t'}. (11)
n=1
Приведем другое выражение для f, используя аппарат тэта-функций Якоби. Точнее, нас интересует функция
6>3(t,x) = 1+2^exp{-n2n2t} cos(2nnx), t>0.
n=1
При x = 0 имеем другое представление [9] 93(t, 0) = 1 + 2VX exp{-n2n2t} =
Vnt'
El exp - . (12)
X
2
Ряд слева быстро сходится при больших значениях ^ а ряд справа — при малых t (которые нас сейчас и интересуют). Если определить = ^ехр{ — пп2^ (п € Ъ, t > 0), то получим соотношение 0(1Д) = у^^), или
лД^^ехр{—пп^} = У~]ехр{— пта^-1} (п € Ъ),
известное как функциональное уравнение для тэта-функции [7]. После технических преобразований получаем
2
f((/) = £ exp
vnt7
m=2n— 1
-m
(n G N). (13)
Для полноты изложения приведем подробные выкладки, которые можно пропустить. Для разбиения на суммы с четными и нечетными номерами слагаемых удобно ввести обозначение ^/ = ^ 1 з 5 . Преобразования:
f (i') = l + 2^(-l)n exp {—n2Л'} =
n=1
= l — 2 E' exp {—n2n2i'}+
те
+ 2 ^ exp {—4nVi'} (12)
n=1
(12)
-2 V' ■ ■ ■ + —^ V exp
—n
Ii7"
—2 V'... + —^ [l+ exp
n=1
E' l Г V-^'
■ ■ ■ + 2 у exp
+ l + 2 Ё
exp
n=1
n l7
—n
Ii7 —n2
+
(12)
(12) l
n2
7 E' expt — 4?!+
+
At7 l
2 E' exp {-n2n2t'} =
П7! l + 2 X] exp{—n2n2i
n=1
' E' exp {- 1} + 2 fC) -
V^F ^ L 4t' J '2 2' - f(t') = тпр Eexp
2 ^' { n2}
- 4FJ.
Для рассматриваемой прикладной задачи с достаточной точностью можно положить (используя лишь частичные суммы из пяти слагаемых):
l. t' G [0, l0-4] 0-4,
3. t' G [l0-1, l]
2. t' G [l0-4, l0-1] - f (i'): (13)
f (t') = 0; ] - f (t'): " f(t'): (11);
4. ^ ^ 1 ^ /(t/) = 1.
Функция /имеет 5-образный вид кривой насыщения (см. врезку на рисунке 1).
Окончание I этапа: с(£*,х) = с1(£ — х)£—1.
II этап (^ ^ ^ =0, начало отсчета t):
С = Бсхх, 0 < х < £, t > 0,
Ф, 0) = С2, с(М)=0, t ^ 0,
с(0, х) = С1(£ — х)£—1, х € [0,£].
Решение: c(t, х) = с2 (£ — х)£—1 —
2E
n=1
С2 - С1 nn
exp
n2n2Dt ] nnx
Г1* —
Проникающий поток:
J2 (t) = D£-1[C1 + (С2 - С1)х
те 2 2
Г-, «v^/ in« Г n2n2Di] х ; + ^(-l)n exp | - }j
n=1
Окончание II этапа: c(t2,x) = с2(£ - ж)£-1.
III этап. Формулы аналогичны с циклической заменой с1 ^ с2, с2 ^ с3.
Результат склейки этапов в единый «экспериментальный» график проникающего потока J(t) (условно t = (t 1, + t2, ¿2 + t3), J = (J1, J2, J3)) показан на рисунке 1.
Рз=70 торр
и
о 3
тЧ
4 2
1 0
p2=50 торр
0.8 ^ 0.6
¡¡i 0.4 0.2 0 -
0.2 0.4 t, 0.6 0.8 1
100
200 t, с
300
Рис. 1. Установление стационарных потоков
7
0
Пересчет давления в поток
В реальном эксперименте измеряется не проникающий поток, а давление газа в выходном объеме V = VOut. Поэтому приведем соответствующую формулу пересчета. Обозначим v скорость откачки вакуумной системы ([v] = м3/с). За основу примем уравнение состояния идеального газа: pV = NkT. Здесь [p] = Па, [V] = м3, N — число частиц (молекул H2), k — константа Больцмана. Дифференцируя по t, получаем pV = NkT.
Подсчитаем баланс частиц за время At:
AN = J SAt - v At N, Je = -SD^ 2 V dx
где [N/V] = 1h2 /м3 — концентрация, J = 1н/(м2 с). Множитель 1/2 появился из-за того, что диффузионный поток атомарный, а частицей в объеме V является молекула H2.
Поделим уравнение материального баланса на At и устремим At ^ 0. В итоге:
NV(t) = Jjf S - £n (t) ^
^ pV = N kT = kT Ji(t)
J S - — N L 2 V .
= kTS- — kT N (t) ^
p(t)V
^ pp(t) = ^J^t) - ~v'p(t), P(0) = 0.
Обозначая 91 = kTS/(2V), 0o = V/v, получаем p = 01J^ - p/0o (po = 0), откуда
/* t t p(t) = 01J exp {T0--} Jt(T) dr.
Формула пересчета давления в поток:
2V 2v
Jl(t) = T^ P(t) + ^F* p(t).
kTS
kTS
что не требует предварительной аппроксимации pp(t). Составной формулы Симпсона достаточно. Следует не забыть, что в текущем пункте используются единицы измерения [kT] = Дж = Н ■ м, [L] = м, [p] = Па = Н/м2 (Па и 7.5 х 10-3 торр), [J] = 1н/(м2 с). В дальнейшем изложении возвращаемся к принятым в данной статье единицам измерения.
Краевая задача водородопроницаемости
По достижении стационара в методе прорыва прекращаем откачку на выходе и поддержание постоянного давления на входе. Комплексный эксперимент переходит в стадию «сообщающихся сосудов»: давление на входе падает, на выходе растет (измеряются p0/(t)). Пока речь идет о прямой задаче моделирования давлений водорода в объемах Vin,out.
Дополнительные входные данные: S = 0.5см2, Vin = 1500см3, Vout = 2200см3, pe(0) = po = p3, С = Гу^. Здесь ориентируемся на экспериментальные условия и данные по сплаву Vg5Ni15, конкретизированные в [5, 11].
Решаем краевую задачу водородопроница-емости, ловушками пренебрегаем (а = 0):
^ = D(T ) ^ dt D(T ) dx2 '
c(0, x) = (i - x)cs£
-1
x=0,£
дс
2в(Т)КТ)ро,е() - Ь(Т)с20/(1) =
= -[2в»РоА*) - Ьс0/№
Уточним, что в зависимостях Б, Ь, в, у от Т берется температура мембраны, а в выражениях для Qin,out — температура газа в емкостях (учитываем поправку к комнатной температуре вследствие нагрева диффузионной ячейки). Модельные кривые давлений молекулярного водорода представлены на рисунке 2.
Измеряемая функция р(Ь) зашумлена, поэтому сначала применяем процедуру сглаживания, а потом вычисляем производную р(£). Заметим, что при достаточно мощной откачке и относительно медленной проницаемости первое слагаемое играет роль малой поправки к приближению 3 « 2vp(í)/(kTS). По мере асимптотического выхода на стационар р ^ 0 ^ 3 = 2ир/(кТБ).
Для вычисления времени запаздывания (см. далее) необходимо интегрирование:
I Мт) ^ = кгёр(*) + кто Iр(т) ^
70
60
50
а
& 40 н
5 30 â * 20
10
0J-.
0 1 000 2 000 3 000 4 000 5 000 6 000 t, с
Рис. 2. Динамика давлений po/(t)
Достаточно стандартных пакетов прикладных программ, если вместо «полной» модели использовать систему ОДУ (9), (10) (в выражения Mq^ подставляем значения температуры водорода в объемах Vin,out). Для этого следует пропустить начальное время to в пределах нескольких минут, пока установится квазистационарный режим (не квазиравновесный). Затем применяем изложенный выше алгоритм для модели быстрой проницаемости. На рисунке 2 отражено, что существенной дополнительной ошибки численного моделирования не возникает (визуально кривые совпадают). В принципе, дождавшись равномерного насыщения образца (выравнивания давлений), можно перейти к ТДС-эксперименту. В данной статье ограничимся «склейкой» только двух типов экспериментов и ссылками [24-26] (ТДС-модели для цилиндрического образца и порошка представлены в публикациях [4, 27, 29]).
Для тестирования излагаемого далее алгоритма решения обратной задачи параметрической идентификации численно генерировались модельные кривые. Затем порождающие их параметры «забывались».
Общая методика идентификации
Определение времени запаздывания
Для полноты изложения кратко опишем метод оценки коэффициента диффузии, предложенный Дайнесом-Бэррером. Асимптотически график потока Ji(t) выходит на стационар Ji = Dc1/£. Соответственно, Q(t) = /о J(т) dr и Ji[t - £2/(6D)] (t ^ t* = ti). Пересечение асимптоты с осью t дает так называемое время запаздывания tq = £2/(6D), по которому оценивается коэффициент диффузии. Аналитически
rt*
tq и t* — / Ji(r)Ji dr, Ji и Ji, t ^ t*. i
Обратим внимание на то, что под интегралом — относительная величина, не требующая абсолютных значений проникающего потока в каких-либо единицах измерения (Ji = sup Ji(t)). Кроме того, значение т0 не зависит от ci. Обычно считается, что на входе быстро устанавливается практически равновесная концентрация с i = ^ v/Pi, что позволяет по соответствующему значению Ji = Dc i/£ дополнительно оценить растворимость Г = ^Wb и проницаемость Ф = Dr. В данной работе считаем такое предположение слишком оптимистическим. Ослабим его с целью повышения точности дальнейших оценок. Полагаем,
что по условиям эксперимента можно считать быстрым установление стационарной входной концентрации с0(^ « с1 < с 1 (t ^ е, е ^ 1) при соблюдении с^) « 0. Само значение с1 еще предстоит уточнить. Таким образом, на данном этапе достоверной считаем только оценку коэффициента диффузии
При новом отсчете времени (^ ^ t0 = 0), интегрируя выражение 7г(^, получаем
Г[^2(т) — 71] Йт » [72 — 71] ■ ^ — ^2(6^)-1], ■)о
где 7г = ^Сг/^, t ^ ^ = t2. Формально, меняя не только начальное время, но и уровень отсчета потока (превышение значения 71), получаем те же выражения для времени запаздывания и оценки Дополнительной информации (об искомых поверхностных параметрах Ь, в) нет, но тройной прорыв позволяет усред-ненно уточнить Для модельного численного эксперимента т01 « т02 « т03 « 20.82 с.
Изотерма. Начальные приближения Ь, в
В экспериментальной практике часто строят и анализируют изотерму — график зависимости стационарного проникающего потока 7 от входного давления р при вакуумной откачке на выходе. Если ориентироваться на закон Сивертса и (квази)равновесную концентрацию с = Г\/р на входе (соответственно 7 = то естественно строить график 7 = 7(у/р). Проанализируем уравнение стационара:
2вдо — Ьс2 = 7 = ^
^ с = а[ — 1 + V! + а-2Г2р], а = Д2£Ь)-1, 7 = ТОГ1 = ^-1а[ — 1 + V + а-2Г2р ].
Асимптотический анализ показывает, что зависимость 7 (^/р) имеет вид параболы (7 а р) при малых давлениях напуска р :
а-2Г2р = х2 (х а л/Р) ^
^ с = а[ — 1 + 1 + х2] = а[0.5х2 + ...]
и вид прямой при относительно больших р :
а-2Г2р » 1 ^ 7 = —^2(2£26)-1 +
По прямолинейному участку изотермы находим (наклон) и, зная оценку определяем начальное приближение коэффициента растворимости Г. По пересечению прямой с осью ординат находим Ь. Зная значения Г = л/2в^/6 и Ь, вычисляем в и Ф = ^Г.
Графическая иллюстрация (используя минимальный набор вычисленных модельных значений 71,2,3) представлена на рисунке 3.
8 7 6
i—I
- 5
w
> ✓ * / *
/ _ / * / *
1■■■■■■■■■ 1 1■■■■■■■■■ 1 ■ ■■■■■■■ 1 ■ ■■■■■■■ 1 1 ■■■■■■■■■ 1 1 ■■■■■■■■■ 1 1 ■■■■■■■■■ 1 .........
3 4 5 6 Vp, VTCpp
Рис. 3. Экстраполяция изотермы стационара J
Отметим, что найденные начальные приближения хорошо согласуются с исходными «забытыми» параметрами:
«истинное» зн-е приближ. зн-е
D 2.00 x 10-5 2.0012 x 10-5
b 5.72 x 10-24 5.8371 x 10-24
s 1.20 x 10-4 1.2221 x 10-4
Г 2.00 x 1020 1.9984 x 1020
Ф 4.00 x 10i5 3.9993 x 10i5
Финальная стадия алгоритма параметрической идентификации
Приведенная таблица отражает лишь вычислительные погрешности при решении прямой и обратной задач. Реальные экспериментальные данные зашумлены. Алгоритм идентификации использует лишь интегральные операторы обработки измерений, что нацелено на определенную помехоустойчивость. В методе прорыва погрешности измерений значительны, требуются данные о проникающем потоке (а это связано с уточнением характеристик вакуумной системы). Модель скачка концентрации растворенного водорода на входе также не отличается высокой точностью. Приходим к выводу, что стадия «сообщающихся сосудов», когда измеряются давления молекулярного водорода в течение длительного времени, характеризуется значительно более высокой точностью измерений.
Поэтому первую стадию комплексного эксперимента воспринимаем как этап начальной оценки коэффициентов Ь, в. Принципиальным моментом является единственность решения обратной задачи, поскольку данные, полученные на тонких мембранах, обобщаются (пересчитываются) на реальные конструкции.
Окончательная «доводка» результата осуществляется локальным варьированием предварительных значений Ь, в уже по модели быстрой водородопроницаемости (система ОДУ). Обычно довольствуются приемлемым визуальным совпадением экспериментальных и модельных графиков давлений. Теоретический анализ сопоставимости уровня погрешностей измерений (а известны лишь примерные оценки) и точности решения обратной задачи параметрической идентификации (математически трудная задача) практически нереален и не является первостепенно актуальным.
Работа выполнена при поддержке РФФИ (грант № 15-01-00744).
Литература
1. Алимов В. Н., Буснюк А. О., Ноткин М. Е., Лившиц А. И. Перенос водорода металлами5-й группы: достижение максимальной плотности потока сквозь ванадиевую мембрану // Письма в ЖТФ. 2014. Т. 40, вып. 5. C. 88-94.
2. Взаимодействие водорода с металлами / Ред. А. П. Захаров. М.: Наука, 1987. 296 с.
3. Водород в металлах / Ред. Г. Алефельд, И.Фелькль. М.: Мир, 1981. Т. 1,506 с.; т. 2,430 с.
4. Заика Ю. В., Родченкова Н. И. Диффузионный пик ТДС-спектра дегидрирования: краевая задача с подвижными границами // Математическое моделирование. 2008. Т. 19, №11. С. 67-79.
5. Заика Ю. В., Родченкова Н. И., Сидоров Н. И. Моделирование водородопроницае-мости сплавов для мембранного газоразделения // Компьютерные исследования и моделирование. 2016. Т. 8, №1. С. 121-135.
6. Изотопы водорода. Фундаментальные и прикладные исследования / Ред. А. А. Юхим-чук. Саров: РФЯЦ-ВНИИЭФ, 2009. 697 с.
7. Ленг С. Эллиптические функции. М.: Наука, 1984. 312 с.
8. Писарев А. А., Цветков И. В., Марен-ков Е. Д., Ярко С. С. Проницаемость водорода через металлы. М.: МИФИ, 2008. 144 с.
9. Уиттекер Э. Т., Ватсон Дж. Н. Курс современного анализа, часть 2. М.: Гл. ред. физ.-мат. лит., 1963. 516 с.
10. Черданцев Ю. П., Чернов И. П., Тюрин Ю. И. Методы исследования систем металл-водород. Томск: ТПУ, 2008. 286 с.
11. Dolan M.D. Non-Pd BCC alloy membranes for industrial hydrogen separation // Journal of Membrane Science. 2010. Vol. 362, iss. 1-2. P. 1228. doi: 10.1016/j.memsci.2010.06.068
56
0
1
2
7
8
9
12. Evard E.A., Gabis I.E., Yartys V.A. Kinetics of hydrogen evolution from MgH2: experimental studies, mechanism and modelling // International Journal of Hydrogen Energy. 2010. Vol.35, iss. 17. P. 9060-9069. doi: 10.1016/j.ijhydene.2010.05.092
13. Gabis I. E. The method of concentration pulses for studying hydrogen transport in solids// Technical Physics. 1999. Vol.44, iss. 1. P. 90-94. doi: 10.1134/1.1259257
14. Handbook of hydrogen storage: new materials for future energy storage / Ed. M. Hirscher. Wiley-VCH, 2010. 353 p.
15. Indeitsev D. A., Semenov B. N. About a model of structure-phase transfomations under hydrogen influence // Acta Mechanica. 2008. Vol. 195, iss. 1. P. 295-304. doi: 10.1007/s00707-007-0568-z
16. Kojakhmetov S., Sidorov N., Piven V., Sipatov I., Gabis I., Arinov B. Alloys based on Group 5 metals for hydrogen purification membranes // Journal of Alloys and Compounds. 2015. Vol.645, supp. 1. P. S36-S40. doi: 10.1016/j.jallcom.2015.01.242
17. Lototskyy M. V., Yartys V. A., Pollet B. G, Bowman R. C. Jr. Metal hydride hydrogen compressors: a review // International Journal of Hydrogen Energy. 2014. Vol.39, iss. 11. P. 58185851. doi: 10.1016/j.ijhydene.2014.01.158
18. Song G, Dolan M.D., Kellam M.E., Liang D., Zambelli S. V-Ni-Ti multi-phase alloy membranes for hydrogen purification // Journal of Alloys and Compounds. 2015. Vol. 509, iss. 38. P. 9322-9328. doi: 10.1016/j.jallcom.2011.07.020
19. Terrani K.A., Balooch M., Wongsawaeng D., Jaiyen S., Olander D. R. The kinetics of hydrogen desorption from and adsorption on zirconium hydride // Journal of Nuclear Materials. 2010. Vol.397, iss. 1-3. P. 61-68. doi: 10.1016/j.jnucmat.2009.12.008
20. The hydrogen economy / Eds M. Ball, M. Wietschel. Cambridge Univ. Press, 2009. 646 p.
21. Varin R.A., Czujko T., Wronski Z.S. Nanomaterials for solid state hydrogen storage. Springer, New York, 2009. 338 p. doi: 10.1007/9780-387-77712-2.
22. Voyt A., Sidorov N., Sipatov I., Dobrotvorskii M., Piven V., Gabis I. Hydrogen
solubility in V85Ni15 alloy // International Journal of Hydrogen Energy. 2017. Vol. 42, iss. 5. P. 30583063. doi: 10.1016/j.ijhydene.2016.10.033
23. Zaika Yu. V., Bormatova E. P. Parametric identification of a hydrogen permeability model by delay times and conjugate equations // International Journal of Hydrogen Energy. 2011. Vol.36, iss. 1. P. 1295-1305. doi: 10.1016/j.ijhydene.2010.07.099
24. Zaika Yu. V., Kostikova E. K. Computer simulation of hydrogen thermal desorption by ODE-approximation // International Journal of Hydrogen Energy. 2017. Vol. 42, iss. 1. P. 405-415. doi: 10.1016/j.ijhydene.2016.10.104
25. Zaika Yu. V., Kostikova E. K. Computer simulation of hydrogen thermodesorption // Advances in Materials Science and Applications. World Acad. Publ. 2014. Vol. 3, iss. 3. P. 120-129. doi: 10.5963/AMSA0303003
26. Zaika Yu. V., Kostikova E. K. Determination of effective recombination coefficient by thermodesorption method // International Journal of Hydrogen Energy. 2014. Vol. 39, iss. 28. P. 15819-15826. doi: 10.1016/j.ijhydene.2014.07.117
27. 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. 2009. Vol.33, iss. 10. P. 3776-3791. doi: 10.1016/j.apm.2008.12.018
28. Zaika Yu. V., Rodchenkova N. I. Hydrogen-solid boundary-value problems with dynamical conditions on surface // Mathematical Modelling. Nova Sci. Publishers. 2013. P. 269-302.
29. Zaika Yu. V., Rodchenkova N. I. Numerical modelling of hydrogen desorption from cylindrical surface // International Journal of Hydrogen Energy. 2011. Vol.36, iss. 1. P. 1239-1247. doi: 10.1016/j.ijhydene.2010.06.121
30. Zhang Y., Maeda R., Komaki M., Nishimura C. Hydrogen permeation and diffusion of metallic composite membranes // Journal of Membrane Science. 2006. Vol. 269, iss. 1-2. P. 6065. doi: 10.1016/j.memsci.2005.06.018
Поступила в редакцию 05.05.2017
References
1. Alimov V.N., Busnyuk A.O., Notkin M.E., Livshits A. I. Hydrogen transport by group 5 metals: Achieving the maximal flux density through a vanadium membrane. Technical Physics Letters. 2014. Vol. 40, iss. 3. P. 228-230. doi: 10.1134/S1063785014030031
2. Vzaimodeistvie vodoroda s metallami [Interactions of hydrogen with metals]. Ed. A. P. Zakharov. Moscow: Nauka, 1987. 296 p.
3. Vodorod v metallakh [Hydrogen in metals]. Eds G. Alefel'd, J.Fel'kl'. Moscow: Mir, 1981. Vol. 1,2. 506,403 p.
4. Zaika Yu. V., Rodchenkova N. I. Diffuzionnyi pik TDS-spektra degidrirovaniya: kraevaya zadacha s podvizhnymi granitsami [Diffusion
©
peak of TDS-spectrum of dehydrogenation: a boundary-value problem with moving bounds]. Matematicheskoe modelirovanie [Mathematical Models and Computer Simulations]. 2008. Vol. 19, iss. 11. P. 67-79.
5. Zaika Yu. V., Rodchenkova N. I., Sidorov N. I. Modelirovanie vodorodopronitsaemosti splavov dlya membrannogo gazorazdeleniya [Modeling of H2-permeability of alloys for gas separation membranes]. Komp'yuternye issledovaniya i modelirovanie [Computer Research and Modeling]. 2016. Vol.8, iss.1. P. 121-135.
6. Izotopy vodoroda. Fundamental'nye i pri-kladnye issledovaniya [Hydrogen isotopes. Fundamental and applied studies]. Ed. A. A. Yukhimchuk. Sarov: RFYaTs-VNIIEF, 2009.
7. Lang S. Elliptic functions. Addison-Wesley Publishing, 1973. 326 p.
8. Pisarev A. A., Tsvetkov I. V., Marenkov E. D., Yarko S. S. Pronitsaemost' vodoroda cherez metally [Hydrogen permeability through metals]. Moscow: MIFI, 2008. 144 p.
9. Whittaker E. T, Watson G. N. A Course of Modern Analysis. Cambridge Univ. Press, 1996. 612 p.
10. Cherdantsev Yu. P., Chernov I. P., Tyurin Yu. I. Metody issledovaniya sistem metall-vodorod [Methods of studying metal-hydrogen systems]. Tomsk: TPU, 2008. 286 p.
11. Dolan M.D. Non-Pd BCC alloy membranes for industrial hydrogen separation. Journal of Membrane Science. 2010. Vol. 362, iss. 1-2. P. 1228. doi: 10.1016/j.memsci.2010.06.068
12. Evard E. A., Gabis I. E., Yartys V. A. Kinetics of hydrogen evolution from MgH2: experimental studies, mechanism and modelling. International Journal of Hydrogen Energy. 2010. Vol. 35, iss. 17. P. 9060-9069. doi: 10.1016/j.ijhydene.2010.05.092
13. Gabis I. E. The method of concentration pulses for studying hydrogen transport in solids. Technical Physics. 1999. Vol.44, iss.1. P. 90-94. doi: 10.1134/1.1259257
14. Handbook of hydrogen storage: new materials for future energy storage. Ed. M. Hirscher. Wiley-VCH, 2010. 353 p.
15. Indeitsev D. A., Semenov B. N. About a model of structure-phase transfomations under hydrogen influence. Acta Mechanica. 2008. Vol. 195, iss. 1. P. 295-304. doi: 10.1007/s00707-007-0568-z
16. Kojakhmetov S., Sidorov N., Piven V., Sipatov I., Gabis I., Arinov B. Alloys based on Group 5 metals for hydrogen purification membranes. Journal of Alloys and Compounds. 2015. Vol.645, supp. 1. P. S36-S40. doi: 10.1016/j.jallcom.2015.01.242
17. Lototskyy M. V., Yartys V. A., Pollet B. G, Bowman R. C. Jr. Metal hydride hydrogen
compressors: a review. International Journal of Hydrogen Energy. 2014. Vol.39, iss. 11. P. 58185851. doi: 10.1016/j.ijhydene.2014.01.158
18. Song G, Dolan M.D., Kellam M.E., Liang D., Zambelli S. V-Ni-Ti multi-phase alloy membranes for hydrogen purification. Journal of Alloys and Compounds. 2015. Vol. 509, iss. 38. P. 9322-9328. doi: 10.1016/j.jallcom.2011.07.020
19. Terrani K.A., Balooch M., Wongsawaeng D., Jaiyen S., Olander D. R. The kinetics of hydrogen desorption from and adsorption on zirconium hydride. Journal of Nuclear Materials. 2010. Vol.397, iss. 1-3. P. 61-68. doi: 10.1016/j.jnucmat.2009.12.008
20. The hydrogen economy. Eds M. Ball, M. Wietschel. Cambridge Univ. Press, 2009. 646 p.
21. Varin R.A., Czujko T., Wronski Z.S. Nanomaterials for solid state hydrogen storage. Springer, New York, 2009. 338 p. doi: 10.1007/9780-387-77712-2
22. Voyt A., Sidorov N., Sipatov I., Dobrotvorskii M., Piven V., Gabis I. Hydrogen solubility in Vg5Ni15 alloy. International Journal of Hydrogen Energy. 2017. Vol.42, iss. 5. P. 30583063. doi: 10.1016/j.ijhydene.2016.10.033
23. Zaika Yu. V., Bormatova E. P. Parametric identification of a hydrogen permeability model by delay times and conjugate equations. International Journal of Hydrogen Energy. 2011. Vol.36, iss.1. P. 1295-1305. doi: 10.1016/j.ijhydene.2010.07.099
24. Zaika Yu. V., Kostikova E. K. Computer simulation of hydrogen thermal desorption by ODE-approximation. International Journal of Hydrogen Energy. 2017. Vol.42, iss. 1. P. 405-415. doi: 10.1016/j.ijhydene.2016.10.104
25. Zaika Yu. V., Kostikova E. K. Computer simulation of hydrogen thermodesorption. Advances in Materials Science and Applications. World Acad. Publ. 2014. Vol. 3, iss. 3. P. 120-129. doi: 10.5963/AMSA0303003
26. Zaika Yu. V., Kostikova E. K. Determination of effective recombination coefficient by thermodesorption method. International Journal of Hydrogen Energy. 2014. Vol. 39, iss. 28. P. 15819-15826. doi: 10.1016/j.ijhydene.2014.07.117
27. 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. 2009. Vol.33, iss. 10. P. 3776-3791. doi: 10.1016/j.apm.2008.12.018.
28. Zaika Yu. V., Rodchenkova N. I. Hydrogen-solid boundary-value problems with dynamical conditions on surface. Mathematical Modelling. Nova Sci. Publishers. 2013. P. 269-302.
29. Zaika Yu. V., Rodchenkova N. I. Numerical modelling of hydrogen desorption from cylindrical surface. International Journal of Hydrogen Energy. 2011. Vol.36, iss.1. P. 1239-1247. doi: 10.1016/j.ijhydene.2010.06.121
30. Zhang Y., Maeda R., Komaki M., Nishimura C. Hydrogen permeation and diffusion of metallic composite membranes. Journal of Membrane Science. 2006. Vol. 269, iss. 1-2. P. 6065. doi: 10.1016/j.memsci.2005.06.018
Received May 5, 2017
СВЕДЕНИЯ ОБ АВТОРАХ:
Заика Юрий Васильевич
рук. лаб. моделирования природно-технических систем, д. ф.-м. н. Институт прикладных математических исследований Карельского научного центра РАН ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: zaika@krc.karelia.ru тел.: (8142) 780059
Родченкова Наталья Ивановна
старший научный сотрудник, к. ф.-м. н. Институт прикладных математических исследований Карельского научного центра РАН ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: nirodchenkova@krc.karelia.ru тел.: (8142) 766312
CONTRIBUTORS:
Zaika, Yury
Institute of Applied Mathematical Research, Karelian Research Centre, Russian Academy of Sciences 11 Pushkinskaya St., 185910 Petrozavodsk, Karelia, Russia
e-mail: zaika@krc.karelia.ru tel.: (8142) 780059
Rodchenkova, Natalia
Institute of Applied Mathematical Research, Karelian Research Centre, Russian Academy of Sciences 11 Pushkinskaya St., 185910 Petrozavodsk, Karelia, Russia
e-mail: nirodchenkova@krc.karelia.ru tel.: (8142) 766312
