CC BY
7
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ФИЗИКО-ХИМИЧЕСКИХ ПРОЦЕССОВ
УДК 533.2, 519.64
ОБ АНАЛИЗЕ ВЛИЯНИЯ СИЛ НА ТРАЕКТОРИЮ ВЗАИМОДЕЙСТВУЮЩИХ МОЛЕКУЛ И КОМПЬЮТЕРНЫЕ МЕТОДЫ ВЫЧИСЛЕНИЯ СРЕДНЕГО СЕЧЕНИЯ УГЛА РАССЕИВАНИЯ
АНИСИМОВА И.В., ИГНАТЬЕВ ВН.
Казанский государственный технический университет им. А.Н.Туполева, 420111, г. Казань, ул. К. Маркса, 10
АННОТАЦИЯ. В работе предлагаются компьютерные технологии вычисления среднего сечения угла рассеивания, значения которого входят в коэффициенты переноса. Компьютерные технологии включают в себя численный алгоритм решения нелинейной системы алгебраических уравнений, алгоритм вычисления несобственного интеграла, описывающего угол рассеивания взаимодействующих молекул и квадратуры с быстро осциллирующей подынтегральной функцией, соответствующей среднему сечению угла рассеивания.
КЛЮЧЕВЫЕ СЛОВА: коэффициенты переноса, взаимодействующие молекулы, угол рассеивания, компьютерные технологии, численный алгоритм, нелинейное уравнение, квадратуры с быстро осциллирующей подынтегральной функцией.
Одной из важных проблем механики сплошной среды является получение теоретически обоснованных формул для коэффициентов переноса. В настоящее время для решения этой проблемы широко используются кинетические уравнения Больцмана [1 - 4]. В зависимости от точности приближения решения получается то или иное выражение для коэффициентов переноса. В эти выражения входят различные виды квадратур и, как правило, со своими особенностями. Поэтому применение известных вычислительных технологий для вычисления той или иной квадратуры из кинетической теории газов требуют тщательного математического обоснования их эффективности применения. Так, для вычисления значения функции угла рассеивания, которая выражается несобственным интегралом первого рода, значение нижнего предела в нём определяется из решения нелинейной системы алгебраических уравнений, алгоритм решения которой предложен авторами [5, 6].
Квадратура, которая описывает среднее сечение угла рассеивания, при некоторых значениях
*
параметра Ь , имеет сильно осциллирующую подынтегральную функцию. Это свойство
подынтегральной функции требует применения обоснованных вычислительных технологий, учитывающих осцилляционный характер подынтегральной функции. Следует отметить важность асимптотического анализа влияния сил на характер поведения функции угла рассеивания. Так, на основе асимптотического анализа обобщенного потенциала
взаимодействия, предложенного автором [5], приведены соотношения, связывающие
*
критическое значение Ь0 с безразмерным числом Л1, характеризующим отношение длины расстояния, на котором функция потенциала меняет знак к расстоянию между центрами взаимодействующих молекул.
Среднее сечение угла рассеивания при столкновениях молекул газа математически описывается с помощью следующей О - квадратуры [1, 2]:
(1)
где
X
(Ь*,g*)=п-2Ь* |
^Т У(г*)
(г* 1„л1 г */ ( g * )
(2)
описывает угол рассеивания взаимодействующих молекул газа. В (1) - (2) использованы следующие обозначения:
Ь* = — - безразмерный прицельный параметр;
У
* г
г =--приведенное межмолекулярное расстояние;
У
(g •)
2 т • g
2
2а
- безразмерная относительная скорость;
^ (г * )= 4
А/
*
г
2(и+3)
- 4
6
+
V г /
^ * • ь * )2
(г* )2 ;
(3)
а (г0 )тт в (2) является минимальным положительным корнем нелинейного уравнения [1, 2]:
\ 2 г
/(г*)=(г;)-() ^
А/12'"*31 ( А/
V г0 J
V го J
+
с g * • ь* 12
V 'о J
= о
(4)
Математическая модель (1) - (3) позволяет оценить величину вклада определенного вида столкновений молекул на переносные характеристики газовой смеси. Этот факт исследуется здесь с помощью потенциала взаимодействия молекул (3), который является обобщением известных и широко используемых потенциалов взаимодействия, таких как Леннарда-Джонса и др. [1 - 4].
Далее проведем асимптотический анализ поведения функции х(Ь*, g *) в зависимости
* *
от параметров Ь , g и подынтегральной функции в несобственном интеграле (2):
1. Если в (3) величина
/ * * \2 g •Ь
очень мала по сравнению с силами притяжения и
отталкивания, то значение величины (г0 ) в этом случае находится в области действия силы отталкивания. Подынтегральное выражение в квадратуре (2) обратно пропорционально
величине
*2 * * / * \ | А/
2 (г ) = 4| —
ч 2(п+3)
Р
где Ре
^6
- 4 —— . Основной вклад в величину угла
.г J V г J
рассеивания (2), соответственно, и переносные характеристики газа здесь вносятся при
* *2 * *2 * / \ значениях г для которых разность g - ре# > 0 будет малой, т.е. g -р*^ ~0(к), где
к << 1. Отметим, что при выполнении неравенства g 2 - р~ > 0 всегда существует такая
точка г е [г0, да), в окрестности которой имеет место 0 < ^ - р^ << 1. В этой точке г ,
*
расположение которой определяется значением g , и ее окрестности значения подынтегральной функции в интеграле (2) принимают большие значения, что существенно
будет сказываться на значениях коэффициентов переноса. Таким образом, основной вклад в
~ *
величину угла рассеивания для этого случая вносят те значения ~ , для которых разность
*2 * г g - р^ мала, т.е. подобласти, где кинетическая энергия взаимодействующих молекул
сравнима с величиной их потенциальной энергии.
*2
со
2
*
*
4
г
2. Если величина
/ * * \2 g • ь
V г у
сравнима со значением потенциальной функции
ч2(и+3) чб
(Г * )= 4|
А/
е. V / *
- 4
А/
*
*2 *2 * * А) Для малых значений величины g условие g - (г ) > 0 будет удовлетворять
г у V г у
то возможны несколько типов соударения молекул:
*2
яе " - ре# У
**
при значениях г > г0, т. е. справа от точки максимума. Причём сила отталкивания при этих
*
же значениях г не играет существенной роли. Поэтому основной вклад на значение угла рассеяния будет вносить сила притяжения. При выполнении соотношения (г ) = g 2
в подобласти действия этой силы существуют точки поворота, в которых относительные
скорости молекул малы. Частицы могут перейти в процесс вращения вокруг общего центра
масс, т. е. в окрестностях этих точек может возникнуть орбитальное движение. Однако
принцип обратимости во времени запрещает устойчивость такого движения и в дальнейшем
происходит разлёт частиц. Таким образом, в этом случае подобласти, где кинетическая
энергия сравнима с потенциальной энергией, существуют справа от точки максимума и
определяются значением силы отталкивания.
*2 * / * \ Б) Если g будет порядка максимального значения величины (г ) , то величина
*2 ~ * — — * * относительного движения g - будет мала в подобласти г ^ г0 . Точки поворота
определяются в области действия силы отталкивания, где значение относительной скорости,
* * / * * \
будет мало. Также как и в случае А) при г < г1 > г0 ) возможно неустойчивое
орбитальное движение молекул вокруг центра масс. При таких значениях параметра g 2 основную роль в потенциале (3) играют силы отталкивания и, соответственно, влияние параметра п будет значительным в этом случае на коэффициенты переноса.
3. Когда величина
/ * * \2 g •ь
*
велика, то член, определяемый центробежными силами
V г у
**
в выражении (3) р. (г ), значительно превосходит потенциал притяжения и поэтому
процесс притяжения молекул для этого случая не играет почти никакой роли и не влияет на
значения коэффициентов переноса. Центробежные силы и силы отталкивания в этом случае
являются определяющими и наибольший вклад в интеграл (2) вносится при значениях * / * * \ * г е ^г0, г1 ). Проанализируем те значения г , при которых могут возникать орбитальные
*2
движения частиц. Как и для всех рассмотренных выше случаев, величина параметра g
* ~*
связана со значением параметра Ь и неравенством
т *2
*2 * Ь
g > (б)
*2
г
Ь *2
ь ** * / * * \
В неравенстве (б) величина при г ^ г0 и малых Ь ^Ь << г0 ) является
г
незначительной. Поэтому усилим (б) его неравенством:
g *2 ^ . (7)
Как следует из этого неравенства, при малых значениях параметра Ь подобласть
*
возможного орбитального движения существует где g 2 = р*., и при любых значениях
*2
Ь
этого параметра, где имеет место g > р ~ +
Прежде чем перейти к обсуждению численных алгоритмов определения значений
несобственных интегралов (1) и (2), рассмотрим алгоритм определения «критических»
**
параметров г0 и Ь0. Для их определения необходимо совместно решить следующую систему нелинейных уравнений:
(г;)-( ь )
g
4
С А/12("+3) С А/16
*
г
V г0 J
*
г
V г0 J
+
С г * • Ь* 12 60 Ь0
= 0;
4
/ \2(«+3)
А/
-4
г \б А/
г
V г0 у
*2
+ •
г
V г0 7
*2
= g 0
*2
(8)
(9)
Как показал асимптотический анализ [5], при значении (г0 ) следует, что
*2 *2 2
Ь0 = g 0 • А/ ,
**
т. е. величина Ь0 связана с параметрами g0 и А/ соотношением:
= А/ из уравнения (9)
Ь = g0 • А/. (10)
**
Из (10) следует, что критические значения Ь0 и g 0 пропорциональны с
коэффициентом пропорциональности А/. Эта зависимость геометрически изображена на рисунке.
Ь
А/ = 4
#0
Рис. Зависимость Ь0 от параметров А/ и g 0
Более точные значения «критических» величин г0 и Ь0 определялись в результате
0 Лшйа/
и
совместного решения системы (8)-(9) методом Ньютона. Начальные значения (г0 ).
(* \ *
Ь0) задавались на основе приближения (10). При заданном g0 критические значения
**
параметров Ь0 и г0 , полученные на основе численного решения системы (8)-(9) согласуются с асимптотическими оценками
(* \ * г ) = А/, Ь = g
0 /тт ' 0 о
А/.
(11)
*
2
г
*
2
г
0
--<
г
0
0
*
г
0
(* * \
Ь , g ),
являющейся аргументом в подынтегральной функции в О — интеграле.
(* * \ * *
Ь , g ) от параметров Ь и g впервые
доказана в работе [7].
*
Проведём анализ поведения этой функции в зависимости от параметра Ь при
**
фиксированных значениях величины g . Так, если значения Ь стремятся к критическому
*
Ь0, при котором взаимодействующие молекулы переходят на неустойчивую орбитальную
траекторию, то подынтегральная функция в квадратуре (1) становится быстро осциллирующей [2]. Компьютерные вычисления значений несобственного интеграла с быстро осциллирующими функциями связаны с определёнными сложностями [8, 9]. На основе анализа о поведении подынтегральной функции в квадратуре (1), исходный интервал интегрирования [о, да) разбивался на сумму трёх квадратур. Первая квадратура это определённый интеграл, длина отрезка интегрирования которого составляет участок монотонного поведения подынтегральной функции в квадратуре (1) (Ь е [о, Ь0 - е]). Вторая
квадратура есть участок осцилляционного поведения подынтегральной функции, т.е.
* * *
Ь е[Ь0 — е, Ь0 + е]. Данный участок интегрирования квадратуры (1) соответствует
неустойчивой «орбитальной» траектории взаимодействующих молекул. Далее происходит
«разлёт» молекул с этой траектории, и дальнейшее их движение описывается квадратурой с
**
монотонной подынтегральной функцией на интервале Ь е[Ь0 + е, да] (третья квадратура). На основе приведённого анализа подынтегральной функции, представим квадратуру (1) в виде суммы, где:
Q=J(1 - cos Xl/) b*db* = 1 J, (12)
0 k=1
где Jk (k = 1,2,3) определялись выражениями:
, * _ ,* b0 ^ b0 + ^ i \ ^ J = J (1 _ cosX) b*db*, J2 = J (1 _ cosx ) b*db*, J3 = J (l _ cos%v) b*db*. (13)
0 , * \ ' » *
bo _s bo +s
**
Значение величин Ь0 и г0 , используемое в квадратурах (2) и (13), определялось из
*
решения системы (8) - (9) численным методом [5 - 6]. После определения значений г0 из
* * *
уравнения (9) определялась величина Ь0 (~0 ) при фиксированном значении величины g0.
После этого осуществлялось вычисление квадратур (13). При этом использовались компьютерные технологии, основанные на ортогональных многочленах Гаусса [8, 9] и кубической сплайн-интерполяции [10 - 11]. Так как на рассматриваемых отрезках разбиения подынтегральные функции у квадратур и /3 является монотонными, то приближённое их
интегрирование основано на идеи минимизации погрешности используемого численного метода [9]. Вычисление их значений осуществлялось адаптивными численными методами [9]. У квадратуры J2 подынтегральная функция является быстро осциллирующей, поэтому
авторами проводился анализ эффективности применения следующих численных алгоритмов: а) алгоритм Гаусса - Лежандра, б) алгоритм Гаусса - Кронрода, в) алгоритм Гаусса -Чебышева, г) алгоритм кубических сплайнов.
Анализ эффективности численных методов для определения значения несобственного интеграла J2 с быстро осциллирующей подынтегральной функцией осуществлялся на примере квадратуры, имеющей подынтегральную функцию близкую по поведению к подынтегральной функции в J2:
. *
Ьо
J2 = | (1 - СОБ0Х) хйХ =
X2 1 , ч COSаX ---( X • SmаX )+--
2 а а
(14)
ьо
Поэтому эффективность того или иного рассматриваемого численного метода заключалась в сравнении результатов вычисления квадратуры (14) по методам а) - г) с точными значениями в зависимости от этого параметра. Из приведенной табл. 1 можно сделать заключение, что среди приближенных методов вычисления интегралов с быстро осциллирующей подынтегральной функцией следует отметить квадратуры Гаусса, основанные на применении ортогональных многочленов. В этих методах узлы разбиения, основанные на минимизации погрешности вычисления, строятся заранее. Адаптивные квадратуры Ньютона-Котеса используют алгоритм половинного деления в процессе вычисления интеграла из условия минимизации погрешности. Поэтому квадратуры Гаусса, в отличие от квадратур Ньютона-Котеса, являются более эффективными для вычисления интегралов с подынтегральными функциями, имеющими большие градиенты. Среди квадратур Гаусса для вычисления интеграла вида (1) можно выделить методы Гаусса-Кронрода и Гаусса-Лежандра, так как в них используются многочлены, ортогональные на полуинтервале [о, да).
Таблица 1
Сравнение результатов вычисления квадратуры J2 численными методами а) - г)
Ь*
о
о
Метод Решение Абсолютная ошибка
а = 1
Точный 499173,5580803917
Гаусса-Лежандра 499173,5580803916 0
Гаусса - Кронрода (адаптивный) 499173,5580803916 0
Гаусса - Чебышева 499173,5580803911 6,40284 •Ю-10
Кубический сплайн 499173,5580789156 1,47615 •Ю-10
а = 5
Точный 500197,6271010171
Гаусса - Лежандра 500197,6271010171 0
Гаусса - Кронрода (адаптивный) 500197,6271010171 0
Гаусса - Чебышева 500197,6271010171 0
Кубический сплайн 500197,6271010072 0 9,8953 •Ю-9
а = 10
Точный 500030,5809604425
Гаусса - Лежандра 500030,5809604425 0
Гаусса - Кронрода (адаптивный) 500030,5809604425 0
Гаусса - Чебышева 500030,5809593803 1,06212 • 10 -6
Кубический сплайн 500030,5809604497 7,21775 • 10 -9
а = 20
Точный 499970,9012289010
Гаусса - Лежандра 499970,9012289010 0
Гаусса - Кронрода (адаптивный) 499970,9012289010 0
Гаусса - Чебышева 499970,9012294314 5,30446 • 10 -7
Кубический сплайн 499970,9012293584 4,57396 • 10 -7
Далее, используя компьютерные технологии, состоящие из численных методов решения нелинейной системы уравнений (8) - (9), квадратуры (2) и приведённых выше вычислительных методов квадратуры (1), осуществлялось вычисление значений среднего угла рассеивания взаимодействующих молекул, т. е. значений О - интеграла. В табл. 2 приведены значения О - квадратуры в зависимости от параметров п и Л1, т. е. вида обобщенного потенциала взаимодействия (3). Анализ таблицы показал, что значения среднего сечения угла рассеивания зависят от данных параметров и это сказывается на значениях коэффициентах переноса.
Таблица 2
А Т *
Значения квадратуры (1) в зависимости от параметров п, Л1, g
А1 * g
1 2 3 4 5
п=2
0,5 100,257 100,258 100,259 100,259 100,259
1,0 100,187 100,238 100,251 100,255 100,256
1,5 100,324 100,296 100,283 100,276 100,271
2,0 100,335 100,342 100,337 100,330 100,324
2,5 100,370 100,366 100,362 100,358 100,355
3,0 100,377 100,389 100,385 100,380 100,377
3,5 100,414 100,404 100,403 100,399 100,396
4,0 100,442 100,426 100,419 100,415 100,412
4,5 100,464 100,445 100,437 100,432 100,427
5,0 100,480 100,462 100,453 100,447 100,442
п=3
0,5 100,257 100,258 100,259 100,259 100,259
1,0 100,199 100,232 100,249 100,254 100,256
1,5 100,333 100,304 100,289 100,281 100,276
2,0 100,344 100,346 100,343 100,338 100,333
2,5 100,368 100,370 100,365 100,363 100,361
3,0 100,385 100,392 100,390 100,386 100,383
3,5 100,418 100,406 100,406 100,404 100,402
4,0 100,450 100,430 100,424 100,420 100,417
4,5 100,472 100,450 100,443 100,437 100,433
5,0 100,487 100,467 100,459 100,453 100,449
п=4
0,5 100,257 100,258 100,259 100,259 100,259
1,0 100,213 100,227 100,248 100,253 100,255
1,5 100,340 100,310 100,294 100,285 100,279
2,0 100,352 100,349 100,348 100,343 100,339
2,5 100,327 100,373 100,368 100,366 100,364
3,0 100,394 100,394 100,393 100,390 100,387
3,5 100,425 100,408 100,408 100,407 100,406
4,0 100,457 100,434 100,427 100,423 100,421
4,5 100,478 100,454 100,446 100,441 100,438
5,0 100,492 100,471 100,463 100,458 100,454
СПИСОК ЛИТЕРАТУРЫ
1. Чепмен С., Каулинг Т. Математическая теория неоднородных газов. М. : ИЛ, 1960. 510 с.
2. Ферцигер Дж., Капер Г. Математическая теория процессов переноса в газах. М. : Мир, 1976. 553 с.
3. Зубарев Д.Н. Неравновесная статистическая термодинамика. М. : Наука, 1975. 440 с.
4. Резибуа П., Ленер М. Классическая кинетическая теория жидкостей и газов. М. : Мир, 1980. 424 с.
5. Анисимова И.В. О нелинейном уравнении в кинетической теории газов и его решениях // Вестник Казанского гос. техн. ун-та им. А.Н.Туполева. 2009. №3. С.74-77.
6. Анисимова И.В., Игнатьев В.Н. Об алгоритме локализации некратных решений нелинейного уравнения и их вычислений // Вестник Казанского гос. техн. ун-та им. А.Н.Туполева. 2010. №2. С.69-72.
7. Анисимова И.В., Гиниятуллина Р.Р. О непрерывности несобственного интеграла от параметра при вычислении угла рассеивания в кинетической теории газов // Вестник Казанского гос. техн. ун-та им. А.Н.Туполева. 2010. №3. С.93-100.
8. Калиткин Н.Н. Численные методы. М. : Наука, 1978. 512 с.
9. Бахвалов Н., Жидков Н., Кобельков Г. Численные методы. М. : Бином, 2003. 632 с.
10. Марчук Г.И. Методы вычислений. М. : Наука, 1980. 608 с.
11. Алберг Дж., Нильсон Э., Уолш Дж. Теория сплайнов и её приложения. М. : Мир, 1972. 319 с.
ABOUT THE ANALYSIS OF FORCES INFLUENCE ON THE TRAJECTORY OF COOPERATING MOLECULES AND COMPUTER METHODS OF THE EVALUATION OF AVERAGE CUT OF DISPERSIONS ANGLE
Anisimova I.V., Ignatiev V.N.
Kazan State Technical University, Kazan, Russia
SUMMARY. Computer process engineerings of evaluation average cut angle of dispersion which values are included into coefficients of transposition are offered in work. Computer process engineerings include numerical algorithm of solution nonlinear system of the algebraic equations, algorithm of evaluation improper integral describing a dispersion angle of cooperating molecules and a quadrature with fast oscillating integrand, corresponding average cut of an angle of dispersion.
KEYWORDS: coefficients of the transposition, cooperating molecules, angle of dispersion, the computer process engineerings, numerical algorithm, the nonlinear equation? quadratures with fast oscillating integrand.
Анисимова Ирина Викторовна, кандидат физико-математических наук, доцент кафедры ВМ КГТУ им. А.Н.Туполева, тел. (843)231-02-25, e-mail: aras/mova/V1@ramWer.rM
Игнатьев Виктор Николаевич, доктор физико-математических наук, профессор, заведующий кафедрой высшей математики КГТУ им. А.Н.Туполева
