Моделирование водородопроницаемости сквозь дефект защитного покрытия Текст научной статьи по специальности «Математика»

Труды Карельского научного центра РАН № 5. 2012. С. 22-32

УДК 519.6:539.2

MОДЕЛИРОВАНИЕ ВОДОРОДОПРОНИЦАЕМОСТИ СКВОЗЬ ДЕФЕКТ ЗАЩИТНОГО ПОКРЫТИЯ

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

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

Рассматривается водородопроницаемость цилиндрической перегородки при наличии дефекта защитного покрытия на входной поверхности конструкционного материала. Лимитирующим фактором является диффузия. На основе неявной разностной схемы разработан вычислительный алгоритм решения соответствующей краевой задачи, приведены результаты численного моделирования проникающего потока.

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

Yu. V. Zaika, E. K. Kostikova. MODELING OF HYDROGEN PERMEABILITY THROUGH PROTECTIVE COATING DEFECT

We consider the hydrogen permeability of a cylindrical barrier with a protective coating defect at the inlet surface of the structural material. Diffusion is the limiting factor. The computing algorithm based on implicit difference scheme was developed for solving the corresponding boundary-value problem. Results of numerical modeling of the permeability flux are presented.

Key words: hydrogen permeability, boundary-value problems, difference schemes, numerical modelling.

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

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

прочность, а нанесенное защитное покрытие должно препятствовать миграции изотопов водорода. Дефекты защитной пленки могут подвергать соответствующую область конструкционного материала прямому воздействию водорода. Исследование влияния на водородопроницаемость геометрических факторов на поверхности (шероховатости, трещин) и в объеме (дефектов кристаллической структуры, пор) представлено в [5]. В статье [12] анализируется модель водородопрони-цаемости цилиндрического образца (радиуса

0

основания Ь и высоты Н), когда диффузия является единственным лимитирующим процессом. На входной поверхности г = 0, покрытой тонкой защитной пленкой, присутствует дефект покрытия малого радиуса г0 («булавочное отверстие»), через который проникает водород. Остальная часть входной поверхности водородонепроницаема, как и боковая поверхность. На выходной стороне г = Н поддерживается вакуум. В начальный момент времени £ = 0 образец обезводорожен. Затем на входной стороне скачкообразно повышается давление молекулярного водорода до уровня р. Если пренебречь относительно быстрым (это зависит от р, материала и размеров образца) переходным процессом, то можно считать, что концентрация растворенного водорода под дефектом поддерживается на постоянном уровне с (находится в равновесии с газообразной фазой по закону Сивертса, с а /р). Растворенный (атомарный) водород диффундирует к выходной поверхности, рекомбинирует в молекулы и десорбируется. С помощью масс-спектрометра регистрируется проникающий поток.

Аналитический анализ соответствующей краевой задачи проведен лишь для случая «бесконечной пластины» (Ь ^ +го) [10, 11]. Основным недостатком является то, что динамика поверхностных процессов, которым в последнее время уделяется повышенное внимание, в модели не учитывается.

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

ДИФФУЗИОННАЯ МОДЕЛЬ

Рассмотрим краевую задачу водородопро-ницаемости цилиндрической перегородки с дефектом в центре круга защитного покрытия входной поверхности, когда лимитирующим

фактором является диффузия [10-12]:

dc / d2c 1 dc d2c \

— = D---------+ - .----1----1

dt

dr2 + r dr + dz2 )

r G (0, L), z G (0,H), t G (0,t*), c(t, r, 0) = Co, r G [0, ro], ro < L,

dc

dz(t, r, 0) = 0, r G (ro, L], t ^ 0, c(t, r, H)=0, r G [0, L], t ^ 0, § (t,L,z)=0, ^ (t +0,z) = 0,

(2)

(3)

(4)

(5)

c(0, r, z) = 0, r G [0, L], z G [0, H]. (6)

