Система свободного доступа Galactica для расчетов взаимодействий N-тел

Размещено на http://www.allbest.ru/

Система свободного доступа Galactica для расчетов взаимодействий N-тел

Смульский И.И.

Аннотация

Система Галактика предназначена для численного решения проблемы взаимодействия N-тел. Она основана на методе высокой точности решения для решения дифференциальных уравнений движения тел, силы взаимодействия которых обратно пропорциональны квадрату расстояния. В статье приведены дифференциальные уравнения движения и метод их решения. Описаны структуры входного файла, файла начальных условий и файлов результатов решения. Представлены примеры этих файлов для разных решаемых задач. Объясняется работа исполняемых файлов. В статье описаны все теоретические и практические вопросы так, что даже начинающий исследователь может использовать систему Galactica в своей работе. галактика тело уравнение движение

Ключевые слова: гравитационные и кулоновские взаимодействия, уравнения движения, компьютерные вычисления, Солнечная система, свободный доступ.

Введение

Для исследования эволюции Солнечной системы за геологические периоды времени была создана программа Galactica. Она основана на высокоточном методе интегрирования дифференциальных уравнений движения материальных точек при их взаимодействии по закону тяготения Ньютона. С помощью программы Galactica был решен ряд задач [1], в том числе: эволюция Солнечной системы за 100 млн. лет, оптимальный полет космического аппарата к Солнцу, эволюция движения астероидов Апофис и 1950DA и превращение их орбит в орбиты спутников Земли; моделирование эволюции оси вращения Земли, воздействие вращающегося Солнца на орбиты планет и др.

Так как программа Galactica и накопленный набор дополнительных средств позволяет решать разнообразные задачи гравитационного взаимодействия тел, то была создана система свободного доступа Galactica. Кроме самой программы система Galactica включает дополнительные компоненты, необходимые для постановки задачи, подготовки ее к решению, решения задачи и анализа ее результатов. Система доступна по адресу: http://www.ikz.ru/~smulski/GalactcW/. В данной статье представлено руководство по ее использованию.

Применяемые в системе Galactica методы во многом отличаются [1] от методов, применяемых в небесной и космической динамике. Результатом этих отличий является более высокая точность решения задач и более широкие возможности по разнообразию постановок задач и их решению. Например, несмотря на то, что метод основан на взаимодействии материальных точек, программа Galactica позволила решить задачу об эволюции вращательного движения Земли и о воздействии вращающегося Солнца. В этих задачах Земля и Солнце уже рассматриваются не как материальные точки.

В программе Galactica решены вопросы высокоточного интегрирования дифференциальных уравнений и выдача информации о вычислительном процессе и его точности. В системе Galactica, в частности, в настоящем руководстве, много внимания уделено вопросам постановки задачи. Мы надеемся, что они позволят исследователям расширить круг задач, решаемых с помощью программы Galactica.

Для тех, кто понимает теоретическую механику, она является очень простой наукой. К сожалению, множество ответвлений в механике в 20^ом веке привели к разному толкованию ее основ. Поэтому в настоящем руководстве мы проводим основные положения механики, необходимые для решения задачи.

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

От первоначальной идеи до окончательного результата задача проходит ряд этапов. Для их успешного решения требуется много разнообразных знаний. Не каждый исследователь ими обладает. Чтобы облегчить работу исследователя, в руководстве большая часть этапов изложена, и приведены конкретные примеры. Одним из важных этапов является анализ полученных результатов. В дальнейшем, мы дополним Руководство этим параграфом с приложением программ по анализу траекторий.

По вышеуказанному адресу система Galactica доступна для работы на персональном компьютере. Для решения больших задач система Galactica создана на суперкомпьютере Сибирского суперкомпьютерного центра (SSCC) Сибирского отделения Российской Академии наук г. Новосибирск. В среде MPI разработан вариант программы Galactica для счета в параллельном режиме. Однако возможности SSCC ограничены, поэтому свободный доступ на суперкомпьютер до настоящего времени не решен.

Так как кулоновское взаимодействие подобно ньютоновскому, то мы разработали вариант программы Galactica для расчета взаимодействия заряженных частиц. С начала 20^ого века в науках о микромире возобладал статистический метод исследования поведения множеств частиц. Во многом этому способствовало то, что в то время невозможно было решать задачи детерминированного взаимодействия многих заряженных частиц. Сейчас решены теоретические задачи [2], [3], которые позволяют задать начальные условия для электронных оболочек атомов. Поэтому мы рассчитываем, что программа Galactica для кулоновского взаимодействия позволит в дальнейшем достигнуть более детерминированного представления об устройстве микромира. После того, как будет хотя бы одна из задач этой области решена, мы дополним систему Galactica кулоновским взаимодействием.

Система Galactica задумана как открытая система. Мы рассчитываем, что с появлением решенных задач другими исследователями, они также в соответствующих папках системы будут размещать свои файлы начальных условий, программы по их вычислению, а также др. разработки. С одной стороны, они позволят другим исследователям повторить решения задачи. Если результат повторится, то решением второго исследователя будет засвидетельствована достоверность выводов, полученных первым исследователем. С другой стороны, основываясь на достижениях каждого исследователя, последующие исследователи могут более успешно ставить и решать новые задачи.

В системе Galactica имеются все необходимые сведения, чтобы даже начинающий исследователь мог использовать ее в своей работе. Имеется неограниченное количество задач, которые могут быть решены с помощью системы Galactica. Приступить к этой работе могут студенты младших курсов, и даже способные ученики старших классов школы. За 2-3 года они могут выйти на высокий уровень и к окончанию университета решить серьезную научную проблему и внести свой вклад в понимание окружающего мира.

Содержание статьи организовано следующим образом.

1. Основные сведения о системе Galactica.

2. Дифференциальные уравнения движения тел как материальных точек.

3. Метод решения.

4. Описание входного файла galacf26.dat.

4.1. Общие сведения.

4.2. Параметры для области взаимодействия.

4.3. Общие параметры решения задачи.

4.4. Параметры для выдачи траекторий определенных тел.

