close

Вход

Забыли?

вход по аккаунту

?

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

код для вставкиСкачать
Труды Карельского научного центра РАН
№ 3. 2010. С. 45–50
УДК 519.6, : 539.2
ОЦЕНКА ПАРАМЕТРОВ ДИФФУЗИИ И ДЕСОРБЦИИ
ВОДОРОДА В КРАЕВОЙ ЗАДАЧЕ ТДС-ДЕГАЗАЦИИ
Ю. В. Заика, Е. К. Костикова
Институт прикладных математических исследований
Карельского научного центра РАН
Рассматривается дегазация пластины, предварительно насыщенной водородом.
Эксперимент проводится методом термодесорбционной спектрометрии (ТДС). В
соответствующей краевой задаче с нелинейными граничными условиями учтены
основные физико-химические процессы: диффузия и десорбция. Предложена
методика оценки параметров модели по результатам измерений.
К л ю ч е в ы е с л о в а : водородопроницаемость, нелинейные краевые задачи,
параметрическая идентификация.
Yu. V. Zaika, E. K. Kostikova. ESTIMATION OF HYDROGEN
DIFFUSION AND DESORPTION PARAMETERS IN THE BOUNDARYVALUE 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).
K e y w o r d s : hydrogen
parameters identification.
permeability,
nonlinear
boundary-value
problems,
___________________________________________________________________
Постановка задачи
Водород рассматривается как один из перспективных экологически чистых энергоносителей. Кроме того, безопасность систем транспортировки и переработки углеводородного сырья
во многом определяется уровнем защиты конструкционных материалов от водородной коррозии. Экспериментальный метод термодесорбци-
онной спектрометрии (ТДС) является одним из
основных при исследовании взаимодействия водорода с твердым телом [Кунин и др., 1972; Водород в металлах, 1981; Взаимодействие водорода.., 1987]. Пластина толщины  из металла
или сплава, нагретая до температуры T = T , находится в камере с газообразным водородом под
давлением p . После насыщения растворенным
атомарным водородом образец быстро охлаждается (отключается ток нагрева), камера вакуумируется, и в условиях медленного нагрева с
помощью масс-спектрометра определяется десорбционный поток. По этой экспериментальной информации судят о характеристиках взаимодействия водорода с исследуемым материалом.
Рассмотрим симметричную по постановке
эксперимента нелинейную краевую задачу ТДСдегазации:
2
∂ c
∂c
= D (T ) 2 , t ∈ ( 0, t* ) , x ∈ ( 0,  ) ,
∂t
∂x
c (0, x ) = ϕ ( x ) = ϕ (  − x ) , x ∈ [ 0,  ],
(1)
(2)
D (T )c x (t ,0) = b(T )c02 (t ) ,
D (T ) c x ( t ,  ) = − b(T ) c2 ( t ) , t ∈ [0, t* ] .
(3)
Здесь c(t , x) − концентрация атомарного водорода (H), растворенного в пластине,
c0 (t ) ≡ c(t ,0) , c (t ) ≡ c(t , ) , c0 ( t ) = c  ( t ) ;
t* − время дегазации; D , b − коэффициенты
диффузии и десорбции;
J (t ) = b(T )c02, (t )
−
плотность десорбционного потока (торцами
пластины пренебрегаем). Параметры модели D
и b зависят от температуры T . Как правило, в
«рабочем диапазоне» достаточно точно выполняется закон Аррениуса:
D(T ) = D0 exp{−ED / [ RT ]} ,
b(T ) = b0 exp{− Eb [RT ]} ,
D0 , E D , b0 , Eb , R = const ( E D , Eb − энергии
активации).
Нагрев
обычно
линейный:
,
.
Сокращенно
будем
обознаv
>
0
T (t ) = T0 + v t
чать D(t ) ≡ D(T (t )) , b(t ) ≡ b(T (t )) .
Что касается начальных данных ϕ (x ) , то в
силу относительной скоротечности подготовительного этапа (охлаждение – вакуумирование – нагрев) можно считать начальное
распределение
практически
равномерным:
ϕ ( x ) = c = const . Несогласованность граничных
условий при этом непринципиальна, поскольку
используем лишь интегральные соотношения
(решение (1)–(3) понимается как обобщенное).
Для тонких мембран следует учесть «начальный
прогиб» концентрации по краям. В статье ограничимся вариантом
ϕ ( x) = c − A0 ( x −  0 ) 2k , k = 1, 0 =  / 2 , A0 > 0 .
Цель работы состоит в разработке вычислительного алгоритма для определения по плотности десорбции (потоку) J (t ) , t ∈ [0, t* ] ( J (t ) ≈ 0 ,
t ≥ t* ), параметров
D0 , E D , b0 , Eb , характеризую-
щих водородопроницаемость конструкционного
материала.
Из-за большого разброса порядков величин
при численном моделировании J (t ) (в соответствии с краевой задачей (1)−(3) и теорией разностных схем [Самарский, Гулин, 1989]) проводилось масштабирование:
x = z , z = [0,1],
ut = Du zz ,
u=c c,
Du z
0,1
D = D 2 ,
= ± bu02,1 ≡ ± J ,
b = bc  ,
2
u(0, z ) = 1 − A(z − 0.5) , A = A 2 c .
Параболическое приближение
Сходимость в нелинейных обратных задачах
параметрической идентификации, как правило,
локальная. «Куполообразный» характер распределения c(t , x) известен. В эксперименте
 < мм . Поэтому целесообразно в качественном