Здесь c(t, r, z) — концентрация атомарного водорода в конструкционном материале (металле или сплаве); c0 = const (для определенности считаем c0 = с); D — коэффициент диффузии. Предполагается, что эксперимент проводится при постоянной температуре T = T, материал практически однороден, так что co, D — константы. Линейные размеры дефекта относительно малы, считаем его круговым. Момент времени t* определяется выходом проникающего потока на стационарное значение. Установление носит асимптотический характер: J(t) « const, t ^ t*. Но t* не следует выбирать слишком большим, чтобы переходные процессы «не потерялись» на фоне стационара. Граничное условие (2) соответствует быстрому (в масштабе t*) достижению локального равновесия. Условие cr (t, +0, z) = 0 следует из симметрии распределения c(t,r, z).

Уточнение постановки задачи. В рамках модели (1)-(6) задача состоит в разработке алгоритма численного моделирования потока водорода с выходной поверхности:

fL dc J(t) = -Dl &

2nr dr.

z=H

Замечание 1. Формально нулевые начальные данные (6) не согласованы с начальным условием (2) при £ = +0 (мгновенный скачок концентрации). Решение краевой задачи следует понимать в рамках теории обобщенных решений. При последующей дискретизации модели на нескольких первых шагах по времени учтен реальный (пусть и относительно быстрый) переходный процесс.

23

Разностная аппроксимация

По стандартной технике [9], введем сетки — кт, к — 0,1,..., К, т — ¿*/К}.

г, — іЛ,г, і —0,...,Жі, Л,г — Ь/Жі 1 ¿7 — і^, і — 0,...,^2, Л.^— Я/^г/

Обозначим через с^- приближенные значения

объемной концентрации с(£&, г,,^-). С целью некоторого технического упрощения изложения при переходе с к-го на (к + 1)-й слой по времени і будем опускать номер слоя:

сі,7 — Сі,7 ,

с.. — *+1/2

с,7 — сі,7 ,

А. . — сй+1 аі7 — сі7 .

Для уравнения (1) рассмотрим неявную разностную схему метода переменных направлений, называемую продольно-поперечной (схемой Писмена-Рэкфорда) [2]. Переход от слоя к к слою (к + 1) осуществляется в два этапа. На первом этапе определяются промежуточные значения Сі,і из соотношений

сі,7 — сі,7 сі,7—1 — 2 сі,7 + сі,7+1

0,5 т^

Л2

+

+ сі—1,7 2 сі,7 + сі+1,7 + сі+1,7 Сі—1,7

^2

2 г,;Л,г

На втором этапе, пользуясь уже известными Сі,і , находим Сі,і из системы уравнений

аі,7 сі,7

0, 5 т^

Аі,І — 1 2 Аі,7 + А

Аі,7+1

Л2

+

(7)

+ сі—1,7 2 сі,7 + Сі+1,7 + Сі+1,7 Сі—1,7

2 пЛг

Данные выражения (для первого и второго этапа) рассматриваются во внутренних узлах (при г = 2,..., N1 — 1, ^ = 1,..., N2 — 1). Погрешность аппроксимации есть 0(т2 + Л>2 + Л^).

Прогонка по радиусу г

Рассмотрим переход с к-го на (к + 1/2)-й слой. В обозначениях к = 2(1 + Л2/[^т]),

А — І-М—1, в, — 1 + |2і1—1

К

К

