АНАЛИТИЧЕСКИЕ И ЧИСЛЕННЫЕ РЕШЕНИЯ ДВУМЕРНЫХ ЗАДАЧ ТЕПЛОПРОВОДНОСТИ И ЭЛЕКТРОННОЙ ОПТИКИ Текст научной статьи по специальности «Математика»

УДК 519.63, 51-73 Вестник СПбГУ. Прикладная математика. Информатика... 2021. Т. 17. Вып. 4 МвС 65М12, 35Л05

Аналитические и численные решения двумерных задач теплопроводности и электронной оптики*

Г. И. Курбатова, Е. М. Виноградова

Санкт-Петербургский государственный университет, Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7^9

Для цитирования: Курбатова Г. И., Виноградова Е. М. Аналитические и численные решения двумерных задач теплопроводности и электронной оптики // Вестник Санкт-Петербургского университета. Прикладная математика. Информатика. Процессы управления. 2021. Т. 17. Вып. 4. С. 345-352. https://doi.org/10.21638/11701/spbul0.2021.403

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

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

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

В работе рассмотрено решение двумерных параболических уравнений в кольцевых областях при отсутствии осевой симметрии. Эти уравнения имеют большое значение в моделировании морских газопроводов, систем электронной оптики, упруго-деформированного состояния тел. Численное решение таких задач получено наиболее популярными аддитивными методами: 1) по локально-одномерной схеме Самарского; 2) по схеме переменных направлений Писмена — Рекфорда; 3) по схеме стабилизирующей поправки Дугласа — Рекфорда; 4) по семейству схем расщепления Яненко.

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 20-07-01086).

© Санкт-Петербургский государственный университет, 2021

Решения задачи Дирихле для стационарного уравнения теплопроводности в полярных координатах для кругового кольца в ряде частных случаев известны. В [2] указанные четыре метода тестировались на следующей задаче:

1 д { дТ\ 1 д /дТ\ _ г дг \ дг ) г2 д(р у д(р) '

г € [Дь 2Дх], р € [0, 2п], Т(ДЬр) = Т1 + (Т2 - Т1)(1+со82р), р € [0, 2п], (2)

Т(2Дьр) = Т2 + (Т2 - Т1)зт2р, р € [0, 2п], (3)

имеющей простое аналитическое решение:

Т(г, ф) = (Т2 - Тг) (§ + - щ) соз2^ + 2Т, - Т2.

Однако граничные условия (2), (3) редко встречаются в практических задачах. В настоящей работе представлено новое аналитическое решение краевой задачи для уравнения Лапласа (1) в кольцевой области при граничных условиях

г = Дь р € [0, 2п] : Т(Дь р, г) = То, (4)

г = Я.2, р € (-ро,ро): Т(К2,р,1)= ТО, (5)

г = Й2, р € (ро, 2п - ро) : Т(Н.2,р,г)= ТО, 0 < ро < п.

Граничные условия (5) позволяют охватить широкий класс задач электронной оптики и установившегося теплообмена морского газопровода при разных вариантах его заглубления в грунт [2, 3].

2. Постановка нестационарной задачи. Рассмотрим выход на установившийся режим распределения температуры Т = Т(г, р, Ь), моделируемого уравнением теплопроводности

дТ а д ( дТ\ а д2Т

г— I +-9-^—9, ¿>¿0, (6)

дЬ г дг \ дг ) г2 др2 '

в кольцевой области О = {(г,р), г € [Д1, Д2), р € [0, 2п]} при краевых условиях (4),

р

1>1о, г € [Й1,Й2), р € [0, 2п] : Т(г,р,Ь) = Т(г,р + 2п,Ь) (7)

и при постоянном начальном условии

Т(г, р, Ьо) = Ть. (8)

Величины То, Т° Т2о, Ть считаются постоянными, а — коэффициент температуропроводности материала. Приведем алгоритмы численного решения нестационарной задачи (4)-(8) аддитивными методами.

3. Локально-одномерная схема Самарского. Эта абсолютно устойчивая разностная схема является аддитивной и экономичной, она допускает обобщение на трехмерные уравнения теплопроводности. Оператор Лапласа в уравнении (6) представляется в виде суммы двух операторов Ьг, Ь^:

— = а(Ьг(Т) + Ь^Т)),