плане за первое приближение взять параболу
c ( t , x ) ≈ c~ ( t , x ) = B ( t ) − A ( t )( x −  0 ) 2 ,
2  0 =  , A ( 0 ) = A0 , B ( 0 ) = c .
Считаем известной равновесную растворимость c = c ( p , T ) ~ p , она определяется давлением и температурой насыщения и пропорциональна корню из давления. Симметрия выполнена, функция B (t ) > 0 аппроксимирует
c(t ,  0 ) , A ( t ) > 0 , t > 0 . Считая, что к моменту t* произошла дегазация образца ( c(t , x ) ≈ 0 ,
t ≥ t* ), определим константу A0 . В силу симметрии и материального баланса
t*
0
0
0
S* ≡  J (τ ) dτ =  [c − A0 ( x −  0 )2 ] dx =
(4)
3
0 0
= c  0 − A  / 3.
Отсюда A0 = 3 (c  0 − S* ) /  30 . Кроме того, условия согласования начальных данных и граничных условий дает D0 / b0 = f ( Eb − E D ) :
D(0) A0 = b(0)[c − A0 20 ]2 
D0
A0 =
b0
(5)
= exp{( E D − Eb ) /[ RT0 ]}[c − A  ] .
Теперь воспользуемся квадратичной аппроксимацией. Из условия материального баланса
выразим B (t ) > 0 и подставим в c~(t , x ) :
0
0
c~ ( t , x ) dx ,
ϕ ( x ) dx − S ( t ) =
2 2
0 0


0
S (t ) ≡

t
0
J (τ ) d τ 
0
Заметим, что параметры начального распределения A0 , c входят в соотношение лишь неявно посредством J (t ) .
Преобразуем теперь Y с учетом связи
D0 / b0 = f ( E X ) (см. (5)):
c  0 − A030 / 3 − S (t ) = B(t ) 0 − A(t )30 / 3,
B(t ) = A(t )20 / 3 + c − A020 / 3 − S (t ) /  0 ,
[
]
 c~(t , x ) = Q (t )  − A(t ) x 2 − x +  2 6 ,
t*
Q (t ) ≡ 2  J (τ ) dτ .
t
Подставим теперь выражение для аппроксив
граничное
условие
мации
c~ (t , x)
2
~
~
выразим
функцию
D (T ) c x (t ,0) = b(T ) c0 (t ) ,
A ( t ) > 0 через параметры модели D , b и известную по экспериментальным данным Q (t ) .
Оба корня квадратного уравнения относительно
A(t ) положительные, выбираем меньший из
них (по физическому смыслу c0 (t ) ≥ 0 ). Из ус-
J = c~0 b получаем соотношение для
оценки D 0 , E D , b0 , E b :
ловия

