Параметрическая идентификация модели водородопроницаемости с динамическими граничными условиями Текст научной статьи по специальности «Математика»

Труды Карельского научного центра РАН № 3. 2010. С. 30-44

УДК 519.6: 539.2

ПАРАМЕТРИЧЕСКАЯ ИДЕНТИФИКАЦИЯ МОДЕЛИ ВОДОРОДОПРОНИЦАЕМОСТИ С ДИНАМИЧЕСКИМИ ГРАНИЧНЫМИ УСЛОВИЯМИ

Ю. В. Заика1, Е. П. Борматова2

1 Институт прикладных математических исследований Карельского научного центра РАН

2Государственное образовательное учреждение высшего профессионального образования Петрозаводский государственный университет

Рассматривается обратная задача определения параметров нелинейной модели водородопроницаемости с динамическими граничными условиями. По известному из эксперимента выходному десорбционному потоку предложен алгоритм оценки параметров переноса: адсорбции, десорбции, растворения и диффузии.

К л ю ч е в ы е с л о в а : водородопроницаемость, математическое моделирование, оценка параметров.

Yu. V. Zaika, E. P. Bormatova. PARAMETRIC IDENTIFICATION OF HYDROGEN PERMEABILITY MODEL WITH DYNAMICAL BOUNDARY CONDITIONS

The inverse problem of parameters determination of nonlinear model of hydrogen permeability is considered. Using the outlet desorption flux known from experiment the algorithm of transition parameters (adsorption, desorption, solution and diffusion) estimation is proposed.

K e y w o r d s : hydrogen permeability, mathematical modelling, parameters estimation.

Введение

Исследования в области водородного материаловедения инициированы перспективами водородной энергетики и проблемой защиты конструкционных материалов от водородной коррозии [Водород в металлах, 1981]. Экспериментальный метод проницаемости является классическим [Кунин и др., 1972]. Вместе с тем при обработке экспериментальных данных имеется разброс в оценках значений параметров моде-

лей. Одна из основных причин: обратные задачи математической физики характеризуются высокой чувствительностью к экспериментальным и вычислительным погрешностям. Кроме того, традиционный «метод подгонки» обычно не гарантирует единственности набора параметров, удовлетворительно аппроксимирующих экспериментальные кривые. Поэтому, во-первых, необходимо определить объем данных, достаточный для однозначной параметрической идентификации модели. Во-вторых, алгоритм оценки

параметров должен быть корректным (в частности, не содержащим дифференцирования зашумленных измерений).

Будем ориентироваться на следующий вариант эксперимента. Пластина толщины i из исследуемого материала (металла или сплава) является перегородкой вакуумной камеры. При заданной температуре T = T образца с входной стороны скачкообразно создается давление p0 = const молекулярного водорода. С выходной стороны производится вакуумирование и с помощью масс-спектрометра определяется де-сорбционный поток J = J(t), t - время. Потоки подразумеваем отнесенными к единице площади, т. е. речь идет о плотностях. По достижении стационара J (t) ~ J = const, t > t*, на входе резко увеличиваем давление до значения p+ > р0 и дожидаемся последующего установления десорбции при t > t > t*. Такой вариант предпочтительнее двух классических экспериментов, поскольку нет необходимости в повторной дегазации образца и «старт» второго этапа происходит не с нулевого начального распределения атомарного водорода в пластине, а с предшествующего стационара (разнообразие повышает информативность). Выход на стационар носит асимптотический характер. Но t* и t* не следует выбирать слишком большими, чтобы переходные процессы «не потерялись» на фоне стационарных.

Давления и концентрацию c(t, х) растворенного водорода считаем относительно малыми, D D(c). Поскольку речь идет об обратных задачах и экспериментальные погрешности в лучшем случае оцениваются в 10-20%, то авторы старались ограничиться минимальным набором коэффициентов, подлежащих определению.

Математическая модель

Рассмотрим следующую модель с поверхностной десорбцией [Габис и др., 1987]:

дс д 2С

?- = D^-2, (t, х)е (0,t.)X(0,i), dt дх

c(0,х) = 0, хе [0,i], P = }Jsp0,

q0(t) = P - bql(t) + Ddc дх

q( (t) = -bq2(t) - D— дх

(1)

(2)

, tе [0,t*],

ф) = ф,0) = gq0(t), се(?) = с(г,I) = gqt(?). (3)

Здесь ц - кинетический коэффициент, s(T) -коэффициент прилипания, Ь(Т) - коэффициент десорбции, g (Т) - параметр соответствия поверхностной и объемной концентраций. Зависимость параметров диффузии и десорбции от температуры Т считаем аррениусовской:

Б = Д, ехр{-Ев /[ЯТ]}, Ь = Ьо ехр{-Еь /[ЯТ]} .

При необходимости допустимы другие зависимости от Т. Температурная зависимость /л формально учтена в коэффициенте s(T). Далее это непринципиально, поскольку в течение эксперимента температура постоянна.

Смысл граничных условий (2): дисбаланс потоков адсорбции, десорбции и диффузии идет на накопление атомов водорода на поверхности. Вакуумная система достаточно мощная, чтобы пренебречь ресорбцией: 1&р1 ~ 0 . Более точная

модель растворения на поверхности:

к + (Т) С0,( ^)[1 - qоЛ ^)qm1ax] -- к-(Т) qoла)[1 - ^ (t)Ст!х] = ±Б(Т) Сх\й1.

Но когда диффузия значительно медленнее растворения и концентрации малы, получаем условие быстрой растворимости с0 , ~ gq0 ,, где

g = к- / к + . Если поверхность изотропна (в смысле Ек- « Ек +), то параметр g слабо зависит

от Т . По измерениям известна функция 3^) = у (t) = Ьq<2(t). Кривые проницаемости

имеют стандартный 8-образный вид кривой насыщения. Для второго этапа эксперимента (t е [£,, t*]) начальное распределение с(и, х) -стационар, а давление р0 заменяем на . Требуется определить значения Б, g, Ь, s при температуре эксперимента Т = Т . В алгоритме идентификации используются интегральные соотношения, устойчивые к помехам.

