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

УДК 519.634

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

В. П. Федотов, О. А. Нефедова

Институт машиноведения УрО РАН,

620219, Екатеринбург, ул. Комсомольская, 34.

E-mails: f edotov_vp@mail .ru; nefedova@imach.uran.ru

Предложен алгоритм для нахождения численно-аналитического решения задач параболического типа (диффузии и теплопроводности). В рамках алгоритма решение задачи осуществляется в три этапа. На первом этапе решается одномерная задача для базового отрезка, которая имеет самостоятельное значение, и в тоже время служит основой для второго этапа. На втором этапе 'рассматривается двумерная задача, ее решение выполняется с использованием модифицированного метода граничных элементов. На третьем этапе применяется приём пошагового интегрирования во времени.

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

Введение. Разработан численно-аналитический алгоритм решения задач параболического типа, в основе которого лежит модифицированный метод граничных элементов (МГЭ). В алгоритме также предложена оригинальная аппроксимация неизвестных граничных функций, когда при решении двумерной задачи в качестве аппроксимирующих функций используются решения одномерной задачи, выраженные через узловые граничные значения и функции влияния. В целях сокращения времени вычислений одномерная задача решается для удобного базового отрезка. Такая аппроксимация, включающая в себя два вложенных друг в друга алгоритма, позволяет получить для двумерной области аналитическое решение, наиболее близкое к действительному.

1. Постановка задачи и её слабая формулировка. В работе рассматриваются однородные нестационарные дифференциальные уравнения параболического типа для задач диффузии и теплопроводности при заданных начальных и граничных условиях:

д Р0М)-±^=°, .6 11;

p(x,to) = Ро(х), xeQ] (1)

p(x0,t) = p*(x0,t), х0 е Si] u(Xo,t) = U*(xo,t), Xq € S2-

Здесь A = d2/dx\ + d2fdx\ — оператор Лапласа; p(x, t) —температура или концентрация примеси в точке х(х\, Х2) в момент времени t ^ to] to — начальный момент времени; Л — не зависящий от координат и времени физический

Владимир Петрович Федотов (д.т.н., проф.), главный научный сотрудник, лаб. прикладной механики. Ольга Анатольевна Нефедова, младший научный сотрудник, лаб. прикладной механики.

параметр; О — исследуемая область; 5 = — гладкая граница области О;

Ро(х), р*(хо,£) и и*(хо, £) — известные функции; и(ж,£) = —<9р(ж, £)/<9п(ж) — тепловой или диффузионный поток через границу с внешней нормалью п(х). Потребуем выполнения (1) в целом по области П с весом О и учтём изменение температуры (концентрации примеси) во времени:

