численным алгоритм расчета сверхзвуковых течений, основанный на квазигазодинамических уравнениях Текст научной статьи по специальности «Физика»

11. Callan C.G., Dashen R., Gross D.J. // Phys. Rev. 1978. D17. P. 2717.

12. DeGrand T., Hasenfratz A., Koisács T. // Prog. Theor. Phys. Suppl. 1998. 131. P. 573.

13. Nambu Y. // Phys. Rev. 1974. DIO. P. 4262.

14. 't Hooft G. // Nucl. Phys. 1978. B138. P.l.

15. Nielsen H.B., Olesen P. // Nucl. Phys. 1979. B160. P. 380.

16. Del Debbio L., Faber M., Greensite J., Olejnik S. // Phys. Rev. 1997. D55. P. 2298.

17. de Forcrand P., D'Elia M. // Phys. Rev. Lett. 1999. 82. P. 4582.

18. Ambjorn /., Olesen P. // Nucl. Phys. 1980. B170 [FS1], P. 265.

19. Olesen P. // Nucl. Phys. 1982. B200 [FS4], P. 381.

20. Cornwall J.M. // hep-th/9901039.

21. Antonov D., Ebert D. // hep-th/9812112.

22. Diakonov D. // Mod. Phys. Lett. 1999. A14. P. 1725, 1909; hep-th/9905084, hep-th/9908069.

23. Жуковский Б.Ч. /У Вестн. Моск. ун-та. Физ. Астрон. 2000. №5. С. 8 (Moscow University Phys. Bull. 2000. N 5. P. 10).

24. Zhukovsky V.C. 11 Int. J. Mod. Phys. 2002. A17. P. 914.

25. Diakonov D., Maul M. 11 hep-lat/0204012.

26. Länge J.D., Engelhardt M., Reinhardt H.// hep-th/0301252.

27. Bordag M. // hep-th/0211080.

28. Banks Т., Casher A. // Nucl. Phys. 1980. B169. P. 103.

29. Ландау Л.Д., Лифшиц E.M. Квантовая механика (нерелятивистская теория). М., 1963.

30. Cooper F., Khare A., Musto R., Wipf A. 11 Ann. Phys. 1988. 187. P. 1.

31. Тернов И.М., Жуковский Б.Ч., Борисов A.B. Квантовые процессы в сильном внешнем поле. М., 1989.

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

УДК 517.958:533.7

ЧИСЛЕННЫЙ АЛГОРИТМ РАСЧЕТА СВЕРХЗВУКОВЫХ ТЕЧЕНИЙ, ОСНОВАННЫЙ НА КВАЗИГАЗОДИНАМИЧЕСКИХ

УРАВНЕНИЯХ

Т. Г. Елизарова*), М. Е. Соколова

(кафедра математики) E-mail: sokolova@afrodita.phys.msu.ru

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

Введение

В работах [1-4] описана система квазигазодинамических (КГД) уравнений, расширяющая возможности классической модели Навье-Стокса (НС) для описания течений вязкого сжимаемого газа. В области применимости уравнений НС дополнительная диссипация, присутствующая в КГД уравнениях, слабо влияет на решение и может использоваться для обеспечения устойчивости численных алгоритмов. Для разреженных течений в микроканалах дополнительные диееипативные слагаемые позволяют получить решение, которое лучше, чем модель НС, описывает данные эксперимента [5].

В статье [6] на основе КГД модели был предложен оригинальный алгоритм расчета нестационарных оеееимметричных течений вязкого газа. В отличие от разработанных ранее численных алгоритмов, опирающихся на КГД уравнения [1, 7], дан-

ный метод базируется на записи КГД уравнений в инвариантном виде, что позволяет естественным образом адаптировать его к различным, в том числе и криволинейным, системам координат. Разностные аппроксимации строятся в потоковой форме непосредственно для векторов плотности потока массы, теплового потока и тензора вязких напряжений, что соответствует записи уравнений газовой динамики в виде законов сохранения и делает алгоритм достаточно компактным и экономичным. НС слагаемые и дополнительные КГД диееипативные слагаемые выписаны отдельно, что позволяет независимо друг от друга варьировать значения коэффициентов вязкости, теплопроводности и релаксационного параметра. В настоящей работе указанный алгоритм реализован для плоской и цилиндрической геометрии и тестируется на примере задач вязкого течения в окрестности цилиндрического торца и невязкого течения в канале с уступом [7, 8].