Численное решение прямой задачи. Вначале входная сторона пластины испытывает ударную нагрузку из-за скачка давления. В связи с этим (и нелинейностью задачи) схема вычислений носит неявный и итерационный характер. Для аппроксимации уравнения (1) используем шеститочечный двухслойный шаблон по схеме Кранка-Николсон. Перейдем к граничным условиям. В стандартных обозначениях

х=i

ч0,1 * Яол(ґк), ск * с(ґк,хі) , здесь ґ0 = О, Х0 = о, хп = І, т = Аґ, А = Ах . Из начальных данных определяем д00і = сО = О, О < і < п . На каждом слое по времени аппроксимируем сх (ґк ,0) * [-3с0 + 4с\ - с\]/2к . Аналогичные выкладки при х = І в изложении опускаем. Заменив производную по ґ конечной разностью (дк0 - дк0-1)/т = (сО - сО-1)/(gт), находим сО = /0(с1, с2) как положительный корень квадратного уравнения. Значения ск, с2 предварительно подсчитываются по явной схеме. С теку-к к

щими со, сп решаем методом прогонки трехдиагональную систему линейных алгебраических уравнений и находим новые приближения ск, с2 (и остальные значения ск, 3 < і < п -1). Снова решаем квадратное уравнение относи-£

тельно с0 и повторяем вычисления до установления (обычно 2-3 итерации). Затем переходим к следующему слою по времени.

Результат счета контролировался выходом распределения с(ґ, х) на линейный стационар при соблюдении материальных балансов:

І ґ ґ I с(ґ, х) &х = -1 Бсх (ґ,0) & +1 Бсх (ґ, І) & ,

0 0 0 ґґ

qt (ґ) = -1 ^ (ґ) & -1 Бсх (ґ, І) & . (4)

00

Для второго этапа (ґ є [ґ*, ґ*]) построения аналогичны. Численные эксперименты проводились в широком диапазоне: Б от 10-9 до 10-3 ст2/5 ; Ь от 10-20 до 10-5 ст2/ 5 ; g от 10-3 до 106 1/ст ; 5 от 10-10 до 10-2 ; р0 от 0,1 до 10 Тогг; /и = 1,46 1021 1/ст25Тогг . Толщина І варьировалась в пределах 2 10-2 - 2 10-1 ст , что соответствует экспериментальной практике. «Внутренние» параметры Б, Ь, g, 5 («внешние» І, р0, р+

задаются), порождающие кривые проницаемости, «забывались» и восстанавливались по излагаемой ниже методике. Затем анализировались качественные возможности алгоритма идентификации.

Вырожденные модели. Наряду с задачей (1)-(3), которую будем обозначать I, рассмотрим еще две. Задачу, в которой граничные условия (2), (3) заменены линейными

с0(0 = с(^0) = с0 > 0, с, (0 = с(^,) = 0, (5)

обозначаем II. Простейшая модель II не учитывает физико-химические процессы на поверхности, диффузия - единственный лимитирующий фактор. Плотность выходного потока атомов водорода определяется как 3и ^) = -Бсх (^ ,). Решение

с(^ х) понимается как обобщенное из-за несогласованности краевых условий при t ^ +0 . Для t > 0 обычно пользуются представлением с(^ х) рядом Фурье. Примем с0 равной установившейся концентрации с0 в задаче I. Это соответствует

ситуации, когда поверхностные процессы значительно быстрее диффузии. Формально задача II -первого рода. Но при (см. (6))

= со = ^Ь-1( Р - 3{)

она связана с нелинейной задачей I.

Рассмотренную в [2а1ка, Вогта1юуа, 2008] модель с граничными условиями

Р - ~с02^) = -Б—

Эх

зш = Ьcїiґ), [Ь ] = -

Ьс2(ґ) = - Б

дх

(с объемной десорбцией) обозначим III. Поскольку в модели I 30, = Ьq0^t = Ьg-2с^,, то переход I ^ III понимаем в следующем смысле: накопление на поверхности несущественно и Ь = Ь / g2. Причина введения моделей II, III в том, что варьирование параметров I в указанных широких пределах может приводить к вырождению.

Равновесие и стационар

Равновесная пара (р, с ) определяется приравниванием к нулю производных:

Лр = Ьц2, с = ^д ^ с = Ут!р , у= g^[лІsЬF1.

Тем самым модель соответствует эксперименту в области с ~ д/р . Определение коэффициента равновесной растворимости у является более простой задачей. Когда в эксперименте насыщения-дегазации пластины нельзя пренебречь поверхностной концентрацией, необходима следующая корректировка. Без ограничения общности поверхности единичны, торцами пренебрегаем. Количество атомов Н после насыщения

х=0

х=І

4

5

численно равно с, + 2Ц = (g, + 2)д - известное число Q после дегазации. Тогда из /лр = Ьц 2 после подстановки д = Q /(g, + 2) определяем комплекс

~ = (g, + 2)7 Л / Ь , Q = у^р .

При g, >> 1 можно считать, что в равновесии поверхностная концентрация несущественна на фоне объемной и ~ , у . Когда известен коэф-

фициент у или ~, исключается один из параметров поверхностных процессов. Но другие эксперименты характеризуются своими погрешностями и требуют затрат. Постараемся обойтись без предварительного определения или ~.

Стационарное распределение водорода линейно. В модели I:

д0, = 0, с( = 0 ^ с(£,,х) = -В 13х + с0,

с» = ^Ь-1(Р - 3).

(6)

Из 3 = Ьд,2 находим значение с, = gqf и приравниваем к с(£«,,):

3-VР-3 + 7,3^В)-1 = 0 (Р = лр>). (7)

Из двух таких уравнений для давлений р0 и р+ исключаем л:

:3/2 Л

Л

J+ J z2 + 2 J+ J z + 2 J+ J

p0+ Pg p0+ P0 p0+ Pg

= 0,

z =-

(8)

gВ В

Корни разного знака, комплекс параметров г определяется однозначно. Затем из (7) находим 5 . Целесообразно предварительно уравнение (8) разделить на 3 и перейти к безразмерной переменной модели III

3-V р - 3+1Ы1/2ю-1 = 0

и по информации {р0, р+, 3, 3 + } однозначно определяются величины 5 и Ь 1/2,В-1. С учетом Ь = Ь / g2 стационарные уровни 3 в моделях I и III одинаковы при фиксированных Р . Параметр

Ь = Ь / g2 в силу 3 0,, = Ьд02, = bg-2с2, имеет

смысл объемного коэффициента десорбции. При фиксированных В, 5 и Р >> 3 из соотно-

шения (7) получаем lg J = -(lg b )/2 + const, т. е. линейную зависимость в логарифмических координатах.

Время запаздывания

Рассмотрим функцию S (t), равную интегралу на отрезке времени [0, t] от плотности десорбции J(т). Это количество атомов H, десорбировавшихся с единичной площадки при х = t за время t. График этой выпуклой функции имеет наклонную асимптоту. Точка пересечения асимптоты с осью t называется временем запаздывания t0. Уравнение касательной к графику S (t) в точке (t., S.):

S = S. + J(t.) • (t -1.) (S = J).

Следовательно, значение t0 e R+ достаточно точно вычисляется как t0 = t. - S. / J , S. = S(t.). Здесь t. - время установления выходного потока: J(t) ~ J = const, t > t.. Для краевой задачи II (c0(t) = ~0, ct (t) = 0, Ju = -Bcx (t, t)) имеется явное выражение t0 = t2/6 B [Кунин и др., 1972]. Важно, что t0 не зависит от ~0.

Замечание. Обычно считают c0(t) = c , где c ~ p01/2 - равновесная с давлением р0 концентрация. В контексте задачи I полагаем ~0 = c0:

быстро устанавливаются поверхностные процессы и лимитирует диффузия. Поскольку в равновесии P = bq2, а в стационаре P - bq02 = J, то c > c0. Кроме того,

J = bq2 ^ c2 = c02 + ct2, J = Jt < Jп .

Когда поток атомов H на входную поверхность значительно превосходит пропускную способность мембраны (P >> J), имеем c0 « c . Наконец, из формулы (6) следует c0 = ct + JtB-1, и если ct << c0, то c0 ~ JtB-1 (в задаче II

~ = 3n tB-1).

Будем говорить, что модель I вырождается в II (I ^ II), если при заданных B, b, s, g и

~0 = c0 = ^b-1(P - J) ( J = Jj ) плотности де-

сорбции Jt (t) = J(t) = bqt2 (t) и Ju (t) = -Bcx (t, t) равны как функции времени. Совпадение подра-

зумевается с погрешностью много меньшей экспериментальной. Аналогично интерпретируем обозначения III ^ II, I ^ III. В силу Ь = Ь / g2 значения с0 и 3 в I такие же, как и в III.

Критерии вырождения. Так как функция £ (^) строго выпукла, но известна приближенно, то в качестве критерия вырождения II можно взять соотношения И ^ = ,2/6В и с0 и ЛВ-1 (с, << с0). Они должны выполняться

совместно. Сравним проникающие потоки. Поскольку в задаче I граничные концентрации равны {с0, с, }, а в П-{ с0,0}, то справедлива

оценка 3г < 3п . Значение 3П = г-1^ //5р0 - 31

определяется из с0 = ^Ь-1 (Р - 31) = 111 ,В-1. Кроме того:

31 = Ьде ^ с, = ^Ь 3 ^ с0/ с, = X,

хЧР3- -1 > 1 (23 <Р).

В обозначениях х, г, I = л/3 уравнение стационара (7) запишется более компактно: 1г +1 = х. Параметр х дополнительно характеризует зону вырождения I ^ II в смысле с, << с0 ~ х >> 1.

Приведем комментарии качественного характера. В модели I концентрации д0, д, могут оказаться сравнимыми (Xй 1). При фиксированных В, g определяющим является произведение ЬР - чем оно меньше, тем «дальше» модели

I и II. Если значительно уменьшить падающий на входную поверхность поток Р = лр0, то проникающий поток 3 будет очень мал и с, » с0. Если уменьшить Ь , то на выходе десорбция будет слабой, что способствует росту с, (с, (?) монотонно растет до уровня с,). Обратно, если Ь велик, то на выходной стороне активная десорбция понижает концентрацию с, и в пределе получаем задачу II с с, (?) = 0. Рост ЬР влечет уменьшение времени установления концентрации: с0 (0) = 0 ^ с0. Когда модели перекрываются (I ^ II), большие (но не слишком) вариации 5, Ь практически не влияют на запаздывание t0 И ,2/6В . В этом и проявляется некорректность обратной задачи. Подчеркнем, что вариа-

ции b, sp0, сохраняющие их произведение, практически не меняют t0j .

Приведенные рассуждения дают основание наряду с временем запаздывания ввести в рассмотрение диффузионное и поверхностное времена: £ = t2 /2B и ц = 2/VbP . Если ц >> £, то лимитирующим фактором являются поверхностные процессы. При г/ << £ лимитирует диффузия (вырождение I ^ II). Разность ~ = tG - tGJ (tGJ = £ / 3) показывает, насколько поверхностные процессы увеличивают запаздывание. Следовательно, t и г) должны быть связаны монотонной зависимостью. Удобнее сравнивать относительные величины: St0 = (tG - tGJ )/£ и п/£. Коэффициент растворения g является связующим между поверхностью и объемом и играет особую роль.

Замечание. Помимо качественных рассуждений и соображений размерности (£~ t2/B,

П ~ 1/VbP ) приведем обоснование множителей.

Рассмотрим стационар в модели II:

c0(t) = ~0, ct (t) = 0, c(t, х) = (t - х)~0 /1, t > t..

Обозначим через v (х) «среднюю скорость стационарного переноса» в сечении х e (0, t). Тогда J = V (х)c(t., х) = B~0/1, dt = ёх / V (х) = (t - х) ёх / B

и после интегрирования на [0, t] получаем характеристическое время т. = t2/2B . Это значение возникает и при вероятностном анализе диффузии [Зельдович, Мышкис, 1973, гл. V]. Обратно. Пусть диффузия медленная на фоне поверхностных процессов. Рассмотрим установление

qo = P - bqG, qo(0) = G ^ q0 = л/PTb

(Bcx(t,0) « 0). Интегрируя, получаем

