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

УДК 550.837

Е.Н. Минаев, В.Ф. Пулин, И.В. Беляев

ОБ ОДНОМ МЕТОДЕ РЕШЕНИЯ СТАЦИОНАРНЫХ КРАЕВЫХ ЗАДАЧ МАТЕМАТИЧЕСКОЙ ФИЗИКИ С НЕЛИНЕЙНЫМИ ГРАНИЧНЫМИ УСЛОВИЯМИ

Представлен численно-аналитический метод, позволяющий свести стационарную краевую задачу к дискретному аналогу интегрального нелинейного уравнения. На конкретном примере показана возможность применения метода при расчете стационарного поля.

E.N. Minaev, V.F. Pulin, I.V. Belyaev

MATHEMATICAL PHYSICS STATIONARY BOUNDARY TASK SOLUTION

METHOD

UNDER NONLINEAR BOUNDARY CONDITION

The method of solution of stationary boundary task of mathematical physics under nonlinear boundary condition is presented in this paper. It is used for calculation of electrical field.

Часто в технических приложениях возникает необходимость расчета стационарного распределения температуры, концентрации растворённого вещества, электрического поля и т.д. С математической точки зрения такие задачи сводятся к решению уравнения Лапласа или Пуассона с соответствующими граничными условиями. При нелинейных граничных условиях, как правило, сначала проводят некоторые аналитические преобразования, а затем применяют численные методы. Например, в работе [1] к краевой задаче применяют метод собственных функций, получают нелинейное интегральное уравнение, которое затем решают численно. Проблема заключается в том, что далеко не всегда возможно решить полученную систему нелинейных уравнений. Например, если система решается итерационными методами, не всегда можно гарантировать сходимость итерационного процесса к решению. В данной работе предложен метод, позволяющий заменить краевую задачу дискретным аналогом интегрального уравнения. В ряде случаев такой подход обеспечивает сходимость итерационного процесса.

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

С учетом осевой симметрии и периодичности расположения анодов, запишем краевую задачу для электростатического потенциала ф = ф(г, х):

д 2ф 1 дф д 2ф _ _ _ „ ,1Ч

—-Н-------— +-т- = 0, 0 < х < xN , 0 < r < R, (1)

dr r dr дх

дф(0, г ) = 0 дф( ^, г ) = 0

дx

дx

_к ддф = 7(Ф) 5

дr

(2)

(3)

где х - продольная координата; г - радиальная координата; xN - половина расстояния между соседними анодами; R - радиус трубы; X - удельная электропроводность раствора.

Выражение (3) представляет собой нелинейную зависимость плотности тока ] от потенциала в приэлектродном слое жидкости. Целью расчета является нахождение распределения плотности тока по координате x в приэлектродном слое г = R.

В отличие от известных методов [2], основанных на эквивалентной замене краевой задачи (1)-(3) интегральным нелинейным уравнением (при помощи функций Грина или метода собственных функций), проведем замену непрерывного граничного распределения плотности тока дискретными значениями с самого начала при постановке задачи ( Л, п = 1,2,3..^ )

дф( x, R)

-к

дr

= л (Фп) ; ^-1 < x < х •

(4)

С учетом граничных условий (2) применим к дифференциальному уравнению (1) и граничному условию (3) интегральное косинус-преобразование Фурье по координате x

хМ

Ф* (Г )= |Ф(Х5 Г )

СОБ

Г *п

х Iёх, * = 1,2,3...

V хк

(5)

и получим для изображения фк (г) обыкновенное дифференциальное уравнение (6) с граничным условием (7)

V хм.

ёг2 г ёг

¿Фк (Я) = Хм м . Г • Г *п

ф* =0;

Б1П

V хм

. Г *п

-Х,

V хм

п_1

ёг Хпк п=1

Решение данного уравнения известно из теории специальных функций [3]

Ф*(г ) = А10

Г * п

----г 1+ ВКК0

V хм )

Г * п

V хм

(6)

(7)

(8)

Г *п

где 10 — г I - модифицированная функция Бесселя нулевого порядка; К0

Г*

V хм

п

-Г I -

V хм

функция Макдональда нулевого порядка; Ак и Вк - константы интегрирования. Учитывая, что в точке г = 0 функции фк и /0 имеют конечные значения, а К0 - бесконечна, получим Вк = 0. Используя далее граничное условие (7), находим Ак.

А* = _‘

кп2 *2 /1(* п Я /хм )

• Г. Г *

Б1П

п

V хм )

_ Б1П

Г *

п

V хм

))

(9)

Возвращаясь далее от изображения к самому потенциалу в соответствии с формулой обратного преобразования Фурье

Ф(х, г ) =

Ф0

+ -

хм *=1

ІФ* (г )со§

* п

V хм

(10)

г

х

х

п_1

п

х

х

N

и учитывая равномерную сходимость ряда (что позволяет поменять местами суммирование по к и по п), получим потенциал

ф(х,Г) = С-|% + У],У Х)

АП п=1 к=11 ( п Я / хМ)

Г . Г к п

Б1П

V хм

. Г к п

к2

СОБ

Г к п

-х

(11)

Рассматривая последнее выражение на границе г = Я и усредняя потенциал на участках хтА< х < хт ,

1 Хт

ф т =^Х 1ф(Х’ Я)ЙХ ’

Хт-1

приходим к системе нелинейных уравнений относительно усредненных потенциалов фт

Фт = #1 -

2 Х 2 N

