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

УДК 5Ї9.6З

О ВЫЧИСЛЕНИИ ТЕМПЕРАТУРНЫХ И ДИФФУЗИОННЫХ ПОЛЕЙ В КУСОЧНО-ПОСТОЯННЫХ АНИЗОТРОПНЫХ СРЕДАХ

© В. Н. Кризский*, А. Р. Бикбаева

Башкирский государственный университет, Стерлитамакский филиал Россия, Республика Башкортостан, 433103 г. Стерлитамак, пр. Ленина, 49.

Тел./факс: +7 (3473) 43 23 SO.

E-mail: [email protected]

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

Ключевые слова: анизотропная среда, температурное поле, диффузионное поле, способ вычисления.

Введение

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

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

Постановка задачи и способ решения

Пусть область исследования О С В может

N

быть представлена в виде О = ^ Ог - объедине-

г=1

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

м

г = и г - внешнюю границу области О, а че-] =1

рез и Гу, соответственно, внутренние и внешние границы областей О..

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

div(Jt ■VUi(P,t))-- aiUi(P,t) - bi

2 dUi (P,t)

dt

= - f (P,t), P є О. , i = IN; Ut (P, t) - Uj (P, t)| ^ =

= О, (j. VU. (P, t), n)-

-(jjVU j(P ^ n)r,nrj = О

. =1 N, j є J =

= {j|j = 1, i -1; Y n Yj * О} aj (P)(jj VUj (P, t), n)-

-pj (P)U j (P, t)| ,єг. = Wj (P, t), \at(P)| + |в(P)| *О,i =

(1)

(2)

(3)

(4)

(3)

] 4 '| у]

= [1,г2,...,гк, к < N; и т (Р, г) ^ О, Р ^ тс, т =

= т1,т2,...,тп, п < N; и. (Р,О) = О, г = .

Здесь fi (Р, г)- функции интенсивности источников/стоков поля в подобластях О.; аг, Ьг -постоянные в о. числовые коэффициенты; и (Р, г) - искомая скалярная функция поля. Переменная г > О — время. Условия (2) - есть условия непрерывности искомой функции и ( Р, г) и плотности потока (аУи (Р, г) , п) на внутренних границах контакта сред с различными тензорами а, здесь п - вектор нормали к границе у. На внешних границах области О будем рассматривать граничные условия третьего рода, которые в частных (в зависимости от функций а(Р) и в( Р )) могут иметь вид граничных условий первого или второго рода. Начальные условия (5) взяты однородными, поскольку без ограничения общности можно считать задачу сформулированной для отклонений от

* автор, ответственный за переписку

некоторого начального состояния значений поля. Условие (4) описывает поведение решения на бесконечно удаленных границах в неограниченных

подобластях области О.

В дальнейшем будем считать, что поле возбуждается точечным источником переменной интенсивности I(г) , расположенным в точке А области О , т.е. f (г) = I(г) Р — А) . Все встречающиеся в задаче функции будем считать достаточно гладкими для использования формул интегральных представлений и интегральных уравнений, а также имеющими необходимый порядок затухания на бесконечности для обеспечения применимости интегрального преобразования Лапласа.

Применим к задаче (1)-(5) способ решения, описанный в работе [1], используя интегральное преобразование Лапласа [3]

тс

и ?(Р) = | и (Р, г)е ?йг , с формулой обраще-

ния

^ «О +їж

и ( р, і ) = — /и т(Р)етйа,

2т з

«О —їж

(6)

Ке « = « > О,

Получим однопараметрическое (по С) семейство краевых задач:

йгу{а, Уи? (Р))— киС (Р) =

= — I ?8(Р — А), Р еОг , г = ;

и?( Р) — и ? ( Р)| ^] = О, (а. УиС( Р), п)—

— (а]Уи?(Р),п)^ = О, г = 1Ж ]е ,г;

а. (Р)(аУи?(Р), п)— Д (Р)иС(Р)| регг =

= ^?(Р), |а (Р)| + |Д(Р)| * О, г = ,г2,...,г*, к < N; и?(Р) ^ О, Р ^тс, ,

