КОСМИЧЕСКИЕ ИССЛЕДОВАНИЯ, 2007, том 45, № 1, с. 93-96
КРАТКИЕ СООБЩЕНИЯ
УДК 502.3:517.956 225
ИСТОЧНИК ДИФФУЗИИ В ЭКСПОНЕНЦИАЛЬНОЙ АТМОСФЕРЕ
© 2007 г. Р. М. Зайдель
г. Москва Поступила в редакцию 20.10.2005 г.
PACS: 96.15.Hy
Вопросы тепло- и массообмена в разреженных и неоднородных средах, таких, как верхние атмосферы планет или звезд, давно привлекали внимание исследователей. В классической работе [1] рассмотрены нестационарные процессы теплообмена в той части высотных слоев земной атмосферы, где "турбулентность и конвекция, а часто и радиация играют малую роль в переносе тепла, и главное значение приобретает молекулярная теплопроводность".
Рассмотренные в [1] нестационарные одномерные задачи с экспоненциальным и степенным распределением плотности могут иметь значение для верхней атмосферы, если горизонтальные размеры области атмосферы, прогреваемой природным источником, достаточно велики по сравнению с высотой над землей, так что можно пренебречь оттоком тепла вбок. Теперь же, в связи с освоением ближнего космоса, появилась необходимость изучить происходящие в верхних слоях атмосферы процессы диффузии опасных для экологии Земли газов от источников техногенного происхождения. В данной статье решена задача о диффузии вещества, выделяемого стационарным точечным источником.
Пусть в начале координат находится постоянный точечный источник с мощностью Q0 частиц/с, диффундирующих в среде, в которой коэффициент диффузии зависит от высоты z по закону
D(z) = D0exp(az) = D0exp(arcos0); a> 0, (1)
где r, 0 - сферические координаты. Принятый выбор знака a означает, что основной поток частиц будет направлен вверх, где коэффициент диффузии больше. Эта ситуация характерна для атмосферы планет.
Обозначим f(r) - плотность частиц, плотность потока частиц
j = - D (r )-Vf( r) t
(2)
= 0 или, с учетом (2), уравнению &у[^(г) • У(г)] = 0. В силу (1) это уравнение приводится к виду
Af (r) + a f = 0.
Удобно сделать замену
a
f( r) = F( r) exp |--zj.
Для функции F(r) получим уравнение
2
A F( r) -a-F (r) = 0.
(3)
(4)
(5)
Здесь уже нет выделенного направления, поэтому решение для (5) будет сферически симметрич-
d2 F (r)
ным, удовлетворяющими уравнению -2"- +
dr2
2dF(r) a2 , n _ +---- F(r) = 0. Решением, исчезающим
при r —»- го, будет
F(r) = ~!exp^-ar); A = const. (6)
Подставляя (6) в (4), получим
f(r) = f(r, 0) = -exp
af _L. \
—(r + z) r
A
— exp
r
(7)
-a r (1 + cos 0)
Запишем радиальную компоненту плотности потока
m3f( r, 0) Do A jr = -D0exp(arcos0) J ' = —2" x
(8)
x
1 + 2 r (1 + cos 0)
exP I -ar(1 + cos0) f.
В стационарных условиях распределение частиц подчиняется уравнению непрерывности <1гу| =
При г —»- 0 ]г = Б А/Г2; поток через сферу малого радиуса равен 4%D0A = Я0; теперь вместо (7) и (8) получим
f (r) = f (r, 0) = во
—в— exp\ r( 1 + cos0) 4nD0r F| 2 v '
Jr
4 nr
1 + a r (1 + cos 0)
X
где R = Vp2 + z2. Проведя вычисления, получим: Jz(z > 0) = Qo; Jz(z = 0) = Qo/2; Jz(z < 0) = 0.
Физический смысл такого результата вполне очевиден: все те частицы, которые были испущены источником в нижнюю полусферу, вследствие уменьшения с глубиной коэффициента диффузии и при отсутствии поглощения, после многократного рассеяния в конечном счете отражаются наверх и обязательно пересекут плоскость z = 0.
Рассмотрим линии тока, форма которых в сферических координатах r, 0 определяется уравне-
dr rd0 тт ...
нием -.-. = ., -. -- „ ч . . Используя (9), получим дf /дг (áf/д Q )/z
форму линий тока в виде
(9)
X exp|-2r( 1 - cos0) j.
Вычислим поток частиц в верхнюю и нижнюю полусферу:
п/2
Jt = J 2пr2Jr(r, 0)sin0d0 =
= во
1-1exp í -a r
Ji =
J 2пr2Jr(r, 0)sin0d0 =
п/2
= | Qoexp r).
Полный поток через сферу произвольного радиуса r, естественно, равен J = J + J4 = const = Q0. Найдем также поток через горизонтальную плоскость на высоте z, записав функцию f(r) из (9) в цилиндрических координатах р, z, полагая r =
= TP
/ 2 , 2 /р + г :
f(p,г) =
во
4 п D0J р2 + z Отсюда находим
exP j-аа (^Р2 + z2 + z)
аг О f
Jz = ^ Tz =
во
4 п(р 2 + z2)
X
Ц
22 р + z
, а i 2. 2. , + 2 Up + z + z)
X
X
X expj-аа(^
Искомая величина равна
■z - z) j.
Jz = J2npdp • jz(p, z) = yexpí az
X
dr
r (0) = 2/a
2/a
sin0
1
1 - cos 0 sin0
ln
1 + cos0
1 - cos 0
1
ln
cos 0о
1 + c o s 0 1 + cos 0о
(10)
< о,
где 0О - направление движения частицы при выходе из источника. Условие г > 0 будет выполнено, если 0 < 0 < 0О; отметим также, что г —»- ^ при 0 —► 0. Из (10) следует, что угол 0 монотонно уменьшается от начального значения 00 до конечного значения г = 0; радиус г также возрастает монотонно от значения г = 0 до г =
Пусть источник мощностью Q1 находится в точке х = у = 0, г = Иъ в которой коэффициент диффузии имеет значение = ^0ехр(аИх). Вместо (9) теперь решением уравнения (3) будет функция
f 1 =
r=
Q1 \ a. , . exp j -2(r1+z1)
(11)
x2+ y2 ■
z1 = z - H j.
Аналогичным образом, плотность частиц от источника мощностью Q2, находящегося в точке х = = у = 0, г = И2 дается формулой
f2 =
в2
a
4 п D2 r2
exp j --(r2 + z2)
r2 =
x2+ y2-
(12)
?2 = z - H2,
D2 = Doexp (a H2).
X
ГР dP
J
-R + a (R + z)
exp(-aR
Рассмотрим случай, когда имеет место соотношение
H1--H2 = H.
(13)
о
п
о
о
ИСТОЧНИК ДИФФУЗИИ в экспоненциальной атмосфере
95
Сумма решений (11) и (12) на плоскости г = 0 будет равна
/1 + /2 =
1 Га,
ехР I -Тг0 х
4 п Я0 г
х &ехр(-(Н) + ^2ехрГ+ аН
г0 = Vх2 + у2 + Н2.
Отсюда видно, что если мощности 61 и 62 подчи нить условию
&ехр | - (н) + ^2ехр Г+(Н) = 0,
то на плоскости г = 0 будет выполняться условие поглощающей границы: при г = 0, / = 0. Таким образом, классический метод изображений, с небольшими изменениями вида (14), оказался в данном случае вполне успешным.
Пусть снова выполнено условие (13). При этом вместо (11) и (12) получим
/1 =
61 Г ацЛ [ а, , ч ехр| -тН ехр! --(г1 + г)
4 п Я0 г1
/2 =
1
62
= *]х2+ / + (г - Н)2,
4 п Я0 г.
ехр 1+°°Н)ехр<! (Г2 + г)
г2 = Vх2 + у2 + (г + Н )2. Составим производные: д/1 _ 61
х
х
Эг 4п Я0
г - Н 1 а Г г - Н
3 Г г1 г1 2 ( Г1
д/2 62
Эг 4п Я 0е
г + Н 1 а Г г + Н
а Н ~2
а
3 г2
г2 2
ехр! -22(Г1 + г)
а "2 Н
ехр! -((Г2 + г)
Полагая здесь г = 0, находим: ] 1г (г = о) = -Яо^.
Эг
г = 0
61 Г а„
= 4ПехрГ~2
-Н?+1. а г-НГ+1
|_ Г3 г 2
а
ехр IГ1,
(15)
72г( г = 0) = -Я0д/2 о г
62 Г, ^П = 44Пехр(+2Н х
х
НГ + 1. (ог
3 г 21 г
1-Г
а
ехр I-2 Г1.
Пусть теперь вместо (14) выполняется условие
61 ехр | -(н) = 62 ехр Г +(Н| = 6. (16)
Отметим, что при а = 0 мы получили бы сразу решение задачи с отражающей границей. Однако в нашем (14) случае, когда а Ф 0, ситуация более сложная. Сложим выражения (15), приняв во внимание условие (16):
Л(г =0) = ] 12(г =0) + ]22(г = 0) = 6 х
х
а (-(Ог), г = Vр2 + Н2; р = л/
(17)
х2+ у2.
Для полного решения задачи с отражающей границей необходимо найти функцию Ф(р, г), удовлетворяющую во всем пространстве уравнению (3), а на плоскости г = 0 дает поток, обратный по знаку величине из (17), при граничном условии Ф(р, г) —► 0 при г —► го. Из формулы (11) следует, что функция
у(р, г) =
а
4Пр2 + (г - Ну
х
(18)
х ехр! -а[7р2 + (г - Н)2 + г - Н]
удовлетворяет уравнению (3), которое в координатах р, г имеет вид
^ + + ^ + а^ = 0
Эр2 РдР Эг2 дг
(19)
Структура уравнения (19) такова, что производные ду/дг, д2у/дг2 и т.д. также удовлетворяют уравнению (19). Рассмотрим с этой точки зрения функцию
Ф(р, г) = |у(р,0dZ.
(20)
Нижний предел интеграла будет определен позднее. Из (20) следует
, . ЭФ ду ЭФ
г у -тг =
г = 0
2
г
0
Запишем уравнение (19) в переменных р,
д2 ¥(Р , С) + 1 Эу(р , С) + д2 ^(р, С) + Эр2 р дР дс2
Эу(р, С) =
ЭС
Предполагаем, что можно менять порядок дифференцирования по р и интегрирования по аргументу Тогда, с учетом (21), имеем:
+ а
= 0.
dp2
= д2 Ф( p , z ) dp2
1 dy(p, Z )• p dP .
dZ =
(22)
1 дФ( p, z) p д p .
Далее получим:
z 2
(•э у(р, z )dz = э^(р,о
dZ2
dZ
= d^(p , z) _ d^(p , Z ) dz dZ
= д2ф( p, z ) _ dy(p , Z )
(23)
Z = ;
dz1
dZ
1
^d^dZ = V(p,Z)|Z0 = V(p, z) - V(p, zo) =
(24)
= дФ( p, z) d z
- V(p, zo).
Собирая выражения (22)-(24), найдем, что Ф(р, z) удовлетворяет уравнению
э^Ф + 1ЭФ + э^Ф + аЭФ = -э¥(р,о +
dp2 pdp dz2
d z
dZ
Z = ;
(25)
Фф, z) = ^
2п D,
dZ
Vp2 + (Z - HУ
■ x
(26)
x exp
(Z - H) 2 + Z - H ]
Преобразуем (26), сделав замену Z - H = psh. Вместо (26) получим
ФФ, z) = ^
2 п D0
exp
а t -2 pe
(27)
to = ln
а
л
p2 + (z - H)2-
z - H
Полагая £ = ^ р^, выразим (27) через интегральную показательную функцию [2]:
Ф<р-2) = -4^ I! ^ = - ^ *(х)'
=а и
p2 + (z - H )2
4 п D0
z - H] > 0.
+ а^(рг0).
Из формулы (18) для у(р, £) следует, что при г0 —- | правая часть уравнения (25) обращается в нуль. Таким образом, для искомого решения получаем формулу
При г —- -I будет 0 < х <1; при этом имеет место разложение £1(х) = -у - 1пх + О(х), где у = 0.577 - постоянная Эйлера. Возрастание плотности частиц на большой глубине можно объяснить запутыванием частиц вследствие уменьшения длины свободного пробега при уменьшении коэффициента диффузии.
Рассмотрим пример. Коэффициент диффузии определяется той же формулой, что и коэффициент температуропроводности. Воспользуемся данными, приведенными в [1]; на высотах ~200 км коэффициент диффузии Б0 ~ 3 • 109 см2/с, а масштаб высот И ~ 50 км. Характерное время выхода на стационарный режим по порядку величины равно т = И2/Б0 = 8 • 103 с = 2 часа. Если за такое время мощность источника мало меняется, то изложенное выше решение может быть использовано.
Как отмечается в [1], на высотах ~120 км роль турбулентных пульсаций, которые способствуют перемешиванию вещества, уменьшается. Поэтому горизонтальную плоскость на этой высоте можно рассматривать как поглощающую границу и, следовательно, применять к данному случаю формулы
(11)—(14).
СПИСОК ЛИТЕРАТУРЫ
1. Голицын Г.С., Романова H.H. О распространении тепла в разреженных и неоднородных атмосферах // Геомагнетизм и аэрономия. 1970. Т. X. № 1.
2. Справочник по специальным функциям / Под ред.
М. Абрамовица и И. Стиган. М.: На
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.