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

УДК 519.63

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

Оценка параметров модели термодесорбции водорода

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

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

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

Е. К. Kostikova, Ju. V. Zaika

Estimation of Parameters of the Hydrogen Thermodesorption Model

Degasification of the plate preliminary sated with hydrogen is considered. The experiment is done with the method of thermal desorption spectrometry (TDS). In a corresponding regional problem with nonlinear boundary conditions the basic physical and chemical processes are taken into account: diffusion and desorption. Is represented the technique of estimating the parameters of the model by results of measuring and results of numerical modeling.

Keywords: hydrogen penetrability, nonlinear regional problems, parametrical identification.

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

Водород рассматривается как один из перспективных экологически чистых энергоносителей. Кроме того, безопасность систем транспортировки и переработки углеводородного сырья во многом определяется уровнем защиты конструкционных материалов от водородной коррозии. Метод термо-десорбционной спектрометрии (ТДС) является одним из основных при исследовании взаимодействия водорода с твердым телом [2, 3, 5]. Пластина толщины I из металла или сплава, нагретая до температуры ' ~ ■' , находится в камере с газообразным водородом под давлением Р. После насыщения растворенным атомарным водородом образец быстро охлаждается (отключается ток нагрева), камера вакуумируется, и в условиях медленного нагрева с помощью масс-спектрометра определяется де-сорбционный поток. По этой информации судят о характеристиках взаимодействия водорода с исследуемым материалом.

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

W= t е c{Q' т) = = ~ т)г v е 10 4

¿XIX&.Q) = мг^еа ытхл«= t * га tл с 2 ■

Здесь С"'—"1--1, - концентрация атомарного водорода (Н), растворенного в пластине, fibCt}" 0*

Ctfil * ifr fif)= ^tft): f- - время дегазации; . ■-■ - коэффициенты диффузии и десорбции;

i*'- ^ ~ ^^ ^St^ - плотность десорбционного потока (торцами пластины пренебрегаем). Квадратич-ность десорбции связана с тем, что водород диффундирует в металле в атомарном состоянии, а поки-

дает поверхность (при

* = 0 Их = t

) в молекулярной форме. Коэффициенты диффузии и десорб-

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

выполняется закон Аррениуса:

ЛС? = ¿?р«ср

Д = contt Въ

энергии активации,

R = S.5L441

Дж/[моль-К] — универсальная газовая по

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

T(t) = Та 4 vi v > О

Сокращенно обозначаем

постоянная). Нагрев обычно линейный: £7(7№ Ь£Я - ШТГЙ

Что касается начальных данных '-т то в силу непродолжительности подготовительного этапа (охлаждение и вакуумирование) обычно считают начальное распределение равномерным:

= Г = СО/ЯП Здесь ■ = ^(Р"' ' ) - равновесная концентрация. Несогласованность начальных и граничных условий при этом непринципиальна, поскольку будем использовать лишь интегральные соотношения (решение задачи (1)-(2) понимается как обобщенное). Для тонких мембран следует учесть «начальный прогиб» концентрации по краям. Ограничимся параболической аппроксимацией

,. = *

= с- яа|> - 5 ^ - 2 5 л» > с-.

Цель работы состоит в составлении вычислительного алгоритма для определения по плотности потока термодесорбции ■ - КЬ ^ ( ДО * 0 : £ Т,) параметров аи, , -V, , характеризующих водородопроницаемость конструкционного материала.

Трудности решения обратных задач хорошо известны [9, 10]. В частности, разработаны градиентные алгоритмы минимизации в пространстве параметров среднеквадратичной невязки экспериментальных и модельных кривых [1]. Но на каждой итерации в общем случае приходится численно решать краевые задачи при текущих приближениях параметров. К тому же, как правило, сходимость лишь локальная. Учет специфики метода термодесорбции позволил разработать алгоритм идентификации, в котором основная вычислительная нагрузка связана с использованием квадратурных формул, а не решением краевых задач.

Для тестирования алгоритма решения обратной задачи сначала численно генерировались модельные кривые, порождающие параметры которых затем «забывались». Из-за большого разброса порядков величин при численном моделировании да проводилось масштабирование: Л — ■- , - —

£7и„ = ±bui i

£> Ъс

°=ё. т.

Лгг*1

Здесь за преобразованными параметрами модели оставляем прежние обозначения ^ , '■■ , ---о. Для определенности авторы ориентировались на данные по вольфраму, являющемуся одним из конструкционных материалов в реакторах [8]: ® = &-034 х Ш-*" 1/смд (. =13 Сги к), о = £00 К, ; =2 К/с, и = аео с { = 0,1 см Ьа = 6х 10"-= см4/с Ей = 3?.5£у Во

=37,629 кДж/моль.

кДж/моль, "D

_'0 = . . х .f = см2/с;

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

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

распределения ■ ■1. Поэтому целесообразно в качественном плане за первое приближение взять параболическую аппроксимацию

dt. * Cit. v> - ffft}- ^(t)i.v - Г^Р,

= f.

AW

Sfo>i

Считаем известной равновесную растворимость ^ = она определяется давлением

молекулярного водорода, температурой насыщения и пропорциональна корню из давления. Симметрия выполнена; функция 0 аппроксимирует срединную концентрацию 7' и0 . ^ ® . Поскольку к моменту окончания эксперимента произошла практически полная дегазация

образца (Гв"«* 1 ^ . ■ ). определим константу -^а в начальных данных ^ = — /(,]'

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

. Известная по итогам эксперимента величина - -

Аш =

Отсюда

равна половине коли-

чества десорбировавшегося водорода (в атомах), отнесенного к ■и" поверхности (■'■

= 0

или

л = г

не важно в силу симметрии). Кроме того, условия согласования ^ ¡Ц при ^ ® началь-

Д

ных данных и граничных условии дает зависимость

:

РШ^ = КоХГ- ЛХУ =* = ^р - 4.

Перейдем к конкретизации функций ^ ' в квадратичной аппроксимации '). Из условия

