close

Вход

Забыли?

вход по аккаунту

?

О вычислении температурных и диффузионных полей в кусочно-постоянных анизотропных средах.

код для вставкиСкачать
Вестник Башкирского университета. 2013. Т. 18. №2
ISSN 1998-4812
313
УДК 519.63
О ВЫЧИСЛЕНИИ ТЕМПЕРАТУРНЫХ И ДИФФУЗИОННЫХ ПОЛЕЙ
В КУСОЧНО-ПОСТОЯННЫХ АНИЗОТРОПНЫХ СРЕДАХ
© В. Н. Кризский*, А. Р. Бикбаева
Башкирский государственный университет, Стерлитамакский филиал
Россия, Республика Башкортостан, 453103 г. Стерлитамак, пр. Ленина, 49.
Тел./факс: +7 (3473) 43 25 80.
E-mail: krizsky@rambler.ru
Рассматривается способ вычисления температурных и диффузионных полей в кусочнопостоянных анизотропных и однородных средах на основе методов интегральных представлений и граничных интегральных уравнений. Приведены результаты сравнительного вычислительного эксперимента, демонстрирующие эффективность предлагаемого подхода.
Ключевые слова: анизотропная среда, температурное поле, диффузионное поле, способ
вычисления.
Введение
Изучение закономерностей распространения в
пространстве нестационарных температурных и
диффузионных полей в областях сложной геометрической формы, с анизотропией физических
свойств сред их наполняющих, тесно связано с решением параболических краевых задач математической физики. Разработка алгоритмов решения
подобного типа задач и программ расчета данных
полей имеет практическое значение для областей
экспериментальной физики, техники и технологии,
требующих высокоточных расчетов для оптимизации процессов, выбора параметров исследуемых
областей, сред и источников поля.
В данной работе, в продолжение результатов
работ [1, 2], обсуждаются результаты вычислительных экспериментов и вычислительные аспекты алгоритма численного расчета нестационарных физических полей применительно к задачам диффузии и
теплопроводности.
Постановка задачи и способ решения
3
Пусть область исследования Ω ⊂ R может
div (σ i ⋅ ∇U i ( P, t )) −
∂U i ( P, t )
=
dt
= − f i ( P, t ), P ∈ Ω i , i = 1, N ;
− aiU i ( P, t ) − bi2
U i ( P, t ) − U j ( P, t )
Ω = U Ω i – объединеi =1
ния N подобластей Ωi , каждая из которых заполнена анизотропным по теплопроводности (диффузии) веществом с постоянным симметричным тензором анизотропии
σ i . Обозначим через
M
Γ = U Γ j – внешнюю границу области Ω , а чеj =1
рез
γi
и Γj , соответственно, внутренние и внеш-
ние границы областей Ω i .
Процессы распространения тепла или диффузионного переноса вещества в области Ω описываются скалярными полями, математические модели которых представляются начально-краевыми
задачами математической физики вида:
* автор, ответственный за переписку
=
= 0, (σ i ∇U i ( P, t ), n ) −
− (σ j ∇U j ( P, t ), n ) γ i ∩γ j = 0,
(2)
i = 1, N , j ∈ J i =
= {j j = 1, i − 1; γ i ∩ γ j ≠ 0};
α j ( P)(σ j ∇U j ( P, t ), n ) −
− β j ( P)U j ( P, t )
P∈Γ j
= ψ j ( P, t ),
(3)
α j ( P) + β j ( P) ≠ 0, i =
= i1 , i 2 ,..., i k , k ≤ N ;
U m ( P, t ) → 0, P → ∞, m =
= m1 , m 2 ,..., mn , n ≤ N ;
N
быть представлена в виде
γ i ∩γ j
(1)
U i ( P,0) = 0, i = 1, N .
(4)
(5)
f i ( P, t ) – функции интенсивности источников/стоков поля в подобластях Ω i ; ai , bi –
постоянные в Ω i числовые коэффициенты;
U ( P , t ) – искомая скалярная функция поля. Переменная t ≥ 0 − время. Условия (2) – есть условия
непрерывности искомой функции U ( P , t ) и плотЗдесь
(
)
ности потока σ ∇U ( P , t ), n на внутренних границах контакта сред с различными тензорами σ ,
здесь n – вектор нормали к границе γ . На внешних границах области Ω будем рассматривать граничные условия третьего рода, которые в частных
(в зависимости от функций α (P ) и β (P ) ) могут
иметь вид граничных условий первого или второго
рода. Начальные условия (5) взяты однородными,
поскольку без ограничения общности можно считать задачу сформулированной для отклонений от
МАТЕМАТИКА
314
некоторого начального состояния значений поля.
Условие (4) описывает поведение решения на бесконечно удаленных границах в неограниченных
подобластях области Ω .
В дальнейшем будем считать, что поле возбуждается точечным источником переменной интенсивности I (t ) , расположенным в точке A области Ω , т.е. f (t ) = I ( t ) δ ( P − A) . Все встречающиеся в задаче функции будем считать достаточно
гладкими для использования формул интегральных
представлений и интегральных уравнений, а также
имеющими необходимый порядок затухания на
бесконечности для обеспечения применимости интегрального преобразования Лапласа.
Применим к задаче (1)-(5) способ решения,
описанный в работе [1], используя интегральное
преобразование
Лапласа
[3]
∞
U ( P ) = ∫ U ( P , t )e
−ωt
dt , с формулой обраще-
ния
Получим однопараметрическое (по
ство краевых задач:
div (σ i ∇U iω ( P) ) − k iU iω ( P ) =
N
i
i = N 1 +1 j∈J i
i
∩
j
(
(
− σ j ∇U ωj ( P ), n
)
ω
P∈Γi
l2
,
≈ ∑ Ak f ( p k ), s > 0,
k =1
(9)
где узлы p k и коэффициенты Ak выбираются
из условий точности формулы (6) для набора функ− am
ций f ( p) = p , m = 0, 2 N L
сильно выполнению равенств
NL
∑A p
k
∑ ∑ ∫ U ω (Q)((σ
i
j
− σ i )∇G (P, Q ), nQ )dγ iQ +
j
ψ lω1 (Q )
(σ l1∇G(P, Q ), nQ )dΓl 2Q + (11)
β l1 (Q )
Γ
+ I ω ⋅ G (P, Q ) + ∑ ∫
l1
ψ lω2 (Q)
G (P, Q )dΓl 2Q .
α l 2 (Q )
Γ
− am
k
k =1
)
N
(13)
NL
(8)
в подобластях Ω i с симметричными эллиптическими
операторами
ω
ω
ω
H U i ( P ) ≡ div σ i ∇U i ( P ) − kiU i ( P ), ki = ai + ωbi2 .
Интегральное представление решения задачи
(7)-(10) имеет вид [1]:
β ⋅ U ω ( P) =
l1
c + i∞
1
e p p − s f ( p )dp ≈
2πi c −∫i∞
(7)
(10)
m = m1 , m2 ,..., mn , n ≤ N
(
l2
где
P ∈ γ l , l ∈ J k , k = N 1 + 1, N , Q ∈ γ j , j ∈ J i , k = N1 + 1, N .
=
= ψ i ( P ), α i ( P ) + β i ( P ) ≠ 0, i = i1 , i2 ,..., ik , k ≤ N ;
]
l 2Q
l 2 Γl 1
ω ) семей-
ω
U mω ( P) → 0, P → ∞,
iQ
l 2Q
Q
l1
)
α i ( P )(σ i ∇U i ( P ), n ) − β i ( P )U i ( P )
ω
Q
l1
= 0, i = 1, N , j ∈ J i ;
γ i ∩γ j
i
j
Обращение преобразования Лапласа (6) программно реализуется с помощью обобщенных
квадратурных формул наивысшей степени точности [4] вида:
U iω ( P ) − U ωj ( P ) γ i ∩γ j = 0, σ i ∇U iω ( P ), n −
l1
ω
Граничные значения функции U i (Q ) находятся как решение системы интегральных уравнений Фредгольма второго рода:
U ω ( P ) − ∑ ∑ ∫ U ω (Q )((σ − σ )∇G (P , Q ), n )dγ = I ω ⋅ G (P, Q ) +
γ γ
(12)
ω
ψ ω (Q )
(σ ∇G(P, Q ), n )dΓ + ∑ ∫ ψα ((QQ))G (P, Q )dΓ ,
+∑ ∫
β (Q )
Q , направленный на внутренних границах γ i в
среду с меньшим, чем i номером.
= − I ω δ ( P − A), P ∈ Ω i , i = 1, N ;
l2
≤ N ).
добластей ( N1
(6)
Re ω = ω0 > 0.
+∑∫
)
В (11) и (12) nQ – вектор внешней нормали в точке
ω + i∞
1 0 ω
U ( P, t ) =
U ( P ) e ωt d ω ,
2πi ω0∫−i∞
i = N 1 +1 j ∈ J i γ i ∩ γ
(
l 1 Γl 1
0
=
ней границы с условиями второго или третьего рода; G P, Q – функция Грина – функция точечного источника поля единичной интенсивности во
вмещающем пространстве, состоящем из N1 по-
l1
ω
[
P ∉Г
 1,
; l1 – номера уча1 / 2, P ∈ Г
стков внешней границы Г , на которых заданы
условия первого рода, l 2 – номера участков внеш-
β =
Здесь
где
узлы
=
− 1 . Это равно-
1
, m = 0, 2 N L − 1 ,
Γ( s + am )
являются
корнями
многочлена
NL
Rn ( x ) = ∏ ( x − pk− a ), определяемые однозначно
k =1
из
условий
ортогональности:
c + i∞
∫e
p
p − s Rn ( p − s ) p − amdp = 0, m = 0, N L − 1 .
c − i∞
Способ решения распространяется и на кусочно-постоянные однородные среды, тензор анизотропии которых имеет вид диагональной матрицы с
одинаковыми значениями на диагонали.
Вычислительный эксперимент
В качестве апробации предлагаемого вычислительного алгоритма рассмотрена задача расчета
Вестник Башкирского университета. 2013. Т. 18. №2
ISSN 1998-4812
температурного поля U ( P , t ) ограниченного од-
U ( P, t ) = U R + (U h − U R ) ×
нородного цилиндра высоты 2h и радиуса R при
r
z

An( 2 ) J 0  µ n  ch µ n
R
 R
×∑
+
µ
k
n
n =1
sh µ n k + ch µ n k
Bih
отсутствии внутренних источников тепла. Математическая модель задачи имеет вид:
1 ∂U ( P, t ) 1 ∂  ∂U ( P, t )  ∂ 2U ( P,
=
r
+
a
∂t
r ∂r 
∂r 
∂z 2 (14)
∞

(U h − U R )µ n2 −


+ ∑∑  An( 2 ) Am(1) 
 2 λ2m
 µ n + 2
(
)
U
−
−
τ
n =1 m =1 
h

k



;
∞
∂U ( P, t )
α
= − 1 U ( P, t ) z = h − U h ,
λ
∂z
z =h
(15)
∂U ( P, t )
α2
= − (U ( P, t ) r = R − U R ) ;
∂r
λ
r=R
∂U ( P, t )
∂U ( P, t )
=0 ,
=0 ;
(16)
∂z
∂r
z =0
r =0
U ( P,0) = U 0 = const.
(17)
λ
– коэффициент температуроЗдесь a =
сρ
проводности, λ , с, ρ – теплопроводность, удель-
(
)
ная теплоемкость и плотность среды, заполняющей
цилиндр, соответственно. Начальное температурное распределение будем считать заданным равномерно с температурой U 0 . В начальный момент
времени поверхности оснований цилиндра начинают обмениваться теплом по закону Ньютона со
средой постоянной температуры U h ≠ U 0 ; также
по закону Ньютона происходит теплообмен боковой
поверхности
со
средой
температуры
U R ≠ U h ≠ U 0 ; α1 , α 2 − соответствующие коэффициенты теплообмена.
Аналитическое решение задачи (14)–(16) известно [5]:
315
∞

λ2
×  µ n2 + m2
k



  ×
 
 
(18)
−1

r
 × J 0  µ n  ×
 R

  µ 2 λ2  
z
× cos λm × exp −  n2 + m2  at ,
h
h  
  R
µ n − корни характеристического уравнеJ (µ ) µ
ния 0 n = n , λm − корни характеристичеJ 1 (µ n ) Bi R
αh
1
ского уравнения ctg λ m =
λ m , Bi h = 1 ,
λ
Bi h
h
α R
2 J 1 (µ n )
,
k = , An( 2 ) =
Bi R = 2 ,
R
λ
µ n [J 02 (µ n ) + J 12 (µ n )]
где
Am(1) =
2 sin λ m
.
λ m + sin λ m cos λ m
Сравнительный анализ результатов аналитического (по формуле (18)) и численного (с помощью программного модуля, реализующего предлагаемый способ вычислений) решений задачи приведен в табл. 1. Вычисления проводились в точке
P (10 м,10 м ) при следующих значениях параметров: U 0 = 293 K , U h = 400 K , U R = 1000 K ,
a = 1.25 × 10−5 м 2 / с , R = 50 м ,
λ = 46,5 Вт /( м ⋅ К ) ).
h = 50 м ,
Таблица 1
Сравнение решений
t,с
1
2
3
4
5
6
7
8
Аналитическое решение, K
293.000296
293.000342
293.000388
293.000435
293.000481
293.000527
293.000573
293.000619
Численное решение, K
293.000016
293.000032
293.000048
293.000064
293.000080
293.000096
293.000112
293.000128
Абсолютная
погрешность, K
2.80·10–4
3.10·10–4
3.40·10–4
3.71·10–4
4.01·10–4
4.31·10–4
4.61·10–4
4.91·10–4
Относительная
погрешность, %
9.56·10–7
1.06·10–6
1.16·10–6
1.27·10–6
1.37·10–6
1.47·10–6
1.57·10–6
1.68·10–6
МАТЕМАТИКА
316
Рис. 1. График зависимости абсолютной погрешности вычислений от значений верхних пределов суммирования.
Вычислительным экспериментом (рис.1) были
определены верхние пределы сумм в формулах
аналитического
решения:
при
значениях
n = m = 18 погрешность вычислений достаточно
мала. Тонкой сплошной линией на графике показана линия экспоненциального тренда абсолютной
погрешности.
Количество верных цифр Kν в мантиссе вычисляемого разработанным программным модулем результата численного обращения преобразования Лапласа по обобщенным формулам наивысшей степени точности (13), в зависимости от
числа узлов N L квадратурной формулы, приведено в табл. 2.
Таблица 2
Зависимость количество верных цифр в мантиссе результата обращения преобразования Лапласа от числа узлов
интегрирования.
NL
Kv
3
8
4
13
Заключение
Предложен способ вычисления нестационарных температурных и диффузионных полей в кусочно-постоянных анизотропных и однородных
средах. Разработано программное средство, реализующее построенный алгоритм решения. Апробация способа решения проведена при расчете температуры нагрева однородного ограниченного кругового цилиндра.
ЛИТЕРАТУРА
1.
2.
3.
4.
5
19
5.
Кризский В. Н. О способе вычисления физических полей в
кусочно-анизотропных средах. Часть I. Стационарные поля // Вестник Башкирского университета. 2009. Т. 14. №3.
С. 726–730.
Кризский В. Н. О способе вычисления физических полей в
кусочно-анизотропных средах. Часть II. Нестационарные
поля // Вестник Башкирского университета. 2009. Т. 14.
№4. С. 1302–1306.
Деч Г. Руководство по практическому применению преобразования Лапласа и Z-преобразования. М.: Наука, 1971.
288 с.
Матвеева Т. А. Некоторые методы обращения преобразования Лапласа и их приложения: дис…. канд. физ.-мат.
наук. С.-П., 2003. 117 с.
Козлов В. П. Двумерные осесимметричные нестационарные задачи теплопроводности. Мн.: Наука и техника,
1986. 392 с.
Поступила в редакцию 31.10.2012 г.
Документ
Категория
Без категории
Просмотров
4
Размер файла
1 746 Кб
Теги
поле, кусочно, вычисления, постоянный, диффузионные, среда, анизотропные, температурных
1/--страниц
Пожаловаться на содержимое документа