4.5. Параметры для изображения результатов на экране дисплея.

5. Файлы начальных условий и результатов интегрирования.

6. Описание файлов траекторий 1t, 2t… traekt.prn.

7. Описание файла сближений dice.dat.

8. Выдача на экран дисплея.

9. Образцы входных файлов galacf26.dat.

10. Файлы начальных условий решенных задач.

11. Файлы исполняемых модулей.

12. Создание файла начальных условий.

13. Порядок решения задачи.

Литература.

Приложение 1. Образец программы в среде MathCad для создания файла начальных условий.

Приложение 2. Моды дисплея в стандартной графике Фортрана.

1. Основные сведения о системе Galactica

Система Galactica предназначена для решения задачи механического взаимодействия тел в виде материальных точек. Взаимодействие происходит по закону тяготения Ньютона. Файлы системы расположены в папке GalactcW. В ней имеется несколько папок: InCndFls, InpFiles и RunFiles и файлы GalDiscrp.pdf и GalDiscrE.pdf с описанием системы Galactica на русском и английском языках, соответственно. В архивированном виде папки можно скачать по адресу http://www.ikz.ru/~smulski/GalactcW/GalacFoldrs.zip.

В состав системы входят исполняемые модули с именем, например, glc3b30 - на суперкомпьютере в операционной системе Unix (glc3b30.exe - в системе DOS); входной файл galacf26.dat и файл исходных данных о взаимодействующих телах и начальных условий их движения с именем, например, sn69JPA2.dat. В дальнейшем с целью сокращения будем его называть как файл начальных условий.

Результаты работы исполняемого модуля выдаются в виде выходных файлов. Форма выдачи определяется ключами, которые задаются во входном файле galacf26.dat. Через определенное количество шагов могут выдаваться файлы с характеристиками движения тел, которые имеют имена: 1, 2, 3… garez28.dat. Файл garez28.dat выдается всегда по завершению. В нем содержатся характеристики тел на последнем шаге счета. Это основной вид выдачи результатов решения задачи.

Кроме того, при задании во входном файле galacf26.dat определенных ключей выдаются параметры траекторий тех тел, номера которых также заданы во входном файле. Имена файлов траекторий: 1t, 2t, 3t … traekt.prn. Размеры файлов определяются количеством данных, которые задаются во входном файле.

При задании определенных ключей выдается файл dice.dat. В нем имеется время сближения, расстояние сближения и номер тела, с которым произошло сближение заданного тела. Каждая строка в файле dice.dat формируется за определенное количество шагов. Их число задано в файле galacf26.dat.

При определенном ключе выдается файл в двоичном коде (или в зависимости от вида машины в неформатированном виде) с именем galgon4. В нем сохранено все состояние задачи в двоичном виде на последнем шаге счета. При повторном запуске исполняемого модуля счет продолжается с этого шага. При этом данные с входного файла и файла начальных условий не считываются. Эта процедура необходима при решении задачи с большими временами счета. Перерывы в счете задачи в этом случае не ухудшают точность ее решения.

В случае появления нештатных ситуаций создается файл ошибок ErrGlctca, в котором приводится их описание. Эти ошибки связаны с открытием файлов, с их считыванием, недостаточной размерностью массивов и невозможностью обеспечить необходимую точность при расчете слияния тел в случае их столкновения.

В программе Galactica также имеется выдача результатов на дисплей персонального компьютера. Эти общие для всех тел параметры задачи позволяют следить за ходом ее решения и контролировать точность счета. При решении задачи на суперкомпьютере они сохраняются в файле выходного потока.

В варианте программы Galactica с графикой на экране дисплея персонального компьютера наряду с общими параметрами задачи выдаются изображения тел с их векторами скорости. Вид изображения определяется ключами и параметрами во входном файле galacf26.dat.

2. Дифференциальные уравнения движения тел как материальных точек

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

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

Согласно закону всемирного тяготения тело с номером k притягивает тело с номером i, и это притяжение выражается силой:

, (1)

где G - гравитационная постоянная;

- радиус-вектор от тела с массой m_k до тела с массой m_i.

Если количество тел равно n, то на i-e тело остальные будут оказывать воздействие с суммарной силой

. (2)

Под воздействием, выраженным силой (2), в соответствии со вторым законом механики , i-ое тело будет двигаться относительно инерциальной (неускоренной) системы отсчета с ускорением

i = 1,2,…,n, (3)

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

Соотношение (3) представляет систему 3n нелинейных дифференциальных уравнений второго порядка, где п может быть любым. Количество тел п задается во входном файле galacf26.dat. Кроме того, массивы исполняемого модуля не должны иметь размерность меньшую п. Например, модуль glc3b30 позволяет рассчитывать взаимодействие включительно до п = 30 тел.

Дифференциальное уравнение (3) представлено в неускоренной (инерциальной) системе координат. Однако все тела движутся с ускорением, и только их центр масс C в системе взаимодействующих тел является неускоренным (при условии, что пренебрегается воздействие других тел на рассматриваемую систему тел). Как правило, из наблюдений получают координаты и скорости относительно какого-либо тела, с которым связана система координат. Например, в Солнечной системе для тела с массой они могут быть относительно Солнца (S): координаты x_Si, y_Si, z_Si и скорости v_xSi, v_ySi, v_zSi. С их помощью определяются координаты и скорости центра масс Солнечной системы, которые в проекции на ось x имеют вид:

X_C x_Si / M_Ss; V_xC v_xSi / M_Ss,(4)

где n - количество тел Солнечной системы; M_Ss - масса Солнечной системы, а в общем случае M_Ss - масса всех взаимодействующих тел.

Следует отметить, что движение выбранного тела, в данном случае Солнца, со временем изменяется, а система отсчета продолжает движение без ускорения. Поэтому систему отсчета связывают с определенным моментом времени в движении избранного тела. Например, систему отсчета, связанную с плоскостью экватора Земли, относят к началу 50-летнего периода: система 1900.0 г., 1950.0 г. или 2000.0 года.

Координаты и скорости тел относительно центра масс C в проекции на ось x будут:

x_Ci = x_Si - X_C; v_xCi = v_xSi - V_xc.(5)

Здесь координаты выражены в метрах, время в - секундах, а массы - в килограммах.

Уравнения (3) при отнесении переменных к характерному размеру системы взаимодействующих тел Am и их массе M_Ss приводятся к безразмерному виду [4-5]. Для Солнечной системы величину Am мы определяем так, чтобы безразмерное время T было в сидерических столетиях. При этом величина Am получилась близкой к полуоси орбиты Плутона. Уравнения в безразмерном виде для проекции x выглядит так:

,(6)

где x_i = x_Ci / Am - безразмерная координата i-того тела;

m_ok = m_k / M_Ss - безразмерная масса k-того тела;

r_ik = - безразмерное расстояние между i-тым и k-тым телом;

v_xi = v_xCik_v - безразмерная скорость i-того тела;

- коэффициент скорости в сек/м;

T = tk_t - безразмерное время;

- коэффициент времени в сид. ст./сек.

В файлах начальных условий в задачах для Солнечной системы безразмерное время T = tk_t - выражено в сидерических столетиях по 36525.636042 дней в столетии.

В безразмерных единицах используются все величины в программе Galactica, в том числе радиусы тел (Ra_o), их количества движения (P_ox), моменты количества движения (M_ox), кинетическая (E_o) и тепловая (E_to) энергии. При задании двух масштабных величин Am и M_Ss могут быть найдены коэффициенты k_v и k_t. После этого результаты решения могут быть пересчитаны к размерным величинам по следующим формулам:

t = T/k_t; x_mi = x_i Am; v_mi = v_i/k_v; Ra_i = Ra_oi Am; P_x = P_ox M_Ss/k_v;

M_x = M_ox M_SsAm/k_v; E = E_o M_Ss/k_v²; E_t= E_to M_Ss/k_v². (7)

В случае, если исходной масштабной величиной было время T, например в сидерических столетиях, то масштабная длина Am из выражения для k_t определяется по формуле

Am = (G M_Ss/k_t²)^1/3. (8)

В режиме Kl1 = 0, т.е. при работе программы Galactica с заданием исходных и начальных условий на основании равномерно распределенной области вещества используется плотность тел _b в кг/м³. По безразмерным радиусу Ra_oi и массе m_oi тел при заданной масштабной величине Am рассчитываются их размерные величины

Ra_i = Ra_o Am; m_i = m_oi3 Ra_i³_b/4. (9)

Затем рассчитывается вся масса тел M_Ss и коэффициенты k_v и k_t. Дальнейший расчет выполняется по формулам (7).

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

3. Метод решения

Дифференциальные уравнения движения тел (6) в программе Galactica интегрируются следующим методом. Значение функции в следующий момент времени T=T₀ + T определяется с помощью ряда Тейлора, который, например, для координаты x имеет вид:

, (10)

где x₀^(k) - производная порядка k в момент T₀.

Значение скорости x' = v_x определяется по аналогичной формуле, а ускорение x₀'' = - по формуле (6). Более высокие производные x₀^(k) определяются по аналитическим выражениям, которые мы вывели в результате дифференцирования уравнений (6). Сейчас используется расчетная схема шестого порядка, т.е. при K=6. При этом шестая производная определяется по разности пятых производных. Добавление каждой производной, за исключением шестой, увеличивало точность решения на 3 порядка. Шестая производная увеличивает точность на 2.5 порядка. Этот метод обеспечивает высокую точность и хорошую ее динамику. Например, при двойной длине числа (17 десятичных знаков) погрешность не увеличивается с продолжением времени счета.

При разработке метода было использовано более десяти различных способов контроля точности [7-8]. Ряд из них используется в программе Galactica для слежения за погрешностями, и сведения о точности выдаются в файлах результатов работы программы. Конкретные данные о точности приведены в [7-8] и др. опубликованных наших работах, а также даны ниже в других параграфах. По результатам решения задачи пользователь может видеть точность полученных результатов и в случае необходимости может ее изменить.

Алгоритм программы Galactica в настоящем Описании не приводится, но некоторые его элементы по мере их привлечения для объяснения ее результатов даны в последующих пунктах.

4. Описание входного файла galacf26.dat

4.1. Общие сведения

В табл.1 приведены параметры входного файла galacf26.dat, причем дополнительно в столбце 1 даны их номера. Все параметры можно подразделить на 4 группы: 1) параметры для области взаимодействия; 2) общие параметры решения задачи; 3) параметры выдачи траекторий определенных тел; 4) параметры для изображения результатов на экране дисплея.

Числами с точками в 4-ом столбце Табл. 1 приведены действительные величины, а без точки - целые. Действительные величины используются в формате F11.4, целые Klmax, Ltk, Ltk2, k2, Kl3, Kl11 и Kli - в формате I12, а остальные - в формате I5. Строковая переменная используется в формате A12.

В первой строке второго столбца Табл. 1 приводится дата варианта создания файла, например, 01.02.2012.

4.2. Параметры для области взаимодействия

В 1-ой группе с номерами строк 1ч20 и 22 приведены параметры для задания исходных данных и начальных условий для взаимодействующих тел, которые получаются из равномерно распределенного вещества в некоторой области пространства. Область задается в строках 1-3 Табл. 1 в виде параллелепипеда со сторонами A, B и C. Размеры их относительны и отнесены к наибольшей стороне A, т.е. A = 1. Начало этой области имеет координаты X0, Y0 и Z0. Вдоль оси X область разделена на Mu промежутков. Размер области B или C разделен на одинаковые по длине промежутки, которые близки к промежутку между делениями стороны A. Масса всей области определяется ее объемом ABC.

Таблица 1. Содержание входного файла