*•' Институт математического моделирования РАН.

Численнш алгоритм

Приведем описание численного алгоритма для оеееимметричных течений согласно [6]. В этом случае система КГД уравнений имеет вид:

др + 1 д{г]тг) + д%

т

= о,

г дг дг

д{риТ) + 1 д{г]тгиТ) + дЦт7ит) др

дг

т

г дг дг

1 д(гЛгт) д(Пгг) _ Пфф г

дг дг г

д(риг) ^ 1 д{г]тгщ) ^ д(1тгиг) др

дг

Ы

г дг дг

= 1_а(гпгг) , д{ъгг)

г

(1) (2)

(3)

дг дг '

дЕ | 1 д(г]тгН) | дЦтгН) | 1 д(гдТ) | % дЬ г

дг

дг

дг

дг

1 д д

= г~дг^Г (ПггИг + + + Пггиг).

(4)

Здесь р — плотность газа, иг и иг — проекции вектора скорости и на оси Ог и Ог соответственно, р — давление, Т — температура, -у — показатель адиабаты, Я — газовая постоянная, Е и Н — полная энергия единицы объема и полная удельная энтальпия, которые вычисляются по формулам:

Е = р

р

и Н =

(Е+р)

I 7-1 р

Компоненты вектора плотности потока массы j вычисляются по формулам:

Згаг = р(иг - ™г), Зтг = Р(Щ ~

где

г

и)г = — р

Т

ПЗХ = -

Р

I д . 2ч д , \ др -^(грит) + ^(рщщ) +—

I д . д 2ч др

Компоненты тензора вязких напряжений П определяются как

1 -1-7*7* — Пгг | И/р И)у.

Н*

ттЯ5 о диг 2 . ПГГ = 2г1 - з'Л и>

П г г = + 11г IV х,

П^5, = П^5 = п ( ^ +

'\дг дг

Пгг = + IIц И}г,

П фф — П фф

К*

П™ = Дуй,

щ IV* +

где П^, П^/, П^, П^/ - компоненты

тензора вязких напряжений НС, а выражения для XV*, го*, Я* и дивергенция вектора скорости сНу и задаются формулами:

ги„ = т

риг-

риг-

дит дг

дщ дг

диг риг-^т- + -7Г

дг ди,

др дг

■ рих-

др дг

Н* = т

иг

др др

дг

и,

дг

1 д(гиг) шу и = - -

дг

диг дг

г дг

Компоненты вектора теплового потока q вычисляются как

д

трщ.

|_7.

1 дг

и.

д

7'

<ь = <£3-

трщ

щ

д

■ риг 7"—

дг

д

1дг д дг

р

1 дг \р д

д

~рщ

дг

1дг д

рщ-

где НС-слагаемые д^5 формулам:

дт

дг '

1?3 =

-к-

№

дг

определяются

по

-к

дт

дг '

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

П = Поо

Т

к =

7 Я

■V,

т =

■Г].

(7 — 1) Рг " ' рБс

Здесь »700 — известное значение г) при температуре Тоо, Рг и Эс - числа Прандтля и Шмидта соответственно, ш — заданный показатель степенной зависимости из промежутка [0.5; 1].

Система уравнений (1)-(4) дополняется начальными и краевыми условиями.

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

с первым порядком. Для численного решения начально-краевой задачи (1)-(4) используется явная по времени разностная схема.

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

Уравнения приводятся к безразмерному виду с помощью замены:

z^z/Rc, r^r/Rc, t У t Cqo/Rc, p^p/poo,

T ^ T JTqq , Ur У Ur j Cqq J Uz У Uz j CQQ J

£7->£7/(poo4,), // -н- llf<:i.

где Rc — линейный размер, c00 — скорость звука, Poo — плотность, Tqo температура в невозмущенном потоке. Обезразмеривание не изменяет вида уравнений.

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

ти и теплопроводности вычисляются как

- Ма I ( Р

Т ~ ReSc р V ~р

r] = rpSc, к =

h

а ■

тр Sc

Рг (7-1)'