т = т1,т2,...,тп, п < N

в подобластях Ог- с симметричными эллипти-

(7)

(8)

(9)

(10)

ческими операторами

Н[и?(Р)] = йу(аУи?(Р))— ки?(Р), кг = а, + ? .

Интегральное представление решения задачи (7)-(10) имеет вид [1]: в- и ?(Р) =

= I I \и?(в)((а] —аг)Уа(Р,б) пе +

1= ^і +1 уе /ї уї пу.

I«

+1«• о(р, а)+ х (^гіу°(р,е), ие )^т 2Й + (11)

и гп Рл(а)

+а ^ 2й ■

12 Гпа12(а)

„I 1, Р * Г

Здесь в = ■< ; /1 - номера уча-

[1/2, Р е Г

стков внешней границы Г, на которых заданы условия первого рода, 12 - номера участков внешней границы с условиями второго или третьего рода; о(р, б) - функция Грина - функция точечного источника поля единичной интенсивности во вмещающем пространстве, состоящем из Nl подобластей (Nl < N).

Граничные значения функции и?(б) находятся как решение системы интегральных уравнений Фредгольма второго рода: и?(Р) — I I |и?(й)([а1 — а)Уе(Р,е),пв'йГб = I?-0(Р,б)+

+11£,‘П ^ (12)

+1е),пе +1\'аС§^(Р,б)йГ'2б,

п гІА,і(Є)'

где

Р е у;, / е ]к, к = N1 +1, N, б е у3, ] е ,, к = N1 +1, N . В (11) и (12) пд - вектор внешней нормали в точке

б, направленный на внутренних границах уг в

среду с меньшим, чем г номером.

Обращение преобразования Лапласа (6) программно реализуется с помощью обобщенных квадратурных формул наивысшей степени точности [4] вида:

С+їж

— / еРР ( р)Ф и

2т •1.

с—їж

и X А^(Рк ^ ^ >О,

(13)

где узлы рк и коэффициенты Ак выбираются из условий точности формулы (6) для набора функций f (р) = р ат, т = О,2NL — 1. Это равносильно выполнению равенств

ЕАр—“ =

1

к=1

, т = О, 2Ыг — 1,

где узлы являются корнями многочлена

Ы,

к=1

из

условий

определяемые однозначно

ортогональности:

| ерр~'Яп (р ~ *) р~атйр = О, т = О, NL — 1.

С—.тс

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

Вычислительный эксперимент В качестве апробации предлагаемого вычислительного алгоритма рассмотрена задача расчета

О

к=1

С + їж

температурного поля и(Р, і) ограниченного однородного цилиндра высоты 2Н и радиуса к при отсутствии внутренних источников тепла, Математическая модель задачи имеет вид:

1 ди(Р,і) 1 д г дттґП .ЛЛ з2г

a

dt

1 _д_І dU(P,t)J d2U(P,

r dr І dr J dz2 (14)

dU (P, t)

dz dU (P, t)

z=h

dr

Uf r=R

dU (P, t)

dz

= -— (U(P, t)|r=R - Ur );

dU (P, t)

z=0

dr

r=0

U (P,0) = U0 = const. Я

(13)

(1б)

(17)

Здесь а =------- коэффициент температуро-

ср

проводности, Л, с, р- теплопроводность, удельная теплоемкость и плотность среды, заполняющей цилиндр, соответственно. Начальное температурное распределение будем считать заданным равномерно с температурой иО . В начальный момент

времени поверхности оснований цилиндра начинают обмениваться теплом по закону Ньютона со

средой постоянной температуры и к Ф иО ; также

по закону Ньютона происходит теплообмен боковой поверхности со средой температуры

ив Ф иь Ф и о; аl, а2 — соответствующие коэффициенты теплообмена.

Аналитическое решение задачи (14)-(16) известно [5]:

U(P, t) = Ur +(Uh - Ur )x

„ a.2) J о I UnR Jch UnR x У---------T 1 1

Unk

+

Bi

sh ц k + ch ц k

h

+УУ

(2) (1) n Am

(Uh - Ur )иП

h R fr^n

f

■(Uh -т)

2 Я2 ^

U + — Un k2

і k

(1S)

x

U + Я Un k2

і k у

x J0І Un

Л z

x cos Яп — x exp mh

r R

f //2 Я2 ^

Un + Яи

vR2 h2 У

x

at

где /Іп — корни характеристического уравне-

^ О (Мп ) Ип

J1 Ц„ ) Bil