4P -4bq0(t) = exp {- 24bPt\([P +4bq0(t)),

1 ~4bP-1/2q0(t) < 2exp{-24bPt}.

При t = т. = N/4bP , N = 2 , левая часть неравенства, характеризующая отклонение от установления на поверхности, меньше 3,7 %. Можно было бы уточнить т., заменив в граничном условии функцию - Bcx (t,0) на J :

90

= Р - Ьд02 - 3 , д0(0) = 0 ^ д0 = V(Р -3)/Ь

^о( t) = - В

дс_

д х

((0 = 3 ),

> 0.

и материального баланса:

312

13В1хСх = |_Р0 -13 Л, =[0 + £,]3 - &

0 0 0 2В

3tо = 0 + t* - 3 = 0 + tо.

Отсюда 0 = 2^. Здесь у 3, &., t.., 10, 0 подразумевается индекс II.

Остановимся на свойствах величины 0 = 01, характеризующей жесткость задачи переноса. Разделим 0 на два слагаемых: 0 = 011 + 0. При t = и из баланса (4) с учетом -Всх (^0) = [.^,(0 - 3] + 3, представления (6) и условия Всх (t,,) = - 3 ^) - д, ^) имеем:

— — _ <• I

3 0 + 3t, - & . - д, = Jо с (^, х) Сх =

т. = 2/^Ь(Р-3) . Но стационарное значение 3 зависит от всех параметров модели, а поскольку 3 < Р / 2, то порядок оценки т. не изменится.

Выясним теперь механизм влияния коэффициента быстрой растворимости g на свойства решений задачи I. Обозначим

= со,

3,2 _

2 В

= соI - 3£.

Здесь & - интеграл от 3(т) = Ьд,2(т), те [0,и]. Подставляя выражения

со = с, + 3/В , с, = ёц, = ,

после деления на 3 для времени запаздывания t0 = и - & / 3 получаем соотношение

0 = 3-1 Гt’[^»(т) - 3] Ст

0

t* - ^ .С ро(т) Ст

Поясним смысл введения всплеск-времени 0 (по аналогии с диффузионным и поверхностным временами). На начальном этапе входной диффузионный поток ^0^) испытывает всплеск и

затем _Р0 ^ 3 . Продолжительность всплеска относительно невелика (<< и), но интеграл может оказаться значительной величиной. Следовательно, 0 - интегральная мера продолжительности этого процесса. Геометрически 0- расстояние от нуля до точки пересечения с осью времени асимптоты графика количества атомов Н, проникших в объем пластины с единичной поверхности. Аналогия с временем запаздывания. Когда модель I вырождается в II, значение 0 (как и ^) можно вычислить аналитически:

0 = 011 = ,2 / 3В = 2^1. Это следует из

~0 = с0 = 3,/В , t0 = ,2 /6В, 3 = -Всх(^,)

^+0=—+4^+1, I = 41 2В 14Ь

(9)

-1/2

?0 +0 = (g -1 + , )[ 3Ь ]

Аналогичное выражение для модели III:

-1/2

tg + 0 = , [ 3 Ь ]

Критерии вырождения £ << т], £>>п асимптотические. Как указать «шкалу соизмеримости» в конкретной задаче, когда диапазон допустимых значений параметров относительно мал (не 5-10 порядков)? Запишем соотношение (9) в переменных

£ = ,2 /2В, п = 2/4ьР : tg + 0 = £ + Оп,

