Оценка параметров диффузии и десорбции водорода в краевой задаче ТДС-дегазации Текст научной статьи по специальности «Математика»

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

УДК 519.6, : 539.2

ОЦЕНКА ПАРАМЕТРОВ ДИФФУЗИИ И ДЕСОРБЦИИ ВОДОРОДА В КРАЕВОЙ ЗАДАЧЕ ТДС-ДЕГАЗАЦИИ

Ю. В. Заика, Е. К. Костикова

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

Рассматривается дегазация пластины, предварительно насыщенной водородом. Эксперимент проводится методом термодесорбционной спектрометрии (ТДС). В соответствующей краевой задаче с нелинейными граничными условиями учтены основные физико-химические процессы: диффузия и десорбция. Предложена методика оценки параметров модели по результатам измерений.

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

Yu. V. Zaika, E. K. Kostikova. ESTIMATION OF HYDROGEN DIFFUSION AND DESORPTION PARAMETERS IN THE BOUNDARY-VALUE PROBLEM OF TDS-DEGASSING

Degassing of a plate saturated with hydrogen is considered. The experiment is based on the thermal desorption spectrometry (TDS) method. The model is a boundary-value problem with nonlinear boundary conditions. The main physical and chemical processes, such as diffusion and desorption, are taken into account. The method of estimating model parameters by measurements is proposed. The work was supported by the Russian Foundation for BasicResearch (grant 09-01-00439).

Key words: hydrogen permeability, nonlinear boundary-value problems,

parameters identification.

Постановка задачи

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

онной спектрометрии (ТДС) является одним из основных при исследовании взаимодействия водорода с твердым телом [Кунин и др., 1972; Водород в металлах, 1981; Взаимодействие водорода.., 1987]. Пластина толщины I из металла или сплава, нагретая до температуры Т = Т, находится в камере с газообразным водородом под давлением р . После насыщения растворенным

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

Рассмотрим симметричную по постановке эксперимента нелинейную краевую задачу ТДС-дегазации:

dc = D(T)^ tе(°,t.), Xе (0,£), (1)

dt дх2

c(°, х) = p(х) = <p(£ - х), х е [0, £], (2)

D(T )cx (t,0) = b(T )c°2(t),

D(T)cx(t,£) = -b(T)c£2(t), tе [0,t,]. (3)

Здесь c(t, x) - концентрация атомарного водорода (H), растворенного в пластине, c0(t) s c(t, 0), ct (t) s c(t, £), c0(t) = c£(t);

t* - время дегазации; D , b - коэффициенты диффузии и десорбции; J(t) = b(T)c)2£(t) -

плотность десорбционного потока (торцами пластины пренебрегаем). Параметры модели D и b зависят от температуры T. Как правило, в «рабочем диапазоне» достаточно точно выполняется закон Аррениуса:

D(T) = D0 expj-^D / [ RT ]},

b(T) = b0exp[- Ebl Mb

D0, Ed , b0, Eb, R = const ( Ed , Eb - энергии активации). Нагрев обычно линейный: T(t) = T0 + vt, v > 0. Сокращенно будем обозначать D(t) s D(T(t)), b(t) s b(T(t)) .

Что касается начальных данных p(х), то в силу относительной скоротечности подготовительного этапа (охлаждение - вакуумиро-вание - нагрев) можно считать начальное распределение практически равномерным: p( х) = c = const. Несогласованность граничных условий при этом непринципиальна, поскольку используем лишь интегральные соотношения (решение (1)-(3) понимается как обобщенное). Для тонких мембран следует учесть «начальный прогиб» концентрации по краям. В статье ограничимся вариантом

p(х) = c - A°(х - £°)2k, k = 1, £0 = £/2, A° > 0.

Цель работы состоит в разработке вычислительного алгоритма для определения по плотности десорбции (потоку) J(t), tе[0,t.] (J(t) ~ 0,

t > t*), параметров D0, ED, Ь0, Eb, характеризующих водородопроницаемость конструкционного материала.

Из-за большого разброса порядков величин при численном моделировании J (I) (в соответствии с краевой задачей (1)-(3) и теорией разностных схем [Самарский, Гулин, 1989]) проводилось масштабирование:

Duz

■■£z, z = [0,1],

= ±bu01 s ±J,

u = c/c D = D/£2,

ut = Duzz, b =bc/£,

(0, z) = 1 - A(z - 0.5)2, A = A£2/c .

Параболическое приближение

Сходимость в нелинейных обратных задачах параметрической идентификации, как правило, локальная. «Куполообразный» характер распределения с(^ х) известен. В эксперименте

* < мм. Поэтому целесообразно в качественном плане за первое приближение взять параболу с^, х) ~ ~^, х) = B^) - A^)(х - *0 )2,

