CC BY
13Статья поступила в редакцию 17.04.11. Ред. рег. № 976
The article has entered in publishing office 17.04.11. Ed. reg. No. 976
УДК 539.219.3:669.788
МОДЕЛИРОВАНИЕ ОДНОМЕРНОЙ ДИФФУЗИИ ВОДОРОДА В МЕТАЛЛАХ. IV. ПРОНИЦАЕМОСТЬ ПРИ КОЭФФИЦИЕНТЕ ДИФФУЗИИ, ЗАВИСЯЩЕМ ОТ КОНЦЕНТРАЦИИ ДИФФУЗАНТА
1 2 В.Н. Лобко , И.Н. Бекман
владимирский государственный университет 600000 Владимир, ул. Горького, д. 87 Тел. (4922)47-98-67, e-mail: lobko_vn@laser-2.vpti.vladimir.ru 2МГУ им. М.В. Ломоносова 119991 Москва, ГСП-1, Ленинские горы, д. 1, стр. 3 Тел. (495)939-32-12
Заключение совета рецензентов: 27.04.11 Заключение совета экспертов: 28.04.11 Принято к публикации: 30.04.11
В первой статье серии были предложены два варианта численных методов моделирования одномерной диффузии водорода в плоскопараллельной металлической пластине, контактирующей с постоянными объемами, при отсутствии информации о краевых условиях на самой пластине. В настоящей работе описанная методика применена для случая проницаемости при коэффициенте диффузии, зависящем от концентрации водорода. Рассмотрены некоторые виды концентрационной зависимости в сравнении с близким по значению постоянным коэффициентом диффузии.
Ключевые слова: диффузия, водород, металлы, моделирование, коэффициент диффузии, проницаемость, метод конечных разностей, концентрационная зависимость коэффициента диффузии, концентрационные профили.
NUMERICAL MODELING OF ONE-DIMENSIONAL DIFFUSION OF HYDROGEN IN METALS. IV. PERMEABILITY AT THE COEFFICIENT OF DIFFUSION DEPENDING ON
DIFFUSANT CONCENTRATION
V.N. Lobko1, I.N. Beckman2
'Vladimir State University 87 Gorky st., Vladimir, 600000, Russia Tel. (4922)47-98-67, e-mail: lobko_vn@laser-2.vpti.vladimir.ru 2Moscow State University 1 Leninskiye Gori, Moscow, GSP-1, 119991, Russia Tel. (495)939-32-12
Referred: 27.04.11 Expertise: 28.04.11 Accepted: 30.04.11
In the first article the mathematical foundations of the method for cases of constant diffusion coefficient and its dependence on the arbitrary concentration of hydrogen were presented. In this work the method for the case of the permeability with the diffusion coefficient, depending on the concentration of hydrogen, was described. Some types of concentration dependence in comparison with the cognate constant diffusion coefficient were considered.
Keywords: diffusion, hydrogen, metals, modeling, coefficient of diffusion, permeability, finite difference method, concentration dependent of diffusion coefficient, profile of concentration of hydrogen.
Введение
В практике инженерных расчетов и теоретического анализа диффузионных задач часто бывает затруднительным использование известных решений для краевых задач из-за отсутствия информации о краевых условиях. Поэтому важным является разработка методов моделирования диффузионных про-
цессов, исходя лишь из геометрических размеров диффузионных систем и табличных констант.
В первых статьях серии были представлены математические основы решения задачи моделирования диффузии в плоскопараллельной пластине при постоянном коэффициенте диффузии Б и при коэффициенте диффузии, зависящем от концентрации водорода Б(с) [1], а также проведена проверка алгоритмов при
International Scientific Journal for Alternative Energy and Ecology № 4 (96) 2011
© Scientific Technical Centre «TATA», 2011
D = const для случаев сорбции-десорбции и проницаемости [2, 3]. В настоящей статье представлена апробация методики на примере газопроницаемости пластины при различных видах функции D(c). Идеология метода основана на использовании известных разностных схем для краевых задач с добавлением уравнений, описывающих параметры диффузионной системы. Алгоритм представлен в двух вариантах: 1-й вариант - на основе дифференцирования и 2-й вариант - на основе интегрирования, которые приводят к двум разным разностным схемам. Сравнительное использование этих схем служит критерием корректности моделирования.
Частным случаем для предложенных схем является D = const, поэтому сначала целесообразно сравнить моделирование по формулам для постоянного коэффициента диффузии D c моделированием диффузии при D, зависящем от концентрации c.
Рис. 1. Сравнение моделирования проницаемости по формулам для D = const и по формулам для D = f(c) (расхождений нет). Выходной поток J в зависимости от времени t. В обоих случаях применен и 1-й, и 2-й варианты моделирования Fig. 1. Comparison of numerical modeling of permeability under formulas for D = const and under formulas for D = f(c) (there aren't divergences). Output flux J depending on time t.
In this cases it is applied both 1st and 2nd versions of numerical modeling
Далее рассмотрено моделирование газопроницаемости для некоторых видов D = fc) по первому и второму вариантам алгоритма для двухатомного (диссоциирующего в материале пластины, например, водорода) газа в сравнении с м2/с. (Аналогичные расчеты были проведены и для одноатомного газа.) Рассмотрены линейная, экспоненциальная, параболическая и гиперболическая зависимости коэффициента диффузии газа от его концентрации. Разработанные программы позволяют обрабатывать практически любые виды D(c) (возможно, за исключением некоторых особых случаев). Поскольку наиболее показательным параметром являются выходные потоки газа в методе проницаемости, приведена именно эта зависимость J(t), хотя в программе, кроме этого, рассчитываются концентрационные профили, временные зависимости концентраций и потоков на поверхностях, а также - давлений в резервуаре и
приемнике. Во всех последующих примерах расчеты проводили при следующих значениях параметров: толщина пластины, 1,5-10-3 м; температура, 680,7 К; площадь пластины, 3,14-10-4 м2; постоянный коэффициент диффузии, 5,66-10-10 м2/с; константа Сивер-тса, 2,45-10-2 моль/м3/Па0,5; объем резервуара, 5-10-2 м3; объем приемника, 5-10-5 м3; давление в резервуаре, 8,95-104 Па; давление в приемнике, = 0 Па; начальная концентрация в пластине, 1-10-12 моль/м3; число разбиений по оси x, 1000; конечное время 10000 с; число разбиений по оси t, 500; точность метода итераций Ньютона, 10-6. Входная концентрация составляла 7,38 моль/м3, выходная в большинстве случаев - в 10 раз меньше.
На рис. 1 приведены результаты моделирования при D = const = 5,66-10-10 м2/с по схемам для постоянного D и для D = fc) (в последнем случае fc) = const). Получено полное совпадение функций J(t) в пределах погрешности расчетов как для 1-го варианта моделирования, так и для 2-го, то есть по трем разным схемам, и это свидетельствует о корректности используемых разностных схем.
Линейная зависимость D(c)
Для моделирования диффузии при линейной зависимости D(c) были взяты следующие функции, изображенные на рис. 2:
1) D = (5,66 + 0,5c)-10-10 м2/с, (рис. 3, кривая 2);
2) D = (-0,5c + 9,41c)-10-10 м2/с, (рис. 4, кривая 2);
3) D = (0,5c + 1,91c)-10-10 м2/с, (рис. 4, кривая 3);
4) D = (5,66 - 0,5c)-10-10 м2/с, (рис. 3, кривая 3).
Рис. 2. Некоторые виды линейной зависимости коэффициента диффузии D от концентрации водорода, использованные для расчета кинетических кривых нестационарной проницаемости: 1 - D = (5,66 + 0,5c)-10-10 м2/с; 2 - D = (-0,5c + 9,41 c)-10-10 м2/с; 3 - D = (0,5c + 1,91 c)-10-10 м2/с; 4 - D = (5,66 - 0,5c)-10-10 м2/с. Пунктир - D = 5,66-10-10 м^ Fig. 2. Some kinds of linear dependence of coefficient of diffusion D on hydrogen concentration, used for calculation of kinetic curves of non-stationary permeability: 1 - D = (5,66 + 0,5c)-10-10 m2/s; 2 - D = (-0,5c + 9,41 c)-10-10 m2/s; 3 - D = (0,5c + 1,91 c)-10-10 m2/s; 4 - D = (5,66 - 0,5c)-10-10 m2/s . Shaped line - D = 5,66-10-10 m2/s
Международный научный журнал «Альтернативная энергетика и экология» № 4 (96) 2011 © Научно-технический центр «TATA», 2011
Рис. 3. Моделирование проницаемости водорода сквозь пластину при линейной функции D(c). Выходной поток J в зависимости от времени t.
1 - D = const = 5,66-10-1° м2/с;
2 - D = (5,66 + 0,5c)-10-1° м2/с;
3 - D = (5,66 - 0,5c)-10-1° м2/с.
В обоих случаях применен и 1-й, и 2-й варианты моделирования (расхождений нет) Fig. 3. Numerical modeling of permeability of hydrogen through a plate at linear function D(c). Output flux J depending on time t: 1 - D = const = 5,66-10-1° m2/s; 2 - D = (5,66 + 0,5c)-10-1° m2/s; 3 - D = (5,66 - 0,5c)-10-1° m2/s. In this cases it is applied both 1st and 2nd versions of numerical modeling (there aren't divergences)
Рис. 4. Моделирование проницаемости при линейной функции D(c). Выходной поток J в зависимости от времени t. 1 - D = const = 5,66-10-1° м2/с; 2 - D = (-0,5c + 9,41 c)-10-1° м2/с; 3 - D = (0,5c + 1,91 c)-10-1° м2/с. В обоих случаях использовали и 1-й, и 2-й варианты моделирования (расхождений нет) Fig. 4. Numerical modeling of permeability at linear function D(c). Output flux J depending on time t. 1 - D = const = 5,66-10-1° m2/s; 2 - D = (-0,5c + 9,41 c)-10-1° m2/s;
3
■ D = (0,5c + 1,91 c)-10"10 m2/s. In this cases it is applied
которым проводится сравнение. При ранних временах ход кривых J(t) почти одинаков и они имеют сравнимый инкубационный период.
На рис. 4 показан такой же случай, но линейные зависимости D(c) равны между собой при больших концентрациях. Зависимость J(t), соответствующая кривой 3, возрастающая, но сама кривая временной зависимости выходного потока имеет большой инкубационный период и лежит значительно ниже кривой для D = const. Кривая 2 (рис. 4) имеет заметно меньший инкубационный период.
Экспоненциальная зависимость В(е)
Рассмотрим экспоненциальную зависимость Б от концентрации газа для функций (рис. 5):
1) D = 5,66exp(0,2c)-10 м2/с;
2) D = 0,3exp(0,5c)-10"10 м2/с;
3) D = 0,3exp(0,7c)-10_1° м2/с.
Рис. 5. Некоторые виды экспоненциальной зависимости коэффициента диффузии й от концентрации водорода, использованные для расчета кинетических кривых проницаемости: 1 - й = 5,66ехр(0,2с)-10"1° м2/с;
2 - D
■■ 0,3exp(0,5c)-10"10 м2/с;
both 1st and 2nd versions of numerical modeling (there aren't divergences)
На рис. 3 представлены результаты расчетов временных зависимостей выходных потоков газа для линейных зависимостей Б(с), возрастающей и падающей, которые при малых концентрациях диффу-занта равны между собой и равны постоянному Б, с
3 - D = 0,3exp(0,7c)-10"10 м2/с. Пунктир 4 - D = const = 5,66-10"10 м2/с Fig. 5. Some kinds of exponential dependences of coefficient of diffusion D on concentration of the hydrogen, used for calculation of kinetic curves of permeability:
1 - D = 5,66exp(0,2c)-10"10 m2/s;
2 - D = 0,3exp(0,5c)-10"1° m2/s;
3 - D = 0,3exp(0,7c)-10"10 m2/s. Shaped line 4 - D = const = 5,66-10"10 m2/s
Результаты моделирования представлены на рис. 6. Обращает на себя внимание очень большой инкубационный период зависимости 2 и ход зависимости 3, когда выходной поток J(t) сначала сильно отстает от кривой для D = const, а затем быстро обгоняет ее. Это объясняется быстрым ростом D с концентрацией газа (см. рис. 5, кривая 3).
International Scientific Journal for Alternative Energy and Ecology № 4 (96) 2011
© Scientific Technical Centre «TATA», 2011
Рис. 6. Моделирование проницаемости при экспоненциальной функции D(c). Выходной поток J в зависимости от времени t. 1 - D = 5,66exp(0,2c)-10-10 м2/с;
2 - D = 0,3exp(0,5c)-10-10 м2/с; 3 - D = 0,3exp(0,7c)-10-10 м2/с; 4 - D = const = 5,66-10-10 м2/c. Во всех трех случаях применен
и 1-й, и 2-й варианты моделирования (расхождений нет) Fig. 6. Numerical modeling of permeability at exponential function D(c). Output flux J depending on time t. 1 - D = 5,66exp(0,2c)-10-10 m2/s; 2 - D = 0,3exp(0,5c)-10"1° m2/s;
3 - D = 0,3exp(0,7c)-10"10 m2/s; 4 -
D = const = 5,66-10"10 m2/s.
Рис. 7. Некоторые виды параболической зависимости коэффициента диффузии D от концентрации газа, использованные для расчета кинетических кривых проницаемости: 1 - D = (1,3c2 - 9,75c + 18,78)-10"10 м2/с;
2 - D = (-1,3c2 + 9,75c + 1,5)-10-10 м2/с.
Пунктир 3 - D = const = 5,66-10-10 м2^ Fig. 7. Some kinds of parabolic dependences of coefficient of diffusion D on concentration of the hydrogen, used for calculation of kinetic curves of permeability: 1 - D = = (1,3c2 - 9,75c + 18,78)-10-10 m2/s; 2 - D = (-1,3c2 + 9,75c + + 1,5)-10-10 m2/s. Shaped line 3 - D = const = 5,66-10-10 m2/s
Результаты показаны на рис. 8. Здесь выходной поток газа, рассчитанный для D(c) вида 2, значительно превышает поток, рассчитанный для постоянного коэффициента диффузии, а поток зависимости 1 - мало от отличается от случая D = const. Зависимость вида 2 соответствует кривой 2 (рис. 7); наблюдаемый эффект объясняется резким повышением коэффициента диффузии D при малых концентрациях газа, его последующее резкое снижение заметно не сказывается на ходе кривой выходного потока J(t).
In all three cases it is applied both 1st and 2nd versions of numerical modeling (there aren't divergences)
Параболическая зависимость D(c)
Параболическая зависимость D, когда коэффициент диффузии с ростом концентрации сначала убывает (возрастает), а затем возрастает (убывает) (см. рис. 7), рассмотрена для вариантов:
1) D = (1,3c2 - 9,75c + 18,78)-10-10 м2/с;
2) D = (-1,3c2 + 9,75c + 1,5)-10-10 м2/с.
Рис. 8. Моделирование проницаемости при параболической функции D(c). Выходной поток J в зависимости от времени t.
1 - D = (1,3c2 - 9,75c + 18,78)-10-10 м2/с;
2 - D = (-1,3c2 + 9,75c + 1,5)-10-10 м2/с;
3 - D = const = 5,66-10-10 м2/а В обоих случаях применен и 1-й, и 2-й варианты моделирования (расхождений нет) Fig. 8. Numerical modeling of permeability at exponential function D(c). Output flux J depending on time t.
1 - D = (1,3c2 - 9,75c + 18,78)-10"10 m2/s;
2 - D = (-1,3c2 + 9,75c + 1,5)-10-10 m2/s;
3 - D = const = 5,66-10"10 m2/s. In this cases it is applied both 1st and 2nd versions of numerical modeling (there aren't divergences)
Гиперболическая зависимость D(c)
Гиперболическая зависимость D взята для гипербол, сдвинутых влево по оси абсцисс таким образом, что вертикальная ветвь пересекает ось ординат (см. рис. 9), по двум вариантам:
1) D = (2/(с + 0,2) + 5,6)-10-10 м2/с;
2) D = (2/(с + 0,08) + 5,6)-10-10 м2/с.
Результаты приведены на рис. 10. Рассчитанные кривые J(t) мало отличаются от тестовой кривой при D = const, заметное отличие - только в малом инкубационном периоде. Причем чем круче возрастает гипербола при приближении концентрации к нулю, тем меньше инкубационный период или отсутствует вообще. Конечно, сказанное справедливо для конкретной толщины пластины.
Международный научный журнал «Альтернативная энергетика и экология» № 4 (96) 2011 © Научно-технический центр «TATA», 2011
Рис. 9. Некоторые виды гиперболических зависимостей коэффициента диффузии D от концентрации газа, использованные для моделирования проницаемости:
1 - D = (2/(c + 0,2) + 5,6)-10-1° м2/с;
2 - D = (2/(c + 0,08) + 5,6)-10-1° м2/с. Пунктир 3 - D = const = 5,66-10-10 м2/c
Fig. 9. Some kinds of hyperbolic dependences of coefficient of diffusion D on concentration of the hydrogen, used for numerical modeling of permeability: 1 - D =
: (2/(c + 0,2) + 5,6)-10"lu m2/s; 2 - D = (2/(c + 0,08) + 5,6)-10"10 m2/s. Shaped line 3 - D = const = 5,6610
Рис. 10. Моделирование проницаемости при гиперболической функции D(c). Выходной поток J в зависимости от времени t. 1 - D = (2/(c + 0,2) + 5,6>10-1° м2/с; 2 - D = (2/(c + 0,08) + 5,6)-10-10 м2/с; 3 - D = const = 5,66-10-10 м2/а В обоих случаях применен и 1-й, и 2-й варианты моделирования (расхождений нет) Fig. 10. Numerical modeling of permeability at hyperbolic function D(c). Output flux J depending on time t. 1 - D = (2/(c + 0,2) + 5,6)-10-10 m2/s;
2
■ D = (2/(c + 0,08) + 5,6)-10"10 m2/s;
3 - D = const = 5,66-10 m2/s. In this cases it is applied both 1st and 2nd versions of numerical modeling (there aren't divergences)
Заключение
Анализ полученных результатов показывает, что различного вида зависимости Б(с) [4, 5] очень сильно влияют на вид кривых выходного потока диффундирующего водорода ДО в режиме нестационарной проницаемости пластины. Следовательно, при решении обратных задач диффузии (в частности, определения коэффициента Б и его возможной зави-
симости от концентрации диффузанта) целесообразно использовать экспериментальные данные по нестационарному потоку, J(t), так как они объективно содержат нужную информацию. При этом требуется разработка математического аппарата решения обратных задач для этого случая.
Для экспоненциальной зависимости двух типов (кривая 2 и 3 рис. 6) получены очень большие инкубационные периоды. Для кривой J(t) с D(c) типа 3 мы имеем ситуацию, когда после большого инкубационного периода выходной поток водорода резко возрастает и затем медленно снижается. Этот эффект может быть использован в технологических целях. При этом моделирование помогает подобрать параметры диффузионной системы для реализации подобных случаев. Отметим, что если экспоненциальная зависимость D(c) возрастает слишком быстро, на смоделированной кривой в зоне резкого возрастания кривой наблюдается разрыв функции.
В случае гиперболической зависимости D(c) обращает на себя внимание практически полное отсутствие инкубационного периода.
Главным выводом является то, что во всех приведенных примерах оба алгоритма показали устойчивость и приводили к практически идентичным результатам. Устойчивость схем математически не доказана, и пока нельзя отдать предпочтение ни 1-му, ни 2-му вариантам моделирования. При практическом использовании схем для моделирования диффузии при D = fc) их устойчивость каждый раз должна проверяться эмпирически.
Поскольку предложенные в работе алгоритмы являются более универсальными, они способны заменить расчетные схемы для D = const. Тем не менее, при D = const лучше использовать схемы для постоянного D, так как расчеты по ним более быстры.
Список литературы
1. Лобко В.Н., Бекман И.Н. Моделирование одномерной диффузии водорода в металлах. I. Метод расчета изотермической диффузии при постоянном и переменном коэффициенте диффузии // Альтернативная энергетика и экология - ISJAEE. 2010. № 10. С. 36-42.
2. Лобко В.Н., Бекман И.Н. Моделирование одномерной диффузии водорода в металлах. II. Сорбция и десорбция в плоскопараллельной пластине, находящейся в замкнутом объеме // Там же. 2010. № 11. С. 38-42.
3. Лобко В.Н., Бекман И.Н. Моделирование одномерной диффузии водорода в металлах. III. Проницаемость через плоскопараллельную пластину, контактирующую с замкнутыми объемами // Там же. 2011. № 4. С. 15-19.
4. Водород в металлах / Под ред. Г. Алефельда и И. Фелькля. В 2-х томах. М.: Мир, 1981.
5. Райченко А.И. Математическая теория диффузии в приложениях. Киев: Наукова думка, 1981.
— TATA — LXJ
International Scientific Journal for Alternative Energy and Ecology № 4 (96) 2011
© Scientific Technical Centre «TATA», 2011