2 b (T )
3 D (T ) 
J (t )
=
− 1  , (6)
 1 + Q (t )
3 D (T )
b (T )
 b (T ) 

,
.
T = T (t ) t ∈ [0, t* ]
Поскольку J (t ) соответствует исходной модели (1)−(3), а на предварительном этапе используется параболическое приближение концентрации H в объеме, то это равенство является приближенным.
График J (t ) растет, потом убывает, причем
на начальном и конечном этапах измерения менее точны. Поэтому нормируем уравнение на
I max = J max (по порядку) и выделим безразмерные переменные:
−1
I ( t ) I max
=
( 1 + 2Q (t )t
I (t ) ≡ J (t ) , X ≡
−1
*
)
−1
J max
X −1 Y ,
t* J max b(T )
3D(T )
,
, Y≡
3D (T )
I max  b(T )
t ∈ [t1 , t2 ] ⊂ (0, t* ) .
Формально допуская
удобно считать новые переменные
X (t ) ≡ X (T (t )) , Y (t ) ≡ Y (T (t )) «аррениусовскими»:
E <0,
X0 ≡
t* J max b0
3 D0
, Y0 ≡
,
3 D0
I max  b0
E X = E b − E D , EY = E D − E b / 2 .
Тогда
−1
f (t; X 0 , E X , Y0 , EY ) ≡ I (t ) I max
−
(
)
−1
− 1 + q(t ) X − 1 Y = 0, q ≡ 2Qt*−1 J max
.
(7)
Y = Y0 exp{− EY /[ RT (t )]} =
= Z 0 exp{− E X /[ RT0 ]}exp{− EY /[ RT (t )]} ,
3 A0 c 2  20 b0
, Z 0 = Z 0 ( b0 ) ↔ b0 .
Z0 ≡
4 I max
Разумеется, I max зависит от всех входных
данных {ϕ , D, b} . Запись Z 0 = Z 0 (b0 ) означает,
что значения c , A0 уже найдены, а J (t ) при решении обратной задачи воспринимается как заданная фиксированная функция времени. Аналогично представим X :
X = X 0 exp{− E X /[RT (t )]} =
t* J max A0
exp{E X /[RT0 ]}exp{− E X /[RT (t )]}.
[c − A020 ]2
Подставляя в уравнение (7), получаем зависимость f = f (t; Z 0 , E X , EY ) . Далее с учетом зашумленности реальных измерений J (t ) и погрешности параболической аппроксимации целесообразно следовать методу наименьших
квадратов (МНК):
t
F ( Z 0 , E X , EY ) =  f 2 (τ ) dτ → min , [t1 , t2 ] ⊂ (0, t* ).
=
2
t1
Производные функции F по указанным аргументам можно выписать явно (подсчет интеграла по квадратурной формуле считаем
элементарной операцией). Задача оптимизации решается стандартными средствами (авторы использовали пакет Scilab) в физически оправданном диапазоне параметров (в пределах
одного-двух порядков). Укажем принятые в
данной работе опорные значения:
c = 1018 1/см3; T0 = 300 К;  = 10 −2 см;
4
b0 = 10 −3 см /c; Eb = 90 кДж/моль;
2
D0 = 10 −2 см /c; E D = 40 кДж/моль.
Результаты вычислительных экспериментов
показали, что без учета погрешностей измерений и модели (т. е. J (t ) определялся решением
прямой задачи (1)−(3)) погрешность оценивания
параметров находится в пределах 20 %. Начальные приближения энергий активации ED , Eb в
диапазоне нескольких десятков кДж/моль можно
указать для конкретного материала из физикохимических соображений. Приближение b0
( Z 0 ) берем в силу
Заменим в скобке
[exp{− μnγ (t,0)}− 1] + 1 :
[]
экспоненту
на
t
c0 (t ) = c −
J (0) = b(T0 )c02, (0) = b0 exp{−Eb / RT0 }[c − A020 ]2 .
Именно как начальное приближение в реальном
ТДС-эксперименте значение J (0) известно с
большой погрешностью.
Замечание. Задача несколько усложняется,
когда равновесную концентрацию c также приходится считать неизвестной. Тогда задача четырехмерная в соответствии с (7). По оценкам
X 0 , E X , Y0 , EY значения c , A0 находятся из
уравнений (4), (5).
Далее переходим к локальному уточнению
оценок D0 , ED , b0 , Eb в соответствии с исходной
моделью (1)−(3). При этом учитываем, что искомых независимых переменных 3 в силу
D0 / b0 = f ( Eb − ED ) .
Функция Грина
Поскольку зависимость от времени J (t ) известна по результатам ТДС-эксперимента, решение краевой задачи удобно представить с помощью функции Грина [Мартинсон, Малов,
2002, гл. 2]. Для ϕ ( x) = c имеем
ν n (t ) ≡
∞
2
4 ∞
J (τ ) dτ + 4 A0 ν n (t ) −  J n (t ),