материального баланса выразим ^^ ^ и подставим в

( ол - 5(0 = Ге&х)1£х. ■ \jmdT

= Щ^-лт- Л: + <?■::.■ 2| Дг)А.

Чтобы найти оставшийся функциональный параметр подставим полученное выражение для в граничное условие ^^ГаГ^К эт0 соотношение позволяет выразить функцию

Л*Т.* О через коэффициенты

модели

D h

и известную по экспериментальным данным информа-

цию Я^. Оба корня квадратного уравнения относительно ^^ положительные, выбираем меньший из них (по физическому смыслу — 0 ) . Из условия У: = - ÍJ = ^ЛО^) 2 ) получаем

соотношение для оценки ^и, , ^в, :

,7ft) _зд<П

г- .■""TV а', - -г -

= Ott-гьц Г 4 3D(T)

^ью тт)

2 = ,1 ''.¡7 6 ГО г . .]■ Поскольку плотность десорбции соответствует исходной модели (1)-(2),

а на предварительном этапе оценки - , --' используется параболическое приближение концентрации Н в объеме, то это равенство является приближенным.

График I ^ имеет характер всплеска с последующим затуханием, причем на начальном и конечном этапах измерения менее точны. Поэтому ограничимся ^ е - нормируем уравнение

на

-/kwäi: fiftJ" ÍA3J

^ \ ' / И ]

выделим безразмерные переменные:

í(ftfíauT(-i}-Cv(i+ m)>/рш Wty-iy.tefaimwmw* V№

Формально допуская значения ^ '"■ ® , удобно считать новые переменные -1'*■''" ^ . ¡ «аррениусовскими»:

An

Er = Вь - &D-

■■ ЗДф ' * = / — '

Тогда, обозначая Ч ■ ■ V1)

, получаем уравнение

Дь^^Л} = (-1) - чфХ) - = о,

Заметим, что параметры --^в, £ входят в соотношение лишь неявно посредством ДЙ.

V т~ = /у 4

Преобразуем величину ' с учетом связи 5'? (см. (4)):

te [^1.^21, (?)

Величина 'шп!!! зависит от всех входных данных * - з}. Запись ХУ = £¿0 ) показывает, что значения ", -^о уже найдены, а .■'" ■ при решении обратной задачи воспринимается как заданная фиксированная функция времени. Аналогично представим ■ '■ •

Х = Г9шхр |

ПИ. ЛП01Ш1 ИЧШ 1 авпм

Ч" ШУ

Ф - Ла'бГ

Подставляя выражения ■* , ' в уравнение (7), получаем зависимость .* = ^Далее с

учетом зашумленности реальных измерений и погрешности параболической аппроксимации целесообразно следовать методу наименьших квадратов (МНК):

Р<Х.г£"Л-ГЕГУ ■ I гсЬ.

'Л

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

Перейдем к изложению результатов численного моделирования. Разностная схема решения задач термодесорбции изложена в [4], здесь на этом не останавливаемся. График плотности потока водорода для параметров, указанных выше, представлен на рис. 1, 2.

ю

500

700

900

1100

1300

ЯТ)х 10"

\Т=1300

/ т = 1200 \

/у 1 - 1000 \\

100

200

300

400

500

Рис. 1. Плотность термодесорбции. Влияние температуры насыщения

300 500 700 900 1100 1300

Рис. 2. ТДС-спектр. Влияние скорости нагрева