No	Описание параметра	Пара-метр	Величина, начиная с 60 позиции
1	Area length along axis X (A = 1) 01.02.2012	A	1.
2	Area height along axis Y (B <=1)	B	1.
3	Area thickness along axis Z (C <=1)	C	1.
4	Initial X coordinate of area	X0	0
5	Initial Y coordinate of area	Y0	0
6	Initial Z coordinate of area	Z0	0
7	Number of divisions of axis X	Mu	2
8	Relative initial density of the area	Ro0	1.
9	Relative initial density of the point body	RoS	2.
10	X - coordinate of the separate body	Xs	1.5
11	Y - coordinate of the separate body	Ys	0.5
12	Z - coordinate of the separate body	Zs	0.5
13	Velocity of the separate body along axis X	Us	-0.4
14	Velocity of the separate body along axis Y	Vs	-0.15
15	Velocity of the separate body along axis Z	Ws	0.5
16	Angular velocity of rotation of area	wu	0.0
17	Relative X-coordinate of the axis of rotation	X01	0.5
18	Relative Y-coordinate of the axis of rotation	Y01	0.5
19	Absolute average density of the area, kg/m^3	Roa	1.E-7
20	Absolute density of bodies, kg/m^3	Rob	5000.
21	Depiction with respect to first body (below)	Kl10	0
22	Step time integration at Kl1=0	dt	5.0e-4
23	Maximum dimension of TR array	Klmax	6000
24	Number of integration steps in one cycle	Ltk	1000
25	Number of cycles	Ltk2	1
26	Key input data from file 4	Kl1	2
27	Key start of counting time	Kl4	1
28	Search of encounters of body No. at Kl3	Nbd	0
29	Key of writing data to file	Kl7	0
30	Number of interacting bodies	k2	12
31	Steps interval of data (pictures) output	Kl3	500
32	Beginning step-by-step of data (pctrs) output	Kl11	0
33	Key of writing of the trajectory parameters	Kl12	0
34	Steps interval of trajectory parmtrs output	Kli	50
35	Display mode	Nre	0
36	Correction of the area size on the screen	Kl14	1
37	Coefficient of area form on the screen	af	0.72
38	Coefficient of mass size of the screen	Rd	0.08
39	Coefficient of velocity on the screen	amv	3.
40	Crrctn factor of step dtk at bodies' mergence	bt	0.8
41	Flat depiction of layer to Z	Ke	0
42	Rotation of area relative to X-axis by angle	al	10.
43	Rotation of area relative to Y-axis by angle	bet	20.
44	Begining of screen area along x	Xe0	300.
45	Begining of screen area along y	Ye0	280.
46	Image size in x	aix	150.
47	Image size in y	biy	70.
48	The initial data file name	FIMA4	sn2010Ap.dat
49	Number issued bodies for their trajectories	Np	2
50	Key of output of the body velocity	Kl15	0
51	Velocity colour of the first body	ipp(1)	3
52	Free key of the first body	jpp(1)	1
54	Number of the first body	kpp(1)	1
55	Velocity colour the second body	ipp(2)	5
56	Free key of the second body	jpp(2)	1
57	Number of the second body	kpp(2)	2
58	Velocity colour of the third body	ipp(3)	8
59	Free key of the third body	jpp(3)	1
60	Number of the third body	kpp(3)	4

Кроме тел, полученных разбиением области на одинаковые объемы, можно задать отдельное тело. Его масса определяется относительной плотностью RoS в строке 9. Величина RoS отнесена к относительной плотности области Ro0, т.е. последняя задается как Ro0 = 1.0. Параметрами Xs, Ys и Zs задаются относительные координаты отдельного тела. Если отдельное тело находится на расстоянии, большем промежутка между уже имеющимися телами, то это тело добавляется к общему количеству тел.

Радиусы тел рассчитываются по их объемам, абсолютной плотности Rob, исходя из их сферической формы. Координаты тел задаются в узлах, а их скорости приравниваются нулю. Скорость отдельного тела задается параметрами Us, Vs и Ws. Массы всех тел приводятся к массе всей системы.

Для задания начальных скоростей тел имеется возможность задать вращение всей области с угловой скоростью wu. Положение оси вращения задается координатами X01 и Y01. Ось вращения параллельна оси Z.

Для интегрирования этих задач шаг интегрирования задается в строке 22 параметром dt. Вся система начальных условий для этих тел создается в системе центра масс. Однако при введении отдельного тела система тел уже не приведена к центру масс. Поэтому в процессе их взаимодействия в системе может появиться направленное движение.

4.3. Общие параметры решения задачи

Эта группа параметров представлена в строках Табл. 1: 21, 23ч32, 40, 48.

Начнем со строки 24, а строку 23 рассмотрим ниже. Ltk - число шагов интегрирования. Оно не может превышать максимальное значение целого числа двойной точности 2147483647. В случае необходимости, большее число шагов задается параметром Ltk2>1. Это число повторения циклов с Ltk шагами.

Ключ Kl1 имеет пять значений: 0, 1, 2, 3 и 4. При Kl1 = 0 не используется файл начальных условий, а все данные задачи определяются параметрами, описанными в п. 4.2. При Kl1 = 1 исходные и начальные условия считываются из файла результатов garez28.dat. В этом случае идет продолжение предыдущего решения задачи. Однако, при этом точность уменьшится за счет округлений при переходе от двоичных чисел к десятеричным в момент создания файла garez28.dat. В большинстве случаев такое уменьшение точности несущественно. Перед запуском задачи в файле garez28.dat можно изменить шаг dT. Это позволяет с большей или меньшей точностью просчитать определенные этапы задачи.

При Kl1 = 2 исходные и начальные условия считываются из файла, имя которого вводится в строке 48 параметром FIMA4. Следующие два значения 3 и 4 ключа Kl1 включают режим продолжения счета без изменения точности. В этом случае после каждого этапа счета, который определяется ключами Kl7 и Kl3, создается неформатированный (двоичный) файл galgon4. При новом запуске задачи игнорируется чтение входного файла и файла начальных условий, а все условия задачи считываются с файла galgon4. При необходимости некоторые условия задачи можно изменить. Для этого разработана дополнительная программа, которая в данный набор программ не входит. Значения ключей отличаются тем, что при первом запуске задачи при Kl1 = 3 решение начинается с файла garez28.dat, а при Kl1 = 4 - с имени файла, которое определено параметром FIMA4.