М ЕАт,п • ]я(ф„); т = 1,2,...#,

Ап3 Ах п=1

(12)

(13)

где коэффициенты Лт,п определяются из числового ряда

Л = У

т,п

Б1П

- Б1П

V хы J

V хм

JJ

Б1П

- Б1П

Vхм J

V хм

ЛЛ 0

Iо( — Я)

х

JJ 11( — Я)

-, (14)

а С1 - неизвестная константа, определяемая из дополнительных условий. Для рассматриваемой задачи расчета электрохимической защиты находим ее из условия равенства потенциала в точке хМ минимальному защитному потенциалу ф(хМ) = фт1П. Система (13) преобразуется к виду

фт = фшт

^ У(ат,п -АЛп)• .)„(фп), т = 2,3,... (М-1),

Ах п3А

п=1

N

к =-У л.

(15)

(16)

В формуле (16)к - плотность тока на аноде, а сама эта формула выражает условие электронейтральности.

В случае плоской задачи система уравнений (15) остается без изменений, но в коэффициентах Лт,п отношение модифицированных функций Бесселя заменяется единицей.

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

1) исходная краевая задача, заданная в области, заменяется на эквивалентную, заданную на границе этой области и требующую дискретизации только на границе;

2) значение искомой функции внутри области определяется её значениями на границе при помощи аналитического выражения. В рассматриваемом случае роль такого аналитического выражения выполняет соотношение (11).

Решение системы (15), (16) проводим методом последовательных приближений, согласно которому

фт+1 = Рп(, фР.... фМ), т = 2,3...(N -1), р = 1,2,3..., (17)

где вектор-функция Гт задается системой (15), (16), ар - номер итерации.

При проведении конкретных расчетов, зависимости плотности тока от потенциала (поляризационные зависимости) были взяты для морской воды при температурах 20-200°С

[4].

На каждой поляризационной зависимости отмечена точка перегиба (0, ф*), ниже которой она аппроксимируется степенной функцией

х

п-1

1

х

х

х

х

т-1

т

п

к

х

N

п=2

] = а

а выше точки перегиба - экспоненциальной функцией

] = Ь0 ехр

{ * \ ф-ф

Коэффициенты этих функций представлены в таблице.

(18)

(19)

/, °с а0, А/м2 а1 Ь0, А/м2 1П м аэ у0, А/м2 Ф* мВ

20 0,013 0,57 0,024 55,65 0,20 138

60 0,0195 0,71 0,037 45,23 0,35 74

120 0,084 0,55 0,070 46,40 0,95 81

200 0,314 0,36 0,097 52,43 1,60 92

Указанные поляризационные зависимости являются нелинейными граничными условиями в выражениях (3), (4), (7), (9), (11), (13), (15).

Проводился расчет электрического поля на границе при температуре 2С°С (А=4,32 Ом-1-м-1), 6С°С (А=7,6С Ом-1-м-1), 12С°С (А=12,52 Ом-1-м-1), 2СС°С (А=19,1С Ом-1-м-1), для плоских и цилиндрических поверхностей. Геометрические параметры расчета были следующими: ширина анода - 0,2 м; ширина защищаемой поверхности между анодами - (1,8-5,8) м; диаметр цилиндрической поверхности - (0,6-2) м. Около 4С вариантов расчета показали, что во всех случаях отмечена приемлемая скорость сходимости (5-2СС итераций), причем в трети случаев сходимость итерационного процесса и единственность решения может быть доказана теоретически еще до начала вычислений, поскольку в этих случаях можно доказать, что вектор-функция ¥т является сжимающим оператором на метрическом пространстве п-мерных векторов {фп} и

выполняется достаточное условие сходимости

N

тах тах £

ф т п=1

с^т

дфг

< 1.

(2С)

При помощи последнего выражения можно объяснить, почему в данной методике сходимость лучше, чем в стандартных. Производные д¥т/дфп при прочих равных условиях будут тем меньше, чем меньше коэффициенты Лтп. Если в данной методике

то в стандартных - по закону к а, где

члены ряда в формуле (14) убывают по закону к -а<2 [1].

Следовательно, Лтп в предложенной методике всегда меньше, чем при классическом подходе.

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

ЛИТЕРАТУРА

1. Гнусин Н.П. Основы теории расчета и моделирования электрических полей в электролитах / Н.П. Гнусин, Н.П. Поддубный, А.И. Маслий. Новосибирск: Наука, 1972. 276 с.

2. Иоссель Ю.А. Математические методы расчета электрохимической коррозиии защиты металлов / Ю.Я. Иоссель, Г.Э. Кленов. М.: Металлургия,1984. 272 с.

3. Янке Е. Специальные функции (формулы, графики, таблицы) / Е. Янке, Ф. Эмде, Ф. Леш. М.: Физматгиз, 1968. 344 с.

с

4. Минаев Е. Н. Контроль коррозии морской техники на основе анализа электрических параметров и полей в водных средах: дис. ... доктора техн. наук / Е.Н. Минаев. Владивосток, 1997. 364 с.

Минаев Евгений Николаевич -

доктор технических наук, профессор кафедры «Общая физика»

Саратовского государственного технического университета

Пулин Виктор Федотович -

кандидат физико-математических наук, доцент кафедры «Общая физика»

Саратовского государственного технического университета

Беляев Илья Викторович -

ассистент кафедры «Общая физика»

Саратовского государственного технического университета.

Статья поступила в редакцию 15.10.07, принята к опубликованию 15.01.08