дТ

1 д ( д\ т 1 д2 Ьг = - — \ г— , Ьи

= аЬр(Т) (10)

дг у дг) ' ^ г2 д(р2

Переход на следующий временной слой выполняется в два этапа. Первый этап:

дТ

— =аЪг(Т) (9)

с начальным условием (8) и граничными условиями (4), (5). Второй этап:

дТ

с условием периодичности (7). В качестве начального распределения температуры на этом этапе используются значения температуры, полученные из решения соответствующих задач первого этапа.

Вводится равномерная пространственная сетка ш = {(гг,р^),гг = Я1 + гНг, р^ = зНр; г = 0,... ,МГ; 3 = 0,..., Мр} с шагами Нг, Нр то перемениым г и р соответственно (Нр = 2п/Мр, Нг = (Я2 — Я1)/Щ^, {Ьп = пт, п = 0,1,... Щ} — равномерная сетка по времени с шагом т. Каждый промежуток времени \Ьп,1п+{\ разбивается моментом времени 1п+1/2 = Ьп + т/2 та две равные части; ТП1+1/2 — температура в узле

('г ,Рз ,^п+1/2)-

На первом этапе для каждого 3 одномерное уравнение теплопроводности (9) аппроксимируется по следующей неявной монотонной разностной схеме:

Тп+1/2 Тп (, , ч Тп+1/2 _ Тп+1/2 , . Тп+1/2 _ Тп+1/2 \

±ЬЗ 1ЬЗ_а I ( , -4+1,3 Чз г \ 1Ъ,3 1г-1,3 \

^ ц ц )'{ }

г = 1,...,щг — 1.

Выражения (11) записываются в виде трехточечных разностных уравнений и дополняются условиями па внутренней г = 0 и внешней г = Щг границах области О:

Т0п+1/2 = То, (12)

рп+1/2 I Т0, 3 =0,..., (по — 1) и 3 = (Щр — по + 1),..., Щр, по = ро/Нр =1 То, з = (по + 1),...,Мр — (по + 1).

Получившаяся система с трехдиагональной матрицей решается методом прогонки.

г

аппроксимируется по неявной разностной схеме

Тп+1 _ тп+1/2 тп+1 _ 2Тп+1 + Тп+1

Ь3 =* ЬЗ+1~ ^ ^ 1 = 0,(14)

1 'г Нр

Уравнения (14) совместно с разностным аналогом условий периодичности (7)

Т^ = Т^1 (15)

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

грп+1 грп+1 . грп+1 _ т-1

— Тг,о + Тг,1 = —рг,о, Вестник СПбГУ. Прикладная математика. Информатика... 2021. Т. 17. Вып. 4 347

3 - Сьо ТП+1 + 3 = , з = - 1, (16)

грп+1 _ ^ ТЛП+1 | тлп+1 _ _ р

V 1

Г

Уг,3

2 + К = Т,

^п+1/2

Матрица этой системы не имеет трехдиагонального вида. Чтобы ее решить, для каждого г использовался метод циклической прогонки, который, как известно [4], является экономичнее итерационных.

4. Метод переменных направлений Писмена — Рекфорда. Он широко применяется для решения двумерных параболических уравнений. Решение задачи (4)-(8) этим методом также проводится в два этапа.

На первом этапе уравнение (6) аппроксимируется неявной разностной схемой по переменной г и явной — то переменной у.

Тп+1/2 _тп

т Ь3 = 0.5 (лг (т^2) +А, (Т™)). (17)

Г

ная — по у.

Т"+1 _ тп+1/2

'" " - = 0.5 (Лг (Т^2) + Л, (Т^1)). (18)

',3 / V \ 1,3

1 ''2_т п

т

Разностные операторы Лг, А^ записываются следующим образом:

Лг [тп+1/2)

■ п+1/2_ грп + 1/1

, г = 0,

ь.1

, п + 1/2_ т п+1/2 , Ч т п+1 /2_ т п + 1/2ч

грп _ г\грп I грп

и' 2 1,2 1г

Уравнения (17) дополняются условиями (8), (12), (13), а уравнения (18) — условиями периодичности (15); в качестве начального распределения на втором этапе применяются значения температуры, полученные из решения задач на первом этапе. На первом этапе уравнения (17) приводятся к трехточечным разностным, кото-

з

г

ме типа (16) с соответствующими коэффициентами Рг,з- Решение этой системы проводилось с помощью циклической прогонки.

Ввиду ограниченности объема статьи алгоритмы метода стабилизирующей поправки Дугласа — Рекфорда и метода расщепления Яненко не приведены, представлены только результаты расчетов относительной погрешности численного решения задачи (4)-(8) на основе полученного аналитического решения.

5. Аналитическое решение стационарного варианта задачи (4)—(8). Трудность аналитического решения этой задачи заключается в наличии разрыва температуры в граничной точке (г = Д2, у = уо)- Построенное аналитическое решение позволяет удовлетворить граничному условию с точкой разрыва. Воспользуемся подходом, предложенном в работе [5]. При граничных условиях (4), (5) распределение

Г2 к2

ат

ат

а

а

температуры симметрично относительно лучей р = пк (к = 0,1), поэтому достаточно найти решение задачи в секторе р € [0, п]. Для решения краевой задачи (1), (4), (5) внутренняя область (г € [Е1,Е2), р € [0,п]) разбивается на такие подобласти: I — р € [0, ро); г € [Й1, Е2); II — р € (ро, п\, г € [Дь Е2). В каждой из них распределение температуры представляется следующим образом:

Функции Т\(г, р), Тц(г,р) являются решениями уравнения Лапласа (1) при однородных граничных условиях в областях I и II:

ТГ1(Е1,р)= ТЕ ,р)=0, Ти(Е1,р)= Ти(Е2,р)=0. (20)

Краевые задачи (20) для функций ТГ1(г,р), ТГ11(г,р) решаются методом разделения переменных. Опуская выкладки, приведем эти решения:

л, л ^ совЬ(Лпу) . г

Щг, р)= Т, ап-Т-ГТ-Г 8т Лп 1п —

П=1 С08П(Л„ро) V Д1

г < \ ^Тч. С0^(Лп(п - р)) . / г \ пп

Тц{г,р)=^Ъп-——-— эш I \п 1п —— 1 , А» =

С08ЦЛ„(п - ро)) \ Д1 ) ЩЕ2/Е1)