~ + 0 = Оп, О = ^[pJ-1(g£ +1)/2.

Переход I ^ II характеризуется Оп ^ 0 (Оп<<£), О = О(В,Ь,g,5). Для оценки влияния поверхностных процессов на водородопро-ницаемость необходимо исследовать окрестность соизмеримости Оп ~ £ .

Для дальнейшего анализа удобно преобразовать множитель О . Введем безразмерный параметр у = -\[&,, характеризующий соизмеримость емкостей поверхности и объема пластины (в равновесии в объемном столбике под единичной площадкой находится атомов водорода). Воспользуемся выражением g, / 3 из уравнения стационара:

1

О = -| у +

21 V

1V 4,Р3 -1

= М| у +

1

М =

л/Р

л/р - 3-3

> 1, Оп=м\у+— \у[£ц

^ t0 + 0 = t0 + 0 - £ = М| у + -1 1Л/£п .

0

В знаменателе дроби в М под знаком корня выполняется л/Р - 3 -43 = 3г > 0,

л/р - 3-43=4Р-43+0(43),

3 = 3(В,Ь,g,5) < Р/2.

Основной вариант. Рассмотрим условие Р >> 3 , когда падающий на вход поток значительно превосходит пропускную способность материала. Этого можно добиться за счет {р0,,}. Уменьшение р0 увеличивает время эксперимента и снижает точность измерений. Поэтому понятно стремление увеличить Р . Но желательно не доводить дело до 4Р >> I, поскольку растет некорректность задачи оценивания Ь, g . Вместе с тем, область эксперимента

4Р >> I (проверяется после определения 5 по двум стационарам) позволяет заметить качественные закономерности, которые в более широком диапазоне подтверждаются численным анализом. Итак, при М « 1 получаем

а = ~ +0 = \у + -'4£п, у = 44~1

V у) (10)

э &т = | = ^(g + 80 = .

где дt0 = ~/£, 80 = 0/£, ~ = tg - tgI,

tgI = I2 /6В, 0 = 0 - 011 = 0 - 2/01. Отсюда ясно,

что величины а и 8(7, характеризующие отклонение модели I от II по входному всплеску и запаздыванию, зависят лишь от седлового параметра g и соотношения времен £ = ,2/2В,

п = 2/л/ЬР ( 7 - от среднего геометрического времени, 87 - от их отношения). Выбором , можно влиять на коэффициент усиления (у + 1/у). Для модели III а = кЛ/£п точка перегиба отсутствует.

Результаты численного моделирования. Для фиксированных значений Р, В (в указанных выше границах) рассматривались значения уе [10-2,102] и варьированием параметра десорбции Ь определялся примерный диапазон соизмеримости Оп/£е [10-2,102]. Величины 0 и 70 характеризуют «вход» и «выход». При

Р >> 3 (4Р >> I) знаем качественное поведение их суммы а. Желательно иметь информацию о них в отдельности.

Рис. 1. Зависимость ~ от 4£ц (перегиб при у = 1)

Рис. 2. Зависимость 0 от у/£/] (монотонность по у)

На рис. 1, 2 представлены зависимости ~, 0 от при различных уе [10 2, 102]. Графики 870, 80 от д/п/£ качественно аналогичны (перегиб по у и его отсутствие). Для определенности фиксированы значения р0 = 0,1, , = 0,01,

В = 10-6 , их вариации не приводят к качественным изменениям графиков. Указанным во врезках значениям у соответствуют g = у2/, . Абсцисса х = 4£п дает значение комплекса Ь5 = 4£2/(х4Лр0), определяющего время запаздывания ^ . Соизмеримость времен Оп ~ £ иллюстрируется таблицей.

Соизмеримость времен Оп ~ £

у 0,01 0,1 1 10 100

Оп/£ 0,3-100 0,3-70 0,1-100 0,3-200 0,3-100

ю1

ю'1

ю'2

«Г5 ю2 ю-1 10° ю’ юг

Рис. 3. «Производная входа по выходу»:

8Э/870 от 7п/£

Рис. 3 демонстрирует почти постоянство «производной входа по выходу», т. е. слабую зависимость 80/870 от поверхностного времени. Управляющим параметром является в основном g . При переходе через значение у = 1 перегиб происходит только у 70 , при

этом 0 растет монотонно с ростом у до значения «насыщения» 0III (для малых у может оказаться 0< 0). Формула (10) получена в предположении 3 << Р , но вычисления показали, что положение точки перегиба не меняется в широком диапазоне параметров (рис. 4, 5). Рис. 6 показывает динамику перехода I ^ III с ростом g .

о’

0°

0-’

10'2 ю'1 10° ю1 ю:

Рис. 4. Сравнительное поведение времен запаздывания, 3 << Р

Рис. 5. Сравнительное поведение времен запаздывания, 3 » Р

х ю4

Рис. 6. Динамика перехода модели I в III с ростом g

Остановимся подробнее на зависимости показателя пропускной способности пластины 3 / Р от параметров моделей I и III. Численно установлено, что при фиксированном £ в модели III приращение времени запаздывания ~ш = 70я - tgI и стационарная водородопроницаемость 3 / Р зависят только от произведения ЬР. Обозначим

Vш = [ЬР]12. Величина УШ имеет размерность скорости и характеризует динамику поверхностных процессов при быстром насыщении поверхности.

Рис. 7. Зависимость 70 от -\1ЬР (модель III)

Рис. 7 демонстрирует уменьшение 7дш с рос-

ътШ -ГГ

том поверхностной скорости V . При VIII = Vя = 2 В / , (VII - «средняя скорость» диффузионного переноса в модели II) с относительно малой погрешностью справедливо ~ «£

(£ = ,2 /2В - диффузионное время).

Рис. 8. Зависимость 3 /Р от \ЬР (модель III)

Рис. 8 показывает, как вместе с 7дш уменьшается проницаемость 3 / Р . При VIII << VИ показатель 3 / Р практически не зависит от VIII и близок к максимальному значению 1/2. С рос-

ътШ 1 I п ■, тШ

том V зависимость 3 / Р от V становится линейной в логарифмических координатах.

В модели I при фиксированном £ отношение

3 / Р зависит только от скорости vI = VЬР / g , поскольку стационарные потоки в I и III совпадают. Характер этой зависимости уточняет рис. 9.

Рис. 9. Зависимость 3 / Р от *]ЬР / g (модель I)

Но время задержки ~ = ґ0 - ґІ0І определяется при

этом не только величиной V1, но и значением g. Поэтому при качественном исследовании модели I использовались не скорости, а времена. С ростом g происходит вырождение 1^ III

и точка с координатами (V11 ,£) попадает на график 70^1 )(рис. 10).

Алгоритм определения параметров

По классической кривой проницаемости нельзя однозначно сделать вывод о соизмеримости поверхностных процессов и диффузии. Если взять 3 и 70 = 7. - & / 3 из задачи I и принять в

II В = Вп = ,2 /67д, ~0 = 3, / Вп , то в II получим то же значение 3 . Близкими могут быть и переходные процессы. Например, при g ~10-2, ЬР ~10-7, В1 = 10-6 концентрация с0 больше

,3 / В1 почти в два раза и 7д » ,2 /2ВI. Однако

3(7) в моделях I и II (при указанных Вг1 и ~0)

совпадают в пределах 5 %. Поэтому и предлагается эксперимент с двумя давлениями и «стартом» на втором этапе не с нуля, а с достигнутого стационара.

Вначале продолжим рассуждения в рамках модели II (~0 = с0). Для второго этапа

(р0 ^ р+ , 7 е [7., 7*]) сделаем замену

сА (7, х) = с(7, х) - с(7*, х) и примем 7* за начало

отсчета времени. Здесь с(7*,х) = 3В-1(,-х),

3 = 3п, 7* = 7* , 7* = 7*. Получим ту же краевую

задачу для сА, только вместо с0,3и (7) будут

Ас> — сд сд и А^ц — 3ц 3ц ( 3ц — Р, ). Новое

время запаздывания вычисляем по формуле

А& А7

А7д = А7* --=3^, А& = Г А3II(т) Ст,

А3II д

А7* = 7 - 7*, А3Я = 3+ - 3П . (11)

В исходном времени интегрирование ведется по те [7*,7*] и А70 определяется длиной отрезка между точкой (7*, &*) и точкой пересечения асимптоты для & + А&(7) (7 > 7*) с горизонтальной прямой & = &*. В модели II А7д = 7д = ,2 /6В и не зависит от значений Ь, g, 5, рд, р+ .

Вырожденный случай. Напомним, что значения 5 и г = ,л/Ь(4В)-1 = Ь1/2, В-1 определяются из анализа стационаров. По известной из эксперимента плотности десорбции 3(7) = Ьд2(7) вычисляем 71 = 7* - &* /3 и

72 =А7* -А&*/ А3, т. е. в (11) вместо А3Я(7)

используем А3(7) = 3(7) - 3, А3 = 3 + - 3 . Если

7 совпадают с высокой точностью, то обоснованно полагаем, что эксперимент находится в зоне II. Ведь в модели I на втором этапе изменилось не только входное давление, но и начальные данные. Хотя теоретически времена запаздывания 7 могут совпасть

и в I при специально подобранных р0, р+. Из ti = ,2/6В находим оценку коэффициента диффузии В . Соотношение (7) дает значение комплекса g / 4Ь . Однозначно определяется коэффициент объемной десорбции Ь = Ь /g2.

Для разделения параметров g, Ь требуется дополнительная информация, например, значение равновесного коэффициента у (у). Если 7 И 72, но Р ~ 3 и (или) Р + ~ 3+, то можно перейти к нелинейной модели, считая полученные оценки В, Ь подлежащими уточнению.

Нелинейная модель. Пусть вычисленные 7

существенно различны в масштабе времени установления 7* . Обратную задачу решаем последовательно в соответствии с усложнением модели: II ^ III ^ I. Известны значения 7., 5 ,

г = Ь1/2, В-1, Р = Лр0 и грубые оценки В, Ь после обработки измерений 3(7) по линейной модели II.

1. В рамках нелинейной модели III справедливы соотношения

7д +0 = £ + ,[~3 ]-1/2 =£(1 + 2/[ г! ])=£[х+1],

7д = 71 = 7* - &*3 1, г.I +1 = X, (12)

г = ,~1/2В-1, х =л/Р3-1 -1 > 1, 0 = 0я(~, В). Имеем нижнюю грань

7д + 0 = ,2/ В = 2[7д0 + 0I ].

Решаем краевую задачу III с оценками В, Ь и вычисляем приближение 0(1). Заменив в уравнении (12) неизвестную величину 0 на 0(1), находим значения £(1) и В(1) = ,2/2£(1), Ь(1) = (В(1) г / ,)2. Повторяем итерации до установления, ориентируясь на невязку

13 - 3 (1,2)|.

2. Переходим к модели I, только когда существенна разница в переходном процессе 3(7) -3Ш(7) (поскольку 3: = 3Ш, см. рис. 6).

Известны 5 , Ь = Ь/g2, В = Ь1/2,/г. Осталось определить параметр g или Ь. Выбирая

0(1) = 0Ш , 0(1) = 2£/3 или даже 0(1) = 0 (когда очень велика разница между 3Ш и 3, g << 1), можно организовать итерационный процесс на основе соотношения (9). Но в конкретной задаче сразу ясно (см. рис. 6, 11) в каком направлении нужно изменять значения g.

Рис. 11. Идентификация моделей НИ по временам запаздывания, Ьш = Ь

Сопряженные уравнения

Выше при уточнении величины 0 использовался метод простой итерации. Можно реализовать и (квази)ньютоновский алгоритм, но это представляется громоздким. Основная проблема для сходимости - начальное приближение. Этот этап вызывает определенные трудности в связи с начальным всплеском входного потока Р0(7) = - Всх (7,0). Попытаемся сгладить эту неприятность, нацелившись на итерации интеграла от входной концентрации, которая монотонна и «глаже» (в силу д0 = Р - Ьд2 - Рд концентрация с0 интегрально зависит от Р0).

Воспользуемся техникой сопряженных уравнений [Марчук, 1992]. Определим сопряженное к (1) уравнение / = -В/хх . Для любого

решения (сопряженной функции) интегрированием по частям получаем

0=Я /[с- Всх] СхС=I (/с)17=д Сх -

(13)

0 0 Смысл такого выбора /(7, х) в том, чтобы после преобразований не появилось двойного интеграла. Остались величины, связанные с краевыми условиями (и измерениями). В простейшем случае (/ = 1) с учетом Л +1 = х получаем уравнение баланса (9):

7д + 0 = Л. +

0 2В л/Ь

+

1

л .х+1 _____

2В х-1 ' Л/Ь'

В тождество (13) можно подставить выражения с(7*, х) согласно (6) и с(7,1) = ^3 / Ь , с(0, х) = 0. Нет информации о с0 (7) и Р0(7) = -Всх (7,0). От одного из соответствующих слагаемых в (13) можно избавиться дополнительным ограничением /(7,0) = 0 или / (7,0) = 0 . Предпочтительнее исключить плотность входного потока Р0 (/(7,0) = 0). Значения 5, г считаем уже известными. Для варианта / = х / , (знаменатель для нормировки) с учетом граничного условия (2) (х = 0) получаем

с, 3,2 I

/ (Ь, В) = ^ - — + -^ + & (7*) +

4В&1/2(7*)

,4ь

В

I

2 3В

| сд(7) с) = 0, (15)