где р = рТ/-у, Ма = («г)оо/соо и Ке= ((«г)оо х хроо-йс)/?7оо числа Маха и Рейнольдса соответственно, К = + К\, Кт, — шаги разностной сетки по пространству, а — численный коэффициент порядка единицы.

Течение в окрестности цилиндрического

торца

Рассматривается задача обтекания цилиндрического торца радиуса Rc однородным сверхзвуковым потоком вязкого сжимаемого газа под нулевым углом атаки.

В качестве граничных условий на входной (правой) границе выбраны условия набегающего потока Роо = 1, (иг)оо = — Ма, («г)оо =0, Роо = 1/7• На оси симметрии ставятся «условия симметрии», на торце цилиндра — «условия скольжения» для скорости, на боковой поверхности — «условия прилипания» для скорости. Поверхность цилиндра предполагается адиабатической. Для аккуратной аппроксимации уравнений вблизи угловой точки (2,1) об-

Таблица 1

Расчет осесимметричного течения

Ма а Сетка h у — h 2, ей Невязка Ps Ps Ts

Теоретические значения, рассчитанные по (5): 2.172 2.280 1.750

1.5 1.0 120 х 120 0.05 Ю-3 10"5 2.112 2.248 1.774

1.5 0.5 120 х 120 0.05 М-3 10"5 2.156 2.280 1.763

1.5 0.5 160 х 120 0.025 5 • М-4 10-4 2.186 2.305 1.757

1.5 0.2 120 х 120 0.05 М-3 10"5 2.189 2.303 1.753

Теоретические значения, рассчитанные по (5): 2.719 3.807 2.333

2 1.0 80 х 60 0.05 10"3 10"5 2.648 3.704 2.331

2 0.5 80 х 60 0.05 М-3 10"5 2.720 3.781 2.317

2 0.2 80 х 60 0.05 М-3 10"5 2.781 3.843 2.303

2 1.0 160 х 120 0.025 5 • М-4 10"5 2.717 3.777 2.317

2 0.5 160 х 120 0.025 5•10"4 10"5 2.760 3.818 2.306

2 0.2 160 х 120 0.025 5 • М-4 10"5 2.791 3.848 2.298

Теоретические значения, рассчитанные по (5): 3.418 8.204 4.000

3 1.0 160 х 120 0.025 10-4 10"5 3.457 8.089 3.900

3 0.5 160 х 120 0.025 10-4 10"5 3.525 8.204 3.879

Теоретические значения, рассчитанные по (5): 3.982 22.330 9.333

5 1.0 160 х 120 0.025 10-4 10"5 4.076 21.914 8.961

5 0.5 160 х 120 0.025 10-4 10"5 4.167 22.287 8.914

Теоретические значения, рассчитанные по (5): 4.402 2203.6 834.33

50 1.0 80 х 60 0.05 10"6 10-4 4.756 2181.1 764.35

50 0.5 80 х 60 0.05 5 • Ю-6 10-4 4.739 2211.8 777.87

ласть расчета разбивается на две подобласти [6]. На левой выходной границе и свободной верхней границе предполагается равенство нулю нормальных производных плотности, давления и компонент скорости.

В качестве начальных условий используются параметры невозмущенного потока. Рассматривается течение вязкого одноатомного газа твердых шаров 7 = 5/3, Рг = 2/3, Sc = 0.77, ш = 0.5, число Рей-нольдса Re = 104. Расчеты проведены на равномерных пространственных сетках для различных значений параметра а и чисел Маха Ма = 1.5, 2, 3, 5, 50. Стационарное решение находится методом установления.

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

Параметры торможения: плотность ps, давление ps, температура Ts (индекс «s» от англ. «stagnation») вычисляются по формулам:

= Р2 1 +

7 — 1 и

TS = T2 1 +

7 — 1 и

z 2

7-1

Ps

(5)

7Ps

где р2, иг2, Т2, С2 — давление, скорость, температура и скорость звука за ударной волной, которые мы получаем из условий Ренкина-Гюгонио [9]:

Р2 = Р1

(7 + 1) Mag 2 + (7 - 1) Ма^

2 '

Р2 =Р 1

27 MaI -7 + 1

7 + 1

Uz2 = Uzi