(21)

Коэффициенты ап, Ьп в рядах (21) находятся из условий непрерывности функций Т1(г, р), Т11(г, р) и их производных по переменной ^и р = ро, г € [Д1, Е2), которые имеют вид

^ 2(-1)п+1 ЫиЪ(Лп(п - ро))

Яп = _ у11) у '__^ гул__/22\

™ 2 1 7гп 1ап11(А„с(г>о) + 1апЬ(Л„(7г — ро))'

Ъп = (Т° - г")2^1)"_1апЬ(Апу0)__з

™ 2 1 пп 1апЬ(Л„((г'о) + 1апЬ(Л„(7г — ро))

Формулы (19), (21)-(23) представляют собой точное аналитическое решение краевой задачи (1), (4), (5).

6. Оценка достаточного числа членов в бесконечных рядах (21), обеспечивающего заданную точность аналитического решения. Рассмотрим, например, подобласть I. Функция Т1(г,р) (21) представлена бесконечным рядом, который имеет, как нетрудно показать, экспоненциальную сходимость. Запишем функцию Т1(г,р) (21) так:

N. со

пп п=1 п=N+1

совЦА^) / 1п Л (24)

пп 1апИ(Лпро)+ 1апИ(Лп (п - ро)) собЩЛп ро) \ Е1)

Воспользуемся мажорирующим рядом и оценим число членов ряда, начиная с ко-

торого сумма бесконечного ряда I ^ Бп I будет меньше заданной малой величи-

V n=N^ + 1 )

ны £о (в К), Опуская выкладки, приведем найденную оценку Ы*(£о) для подобласти I: Вестник СПбГУ. Прикладная математика. Информатика... 2021. Т. 17. Вып. 4 349

N.(eo)>U In . ° Д e =

в\ £o(exp (в) - 1) J П ln(R2/Ri)

N

Таким образом, функция T\(r, ф, N*) = Y1 Sn с точностью до е0 является решением

n=1

уравнения (1) с граничными условиями T\(Ri, ф) = T\(R.2, ф) = 0 в подобласти I. При ф = ф0 = п/2 выражение для Sn (24) упрощается и бесконечный ряд для T\(r, ф) (21) суммируется, функция T\(r,n/2) имеет вид

Если ф0 = п/2, то та луче ф = ф0, исключая особую точку r = R2, оценка N*(e0) усложняется, но не представляет принципиальной трудности.

Аналогичным образом получаются оценки для функции Тц(г, ф) (21) в подобласти II.

7. Результаты расчетов. Для каждого набора параметров задачи (4)-(8) рассчитывался момент времени t*, начиная с которого с заданной точностью е устанавливается стационарное распределение температуры Тф.

Для сравнения эффективности рассматриваемых численных методов вычислялась величина

— средняя по области G относительная погрешность расчета температуры по k-му методу для момента времени, большего момента t*:

- _!_ Y у /E^bSl) (25)