[ [ ~Ар(х, £) - X, гк^)(1^{х)(И = О, (2)

Л дЬ

где Ьк — момент наблюдения; £(£1,^2)—произвольная точка области. В соответствии с методом граничных элементов [1] в качестве весовой функции С(£, ж, выберем фундаментальное решение исходного уравнения. Фун-

даментальное решение С(^,х,1к^) описывает реакцию в точке ж в момент времени Ьк на действие единичного точечного источника вида дельта-функции Дирака 5(^,х,1к^), помещенного в точку £ неограниченной области в момент времени £, и определяется из уравнения

ДС(£, ж, гк, £) - = §((,х, гк^). (3)

Решение уравнения (3) известно:

Здесь Т = гК - Г2 = ЗД; ^ = Ж; - &, * = 1,2.

Направленный поток , обусловленный С*, задаётся выражением

те Ж* *) = -^ = -^п- = ^ = -^-ехрГ--^

; дхг г 2Лг 8тт\2т2 V

где Пг{ж)—компоненты единичного вектора внешней нормали к линейному элементу, проходящему через точку Ж(Ж1,Ж2); (I = (х\ — {1)^1 + (Ж2 — £2)^2-В общем случае т-мерного пространства зависящее от времени фундаментальное решение имеет вид

а'^Лк.О = (4ТДт/2ехр(-^)н(г)’

где Н (г) — функция Хевисайда.

Проинтегрируем выражение (2) по частям дважды по ж, а затем выполним интегрирование по £. Полученное выражение для температуры (концентрации примеси) р(£, Ьк) в произвольный момент времени Ьк в произвольной внутренней точке £ области П имеет вид

р(£^к) = — А [ [ и(х,£)С*(£,х,1к^)с13(х)сИ+

Ло

+ А / / р(ж, £)_Р*({, ж, ^с13(х)сИ + / Ро{х)С*{^х^к^о)сЮ,{х). (4)

./ £о П

Граничное интегральное уравнение для задачи (1) запишется так:

1 ftK Г

-р(хо, tx) = — А / / u{x,t)G*{xo,x,tK,t)dS{x)dt+

2 Jto Js

2. Одномерный случай. В качестве области Q рассмотрим базовый отрезок длины L, который лежит на оси абсцисс и один из концов которого совпадает с началом координат. Тогда граница S области Q вырождается до двух точек: 0 и L. Перепишем интегральное уравнение (4) в предположении, что Ро = const:

Для получения аналитического решения уравнения (5) используем алгоритм, предложенный в [2]. Рассмотрим плоскость {x,t} и одномерную однородную область, простирающуюся от х = 0 до х = L, с заданными вдоль неё начальными условиями р$ = const. Другими границами являются прямые, параллельные оси времени, вдоль которых заданы две постоянные граничные величины. Принимая время t в качестве параметра, две неизвестные граничные величины также предполагаем постоянными во времени. Соотношение (5) запишется в виде

где «1 = и(0+), и2 = и(Ь~), рх = р(0+), Р2 = р(Ь~).

Интегралы по времени и пространственной переменной в правой части выражения (6) были вычислены аналитически, затем точка £ поочередно устремлялась к точкам 0 и Ь. В результате была получена система двух уравнений для нахождения двух неизвестных граничных значений.

Функциональная зависимость р (£, Ьк) для произвольной точки отрезка [О, Ь] принимает вид

+ Ро [ G*(£,x,tK,to)dx. (5)

+ Ро [ G*(i,x,tK,h)dx, (6)

p(i,tK)=-u2

-(L-Q2

4X(tK-to)

+

+ U\

X(tR — to)

expf ^---------) - ^erfcf— ^

\4X(tK-t0)J 2 \2v/\(tK-t0)

IT

+

+ 'Pl.rU (L~^) \ + £W

+ 2 \2^X(tK — to)' ^ V

2д/Л(tx ~ to)^

+

i * + 2»

erf

(£-0

+

erf (-

2 л/Щк - to)' ^2 л/Щк - to)

(7)

ry

Здесь erfc(y) = 1—erf (у), erf (у) = (2/у/ж) / exp (—-г2) dz — функция ошибок.

Jo

Для нахождения функциональной зависимости и(£, tx) используются формула для определения потока u(x,t) = —dp(x,t)/dn(x) и свойство симметрии функции G* относительно х и

и(.£,Ък) = А

'£о

«(ж)

dG*(£,x,tK,t)

— р(х)

dF*(£,x,tK,t)

дС

rL

dt—

о

Ро

dG*(i,x,tK,to)

■)о д(

или, после аналитического интегрирования правой части выражения (8

dx, (8)

u(£,tK) = —erfc (—^

^ K> 2 \2y/X(tK-t0)

, Ul f + —erfc

£

2л/Щк - to) ex Pl (

2л/тт\(tx — to) \4A(ix —^o)/ 2y/nX(tK — to) \4A(ix — £o)/

+

2\/irX(tK ~ to)

exp -

(L-Q2 4X(tK ~ to)

- exp -

£2

4Л(tK ~ to)

+

(9)

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

3. Двумерный случай и модифицированный метод граничных элементов.

Разобьём границу 5 = 51 и 62 двумерной области О на конечное число отрезков а/\, ] = 1, 2,... , N + М, в предположении, что на одной части

границы 51 = {[aj-l, а/\, ] = 1,2,... , N} задана температура или концентрация примеси, а на другой части 5г = {[%-1, %], 3 = N + 1, N + 2,..., N + М} задан тепловой или диффузионный поток. Тогда интегральное уравнение (4) перепишется в виде

rtK N а.

p(t,tK) = а / 1

j_ 1 Jaj-i

'to

—и^\ж, t)G*{£, ж, tx, t)+

+ р^*(ж, t)F*(£, ж, tx,t))dxdt+

N-\-M ('O’j

+ A / ^ / (—u^')*(x,t)G*(£i,x,tK,t)+p^\x,t)F*(£i,x,tK,t))dxdt+

j=iV+l “j- 1

+ / Рв(х)С*((,х,гк,Ь)с1£}(х). (10)

Здесь звездочкой обозначены известные величины и введены следующие обозначения: и^\х,Ь) = и(х,і), р(її*(х,і) = р*(х,і) при х Є [aj-l,aj\ Є 5і; и^*(ж,і) = -и*(ж,і), р(її(ж,і) = р(ж,і) при ж Є [aj-l,aj] Є 5г-

Интеграл по О в правой части уравнения (10) преобразуем с помощью второго тождества Грина к граничному интегралу

1^р*0(х)С*((,х,ік,іо)(ІП(х) = ^~^|^ехр(~4Л(^ _^0))Р°^ +

1 „ / г2

+ 2Е\тк-и)Г(х)Г(х)' (и)

ЧТО ВОЗМОЖНО ДЛЯ гармонической функции Ро(х). Здесь

/У

(ехр (г)/г)йг, и$(х) = — др$(х) / дп(х).

-ОС

Для аппроксимации неизвестных граничных условий в уравнении (10) используются функциональные зависимости (7) и (9).

Модифицированный метод граничных элементов [3] позволяет заменить произвольное граничное разбиение \dj-i, а?] на базовый отрезок [0, Ь\. Выполним преобразование координат для разбиения [aj-l,aj\, при котором длина отрезка не изменяется, точка (lj-l отображается в начало координат, а точка отображается в точку Ь(Ь, 0). Такое преобразование является комбинацией параллельного переноса и поворота отрезка на угол ср (угол наклона отрезка к оси абсцисс). Тогда произвольная точка на плоскости ж(ж1, Жг) отображается в точку ж(ж1,ж2), связанную с ней соотношением

ж = Бж + С, х = В~1(х — С),

где В — матрица поворота; С — вектор координат той точки граничного отрезка, которая при преобразовании отображается в начало координат. При таком преобразовании исследуемая система объектов жестко перемещается как единое целое. Поскольку величина р — скалярное поле в задаче (1), то справедливы равенства

го,, гЬ ^

/ и(х, £)С*(£, ж, Ьк, £)(1х = / и(х,1)0*(^,х,1к,^с1х,

7»,1 -)°[ь_ _ (12)

/ р(х, Ь)Р*(£, ж, Ьк, £)(1х = / р(х,1)Р*(^,х,1к,^с1х.

«/ а ^ ]_ 0

Здесь { = В-1^ ~ С).

Запишем вид функции влияния и ее нормальной производной для базового отрезка [0,£], на котором Ж2 = 0, П1(ж) = 0, Пг(ж) = —1:

С*(^,ж,^,^) = тД-ехрГ-^Ж ^ +^2

АжХт \ 4 Ат /

^*(С,ж,ік,і) = 0 ^2 9ехр(-(Х ^ +^2

8тгА2т2 4 Ат

Подставим выражения (7), (9), (13) в уравнение (10) и учтём равенства (11) и (12). В результате получим

ҐІК ^ - — — —

Р(&К) = / ^2\~Сі(иЗ-ІІІС) +иЗІ2)) -С2(^*_і4?) -Р^і4)) +

*° 3 = 1

/•*ЛГ м+м

+ 2сзР*_14° М + / Е [Сз + ^■/?)) + сз {и*з-1 (^4° - 4°) “

*° І=ЛГ+1

-и*(С44?) -4°)) -2сіи*_14?)

(і£~Ь

г N _ _ дг+м

+ /

Ло

(4й - 4") + £

■ ' ^=м+1

лг+м _

+ ^ Е (ро,--14° + 2<-14°)- (!4)

3 = 1

Здесь и*, р./, р* — узловые значения соответствующих функций; звёздочкой, как и прежде, обозначены заданные величины;

1 сі о /М^-^о)

Сі =----;-------С2 = —. , Сз =-к С4 = 2\ -:

8тт (Ьк ~ £) л/ттМТ^То) 1б7гЛ (^ - І) V 7Г

/<« = Л-Ч, ..У-гГ (Ж.Т.!1)2

/о 4 2 д/ Л (і - £ог V 4 Л(іаг-і)

^ “Л<г^(ГйУ

IР = Гехр(-7т7^)ехр('-^#£±Й )<Ь,

3 Уо V 4А(і — іо) / 4Л(^-і) 1

г(0 _ А*'™™/' (1 ~ ж)2 > — /" (^~Сі)2+С2

= 7о ехЧ-ЩГ^) і ехЧ~ 4А(І*-І) 1Аг-

1р=£ Н~%^)'2Ух-

(О - (^-Єі)2+Й

^ =/„ таг1Н2^л«-го)^хрг 4л«У-о г-

4° = [{1~ Х)ак (2л/лс7-І0)) еХР (~ (^к -+<)2 ) ёХ'

1{Р = I ----^--=2 ехр (Л-^2+Н <&>

Л (ж-^З+^з V Щік-іо) )

7(0

19

ГЬГ ((X -її)2 +Ґ2\ІТ і ЛЩік-іо)) ‘

Соответствующее уравнению (14) граничное интегральное уравнение для произвольной граничной точки Жо(жо1, Ж02) имеет такой вид:

1 ft К ^

-p(x0,tK) = / J2[~Cl{U3-lJ<iXO) +UiJ2°]) - C2{p*j-lAX0) -P*jjt°]) +

По j=1

— 1 r _ _

+2c3p*_14xo)]dt+ \ф3-^1о)+p3j{2o))+

to j=N+l'

+ C3 [u*_! (c4 jf0) - 4Xo)) - U* (c4 Jf0) - 4X0))] - 2c1U*_14X0)

,tKr N _ _ N+M

dt-h

dt-h

+ Г Ec2^_i(jfo)-jfo))+ £ c3^_1(jfo) + jfo)+2jfo))

4=1 j=N+1

1 W+M _

+ 2^ E (poj-i4:CO'> + 2'u0j-i4:C0'>)> (15) j=l

где Xo = B~l(xо — С). Интегралы Jn ^ = Jn(xo,x,tK,t) в уравнении (15) соответствуют интегралам iffi = In(£,x,tK,t) из уравнения (14) с заменой точки C(Ci,C2) на точку Ж0(Ж01,Ж02).

Уравнение (15) интегрируется по пространственной переменной. Для получения результата интегрирования в аналитическом виде специальная функция ошибок erf (у) была приближена кубическим сплайном.

4. Интегрирование по времени. Рассмотрим заключительный, третий этап. Разобьём интервал интегрирования [to,tK] на К шагов [tk~i,tk\, считая граничные условия постоянными на каждом временном шаге, и будем выполнять интегрирование по времени.

Использование модифицированного МГЭ позволяет ограничиться в конце каждого временного шага подсчётом значений функций лишь в граничных узлах, для того чтобы использовать их при подсчёте значений функции температуры (концентрации примеси) внутри и на границе области Q на следующем шаге по времени.

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

к к

Е +Е RkKuk+ад* = о. (16)

к—1 к— 1

(*tk

гч

Здесь Я*кк-> Щк ~ матрицы коэффициентов / д*(хо,tк,t)dt,

Г* к

/ г*(хо,ЬК, t)dt при к = 1, 2,... , К. Подынтегральные функции д*(хо, £) и г*(ж0,£к,£)—аналитические выражения от различных комбинаций произ-

ведений функции ошибок ег((Хо,Ьк, £) и функции ехр(ж0, Коэффициен-

ты матриц Рк и 11к — узловые значения температуры (концентрации примеси) и теплового (диффузионного) потока соответственно для каждого момента

времени Коэффициенты матрицы обусловлены интегрированием ■1^°0\ когда Жо € [aj-l,aj\ € 62; Р0* — вектор коэффициентов Ро^_1-

Система уравнений (16) решается численно для момента времени ^ = при уже найденных на предыдущем шаге значениях Р^ и II^ при к = 1,2,, К — 1. После решения системы уравнений (16) и определения всех граничных величин в каждый момент времени рассчитывается температура (концентрация примеси) в произвольной внутренней точке £ области П.

5. Пример. В качестве примера применения предложенного алгоритма была рассмотрена задача о распределении температуры в квадратной области размером 3 х 3 м с начальной температурой р$(х) = 0°С и коэффициентом теплопроводности Л = = 0,204 Вт/(м-°С), для которой задано граничное условие р*(жо,£) = 20°С. На рисунке приведены значения температуры, полученные с помощью предложенного метода (сплошная линия), с помощью пакета РгееРет++-сз (штриховая линия) и для аналитического решения [4] (штрих-пунктирная линия). Все данные приведены для точки, расположенной в центре пластины.

Отклонение результатов численных расчётов от аналитического решения вблизи начального момента времени обусловлено тем, что принятый в качестве аналитического решения ряд Фурье даёт на этом интервале времени заведомо неточное приближение, что объясняется видом ряда Фурье, а также разрывом между начальными и граничными условиями. Однако в данной постановке задачи поведение решения на начальном интервале времени является несущественным, поскольку основной интерес представляет развитый процесс теплопроводности. Начиная с некоторого момента времени все три решения задачи стремятся к асимптотическому значению.

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

Работа выполнена при поддержке ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 годы (госконтракт № 02.740.11.0202) и Программы Президиума РАН № 14 «Интеллектуальные информационные технологии, математическое моделирование, системный анализ и автоматизация».

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Бреббия К., Теллес Ж., Вроубел Л. Методы граничных элементов. М.: Мир, 1987. 526 с. [Brebbiya К., Telles Zh., Vroubel L. Boundary element techniques. Moscow: Mir, 1987. 526 pp.]

2. Бенерджи П., Баттерфилд В. Методы граничных элементов в прикладных науках. М.: Мир, 1984. 494 с. [Benerdzhi Р. К., Baiterfild R. The boundary elements methods in applied sciences. Moscow: Mir, 1984. 494 pp.]

3. Федотов В. П., Спевак Л. Ф. Модифицированный метод граничных элементов в задачах механики, теплопроводности и диффузии. Екатеринбург: УрО РАН, 2009. 161 с. [Fedotov V. P., Spevak L. F. A modified boundary element method in problems ol mechanics, heat transler, and diffusion. Ekaterinburg: UrO RAN, 2009. 161 pp.]

4. Лыков А. В. Теория теплопроводности. М.: Высш. шк., 1967. 600 с. [Lykov А. V. Theory ol Heat Conduction. Moscow: Vyssh. shk., 1967. 600 pp.]

Поступила в редакцию 04/XI/2010; в окончательном варианте — 27/IX/2011.

MSC: 65М38

APPLICATION OF THE MODIFIED BOUNDARY ELEMENT METHOD FOR THE SOLUTION OF PARABOLIC PROBLEMS

V. P. Fedotov, O. A. Nefedova

Institute of Teoretical Engineering, Ural Branch of RAS,

34, Komsomolskaya St., Yekaterinburg, 620083, Russia.

E-mails: f edotov_vp@mail .ru; nefedova@imach.uran.ru

An algorithm for finding numerically-analytical solution of parabolic problems (diffusion and heat conduction) is proposed. The problem is solved by the proposed algorithm in three steps. At the first step the one-dimensional problem is solved for a base interval of integration. This problem is of independent significance as well as the basis for the second step. At the second step the two-dimensional parabolic problem is considered.

Its solution is performed using the modified boundary elements method. At the third step, the method of step-by-step integration over time is used.

Key words: linear parabolic equation, modified boundary element method, numerical-analytical integration.

Original article submitted 04/XI/2010; revision submitted 27/IX/2011.

Vladimir P. Fedotov (Dr. Sci. (Techn.)), Chiel Research Scientist, Lab. ol Applied Mechanics. Olga A. Nefedova, Associate Researcher, Lab. ol Applied Mechanics.