При первоначальном счете задачи с ключами Kl1 = 2 или 4 первый шаг уменьшается в 1000 раз, т.е. равен 0.001dT.

При значении ключа Kl4 = 1 в строке 27 продолжается счет времени T со значения, данного в файле начальных условий, а при Kl4 = 0 счет времени начинается с T = 0.

Параметр Nbd указывает номер тела, для которого будут запоминаться сближения с телами на минимальные расстояния на интервале шагов Kl3. В этом случае создается файл dice.dat, в который после Kl3 шагов заносится характеристики сближения: время, расстояние и номер тела. При Nbd = 0 сближения не ищутся, и файл dice.dat не создается.

При значении ключа Kl7 = 1 в строке 29 через Kl3 шага происходит выдача результатов интегрирования задачи на промежуточных интервалах. Они формируются в файлах 1, 2, 3 … garez28.dat. Общее число файлов определяется как LtkLtk2/Kl3. При запуске задачи на счет необходимо удостовериться, что число файлов не будет запредельно большим. Промежуточной выдачи результатов не происходит при ключе Kl7 = 0.

Параметр k2 в строке 30 определяет число взаимодействующих тел. Если массивы исполняемого модуля имеют размерность меньшую k2, на экран дисплея и в файл ошибок ErrGlctca будет выдано сообщение «*08G Change the dimension of arrays *** Change l2 >= k2», где k2 - заданное число. В этом случае нужно перейти на исполняемый модуль с большим числом тел. Следует также иметь ввиду, что из файла начальных условий будет считываться количество параметров для k2 тел, независимо от того, для скольких тел они приведены в файле.

В строке 31 параметр Kl3 определяет число шагов, через которое происходит выдача результатов интегрирования в файлы: 1, 2, 3 …, а также на экран дисплея персонального компьютера (в файл выходного потока суперкомпьютера). Необходимо иметь ввиду, что большая частота выдачи результатов уменьшает скорость счета программы Galactica.

Ключ Kl11 в строке 32 задает номер шага, начиная с которого в течение 10 шагов dT будет идти пошаговая выдача результатов. С помощью этого ключа предоставляется возможность детально изучить ситуацию, которая обнаружилась в отделенном от начала счета месте. При Kl11 = 0 эта опция отключается.

В строке 40 задается коэффициент коррекции bt шага dT при слиянии тел. Его целесообразно варьировать в ситуациях, когда при сближении тел расчет их слияния затрудняется.

В строке 48 параметр FIMA4 определяет имя файла исходных данных и начальных условий. Имя файла может быть задано строковой переменной с символами не более 12.

4.4. Параметры выдачи траекторий определенных тел

Для исследования траекторий определенного тела необходимо иметь его кинематические характеристики с высоким разрешением по времени. Для этой цели предназначены параметры 33-34 и 49-60. Ключ Kl12 = 1 в строке 33 включает опцию записи параметров траектории в файл. При этом выбираются кинематические данные определенных тел через заданное число шагов и после накопления их количества, определенного параметром Klmax, выводятся в файлы траекторий 1t, 2t, 3t … traekt.prn. Параметром Kli определяется число шагов по времени dT выборки кинематических характеристик тел.

Параметр Np в строке 49 определяет количество тел, для которых выводятся характеристики траекторий. А ключ Kl15 = 1 в строке 50 включает опцию выдачи компонент скорости. При Kl15 = 0 выдаются только координаты тел.

В последующих строках тройками параметров ipp(), jpp(), kpp() задается каждое тело. Номер тела определяет параметр kpp(). Параметр ipp() определяет цвет вектора скорости при использовании исполняемого модуля программы Galactica с графическим режимом. Параметр jpp() - не используется в настоящей версии программы. Цифрами в скобках при вышеупомянутых параметрах даны номера последовательности тел в файлах траекторий 1t, 2t, 3t … traekt.prn. Величина этих файлов определяется Klmax. В представленных компиляциях исполняемых модулей величина Klmax не должна превышать 15000. Для этих же модулей число тел Np не должно превышать 3 при Kl15 = 1 и не должно превышать шести тел при Kl15 = 0. Для задания дополнительных трех тел в файле galacf26.dat добавляются три строки с параметрами ipp(), jpp() и kpp()для каждого тела.

В случае неверного задания вышеупомянутых параметров в файл ошибок ErrGlctca и на экран дисплея выдается сообщение «*05G Change the dimension of TR *** Kl12d0 = 750, Kl12e0 = 20.» Здесь 750 и 20 - заданные размеры двумерного массива для параметров траекторий TR(Kl12d0, Kl12e0), где Kl12d0 - количество записей, а Kl12e0 - количество чисел в одной записи.

4.5. Параметры для изображения результатов на экране дисплея

Строками 21, 35-39, 41-47 в Табл. 1 задаются параметры для изображения тел и векторов их скоростей (см. рис. 1) на экране. Эти параметры используются при работе исполняемых модулей с графикой.

Используется система координат та же, что и в файле начальных условий, т.е. связанная с центром масс всей системы. Если необходимо представить движение тел относительно какого-либо тела, то в строке 21 для ключа Kl10 задается значение 1. В этом случае изображение координат и скоростей всех тел выдается относительно тела, номер которого задан параметром kpp(1) в строке 54.

Интервал шагов выдачи изображений задается, как уже упоминалось, ключом Kl3 в строке 31. Параметр Nre в строке 35 используется для учета моды дисплея. Моды дисплея в стандартной графике Фортрана приведены в Приложении 2. Влияние моды зависит от особенностей компьютера и исполняемого модуля. Необходимое значение может быть выбрано из представленных в Приложении 2. Значение Nre = -3 является, по-видимому, наиболее приемлемым.

Параметр Kl14 в строке 36 для обеспечения видимости тел при неограниченном расширении системы в данный версии программы Galactica не используется.