о

ад = & = |д 3(т) Ст, &1/2(7*) = &1/2 = |д I(т) Ст,

I ^л/3, сд = ^Ь-1( Р - 3).

По аналогии с 70 введем времена запаздывания для с0(7) и I(7):

__1 р7« — -^ <• 7*

7с0 = 7* - с0 |0 с0(т) Ст, 71/2 = 7* -7 . /(т) Ст.

Тогда уравнение (15) с подстановкой г7 +1 = х перепишется в форме

,2 7

[х + 2] + ^Ь = [х- 1])0 + 71/2 - ^>х. (16)

Подставим выражение л/Ь из уравнения (14):

£= 0[х- 1] + 7с>х 71/2

(£ = £2/2В, 0 = 0-0я =0-,2/3В)

(17)

Итак, диффузионное время £ представлено в виде £ = 71 + 72 - 73 (рис. 12). Рассмотрим асимптотику 72 - 73 с ростом g . Стационарные концентрации сд, (с0 > с,) удовлетворяют соотношениям

Р - с>2Ь / g2 = 3, с,2Ь / g2 = 3 э

сд/с, = х = VР3-1 - 1 > 1, 3 < Р/2.

В силу уравнения стационара х = гI +1 формально при g — +^ (г = ,4Ь /(4В) — +0) имеем х —— +1, с0 - с, —— +0, сд2,Ь /g2 — Р/2. Кон-

центрации сд, растут согласованно с g. Далее

большая входная концентрация на начальном этапе приводит к быстрому переносу к выходной стороне, но из-за большого g относительно малы д, и соответствующий десорбционный

отток. Наблюдается тенденция к выравниванию интегралов от сд, (7), а не только

