CC BY
25
Труды арельского нау ного центра
8. 2016. С. 11-23 DOI: 10.17076/mat395
519.6:539.2
АППРОКСИМАЦИЯ КРАЕВОЙ ЗАДАЧИ ТЕРМОДЕСОРБЦИИ СИСТЕМОЙ ОДУ: УСКОРЕНИЕ СХОДИМОСТИ
Ю. В. Заика, Е. К. Костикова
Ин тит т прикладны ате ати е ки и ледовани арел кого на ного ентра Р
рамка те нологи ески ада водородного материаловедения (вкл ая проект ITER) ведется интенсивны поиск ра ли ны по на на ени конструкционны материалов с аранее аданными пределами водородопроницаемости.
дним и кспериментальны методов является термодесорбционная спектрометрия (ТДС). браец, насы енны водородом, дегаируется в условия ваку-умирования и монотонного нагрева. С помо ь масс-спектрометра регистрируется десорбционны поток, по воля и судить о арактере в аимоде ствия и отопов водорода с твердым телом. нтерес представля т такие параметры переноса, как ко ициенты ди у ии, растворения, десорбции. статье представлены распределенная краевая ада а термодесорбции и исленны метод моделирования ТДС-спектра, требу и лишь интегрирования нелине -но системы обыкновенны ди еренциальны уравнени ( Д ) невысокого порядка. ло ен метод ускорения с одимости на основе выделения интегрируемо слабо особенности.
л евые слова: водородопроницаемость; термодесорбция; нелине ные краевые ада и; динами еские грани ные условия; исленное моделирование.
Yu.V. Zaika, E. K. Kostikova. APPROXIMATION OF THE BOUNDARY-VALUE PROBLEM OF HYDROGEN THERMAL DESORPTION BY ODE SYSTEM: CONVERGENCE ACCELERATION
One of the technological challenges for hydrogen materials science (including ITER project) is the currently active search for structural materials with various potential applications that will have predetermined limits of hydrogen permeability. One of the experimental methods is thermodesorption spectrometry (TDS). A hydrogen-saturated sample is degassed under vacuum and monotone heating. The desorption ux is measured by mass spectrometer to determine the character of interactions of hydrogen isotopes with the solid. We are interested in such transfer parameters as the coe cients of di usion, dissolution, desorption. The paper presents a distributed boundary value problem of thermal desorption and a numerical method for TDS-spectrum simulation, where only integration of a non-linear system of low order ordinary di erential equations (ODE) is required. The method of convergence acceleration based on the allocation of integrable weak singularity is explained. This work is supported by the Russian Foundation for Basic Research (project 15-01-00744).
Keywords: hydrogen permeability; thermal desorption; nonlinear boundary-value problems; dynamical boundary conditions; numerical simulation.
Е ЕНИЕ
нтерес к взаимодействию изотопов водорода с различными материалами носит многоплановый характер [1, 2, 5, 6, 10, 12 20].
нтузиасты говорят не только о водородной энергетике, но и о водородной экономике [19]. Некоторые частные задачи водородного материаловедения, связанные с темой статьи, исследованы в [21 24]. кс-периментальный опыт показывает, что лимитирующими являются не только диффузия, но и физико-химические процессы на поверхности [1, 2]. Остановимся на моделировании термодесорбции, учитывая лишь основные факторы и информативность рассматриваемого ТДС-эксперимента. Статья является продолжением [4], поэтому представим постановку задачи и модель компактно.
А ЕА ИЧЕ КА ЕЛ ЕЕНА
Рассмотрим перенос водорода сквозь образец тестируемого материала (пластину толщиной £). Нагрев медленный, практически равномерный, так что градиент температуры пренебрежимо мал и диффузионный поток можно считать пропорциональным градиенту концентрации. Концентрация растворенного водорода мала. Материал достаточно однороден, чтобы пренебречь взаимодействием H с ловушками (микродефектами структуры, которые могут удерживать водород). Для диффузии примем стандартное уравнение:
ct(t,x) = D(T )cxx(t,x), (t,x) Qtt, (1)
где t время, Qtf =(0,t) (0,£); c(t,x) концентрация диффундирующего водорода (атомарного). Зависимость коэффициента диффузии D от температуры T(t) соответствует закону ррениуса: D = D0 exp{ ED/[RT(t)]}.
Трудности численного анализа связаны не с уравнением (1), а с динамическими нелинейными граничными условиями. Пусть пластина контактирует с газообразным H2 и поверхность является потенциальным барьером (см. [1, с. 177 206; Габис, Компаниец, Кур-дюмов]). Тогда с учетом (де)сорбционных процессов краевые условия следуюшие:
c(0,x) = c(x), x [0,4 t [0,t ], (2)
co(t) = g(T)qo(t), ci(t) = g(T)q£(t), (3)
qo(t)= MT)po(t) b(T)q2(t) + Dcx(t, 0), (4)
Ш= №(T)pe(t) b(T )q|(t) Dcx(t,£), (5)
b(T) = bo exp{ Eb[RT] , s(T) = ..., где c0(t) c(t, 0), ci(t) c(t,£) граничные об -емные концентрации диффундирующего водорода; q0(t), qi(t) концентрации на поверхностях (x = 0,£); g(T) параметр локального равновесия между концентрациями на поверхности и в приповерхностном об еме; ß кинетический коэффициент; s(T) параметр, отражающий тот факт, что только малая часть «налетающего» водорода окажется в форме атомов на поверхности; p0(t), pi(t) давления газа (H2); b(T) коэффициент десорбции. Говоря о потоках (в частности, десорбционных J0,i = bq0i), подразумеваем их плотность.
В модели фигурирует как молекулярный, так и атомарный водород. Для единообразия подсчет будем вести в атомах: [c] = 1/см3, [q] = 1/см2, [Dcx] = [J] = 1/см2с (J = bq2). В кинетической теории газов величина ßp определяет число частиц (в данном случае молекул H2), соударяющихся с единичной площадкой поверхности в единицу времени. Но за счет множителя s удобно воспринимать слагаемое ßsp как плотность потока атомов, оседающих на поверхности. то интегральный показатель, без разделения процесса на стадии. Вместо s можно написать 2s и интерпретировать s как долю адсорбируемых атомов H.
В масштабе рассматриваемого далее ТДС-эксперимента удобно выбрать [p] = Topp, откуда ß(T) = (2nmkT) 1/2 2, 484 1022/ T. Здесь [ß] = 1h2/(Toppсм2с), [T] = K, k постоянная Больцмана, m масса молекулы водорода. Зависимостью кинетической константы от температуры (ß 1/ T) часто пренебрегают на фоне экспоненты в s(T).
аеа е. В случае высокой степени пористости разумно уравнения (3) (5) заменить на
ßs(T)po,i(t) b(T)c0,i(t)= D(T)cx|o,i. (6)
Плотность потока абсорбции ßsp (в приповерхностный слой) пропорциональна давлению H2 снаружи. Формально (6) получается из уравнений (4), (5) при малой скорости накопления на поверхности (q 0). Коэффициент десорбции обозначается одной буквой, хотя в (6) и (4), (5) это разные величины. х согласованность понимаем в смысле bsurface = g2bvolume (см. (3)). В [10] приводится подробный анализ различных стадий проникновения водорода из газа в металл и обратно. Коэффициент b в условии (6) является эффективным коэффициентом рекомбинации.
е т
оде
е е
В камеру c нагретой пластиной из исследуемого металла или сплава подается водород в газовой фазе при сравнительно большом давлении. После того как образец поглотит достаточное количество водорода (до состояния равновесного насыщения), он быстро охлаждается (отключается ток нагрева). В режиме последующего постоянного вакуумирова-ния камеры лента снова нагревается. С помощью масс-спектрометра измеряется давление молекулярного водорода в вакуумной камере, обусловленное десорбцией J(t) = b(t)q2(t). «Перерасчет» p() J() определяется характеристиками установки, выходными данными считают функцию J(t). Для упрощения обозначений примем сокращенную запись:
b(t) b(T(t)), D(t) D(T(t)), s(t) s(T(t)).
Коэффициент g быстрого растворения (локального равновесия «поверхность об ем») считаем константой (Ek- = Ek+) в пределах рассматриваемого ТДС-пика десорбции. Для метода ТДС выполнена симметрия:
p(t) = po,e(t), q(t) = qo/(t), co(t) = ce(t), (7)
D(t)cx(t, 0) = D(t)cx(t,£), c(x) = c = const.
Время t окончания эксперимента определим условием практически полной дегазации: p(t) 0,t ^ t , c(t ,x) 0,x [0,£].
HE -И EEH ИАЛ H Е А НЕНИЕ Е Е ИИ
Принятая модель ТДС-дегазации: ct = D(t)cxx , c(0,x) = c, co,e(t) = gq(t), q(t) = r(t) b(t)q2(t) + D(t)cx(t, 0), r(t) ps(t)p(t), J(t) b(t)q2(t).
Вакуумную систему считаем достаточно мощной, чтобы после насыщения, на этапе дегазации пренебречь ресорбцией (r(t) = 0) в динамическом граничном условии (q = ...). Разрешимость краевой задачи (и в более общем случае учета обратимого захвата в об еме) обоснована в [3]. Ограничимся прямой задачей численного моделирования ТДС-спектра J = J (T). тот этап необходим как итерационная составляющая алгоритма идентификации. Нагрев T(t) обычно реализуют линейным (T (t) = To + vt). Скорость нагрева v невелика (< K/c). По достижении максимальной температуры (если дегазация еще не завершилась) нагрев прекращается: T(t) = Tmax.
Теперь мы в состоянии математически сформулировать задачу. В [3, 22] представлены разностные схемы решения краевой задачи (в том числе и с учетом ловушек различных типов). Но чтобы сравнивать модельные и экспериментальные ТДС-спектры, нужна только поверхностная концентрация (1 = Ьд2). Естественно попытаться избежать итерационного решения краевой задачи при текущих приближениях параметров модели Ер, Ь0, Еь, в0, Е3, д. С этой целью проведем преобразования, придя в итоге к необходимости интегрировать лишь систему ОДУ невысокого порядка.
Оставив прежнее обозначение Ь, выполним замену времени Ь = ^^йв :
Сг{Ь, х) = ехх(Ь, х), с(0, х) = с, Со/ = дд(Ь), (8)
Схо = Сх I = 4(1) + [1 (Ь) т(Ь)]Б 1(Ь). (9)
Считаем д(Ь) функциональным параметром, а (9) дополнительным соотношением к линейной задаче (8). Сделаем замену, приводящую краевые условия в (8) к однородным:
с = с(Ь, х) дд(Ь), Сг(Ь, х) = Схх(Ь, х) + /(Ь), /(Ь) = дд(Ь), С(0, х) = ф(х) = 0, С о,е = 0.
Запишем решение с помощью функции мгновенного источника (функции Грина) [8, гл.2]:
С(Ь,х)= Г [£С(х,С,Ь т)/(т) й£йт,
>0 J о 2
J
. . 2 v-^ г n2n2 л
G(x,i,t) = jl^exP{ ~jrt)
sin ■
nnx .
n=1
£
sin ■
£
В динамические граничные условия входит производная cx(t, 0):
n2 n2
cx 0 = y Jo q(r£2
(r t^ dr,
Е Еп=1,з,5.... При т = Ь ряд расходится, так что подразумевается почленное интегрирование. В исходном времени Ь имеем
Сх(Ь, 0) = Сх(Ь, 0) = Сх(Ь,£) = Сх(Ь,£), Сх(Ь, 0) =
Г q
t - п2п2 *
£ q(r)exp| \ D(s) ds}dr.
Окончательно динамическое граничное условие запишется в форме
q(t) = r(t) b(t)q2(t)
4gD ST ft •
(10)
( n2n2
q(r)exp| J D(s) dsjdr.
Jo
Уравнение эквивалентно исходной краевой задаче: решение q(t) однозначно определяет c(t,x). Нас интересует отрезок времени [t 1, (0,t ), соответствующий выражен-
ному пику термодесорбции (измерения при t близких к 0, t малоинформативны). Для мощной вакуумной системы ресорбцией пренебрегают (r(t) = 0), что упрощает уравнение.
еаеа оа ада
Для численного моделирования удобно перейти к безразмерным переменным, выполнив замены: t = J^D^ds/^2, x = x/£, v = q/с (с = gc). Не меняя обозначения t, получаем:
b(t)£2ç
b(t)
V(t) = b(t)v2(t)
(11)
4gl
EУ vv(r) exP { n2n2[t t]}dr,
г>(0) = 1 (начальные данные). Ограничимся первыми к слагаемыми ряда. Обозначим:
**(*)= / ¿(т) ехр { (2г 1)2п2[^ т]} ^т.
Jo
Продифференцируем по Ь и подставим вместо V выражение (11) (сумма до п = 2к 1):
Zi = bv
[(2г 1)2п2 + 4g^Zi
4дф1 + ... + 1 + ¿¿+1 + ... + гк).
В итоге получаем следующую систему обыкновенных дифференциальных уравнений:
d_ dt
fv\
Zi \zk)
= b(t)v2 (t)
1 1
(12)
4gl
1
11 ai 1
\1 1
1 1
O-kj
/Zi\
Z2
Zk
ai 1 + (2i 1)2n2/(4gl), v(0) = 1, Zi(0) = 0 (размерность матрицы (k + 1) k).
И ЛЕНН E
ЕЛИ
АНИЕ
Для численных экспериментов использовались данные по никелю и нержавеющей стали марки 12Х18Н10Т. Варьирование параметров позволяет оценить степень чувствительности («производные») ТДС-спектра к выделенным лимитирующим факторам. Для совпадения результатов с решением исходной краевой задачи разностным методом [22] потребовалось
интегрирование системы (12) порядка 10 15 (авторы пользовались пакетом Scilab 5.5). Результаты моделирования представлены в [4].
Данные по никелю [5]: [E] = Дж/моль; D = 7, 5 10 3 exp 40 000/RT [см2/с]; bvol = 1, 53 10 18 exp 43 200/RT [см4/с]; s = 1, 8 10 2 exp 61400/RT ; p = 37,4[Торр]; l = 0, 02 [cm]; T = 770 [K], T(0) = 295 [K]; T = 0,5[K/c]. Для определенности зафиксируем bo surf = 1, 53 10 14 [см2/с]; g = 100[1/с].
Данные по стали марки 12Х18Н10Т [5, 18]: D = 3,09 10 4 exp 27 780/RT ; bvol =
5, 05
10
13
exp 97140/RT ; s = 7,05
10 2 exp 59 510/RT ; T(0) = 300 К. Для определенности b0surf = 5, 05 10 9; l = 0,1; g = 100; p = 100; T = 850, T = 0, 5 [K/c].
Л
ЕЛЕНИЕ К ЕНИ
ЕНН И
И И
Проанализируем, почему при достаточно точной аппроксимации интегро-дифференци-ального уравнения (11) системой ОДУ (12) требуется ввести относительно большое количество переменных ¿¿(¿). Функция
Ö(s) = 4^ exp {
n2n2s
} E E
«.= 1,3,5,...
имеет конечные значения при в > 0. Ряд быстро сходится при больших в, а при формальной подстановке в = 0 (когда переменная интегрирования т достигает верхнего предела Ь) получаем расходящийся ряд. «Спасает» почленное интегрирование. При этом тах ¿¿(¿) = 0(п 2) (п = 2г 1, г ^ 1), что приводит к медленной сходимости. Десяток-другой ОДУ указанного выше типа не является вычислительной проблемой, но является препятствием для использования метода специалистами в области водородного материаловедения. Поставим задачу найти компромисс: ОДУ должно быть немного (4 5), чтобы можно было использовать стандартные пакеты программ (например, свободно распространяемый БсПаЪ). Выделим интегрируемую особенность и соответствующим образом «подправим» параметры в уравнениях. Условно будем вести речь о методе эффективных коэффициентов.
Проведем преобразования, используя тэта-функции коби. Точнее, нас интересует
#3(Ь,х) = 1 + 2^ ехр п2п2Ь ео8(2ппх).
При х = 0, Ь > 0 имеем альтернативное представление 03(Ь, 0) (см. [6, с. 179]; [11, с. 353]):
1 + 2 E exp n2n2t = —= E exp Г П
n=1
t
14
Ряд слева быстро сходится при больших Ь, а ряд справа при малых значениях Ь (которые нас сейчас и интересуют). Если определить 0(Ь) = ^ехР пп2Ь (п Z, Ь > 0), то получим функциональное уравнение для тэта-функции [7, с. 261] 0(1/Ь) = 10(Ь), или ехр пп2Ь = ^ехр пп2/Ь (п Проведем вспомогательные преобразования, разбивая промежуточную сумму на нечетные и четные слагаемые:
1
= У^ехр
у} = 1 + 2^ехр п2п2Ь =
п=1
1 + 2^ + 2^ехр к2п24Ь =
к=1
2У + —=Уехр
п2 4
Отсюда (после вычитания из первого ряда последнего и удвоения) получаем искомое разложение при в > 0:
в(в) ехр {
22
1 + 5
п"п"в\ = —=-, пв
Б (в)
2
п=1
дп , д
ех^ 4в
Ряд 8 очень быстро сходится при малых в. При в +0 имеем 8 0 и интегрируемую особенность в 1/ пв. Ряд еще и знакопеременный, что влечет неравенства 8 < 0 (точнее Б ( 1,0)) и Б ^ 2 д . Последняя оценка грубая (по сравнению с Б2т ^ Б ^ Б2т+\ в терминах частичных сумм), но нам ее будет достаточно при относительно малых в. График функции Б (в) имеет Б-образный вид кривой насыщения и уже Б(1) 0, 999. В получаемом с помощью описанного представления в(в) уравнении {р ^ 0, )(0) = 1) )(Ь) =
= Ь(ь)у2(ь) к/ * 1+1Б ( Т) )(т) йт (13)
дробь (слабая особенность под интегралом) быстро убывает от бесконечности (т = Ь) до практически нуля (т = Ь 1). Но не учитывать более отдаленную предысторию опасно, поскольку в общем случае (Ь 1, к 1) функция )(Ь) быстро убывает с ростом Ь.
Оценим, насколько приемлемым (в контексте ТДС-эксперимента) является приближение в(в) 1/ пв, которое элементарной формулой выделяет интересующую нас особенность ( Б 1 при малых в). Для определенности зададимся уровнем 2% ошибки в числителе: 1 + Б 1, Б < 2ехр 1/4в < 2 10 2.
Округленно получаем в ^ 1/20. В экспериментальных работах точность измерений в ТДС-эксперименте априорно оценивается в пределах 10 20 %. то позволяет в рамках принятой модели уравнение (11) переписать в форме
ЦЬ) = ъ(г)ь2(г)
ГЬ к
4кЕ / ■■■ 4кЕ
(14)
Ч к
(Ь ^ Н, Ь дЬ£2/Б, к д£) и в последней сумме («главной части») воспользоваться аппроксимацией в(в) 1/ пв (в = Ь т) с приемлемым ограничением 0 ^ в ^ Н = 1/20. та оценка в безразмерном («штрихованном») времени. Возвращаясь к исходному физическому времени, имеем
Н // н°й^/£2. Типичные значения £ = 0, 02 см, Б = 10 7(10 6) см2/с, что приводит к грубой оценке Н 200(20) с. Время Н не является исчезающе малым в масштабе часов ТДС-эксперимента. Таким образом, изложенный способ выделения главной части интеграла не приводит к дополнительным особенностям математического характера при обработке экспериментальных данных.
Метод етв о е тов
В исходном уравнении (11) выделим в сумме небольшое число слагаемых (для определенности три с учетом быстрого роста показателя в экспоненте п2 = 1, 9, 25):
1 + 2,
Е
п=7,9,
Для 1 введем, как и ранее, переменные гг(Ь) = )(т) ехр { п2п2[Ь т]} йт, п = 2г 1,
Jo
г = 1, 2, 3; г^0) = 0, Ь > 0 Zi < 0. Сумму 2 аппроксимируем, воспринимая ее как возмущение. При этом выделим интегрируемую особенность в(в) 1/ пв (в = Ь т 1) и будем чередовать завышающие оценки с занижающими. ель: остаться в классе ОДУ невысокого порядка (ценою соответствующей корректировки коэффициентов уравнений).
Фиксируем параметр Н «опасного сближения» т и Ь (Ь т ^ Н) и рассмотрим начальный отрезок времени Ь [0, Н]. Оценка 1/ пв немного завышена (Б < 0, Б 1), поверхностная концентрация монотонно убывает при вакуумной дегазации () < 0). Поэтому
4к^ ^ 4к^
1
^ к у(т) -
Ь \/п[Ь т ]
йт. (15)
г
2
1
Считая погрешность оценки малой, вместо уравнения (11) получаем
v
v(t) = b(t)v2(t) 4 к Е Zi(t)
i=1,2,3 v(r )
ft
О л/nit r]
dr, t [0,h]. (16)
wfc+i(t) = b(t)v2(i) к
f* wfc+i(r) О л/П^ r]
dr,
соответствует физическому смыслу: в исходном уравнении (11) мы убрали в правой части положительные слагаемые, что уменьшает ), но множитель [...] это компенсирует.
С учетом ¿¿(¿) = п2п2£(Ь) + )(£} получаем следующую систему ОДУ () = )(Ь), £ = ¿¿(Ь)):
Переменные г»(Ь) введены для учета «эффекта накопления» к последующему интегрированию на интервале Ь > Н. При малом Н и начальных данных £¿(0) = 0 ими можно пренебречь (даже уменьшив погрешность, используя сразу правое неравенство в (15)):
)(Ь)= Ь(Ь))2(Ь) к Г /(Т) = йт, (17)
() () () Ус л/П[Ь ту , у 7
)(0) = 1, Ь [0, Н]. Для переменной -ш(Ь) = )(Ь) ()(Ь) = 1 + /С эд^т) получили нелинейное интегральное уравнение с так называемой слабой особенностью. Если иметь в виду итерации )с(Ь) = 5(Ь))§(Ь) ()с(0) = 1), адс = Ь)2,
с v = bv2
¿1 = bv2
¿2 = bv2
-¿3 = bv2
z1 = bv2 [4 ä + п2] z1 4 к z2 4 ä z3,
то можно воспользоваться развитой теорией линейных интегральных уравнений со слабой особенностью [9] или операционным исчислением (справа свертка 1/ пЬ). Такие итерации, излишне усложняющие общий вычислительный алгоритм, трудно интерпретировать физически. К тому же использовалась завышенная оценка В 1/ пз в (15).
Заметим, что за счет особенности знаменателя «наибольший вклад» в свертку функций )(Ь) 1/ пЬ вносят значения )(т), т Ь. Заменим под интегралом )(т) на )(Ь) (с погрешностью дроби 0( Ь т)). С учетом убывания абсолютной величины скорости изменения поверхностной концентрации это приводит к занижению интегрального слагаемого, что является определенной компенсацией предшествующего завышения. Здесь напомним, что плотность термодесорбции определяется как 7 = Ьд2, так что монотонное уменьшение д(Ь) () = д/д) с одновременным более интенсивным ростом Ь(Т(Ь)) по мере нагрева приводит к немонотонному ТДС-пику. В итоге, интегрируя особенность, вместо (16) приходим к ОДУ
)(Ь)[1 + 2к/Тп] = Ь(Ь))2(Ь) 4к^ ^¿(Ь).
В случае (17) сумма £ будет отсутствовать. Отметим, что коррекция коэффициента при )
(18)
для которой v(0) = 1, Zi(0) = 0, b = b(T(t)) b(T(t))/[1 + 2K^t/п], к к/[1 + 2K^iTn] (когда h относительно мало [...] 1). По сравнению с исходной системой произошла лишь корректировка коэффициентов за счет анализа особенности при т t в интегралах.
При переходе к интервалу t > h рассмотрим уравнение (14) при t ^ 0, для единообразия формально полагая v(r) =0, r < 0. В силу относительной малости h во второй сумме воспользуемся аппроксимаций 0(s) 1/ ns. В первой сумме оставим «главные» первые слагаемые. При этом отбрасываем положительную величину. В порядке компенсации верхний предел интегрирования t h заменим на близкий t. Приходим к уравнению
v(t) = b(t)v2(t) 4 к E zi(t)
¿=1,2,3
' * _v(r) Ith'
ж
Vn[t r ]
dr, t ^ 0.
Далее с учетом подынтегральной особенности заменяем ))(т) на ))(Ь) (памятуя о )(т) = 0, т < 0) и получаем систему (18). Только при Ь > Н заменяем [1+2 кл/Ь/п] на «установившееся» значение [1+2 кл/Н/п]. При использовании пакетов программ (например, БсПаЪ) можно использовать универсальную запись системы в форме (18) для Ь ^ 0 с применением функции насыщения ва^ [...] = 1+2 ку^а^/^Н/п.
Приведенные рассуждения в совокупности с исходными физико-химическими упрощениями можно подытожить в совокупности: предложена простейшая модель (18) термодесорбции с учетом динамического взаимодействия поверхностных процессов с диффузией. По функции )(Ь) = д(Ь)/<? (Ь = Ь ) после возвращения к физическому времени получаем модельный ТДС-спектр 7(Т) = Ь(Т)д2(Ь) (Т(Ь) Ь).
Прокомментируем рисунок 1. Увеличение безразмерного параметра Н ухудшает аппроксимацию особенности, но улучшает возможно-
2
2
сти приближения «хвоста» (Ь — т > К) небольшим числом переменных г% (Ь). Восходящий фронт ТДС-спектра аппроксимируется с высокой точностью в широком диапазоне К, погрешность накапливается к окрестности пика. Динамика аппроксимации пика немонотонна по К (значения указаны по убыванию тах). При этом, например, К = 0,1 лучше «выглядит на пике», а К = 1 точнее на спаде (Т > 600К). С учетом значительной погрешности определения плотности потока термодесорбции минимизация какой-либо невязки непринципиальна. Вычислительные эксперименты приводят к приемлемому диапазону 0,1 < К < 0,8. Рисунок 2 демонстрирует искомую «седловую точку»: и К = 0, 3 не мал, и г ^ 3 невелико. Более того, система (18) проявляет устойчивость параметрической идентификации в следующем смысле. Возьмем завышенное значение К = 10 (в том
числе и при пересчете на реальное время ТДС-эксперимента). Судя по рисунку 1, для компенсации возникшей существенной погрешности формально следует уменьшить запас водорода в образце (уменьшить с). Достаточно в приемлемых пределах изменить исходные значения параметров (увеличить предэкспоненту Ьо на 12% или уменьшить д на 11%), чтобы получить удовлетворительное приближение ТДС-спектра. За «эталон» принят график, полученный решением исходной краевой задачи разностным методом [22].
Рисунки 3, 4 демонстрируют возможности аппроксимации в более сложном случае, когда наблюдается двухпиковая структура ТДС-спектра. Отметим, что обычно этот эффект интерпретируется как влияние ловушек с различными энергиями связи (в модели этого нет, хотя технически нетрудно добавить соответствующие слагаемые в уравнение диффузии).
Метод шагов
Вернемся к уравнению (13), полагая (для единообразия) формально ¿(т) = 0 при т < 0. Поскольку Б(1) -0, 999, то фактически предыстория «длины» более чем Ло = 1 забывается и можно перейти к уравнению
/Щ) 10"13 200
400
800
1000
¿(ь) = —ь(ф2(ь)-к
2,+. у 1+Б(Ь —т)
л/п\ь—т\
¿(т) ¿т. (19)
Ограничимся сначала отрезком Ь [0,Л], Л ^ Л0. Особенность при т Ь (£(+0) = 0) позволяет в первом приближении заменить ¿(т) ¿(Ь). Получаем ОДУ:
-о-о- к = 1
-А-А- к = 0.1
—О- к = 0.05
Рис. 6. ТДС-спектр (никель), метод шагов
¿(Ь)[1+ ст(Ь)+2 Ь к/ П = -Ь(Ь)-2(Ь), (20)
Г* Б(в) , £(+0) Уо в
7(Т(0) 10"13_^оо
+0
В силу £ (в) = 2 £га=1 (д = - ехр —1/(4в) ) особенность в нуле подынтегральной функции несущественна. Можно численно интегрировать по в [е, Ь], заменив нижний предел на малое е > 0 (скажем, 10-4). Графики функций Б (в) и Б (в)/ в изображены на рисунке 5. Ряд Б быстро сходится при в [0,1], достаточно десятка слагаемых. При переходе к Ь > Л, учитывая в (19) только интеграл по т [Ь — Л, Ь], следует только «заморозить» коэффициент [...] при ¿(Ь): а(Ь) о"(Л), Ь Л. Результат моделирования иллюстрируют рисунки 6, 7. Если обрабатывать данные, соответствующие восходящему фронту ТДС-пика, то уже одного скалярного уравнения (20) достаточно для оценки коэффициентов модели. Для всей кривой нужно пользоваться более точным приближением, чем ¿(т) ¿(Ь).
\ N 2
\ ЭД-'^?". % = - ехр{-1/(4Л)}, N > 10
\ "="
V... -{ (Л) + 2 Л }
0 \ 0.2 0.4 0.6 0.8^--- 1 ^ Л
400 500 600 700
900 1 000
Рис. 5. Функции Б (в) и Б(в)/ в
Рис. 7. ТДС-спектр (сталь), метод шагов
Уточнение подынтегральной аппроксимации можно провести следующим образом. Фиксируем шаг Л, например Л = Л0. В уравнении (19) воспользуемся под интегралом линейной интерполяцией ¿(т) ¿(Ь) + А[Ь — т]. Параметр А является константой по отношению к переменной интегрирования т. На отрезке [кЛ, (к + 1)Л] (к ^ 0) полагаем А = —(¿(Ь) — ¿(кЛ))/(Ь — кЛ), чтобы при т = кЛ получить «предыдущее» значение ¿(кЛ).
Вместо линейной интерполяции можно воспользоваться квадратичным приближением ¿(т) ¿(Ь) + А[Ь — т] + — т]2. Условия т = Ь ¿(Ь), т = кЛ ¿(кЛ) приводят к уточнению: Ь [кЛ, (к + 1)Л], т [кЛ,Ь],
¿(Ь) — ¿(кЛ) г
»<т> «м- ,—к,.'[Ь—т]—
— кЛ][Ь — т ] + В[Ь — т ]2.
г
т
1200
1400
Т
Л
®
Выбор B < 0 позволяет учесть вогнутость функции v (v < 0). Подчиним B условию гладкой склейки уравнений: ¿((k + 1)h 0) = ¿((k + 1)h + 0). Здесь не будем приводить уточняющие математические выкладки.
В заключение приведем реализацию метода шагов при малом h, когда можно пользоваться приближением ©(s) 1/ ns во втором интеграле в правой части уравнения (14). Ограничимся линейной интерполяцией ¿(т)
¿(t) + A[t т], A = (¿(t) ¿(kh))/(t kh).
Фиксируем шаг h безразмерного времени (t = t ) для удовлетворительной аппроксимации ©(s) 4 J2 ехР { n2n2s} 1/ ns (s (0, h], s = t т). Воспользуемся этой оценкой на начальном отрезке t [0, h] в исходном интегро-дифференциальном уравнении (11):
V =
bv2 4^ £/.••
V(t) = b(t)v2(t) к
v(r)
о л/nït Т]
dr.
В свертке ¿(t) 1/ nt вторая функция является оценкой, так что разумно вместо «точного» решения аппроксимировать и первую функцию. Поскольку из-за особенности наибольший вес имеют значения ¿(т) при т t, то после замены ¿(т) на ¿(t) приходим к задаче ¿(t) = b(t)"U2(t) (¿(0) = 1), где b = b/ [1 + 2 кл/t/n]. В квадратурах получаем конкретные функции ¿1,г) 1. Подставляя ¿1 под интеграл, получаем уравнение Риккати г? = ¿¿2 + /1 (t). терации можно повторять. Но в итоге получим алгоритм, сравнимый по «энергоемкости» с решением исходной краевой задачи разностными методами. В целях упрощения воспользуемся под интегралом линейной интерполяцией ¿(т) ¿(t) + A[t т]. С учетом W(0) = b(0) получаем значение A = [¿(t) + b(0)]/t и уравнение Риккати
V(t)
_ 4 к t
1 + 4Tn-J
= b(t)v2(t) +
2 кЬ(0) t
3 п '
г(0) = 1, £ [0, Л]. В физических терминах уравнение г = Ьг2 соответствует быстрой десорбции водорода с поверхности без диффузионной «подкачки» из об ема. Дроби корректируют уравнение с учетом растворения и диффузии ( к = £ = £ (^)): формально уменьшается коэффициент десорбции и появляется положительное слагаемое в правой части.
Переходим к отрезку t [kh, (k + 1)h] (k ^ 1), выделяя при r t особенность:
V(t)= b(t)v2(t) к
v(r )
kh
vMt r]
dr 4 K$fc(t), (21)
¡■kh
k 1 f (j+1)h E E / V(r)exp{ n2n2[t r]} dr.
j=o
' jh
В уравнении (21) в рамках линейной интерполяции заменяем (t [kh, (k + 1)h])
) ад ißi-kh^ [t Т]. (22)
Заметим, что в силу монотонного стремления к нулю v < 0 правая часть меньше левой. то-бы компенсировать завышенные оценки 0(s) и V(t) в (21), при аппроксимации функции Фк (t) будем пользоваться для интегралов формулой правых прямоугольников:
Ek
V(jh)2(t jh), t ^ kh,
j=1
2(o >; ~ e' e > 0.
_ 1 „ n2 n2h
E 1 e_e n2n2Ç
1,3,...,N
n2n2
Функция 2 табулируется предварительно, N 1 (формально N = + ). Влия-
ние предыстории быстро снижается, так что достаточно удерживать несколько последних слагаемых по ] (^ = к, ] = к 1,...). Применение формулы трапеций не приводит к существенному уменьшению погрешности.
В итоге в качестве упрощенной модели в классе ОДУ предлагается
v(t)
1 +
4 к t kh
= b(t)v2(t)
2 t kh
* - Wk 4 кФк(t), k ^ О, (23) 3 п
ык г(кН), £ [кН, (к + 1)Н], г(0) = 1
(Фо()=0, аде = Ь(0)).
Рисунки 8,9 показывают, что параметр Н можно значительно увеличить при явном учете интегральной предыстории функциями Фк(£). На рисунках 10,11 формальное значение N = 0 соответствует отсутствию учета предыстории.
t
t
400 450 500 550 600 650
Рис. 8. Модель (23), влияние Л
Лги)) ■10-14 200 400
300 400
700 800
Рис. 10. Модель (23), влияние N
У(Г(г)) 10"14 400
Рис. 9. Модель (23), влияние Л
10-14 400
та1епа! 12Х18Н10Т
- N = 39
---N = 19
- N = 1
N = 0
Рис. 11. Модель (23), влияние N
вадрати ное глажмвание. Воспользуемся в (21) вместо линейной интерполяции (22) квадратичной: ¿(т) ¿(Ь) + А[Ь т]+ В[Ь т]2. В силу т = Ь ¿(Ь), т = к Л ¿(кЛ) получаем
¿(Ь) ¿(кЛ) г
В[Ь кЛ][Ь т ]+ В[Ь т ]2, (24)
Ь [кЛ, (к + 1)Л], т [кЛ, Ь]. Выбор В < 0 позволит учесть вогнутость функции у (й < 0). Подберем значение параметра В < 0 из условия «гладкой склейки» уравнений (23).
Подстановка вместо (22) квадратичной аппроксимации (24) в (21) приводит к появлению в правой части уравнения (23) дополнительного слагаемого (Ь) 4Вк[Ь кЛ]5/2/[15 п]:
¿(Ь)
1+
4к Ь кЛ
3 п
= 6(ф2(Ь)
2к Ь кЛ
-—=-4кФй (Ь)+ ^ (Ь). (25)
3 п
Сравним теперь к-е уравнение при значении Ь = (к + 1)Л = (к + 1)Л 0 с (к + 1)-м уравнением при Ь = (к + 1)Л = (к + 1)Л + 0:
_ 4к Л
=
2 2к Л
2 ^^ ^+
3 п
+
4Вк
= Л5/2 4кФк ((к + 1)Л),
(26)
15 п
¿((к + 1)Л) = 6((к + 1)Л)^2((к + 1)Л)
4кФй+1((к + 1)Л). (27 Замечаем, что Фк+1((к + 1)Л) =
= Фк ((к + 1)Л) + ¿((к + 1)Л)7(Л),
7 (Л) £
1 е
и2ж2Ь,
П2П2
Здесь для ^(0) полагаем N + . Подчиним выбор Л и В условиям: 3 П7(Л) = Л, 2ВЛ2 = 5^. Тогда уравнения (26) и (27) совпадут. то означает, что при переходе к по-
следующему Н-слою по времени будет непрерывной не только траектория г(£) (по построению), но и производная г(£). Кроме того, полученное значение Вк В отрицательное в силу Wk < 0, т.е. учитывается вогнутость г(£).
Осталось требование 7(Н)/ Н = 1/[3 п]. При Н + имеем
Е
1- e
—n2n2h
П2П2
Е
1
n2n2
Y (h) h
0.
отрезке безразмерного времени £ = £ [кН, (к + 1)Н], эффективно воспроизводит зависимости, представленные в [22], которые получены «точным» решением исходной краевой задачи разностным методом (или системы (12) порядка 10).
Если h +0, то по правилу Лопиталя
lim y(h)/ h = lim 2 e-n"n"h = lim 2 h 1/[4 nhh] = 1/[2 П] > 1/[3 П]. Подходящее значение h > 0 определяется численно. Оценка h < 1 для возможности сглаживания «на стыке шагов» согласуется с приведенными выше физическими соображениями.
Рисунок 12 демонстрирует однозначность выбора h из условия сглаживания Y(h) = h/[3 п]. Рисунки 13 и 14 показывают, что модель (25), требующая лишь численного интегрирования скалярного ОДУ на текущем
0.1 0.2 0.3 0.4
Рис. 12. Модель (25), выбор h
0.5 0.6
jcm-io
J(T(t))-W
400 500 600 700
900 1 000 1 100
Рис. 13. Модель (25), влияние 6q, Еь
Рис. 14. Модель (25), влияние 6q, Еь
Заключение
В статье представлена краевая задача термодесорбции с нелинейными динамическими граничными условиями для моделирования ТДС-спектра дегазации конструкционного материала, предварительно насыщенного водородом. Решение прямой задачи необходимо для качественного сопоставления с экспериментальными данными и как итерационная составляющая алгоритма параметрической идентификации. Вариации параметров позволяют оценить чувствительность ТДС-спектра к изменению выделенных лимитирующих факторов. Для моделирования термо-
десорбционного потока предложен вычислительный алгоритм, требующий (вместо приближенного решения распределенной краевой задачи с текущими приближениями параметров) лишь интегрирования нелинейной системы ОДУ невысокого порядка. Представлены различные модификации с целью понижения порядка системы на основе выделения интегрируемой особенности. Приведены результаты вычислительных экспериментов, использующих экспериментальные данные по никелю и стали 12Х18Н10Т.
Работа выполнена при поддержке РФФИ (грант №15-01-00744).
ИЕ А А
1. и си водорода с металлами / ед. . . а аров. М.: аука, 1987. 296 с.
2. в металла / ед. . ле ельд и . Фелькль. М.: Мир, 1981. Т. 1, 506 с.; т. 2,430 с.
3. и . . нтегральные операторы прогно ирования и иденти икация моделе водородопроницаемости. етро аводск: ар
, 2013. 505 с.
4. и . ., си . . ппрокси-мация краево ада и термодесорбции водорода системо Д // Труды арельского нау ного центра . 2015. 10. C. 42-53. doi: 10.17076/mat147
5. водорода. Фундаментальные и прикладные исследования / ед. . . им-
ук. Саров: Ф - Ф, 2009. 697 с.
6. и . ., и . ., . ., и . . роблемы дега ации металлов.
М.: аука, 1972, 324 с.
7. . ллипти еские ункции. М.: ау-ка, 1984. 312 с.
8. ис . ., . . Ди ерен-циальные уравнения математи еско и ики. М.: МТ , 2002. 368 с.
9. ии . ., ии . . рибли ен-ные методы решения ди еренциальны и интегральны уравнени . М.: аука, 1965, 384 с.
10. ис . ., . ., -
. ., . . роницаемость водорода
ере металлы. М.: МФ, 2008. 144 с.
11. и . ., с . . урс современного анали а, асть 2. М.: л. ред. и .мат. лит., 1963. 516 с.
12. . ., . ., -и . . Методы исследования систем металл-водород. Томск: Т , 2008. 286 с.
13. Evard E. A., Gabis I. E, Yartys V. A. Kinetks of hydrogen evolution from MgH2: experimental studies, melanism and modelling // International Journal of Hydrogen Energy. 2010. Vol. 35. P. 9060-9069. doi: 10.1016/j.ijhydene.2010.05.092
14. Gabis I. E. The method of œnœntration pulses for studying hydrogen transport in solids // Te^nkal Physks. 1999. Vol. 44, no. 1. P. 90-94. doi: 10.1134/1.1259257
15. Handbook of hydrogen storage: new materials for future energy storage / Ed. M. Hirscher. Wiley-VCH, 2010. 353 p.
16. Indeitsev D. A., Semenov B. N. About a model of structure-phase transfomations under hydrogen in uence // Acta Mechanica. 2008. Vol. 195. P. 295-304. doi: 10.1007/s00707-007-0568-z
17. 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
18. Popov V. V., Denisov E. A. Inhibition of hydrogen permeability by TiN: evaluation of kinetic parameters / Hydrogen Materials Science and Chemistry of Carbon Nanomaterials / Eds T. N. Veziroglu et al. Springer. 2007. P. 671-680. doi: 10.1007/978-1-4020-5514-0
19. The hydrogen economy / Eds M. Ball, M. Wietschel. Cambridge Univ. Press, 2009. 646 p.
20. Varin R.A., Czujko T., Wronski Z.S. Nano-materials for solid state hydrogen storage. Springer, New York, 2009. 338 p. doi: 10.1007/9780-387-77712-2
21. Zaika Yu.V., Bormatova E.P. Parametric identi cation of a hydrogen permeability model by delay times and conjugate equations // International Journal of Hydrogen Energy. 2011. Vol. 36, no1. P. 1295-1305. doi: 10.1016/j.ijhydene.2010.07.099
22. Zaika Yu. V., Kostikova E. K. Computer simulation of hydrogen thermodesorption // Advances in Materials Science and Applications. World Acad. Publ. 2014. Vol. 3, iss. 3. P. 120-129. doi: 10.5963/AMSA0303003
23. Zaika Yu. V., Rodchenkova N. I. Boundary-value problem with moving bounds and dynamic boundary conditions: di usion peak of TDS-spectrum of dehydriding // Applied Mathematical Modelling. Elsevier. 2009. Vol. 33, no 10. P. 37763791. doi: 10.1016/j.apm.2008.12.018
24. 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.
Поступила в редакцию 30.05.2016
References
1. Vzaimodeistvie vodoroda s metallami [Interaction of hydrogen with metals]. Ed. A. P. Zakharov. Moscow: Nauka, 1987. 296 p.
2. Vodorod v metallakh [Hydrogen in metals]. Eds G. Alefeld, J. Volkl. Moscow: Mir, 1981. Vol. 1, 506 p.; vol. 2, 430 p.
3. Zaika Yu. V. Integral nye operatory prognozirovaniya i identi katsiya modelei vodorodopronitsaemosti [Integral forecasting
operators and identi cation of hydrogen permeability models]. Petrozavodsk: KarRC of RAS, 2013. 505 p.
4. Zaika Yu. V., Kostikova E. K. Approksimacija kraevoj zadachi termodesorbcii vodoroda sistemoj ODU [Approximation of the boundary-value problem of hydrogen thermal desorption by ODE system]. Trudy Karel skogo nauchnogo centra RAN [Transactions of KarRC of RAS]. 2015. No 10. P. 42-53. doi: 10.17076/mat147
5. Izotopy vodoroda. Fundamental nye i prikladnye issledovaniya [Hydrogen isotopes. Fundamental and applied studies]. Ed. A. A. Yuk-himchuk. Sarov: RFYaTs-VNIIEF, 2009. 697 p.
6. Kunin L. L., Golovin A. I., Surovoi Iu. I., Khokhrin V. M. Problemy degazatsii metallov [Problem of metal degassing]. Moscow: Nauka, 1972. 324 p.
7. Leng S. Elliptic functions. Addison-Wesley publishing, 1973.
8. Martinson L. K., Malov Yu. I. Di erentsial nye uravneniya matematicheskoi ziki [Di erential equations of mathematical physics]. Moscow: MGTU, 2002. 368 p.
9. Mikhlin S. G, Smolitskii Kh. L. Priblizhennye metody resheniia di erentsialnykh i integralnykh uravnenii [Approximate di erential and integral equation solving]. Moscow: Nauka, 1965. 384 p.
10. 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.
11. Whittaker E.T., Watson G.N. A Course of Modern Analysis. Cambridge University Press, 1996, 612 p.
12. 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.
13. 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. 9060-9069. doi: 10.1016/j.ijhydene.2010.05.092
14. Gabis I. E. The method of concentration pulses for studying hydrogen transport in solids. Technical Physics. 1999. Vol. 44, no. 1. P. 90-94. doi: 10.1134/1.1259257
15. Handbook of hydrogen storage: new materials for future energy storage. Ed. M. Hirscher. Wiley-VCH, 2010. 353 p.
Заика Юрий Васильевич
руководитель лаб., д. ф.-м. н. Институт прикладных математических исследований Карельского научного центра РАН ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: zaika@krc.karelia.ru тел.: (8142) 766312
Костикова Екатерина Константиновна
научный сотрудник, к. ф.-м. н. Институт прикладных математических исследований Карельского научного центра РАН ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: kostikova@krc.karelia.ru тел.: (8142) 766312
16. Indeitsev D. A., Semenov B. N. About a model of structure-phase transfomations under hydrogen in uence. Acta Mechanica. 2008. Vol. 195. P. 295304. doi: 10.1007/s00707-007-0568-z
17. 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.
18. Popov V. V., Denisiv E. A. Inhibition of hydrogen permeability by TiN: evaluation of kinetic parameters. Hydrogen Materials Science and Chemistry of Carbon Nanomaterials Eds T. N. Veziroglu et al. Springer. 2007. P. 671-680. doi: 10.1007/978-1-4020-5514-0
19. The hydrogen economy. Eds M. Ball, M. Wietschel. Cambridge Univ. Press, 2009. 646 p.
20. Varin R.A., Czujko T., Wronski Z.S. Nano-materials for solid state hydrogen storage. Springer, New York, 2009. 338 p. doi: 10.1007/9780-387-77712-2
21. Zaika Yu. V., Bormatova E.P. Parametric identi cation of a hydrogen permeability model by delay times and conjugate equations. International Journal of Hydrogen Energy. 2011. Vol. 36, no. 1. P. 1295-1305. doi: 10.1016/j.ijhydene.2010.07.099
22. Zaika Yu. V., Kostikova E. K. Computer simulation of hydrogen thermodesorption. Advances in Materials Science and Applications. World Acad. Publ. 2014. Vol, 3, iss. 3. P. 120-129. doi: 10.5963/AMSA0303003
23. Zaika Yu. V., Rodchenkova N. I. Boundary-value problem with moving bounds and dynamic boundary conditions: di usion peak of TDS-spectrum of dehydriding. Applied Mathematical Modelling. Elsevier. 2009. Vol. 33, no. 10. P. 37763791. doi: 10.1016/j.apm.2008.12.018
24. 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.
Received May 30, 2016
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
Kostikova, Ekaterina
Institute of Applied Mathematical Research, Karelian Research Centre, Russian Academy of Sciences 11 Pushkinskaya St., 185910 Petrozavodsk, Karelia, Russia
e-mail: kostikova@krc.karelia.ru tel.: (8142) 766312