Для оценки значений предэкспонент и энергий активации диффузии и десорбции £ ii.ii А £ использовались метод наименьших квадратов (МНК) и метод моментов (ММ) применительно к уравнению (7) ( / = 0 ), в которое подставлены выражения ' согласно форму-

£>*Г} = Da «го ■

. ГЮ= Т„ 4*

и выражение

лам (6), " =

Дф — Вш т-с-г из соотношения (4). На рис. 3, 4 показано, что задача Ья * г-~и1- () хорошо

обусловлена по каждому из коэффициентов ^ '■" (один из них фиксировался равным «истинному»). При этом дополнительное соотношение (4) не учитывалось при построении поверхности на рис. 4, но для рис. 3 оно оказалось необходимым, иначе отсутствует экстремум в физически оправданном диапазоне.

Более подробно остановимся на уравнении = 0 полученном после подстановки в

(7) выражений : ^о в соответствии с формулами (4), (6). «Истинные» значения параметров:

= £и= 1.929 ¿Ту = 17.8+9.

На рис. 5, 7, 9 варьируем лишь один из параметров

-а- Решение задачи И И-: по ^я (см. рис.5) дает относительную погрешность

. Это приводит в соответствии с формулой (8) к ^^о' = ^-77 и

бШо) = 0.17

в силу (4). Постоянство & у означает постоянство

Обратимся к методу моментов. Обозначим

М

, = М^^Ец йу)= | %<t)f(t;ZQiEXtEY)<lt.

чда = -. - ^ &}=

ИГ1«.

t е [ti,t:l,it = SO

Ограничимся функциями Г- " И - ,

= -т£0 (1^*= ЭСО) Эксперименты показали, что использование метода моментов дает точность в среднем такую же, как и метод наименьших квадратов (см. рис. 6, 8, 10). Результаты применения параболического приближения позволили решить обратную задачу для исходной распределенной модели с погрешностями, указанными в таблице 1. Подчеркнем, что параболическое приближение является грубым для краевой задачи (1)-(2). Его цель - «попасть в порядки» оцениваемых коэффициентов -■' , . Значения предэкспоненты '^о определяется заметно хуже, что объясняется его малым абсолютным значением (коэффициент при квадрате концентрации). Можно было бы перейти к двойной точности вычислений на всех этапах моделирования, но это излишне при зашумленных входных данных обратной задачи.

Таблица 1

Параметр Исходные данные Полученные значения Относительная погрешность

К е* io-ls 1,514 МО"11 152.3»

Еъ # 39.&S9 45.100 lflJBfi

D. 4ix 10"я 2,330 х Iß"3

37,629 36.545 7 1%

Hf.

Ä'fZJ - 52 %

Й(Ь„) = 77%

<S(D0)-I7%

г0 ист.

52%

Siь„)- 77%

17%

7,„ ист.

0 0.04 ^чДОВ Ü12 0.16 0.2

м\

\ z X^MJ

\ *

Рис. 5.МНК. Поиск :0

Рис. 6. ММ. Поиск - о■

l|f||,

— 138%

ä(Ki;)=7%

Ех ис!.

Ех нет.

0 2 Гр К 1(1

// / f /

y7 s; 142 %

А /м, 13%

/ V <i{E0H7%

Рис. 7. МНК. Поиск аХ

Рис. 8. ММ. Поиск -X

imi,

¿(Е[,)=25%

Гуист.

20 21 Рис. 9. МНК. Поиск Zy

6( Еу) S 29 % i(Eb) S 26% <S(E„) i 2K %

Рис. 10. MM. Поиск Еу

Начальные приближения энергий активации -'V. £ в диапазоне нескольких десятков кДж/моль можно указать для конкретного материала из физико-химических соображений. Приближение о^

берем в силу * I 01* . Именно как приближение, в ТДС-

эксперименте /" '■'' известно с большой погрешностью.

Задача несколько усложняется, когда равновесную концентрацию с также приходится считать не-

Ул Ev

Y значения с,

известной. Тогда задача четырехмерная в соответствии с (7). По оценкам ■'■ а, ■■>о находятся из уравнений (3), (4). Далее переходим к локальному уточнению оценок ¿-'о.

в соответствии с исходной моделью (1)-(2), при этом учитываем, что искомых независимых пе-

т- " - ■

ременных 3 в силу .

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

Для уточнения оценок ^ , ^ необходимо иметь в распоряжении семейство уравнений, связывающих параметры с экспериментальной информацией. Укажем один из возможных подходов, основанный на общей идее использования сопряженных уравнений [7].

Интегрированием по частям для достаточно гладкой функции V'"1 - г Л'} получим:

♦ I Bttüim^it- о sft - \ пик/т^ЮФхЪи} ctt-

iv-

- С ^ А ) «я + г <£ж.