предельных значений сд,. Следовательно, разность 72 - 73 = х)с0 - 71/2 убывает с ростом g (рис. 13), поскольку

I(7) = с, (7)g-14Ь, 71/2 = 7* -1-1 Г I(т) Ст =

* 0

= 7* - с, с, (т) Ст = ^, .

Рис. 12. Зависимость слагаемых в формуле (17) от v

Итак, при g << 1 0(о) ~ 0 , а при больших g следует брать 0(о) из £2 / 2D ~ [0 -12 / 3D][x -1] •

Идентификация модели III. Вырождение I ^ III происходит с ростом g при

b = b / g2 = const. Уравнения вида (14), (16) для модели III получаются формальным предельным переходом 1/Vb = 1/(gb1/2) ^ 0 :

to + 0 :

ІІ. х+1 2D х-1

,2

6D

[х+ 2] = [х- 1]t0 + t1/2 -х)«

(г = іЬ1/2В-1). Полагаем из ґ1/2 - Хс0 « 0, находим из второго уравнения приближение В(0) и значение Ь(0) из г . Решая краевую задачу численно, находим ^ и так до установления

(несколько итераций). Если рассматриваем переход III, то, совместив таким образом и

3Ш, растворимость g рациональнее подбирать по монотонности (см. рис. 6, 11).

В общей модели I выбор начального приближения - проблема, которую будем решать

аппроксимацией интеграла К = J с0 (ґ) Жґ в соотношениях (16), (17) (ґ(°) « ?).

Рис. 13. Зависимость 72 и 73 от g

Вместе с тем, 0 быстро растет (в том числе из-за роста 7*), так что при g >> 1 практически ,2 /(2В) и [х -1][0 - ,2 /(3В)] (установлено численно). Это дает хорошее начальное приближение 0 и 0(д) в форме линейной функции от £.

Аппроксимация интеграла К. Когда входная поверхность насыщается очень быстро, можно считать

К и ис0 = и^Ь-](Р -3) = и.,4Р -3 /(гВ).

Следующий шаг - учтем интегрально входной поток значением 3 :

до = -Ьдо + Р -3, до(0) = 0 э с0()) = = 490()) = сд1апЬ{ш72}, а = 2^Ь(Р-3).

Такое приближение функции с0 (7) качественно

отражает реальную картину: график является 8-образной кривой насыщения с характерными показателем экспоненты -а и временем установления т*. Поскольку Рд (7) — 3 , то стационарная концентрация не изменится (1апЬ — 1). Вычисляем интеграл от с0(7) (7*: ехр{-а7*} » 0):

К и сд [7* - [Ь(Р - 3)]-1/21п 2] =

= с>[)* - g-1[~(Р -3)]-121п2] '

Эта оценка неэффективна при I — III с ростом g , вырождаясь в К и с0)* , 7с0 и 0 (в то время как 7с0 растет, см. рис. 4, 5, 14).

При относительно малых g полагаем в (14) 0 и 0 на фоне большого 70 и принимаем

7с0 = 7* - с0-1К и 1п2/д/Ь(Р - 3) в (16). Тогда получаем систему линейных уравнений для оценки параметров £,п (а значит и В, Ь):

г х+1 а/Р г х+2 £^—+ п^ = 70, £-— + х-1 21 0 3

+ ^"2-^ (х+ 1п2 - 1) = (х- 1)7о + 71/2 .

Перейдем к поиску приближения К в более широком диапазоне возможных значений растворимости g . В уравнении насыщения входной поверхности представим один сомножитель д0 (7) полученной аппроксимацией

~о()) = g 1с01апЬ{7а/2}: д0(0) = 0,

(18)

~0() ) = ^ Ь 1(Р - 3) гапЬ{^Ь( Р - 3)}

к(т) Ст = |Рд(т) Ст = 0 = £X——+п^. (19)

д о х-1 2

Из физических соображений примем Рд(7) = 3 [1 + А ехр{-Д}]. При выборе показателя в ориентируемся на величину В(п/,)2, характеризующую скорость убывания входного диффузионного потока в аппроксимирующей линейной краевой задаче II. Выше принята верхняя оценка показателя скорости насыщения поверхности а = 2д/Ь(Р - 3) . Значение в уточним позже. Предэкспонента А фиксируется в силу интегрального соотношения (19) (да и интересует нас именно интеграл К ):

7* 7*

| Р0 Ст=3[7* + А/в] =| Рд Ст=3[0+7*] э А/в=0.

д о

Уточнение показателя в. Интегрируя линейное уравнение

дд =-дд1апЬ{7а/2}а/2 + (Р - 3) - 3Аехр{-в7}, получаем

дд(7) = д/ Ь-1 (Р - 3) 1апЬ{7а/2} -__________3а

ехр{7а /2}+ехр{-7а/2}

-х

X

ехр{(-в- а/ 2)7} -1 + ехр{(-в+а/2)7} -1

- в - а/2

К = | с0(7) С7 и с

о

^ р1 и7 -1

в+а/2

7* —

1п2

л/Ь( Р - 3)

+

(20)

+ 4 4А3

2 ■'о иг +1

л п Л

Си +------1

Г_____

= а Г=~в'

90 (7) = -д0 (7) 1апЬ{7а / 2}а /2 + Р - Р0 (7).

По построению, если заменить входной поток Р0 (7) на 3 , то функция

удовлетворяет уравнению (18). Ищем такую аппроксимацию Рд(7) — 3 , которая бы привела к

оценке К И К , пригодной и при больших g .

В силу условия (2) ориентируемся на Рд < Р .

Кроме того, требуем соблюдения материального баланса (14):

в [/-2][Г+ 2]

При интегрировании выбиралось и = ехр{-в7} , ехр{-в7*} и 0 . Рассмотрим показатель в = В(п/,)2, характерный для краевой задачи II (решение представимо рядом по ехр{-В(лк / ,)2}). Тогда грубо

д/Ь(Р-3) иТЬР, 7 = а/в и £/п.

Число 7 имеет важный физический смысл - оно характеризует соизмеримость диффузионного и поверхностного времен. Поскольку речь идет о приближении, то разумно остановиться на рациональных значениях (чтобы интеграл вычислялся в элементарных функциях), например, 7 = 1/3,1/2,1, 2,3.