0
 n =1
n =1
1 − exp{− μ n γ (t ,0)} ,
μn
t
J n (t ) ≡  exp{− μnγ (t ,τ )}J (τ ) dτ .
0
Соотношение
J = bc0 (t ) имеет форму
семейства уравнений для оценки параметров:
Φ = Φ(t; D0 , E D , b0 , Eb ) = 0 ,
t ∈ [t1 , t 2 ] ⊂ (0, t* ). При численной реализации ряды заменялись частичными суммами:
(8)
Φ = J (t ) − b ×
t
N1


2
4 N2
ϕ (0) −  J dτ + 4 A0 ν n (t ) −  J n (t ) = 0.
0
 n =1
n =1


Используя пакет Scilab, численно решалась за-
дача

t2
t1
Φ 2 dτ → min при N1 = N 2 = 5 .
× {G ( x, t ,0,τ ) + G ( x, t , ,τ )}dτ ,
За начальное приближение принимались значения, полученные в предыдущем пункте. Уровень ошибок оценивания в среднем понизился
до нескольких процентов (на «идеальном» сигнале J (t ) из (1)−(3)). Численные эксперименты
показали, что при N i ≥ 10 невязка резко возрастает, возможно, из-за накопления вычислительных погрешностей.
G ( x, t , y ,τ ) =
Сопряженные уравнения
t
c(t , x ) = c −  D −1 (τ )J (τ ) ×
0
=
n π
1 2
+  exp  2
  n =1
 
∞
2
2

(τ − t ) cos nπx cos nπy .



Для уточнения оценок параметров модели
используем соотношение J = bc02, , поэтому
далее нас интересует представление граничной концентрации c0 (t ) . Для варианта
ϕ ( x ) = c − A( x −  0 ) получаем:
t
2
c0 (t ) = c −  J (τ ) dτ −
0
 4 ∞ 1

− A0   +  exp{− μnγ (t ,0)} −
12  n =1 μn

Если неопределенность значений D, b очень
велика (например, 3–5 порядков), то нелишне
иметь в распоряжении дополнительное семейство уравнений, поскольку на универсальный алгоритм параметрической идентификации нелинейной распределенной модели рассчитывать не
приходится. Укажем один из возможных подходов, основанный на общей идее использования
сопряженных уравнений [Марчук, 1992].
Используя интегрирование по частям, для
произвольной достаточно гладкой функции
ψ (t , x) получим:
t* 
t
−
4 ∞
 exp{− μnγ (t,τ )}J (τ ) dτ ,
 n =1 0
t
τ
γ (t ,τ ) ≡ D ( s ) ds , μ n ≡ (nπ  0 )2 .
0 =  ψ (t , x )[ct − Dc xx ]dxdt =
0 0
t*
t*
0
0
=  J ( t )ψ ( t ,0) dt +  J ( t )ψ (t , ) dt +
t*
+  D (t ) J (t )b −1 (t )ψ x (t , ) dt −
0
t*
−  D (t ) J (t )b −1 ( t )ψ x (t ,0) dt −
0


0
0
− c ψ (0, x ) dx + A ( x −  0 ) ψ (0, x ) dx .
2
(9)
Здесь опущен двойной интеграл, поскольку
дальнейшем изложении считаем функцию
ψ (t , x) подчиненной сопряженному уравнению
в
∂ψ ∂t = − D ∂ 2ψ ∂x 2. Подчеркнем, что краевые условия не ставятся, «пробных» функций ψ бесконечно
много. Выберем, например, ψ (t , x ) = β (t ) expσx .
При нормировке β (t* ) = 1 получаем
ψ ( t , x ) = exp{σ 2γ ( t* , t )} exp σx ,
γ (t,τ ) ≡

