CC BY
22
Труды Карельского научного центра РАН № 5. 2011. С. 97-102
УДК 519.6:539.2
РАЗНОСТНАЯ СХЕМА ДЛЯ КРАЕВОЙ ЗАДАЧИ ВОДОРОДОПРОНИЦАЕМОСТИ ПРИ НАЛИЧИИ ДЕФЕКТА ЗАЩИТНОГО ПОКРЫТИЯ
Н. И. Родченкова, Е. К. Костикова
Институт прикладных математических исследовании Карельского научного центра РАН
Рассматривается водородопроницаемость цилиндрической перегородки при наличии на входной поверхности дефекта защитного покрытия. Для различных конструкционных материалов, когда лимитирующими являются процессы диффузии и десорбции, представлены соответствующие математические модели в форме краевых задач с условиями I рода и нелинейными динамическими граничными условиями. Разработана разностная схема для численного моделирования проникающего потока.
Ключевые слова: водородопроницаемость, нелинейные краевые задачи, разностные схемы.
N. I. Rodchenkova, E. K. Kostikova. DIFFERENCE SCHEME FOR THE BOUNDARY-VALUE PROBLEM OF HYDROGEN PERMEABILITY IN THE PRESENCE OF PROTECTIVE COATING DEFECT
We consider the permeability of cylindrical barrier to hydrogen in the presence of protective coating defect on the inlet surface. For different structural materials, when diffusion and desorption are the limiting processes, corresponding mathematical models in the form of boundary-value problems with first-type conditions and nonlinear dynamic boundary conditions are presented. The difference scheme for numerical modelling of the permeation flux has been developed.
Key words: hydrogen permeability, nonlinear boundary-value problems, difference schemes.
Постановка задачи
Снижение проникновения водорода и его изотопов сквозь стенки из конструкционных материалов является важнейшей задачей при решении комплексных проблем хранения и транспортировки водорода, защиты от водородного охрупчивания, контроля содержания трития в защитных системах будущих термоядерных реакторов (проект ITER). Конструкция из металла или сплава обеспечивает
необходимую механическую прочность перегородки, а нанесенное защитное покрытие должно препятствовать миграции изотопов водорода. Дефекты защитной пленки могут подвергать соответствующую область конструкционного материала прямому воздействию водорода [Писарев и др., 2008]. В статье [Zajec, 2011] поставлена задача математического моделирования водородопроницаемости цилиндрического образца радиуса Ь и высоты
■©
H в случае, когда диффузия является единственным лимитирующим процессом. На входной поверхности z = 0, покрытой тонкой защитной пленкой, присутствует дефект малого радиуса r0 (булавочное отверстие), через который проникает водород. Остальная часть входной поверхности водородонепроницаема, как и боковая поверхность. На выходной стороне z = H поддерживается вакуум. В начальный момент времени t = 0 образец обезводоро-жен. Затем на входной стороне скачкообразно повышается давление молекулярного водорода до уровня р. Если пренебречь относительно быстрым (это зависит от р, материала и размеров образца) переходным процессом, то можно считать, что концентрация растворенного водорода под дефектом поддерживается на постоянном уровне Со (находится в равновесии с газообразной фазой по закону Сиверт-са, с0 а /р). Растворенный (атомарный) водород диффундирует к выходной поверхности. С помощью масс-спектрометра регистрируется проникающий поток.
Аналитический анализ краевой задачи без учета поверхностных процессов проведен лишь для случая полупространства (r ^ +го) [Warrick et al., 1992; Rajendran, Sangaranarayanan, 1995]. Основным недостатком такой постановки задачи является то, что поверхностные процессы, которым в последнее время уделяется повышенное внимание, в модели не учитываются.
Цель работы — построить разностную схему для моделей водородопроницаемости цилиндрического образца заданных размеров при наличии дефекта как без учета, так и с учетом влияния поверхности.
Диффузионная модель
Рассмотрим краевую задачу водородо-проницаемости цилиндрической перегородки только с учетом диффузии в объеме:
І = D
д2с І дс д2 с дг2 + r дr + дz2
r Є (0, L), z Є (0, H), t Є (0, t*),
c(t,r,0) = Со, r Є [0,го], го < L, (2)
дС
— (t,r, 0) = 0, r Є (го, L], t ^ 0, (3)
дz
c(t, r, H) = 0, r Є [0, L], t ^ 0,
(4)
Iе (^£,г) = 0, Iе ^ +М=0, (5)
с(0,г, г) = 0, г е [0, £], г е [0,Я].
Здесь с(£, г, г) - концентрация атомарного водорода в конструкционном материале (метал-
ле или сплаве); Б - коэффициент диффузии. Момент времени £* определяется выходом проникающего потока на стационарное значение. Следует отметить, что установление носит асимптотический характер. Но £* не следует выбирать слишком большим, чтобы переходные процессы «не потерялись» на фоне стационара. Условие сг (£, +0, г) следует из симметрии распределения с(£, г, г).
Уточнение поставки задачи. Цель состоит в разработке разностной схемы для моделирования потока водорода с выходной поверхности цилиндрического образца:
2nr dr.
z=H
Изменения выражения .] (£) с учетом рекомбинации атомов водорода в молекулы на поверхности (в приповерхностном объеме) приведены в последующих разделах статьи.
Численное моделирование позволяет выявить особенности кинетики проникновения водорода и оценить влияние параметров модели, включая размеры образца и дефекта.
Разностная аппроксимация краевой задачи
Будем следовать методике, разработанной для задачи термодесорбции [Родченкова, Заика, 2010]. Следуя стандартной технике [Самарский, 1971], введем пространственную сетку
ri = i hr, i =0, І, zj = jhz, j = 0, І,
,N1 =[L/hr ]
,N2 = [H/hz]
и сетку по времени мт={ = кт, к = 0,1,..., К = [і*/т]}.
Обозначим через сй- приближенные значения объемной концентрации с(І&, Гг,2-). При переходе с к на к + 1 слой по времени будем использовать следующие обозначения, опуская номер слоя по времени: сг,- = С-, сг,- =
С^+1/2, Сг,і = С—\ где (Гг,2-) Є ^^ Є .
i,j
Для уравнения (1) рассмотрим неявную разностную схему метода переменных направлений, называемую продольно-поперечной (схемой Писмена-Рэкфорда) [Косарев, 2000]. Переход от слоя к к слою к + 1 осуществляется в два этапа. На первом этапе определяются промежуточные значения с^,- из соотношений
Ci,j с
i,j
0,5 tD
ci,j—1 2 ci,j + ci,j+1
h2
+ ^i—1,j 2 Ci,j + ^i+1,j + ^i+1,j ^i—1,j
h2
2 rihr
На втором этапе, пользуясь найденными Сг,-, находим Сг,- из соотношений
Сг,- Сг,- Сг,-—1 2 Сг,- + Сг,-+1
+
0, 5 ті Н2
Сг—1,- 2 Сг,- + Сг+1,- Сг+1,- Сг—1,-
Н2
+
2 ГгНг
. (6)
Данные соотношения рассматриваются во внутренних узлах сетки (г = 2,..., N1 — 1, 3 = 1,..., N2 — 1). Суммарная погрешность аппроксимации 0(т2 + ^2 + Л^).
Прогонка по радиусу г
Рассмотрим переход с к-го на (к + 1/2)-й слой. В обозначениях
А =
1 — [2І]
1
К
В =
1 + [2і]
1
к = 2
К
Нг \2 с г,- — 1 2сг,- + сг,-+1
+
К
Н2 1+і? 1-2
К
"г,-
при каждом фиксированном 3 = 1, 2,..., N2—1 получаем при к ^ 0
АгСг—1,- Сг,- + ВгСг+1,- + ^г,- — °.
Значения в начальный момент времени (на нулевом слое) известны: с0,- = 0. Следуя методу прогонки, ищем приближенные значения концентрации в узлах сетки на (к + 1/2)-м слое (к ^ 0) по времени в виде
По формуле (7) определяем искомые значения Сг,-, г = 0,..., N1 — 1, 3 = 1,..., N2 — 1.
Теперь найдем оставшиеся приближения сг,- при 3 = 0 и 3 = ^, г = 0,1,..., ^. Значения = 0 согласно граничному условию (4) (г = Н). Обозначим го = тах{г : гг ^ го}. Тогда из граничных условий (2), (3) (г = 0) получаем соотношения сг,0 = с0, г = 0,1,...,г0 и Сг,0 = (4Сг,1 — Сг,2)/3, г = го + 1, . . . , N1.
Прогонка по переменной г
Поскольку в цилиндрических координатах возникает особенность при г ^ +0, и на границе г = 0, задано смешанное краевое условие, то переход с (к + 1/2)-го слоя на (к + 1)-й совершается в два шага.
Первый шаг: г = 1, г ^ +0 (аппроксимируем сг/г к сгг), реализуется алгоритм прогонки для уравнения с* = Б(2сгг + с^).
Второй шаг: г = 2,..., N — 1, г > 0, прогонка для уравнения с* = Б(сгг + сг/г + с^). Аппроксимация уравнения с* = Б(2сгг + с,^):
С1,- с1,- С1— 1 2 С1- + С1,—+1
0,5 ті
+2
Н2
С 0,- — 2 С 1,- + С 2,
Н2
В обозначениях О = 2(1 + Н2/[1т]), Нл2/
— / Н \2
Сг,- = аг+1,- Сг+1,- + вг+1,-, І = 0,...,^ — 1.(7) І?1,І = 2((С0- — 2С 1- + ^ +(О — 2)С 1
Прогоночные коэффициенты: І = 2, 3,..., N1, получаем
Вг_
аг,-:
г1
1 Аг—1 аг—1,-
А,- = '
^г—1,- + Аг—1 Дг— 1,1 Аг— 1аг— 1,-
При г ^ +0 имеем
сг/г = (сг(£, г, г) — сг(£, 0, г))/г к сгг.
Начальные коэффициенты а1,-, ,01,- находим из аппроксимации уравнения с* = Б(2сгг+с^) на (к + 1/2)-м слое, г = 1, и условия сг|+0 = 0: “1,- = (6 — к)/4, в1,- = ^1,- к/4.
Ближайшая цель — найти значение с щ -, необходимое для реализации обратного хода прогонки. Запишем аппроксимацию первого из граничных условий (5) (г = Ь). В граничном узле с точностью до 0(^2) имеем:
2^гсг(^+1/2, Ь, ) ~см1-2,- — 4 + 3с^1,-.
Подставляя значения с щ_2,-, с щ _1,- из соотношения (7), имеем
С1,-—1 — О С1,- + С1-+1 + ^1,- — 0, к ^ 0. (8)
Ищем приближенные значения концентрации на (к + 1)-м слое по времени в виде: к ^ 0,
С1,- = «1,-+1 С1-+1 + в1,-+1, 3 = 0,1, ...,N2 — 1. (9) Прогоночные коэффициенты: 3 = 2, 3,..., N2,
а1,-
1
О — «11
, А,- =
^1,-—1 + в1,-—1 О — «11
Начальные коэффициенты «1,1, ,01,1 находим из формулы (9) при 3 = 0 и условия (2) (С1,0 = с0): а1,1 = ° ^1,1 = с0.
Для разностного соотношения (6) в обозначениях
1 2іі Сг—1,-
в^!, - (4 а^1-1, - ) в^1-1,І 3 + а^1,і (а^1-1,І — 4)
21
—іт)+(1+2і) Сг+1-
при каждом г = 2, 3,..., N1 — 1 получаем
Сг,-—1 ОСг,- + Сг,-+1 + -^г,- — 0, к ^ 0. (10)
Ищем значения концентрации в узлах сетки на (к + 1)-м слое по времени в виде: к ^ 0,
Сг,- — аг,-+1 Сг,-+1 + вг,-+1, 3 —0, 1,...,^2 1. (11)
Прогоночные коэффициенты: 3 = 2, 3,..., N2,
Прогоночные коэффициенты:
1
О «г,?—1
вг,- =
^г,;— + Дг,- —1
О аг,-—1
д с
из»—Ь с2(1 г, 0) = — і — дг
д с
Ьс 2(і,г,Я ) = — і— дг
^=0
-=я
г є [0, Го], (12) г є [0,£], (13)
■](£) = [ Ь с2(£, г, Н) 2пг ^г.
0
Здесь Ь — коэффициент объемной десорбции (эффективной рекомбинации [Писарев и др., 2008]), ^ — кинетическая константа, р — давление молекулярного водорода, в — коэффициент прилипания водорода к поверхности. При необходимости можно учесть различие входной и выходной поверхностей: Ь = Ь1 при г = 0 и Ь = Ь2 при г = Н.
1. Метод встречных прогонок. Рассмотрим переход с (к + 1/2)-го слоя на (к + 1)-й. Для соотношений (8), (10), используя метод встречных прогонок [Самарский, Николаев, 1978], будем искать приближенные значения концентрации в следующем виде (г = 1,..., г0, к ^ 0):
Сг,- — аг,-+1 Сг,-+1 + вг,-+1 + 0г,-+1 Сг^
3 =0,1,...,^ — 1;
Сг,- = ССг,-—1 Сг,-—1 + /5г,-—1 + Ог,-—1 Сг,М2 ,
3 = 1, 2,...,N2.
с^г,-—1 = аг,^2-(3--1)>
Начальные коэффициенты при г ^ г0 (г ^ г0) находим из (11) при 3 = 0 и условия (2): аг,1 = 0, вг,1 = с0. Начальные коэффициенты при г > г0 (г > г0) определяем из (10) при 3 = 1 и условия (3): аг,1 = 2—О/2, Дд = —д/2.
Граничные значения, необходимые для реализации обратного хода метода прогонки, находим из условия (4): Сг,щ = 0. По формулам (9), (11) вычисляем Сг,-, г = 1,.. .,N1 — 1, 3 = 0,... ,N2 — 1.
Найдем оставшиеся приближения Сг,- при г = 0 и г = N1, 3 = 0,1,..., N2. Используя второе из граничных условий (5) (сг|+0 = 0), получаем С0,- = (4с1,- — с2,-)/3. Согласно первому условию (5) Сщ1,- = (4СЩ1_1,- — Сщ1_2,-)/3. Модификация модели с учетом объемной десорбции
Вместо краевых условий (2), (4) используем
аг,-+1—(О — аг,- ) \
вг,-+1=(вг,-+-^г,-)аг,-+1, Дг,-— 1=(/Сг,-+-^г,-),С:г,-— 1,
7г,-+1=7г,-аг,-+1, 7г,-—1 = 7г, щ2_0'_1).
Отметим, что выражения для величин -1,- и -г,-, г > 1, различны. Начальные прогоночные коэффициенты: аг,1 = вг,1 = /5г, щ _1 = 0, 7г,1 = 1. 2
В обозначениях А = аг,щ _1 — 4,
А = 3 + аг,щ2 А, В = вг,щ2_1 + вг,щ2 А,
В = Дд + /?г,0, ^ = 7г,^2_1 + 7г,щ А
получаем аппроксимацию граничных условий (12), (13) с точностью до О(^) в виде
'2Н- , л2
— Ь с2,0 + Асг,о +
2Н-
Ь с2^2+ Асг,^2+^Сг,о + В] = 0.
Система уравнений имеет единственное решение Сг,0 > 0, Сг,щ > 0, по крайней мере, при малых Л,г ~ ^ (сравнимых по порядку).
Осталось получить значения концентрации Сг,- при г > г0. Изменения коснутся лишь границы г = Н .В объеме выполняется соотношение (10), прогоночные коэффициенты остаются без изменений. Модификация появляется после выполнения прямого хода прогонки. Для нахождения Сг,щ аппроксимируем (13) с точностью до 0(^2) и используем (11):
2^ Ь^—1 С2,щ+АСг,щ2+ В = 0,
где г > г0, А = 3 + аг,щ А, В = вг,щ2_ + вг,щ2 А, А = аг,щ_1 — 4. При малых Л,г ~ ^ кор)ни квадратного уравнения имеют разные знаки, по физическому смыслу задачи выбираем положительный.
При переходе с к-го на (к + 1/2)-й слой оставшиеся граничные значения сг,0, сг,щ для г = 0,..., N1 определяются как положительные корни квадратных уравнений, полученных после подстановки в (12), (13) выражений
С-(£,Гг, 0)
—сг,2 + 4 Сг,1 — 3с
сг,0
2Н-
(14)
- /. ТТ\ Сг,^2-2 4 Сг,^2-1 + 3сг,^2
С-(‘-Гг-Я> « -------------------2Н-----------
где Сг,1, Сг,2, Сг,^2-1, Сг,^2- уже известны по результатам прогонки по Г.
«г,7 =
0
(ш)
2. Итерационный метод. На (к + 1)-м слое по времени аппроксимируем с^ (^^+1, г, 0) ~ [—3сг,о + 4йг,1 — сг,2]/2Л,.г. Подставляя в граничное условие (12) при г = 0, находим Сг,0 = /о(Сг,1, Сг,2) как положительный корень квадратного уравнения. Аналогичным образом определяем с^2 = /^2 (<Ч^2-1, <4^-2) как положительный корень из (13) при г = Н. Значения Сг,1, Сг,2, с*,^-, предваритель-
но подсчитываются по явной разностной схеме, аппроксимирующей уравнение диффузии (1). С текущими сг,о, сг,^2 решаем методом прогонки по переменной г трехдиагональную систему линейных алгебраических уравнений и находим новые приближения концентраций Сг,1, Сг,2, С*,^-, — (и остальные значения
сг,- для ] = 3,..., N — 3, г = 0,..., го). Снова решаем квадратные уравнения относительно Сг,о, Сг,^ и повторяем вычисления до установления граничных значений (обычно 2-3 итерации). Нахождение приближений сг,- при г > го описано в предыдущем пункте. Затем переходим к следующему слою по времени.
Модификация модели с учетом поверхностной десорбции
Вместо условий (2), (4) используем
^ = ^р — Ьдо2(£,г) + , г е [0, го], (15)
д£ 0 дг о
ilr = (*,r) - D|
r G [0,L], (16)
Ci,0 ci,0
gT
=o
^sp - b( + Ddz Cj,o
C*,0\ 2
+(1 - o) ^sp -+ Ddz Ci,o
Ci,N2 Ci,N2
gT
=o
- b(^)2 - DdzC*„
+(1 - o)
— b
- Ddz Ci
N2
Здесь дг Сг,о, дг сг,^ определяются выражениями, аналогичными (14). В дальнейшем полагаем а = 1/2. Как и в модификации модели для объемной десорбции (сохраняя обозначения), нахождение сг,о, сг,^2 сводится к решению системы квадратных уравнений:
' Ь л2 /ДА 2\
~2 Сг о + ( ТГй 1--) Сг,о + Сг,^2+ С1 = 0
g2 ",0 \2hz gT
b л0 /DA 2' __л
72 C2 n + ^Т“ + ~ ^i,N2+ 777“1Ci,0 + C2 = 0,
2hz
DG
I g2 i,n2 \^2hz gV i,N^ 2hz
DB b 2 2
C1 = 2^ + g2 C2,0 - gTCi,0 - Ddzci,0 - 2^sp,
2hz g gT
b
2
С2 = ХТ I 2 с2^-------с»,^2+ Дд-гсг,^2 .
2яг g2 ь^2 gr 2 2
Все значения сг,- на к-м слое известны.
2. Итерационный метод. Для краевых
условий (15), (16) используется схема с весами
(при £ = ^^+1/2 + т/4 порядок аппроксимации
0(т2 + Ь2)):
Ci,0 Ci,0
0,5 Tg
C*,N2 C*,N2
0,5 Tg
c(t, r, 0) = gq0 (t, r), c(t, r, H) = gqH (t, r),
J(t) = [ bq^ (t, r) 2nrdr.
J0
Здесь q0, qH — поверхностные концентрации на входной и выходной поверхности, g = const — коэффициент соответствия концентраций атомов водорода в объеме и на поверхности (коэффициент быстрого растворения).
1. Метод встречных прогонок. Остановимся на вычислении граничных концентраций Сг,о, Ci,N , i ^ i0. Чтобы сохранить порядок аппроксимации 0(т2 + h^), для краевых условий (15), (16), будем использовать схему с весами:
=0,5
+0,5
=0,5
+0,5
^sp - Ь(—°) + Ddz Ci,0 ^sp - + Ddz Cj,0
- b( g ^ - Ddz Ct,N2
- MCi^l2 - Ddz Ci.
g
N2
Граничные значения сг,о, сг,^ (г ^ го) определяются как положительные корни квадратных уравнений:
b
/ 3D 4
2 C2,0 + ( ^7Т“ + Ci,0 + C1 = 0,
g2 . \2hz gT
b 2 /3D 4 ,
g2 C2,N2+(v2h; + gr)Ci,N2+ C2=0
D / \ b _2
Cl = 2h;(Ci,2 - 4Ci,^ + g2^0 4
-----Ci,0 - Ddz Ci,0 - 2^sp,
gT
D
C2 = 2h ^*.N2-2 - 4Ci,N2-l)
b 2 4
+ g2 C^,N2 + DdZCi,N2 - gTCi,N2 .
g
В выражениях, аналогичных (14), значения Сг,{1,2}, Сг, №-1,^2-2} предварительно подсчитываются по явной разностной схеме, аппроксимирующей уравнение диффузии (1). С текущими Сг,о, Сг,^ решаем методом прогонки трехдиагональную систему линейных алгебраических уравнений и находим новые приближения Сг,{1,2} , Сг, №-1^2-2} (и остальные значения Сг,-, ,] = 3,..., N2 — 3, г = 0,..., го). Снова решаем квадратные уравнения относительно Сг,о, Сг,^ и повторяем вычисления до установления граничных значений (обычно 2-3 итерации). Для нахождения граничной концентрации Сг,^ при г > го аппроксимируем (16) С точностью до 0(т2 + ) и используем про-
гоночные коэффициенты, найденные прямой прогонкой при использовании условия (3):
%-2 СЦ+ [Д(2^)-1 А + 2(^)-1]Сг,^+ ® = 0
где А = 3 + аг,^2 А, А = «г,^2-1 — 4, ® = —2сг,^2/(^ + д(вг,^2-1 + А,*2А)/(2Л-.г). При
малыше Л,,г, т корни кв£1дратного уравнения разных знаков, по физическому смыслу задачи выбираем положительный. Вычисленные аппроксимации сг(£, г, Н) в модели (1)—(5); концентрация с(£, г, Н) в случае граничных условий (12), (13); поверхностная концентрация (£, г) для модификации модели (15), (16) дают возможность приближенно вычислить .] (£).
Заключение
Таким образом, представленный вычислительный алгоритм позволяет по входным данным краевых задач водородопроницаемости
СВЕДЕНИЯ ОБ АВТОРАХ:
Родченкова Наталья Ивановна
научный сотрудник, к. ф.-м. н.
Институт прикладных математических исследований КарНЦ РАН
ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: [email protected] тел.: (8142) 766312
Костикова Екатерина Константиновна
младший научный сотрудник
Институт прикладных математических исследований КарНЦ РАН
ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: [email protected] тел.: (8142) 766312
(с учетом различных лимитирующих факторов) численно моделировать проникающий поток водорода сквозь перегородку из конструкционного материала при наличии дефекта защитного покрытия.
Работа выполнена при поддержке РФФИ (грант № 09-01-00439).
Литература
Косарев В. И. 12 лекций по вычислительной математике. М.: МФТИ, 2000. 224 с.
Писарев А. А., Цветков И. В., Маренков Е. Д., Ярко С. С. Проницаемость водорода через металлы. М.: МИФИ, 2008. 144 с.
Родченкова Н. И., Заика Ю. В. Численное моделирование десорбции водорода с цилиндрической поверхности // Труды Карельского научного центра РАН. 2010. № 3. С. 72-82.
Самарский А. А. Введение в теорию разностных схем. М.: Наука, 1971. 553 с.
Самарский А. А., Николаев Е. С. Методы решения сеточных уравнений. М.: Наука, 1978.592 с.
Rajendran L., Sangaranarayanan M. V. A two-point Pade approximation for the non-steady-state chronoamperometric current at ultramicrodisc electrodes // Journal of Electroanalytical Chemistry. Elsevier, 1995. Vol. 392. P. 75-78.
Warrick A. W., Broadbridge P., Lomen D. O. Approximations for diffusion from a disc source // Applied Mathematical Modelling. Elsevier, 1992. Vol. 16. P. 155-161.
Zajec B. Hydrogen permeation barrier - recognition of defective barrier film from transient permeation rate // International Journal of Hydrogen Energy. Elsevier, 2011. Vol. 36. P. 7353-7361.
Rodchenkova, Natalia
Institute of Applied Mathematical Research, Karelian Research Centre, Russian Academy of Science 11 Pushkinskaya St., 185910 Petrozavodsk, Karelia, Russia
e-mail: [email protected] tel.: (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: [email protected] tel.: (8142) 766312