Пусть у = 1/2 (поверхностные процессы медленнее, относительно низкие температуры). Тогда получаем оценку

7* —

1п2

л/ь( р - 3)

8п + 41п 2 - 7

Ь( Р - 3)

60

А =0= ,2

—I _= = 0 =-+

44^-3 В

7 0 = -

с0

1п2

л/Ь( Р - 3)

+ 0-

«£+1 - 7 л/Ь tо,

3 8п + 41п 2 - 7

g

Р-3

60

а

(• 7* (• 7* -~

|о Р(7) С7 =|д Р(7) С7 Э 0 = --

А^ г11-а Го 1-

Вычисляя как и выше соответствующий интеграл К, получаем

+ Ау г 11_и^Си ( = 2лjb(P-J = ув ).

а -Ч + иу ' у У ’

Р21п2 - 3п + 23Ай(^)

(Р - 3 )а

И(у) =

Г

(Р - 3 )а

-1 1 - иУ

[у-2][Г+ 2]

+и

и' у п

—Си +----------

г 2 2

Отсюда следует, что при малых g имеем прежнюю оценку 7с0 (0 — 0). С ростом g при вырождении I — III первое слагаемое убывает и получаем быстрый рост 7с0 ~ 0. При у = 1, 2 изменится лишь константа в форме дроби: [2п-21п2-1]/3, [п-3]/2. Правда, в случае у = 2 (поверхностные процессы быстрее, относительно высокие температуры) приходится раскрывать неопределенность дроби в (20). Это делается стандартно. Обозначим числитель через /(Г). Тогда /(2) = (,

/'(2) = [п- 3]/2, /(у) = /'(2)(у-2) + о(| 7-21). Сокращаем на у-2 и переходим к пределу у — 2 . Не останавливаясь на арифметических преобразованиях, отметим лишь, что уточненное представление запаздывания 7 0 по-прежнему линейно по искомым переменным £,п .

Если нужно более точно учесть начальный быстрый рост входного потока от нуля, то полагаем

Рд(7) и Рд(7) = 1апЬ{7а/2}3[1 + Аехр{-в}]. Первый множитель в масштабе времени насыщения поверхности растет от нуля до единицы (8-образно), появляется начальный всплеск. Величина А по-прежнему определяется материальным балансом:

21п2

Линейность запаздывания 7с0 по п, 0 (а значит и по п, £) сохраняется.

Заключение

Параметры поверхностных и диффузионных процессов определяются в одном эксперименте. Это представляется более корректным, чем использование «табличных данных» по каждому коэффициенту из различных источников. К тому же свойства поверхности существенно зависят от условий производства материала. Информация на входе «черного ящика» (о концентрации или потоке при х = 0 ) существенно улучшила бы обусловленность обратной задачи идентификации. Желательно, чтобы скорости диффузии и поверхностных процессов были соизмеримыми по порядкам ([~1/2л/Р ] = [ В / ,] = с т / 5). Следовательно, речь идет о не слишком высоких температурах: химические процессы активируются с ростом Т значительно быстрее диффузии, модель I вырождается в II. В прикладном плане модель I нацелена на исследование тонких пленок или стенок, когда концентрации малы и требуется оценить медленное (но на значительном интервале времени) накопление водорода или его изотопов с внешней стороны.

В алгоритме используются только интегральные операторы обработки информации 3(7) , что нацелено на помехоустойчивость оценивания. Рассмотренный диапазон порядков величин очень широк. В конкретном случае следует переходить к безразмерным переменным х1:

В = Вх1, Ь = Ьх2, 5 = 5х3, g = «х4 (чертой фиксируются характерные порядки) с учетом общей априорной информации о материале (можно ли считать коэффициент растворимости малым/большим; достаточно ли велика/мала температура, чтобы считать поверхностные процессы сравнимыми с диффузией).

Работа выполнена при поддержке РФФИ (грант 09-01-00439).

Литература

Водород в металлах / ред. Г. Алефельд,

B. Фёлькль. М.: Мир, 1981. Т. 1. 506 с., Т. 2. 430 с. Габис И. Е., Компаниец Т. Н., Курдюмов А. А.

Поверхностные процессы и проникновение водорода сквозь металлы // Взаимодействие водорода с металлами / ред. А. П. Захаров. М.: Наука, 1987.

C. 177-206.

Зельдович Я. Б., Мышкис А. Д. Элементы математической физики. М.: Наука, 1973. 352 с.

СВЕДЕНИЯ ОБ АВТОРАХ:

Заика Юрий Васильевич

зав. лаб. моделирования природно-технических систем, д. ф.-м. н.

Институт прикладных математических исследований КарНЦ РАН

ул. Пушкинская, 11, Петрозаводск, Республика Карелия,

Россия, 185910

эл. почта: 7а1ка@кгс.кагеНа.ги

тел.: (8142) 766312

Борматова Елена Павловна

доцент, к. т. н.

ГОУ ВПО Петрозаводский государственный университет, кафедра математического моделирования систем управления

ул. Ленина, 33, Петрозаводск, Республика Карелия,

Россия, 185910

эл. почта: Ьогша1@ре1геи.ги

тел.: (8142) 719606

Кунин Л. Л., Головин А. И., Суровой Ю. И., Хохрин В. М. Проблемы дегазации металлов. М.: Наука, 1972. 324 с.

Марчук Г. И. Сопряженные уравнения и анализ сложных систем. М.: Наука, 1992. 336 с.

Zaika Yu. V., Bormatova E. P. Algorithms of parameters estimation of hydrogen permeability model // NATO Science for Peace and Security, Series C, Carbon Nanomaterials in Clean Energy Hydrogen Systems, Springer, 2008. P. 403-414.

Zaika, Yury

Institute of Applied Mathematical Research, Karelian Research

Centre, Russian Academy of Science

11 Pushkinskaya St., 185910 Petrozavodsk, Karelia, Russia

e-mail: zaika@krc.karelia.ru

tel.: (8142) 766312

Bormatova, Elena

Petrozavodsk State University

33 Lenin St., 185910 Petrozavodsk, Karelia, Russia

e-mail: bormat@petrsu.ru

tel.: (8142) 719606