В (25) Ta(ri, pj) — температура, определенная по найденному аналитическому решению (21)-(23), — температура в узле (i, j) в момент времени tn > t*, рассчи-k

hr, hp, т (при параметре а = 0.9 k = 1 соответствует расчету по локально-одномерной k=2

k = 3 — схеме стабилизирующей поправки Дугласа — Рекфорда, k = 4 — по семейству схем Япепко).

В качестве иллюстрации расчетов величин 5(к) представим результаты вычислений при следующем варианте параметров задачи:

R1 = 0.5 м, R2 = 0.54 м, ф0 = п/2, Ть = 283.15 К, T0 = 293.15 К, T0 = 278.15 К,

T20 = 283.15 К, a = 5.33 • 10"6 м2/с, Np = 240, Nr =41, т = 0.001 с. (26) t*

вившийся режим с точностью до сотых градуса t* « 9.6 мин. Количество N* членов ряда, обеспечивающее равномерно по всем r и ф из расчетной области точность аналитического решения порядка 10"8 N* = 19. Приведем результаты расчетов средней относительной погрешности 5(к) указанными в п. 1 методами для варианта параметров (26):

¿(1) = 7.71 • 10"4; 5(2) = 6.57 • 10"4; 5(3) = 7.75 • 10"4; ¿(4) = 7.83 • 10~4.

7. Заключение. Проведенные расчеты по локально-одномерной схеме Самарского, по схеме переменных направлений Писмена — Рекфорда, по схеме стабилизирующей поправки Дугласа — Рекфорда и по семейству схем Яненко двумерного уравнения теплопроводности в кольце с неоднородными граничными условиями по p для разных вариантов параметров подтвердили вывод, сделанный в работе [2]: наиболее эффективным для решения таких задач является метод Писмена — Рекфорда. При одних и тех же значениях шагов он точнее в расчетах временных характеристик и не уступает другим рассмотренным методам в точности расчета пространственного распределения температуры. Однако обобщение метода Писмена — Рекфорда на трехмерные задачи сопряжено с большими трудностями [1]. Этого недостатка лишена локально-одномерная схема Самарского, точность расчета пространственного распределения температуры которой не меньше, чем у метода Писмена — Рекфорда.

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

Литература

1. Самарский А. А., Вабищевич П. Н. Аддитивные схемы для задач математической физики. М.: Наука, 2001. 319 с.