Параметр af в строке 37 позволяет сохранить форму изображения, например, квадрат представить как квадрат, а не как прямоугольник. При изображении тела и скорости их размеры на экране определяются с помощью отнесения к максимальным значениям массы и скорости. Это позволяет их видеть на экране при изменении масс и скоростей тел в широких диапазонах при решении задачи. Для первоначального изображения тела на дисплеях разной конфигурации используются корректирующие коэффициенты Rd и amv в строках 38 и 39.

Если система тел представлена многослойной областью вдоль оси Z, то с помощью параметра Ke можно задать номер слоя, тела которого будут выведены на экран. При Ke = 0 эта опция отключена.

В строках 42 и 43 задаются повороты системы тел при изображении относительно осей X и Y с помощью параметров al и bet, соответственно. Эти величины задаются в градусах.

В строках 44 и 45 с помощью параметров Xe0 и Ye0 задается верхнее левое начало области изображения на экране. Для изменения масштаба изображения вдоль осей X и Y используются параметры aix и biy, соответственно. Чтобы увеличить изображение, необходимо увеличить aix и biy.

Как уже ранее отмечалось, вектора скоростей отдельных тел можно выделить определенным цветом. Номер тела задается параметрами kpp() в строках 54, 57 и т.д., номер цвета - параметрами ipp() в строках 51, 55 и т.д. В стандартных графических операторах Фортрана соответствие между видом цвета и его номером представлено файлом FrtnColr.bmp в папке InpFiles.

Вид изображения зависит также от параметра Mu. Поэтому желательно задавать значение этого параметра, а также параметры B и C в 2-ой и 3-ей строках такими, какими они имеются в файле начальных условий. Позиции этих параметров описаны в п. 5.

Форма изображения результатов интегрирования зависит от многих факторов, в том числе от версии Фортрана и компилятора, от операционной системы компьютера и от вида экрана и его разрешения. В папке InCndFls с образцами файлов начальных условий решенных нами задач приведены изображения тел в файлах с расширением «.bmp». Имена этих файлов совпадают с именами файлов начальных условий. Например, для файла начальных условий sn49jplc.dat изображения тел приведены в файле Sn49jplc.bmp.

Изображения тел также имеются в папке InpFiles\Area. Здесь представлено 6 образцов входного файла galacf26.dat, в именах которых вместо расширения dat приведены другие символы. Они означают тип задачи, в которой рассматривается взаимодействие тел, образованных из области вещества. Изображение тел этой задачи в начальный момент времени представлено файлами с расширением bmp, а их имена представлены символами из имени входного файла. Например, входной файл galacf26.64r представляет взаимодействие 64 радиально движущихся тел, а в файле 64r.bmp представлено их изображение.

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

5. Файлы начальных условий и результатов интегрирования

Файлы начальных условий, например: sn49f.dat, sn49jplc.dat и др. в папке InCndFls, а также выходные файлы: 1, 2, 3 и т.д. и файл garez28.dat имеют одинаковую структуру. При описании файлов будем использовать два вида обозначения данных: 1) прямым шрифтом без индексов, и 2) курсивом. Первый вид применяется в программах, а второй вид - в печатных материалах.

Каждый из вышеупомянутых файлов начинается из набора 24 чисел: T0, omm, Um, dTp, Px, Py, Pz, AMx, AMy, AMz, Spsx, Spsy, Spsz, E, Em, Ett, dT, i2b, j2b, k2b, A, B, C, Mu. Эи числа являются общими для всех взаимодействующих тел. Затем следуют пятнадцать чисел для каждого k-ого тела: om(k), (X(k,q),q=1,3), (U(k,q),q=1,3), (dUp(k,q),q=1,3), (Sp(k,q),q=1,3), Ra(k), Et(k). Здесь q - индекс перечисления переменных по трем координатам x, y, z.

Общие параметры представляют собой следующее:

Т0 - время в относительных единицах, например, в файлах для Солнечной системы в сидерических столетиях по 36525.636042 дней в столетии;

omm - максимальная масса тела, отнесенная к массе всех тел;

Um - максимальная скорость в относительных единицах;

dTp - шаг по времени в относительных единицах на предыдущем шаге счета;

Px, Py, Pz - проекции количества движения всей системы тел на оси координат;

AMx, AMy, AMz - проекции момента количества движения всей системы тел;

Spsx, Spsy, Spsz - проекции суммы моментов количеств движений от собственного вращения тел вокруг своих осей (проекции спинов);

Е - суммарная кинетическая энергия всех тел в данный момент Т0;

Em - максимальная кинетическая энергия всех тел от начала интегрирования;

Ett -тепловая энергия всех тел, которую они приобрели при их слиянии в результате столкновения;

dT - текущий шаг по времени в относительных единицах;

i2b, j2b, k2b - количества узлов в результате деления массива с веществом по осям x, y и z, соответственно;

A, B, C - размеры этого массива по осям x, y, z;

Mu - число промежутков деления массива по оси x.

Последние семь параметров используются при задании взаимодействующих тел в виде области равномерно распределенного в пространстве вещества с плотностью Ro0.

Далее следуют 15 параметров для каждого тела:

om(k) m_ok - безразмерная (относительная) масса k = i-того тела (k и i - номера тел в разных документах);

(X(k,q),q=1,3) x_i, y_i, z_i - безразмерные барицентрические координаты k = i-тела; для файлов Солнечной системы в экваториальной системе координат;

(U(k,q),q=1,3) v_xi, v_yi, v_zi - безразмерные скорости k = i-того тела;

(dUp(k,q),q=1,3) x_i⁽⁵⁾, y_i⁽⁵⁾, z_i⁽⁵⁾ - безразмерные производные пятого порядка k = i-того тела;

(Sp(k,q),q=1,3) S_pxi, S_pyi, S_pzi - безразмерные моменты количества движения (спины) от собственного вращения k = i-того тела. Спины появились у тела, если оно образовалось слиянием других тел в результате их сближения на расстояние, меньшее суммы их радиусов. Изначально вращение тел можно не учитывать и их спины задавать S_pxi=S_pyi=S_pzi=0. Так сделано в приведенных ниже образцах файлов начальных условий.