2 + (7 - 1) Mag (7 + 1) Ma I '

T2= 7

P2 P2

C2

здесь р1, иг 1, р\, Тх — параметры набегающего потока.

Для всех чисел Маха, включая Ма = 50, полученные в расчете значения параметров торможения соответствуют теоретическим значениям (5). Наиболее точные значения получаются при а = 0.5. Для чисел Маха 1.5 и 2 уменьшение шага сетки в два раза мало влияет на точность решения.

Изолинии плотности и линии тока для чисел Маха 1.5, 5 и 50 (а = 0.5) для установившегося течения приведены на рис. 1-3. Наглядно видно изменение формы и положения ударной волны, а также структуры течения с ростом скорости набегающего потока. На рис. 4 представлены профили плотности в ударной волне. Приведенные графики демонстрируют высокую точность расчета скачков уплотнения и отсутствие осцилляций решения при больших числах Ма.

Рис. 1. Изолинии плотности и поле течения (Ма= 1.5) 3

Рис. 2. Изолинии плотности и поле течения (Ма = 5)

0 12 3 4

2

Рис. 3. Изолинии плотности и поле течения (Ма = 50)

Течение в канале с уступом

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

Р

5,0-

0,5-1-1-1-1-.-1-1-1-1-1—

2,0 2,5 3,0 3,5 4,0

I

Рис. 4. Профили плотности в осесимметричном течении

Таблица 2 Расчет течения в канале с уступом

V Сетка Нх — Ьу ей Число шагов г ргшп ртах

1 120 х 40 0.025 0.001 4000 4 0.551 4.464

2 240 х 80 0.0125 0.0005 8000 4 0.377 4.553

3 480 х 160 0.00625 0.0001 40000 4 0.247 4.595

течения в канале е уступом [7, 8]. Сложная конфигурация ударных волн, формирующихся с течением времени в канале, служит известным тестом для оценки работоспособности методов высокого порядка точности для решения уравнений Эйлера и Навье-Стокса. В соответствии с работами [7, 8] задача решается в следующей постановке: длина канала — 3, ширина — 1, высота уступа, расположенного на расстоянии 0.6 от начала канала, равна 0.2. Рассматривается невязкий нетеплопроводный газ с показателем адиабаты 7 = 1.4. На входной (правой) границе задается однородный поток с числом Маха Ма = 3, на выходной (левой) предполагается равен-

1

^0.5

Рис. 5. Изолинии плотности для

ство нулю производных от всех параметров потока. На стенках канала и уступа задаются отражающие граничные условия. Расчеты проводились на трех равномерных пространственных сетках (табл. 2). При вычислении релаксационного параметра г коэффициент а = 0.3.

Для сопоставления с результатами работ [7, 8] на рис. 5 приведено распределение плотности на момент времени £ = 4 на сетке 240 х 80 (приведено 50 изолиний, расположенных эквидистантно). Картина течения соответствует данным [8], полученным с использованием кусочно-параболических схем третьего порядка точности по пространству, и результатам [7], где применяются методы расщепления первого, второго и третьего порядка точности.