—(.г

\2 сі,7— 1 2сі,7 + сі,і + 1

+

К

-іі

при каждом фиксированном і — 1, 2,..., N2—1 получаем при к ^ 0

А,с,—1 ,7 Сі,і+ВіСі+1,і+^і,і — 0, (0<і<^). (8)

Значения в начальный момент времени известны: с0і — 0. Следуя методу прогонки, ищем приближенные значения концентрации в узлах сетки на (к + 1/2)-м слое (к ^ 0) в виде

сі,7 — аі+1,7 сі+1,7 + ві+1,і, і — 0,...,Ж1 1. (9)

Прогоночные коэффициенты: і — 2, 3,..., N1,

а

і- 1

і,7

в,7 —

1 — Аі—1 аі— 1,7 ’ Я—1,7 + Аі—1 ^і—1,

1 АІ — 1аІ— 1,7

(10)

Варианты начальных коэффициентов

• Запишем соотношение (8) для і — 1:

13

2К с0,7 — с1,7 + 2К С2,7 + ^1,7 — 0. (11)

Используя второе из граничных условий (5) (сг|+о — 0), с точностью до 0(^2)

— С2,7 + 4с1,7 — 3со,7 ~

~ 2^гсг(tfc+1/2, +0, г7) — °. (12)

Учитывая (11) и (12), находим а1,7 — (6 — к)/4, ^1,7 — ^1,7 к/4. (13)

• При г ^ +0 имеем

сг/г — (сг(¿, г, г) — сг(¿, +0, ¿))/г « сгг.

Начальные прогоночные коэффициенты 01,7, ,01,7 находим из аппроксимации уравнения с* — ^(2сгг+) на (к+1/2)-м слое при і — 1 и условия сг |+о — 0:

2(с0,7+с2,7)/К—с1,7(1+2/к)+-^1,7 —0, (14) с2,7 — 4с1,7 — 3с0,7, откуда “1,7 — (6 — к)/4, в1,7 — ^1,7 к/4. (15) Получили те же выражения, что и (13).

• Для нахождения начальных прогоноч-ных коэффициентов воспользуемся осевой симметрией. Значению г — 0 соответствует ось цилиндра. Формально пространственную область [0, Ь] х [0, Я] можно рассматривать как половину осевого сечения цилиндра. Ориентируясь на общую идею метода фиктивных областей [8], формально добавим симметричный узел (і — —1) и положим с—1,7 — 61,7,

С—1,7 ~ c(tfc+1/2, , г7 ) — c(tfc+1/2, , г7).

г

Продолжая действовать формально, запишем уравнение (8) для г = 0 (г « 0), учитывая условие симметрии с-1,^ = 61,7:

^ (1 — I + 17 (1 + | )+«* = 0

Особенность (деление на г « 0) при условии с-1,7 = С1,7 «сокращается», получаем «1,7 = 2/к, ,01,7 = ^о,7. Отметим, что при Лг/^т] ^ 1 начальные прогоночные коэффициенты такие же, как и в (13), (15). Сравнимости по порядку Л2/[^т] ~ 1 целесообразно придерживаться из соображений соизмеримости ошибок аппроксимации производных по времени и пространственным переменным.

Продолжим функцию симметрично: с(^г, я) = с(£, —г, г) (г < 0), т. е. допустим г € [—Ь,Ь], оставив то же уравнение диффузии. Непосредственно при г = 0 имеем за счет сг |+0 = 0 разрешимую особенность: с* = ^(2сгг + с,^), г = ±0. Формально «склеим» уравнение диффузии на [—Ь,Ь]. Аналогично соотношению (14) получаем

2(с-1,7 + с1,7)/к — со,7(1 + 2/к) + ^о,7 = °. Положив с-1,7 = С1,7, находим

Численные эксперименты показали, что все варианты нахождения начальных коэффициентов при расчете на достаточно мелкой сетке приводят к одинаковым результатам.

После выполнения прямого хода прогонки (вычисления «¿,7, вг,7) ближайшая цель — найти значение 7, необходимое для реализации обратного хода прогонки. Запишем аппроксимацию первого из граничных условий (5) при г = Ь (сг (¿, Ь, я) = 0) .В граничном узле с точностью до 0(Л_2) [1] имеем

2Лгсг(¿^+1/2, Ь, я7) ~см1_2,7 — 4с^1_1,7 + Зс^1,7.

Подставляя значения _2,7, _1,7- из соотно-

шения (9), находим

вЫ1, 7 (4 а^!-1,7 ) в^!-1,7

3 + а^1,7 (а^1-1,7 — 4)

По формуле (9) (с,,7 — аі+1,7 сі+1,7 + ^,+1,7, і — 0, ...,N1 — 1) определяем искомые значе-

Теперь найдем оставшиеся приближения с,,7 при і — 0 и і — Ж2, і — 0,1,..., N. Значения с,,^ — 0 согласно граничному условию (4) при г — Я (с(£, г, Я) — 0). Введем обозначение і0 — тах{і: г, ^ го}. Тогда из граничных условий (2), (3) (г — 0) получаем приближения с,,о — со, і — 0,1,..., іо и с,,о — (4с,,1 — с,,2)/3, і — іо + 1,..., N1.

Прогонка по переменной z

Поскольку в цилиндрических координатах возникает особенность при г ^ +0, и на границе г — 0 задано смешанное краевое условие, то переход с (к + 1/2)-го на (к + 1)-й слой по времени совершается в два шага.

Первый шаг: і — 1, г ^ +0 (аппроксимируем сг/г ~ сгг), реализуется алгоритм прогонки для уравнения с* — ^(2сгг + с,^). Аппроксимация уравнения с* — ^(2сгг + с,^):

¿1,7 с1,7 А1,7—1 2 ¿1,7 + А1,7+1

0,5 т^

с0,7 — 2 с1,7 + с2,7

Л2 .

+

+2

В обозначениях С — 2(1 + /|^т]),