-в -о

Здесь опущен двойной интеграл, поскольку в дальнейшем изложении считаем функцию " ' ^

££

подчиненной сопряженному уравнению д-£ да1 . Кроме того, пренебрегаем интегралом от V &и ■ ' по ' с учетом С" *" ® . Косвенно ограничиваемся не слишком быстро растущи-

по

ми по 1 функциями V■»■ ■ V Подчеркнем, что краевые условия не ставятся, «пробных» функций ^ бесконечно много. Простые варианты = Л приводят к уравнению материального баланса, которое уже использовалось для оценки константы Выберем, например, ^йкрол При

нормировке * получаем

л

ijbtC, jt> = ехр( а" у* t„r t <7* f, 1. тj ■ [ Di i■ rfs.

А

Перепишем соотношение (10) в обозначениях

(?f - 1

Л" = i Jß dt,

Y = I D-ifb'1 ß dt,

^ -1

Получаем семейство уравнений. Параметр ^ целесообразно варьировать в пределах ^! " 1 .В таблице 2 приведены значения параметров, полученные решением системы уравнений (11) для

'-■■'- 1 ■ 1 . Энергии активации восстанавливаются с большой точностью (их влияние на кинетику дегазации очень велико). Предэкспонента я определяется хуже в силу ее малого абсолютного значения.

Таблица 2

Параметр Исходные данные Полученные значения Относительная погрешность

Ьа 6x10-« S.ißgx Ю-1* G Г*

fö SftSSfl 39.559 096

Do 41 х IQ"* 4104* 10"3 l.^SS

Еэ 3?<S29 036

Рис. 11 Экстремум Gl v

Рис. 12. Экстремум

На рисунках 11,12 представлены поверхности £ = 'Ч-11" 5 при фиксированном П — О шп а ■ з соответственно. Из рис. 11 видно, что для коэффициента '■■ важно найти хорошее начальное приближение. На рисунках 13—16 кривым 1, 2, 3 соответствуют ^ = * 9 ¿А . указано максимальное значение относительной погрешности при выбранных & . Один параметр варьируем, остальные подставляем «истинные». Видно, что уравнение = ® хорошо обусловлено по каждому из параметров модели. Более «удачным» является ^ = 9 .

F х 10

Fx 10

10

-10

5(Ъо) = оз %

1

5.975е-12 <Ш5?Т?-

Ь0 ист.

^чЛ

\з

10

5(ЕЬ) = 0.04%

39.5

39>Г ист .у

39.6

Рис. 13. Вариация

F х 10

Рис. 15. Вариация ^'о

Рис. 14. Вариация —Ь

Fx 10

10

-10

<S(ed) = o.ooi %

1___—/

■ЗТ6285~ ^^ЪТ.ЫЧ / 37.6295

F., ист./

3/

Рис. 16. Вариация -Г1

Таким образом, последовательность представленных этапов параметрической идентификации позволяет восстановить параметры модели с относительной погрешностью, которая с запасом «поглощается» точностью ТДС-эксперимента. Входные данные А^ в уравнении стоят под знаком интеграла, что обеспечивает определенную помехоустойчивость. Если мембрану нельзя считать тонкой (это зависит от материала и условий эксперимента), следует использовать вариант

(в середине пластины равновесная концентрация, и только у самого края концентрация немного меньше).

Библиографический список

1. Алифанов О. М., Артюхин Е. А., Румянцев С. В. Экстремальные методы решения некорректных задач [Текст] / О. М. Алифанов, Е. А. Артюхин, С. В. Румянцев. - М. : Наука, 1988. - 288 с.

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

3. Водород в металлах [Текст] : [пер. с англ.]. В 2 т. / под ред. Г. Алефельда, В. Фёлькля. - М. : Мир, 1981.

4. Заика, Ю. В., Костикова, Е. К. Разностная схема для краевой задачи ТДС-дегазации с динамическими граничными условиями [Текст] // Ученые записки ПетрГУ. Серия «Естественные и технические науки». - № 7 (101). -2009. - С. 65-70.

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

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

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

8. Писарев, А. А., Цветков, И. В., Маренков, Е. Д., Ярко, С. С. Проницаемость водорода через металлы [Текст] / А. А. Писарев, И. В. Цветков, Е. Д. Маренков, С. С. Ярко. - М. : МИФИ, 2008. - 144 с.

9. Самарский, А. А., Вабищевич, П. Н. Численные методы решения обратных задач математической физики [Текст] / А. А. Самарский, П. Н. Ващибевич. - М. : Едитоиал, УРСС, 2004. - 480 с.