2*0 = * , А(0) = А0, В(0) = с .

Считаем известной равновесную растворимость с = С( р, Т) ~ у[р, она определяется давлением и температурой насыщения и пропорциональна корню из давления. Симметрия выполнена, функция В(^ > 0 аппроксимирует с(^*0), А(t) > 0 , t > 0 . Считая, что к моменту t» произошла дегазация образца (с(^ х) ~ 0, t > t*), определим константу А0. В силу симметрии и материального баланса

S, s j J(t)dT = j °[c-A0(x-£°)2]dr =

= c£°- 4)£°/3.

(4)

Отсюда А0 = 3(с*0 - )/£\. Кроме того, ус-

ловия согласования начальных данных и граничных условий дает D0 / Ь0 = /(Еь - ED):

Д0)А/ = Ь(0)[с - АД]2 ^

Ь0 (5)

= ехр{( ED - Еь )/№]} - А^]2.

Теперь воспользуемся квадратичной аппроксимацией. Из условия материального баланса выразим В(^ > 0 и подставим в ~(^ х) :

Р * 0 Р * 0 ~

[ <р(х)dх - 5(0 = [ с(^х)dх,

J 0 J 0

5(t) = [ J(т) dт ^

0

х

°,1

и

с*0 - А-* /3 - 8(?) = В^)*0 - A(t)*31 /3,

B(t) = A(t)*20 /3 + с - А0*0/3 -5(?)/*0,

^ ~(?, х) = Q(t)/* - А(?)[х2 - *х + *2/б],

р?*

Q(t) = 2| J(т) dт.

Подставим теперь выражение для аппроксимации с(?, х) в граничное условие

D(T )~х (?,0) = Ь(Т )~02(?), выразим функцию

А (?) > 0 через параметры модели D, Ь и известную по экспериментальным данным Q(t) . Оба корня квадратного уравнения относительно А(?) положительные, выбираем меньший из

них (по физическому смыслу с0(?) > 0). Из условия ^[7 = ~0 л/Ь получаем соотношение для оценки D0, ЕГ), Ь0, Еь:

J (t) 3D (T )

1 + Q(,)Ш± - 1

3 D (T)

, (б)

I

(по порядку) и выделим безраз-

мерные переменные:

I (t) I™ = (V1 + 2Q (t К1 J max x -1) Y,

I(t) = J), x =

t* J mxb(T)

Y = -

3D(T)

3D(T ) ’ I max ^Vb(T) ’

t є [t1, t2] с (0, t*) . Формально допуская E < 0 , удобно считать новые переменные x(t) = x(T(t)) , Y (t) = Y (T(t)) «аррениусовски-ми»:

3Do

x 0

t* J max b0

3D0

Y0

Ex = Eb - Ed , Ey = Ed - Eb / 2 Тогда

f (t; x 0, Ex ,Y0, Ey ) = I (t) I max --(1 + q(t)x -1) Y = 0, q = 2Qt*-1-J^.

(7)

Заметим, что параметры начального распределения Д, С входят в соотношение лишь неявно посредством J (і) .

Преобразуем теперь 7 с учетом связи

А)/ *0 = / (Ех ) (см. (5)):

7 = 7) ехр{- { /[ЯТ(і)]}=

= ^ ехр{- Ех /[Д70]}ехр{-Е7 /[ЯТ(і)]},

_3 АоС 2 і 2оУ*7

Z 0 =

, Z0 = Zo(bo) ^ bo.

41

max

Разумеется, /max зависит от всех входных данных {p, D, b}. Запись Z° = Z°(b°) означает, что значения c, A° уже найдены, а J(t) при решении обратной задачи воспринимается как заданная фиксированная функция времени. Аналогично представим X:

X = X°exp{- Ex /[RT(t)]} =

Ь (Т) * Ь (Т)

Т = Т(?) , ? е [0, ?*

Поскольку 7(?) соответствует исходной модели (1)—(3), а на предварительном этапе используется параболическое приближение концентрации Н в объеме, то это равенство является приближенным.

График 7 (?) растет, потом убывает, причем на начальном и конечном этапах измерения менее точны. Поэтому нормируем уравнение на

[c - A/o]2

^expEx /[R70]}exp{-Ex /[RT(t)]}.

Подставляя в уравнение (7), получаем зависимость / = /(?; 20,Ех,Ег ) . Далее с учетом зашумленности реальных измерений 7(?) и погрешности параболической аппроксимации целесообразно следовать методу наименьших квадратов (МНК):

Е(2й,Ех,Ег) = [2 /2(т)dт^ min, [t1,?2] ^ (0,?*)-•’<1

Производные функции Е по указанным аргументам можно выписать явно (подсчет интеграла по квадратурной формуле считаем элементарной операцией). Задача оптимизации решается стандартными средствами (авторы использовали пакет 8сПаЬ) в физически оправданном диапазоне параметров (в пределах одного-двух порядков). Укажем принятые в данной работе опорные значения:

с = 1018 1/см3; Т0 = 300 К; * = 10 “2 см;

Ь0 = 10-3 см4/с; Еь = 90 кДж/моль;

D0 = 10-2 см2/с; ED = 40 кДж/моль.

Результаты вычислительных экспериментов показали, что без учета погрешностей измерений и модели (т. е. 7 (?) определялся решением

прямой задачи (1)-(3)) погрешность оценивания параметров находится в пределах 20 %. Начальные приближения энергий активации ED, Еь в диапазоне нескольких десятков кДж/моль можно

указать для конкретного материала из физикохимических соображений. Приближение Ь0

(Ъ0) берем в силу

7 (0) = Ь(Т0 )с02,*(0) = Ь0 ехРЬЕь / Щ}\° - А0*20]2.

Именно как начальное приближение в реальном ТДС-эксперименте значение 7(0) известно с большой погрешностью.

Замечание. Задача несколько усложняется, когда равновесную концентрацию с также приходится считать неизвестной. Тогда задача четырехмерная в соответствии с (7). По оценкам Х0, Ех, У0, Ег значения с, А0 находятся из

уравнений (4), (5).

Далее переходим к локальному уточнению оценок D0, Е1), Ь0, Еь в соответствии с исходной

моделью (1)-(3). При этом учитываем, что искомых независимых переменных 3 в силу DJЬ = I(Еь -ED).

Функция Грина

Поскольку зависимость от времени 7(?) известна по результатам ТДС-эксперимента, решение краевой задачи удобно представить с помощью функции Грина [Мартинсон, Малов, 2002, гл. 2]. Для (р(х) = с имеем

х

c(t,x) = c - jD 1(t)J(t):

o

х {G( x, t,0, t) + G (x, t, і, т)},

nnx nny

G( х, t, y,T) =

1 2 I n П / \ I плх i

= — + — > exp^—-— (T-1)!>cos---------cos-

£ £ [ £2 J £ £

Для уточнения оценок параметров модели используем соотношение J = bc° £ , поэтому далее нас интересует представление граничной концентрации c°(t). Для варианта р(х) = c - A(х - £°) получаем:

2 1

c°(t) = c -~ jJ(t)dT-

- 4і

4 ^ і i exp{-inY(t, t)}J (t) dr,

і n =1 0

I

Y(t, T) = j D( s) ds, ln = (п/ і o )2.

Заменим в скобке [...] экспоненту на [exp[„Y(t,0)}-1]+1:

2 1 ю 4 ю

c°(t) = c - - j J (T) dT + 4 A° (t) - - > Jn (t),

£ A n=1 £ n=1

(t) =

1 - exp{-i„Y(t,0)},

ln ’

1

7п (?) = | ехр{-^„г(?,т)}7 (т)dт.

0

Соотношение -\[Г = 4Ъс0(?) имеет форму семейства уравнений для оценки параметров:

Ф = Ф(?; Do,ED, Ь0,Еь) = 0,

? е [?1, ?2] с (0, ?*). При численной реализации ряды заменялись частичными суммами:

Ф = Л-л(Ь X (8)

2 ‘ N1 4 N2

p(0) - -j j J dT + 4 Ao УГ Vn (t) - 7I Jn (t)

Vn (t) - і I n>

I

= 0.

Используя пакет Scilab, численно решалась за-

дача

j^2 Ф 2dT ^ min при N1 = N2 = 5 .

За начальное приближение принимались значения, полученные в предыдущем пункте. Уровень ошибок оценивания в среднем понизился до нескольких процентов (на «идеальном» сигнале 7 (?) из (1)-(3)). Численные эксперименты показали, что при Ni > 10 невязка резко возрастает, возможно, из-за накопления вычислительных погрешностей.

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

Если неопределенность значений D, Ь очень велика (например, 3-5 порядков), то нелишне иметь в распоряжении дополнительное семейство уравнений, поскольку на универсальный алгоритм параметрической идентификации нелинейной распределенной модели рассчитывать не приходится. Укажем один из возможных подходов, основанный на общей идее использования сопряженных уравнений [Марчук, 1992].

Используя интегрирование по частям, для произвольной достаточно гладкой функции у/(<, х) получим:

<* *

0 = Ц^(?, х)Х - Dcхх ]dхdt =

0 0

<* <*

= 17 (? )^(?,0) dt +17 (? )^(?, *) dt +

Г

/*

- jD(t)д/J(t)b-1(t)щх (t,£ ) dt -

-1D(t^7(?)Ь (?)/х(?,0) ^ -

0

£ £

- с| у/(0, х) dх + а| (х - *0 )2/(0, х) dх . (9)

0 0

Здесь опущен двойной интеграл, поскольку в дальнейшем изложении считаем функцию /((?, х) подчиненной сопряженному уравнению

Э //Э? = —D д 2//дх2. Подчеркнем, что краевые условия не ставятся, «пробных» функций / бесконечно много. Выберем, например, /(?, х) = в(?) ехр ох . При нормировке в(?*) = 1 получаем

/(?, х) = ехр{о2^(?*, ?)}ехр ох,

г(?,т) = 10^ (*) ^ .

Перепишем (9) в обозначениях X = |<* Jpdt, У = |<* D^^Jb—1вdt:

г=ехР°+1, ох- + о'-у + ^01 X

ехр о* -1 о

х{2А, - о2 (х - А0[*20 -*г/о])}= 0. Аналогично рассматривается вариант

/(?,х) = в(?)этох (Д(()со8ох/: /(<, х) = ехр{-ст2у(/*, /)}ш ох

(/(<, х) = ехр{-о2^(<», <)}с08ох/,

каХ + а2Г +

Ж°).

а

х{2A° + а2(c -(°[£° -£к/а])}= 0,

Ks-

sina£ cos а-1

Ks-

cosa+1

- sm<

Получаем дополнительное семейство уравнений. Параметр о целесообразно варьировать в пределах о* ~1.

Заключительные замечания

В предлагаемой итерационной процедуре оценивания используется подсчет интегралов по

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

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

зав. лаб. моделирования природно-технических систем, д. ф.-м. н. Институт прикладных математических исследований КарНЦ РАН

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

Россия, 185910

эл. почта: zaika@krc.karelia.ru

тел.: (8142) 766312

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

(р(х) = с - А0(х - *0) , к > 1 (в середине пластины равновесная концентрация и только у самого края концентрация немного меньше).

При обработке зашумленных экспериментальных кривых следует, по-видимому, отказаться от использования условия согласования краевых условий (5). В интегральном смысле (рассматриваем решение краевой задачи (1)-(3) как обобщенное) считаем ф(х) = с (А0 = 0). В силу (4) имеем оценку равновесной концентрации с . Далее считаем параметры X 0, Ех ,У0, ЕУ

(однозначно определяющие D0, Ев, Ь0, Еь) независимыми в соотношении (7). Последующие преобразования X, У не проводим, уравнение (5) и соответствующее соотношение D0/Ь0 = I(Еь - Ев) исключаем из модели. Задача оптимизации становится четырехмерной.

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

Литература

Взаимодействие водорода с металлами / ред. А. П. Захаров. М.: Наука, 1987. С. 177-206.

Водород в металлах / ред. Г. Алефельд, В. Фёлькль. М.: Мир, 1981. Т. 1. 506 с. Т. 2. 430 с.

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

Мартинсон Л. К., Малов Ю. И. Дифференциальные уравнения математической физики. М.: МГТУ им. Н. Э. Баумана, 2002. 368 с.

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

Самарский А. А., Гулин А. В. Численные методы. М.: Наука, 1989. 432 с.

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

Костикова Екатерина Константиновна

стажер-исследователь

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

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

эл. почта: fedorova@krc.karelia.ru тел.: (8142) 766312

Kostikova, Ekaterina

Institute of Applied Mathematical Research, Karelian Research

Centre, Russian Academy of Science

11 Pushkinskaya St., 185910 Petrozavodsk, Karelia, Russia

e-mail: fedorova@krc.karelia.ru

tel.: (8142) 766312