— / Л \2

^1,7 — 2((с0,7 — 2 ¿1,7 + ¿2,7) +(С — 2)С1„

а1,7 — 4/(2 + к), в1,7 — ^0,7к/(2 + к). (16) получаем

¿1,7—1 — С ¿1,7 + ¿1,7+1 + ^^1,7 — 0, к ^ 0. (17)

Ищем приближенные значения концентрации на (к + 1)-м слое по времени в виде

¿1,7 — “1,7+1 ¿1,7+1 +^1,7+1, і — ^..^^ — 1. (18)

Подставляя в (17), находим прогоночные коэффициенты: і — 2, 3,..., N2,

1

“1,7 —

С — “1 7_

7—1

і^1,7—1 + в1,7—1 (19)

С — 01,7—1 . ( )

Начальные коэффициенты «1,1, Дд находим из формулы (18) при j = 0 и условия (2) (С1,0 = с0): а1,1 = 0, ^1,1 = со.

Второй шаг: г = 2,..., N1 — 1, г > 0, прогонка для уравнения (1) (с* = ^(сгг + сг/г + )).

Из разностного соотношения (7), обозначив

■^,7 —

1—2і^сі—1,7—

ния с

і — 0,..., N — 1, і — 1,..., N — 1.

■ 21 — £ )с

і,7

+ ( 1 + )сі+1,7

25

2

при каждом г = 2, 3,..., N1 — 1 получаем ¿¿,7-1 — ^¿¿,7 + ¿¿,7+1 + ^^¿,7 = 0, к ^ 0. (20)

Ищем значения концентрации в узлах сетки на (к + 1)-м слое по времени в виде

¿¿,7 — а.,7+1 ¿¿,7+1 + вг,7+1, ] — 0, ..., -^2. (21)

Прогоночные коэффициенты: ^ — 2, 3,..., N2,

а. . =_____1____ в- ■ = ^■?'-1 + в^,^-1 (22)

а-,7 = С _ а , в-,7 = С _ а . (22)

С а.,. —1 С а.,. —1

Начальные коэффициенты при г ^ г0 (г ^ го) находим из (21) (_?’ = 0) и условия (2) ¿¿,0 = с0: а.,1 = 0, в.,1 = Со. Начальные коэффициенты при г > г0 (г > г0) определяем из (20) при ] = 1 и условия (3) сг(¿,г, 0) = 0,г > г0:

а.,1 = 2 — С/2, в.,1 = ^^¿,1/2.

Граничные значения, необходимые для реализации обратного хода прогонки, находим из условия (4): ¿¿, = 0. По формулам (18), (21)

вычисляем значения ¿¿,. для г = 1,..., N1 — 1, ] = 0,... ,N2 — 1.

Найдем оставшиеся значения ¿¿,. при г = 0 и г = N1, ^ =0,1,..., N2. Используя второе из условий (5) (сг|+0 = 0), из аппроксимации (12) получаем ¿0. = ^¿^ — ¿2,.)/3. Аналогично, согласно первому условию (5) сг(¿,Ь,г) = 0 и

¿N^7 (4¿^V1 —1,. ¿^!-2,. )/3.

Как и на этапе прогонки по г можно предложить другой вариант вычислений с использованием дополнительного (фиктивного) узла. В этом случае значения концентрации на границе (г = 0) можно искать в цикле прогонки по переменной г, а при нахождении ^0. формально считать, что имеются узлы ¿—1,. (¿-1,:? = ¿1,.).