Ra(k) - безразмерный эквивалентный радиус k = i- того тела; он рассчитывается по массе и средней плотности тела при его сферической форме;

Et(k) - E_ti - безразмерная тепловая энергия k = i-того тела. Она возникла при слиянии нескольких тел в одно за счет перехода кинетической энергии в тепловую.

Изначально тепловая энергия тел может не учитываться, т.е. E_ti=0, например, в приведенных образцах файлов начальных условий E_ti=0. При слиянии тел с массами m_oi и m_ok тепловая энергия слившегося тепла определяется по формуле:

, (11)

где m_oi m_ok , а v_rki - безразмерная радиальная скорость тела с массой m_ok относительно тела с массой m_oi.

После слияния масса i-того тела увеличивается: m_oi + m_ok, а масса k-того тела приравнивается нулю. Изменения координат, проекций скорости и спина, радиуса тела, которые вызваны слиянием, пересчитываются по алгоритму (см. параграф 11.3 книг [4-5]).

В последней строке файла начальных условий могут быть даны пять чисел: дата создания начальных условий, количество тел, масса системы, масштабный размер в метрах и коэффициент скорости: дата, k2, M_Ss, Am, k_v. Причем дата приведена целым числом из восьми чисел: год, месяц и день. Эта строка не считывается программой Galactica и отсутствует в файлах результатов.

Образец программы в среде MathCad для создания файла начальных условий приведен в Приложении 1.

6. Описание файлов траекторий 1t, 2t… traekt.prn

Как уже упоминалось в п. 4.4 имена файлов траекторий могут быть: 1t, 2t, 3t … traekt.prn. Все они имеют одинаковый состав, за исключением завершающего файла traekt.prn. Каждая запись относится к одному моменту времени и содержит Kl12e чисел Их количество определяется выражением

Kl12e = 3(Kl15 + 1)Np + 2, (12)

где Kl15 и Np - параметры, заданные во входном файле galacf26.dat и описанные в п. 4.4.

Начинается запись с момента времени T в безразмерных единицах, затем следует E - кинетическая энергия всей системы. Затем для каждого тела выводятся три координаты: x, y, z и при ключе Kl15 = 1 проекции скорости vx, vy, vz. Все эти параметры, как и время, даны в безразмерных единицах (см. п. 2).

Количество записей в файлах и их объем определяется параметрами, описанными в п. 4.4. При небольшом количестве записей создается только один файл траекторий traekt.prn. В его конце производится запись общих 10 параметров задачи: T, Tem, Em, Ett, AMz1, Nb, Lsls, Um, omm, Amz,

где T - время конца счета задачи;

Tem - момент времени с наибольшей кинетической энергией системы;

Em - максимальная кинетическая энергия системы;

Ett - тепловая энергия системы в конце счета задачи;

AMz1 - момент количества движения всей системы относительно оси z на первом шаге счета;

Nb - количество тел в системе в начале счета;

Lsls - количество слияний, т.е. количество тел, которые при сближении с другим телом были объединены в одно тело;

Um - наибольшая скорость за все время счета;

omm - масса наибольшего тела;

AMz - момент количества движения всей системы относительно оси z в конце счета.

В случае выдачи предшествующих файлов, например 1t, файл traekt.prn имеет такой же объем, и после завершающей записи в нем содержатся записи, оставшиеся от предшествующего файла.

Кинетическая энергия системы в программе Galactica рассчитывается по формуле

. (13)

Как уже отмечалось, расчет тепловой энергии при слиянии двух тел проводится по формуле (11). Как будет показано далее в п. 8 момент количества движения AMz = M_z определяется по формуле (14).

7. Описание файла сближений dice.dat

Как уже упоминалось в п. 4.3, при ненулевом параметре Nbd происходит расчет расстояний сближений тела с номером Nbd с другими телами. После завершения Kl3 шагов счета в файл dice.dat записывается три числа: время сближения на наименьшее расстояние с определенным телом, номер этого тела и величина расстояния. Количество таких записей в этом файле определяется как LtkLtk2/Kl3.

Эта опция использовалась нами для исследования сближений астероидов с планетами и Луной [9-10]. На рис. 1 и рис. 7 работы [10] представлены результаты, которые получены на основании данных файла dice.dat.

8. Выдача на экран дисплея

При решении задачи на персональном компьютере программой Galactica с графической опцией на экран дисплея через определенное число шагов интегрирования уравнений выдается численная информация о процессе интегрирования и картина положения планет с векторами скорости (см. рис. 1). Эти данные позволяют наблюдать процесс решения и качественно определять его достоверность. При решении программой без графики на экран выдается только числовая информация, а при решении на суперкомпьютере эта информация поступает в файл выходного потока.

Первые 16 параметров повторяют значения общих параметров:

TО, omm, Um, dTp, Px, Py, Pz, AMx, AMy, AMz, Spsx, Spsy, Spsz, E, Em, Ett,

файлов начальных условий и результатов интегрирования, описанных в п. 5. Оставшиеся три параметра представляют номер шага интегрирования L_t, его величину dT и изменение дM_z момента количества движения в проекции на ось z (следует отметить, что полный номер шага интегрирования состоит из произведения Ltk·Ltk2).

Рис. 1. Представление на экране дисплея (вид с южного полюса) результатов численного интегрирования системы уравнений (6). Курсивом даны пояснения и расшифровка числовых параметров.

Эти величины позволяют легко анализировать эволюцию системы и процесс интегрирования уравнений (6). Каждая из этих величин: наибольшая масса omm = m_max, количество движения P, момент количества движения M, сумма собственных моментов количества движения S, кинетическая E, тепловая энергия E_t всей совокупности тел и относительное изменение (дM_z) проекции момента М на ось z, постановщику задачи позволяет получить представление о ходе решения, а в случае ошибки - локализовать ее место нахождения.

Одним из важных факторов контроля является относительное изменение момента дM_z. При отсутствии внешних воздействий момент количества движения системы, например в проекции на ось z,

, (14)

не должен изменяться. Поэтому относительное изменение момента

дM_z = (M_z - M_z0)/M_z0, (15)