Отчетливо прослеживается образование вторичных волн отражения от верхней стенки канала и верхней поверхности уступа. За волной разрежения над углом уступа плотность газа минимальна, вблизи контактного разрыва за тройной точкой над уступом плотность газа максимальна. Полученные в расчетах минимальное и максимальное значения плотности на разных сетках приведены в табл. 2. Согласно [81, ртах/р00 =4.33, рШт/Роо = 0.18. Видно хорошее согласие максимального значения плотности и уточнение величины при сгущении пространственной сетки.

Таким образом, на примере расчета двух типичных течений в работе протестирован предложенный в [6] новый численный алгоритм, основанный на КГД уравнениях. Алгоритм представляет собой явную по времени разностную схему первого порядка точности. Достоинствами метода являются простота реализации и высокая точность, которая, по мнению авторов, обеспечивается использованием регуляри-заторов специального вида, присутствующих в КГД уравнениях. Алгоритм может использоваться для расчета вязких и невязких сверхзвуковых нестационарных течений с высокими числами Маха, возника-

в канале (Ма = 3, 1 = 4)

ющих, например, при входе летательных аппаратов в атмосферу планет.

Работа выполнена при финансовой поддержке РФФИ (грант 01-01-00061).

Авторы благодарны Ю. В. Шеретову за сотрудничество.

Литература

1. Елизарова Т.Г., Четверушкин Б.Н. 11 ЖВМ и МФ. 1985. 25, № 10. С. 1526.

2. Шеретов Ю.В. Математическое моделирование течений жидкости и газа на основе квазигидродинамических и квазигазодинамических уравнений. Тверь, 2000.

3. Елизарова Т.Г., Шеретов Ю.В. // ЖВМ и МФ. 2001. 41, №2. С. 239.

4. Елизарова Т.Г., Соколова М.Е. // Вестн. Моск. ун-та. Физ. Астрон. 2001. №5. С. 19 (Moscow University Phys. Bull. 2001. N 5. P. 24).

5. Elizaroisa Т.О., Sheretov Yu.V. // J. La Houil. Blan. 2003. N 5. P. 66.

6. Шеретов Ю.В. // Применение функционального анализа в теории приближений. Тверь, 2001. С. 191.

7. Граур И.А. // ЖВМ и МФ. 2002. 42, №5. С. 698.

8. Woodward P., Colleta Р. // J. Comput. Phys. 1984. N 54. P. 115.

9. Ландау Jl.Д., Лившиц Е.М. Гидродинамика. М., 1986.

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

УДК 539.124.17

ТЯЖЕЛЫЕ МАЙОРАНОВСКИЕ НЕЙТРИНО В РОЖДЕНИИ ДИЛЕПТОНОВ НА ЛЕПТОН-ПРОТОННЫХ КОЛЛАЙДЕРАХ

А.Али*), A.B. Борисов, Д.В. Журидов

(.кафедра теоретической физики) E-mail: borisov@ave.phys.msu.su

Вычислены сечения глубоконеупругих процессов рождения пар лептонов одного знака заряда е+р —>• Ре£+£'+Х и v,;p —>• е£+£'+Х (£,£' = е. //,. т). обусловленных обменом тяжелым майорановским нейтрино. Рассмотрены возможности наблюдения этих процессов на существующем ер-коллайдере HERA и проектируемом суперколлайдере VLHC.

1. В настоящее время экспериментально установлено, что нейтрино обладают массой (все аспекты физики нейтрино, в том числе проблема массы и соответствующие экспериментальные данные обсуждаются в обзоре PDG [1]). Наиболее убедительными являются наблюдения тремя независимыми группами оецилляций нейтринных ароматов [2-4] (ссылки на более ранние работы см. в обзоре [5], где анализируются современные данные по солнечным, атмосферным, реакторным и ускорительным нейтрино). Указанные осцилляции объясняются тем, что нейтрино щ определенного аромата £ = е, р, т представляет собой когерентную суперпозицию состояний щ с определенными массами пц (см., напр.,

[6,7]):

'•'<• Ц1 (!)

i

где коэффициенты 17ц образуют матрицу лептонного смешивания.

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

^ А. Ali (DESY, Hamburg).

кварков и заряженных лептонов в стандартной модели [6, 7]. Майорановекое же нейтрино — истинно нейтральная частица, тождественная своей античастице. Майорановский массовый член в лагранжиане не сохраняет лептонное число, изменяя его на две единицы. Поэтому майорановские нейтрино могут приводить к ряду процессов с несохранением лептонного числа, поиск которых представляет одно из важных направлений современной экспериментальной физики элементарных частиц. Заметим, что нейтринные осцилляции нечувствительны к типу массы нейтрино. В течение многих лет ведется поиск безнейтринного двойного бета-распада ядер [1, 6], обусловленного обменом виртуальными май-орановскими нейтрино. Недавно появилось первое экспериментальное указание на его обнаружение [8], пока не получившее независимое подтверждение другими группами. К числу процессов, в которых рождаются пары лептонов одного знака заряда и тем самым проявляется майорановская природа нейтрино, относятся редкие распады мезонов типа К+ —> 7г^ р+р+ (см., напр., [9]) и глубоконеупругие адрон-адронные и лептон-адронные столкновения: рр е^Х [9, 10], е+р (■' X [11, 12],

^М^р-р+р+Х [13].