Вычислительный Алгоритм

Для контроля вычислений используем критерий материального баланса в объеме:

Г** г [• г0

=

—1 о д

- О — 2пг ^г

- .'о дг о

СГч д

- О — 2пг^г

- .'о дг В

гВ гЬ

Слева — разность количеств атомов водорода, растворившихся в образце на входной стороне сквозь дефект покрытия и покинувших выходную сторону, а справа — количество Н, находящихся в материале. Вместо можно использовать и другие «контрольные» значения ¿. Обозначим входной и выходной потоки через

Г Г

^о(^) = — О / Jо ■]В (і) = — О [

о

0 дс(і, г, 0) дг

дс(і, г, Я) дг

2пг^г,

2пг^г.

В программной реализации предусмотрен контроль дисбаланса относительно входного потока на уровне процента:

і, і, ВЬ

/іо(і) ¿і-/ ^я(і) ¿і —//с(і*, г, г)2пг гіг^г

о________________о__________________00__________________________________

і,

/ 70 (і)

о

с(і*, г, г) 2пг гіг^г.

оо

Для повышения точности результата вычислений использовалось правило Рунге-Ромберга [1]. Кратко приведем его реализацию в контексте рассматриваемой задачи. Первая и вторая производная по переменной г заменялись центральными разностными производными, поэтому, выполнив дополнительную серию расчетов с шагом по радиусу (Л,г/2), можно улучшить аппроксимацию в объеме на два порядка. Уточненные значения концентрации находим по формуле

С(*, г, г) = (4с^/2^, г, г) — (¿, г, г))/3. (23)

Чтобы сохранить общий порядок аппроксимации 0(^4 + Л>2 + т2), в граничных узлах сетки с точностью до О(Л^) (используем метод неопределенных коэффициентов [7]) полагаем для Сг

о

С

о

/, т ч 3с^ -4,7 - 16с -3,7 + 36с^1 -2,7 - 48с^1 -1,7 + 25с^,7

Сг& ^ ^) и —1-------------------1-----------24^-------------1-------------—. (24)

С26)

Приведем в компактной форме алгоритм вычислений. Переход от к-го слоя по времени к (к + 1)-му осуществляется в два этапа.

Х-й этап: на (к + 1/2)-м слое (к ^ 0) проводим следующие вычисления.

1. Вычисляем прогоночные коэффициенты а^-, ві,7 по формулам (13). Находим значение концентрации с^1,7' для граничного узла из аппроксимации сг (і, Ь, г) = 0 с точностью до 0(^4) по формуле (24). Обратным ходом прогонки восстанавливаем значения концентрации в оставшихся узлах.

2. В соответствии с граничными условиями (2), (3) и (4) определяем значения концентрации с^- для і = 0 и

І = N2, і = 0,1,..., Ж1.

3. Выполняем повторно пункты 1, 2 с шагом по радиусу Л,г/2. По правилу Рунге-Ромберга (формула (23)) находим уточненные значения концентрации во всех узлах слоя.

ХХ -й этап — переход на (к + 1)-й слой.

1. Для і = 1, 2,..., N1 — 1 вычисляем наборы прогоночных коэффициентов а:^, ві,і, і = 1, 2,..., N — 1 по формулам (19), (22). В силу граничного условия Сі,^2 = 0. Обратным ходом прогонки восстанавливаем значения концентрации в узлах сетки по формулам (18), (21).

2. В соответствии с граничными условиями сг |г=+о,ь = 0 определяем значения концентрации с^- (і = 0 и

і = N1, і = 0,1,..., N2) с точностью до 0(^4) согласно формуле (24).

3. Выполняем повторно пункты 1, 2 с половинным шагом по радиусу (на (к + 1/2)-м слое по времени). По правилу Рунге-Ромберга (см. (23)) находим уточненные значения концентрации во всех узлах слоя.

4. Проверяем выполнение материального баланса с заданной точностью.

Результаты моделирования

Безразмерная форма краевой задачи