2. Kurbatova G. /., Klemeshev V. A. Effective numerical methods for calculating non-stationary heat and glaciation dynamic processes for offshore gas pipelines // Energy. 2020. Vol. 205. N 117995 (7 p).

3. Курбатова Г. И., Ермолаева Н. Н., Филиппов В. В., Филиппов К. Б. Проектирование газопроводов в северных морях. СПб.: Лань, 2020. 352 с.

4. Самарский А. А., Николаев Е. С. Методы решения сеточных уравнений. М.: Наука, 1978. 561 с.

5. Vinogradovа Е. М., Egorov N. V., Starikova А. V., Varayun М. I. Calculating a multipole cylindrical electrostatic system // Technical Physics. 2017. Vol. 87. Iss. 5. P. 791-794.

Статья поступила в редакцию 6 августа 2021 г. Статья принята к печати 13 октября 2021 г.

Контактная информация:

Курбатова Галина Ибрагимовна — д-р физ.-мат. наук, проф.; [email protected] Виноградова Екатерина Михайловна — д-р физ.-мат. наук, проф.; [email protected]

Analytical and numerical solutions of two-dimensional problems of heat conduction and electronic optics*

G. I. Kurbatova, E. M. Vinogradova

St. Petersburg State University, 7-9, IJniversitetskaya nab., St. Petersburg, 199034, Russian Federation

For citation: Kurbatova G. I., Vinogradova E. M. Analytical and numerical solutions of two-dimensional problems of heat conduction and electronic optics. Vestnik of Saint Peters-bury University. Applied Mathematics. Computer Science. Control Processes, 2021, vol. 17, iss. 4, pp. 345-352. https://doi.org/10.21638/11701/spbul0.2021.403 (In Russian)

This work is devoted to a comparative analysis of the effectiveness of various methods for the numerical solution of a non-stationary two-dimensional problem of heat conduction in

* This work was supported by the Russian Foundation for Basic Research (grant N 20-07-01086).

a cylindrical coordinate system in the absence of axial symmetry. The analysis is carried out on the basis of the obtained analytical solution of the two-dimensional problem of heat conduction in a steady state. The results of this work are of great applied importance in calculating the processes of gas transportation through offshore pipelines, as well as in calculating the electrostatic field of focusing systems during transportation of charged particle beams in problems of electronic optics.

Keywords: offshore gas pipelines, frontal interval method, heat exchange model, numerical method with increased accuracy, electronic optics.

References

1. Samarskii A. A., Vabishchevich P. N. Additivnyie skhemy dlia zadach matematicheskoy fiziki [Additive schemes for problems of mathematical physics]. Moscow, Nauka Pabl., 2001, 319 p. (In Russian)

2. Kurbatova G. I., Klemeshev V. A. Effective numerical methods for calculating non-stationary heat and glaciation dynamic processes for offshore gas pipelines. Energy, 2020, vol. 205, no. 117995 (7 p).

3. Kurbatova G. I., Yermolayeva N. N., Filippov V. B., Filippov K. B. Proyektirovaniye gazoprovodov v severnykh moryakh [Design of gas pipelines in the northern seas]. St. Petersburg, Lan' Publ., 2020, 352 p. (In Russian)

4. Samarskii A. A., Nikolaev E. S. Metody resheniya setochnykh uravneniy [Numerical methods for grid equations}. Moscow, Nauka Publ., 1978, 561 p. (In Russian)

5. Vinogradova E. M., Egorov N. V., Starikova A. V., Varayun M. I. Calculating a multipole cylindrical electrostatic system. Technical Physics, 2017, vol. 87, iss. 5, pp. 791-794.

Received: August 06, 2021. Accepted: October 13, 2021.

Authors' information:

Galina I. Kurbatova — Dr. Sci. in Physics and Mathematics, Professor; [email protected]

Ekaterina M. Vinogradova — Dr. Sci. in Physics and Mathematics, Professor; [email protected]