-, Я

корни характеристиче-

і Ї 1 • -lh

ского уравнения ctg Яи =---Яи , Bih = ——,

BL Я

| 2 Ji (цп) ,

Un [jО2 (Un )+ J12 (Un )j

и, = М, к = А, А» ■

В Л Я п

А<1) _ ______2 в*П Лт_______

т л л

Лт + вШ Лт ^ Лт

Сравнительный анализ результатов аналитического (по формуле (18)) и численного (с помощью программного модуля, реализующего предлагаемый способ вычислений) решений задачи приведен в табл. 1. Вычисления проводились в точке Р(1О м,1О м) при следующих значениях параметров: ио = 293 К, иА = 4ОО К , и, = 1ООО К,

а = 1.25 х 1О—5 м2/ с, В = 50 м, к = 50 м,

Л = 46,5 Вт /(м - К)).

Таблица 1

Сравнение решений

t, с Аналитическое решение, К Численное решение, К Абсолютная ногрешность, К Относительная погрешность, %

1 293,ООО296 293.ОООО16 2.80^ 10-4 9.56-10-7

2 293.ООО342 293.ОООО32 З.1040-4 1.06^ 10-6

З 293,ООО388 293.ОООО48 З.4040-4 1.16^ 10-6

4 293.ООО435 293.ОООО64 З.7110-4 1. 27 ■ 1 О-6

5 293.ООО481 293.ОООО8О 4.01 ■ 10-4 1.37 ■ 10-6

6 293.ООО527 293.ОООО96 4.З1 ■ 10-4 1.47 ■ 10-6

7 293.ООО573 293.ООО112 4.61 ■ 10-4 1.57 ■ 10-6

8 293.ООО619 293.ООО128 4.91 ■ 10-4 1.68^ 10-6

n=1 m=1

-Ї

0

0

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

Вычислительным экспериментом (рис.1) были определены верхние пределы сумм в формулах аналитического решения: при значениях

п = т = 18 погрешность вычислений достаточно мала. Тонкой сплошной линией на графике показана линия экспоненциального тренда абсолютной погрешности.

Количество верных цифр Ку в мантиссе вычисляемого разработанным программным модулем результата численного обращения преобразования Лапласа по обобщенным формулам наивысшей степени точности (13), в зависимости от

числа узлов NL квадратурной формулы, приведено в табл. 2.

Таблица 2

Зависимость количество верных цифр в мантиссе результата обращения преобразования Лапласа от числа узлов интегрирования.

К

4

13

5

19

Заключение

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

ЛИТЕРАТУРА

1. Кризский В. Н. О способе вычисления физических полей в кусочно-анизотропных средах. Часть I. Стационарные поля // Вестник Башкирского университета. 2009. Т. 14. №3.

С.726-730.

2. Кризский В. Н. О способе вычисления физических полей в кусочно-анизотропных средах. Часть II. Нестационарные поля // Вестник Башкирского университета. 2ОО9. Т. 14. №4. С. 13О2-13О6.

3. Деч Г. Руководство по практическому применению преобразования Лапласа и X-преобразования. М.: Наука, 1971. 288 с.

4. Матвеева Т. А. Некоторые методы обращения преобразования Лапласа и их приложения: дис.... канд. физ.-мат. наук. С.-П., 2003. 117 с.

5. Козлов В. П. Двумерные осесимметричные нестационарные задачи теплопроводности. Мн.: Наука и техника, 1986. 392 с.

Поступила в редакцию 31.10.2012 г.