где M_z0 - величина момента количества движения в начальный момент времени, должно быть равно нулю. Если его значение не равно нулю, то это свидетельствует о погрешностях при численном интегрировании задачи. О показателе точности дM_z и связи его с погрешностью координат и скоростей более детально дано в работах [7-8].

Отличие от нуля составляющих количества движения всей системы тел Px, Py, Pz может информировать о появлении погрешностей в каком-то направлении. Необычно высокие значения максимальных массы omm и скорости Um постановщику задачи также подскажут вероятную причину погрешности.

9. Образцы входных файлов galacf26.dat

В папке InpFiles представлены образцы входных файлов: в самой папке и входящей в нее папке Area. Входной файл galacf26.dat предназначен для решения осесимметричной проблем 3000 тел с файлом начальных условий axsymb3t.dat. Пример входного файла galacf26.v01 использовался в задаче о многослойных кольцевых структурах [3]. В файле FrtnColr.bmp представлены номера цветов от 1 до 12 в стандартной графике Фортрана. Эти цвета можно задавать для изображения векторов скорости определенных частиц. Как уже отмечалось ранее, цвета задаются параметрами ipp() во входном файле.

В папке Area представлено 6 примеров входного файла galacf26.dat при задании начальных условий для взаимодействующих тел, полученных разделением области вещества на определенное число тел. Вместо расширения «dat» имена этих файлов содержат символы, которые идентифицируют получаемую систему тел. Для каждого входного файла приведен графический файл с расширением bmp, имя которого начинается с этих символов. Ниже перечислены задачи, которым соответствуют эти файлы:

1) galacf26.3br - радиальное взаимодействие трех тел;

2) galacf26.3b1 - три тела при задании угловой скорости вращения первоначальной области;

3) galacf26.4bo - четыре тела при задании угловой скорости вращения;

4) galacf26.64r - радиальное взаимодействие 64 тел;

5) galacf26.64w - взаимодействие 64 тел при задании угловой скорости вращения области;

6) galacf26.65b - взаимодействие 65 тел, в том числе 64 тела области и отдельного тела, расположенного на расстоянии от области.

Образцы входных файлов galacf26.dat имеются также в папке InCndFls\Preprtn для запуска на счет трех образцов файлов начальных условий: axsymb30.dat, axsymb3c.dat и axsymb3t.dat в графическом режиме счета.

10. Файлы начальных условий решенных задач

В папке InCndFls представлены файлы исходных данных и начальных условий решенных задач. Структура данных в этих файлах описана в п. 5. Здесь рассмотрим некоторые сведения об их использовании.

В этой папке находятся три вида файлов: 1) файлы начальных условий; 2) файлы изображений тел, которые имеют имена такие же, как и имена файлов первого вида; 3) входные файлы galacf26.dat, имена которых состоят из g26 и имен файлов первого вида. Файлы последнего вида приведены не для всех, а только для типичных задач.

Файлы sn49f.dat и sn49jplc.dat применялись при исследовании эволюции и Солнечной системы на интервале времени в 100 млн. лет [7-8]. В них представлены массы, координаты и скорости планет, Луны и Солнца на дату 30.12.1949 г. с юлианским днем JD = 2433280.5. Данные в файле sn49jplc.dat получены по более поздним эфемеридам Jet Propulsion Laboratory (JPL), чем файл sn49f.dat, а именно по DE406. Поэтому данные в этом файле точнее отражают свойства Солнечной системы. В этих файлах приводятся данные для 11 тел. Кроме того, зарезервирована строка для 12-ого тела. В эту строку можно занести данные космического аппарата или астероида, чтобы исследовать его движение в Солнечной системе.

В файле sn69JPA2.dat приведены начальные условия на дату 28 июня 1969 г. с JD = 2440400.5. В нем кроме вышеупомянутых тел имеются параметры трех астероидов: Церера, Паллада и Веста, т.е. всего 14 тел. Массы, координаты и скорости в этом файле получены из опорных данных эфемерид DE405 [11]. По сравнению с файлами sn49f.dat и sn49jplc.dat расчет по этому файлу дает более высокую точность совпадения результатов расчета с эфемеридами DE406 для планет от Юпитера до Плутона.

В файлах sn49e5n.dat, sn49e5n2.dat и sn49e5n3.dat представлены данные о составных моделях (1-ой, 2-ой и 3-ей, соответственно) вращения Земли [7, 12]. Данные об одиннадцати телах Солнечной системы предшествуют данным о телах модели вращения Земли. Их параметры основаны на тех же данных, что и параметры файла sn49f.dat, за исключением массы Земли. В этих файлах масса центрального тела составной Земли уменьшена на массу периферийных тел модели.

В файлах sn49S6J4.dat и sn49S6J5.dat представлены данные о 4-ой и 5-ой, соответственно, составных моделях Солнца [13-17]. Вначале идут данные об 11 телах Солнечной системы. Они основаны на тех же данных, что и параметры файла sn49jplc.dat, за исключением массы Солнца. Масса центрального тела модели Солнца уменьшена на массу периферийных тел его модели. Затем приводятся параметры периферийных тел модели. В файле sn49S6J5.dat модель состоит из 10 тел, а в файле sn49S6J4.dat - из 5 тел.

В файле sn49S6J6.dat представлены исходные данные и начальные условия задачи о воздействии составной модели Солнца только на Меркурий [17]. В этом файле представлены два объекта: Меркурий и составная модель Солнца, которая состоит из центрального тела в окружении 10 периферийных тел.

В файлах sn08AsAp.dat и sn08AsDA.dat [9, 10, 18-21] приведены данные об 11 телах Солнечной системы и астероидах: Апофис и 1950 DA, соответственно, на дату 30.11.08 г. с юлианским днем JD = 2454800.5. Массы 11 тел Солнечной системы такие же, как и в файле sn49jplc.dat. Их координаты и скорости основаны на эфемеридах DE406. Координаты и скорости астероидов рассчитывались по нашей методике [10] на основании их орбитальных элементов, взятых из базы малых тел NASA [22-23].

...