CC BY
61
Труды Карельского научного центра РАН № 10. 2015. С. 54-68 DOI: 10.17076/mat137
УДК 519.6:539.2
КРАЕВАЯ ЗАДАЧА ВОДОРОДОПРОНИЦАЕМОСТИ МЕМБРАН ГАЗОРАЗДЕЛЕНИЯ
Ю. В. Заика, Н. И. Родченкова
Институт прикладных математических исследовании Карельского научного центра РАН
Производство высокочистого водорода необходимо для экологически чистой энергетики и различных химико-технологических процессов. Значительная часть водорода будет производиться за счет конверсии метана, а также его выделения из других углеводородных газов, не вовлеченных в процесс производства энергии. Методом измерения удельной водородопроницаемости исследуются различные сплавы, перспективные для использования в газоразделительных установках. Требуется оценить параметры диффузии и сорбции с тем, чтобы иметь возможность численно моделировать различные сценарии и условия эксплуатации материала (включая экстремальные), выделять лимитирующие факторы. В статье представлены нелинейная модель водородопроницаемо-сти и ее модификации в соответствии со спецификой эксперимента, разностная схема решения краевой задачи и результаты численного моделирования.
Ключевые слова: водородопроницаемость, нелинейные краевые задачи, разностные схемы, численное моделирование.
Yu. V. Zaika, N. I. Rodchenkova. BOUNDARY-VALUE PROBLEM OF HYDROGEN PERMEABILITY OF GAS SEPARATION MEMBRANES
High-purity hydrogen is required for clean energy and a variety of chemical technology processes. A considerable part of hydrogen is to be obtained by methane conversion and its separation from other hydrocarbon gases not involved in energy production. Different alloys, which may be well-suited for use in gas-separation 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 experimental conditions of the material usage (including extreme ones), and identify the limiting factors. This paper presents a nonlinear model of hydrogen permeability and its modifications in accordance with the specifics of the experiment, the difference scheme for the solution of the boundary-value problem, and the results of numerical modelling. This work is supported by the Russian Foundation for Basic Research (Project No. 15-01-00744).
Key words: hydrogen permeability, nonlinear boundary-value problems, difference schemes, numerical simulation.
Введение
Интерес к взаимодействию водорода с различными материалами носит многоплановый характер [1-10]. Достаточно упомянуть задачи энергетики, защиты металлов от водородной коррозии, проектирования химических реакторов, ракетостроения. Гидриды позволяют удерживать большое количество водорода. С этим связаны перспективы водородных аккумуляторов и двигателей с высоким уровнем безопасности: без высоких давлений и низких температур. На обратимом легировании металлов водородом основаны пластифицирование и термоводородная обработка титановых сплавов. Некоторые частные задачи материаловедения исследованы в [11-15]. Энтузиасты говорят не только о водородной энергетике, но и о водородной экономике [7].
Экспериментальный опыт показывает, что лимитирующими являются не только диффузионные процессы, но и физико-химические явления на поверхности [1, 2]. Параметры переноса зависят и от технологических особенностей получения партии материала, поэтому вряд ли следует ориентироваться на получение «табличных данных», нужны эффективные алгоритмы обработки экспериментальных кривых. В статье остановимся на методе проницаемости, учитывая лишь основные факторы и информационные возможности рассматриваемого эксперимента. Основой для проведенных исследований послужили работа [16] и данные по водородопроницаемости некоторых перспективных сплавов [20].
Модель водородопроницаемости
Вначале кратко опишем эксперимент. Образец конструкционного материала, нагретого до фиксированной температуры, является перегородкой вакуумной камеры. Предварительно проведена дегазация. В начальный момент времени на входной стороне создается давление молекулярного водорода (впрыскивается порция газа). Измеряется падающее давление на входе и растущее давление водорода в выходной емкости. Информационные возможности эксперимента ограничены, поэтому в модели водородопроницаемости будем учитывать только основные факторы для прикладной задачи мембранной фильтрации.
Рассмотрим перенос водорода сквозь образец материала (пластину толщины £ и площади 5). Температура Т постоянна в течение одного эксперимента. Концентрация растворенного водорода (в атомарном состоянии) относительно мала, и диффузионный поток можно считать пропорциональным градиенту кон-
центрации. Часть атомов Н взаимодействует с ловушками (микродефекты различной природы, включая микрополости), которые могут удерживать водород. Ориентируясь на прикладной смысл задачи и возможности метода проницаемости, ограничимся представлением об «ограниченном стоке» без дополнительной детализации. В качестве модели диффузии с ограниченным захватом в объеме примем нелинейную систему уравнений
dc d 2 c
Ж = DT>8X2 - fT.z.c).
I = f ^ aT)
1 -
z(t, x)
(1)
c - flout(T)z, (2)
где c(t, x) - концентрация диффундирующего водорода (атомарного); z(t,x) - концентрация захваченного диффузанта; D - коэффициент диффузии; a = ain и aout - коэффициенты поглощения и высвобождения атомов H ловушками. Знак тождества часто используем в смысле равенства по определению. Величину zmax считаем малой, захват носит характер поправки и не требует более детального моделирования. Для конструкционных материалов (рассматривается металлический сплав) в рабочем диапазоне температур T е [500, 800] K, как правило, aout ^ ain ив процессе насыщения (проницаемости) нет необходимости усложнять модель (далее полагаем aout = 0). Величины D, а зависят от температуры T образца по закону Аррениуса с предэкспонен-циальными множителями Do, ao и энергиями активации Ed, Ea (R - универсальная газовая постоянная): D = D0exp{—ED/[RT(t)]}, a = a0 exp{—Ea/[RT(t)]}. Начальные данные: в силу предварительной дегазации
c(0,X = 0, z(0,x) = 0, x е [0,4 (3)
Из баланса потоков получаем следующие нелинейные граничные условия:
dQir dt
= [^(T)s(T)po(t) —
dc
— b(T )c2(t)]S = —SD-
ж=0
dQout dt
= [^(T)s(TMt) —
dc
— b(T )c2(t)]S = SD-
x=t
(4)
(5)
Здесь фп(*), фои^) - количество атомов водорода во входной емкости объема Ип и выходной емкости объема РО^, с0(£) = с(£, 0), = с(£,£). Газообразный водород в рассматриваемом «рабочем» диапазоне температур находится в молекулярной форме, но для
55
z
max
единообразия, поскольку сквозь металлическую мембрану диффундирует атомарный водород, подсчет ведем в атомах. Согласно кинетической теории газов, плотность Jp падающего на поверхность потока частиц связана с давлением р по формуле Герца-Кнудсена: Jp = рД/2nmkT (k - постоянная Больцма-на, m - масса молекулы водорода). В контексте рассматриваемой методики эксперимента удобно в качестве единиц измерений выбрать [£] = cm, [р] = Torr. Тогда численно получаем зависимость Jp = ßp, ß(T) и 2.474 ■ 1022/vT ([ß] = 1h2/(Torrcm2 s), [T] = K). На поверхности происходят процессы физической адсорбции, хемосорбции, диссоциации молекул на атомы, растворения. Лишь малая часть «налетающих» атомов H окажется в абсорбированном состоянии в объеме. Это отражается множителем s. Итак, ßsp - результирующий поток атомов в объем сквозь поверхность без разделения на более элементарные стадии. По контексту слово «плотность» обычно опускаем. Можно вместо s написать 2s и интерпретировать безразмерный вероятностный множитель s как долю абсорбируемых атомов H.
Далее, J0,e = bc^ ([J] = 1н/(cm2 s)) - это плотности десорбции из образца (квадратич-ность является следствием соединения двух атомов водорода в молекулу), b - коэффициент десорбции. Для s и b также предполагаем аррениусовскую зависимость от температуры. По крайней мере формально: в экспоненте «энергия активации» Es может оказаться и отрицательной величиной как линейная комбинация энергий активаций и теплот поверхностных процессов на пути «из газа в раствор».
Отметим следующее. Если с обеих сторон мембраны поддерживать постоянное давление насыщения р молекулярного водорода при постоянной температуре T, то через некоторое время установится равновесная концентрация с растворенного атомарного диффузионно подвижного водорода. Из модели (4), (5), приравнивая производные к нулю, получаем с = rVP, Г = л/ßs/b. Таким образом, модель соответствует диапазону адекватности закона Сивертса (с гс д/р).
Уточним экспериментальные условия. Объемы Vin,out - несколько литров, толщина мембраны £ меньше mm, площадь S - около cm2, давление напуска p0(0) - несколько десятков Torr. Диапазон [pmin,pmax] невелик, ограничимся zmax = ис, и ^ 0.1. Это не приведет к нарушению закона Сивертса (с + zmax гс д/р), причем с + zmax и с = Г^/р в пределах экспериментальной точности.
Остается определить величины Qin, Qout. В масштабе времени установления диффузии газ находится в термодинамическом квазиравновесии с поверхностью, поэтому воспользуемся формулой N = рУ/(kT). Здесь N - количество частиц газа, занимаемого объем У при температуре T и давлении р (в системе СИ [р] = Pa, [У] = m3, [k] = J/K). С учетом соотношений Torr = 133.322 Pa, Pa = J/m3 (формально), получаем для соответствующих давлений и объемов в граничных условиях (4), (5) Q = 2N = арУ/T, а и 1.931 ■ 1019. Здесь р, У, T означают численные значения в выбранных единицах (Torr, cm3, K).
Функция стока
Остановимся более подробно на функции стока f(T,z,c) = a[1 — zz—lx]c ([a] = 1/s). Скорость поглощения атомов H увеличивается с ростом концентрации диффузанта и уменьшается по мере заполнения ловушек. Уравнение диффузии обладает неограниченной скоростью передачи возмущения, но переход к модификации уравнения «с конечной скоростью» не представляется целесообразным. Более разумно учесть захват части атомов неоднородностями структуры материала. Косвенным подтверждением являются трудности полной дегазации - процесс требует высоких температур, глубокого вакуумирования и длительного времени. Диффундирующий растворенный водород относительно быстро десорбируется, а «оторвать» захваченные атомы H значительно труднее. Как уже упоминалось, учесть возможность полной дегазации можно введением дополнительного параметра: f = ain [1 — zz—¿У c — aoutz. Захват ограниченной емкости можно интерпретировать и следующим образом. Часть атомов водорода идет на «расширение решетки» металла, причем эта связь приоритетна по сравнению с диффузионным перемещением вглубь. Тогда при достаточно большом коэффициенте a сначала будет преимущественно заполняться эта «ловушка», обеспечивающая и определенное запаздывание переноса водорода.
Величину zmax можно считать независимым параметром модели. Связь zmax = ис принята по следующим соображениям. Значение zmax зависит от материала и внешних условий, что «коррелирует» с равновесной растворимостью диффузионно подвижного водорода. В рамках модели с = ^s^^T). Полагая в указанном диапазоне температур и давлений эту связь монотонной, принимаем пропорцию zmax = и с (и ^ 1) с учетом того, что материал однороден и захват невелик. Общая растворимость определяется кон-
®
центрацией с + zmax a p. Это согласуется с известным мнением [17, с. 512]: «Правильнее, по всей вероятности, считать, что часть водорода прочно удерживается дефектами решетки и поэтому не должна учитываться при определении концентрации, градиент которой входит в уравнение закона Фика, т. е. часть водорода является полуинертной примесью. Закон Фи-ка применим только к оставшейся части».
Итак, эффект захвата (включая «расширение решетки») учтен без дополнительных параметров в соответствии с ограниченной информативностью метода проницаемости (на входе давление монотонно падает, на выходе монотонно растет до установления). Другое дело, например, пористый вольфрам [3], когда микропоры достаточно велики для рекомбинации атомов водорода в молекулы. Это уже требует более детального моделирования захвата.
Возможные модификации
Представленную модель (1)-(5) будем считать базовой. При необходимости можно учесть емкость поверхности как дополнительной «удерживающей силы». Кроме того, на поверхности могут образовываться заметные оксидная пленка и гидридная фаза (по существу слой другого материала), что существенно сказывается на водородопроницаемости. За счет химических связей в гидриде значительно больше водорода, чем в том же объеме жидкого водорода. В принципе можно считать мембрану трехслойной и учесть динамическое накопление в приповерхностных слоях, в том числе и за счет гидридообразования. Появятся дополнительные коэффициент диффузии и коэффициенты в условиях сопряжения на стыках слоев (равенство диффузионных потоков, но скачки концентраций) . . . Далее, при большом перепаде давлений мембрана испытывает изгиб. Простейший вариант (без учета напряжений и деформаций) - считать значения s, b различными при x = 0,£. При такой детализации образуется «снежный ком» параметров с неизвестными априори значениями.
Для примера остановимся на «интегральном» учете влияния включений гидридной фазы. Выделим слой малой толщины Л ^ £, который будет играть роль «физической» поверхности. Вместо Zmax = const полагаем Zmax(x): Zmax = аС (x G [Л, £ — Л], a < 0.1) и
Zmax = (С — 2(Z — a^-1x+
+ (Z — a^-2x2, Z> 1, x е [0, Л],
Zmax = Сс + 2(С — a^-1(x — £) + (6)
+ (Z — a^-2(x — £)2, x е [£ — Л,4
Прокомментируем эти формулы. Внутри мембраны (ж € [Л, £ — А]) учитываем только малый захват дефектами структуры материала. Образование гидрида связано с изменением объема (как правило, с увеличением). Если это происходит внутри в заметном масштабе, то материал растрескивается и непригоден для рассматриваемых целей. На поверхности возникнуть гидридному «вздутию» значительно легче. Растет зародыш гидрида за счет растворенного атомарного водорода, из газовой фазы поглощение значительно медленнее. Обозначим (с (£ > 1) максимальную концентрацию водорода в приповерхностном слое с учетом частичной заполненности гидридом. Формулы (6) означают, что максимальная емкость ловушек быстро возрастает в приповерхностных слоях (используем квадратичную зависимость). Параболы гладко «стартуют» (с нулевой производной) с уровня ас при ж = Л, ж = £ — Л и стремятся к ("с по мере приближения к геометрическим входной и выходной поверхностям. Различные «емкости» считаем пропорциональными равновесной растворимости диффундирующего водорода, так что при полной дегазации насыщенного водородом образца останется пропорциональность общей концентрации корню из давления насыщения. Отметим также, что в случае ¿тах(ж) увеличивается интенсивность захвата /(Т, г, с) именно у поверхности (при равных г, с): там нет «симметрии сил» и дефектов больше.
Модифицируем изменение «свободных» поверхностей, через которые идет процесс растворения и десорбции. Обозначим их 50/ (геометрические объекты или площади - различаем по контексту). При частичной занятости начальной 5 «непроводящими» гидридными включениями имеем 50/ < 5. Примем модель
50 = 5, £ < г [г(г, 0) < ас, 0) = ас], 50 = + (1 — [Сс — г(г, 0)][(с — ас]-1,
£ € (0,1), г ^ £ По постановке эксперимента и модели динамика концентраций носит монотонный характер. Водород начинает поступать сквозь поверхность 50 = 5 с наибольшей локальной интенсивностью захвата, так что быстрее всего уровень насыщения дефектов ас достигнет концентрация 0) (если вообще это произойдет, иначе 50 = 5 и ничего не меняем). Если наступит ситуация 0) = ас (с дальнейшим ростом 0)), то это свидетельствует о начале образования и роста зародышей гидридной фазы. С этого момента 50 начинает плавно уменьшаться. В случае до-
стижения максимума ("с уменьшение Бо останавливается на уровне < Б.
На выходе в приповерхностном слое концентрация диффундирующего водорода снижается, но быстрее растет гтах(ж). Примем аналогичную модель для Бе (с заменой г(£, 0) на г(£, •£)). Имеется определенная асимметрия: момент начала уменьшения Бе наступит позже. Добавим также, что неявно предполагаем гидридные включения достаточно дискретными, чтобы в целом перенос осуществлялся в перпендикулярном к поверхности направлении. Пренебрегаем и скоростью распада гидрида, что возможно при значительном снижении давления на входе.
Отметим, что задействованные дополнительные параметры Л (толщина «физической» поверхности), ("с (максимальная концентрация в ловушках с учетом гидридной фазы), £ (минимальная доля свободной поверхности) можно грубо оценить из физических соображений или даже визуально по итогам эксперимента. Значения зависят не только от материала, но и от условий эксперимента (в частности, давления напуска и температуры). Предложенная модификация (зависимости гтах(ж), 5о,е(£)) приводит лишь к некоторым техническим более громоздким выкладкам в излагаемом далее вычислительном алгоритме.
Безразмерная форма краевой задачи
Представим компактно базовую модель:
I = ^)£ - )
1 -
г(£, ж)
с(£, ж),
дг
= а
1
г(£, ж)
с, с(0, ж) = 0, г(0,ж) = 0,
МТЖТ)ро,е(£) - Ь(т)с0,е(^) = ТДТ) —
ж=о,е
т,ои1
= -[МТ)з(Т)ро,е(£) - Ь(т)с0,е(*)]Б,
-1
ди
дТ
д 2и ду2
— а
1 - ^ «т
и,
дТ
= аа
1 - ^ «т
и,
и,« е [0, 1], у е (0, 1), Т > 0, «т = «т а(Т) = а^-1, ио,1(т) = и(т, у) |у=о,1, и(0,у) = 0, «(0,у) = 0, у е [0,1],
род(т) = Ро,е(£)Ро-1, Р о = ^Ро, Ж (Т) = ЬсШ-1 = Р^-Ос]-1,
ди
Ж[ро - ио] = - ду
у=о
ди
- «Я = ^
= = Vmс,
У=1
^ = -^т^ [ро - и2], Ро(0) = 1,
^ = -^т^ [Р1 - и2], Й(0)=0,
= РоТ 1,
= РоТ 1.
= аРо,е(^)Ип,ои1Т
Выберем соответствующие нормировки. Температура Т фиксирована. По максимальному давлению ро = ро(0) определим соответствующие равновесную концентрацию диффундирующего водорода с = Гу'ро и количества атомов ^п = а^оИп/Т, = «РоКи /Т. Перейдем к безразмерным координате у = ж/^, концентрациям и = с/с, у = г/с (гтах = и с ^ с) и времени т: £ = ([^] = ст2Д). Получаем следующую краевую задачу:
«Емкость» равна количеству атомов водорода в образце в режиме равновесного насыщения при давлении с = Ро(0) и температуре Т. Отметим, что транспортный параметр Ж [3] играет определяющую роль при анализе варианта метода проницаемости, когда на выходе производится постоянное вакуумиро-вание (метод прорыва).
Замечание. Формально нулевые начальные и граничное условие Ш[ро - и2] = -дуи|у=о не согласованы при £ ^ +0 (ро(0) = 1, ио(0) = 0, дуи(0, 0) = 0). На самом деле «мгновенный» напуск водорода на входе длится некоторое время, пусть и пренебрежимо малое. В вычислительном алгоритме решения краевой задачи это фактически учтено. В принципе дискретную аппроксимацию можно считать моделью «в атомах», а краевую задачу - ее компактным представлением, «континуальным замыканием». Иначе следует вести речь в терминах теории обобщенных решений.
Модель квазистационарной проницаемости
Для приемлемого времени установления равновесия или стационара в экспериментальной практике используются тонкие мембраны (доли миллиметра). Материал для газоразделения подбирается с высокой водородопрони-цаемостью. По этим причинам быстро устанавливается квазистационарный режим, когда распределение атомов водорода практически линейное с относительно медленным дрейфом. В такой ситуации нерационально решать
г
тах
г
тах
на каждой итерации идентификации модели нелинейные краевые задачи при текущих приближениях параметров. Изложим квазистационарную модель водородопроницаемости.
Экспериментально фиксируется время ¿0 ^ , когда на выходе начинается заметный рост давления молекулярного водорода, ¿¿¿(¿0) > 0. Подразумеваются масштабы давления превышения шума и времени установления проницаемости ¿*. За время ¿о заполняются ловушки малой емкости (¿тах ^ с) и устанавливается линейный квазистационар: с(£0,х) и с(£0,х) = А(^)х + В(¿0) (В и С0, А < 0). Выходной поток еще незаметен: ¿(¿0,£) = А(^)£ + В(¿0) = 0. По падению входного давления вычисляем количество «пропавших» атомов ^¡п — аУтр(£0)/Т. Захвачено ££.гтах атомов Н, остаток равен
5/0 [А(*0)ж + В(^)]—х = 5[В(¿0)£ + А(^0)^2/2]. Это позволяет при фиксированном ¿тах определить значения А(£0) и В(^). При £ = ¿0 приближенно имеем аРО^Т-1°ре(10) =
dQout _ _SDdc dt dx
x=l
_ —SDA(to) > 0.
Это соотношение дает принципиальную возможность оценки коэффициента диффузии В. Измерения давления зашумлены, операция дифференцирования вычислительно некорректна, момент времени ¿0 и наклон графика Рг(Ь) достаточно условны, так что имеется в виду лишь начальная грубая оценка.
Дальнейшие выкладки проведем уже для безразмерной модели. Поскольку ловушки уже заполнены, то / = 0 и для уравнения диффузии квазистационаром (дти и 0) является линейное распределение
и(т, у) и А(т)у + В(т), т ^ Т0 = ВГ2^, а4(Г0) = А^с-1, В(т0) = В (¿0 )с-1.
Остается определить А(т), В(т) (т > т0). Воспользуемся граничными условиями с учетом
ду и = А:
W [р0 — В2] = — А, Ж [р>1 — [А + В]2] = А,
в0 0 = А в "—Г = —А4, в0,1 = ^1п,си1 ФтЛ
Из последних двух уравнений получаем 7<-р0/-т = —йр1 /—т, где 7 = = Ип/Кш^
Условие согласования наклонов 7^0^) = —¿¿¿(¿) (£ ^ ¿0) означает равенство потоков на входе и выходе мембраны и может считаться критерием выхода на квазистационарный режим проницаемости. Разумеется, по-
добные асимптотические равенства проверяются с большой погрешностью и носят скорее вспомогательный качественный характер. Интегрируя с учетом ¿51 (то) _ 0, имеем простую связь давлений р\(т) _ 7[ро(т0) — ¿50(т)], т ^ т0. Остается определить давление ¿50.
Заменяем в первом уравнении A _ во5о (¿50 = dp0/d-т) и выражаем
B_ {¿5о + eoW-15о}1/2.
В силу ро < 0 и eoW-1р5о _ —¿50 + u2 > —¿50 под корнем положительное число. Заменяя во втором уравнении A, B полученными выражениями, приходим к соотношению
y[5o(to) — ¿5о ] — [во5о + {5о + во W-15о }1/2]2
_ eoW-1р5о, ¿50 _ ¿50(т), т ^ то. (7)
Это обыкновенное дифференциальное уравнение первого порядка, неразрешенное относительно производной. Не будем искать выражение dpo/dT _ F(¿50) с целью последующего интегрирования, а рассмотрим (7) как алгебраическое уравнение относительно производной.
Удобно ввести нормированную переменную Z(т) _ —e0W-1dj50/dT £ (0,1), которая монотонно убывает и удовлетворяет соотношению вида Z _ F(Z):
y[5o(to) — ¿50 (т)] — [V5o(t ) — Z (т) — WZ (т )]2
_ —Z(т), т ^ то.
(8
Начинаем с т = тз. Начальное приближение Z(т0) = — р0(£0)в0£2/(Жр0В) можно найти грубо по начальному наклону графика (¿¿(¿0) = —7^0(^0)). Для поиска неподвижной точки Z = ^(Z) применяем метод последовательных приближений Zfc+l = ^ (Zfc). По Z(т0) вычисляем —¿50/—т = ¿50 = —Z(т0)Ж/в0 и ¿50(т0 + ¿т) и ¿50(т0) + ¿50£т с малым шагом ¿т. Далее смещаемся в точку т = т0 + ¿т и снова решаем уравнение (8) относительно Z(т0 + ¿т). Итерационный процесс описан. Он не связан с решением краевых задач и требует незначительных вычислительных ресурсов.
Отметим, что формально в уравнение (8) входит только параметр Ж. Напомним, что
W _ bc£D-1 _ Po£[DC]
1
P о
с = Г/С0, г = [^Г1]1/2, £ = В-1£2т.
Поэтому в исходных переменных расчет модельных давлений ¿0(£), потребует текущих приближений В, Ь, 8 (помимо входных данных 5, Ут,ои1, Т, ¿0 = ¿0(0), ¿0).
59
Представляет интерес сравнение квазистационарного приближения с квазиравновесным, когда плотность потока водородопрони-цаемости сквозь мембрану моделируется формулой Ричардсона
1
^ = -яг /Ш - л/рос£)]^
Она получается, если в выражении =
- градиент концентрации заменить разностным отношением (се - со) /^ с подстановкой равновесной зависимости с = Г^р. Учет общей концентрации с + гтах не изменит .
Насколько такие приближения правомерны, может показать (в рамках принятой модели) решение исходной краевой задачи, включающей и начальный (во многом определяющий) этап водородопроницаемости.
ВЫЧИСЛИТЕЛЬНЫЙ Алгоритм
Следуя технике разностных схем, введем сетку О = {ут = тЪу, т = 0,1,..., М} (Ъу = 1/М) по пространственной переменной и сетку по времени ш = {тп = пЪт, п = 0, 1,...}. Обозначим через {ит}, {Кт} приближенные значения концентраций в объеме (и(тп,ут)) и ловушках («(тп, ут)). Для уравнения диффузии в безразмерной форме рассмотрим неявную разностную схему
ЦП+1 - ип ип+1 - 2ип+1 + ип+1
ът
ъу
- а [1 - Кп+1 ]ип+1,
т/п = Кп
г т "та^ т
(9)
Значения в начальный момент известны: и° = = 0 (0 ^ т ^ М). Следуя методу прогонки, ищем приближенные значения концентрации в узлах сетки на (п + 1)-м слое в виде
иП+1 = ат+1 иП+1 + вт+1, т = 0,..., М -1. Прогоночные коэффициенты (т: 1,..., М-1):
а т+1 =
Ат +1 =
2 + в - ат + а^у [1 - Ктп+1]'
вт + в ит
2 + в - ат + а М [1 - КП+1].
(11)
Для нахождения начальных коэффициентов а1, в1 воспользуемся следующими соображениями. Подсчитаем предварительно ип+1 по явной разностной схеме (в равенстве (9) справа заменяем п + 1 на п). На (п + 1)-м слое по времени аппроксимируем производную дуи|у=о - [-Эи^1 + 4Цп+1 - и^1 ]/2Ъу и подставим в условие Ж[ро - и2] = -дуи|у=о. Выражение величины РП+1 для подстановки в левую часть берем из аппроксимации
От = -^тЖ [Ро - ио] ^
рп+1 _ рп
Оп= (-) -ОтЖ[РП+1 - 0ип+1)2]
^+1 = РП + А1(иоп+1)
^ ро = "
2
_ ЪтОтЖ
1 + А
1
А1 — ^
Здесь, на начальном этапе, чтобы иметь возможность использовать алгоритм прогонки, величину КП+1 предварительно вычисляем из неявной схемы по К, сохраняя значение и,
п :
т/
Кп+1 _ т/п
ът
= а«-ах[1 - КП+1 ]ц
т/п + ъ а«-1 ип
т/п+1 = тт + ът а^тахит
1 + ът а«
ип тахит
(10)
ип-1 - [в + 2 + аъ2(1 - ТС+1)] ип+1+ + ип+1 + вип = 0,
в — Ъу Ъ— .
В итоге получаем и^1 = /о^^1, и2п+1) (положительный корень квадратного уравнения):
В — , и^1 =0.5
1 2Ъу Ж' о
- ЭВ1+
Это соответствует и последовательности во времени захвата диффундирующего атома. Данные соотношения рассматриваются во внутренних узлах сетки (т = 1, . . . , М - 1).
Рассмотрим уравнения перехода с п-го на (п + 1)-й слой по т (п ^ 0, 0 <т< М):
+ л/9В2 - 4.В(ип+1 - 4и1п+1) + 4/Р(
Зная выражение ип+1 = а1ип+1 + в1 и значение ип+1, получаем а1 =0, в1 = ип+1.
Поскольку при £ = 0 в континуальной модели имеем формально скачок на входе, то проверим корректность вычисления ио1 на начальном этапе: и/ = и1 =0 ^ и1 е (0,1) УВ > 0.
Ближайшая цель - найти значение иМ+1, необходимое для реализации прогонки. Предварительно подсчитываем иМ^м-1 по явной разностной схемео. Далее опре)деляем зависимость иМ+1 = /м (и]М+11, и_п+_12) из граничных условий при у = 1, используя аппроксимацию
1
dyu|y=i « [UM+2 - 4UM+-1 + 3UM+1]/2h
pin. pout [Torr]
1 + A2 2hy W
A2 ^ hrgmW, UM+1 = 0.5x
Q.
out
B2 =
- 3B2 W9B2 - 4B*(UM+-2 - 4UM-\) +4pn
Следующий этап: с текущими приближениями ЦП+1, ЦМ+1 решаем обратным ходом прогонки трехдиагональную систему линейных уравнений и находим новые приближения концентраций ЦП^1, Ц^й м- 1 (и остальные значения ЦП+1 для т = 3,..., М — 3).
После этого корректируем значения V"+1, заменяя ЦП в (10) на текущие приближения ит+1, и значения коэффициентов (11) (выражения а1 = 0, в1 = Ц0П+1 остаются). Далее снова пользуемся формулами
и0га+1 =/о(иГ+\иП+1),
ипм+1 = /м №Л,Ц5-2)
и повторяем вычисления, возвращаясь к предыдущему абзацу, до установления граничных значений ЦПм1 (обычно 2-3 итерации).
Заключительный этап: переходим к следующему слою по т, вычислив род по формулам
-+1 = РП + U0ra+1)
pQ =
2
1 + A
1
irn+1
_pn + A^ UM+1) p1 =
2
1 + A'
2
Результаты моделирования
Для определенности ориентируемся на экспериментальные возможности [16] и данные по сплаву V85NI15 [20]. По результатам вычислительных экспериментов а = 0.1, причем параметр скорости захвата в широком диапазоне a ^ 0.1 практически не влияет на проникающий поток, поскольку ловушки быстро насыщаются. Ориентировочно i œ 0.02 cm, S œ 0.6cm2, Vin œ Vu œ 1500cm3.
Модельные кривые давлений молекулярного водорода представлены на рис. 1-3. Скорость убывания «просвета» между кривыми (pin(t) — pout(t)) указывает на рост проницаемости с увеличением температуры.
400 C
D = 1.27e-4 Pin
b = 7e-23
s = 9e-5
\bar p0 = 32.87
\mu = 9.5e+20
a = 0.01, sigma = 0.1
zmax = 0.1\bar c = 2e+19 Pout
Г = 3.5e+19
50 55 60 time, min
Рис. 1. Динамика давлений, 400 °C
Pin. Pout [Torr]
450 C
D = 1.36e-4
b = 8.55e-23 Pin
s = 1.93e-4
\bar p0 = 32.09
\mu = 9.2e+20
a = 0.01, sigma = 0.1
zmax = 0.1\bar c = 2.58e+19 Pout
Г = 4.56e+19
0 4 8 12 16 20 24 28 32 36 40
time [min]
Рис. 2. Динамика давлений, 450 °C
pin. pout [Ton"]
500 C
D = 1.43e-4
b = 8.76e-23
s = 2.01e-4 Pin
\bar p0 = 32.55
\mu = 8.9e+20
a = 0.01, sigma = 0.1
zmax = 0.1\bar c = 2.58e+19 Pout
Г = 4.52e+19
0 4 8 12 16 20 24 28 32 36 40
time [min]
Рис. 3. Динамика давлений, 500 °C
61
32
28
24
20
16
12
8
4
0
0
5
10
15
20
25
30
35
40
45
32
28
24
20
16
12
8
4
0
32
28
24
20
16
12
8
4
0
Что касается обратной задачи параметрической идентификации, то следует учитывать ограничение монотонности с ростом температуры, что соответствует росту водородопрони-цаемости (рис. 1-3) и ориентировано на арре-ниусовскую зависимость физико-химических параметров переноса.
Тут возникает некоторая опасность. Для определенности остановимся на коэффициенте объемной десорбции (эффективном коэффициенте рекомбинации) b = bvoi. Рассмотрим баланс потоков, учитывая поверхность как потенциальный барьер (см. [2, с. 177-206; Габис, Компаниец, Курдюмов]):
<Mi) = MT )s(T (t)[1 - Mi)]2-
- bsurf (T)9o2^(i) ± D(T)dxc|o,e ,
k-(T)[1 - co;e(i)cmix]9o;^(i)--k+(T) [1 - 0o,*(i)] co,*(t) = T D(T)Сл |x=0>£ ,
где q [1H/(cm2)] - поверхностная концентрация, 0(t) = q(i)/qmax - степень заполнения поверхности, множитель [1 -0]2 учитывает необходимость двух сорбционных центров для диссоциации молекулы водорода на атомы. Величина qmax ~ 1015 лимитируется монослоем атомов водорода на поверхности. Первое уравнение читается так: рассогласование потоков адсорбции, десорбции и диффузии (справа) идет на накопление атомов H на поверхности. Второе уравнение согласовывает потоки растворения и выхода из объема на поверхность. При относительно малых концентрациях и накоплении (0 ^ 1, cq ^ cmax, q « 0), медленной диффузии по сравнению с растворением на поверхности, имеем в пределе c = gq, g = k-/k+ и b = bvoi = bsurf/g2. Косвенно эти рассуждения фиксируют достаточно высокие температуры адекватности исходной базовой модели. Предполагая bsurf, k-, k+ аррениусов-скими, можем в принципе получить «отрицательную энергию активации» у коэффициента эффективной рекомбинации b. Ориентируясь, однако, на широкий спектр исследований конструкционных материалов [3, 5], мы постулировали в процессе параметрической идентификации монотонный рост b по температуре.
Еще одну опасность продемонстрируем на рис. 4. Согласно принятой модели, равновесная концентрация с при условиях насыщения р = const, T = const определяется (после приравнивания к нулю всех производных) как с = Г^р, где Г = л/^s/b - коэффициент растворимости. Но следует помнить, что учтен только диффузионно подвижный растворен-
ный водород. При непродолжительной дегазации выделится именно он. Если нацеливаться на полную дегазацию, то растворимость другая: с + zmax = rmaxyp. В принятой модели имеем Гтах = Г[1 + а] при сохранении закона Сивертса с + zmax a ур. По сравнению с рис. 1 на рис. 4 с гипотетическим значением zmax = 10 с фиксирована общая растворимость на порядок больше. Если фиксировать установившийся поток водородопроница-емости J = —Ddxc в эксперименте прорыва, когда p0(t) = р = const, а на выходе вакууми-рование, то в предположениях c0 = с0 = Гур,
ce = 0 имеем J = Фур, Ф = Dr = Dy^s/b. Если же коэффициент проницаемости Ф вычислять по имеющимся значениям коэффициентов диффузии и растворимости, то можем получить Drmax. По-видимому, это одна из причин разброса литературных данных по растворимости и проницаемости. Проблема в том, что на очень тонких мембранах в режиме проницаемости трудно обнаружить «пропажу» водорода в ловушках даже при большой их емкости, существенно влияющей на общую растворимость в материале (при пересчете на кубические cm и m в условиях S ^ 1). С целью повышения точности параметрической идентификации следует дополнительно провести эксперименты насыщения-дегазации или проницаемости более массивного образца (^ ~ mm), чтобы заметнее проявились объемные характеристики материала на фоне поверхностных. Визуальным критерием сопоставимости влияния поверхности и объема можно считать заметную несимметричность графиков давлений p0,e(t) (с поправкой на различие объемов Vin)0ut).
Pin, Pout [Torr]
400 C
D = 1.27e-4
b = 8.15e-23 P,n
s = 9.4e-5
\bar p0 = 32.87
\mu = 9.54e+20
sigma = 10, a = 0.01
Г = 3.32e+19
rmax = r(1+sigma)=3.65e+20 ^^^ pout
zmax = 10\bar c = ______________
0 5 10 15 20 25 30 35 40 45 50 55 60
time [min]
Рис. 4- Давления при zmax = 10с, 400 °C
32
28
24
20
16
12
8
4
0
Перейдем к анализу динамики объемных концентраций в рассматриваемом эксперименте «сообщающихся сосудов». Для определенности далее ограничимся T = 400 °C, поскольку качественная картина остается неизменной и при других температурах в рассматриваемом диапазоне. Начнем с приповерхностных (ж = 0,£) - рис. 5,6. В течение минуты происходит практически полное заполнение ловушек: z0 ~ zg & zmax. Концентрации c0(t) = c(t, 0), cg(t) = c(t,£) быстро стабилизируются, но, как видно в масштабе часа, это локально: происходит смена переходного режима всплеска на медленный тренд. Сравнимое время стабилизации cg (по отношению к c0, когда на входе скачком создается давление в десятки Torr) объясняется большой водоро-допроницаемостью сплава и £ ^ 1. Высокая скорость переходных процессов ожидаема, поскольку характеристическое время диффузии £2/D (с коэффициентом D & 10-4 на расстояние £ ~ 0.02) составляет всего несколько секунд. Общая концентрация определяется суммой c + z. Квазиравновесные (сивертсовские) концентрации растворенного диффундирующего водорода c0;g(t) а ^p0,g (графики barc0,g) определяются по давлениям молекулярного водорода соотношениями ^sp0,g = (&) bc0 g: C0,g(t) = Г^р0, g. С учетом захвата следует оперировать суммами с0 ,g(t)+zmax и c0 ,g(t)+z0 ,g(t). По рис. 6 можно оценить, насколько рассогласование концентраций cg — Cg на выходе существенно больше входного C0 — c0. Лишь асимптотически (по мере приближения к истинному равновесию) происходит сближение этих величин. Конечно, с точки зрения газоразделения нас в первую очередь интересует проникающий поток. Но аппроксимация градиента концентрации dxc разностным отношением [cg(t) — c0(t)]/£ ухудшается переходом к квазиравновесным оценкам как на входе (60 > c0), так и на выходе (Cg < cg).
Проиллюстрируем последние рассуждения. Примем базовую модель за «начало отсчета» и упростим ее в предположении квазиравновесности приповерхностных концентраций и линейного распределения в объеме (приближение Ричардсона для потока проницаемости):
dQinout = ±SD(T) Iе
dt v ! dx
= ±SD(T)
x=0,g
cg(t) - co(t)
í
= ±SD(T)r
VPg(t) - VPo(t)
^ ^-1 ^ Фо/
= «РоДО Цп^Т ^ =
Решая численно систему двух дифференциальных уравнений для Ро/(*) с начальными данными ро(0) = рс0, р^(0) = 0, получаем приближение (пунктир на рис. 7). Достаточно интегрировать одно уравнение в силу
&(*) = — ИпК-и1 рРо(^) ^
= ИпК—1 [ро — Ро(*)] .
Результат неудовлетворителен, если сравнить с базовой моделью (рис. 1). Но ценою существенной вариации «истинных» значений параметров можно добиться хорошего приближения (сплошные модельные линии). Это одна из основных причин разброса оценок: читая в литературе об «измеренном» коэффициенте X, нужно следить, по какой модели его значение на самом деле вычислялось. Разумеется, каждая модель - всего лишь модель, но при сравнении мы за «начало отсчета» берем базовую модель, поскольку квазистационарное приближение является ее дальнейшим и существенным упрощением. «Истинный» коэффициент проницаемости (данные приведены на рис. 7) равен Ф = ^Г ^ 4.438 ■ 1015. При подгонке по квазиравновесностационар-ной модели получается 2.453 ■ 1015, т. е. почти в два раза заниженное значение. Сравнивая с отношением [со — с^]/[со — ~ 2 (рис. 6), заключаем, что значительно более точной была бы квазистационарная модель
роД*) = [а^1П'О;1£]-1Т — со(*)]. Но информации о граничных концентрациях нет, и «вынужденная сивертсовская» подстановка c0'g = Г^ занижает коэффициент проницаемости. Формально можно принять модель
dpo,g dt
±SDr
aVi
in,out
T n
Vpg(t) - Vpo(t)
Vi(t) = VinVo-1 [po - po(t)], n e (0,1), to ^ 0,
но об априорном множителе n нет информации (кроме n ^ 1 с ростом момента to).
Ясно, что по мере стремления к равновесию (не обязательно реализовывать подгонку с момента времени to = 0, лучше пропустить начальный переходный процесс) и с повышением температуры поправочный коэффициент П будет расти (n ^ 1). Эта динамика проиллюстрирована на рис. 8 (t0 = 5min). Пунктиру соответствует квазиравновесная модель с исходными данными на рис. 3, а сплошные линии получены подгонкой за счет искажения коэффициента проницаемости Ф.
concentration x 1e-20
concentration x 1e-20
bar c0
c0 400 C
q D = 1.27e-4 b = 7e-23 s = 9e-5 \bar p0 = 32.87 \mu = 9.5e+20 a = 0.01 zmax = 0.1\bar c = 2e+19 Г = 3.5e+19 rmax=r(1+0.1)=3.85e+19
Zi bar ci
1 z0
2.2 2
1. 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
time, min
10 15 20 25 30 35 40 45 50 55 60
time [min]
Рис. 5. Приповерхностные концентрации, 1 min Рис. 6. Приповерхностные концентрации, 60 min
pin, pout [Torr] Quasi-equilibrium model
pin, pout [Torr] Quasi-equilibrium model
0
05
D = 1.27e-4 b = 7e-23 s = 9e-5 \bar p0 = 32.87 \mu = 9.5e+20 Г = 3.5e+19
D = 1.32e-4, D0 = 3.25e-4 b = 11e-23 s = 4e-5 Г = 1.86e+19
0 4 8 12 16 20 24 28 32 36 40 44 43 52 56 60
time [min]
500 C
"---.
Ф = Dr = 6.47e+15 ' ' —1 — ---
\bar p0 = 28.5 ___ —' _________
"""""__
---
Ф = Dr = 4.32e+15
0 4 8 12 16 20 24 28 32 36 40
time [min]
Рис. 7. Квазиравновесная модель, 400 °C
Рис. 8. Квазиравновесная модель, 500 ° C
400 C
32
28
20
6
2
4
Отметим также следующее обстоятельство. В квазиравновесную модель входит только комплекс параметров переноса Ф = ^Г = ^л/^в/6. При вариациях значений Ь, в, сохраняющих коэффициент проницаемости Ф, модельные графики давлений не изменятся. Для тонких мембран с большой водородопро-ницаемостью ситуация практически квазиста-ционарна (линейные распределения концентрации), но отлична от квазиравновесной (из-за отличий со,е от Со,е). С ростом времени квазистационары все ближе к квазиравновесным распределениям. Это означает, что обратная задача параметрической идентификации исходной базовой модели будет плохообусловле-на (в смысле слабой чувствительности модельных давлений к вариациям параметров, сохраняющих значение Ф). Как показано на рис. 9
(параметры указаны в соответствии с аппроксимацией при £ = 0.02 ст, рис. 1), увеличение толщины мембраны в разы кардинально не решит проблему: качественно картина сохраняется, а изменение уровней давлений происходит в значительно меньшем масштабе. Для «раскачки» комплекса ^Г целесообразно провести эксперименты насыщения-дегазации (хотя бы при нескольких температурах) более массивного образца. Это дало бы информацию о коэффициенте растворимости Г, что позволит «отделить» поверхностные параметры Ь, в от объемного Если емкость ловушек значительна, то следует фиксировать начальный всплеск десорбции диффузионно подвижного водорода и последующий длительный этап высвобождения из ловушек (для оценки Г и
Гтах (-^шах)) .
0 10 20 30 40 50 60 70 80
time [min]
Рис. 9. Зависимость давлений от толщины t
c(t,x) x 1e-20, [c(t,x) + z(t,x)] x 1e-20
thickness [cm]
Рис. 11. Распределение концентрации ~ 60 шш
Перейдем к распределениям концентрации в объеме, по толщине мембраны (ж € [0,£]) -рис. 10-12. Очень быстро устанавливается квазистационарный режим проницаемости, характеризуемый практически линейным распределением с(*, ж). Характеристическое время диффузии - несколько секунд. Учет захвата приводит к вертикальному сдвигу графиков, спадающему к ж = £ лишь на этапе насыщения ловушек (при г = гтах сдвиг параллельный). В течение полминуты градиент растет (по абсолютному значению), а затем концентрация на входе начинает падать, на выходе - продолжает расти. В середине пластины формируется точка перегиба, что отчетливо видно при относительно больших временах (рис. 11). Рис. 12 иллюстрирует динамику различий квазистационарных и квазиравновесных (пунктиром) распределений. Под
thickness [cm]
Рис. 10. Распределение концентрации ~ 30 шш
c(t,x) x 1e-20, cR(t,x) x 1e-20
thickness [cm]
Рис. 12. Сравнение с приближением Ричардсона
квазистационаром мы понимаем практически линейный по ж профиль концентрации с(*, ж) с насыщенными ловушками (как при установлении стационарной проницаемости), который относительно медленно меняется со временем. Квазиравновесным считаем такой квазистационар, который характеризуется заменой Со/(*) на «сивертсовские» концентрации
Со/ = Гу/роД*) (г = 2тах) .
Завершим обсуждение численного моделирования анализом динамики потоков (их плотности). Вследствие скачкообразного напуска водорода (под достаточно большим давлением) на входе происходит быстрый переходный процесс (рис. 13). Резкий всплеск (плотности) входного диффузионного потока 70(*) = —0) (материал вначале «пуст») сменяется спадом и стабилизацией. На выходе по-
ток диффузии J¿(í) = —Ddxc(í,¿) практически совпадает с десорбционным, поскольку давление в выходном объеме еще пренебрежимо мало для регистрации заметной ресорбции ^(T)s(T)pe(í)- «Слипание» J0(t), J*(í) говорит о том, что наступил квазистационарный режим: сколько вошло атомов водорода в единицу времени в образец, столько и вышло, распределение в объеме линейное. Поток водоро-допроницаемости в приближении Ричардсона —D [c,(í) — co(t)] Д = —Dr [^¡(í) — Vp0(t) обозначен как Jr. Локальная стабилизация потоков в масштабе часа (рис. 14) выглядит лишь переходом к длительным монотонным трендам. Потоки не отражены в силу
баланса — bc0/ = ±J0/. Заметна ошиб-
ка превышения Jr над «истинным» уровнем проницаемости J0 = J¿. Качественные рассуждения понятны из физических соображений, здесь важны количественные характеристики.
flux [cm-2 sec-1] x 1e-18
6 7
time [sec]
Рис. 13. Динамика потоков ~ 10 sec
flux [cm-2 sec-1] x 1e-1í
400 C D = 1.27e-4 J0
b = 7e-23 _ _ _ Jl
s = 9e-5 Jr
\bar po = 32.87 bc02
bco2 \mu = 9.5e+20 — — bcl2
a = 0.01, l = 0.019
zmax = 0.1\bar c 2e+19
Г = 3.5e+19
rmaX=r(1+0.1)=3.85e+19
Jr
be2
J0. Jl -----------------
0 8 16 24 32 40 48 56
time [min]
Рис. 14. Динамика потоков ~ 60 min
Заключение
Модель ориентирована на прикладную задачу выбора материалов для мембранной технологии выделения особо чистого водорода. Физико-технический характер задачи предполагает оценку основных интегральных показателей водородопроницаемости. Но параметров не должно быть слишком много, учитывая ограниченную информативность эксперимента. Введен параметр максимальной емкости стока -тах без детализации многообразия ловушек. Мембраны в установках газоразделения тонкие, материал достаточно однороден, с высокой водородопроницаемостью, так что захват носит характер малой поправки.
Результаты вычислительных экспериментов демонстрируют способность принятой модели удовлетворительно аппроксимировать экспериментальные кривые (см. [16]). Результаты моделирования соответствуют физическим представлениям качественного характера, но позволяют дополнить их информацией о «производных» выходных данных по отношению к вариациям параметров водородопро-ницаемости материала. Однозначность определения значений параметров нельзя гарантировать без дополнительных физически оправданных предположений о монотонном росте с увеличением температуры в рассматриваемом «рабочем» диапазоне.
При относительно высоких давлениях напуска комплекс параметров ВГ = В^ф является определяющим для тонких мембран с высокой водородопроницаемостью, что требует разработки специального математического аппарата при решении обратной задачи параметрической идентификации. В частности, с высокой точностью невозможно оценить емкость ловушек, чей вклад в общую растворимость в объеме материала может оказаться значительным. Целесообразно, наряду с проницаемостью, проводить оценку растворимости методом насыщения-дегазации. Это позволит «разделить» влияние объемных и поверхностных процессов.
Показано, что для очень тонких мембран с высокой водородопроницаемостью пользоваться приближением Ричардсона («сиверт-совские» концентрации в приповерхностном объеме) можно лишь с целью оценки порядка коэффициента проницаемости. Представляется более целесообразной разработка предложенной квазистационарной модели с контролем точности вычислений в рамках базовой модели (требующей значительно больших вычислительных ресурсов).
0
2
3
4
5
Работа выполнена при поддержке РФФИ (грант № 15-01-00744).
Литература
1. Водород в металлах/ Ред. Г. Алефельд, И. Фелькль. М.: Мир, 1981. Т. 1, 506 c. T. 2, 430 с.
2. Взаимодействие водорода с металлами / Ред. А. П. Захаров. М.: Наука, 1987. 296 с.
3. Писарев А. А., Цветков И. В., Марен-ков Е. Д., Ярко С. С. Проницаемость водорода через металлы. М.: МИФИ, 2008. 144 с.
4. Черданцев Ю. П., Чернов И. П., Тюрин Ю. И. Методы исследования систем металл - водород. Томск: ТПУ, 2008. 286 с.
5. Изотопы водорода. Фундаментальные и прикладные исследования/ Ред. А.А. Юхимчук. Саров: РФЯЦ-ВНИИЭФ, 2009. 697 с.
6. Varin R.A., Czujko T., Wronski Z.S. Nanomaterials for solid state hydrogen storage. Springer, New York, 2009. 338 p. doi:10.1007/978-0-387-77712-2
7. The hydrogen economy / Eds. M. Ball, M. Wietschel. Cambridge Univ. Press, 2009. 646 p.
8. Handbook of hydrogen storage: new materials for future energy storage / Ed. M. Hirscher. Wiley-VCH, 2010. 353 p.
9. Gabis I. E. The method of concentration pulses for studying hydrogen transport in solids // Technical Physics. 1999. Vol. 44, N 1. P. 90-94. doi:10.1134/1.1259257
10. 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. P. 5818-5851. doi:10.1016/j.ijhydene.2014.01.158
11. Indeitsev D. A., Semenov B. N. About a model of structure-phase transfomations under hydrogen influence // Acta Mechanica. 2008. Vol. 195. P. 295-304. doi:10.1007/s00707-007-0568-z
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. P. 90609069. doi:10.1016/j.ijhydene.2010.05.092
13. 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. 37763791. doi:10.1016/j.apm.2008.12.018
14. 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, N 1. P. 1295-1305. doi:10.1016/j.ijhydene.2010.07.099
15. 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.
16. Kojakhmetov S., Sidorov N., Piven V., Sipatov I., Gabis I. and Arinov B. Alloys based on group 5 metals for hydrogen purification membranes // Journal of Alloys and Compounds. 2015. doi:10.1016/j.jallcom.2015.01.242
17. Даркен Л. С., Гурри Р. В. Физическая химия металлов. М.: Металлург-издат, 1960. 585 c. doi:10.1002/ange.19540660409
18. Бокштейн Б. С. Диффузия в металлах. М.: Металлургия, 1978. 248 c.
19. Алимов В. Н., Буснюк А. О., Ноткин М. Е., Лившиц А. И. Перенос водорода металлами 5-й группы: достижение максимальной плотности потока сквозь ванадиевую мембрану // Письма в ЖТФ. 2014. Т. 40, вып. 5. C. 88-94.
20. Dolan M. D. Non-Pd BCC alloy membranes for industrial hydrogen separation // Journal of Membrane Science. 2010. Vol. 362. P. 12-28. doi:10.1016/j.memsci.2010.06.068
Поступила в редакцию 31.03.2015
References
1. Vodorod v metallakh [Hydrogen in metals] / Eds. G. Alefel'd, J. Fel'kl'. Moscow: Mir, 1981. Vol. 1, 506 p. Vol. 2, 403 p.
2. Vzaimodeistvie vodoroda s metallami [Interactions of hydrogen with metals] / Ed. A. P. Zakharov. Moscow: Nauka, 1987. 296 p.
3. 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.
4. 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.
5. Izotopy vodoroda. Fundamental'nye i prikladnye issledovaniya [Hydrogen isotopes. Fundamental and applied studies] / Ed. A.A. Yukhimchuk. Sarov: RFYaTs-VNIIEF, 2009. 697 p.
6. Varin R.A., Czujko T., Wronski Z.S. Nanomaterials for solid state hydrogen storage. Springer, New York, 2009. 338 p. doi:10.1007/978-0-387-77712-2
®
7. The hydrogen economy. Eds. M. Ball, M. Wietschel. Cambridge Univ. Press, 2009. 646 p.
8. Handbook of hydrogen storage: new materials for future energy storage / Ed. M. Hirscher. Wiley-VCH, 2010. 353 p.
9. Gabis I. E. The method of concentration pulses for studying hydrogen transport in solids. Technical Physics. 1999. Vol. 44, N 1. P. 90-94. doi:10.1134/1.1259257
10. 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. P. 5818-5851. doi:10.1016/j.ijhydene.2014.01.158
11. Indeitsev D. A., Semenov B. N. About a model of structure-phase transfomations under hydrogen influence. Acta Mechanica. 2008. Vol. 195. P. 295304. doi:10.1007/s00707-007-0568-z
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. P. 90609069. doi:10.1016/j.ijhydene.2010.05.092
13. 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. 37763791. doi:10.1016/j.apm.2008.12.018
14. Zaika Yu. V., Bormatova E. P. Parametric identification of a hydrogen permeability model by
СВЕДЕНИЯ ОБ АВТОРАХ:
Заика Юрий Васильевич
зав. лаб. моделирования природно-технических систем, д. ф.-м. н. Институт прикладных математических исследований Карельского научного центра РАН ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: zaika@krc.karelia.ru тел.: (8142) 766312
Родченкова Наталья Ивановна
научный сотрудник, к. ф.-м. н. Институт прикладных математических исследований Карельского научного центра РАН ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: nirodchenkova@krc.karelia.ru тел.: (8142) 766312
delay times and conjugate equations. International Journal of Hydrogen Energy. 2011. Vol. 36, N 1. P. 1295-1305. doi:10.1016/j.ijhydene.2010.07.099
15. 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.
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. doi:10.1016/j.jallcom.2015.01.242
17. Darken L. S., Gurri R. V. Fizicheskaya khimiya metallov [Physical chemistry of metals]. Moscow: Metallurgizdat, 1960. 585 p. doi:10.1002/ange.19540660409
18. Bokshtein B. S. Diffuziya v metallakh [Diffusion in metals]. Moscow: Metallurgiya, 1978. 248 p.
19. 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, N 3. P. 228-230. doi:10.1134/S1063785014030031
20. Dolan M. D. Non-Pd BCC alloy membranes for industrial hydrogen separation. Journal of Membrane Science. 2010. Vol. 362. P. 12-28. doi:10.1016/j.memsci.2010.06.068
Received March 31, 2015
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) 766312
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