t
0
D ( s ) ds .
Перепишем (9) в обозначениях
t*
t*
0
0
X ≡  Jβ dt , Y ≡  D Jb −1 β dt :
exp σ + 1
β (0)
, κσX + σ 2Y + 2 ×
κ≡
exp σ − 1
σ
× {2 A0 − σ 2 (c − A0 [ 20 − κ σ ])}= 0.
Аналогично
рассматривается
вариант
2
:
ψ (t , x) = β (t ) sin σx (β (t ) cosσx ) ψ (t , x) = exp{− σ γ (t* , t )}sin σx
( ψ (t, x) = exp{− σ 2γ (t* , t )}cosσx ) ,
κσX + σ 2Y +
β ( 0)
×
σ2
× {2 A0 + σ 2 (c − A0 [  20 − κ σ ])}= 0,
sin σ 
cos σ + 1  .
κ≡
κ ≡

cos σ − 1 
− sin σ 
Получаем дополнительное семейство уравнений. Параметр σ целесообразно варьировать в
пределах σ ~ 1 .
Заключительные замечания
В предлагаемой итерационной процедуре оценивания используется подсчет интегралов по
времени, а не численное решение краевой задачи
при текущих приближениях параметров (как того
требует, например, «стандартный» градиентный
метод минимизации невязки). Если мембрана не
очень тонкая, следует использовать вариант
ϕ ( x ) = c − A0 ( x −  0 ) 2 k , k > 1
(в середине пластины равновесная концентрация и только у самого края концентрация немного меньше).
При обработке зашумленных экспериментальных кривых следует, по-видимому, отказаться от использования условия согласования
краевых условий (5). В интегральном смысле
(рассматриваем решение краевой задачи (1)−(3)
как обобщенное) считаем ϕ ( x) = c ( A0 = 0 ). В
силу (4) имеем оценку равновесной концентрации c . Далее считаем параметры X 0 , E X , Y0 , EY
(однозначно определяющие D0 , E D , b0 , Eb ) независимыми в соотношении (7). Последующие
преобразования X , Y не проводим, уравнение
(5)
и
соответствующее
соотношение
D0 / b0 = f ( Eb − ED ) исключаем из модели. Задача оптимизации становится четырехмерной.
Работа выполнена при поддержке РФФИ
(грант 09-01-00439).
Литература
Взаимодействие водорода с металлами / ред.
А. П. Захаров. М.: Наука, 1987. С. 177–206.
Водород в металлах / ред. Г. Алефельд, В.
Фёлькль. М.: Мир, 1981. Т. 1. 506 c. T. 2. 430 c.
Кунин Л. Л., Головин А. И., Суровой Ю. И.,
Хохрин В. М. Проблемы дегазации металлов. М.:
Наука, 1972. 324 c.
Мартинсон Л. К., Малов Ю. И. Дифференциальные уравнения математической физики. М.: МГТУ
им. Н. Э. Баумана, 2002. 368 c.
Марчук Г.И. Сопряженные уравнения и анализ
сложных систем. М.: Наука, 1992. 336 с.
Самарский А. А., Гулин А. В. Численные методы.
М.: Наука, 1989. 432 с.
СВЕДЕНИЯ ОБ АВТОРАХ:
Заика Юрий Васильевич
зав. лаб. моделирования природно-технических систем, д. ф.-м. н.
Институт прикладных математических исследований
КарНЦ РАН
ул. Пушкинская, 11, Петрозаводск, Республика Карелия,
Россия, 185910
эл. почта: zaika@krc.karelia.ru
тел.: (8142) 766312
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
Документ
Категория
Без категории
Просмотров
7
Размер файла
303 Кб
Теги
диффузия, оценки, дегазации, водорода, тдс, краевой, десорбции, задачи, параметры
1/--страниц
Пожаловаться на содержимое документа