Для численного моделирования перейдем к безразмерным переменным, используя естественные для данной задачи нормировки: и = с/с0, р = г/Ь, £ = г/Н, т = (^/Ь2)£. Значение Ь2/^ интерпретируется как характерное время диффузионного выравнивания концентрации [4] в области с линейными размерами порядка Ь (в вертикальном направлении, по оси г, диффузия при сравнимых Н быстрее вследствие вакуумирования с выходной стороны). Обозначив р0 = г0/Ь, К = Ь/Н, приходим к безразмерной краевой задаче:

ди 1 д / ди\ 2д2и

дт р др\др/ + д(2 ,

р е (0,1), ( е (0,1), т > 0, и(т,р,0) = 1, р е [0,р0], р0 < 1, ди

—(т,р, 0) = 0, р е (р0,1], т ^ 0, и(т, р, 1) = 0, р е [0,1], т ^ 0,

|(т,1.0=0, |(т,+о.с)=0.Се[о. 1]. и(0,р.С) = 0. р е [0,1], с е [0,1].

Варьируемыми являются К > 0 и р0 е (0, 1). Формально получили уравнение диффузии в анизотропной среде: = 1, ^ = К2. Чем

меньше Н (тоньше перегородка), тем быстрее диффузия в направлении £ (оси г). Коэффициент ^ отражает соотношение геометрических характеристик образца. За «время» т характерной «длиной» по р является величина у/т, а по £ — Ку7^. Установление (ит « 0) определяется значением т* = (^/Ь2)£*.

Введем безразмерный усредненный (по основанию площади п) выходной поток

д2 ди

д<

С =1

2пр^р.

Функция и (т) монотонно растет, выходя асимптотически на стационарное значение и* = и (т*). В отсутствие защитного покрытия (формально р0 = 1) имеем максимум и* = К2.

Исходный усредненный поток имеет вид

^ (і) = Осо пЬ2

Я

V (т), (т = ОЬ-2і),

1

1

ди

С=1

2прф.

27

Величина Dco/Н равна плотности проникающего стационарного потока в отсутствие защитного покрытия (ро = 1). Следовательно, число V(DL-2t) є (0,1) является долей J(t)/(nL2) в максимально возможной плотности выходного потока Dco/H. В пределе имеем максимум V* = R-2U* = 1 (р0 = 1).

Для определенности фиксируем диапазон коэффициента диффузии D: 10-7 — 10-5 см2/с, и будем рассматривать цилиндрическую перегородку, соизмеримую с мелкой монетой. В безразмерной модели остановимся на следующих значениях параметров:

ро = {1/50, 2/50,..., 1/5},

R = {1/3,1/2, 7/10,1, 3/2, 3, 5, 7,10}.

На начальных шагах по времени, когда в образце еще пренебрежимо мало водорода, но в приповерхностном объеме под дефектом покрытия уже устанавливается концентрация и(р, 0) = 1 (р ^ р0), происходит резкий скачок концентрации в узлах, прилегающих к дефекту, что приводит к росту погрешности по критерию материального баланса. Чтобы сгладить этот начальный всплеск, при решении безразмерной задачи концентрацию водорода под дефектом защитного покрытия наращиваем постепенно, но относительно быстро. Это соответствует физическим предположениям и лишь формализовано в модели начальным скачком концентрации. Для относительно малых шагов Ат (порядка секунды в исходном времени t — за это время ничего значительного в «копейке» не произойдет) определим и(Ат, р, 0) = 10-3, и(2Ат, р, 0) = 10-2, и(3Ат,р, 0) = 10-1, и(4Ат,р, 0) = 0,5, и(5Ат, р, 0) = 1, р ^ р0. За пять шагов выходим на уровень и = 1 под дефектом и за счет этого сглаживаем начальный модельный всплеск входного потока. Формально при контроле материального баланса сдвигаем время начала эксперимента: 0 ^ 5Ат.

Выходными параметрами эксперимента являются т*, т0, U*, V*, где т0 — время запаздывания, вычисляемое по формуле

Т0 = т* — S(t*)/S(t*) = т* — S (t*)/U (т*),

S(т) =/ U(т) dT (S = dS/dT).

0

В исходном времени t0 = t* — Jg* J(t) dt/J (t*). Геометрически это точка пересечения асимптоты графика количества проникающего

сквозь перегородку водорода с осью ¿. С учетом асимптотического характера выхода на стационарный режим проницаемости точность вычислений т0(£0) возрастает с ростом т*(£*). Время выхода на стационар, значения установившегося потока и время запаздывания — регистрируемые величины, являющиеся входной информацией при решении обратной задачи оценки кинетических параметров проницаемости по экспериментальным данным.

Критерий останова и тестирование

При выборе критерия останова использовались два условия:

{и(т„) > а(р0К)2} Л { и(т"> п1 < е}.

Здесь первое условие необходимо, чтобы исключить останов в начале эксперимента, когда выходной поток очень мал. Ориентиром является тахи(т*) = К2 (р0 = 1) и пропорциональность (в линейном приближении) потока площади дефекта (~ р0). Коэффициент а является эмпирическим, в рассматриваемом диапазоне параметров модели использовалось значение а = 1/30. Второе условие срабатывает, когда изменение потока относительно мало в течение «заметного» времени реального эксперимента (например, — ¿га-п ~ 1 мин.).

Второй вариант критерия останова:

ГГтг-\ ^ „Ь Г и(тП — и(т«-п^ _} {(./(т„) < ^ Л {-------С7д2-----< 4

В этом случае первое условие учитывает, что графики усредненного потока имеют точку перегиба: сначала кривая выпуклая, а по мере насыщения образца становится вогнутой. Цель первого условия в этом критерии та же — не остановить вычисления преждевременно.

В данной работе решается прямая задача с целью выявления качественных зависимостей величин и*, т*, т0 от параметров модели.

На рис. 1-3 представлены результаты численных экспериментов для различных С ([О] = ем2/е): С1 = р0 (моменты останова для соответствующего значения е обозначены О), С2 = и(тга-п)(■), Сз = тга — тга-п (А). Результаты для различных К, р0, О показали, что на модельных кривых наиболее приемлемо ограничение «производной» выходного потока (А). В последующих расчетах использовался второй критерий при С = Сз и е = 10-з.

1.01

о ю-1доЛ—,10"* ■ ю^до^.-ло'6

А 10'1,10'2,...,10‘6

0.00 0.05 0.10 0.15 0.20 0.25

Рис. 1. Останов: О = 10-5, К =5, р0 = 0, 2

1.00

0.99

0.33 1 Ь5 3 5 7

Рис. 4. Контроль: то, р0 = 1

и

и

Рис. 2. Останов: О = 10 7, К =5, р0 = 0, 08

Рис. 5. Выбор т*, К < 1, 5

и

и

Рис. 3. Останов: О = 10 5, К = 0, 8, р0 = 0, 2

Дополнительно проведены контрольные расчеты при р0 = 1 — это соответствует пластине без защитного покрытия. Расчеты подтвердили, что выполняется асимптотический

Рис. 6. Выбор т*, К ^ 1,5

выход на линейное стационарное распределение и(т, р, С) = 1 — £ (т ^ т*). Кроме того, как известно, время запаздывания проникающего потока сквозь пластину толщины £

■©

равно £2/[60] [3], в нашем случае 1/[6К2} (рис. 4). Также проверена асимптотика и* ^ К2 при ро ^ 1. Правильность автоматического выбора момента окончания эксперимента для различных К проиллюстрирована на рис. 5-6.

На рис. 7-16 представлены результаты численного моделирования при варьировании параметров К и ро. По итогам эксперимента в безразмерной форме рассматриваем три основные характеристики т*, т0, и* (V* = К-2и*). Каждый из выходных параметров эксперимента представлен на двух рисунках: отдельно для К ^ 1,5 и К ^ 1,5, так как полученные значения на этих промежутках отличаются на порядок. По умолчанию очередность перечисления значений соответствует убыванию максимумов. «Водораздел» К = 1, 5, разумеется, достаточно условен. В широком диапазоне параметров (в том числе и П) на качественном уровне следует говорить об относительно малых и больших К и ро.

На рис. 7-10 заметна слабая зависимость (безразмерных) временных характристик т*, то от радиуса дефекта. Пунктиром обозначены графики с повторяющимися значениями К на смежных рисунках. В пределах ро < 0,2 (при ро ^ 1 уже трудно говорить о наличии защитного покрытия) наблюдается экстремум функций т*(ро), то(ро) при К < 1, 5 (смена знака производных на рис. 11, 12). При К ^ 1, 5 сохраняется монотонность. По рисункам 13, 14 можно оценить скорость изменения стационарного уровня и* при вариациях К и ро.

В целом, на качественном уровне (ориентируясь на экспериментальные погрешности в несколько десятков процентов) можно утверждать, что при К = Ь/Н > 3 влияние граничных условий при г = Ь незначительно и можно практически считать цилиндр пластиной (Ь ^ +го). Это согласуется с результатом работы [12] и позволяет использовать аналитические методы исследования [10].

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

уровня стационарного режима водородопро-ницаемости. Разностная схема для случая, когда существенно влияние поверхностных процессов, представлена в [6].

10

10

10° _

10

Я = 1/3, 0.5, 0.7, 1, 1.5

—"

0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20

Рис. 7. Время эксперимента, К ^ 1, 5

ро

10

10° I

10

Я = 1.5, 3, 5, 7, 10

0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20

Рис. 8. Время эксперимента, К ^ 1, 5

ю

10

10

Я = 1/3, 0.5, 0.7, 1, 1.5

0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20

Рис. 9. Время запаздывания, К ^ 1, 5

Л

10

10

10"

Д = 1.5, 3, 5, 7, 10

0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20

Рис. 10. Время запаздывания, К ^ 1, 5

ро

0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.201

Рис. 13. Стационарный поток, К ^ 1, 5

(к*/с1р{

. \\ Я = 1/3, 0.5, 0.7, 1, 1.5 6_ ^ = 10, 7, 5, 3, 1.5 /

5-

4-

о.о8 о.1о \ о.о\ о^к. олб Тт Ро 3-

\ \ 2-

\ 1-

\ 0_

Рис. 11. Скорость изменения т*, К ^ 1, 5

0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20

Рис. 14. Стационарный поток, К ^ 1, 5

Ро

т0/т*

Рис. 15. Доля времени запаздывания, К ^ 1, 5

То/Т*

Рис. 16. Доля времени запаздывания, R ^ 1, 5

Литература

1. Волков Е. А. Численные методы. М.: Наука, 1987. 248 с.

2. Косарев В. И. 12 лекций по вычислительной математике. М.: МФТИ, 2000. 224 с.

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

4. Ландау Л. Д., Ахиезер А. И., Лифшиц Е. М. Курс общей физики. Механика и молекулярная физика. М.: Добросвет; КДУ, 2011. 338 с.

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

6. Родченкова Н. И., Костикова Е. К. Разностная схема для краевой задачи водородопрони-цаемости при наличии дефекта защитного покрытия // Труды Карельского научного центра РАН. Сер. Математическое моделирование и информационные технологии. Вып. 2. 2011. № 5. С. 97-102.

7. Рябенький В. С. Введение в вычислительную математику. М.: Физматлит, 2000. 296 с.

8. Самарский А. А., Вабищевич П. Н. Вычислительная теплопередача. М.: Едиториал УРСС, 2003. 784 с.

9. Самарский А. А. Введение в теорию разностных схем. М.: Наука, 1971. 553 с.

10. 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.

11. 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.

12. Zajec B. Hydrogen permeation barrier-recog-nition of defective barrier film from transient permeation rate // International Journal of Hydrogen Energy. Elsevier, 2011. Vol. 36. P. 73537361.

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

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

зав. лаб. моделирования природно-технических систем, д. ф.-м. н.

Институт прикладных математических исследований Карельского научного центра РАН ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: [email protected] тел.: (8142) 766312

Zaika, Yury

Institute of Applied Mathematical Research, Karelian Research Centre, Russian Academy of Sciences 11 Pushkinskaya St., 185910 Petrozavodsk, Karelia, Russia

e-mail: [email protected] tel.: (8142) 766312

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

младший научный сотрудник

Институт прикладных математических исследований Карельского научного центра РАН ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: [email protected] тел.: (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: [email protected] tel.: (8142) 766312