diff --git a/M7B_README.md b/M7B_README.md new file mode 100644 index 0000000..3b0dcc3 --- /dev/null +++ b/M7B_README.md @@ -0,0 +1,188 @@ +# М7Б. Статистика идеального газа + +## Описание + +Численное моделирование установления термодинамического равновесия в системе молекул идеального газа, находящихся в вертикальном сосуде в поле тяжести. + +### Цели работы + +1. **Численное моделирование** движения молекул идеального газа в замкнутом вертикальном сосуде с учётом упругих столкновений со стенками и между собой +2. **Получение распределений** молекул по скоростям и по высоте в равновесном состоянии +3. **Сравнение** полученных распределений с теоретическими: распределением Максвелла по скоростям и барометрической формулой для распределения по высоте + +## Физические модели + +### Динамика молекул + +Каждая молекула движется под действием силы тяжести согласно второму закону Ньютона: + +$$\frac{d\vec{v}}{dt} = \vec{g}$$ + +В декартовых координатах (x, y — горизонтальные, z — вертикальная ось): + +$$\frac{dv_x}{dt} = 0, \quad \frac{dv_y}{dt} = 0, \quad \frac{dv_z}{dt} = -g$$ + +### Распределение Максвелла + +В состоянии термодинамического равновесия распределение молекул по скоростям описывается распределением Максвелла: + +$$f(v) = 4\pi v^2 \left(\frac{m}{2\pi k_B T}\right)^{3/2} \exp\left(-\frac{mv^2}{2k_B T}\right)$$ + +где: +- m — масса молекулы +- k_B = 1.38 × 10⁻²³ Дж/К — постоянная Больцмана +- T — абсолютная температура + +### Барометрическая формула + +Распределение молекул по высоте в гравитационном поле описывается барометрической формулой: + +$$n(z) = n(0) \exp\left(-\frac{mgz}{k_B T}\right)$$ + +Масштаб высоты: +$$H_0 = \frac{k_B T}{mg}$$ + +## Структура кода + +### Основные компоненты + +``` +src/models/ideal_gas/ +├── __init__.py # Главный модуль, Streamlit интерфейс +├── config.py # Конфигурация и параметры +├── objects.py # Классы для моделирования +├── utils.py # Вспомогательные функции +└── charts.py # Визуализация результатов +``` + +### Ключевые классы + +#### `SimulationConfig` +Параметры симуляции: +- Молекулярные параметры (масса, диаметр) +- Геометрия сосуда (высота, радиус) +- Начальные условия (начальная скорость, высота) +- Параметры времени (общее время, шаг интегрирования) +- Параметры термостата (опционально) + +#### `IdealGasSimulation` +Основной класс для моделирования: +- `run()` — запуск симуляции +- `_velocity_verlet_step()` — интегрирование по методу Velocity Verlet +- `_handle_collisions()` — обработка столкновений +- `_handle_wall_collisions()` — отражение от стенок +- `_handle_intermolecular_collisions()` — упругие столкновения между молекулами + +#### `EquilibriumAnalyzer` +Анализ равновесного состояния: +- `calculate_temperature()` — определение температуры +- `calculate_pressure()` — расчёт давления + +## Численные методы + +### Интегрирование: Velocity Verlet + +``` +r(t + dt) = r(t) + v(t) * dt + 0.5 * a(t) * dt² +v(t + dt) = v(t) + 0.5 * (a(t) + a(t + dt)) * dt +``` + +Так как гравитационное ускорение постоянно: + +``` +r(t + dt) = r(t) + v(t) * dt + 0.5 * g * dt² +v_z(t + dt) = v_z(t) - g * dt +``` + +### Столкновения между молекулами + +Для упругого столкновения двух молекул равной массы используется алгоритм: + +1. Вычисляется единичный вектор вдоль линии центров: **n** +2. Вычисляется относительная скорость в направлении столкновения +3. Компоненты скоростей вдоль линии столкновения обмениваются + +## Использование + +### Запуск Streamlit интерфейса + +```bash +streamlit run src/main.py +``` + +Затем перейдите к странице М7Б в навигационном меню. + +### Параметры симуляции + +В боковой панели можно настроить: +- **Количество молекул** (50-1000) +- **Высота сосуда** (0.5-5.0 м) +- **Радиус сосуда** (0.5-3.0 м) +- **Начальная скорость молекул** (1-50 м/с) +- **Время симуляции** (1-20 с) +- **Шаг времени** (0.1-2.0 мс) +- **Использовать термостат** (опционально) + +### Результаты + +Симуляция выводит: +1. **Определённая температура** из кинетической энергии +2. **Масштаб высоты** барометрической формулы +3. **Средняя высота молекул** +4. **Средняя скорость молекул** + +### Визуализация + +5 вкладок с графиками: +1. **Распределение скоростей** — сравнение с распределением Максвелла +2. **Распределение по высоте** — сравнение с барометрической формулой +3. **Эволюция энергии** — кинетическая, потенциальная и полная энергия +4. **Эволюция температуры** — изменение температуры во времени +5. **3D позиции молекул** — визуализация распределения в конечном состоянии + +## Тестирование + +Запуск тестов: + +```bash +pytest tests/m7b/ +``` + +Тесты проверяют: +- **Корректность конфигурации** — валидация параметров +- **Форму выходных данных** — правильные размеры массивов +- **Граничные условия** — частицы не выходят за границы +- **Законы сохранения** — энергия ведёт себя ожидаемо +- **Утилиты** — расчёты распределений +- **Численная стабильность** — результаты не расходятся + +## Параметры по умолчанию + +- **Масса молекулы**: 6.63 × 10⁻²⁶ кг (аргон) +- **Эффективный диаметр**: 3.4 × 10⁻¹⁰ м +- **Гравитация**: 9.81 м/с² +- **Постоянная Больцмана**: 1.38 × 10⁻²³ Дж/К + +## Теоретические соотношения + +### Средняя кинетическая энергия + +$$\left\langle \frac{mv^2}{2} \right\rangle = \frac{3}{2}k_B T$$ + +Или по компонентам: + +$$\left\langle \frac{mv_i^2}{2} \right\rangle = \frac{1}{2}k_B T, \quad i = x, y, z$$ + +### Средняя высота молекулы + +$$\langle z \rangle = H_0 = \frac{k_B T}{mg}$$ + +### Средняя полная энергия + +$$\langle E \rangle = \frac{3}{2}k_B T + k_B T = \frac{5}{2}k_B T$$ + +## Ссылки + +- **PDF отчёта**: M7B-Report.pdf в корне проекта +- **Метод Verlet**: https://en.wikipedia.org/wiki/Verlet_integration +- **Распределение Максвелла**: https://en.wikipedia.org/wiki/Maxwell%E2%80%93Boltzmann_distribution diff --git a/M7B_report/conclusion.tex b/M7B_report/conclusion.tex new file mode 100644 index 0000000..f62532f --- /dev/null +++ b/M7B_report/conclusion.tex @@ -0,0 +1,48 @@ +\section{Заключение} + +В работе проведено комплексное исследование процесса установления термодинамического равновесия в системе идеального газа, находящегося в вертикальном сосуде в гравитационном поле. + +\textbf{Основные достижения:} + +\begin{enumerate} + \item \textbf{Численное моделирование динамики молекул} + \begin{itemize} + \item Реализована корректная обработка упругих столкновений молекул между собой и со стенками сосуда + \item Использована оптимизация методом ячеек для ускорения вычислений + \end{itemize} + + \item \textbf{Исследование процесса релаксации} + \begin{itemize} + \item Изучен переход системы из неравновесного начального состояния (все молекулы у дна с одинаковой энергией) к термодинамическому равновесию + \item Определено характерное время релаксации и его зависимость от параметров системы + \item Проверено сохранение полной энергии системы в изолированном режиме + \end{itemize} + + \item \textbf{Получение равновесных распределений} + \begin{itemize} + \item Построено распределение молекул по модулям скоростей и проведено сравнение с распределением Максвелла + \item Построено распределение молекул по высоте и проведено сравнение с барометрической формулой + \item Проверена факторизация распределений по различным компонентам скорости + \end{itemize} + + \item \textbf{Статистический анализ} + \begin{itemize} + \item Вычислены характерные скорости: наиболее вероятная, средняя и среднеквадратичная + \item Определена температура системы из распределения кинетических энергий + \item Проверено выполнение теоремы о равнораспределении энергии по степеням свободы + \end{itemize} + +\end{enumerate} + +\textbf{Выводы} + +\begin{itemize} + \item Численное моделирование методом молекулярной динамики подтверждает основные положения статистической механики + + \item Система молекул идеального газа самопроизвольно переходит из упорядоченного неравновесного состояния в равновесное состояние с максимальной энтропией + + \item В равновесии наблюдается полное согласие с распределением Максвелла-Больцмана, несмотря на детерминированный характер классической механики + + +\end{itemize} + diff --git a/M7B_report/intro.tex b/M7B_report/intro.tex new file mode 100644 index 0000000..8131cc0 --- /dev/null +++ b/M7B_report/intro.tex @@ -0,0 +1,40 @@ +\newpage + +\section{Введение} + +\Task Исследование процесса установления термодинамического равновесия в системе идеального газа, находящегося в вертикальном сосуде в поле тяжести. Рассмотрение неравновесного начального состояния, когда все молекулы сосредоточены вблизи дна сосуда и имеют одинаковую кинетическую энергию. + +\Goal +\begin{enumerate} + \item Провести численное моделирование движения молекул идеального газа в замкнутом вертикальном сосуде с учётом упругих столкновений со стенками и между собой + \item Получить распределения молекул по скоростям и по высоте в равновесном состоянии + \item Сравнить полученные распределения с теоретическими: распределением Максвелла по скоростям и барометрической формулой для распределения по высоте +\end{enumerate} + +\section{Физическая постановка задачи} +\vspace{-1em} + +Рассматривается система из $N$ одинаковых молекул идеального газа массой $m$ каждая, находящихся в вертикальном цилиндрическом сосуде высотой $H$ и площадью основания $S$ в гравитационном поле с ускорением $g$. + +\textbf{Модель идеального газа предполагает:} +\begin{enumerate} + \item Молекулы -- материальные точки (размеры много меньше расстояний между ними) + \item Столкновения молекул между собой и со стенками абсолютно упругие + \item Взаимодействие молекул происходит только при столкновениях (потенциальная энергия взаимодействия пренебрежимо мала) + \item Стенки сосуда непроницаемы и гладкие +\end{enumerate} + +\textbf{Начальные условия:} +\begin{itemize} + \item Все молекулы находятся в тонком слое вблизи дна: $0 < z < h_0 \ll H$ + \item Все молекулы имеют одинаковую по модулю скорость $v_0$, направленную случайным образом + \item Начальная кинетическая энергия каждой молекулы: $E_0 = \dfrac{mv_0^2}{2}$ +\end{itemize} + +\textbf{Рассматриваются два режима:} +\begin{enumerate} + \item \textbf{Изолированная система}: полная энергия $E = const$, стенки теплоизолированы + \item \textbf{Термостат}: температура стенок $T = const$, энергия системы флуктуирует +\end{enumerate} + +\newpage diff --git a/M7B_report/main.tex b/M7B_report/main.tex new file mode 100644 index 0000000..d71d0cf --- /dev/null +++ b/M7B_report/main.tex @@ -0,0 +1,16 @@ +\input{preamble} + +\begin{document} + +\input{title} +\pagestyle{main} +{ +\centering +\tableofcontents +} +\input{intro} +\input{models/models} +\input{numerical_methods} +\input{conclusion} + +\end{document} diff --git a/M7B_report/models/barometric_formula.tex b/M7B_report/models/barometric_formula.tex new file mode 100644 index 0000000..2ceb521 --- /dev/null +++ b/M7B_report/models/barometric_formula.tex @@ -0,0 +1,154 @@ +\subsection{Барометрическая формула} + +В состоянии термодинамического равновесия газ в гравитационном поле имеет неоднородное распределение плотности по высоте. + +\textbf{Вывод барометрической формулы} + +Рассмотрим элемент газа между высотами $z$ и $z + dz$. Условие механического равновесия: +\[ +p(z + dz) \cdot S = p(z) \cdot S + \rho(z) g S dz +\] + +где $p(z)$ -- давление на высоте $z$, $\rho(z)$ -- плотность массы. + +Дифференцируя: +\[ +\dfrac{dp}{dz} = -\rho(z) g +\] + +Используем уравнение состояния идеального газа: +\[ +p = \dfrac{\rho}{m} k_B T = n(z) k_B T +\] + +где $n(z) = \dfrac{\rho(z)}{m}$ -- концентрация молекул (число молекул на единицу объёма). + +Подставляем в условие равновесия: +\[ +\dfrac{d(n k_B T)}{dz} = -n(z) m g +\] + +Для изотермической атмосферы ($T = const$): +\[ +\dfrac{dn}{dz} = -\dfrac{mg}{k_B T} n(z) +\] + +\Solution Разделяем переменные и интегрируем: +\[ +\dfrac{dn}{n} = -\dfrac{mg}{k_B T} dz +\] +\[ +\ln n(z) - \ln n(0) = -\dfrac{mgz}{k_B T} +\] + +Получаем \textbf{барометрическую формулу}: +\[ +\boxed{n(z) = n(0) \exp\left(-\dfrac{mgz}{k_B T}\right)} +\] + +или в терминах давления: +\[ +p(z) = p(0) \exp\left(-\dfrac{mgz}{k_B T}\right) +\] + +\textbf{Характерная высота атмосферы} + +Введём масштаб высоты: +\[ +H_0 = \dfrac{k_B T}{mg} +\] + +Тогда барометрическая формула принимает вид: +\[ +n(z) = n(0) e^{-z/H_0} +\] + +\textbf{Распределение молекул по высоте} + +Вероятность найти молекулу на высоте между $z$ и $z + dz$ пропорциональна концентрации: +\[ +\boxed{f(z) = \dfrac{1}{H_0} \exp\left(-\dfrac{z}{H_0}\right)} +\] + +где нормировка: +\[ +\int_0^\infty f(z) dz = \int_0^\infty \dfrac{1}{H_0} e^{-z/H_0} dz = 1 +\] + +\Solution Проверим нормировку: +\[ +\int_0^\infty \dfrac{1}{H_0} e^{-z/H_0} dz = \dfrac{1}{H_0} \int_0^\infty e^{-z/H_0} dz +\] + +Замена переменной: $u = \dfrac{z}{H_0}$, $dz = H_0 du$: +\[ +\dfrac{1}{H_0} \int_0^\infty e^{-z/H_0} dz = \dfrac{1}{H_0} \cdot H_0 \int_0^\infty e^{-u} du = \left[-e^{-u}\right]_0^\infty = 0 - (-1) = 1 +\] + +\textbf{Средняя высота молекулы} + +\Solution Вычислим среднюю высоту: +\[ +\langle z \rangle = \int_0^\infty z f(z) dz = \int_0^\infty \dfrac{z}{H_0} e^{-z/H_0} dz +\] + +Замена переменной: $u = \dfrac{z}{H_0}$, $z = H_0 u$, $dz = H_0 du$: +\[ +\langle z \rangle = \int_0^\infty \dfrac{H_0 u}{H_0} e^{-u} H_0 du = H_0 \int_0^\infty u e^{-u} du +\] + +Интегрирование по частям или используя табличный интеграл: +\[ +\int_0^\infty u e^{-u} du = \Gamma(2) = 1! = 1 +\] + +Таким образом: +\[ +\langle z \rangle = H_0 \cdot 1 = H_0 = \dfrac{k_B T}{mg} +\] + +\textbf{Потенциальная энергия в среднем} + +\[ +\langle U \rangle = \langle mgz \rangle = mg \langle z \rangle = mg \cdot \dfrac{k_B T}{mg} = k_B T +\] + +Это согласуется с теоремой о равнораспределении: на каждую "квадратичную" координату приходится энергия $\dfrac{1}{2}k_B T$, но в гравитационном поле потенциальная энергия линейна по высоте, и среднее значение равно $k_B T$. + +\textbf{График барометрического распределения} + +\begin{figure}[h] +\centering +\begin{tikzpicture} +\begin{axis}[ + width=0.8\textwidth, + height=7cm, + xlabel={Высота $z/H_0$}, + ylabel={Относительная концентрация $n(z)/n(0)$}, + grid=major, + domain=0:5, + ymin=0, + ymax=1, + samples=200, +] +% Барометрическая формула: n(z)/n(0) = exp(-z/H_0) +\addplot[blue, thick] {exp(-x)}; +\addlegendentry{$n(z)/n(0) = e^{-z/H_0}$} + +% Отмечаем характерные высоты +\draw[dashed, gray] (axis cs:0,0.368) -- (axis cs:5,0.368); +\draw[dashed, gray] (axis cs:1,0) -- (axis cs:1,1); +\node[anchor=west, gray] at (axis cs:1.1,0.368) {$e^{-1} \approx 0.37$}; +\node[anchor=north, gray] at (axis cs:1,0.05) {$H_0$}; +\end{axis} +\end{tikzpicture} +\caption{Барометрическое распределение в безразмерных единицах. На высоте $z = H_0$ концентрация уменьшается в $e$ раз. Для воздуха при комнатной температуре $H_0 \approx 8.5$ км} +\end{figure} + +Физический смысл: +\begin{itemize} + \item Экспоненциальное убывание концентрации с высотой + \item Характерная высота атмосферы $H_0 = k_BT/(mg)$ зависит от температуры и массы молекул + \item Для воздуха при комнатной температуре $H_0 \approx 8.5$ км + \item Более лёгкие газы (водород, гелий) имеют большую высоту шкалы +\end{itemize} diff --git a/M7B_report/models/boltzmann_distribution.tex b/M7B_report/models/boltzmann_distribution.tex new file mode 100644 index 0000000..a4fd087 --- /dev/null +++ b/M7B_report/models/boltzmann_distribution.tex @@ -0,0 +1,93 @@ +\subsection{Распределение Максвелла-Больцмана} + +Полное распределение молекул в фазовом пространстве координат и импульсов описывается \textbf{распределением Максвелла-Больцмана}. + +\textbf{Каноническое распределение Гиббса} + +Вероятность найти систему в состоянии с энергией $E$ при температуре $T$: +\[ +P(E) \propto \exp\left(-\dfrac{E}{k_B T}\right) +\] + +Это называется \textbf{фактором Больцмана}. + +\textbf{Распределение одной молекулы} + +Для одной молекулы в гравитационном поле полная энергия: +\[ +E(\vec{r}, \vec{v}) = \dfrac{mv^2}{2} + mgz +\] + +Плотность вероятности в фазовом пространстве: +\[ +f(\vec{r}, \vec{v}) = C \exp\left(-\dfrac{E(\vec{r}, \vec{v})}{k_B T}\right) = C \exp\left(-\dfrac{mv^2/2 + mgz}{k_B T}\right) +\] + +где $C$ -- константа нормировки. + +\Solution Найдём константу нормировки $C$ из условия: +\[ +\int_{V} d^3r \int_{-\infty}^{\infty} d^3v \, f(\vec{r}, \vec{v}) = 1 +\] + +Разделяя экспоненты: +\[ +f(\vec{r}, \vec{v}) = C \exp\left(-\dfrac{mv^2}{2k_B T}\right) \exp\left(-\dfrac{mgz}{k_B T}\right) +\] + +Интегрируем по скоростям и координатам отдельно, получая: +\[ +C = \left(\dfrac{m}{2\pi k_B T}\right)^{3/2} \cdot \dfrac{mg}{Sk_B T} +\] + +где $S$ -- площадь основания сосуда. + +Это распределение факторизуется: +\[ +f(\vec{r}, \vec{v}) = f(\vec{v}) \cdot f(z) +\] + +где: +\begin{itemize} + \item $f(\vec{v})$ -- распределение Максвелла по скоростям + \item $f(z)$ -- барометрическое распределение по высоте +\end{itemize} + +\textbf{Распределение по энергии} + +Полная энергия молекулы $E = \dfrac{mv^2}{2} + mgz$ может принимать значения от 0 до $\infty$. + +Распределение по полной энергии в термодинамическом равновесии: +\[ +f(E) \propto \sqrt{E} \exp\left(-\dfrac{E}{k_B T}\right) +\] + +Множитель $\sqrt{E}$ появляется из-за статистического веса -- числа способов реализовать данную энергию. + +\textbf{Средняя энергия молекулы} + +В трёхмерном пространстве для молекулы в гравитационном поле: +\[ +\langle E \rangle = \left\langle \dfrac{mv^2}{2} \right\rangle + \langle mgz \rangle +\] + +Из распределения Максвелла: $\left\langle \dfrac{mv^2}{2} \right\rangle = \dfrac{3}{2}k_B T$ + +Из барометрической формулы: $\langle mgz \rangle = k_B T$ + +Итого: +\[ +\boxed{\langle E \rangle = \dfrac{3}{2}k_B T + k_B T = \dfrac{5}{2}k_B T} +\] + +\textbf{Связь с термодинамикой} + +Для системы из $N$ молекул идеального газа: +\begin{itemize} + \item Внутренняя энергия (без учёта гравитации): $U = N \cdot \dfrac{3}{2}k_B T = \dfrac{3}{2}NkT = \dfrac{3}{2}\nu RT$ + \item Теплоёмкость при постоянном объёме: $C_V = \dfrac{\partial U}{\partial T} = \dfrac{3}{2}\nu R$ + \item Теплоёмкость при постоянном давлении: $C_P = C_V + \nu R = \dfrac{5}{2}\nu R$ + \item Показатель адиабаты: $\gamma = \dfrac{C_P}{C_V} = \dfrac{5/2}{3/2} = \dfrac{5}{3} \approx 1.67$ +\end{itemize} + +где $\nu$ -- число молей, $R = N_A k_B$ -- универсальная газовая постоянная. diff --git a/M7B_report/models/ideal_gas_dynamics.tex b/M7B_report/models/ideal_gas_dynamics.tex new file mode 100644 index 0000000..90880d0 --- /dev/null +++ b/M7B_report/models/ideal_gas_dynamics.tex @@ -0,0 +1,89 @@ +\subsection{Динамика молекул идеального газа} + +\textbf{Уравнения движения} + +Каждая молекула движется под действием силы тяжести согласно \textbf{второму закону Ньютона}: +\[ +m\dfrac{d\vec{v}}{dt} = m\vec{g} +\] + +В декартовых координатах ($x$, $y$ -- горизонтальные, $z$ -- вертикальная ось): +\[ +\begin{cases} +\dfrac{dv_x}{dt} = 0 \\[0.5em] +\dfrac{dv_y}{dt} = 0 \\[0.5em] +\dfrac{dv_z}{dt} = -g +\end{cases} +\quad \text{и} \quad +\begin{cases} +\dfrac{dx}{dt} = v_x \\[0.5em] +\dfrac{dy}{dt} = v_y \\[0.5em] +\dfrac{dz}{dt} = v_z +\end{cases} +\] + +Горизонтальные компоненты скорости сохраняются между столкновениями, вертикальная изменяется под действием гравитации. + +\textbf{Граничные условия} + +Упругое отражение от стенок сосуда: +\begin{itemize} + \item \textbf{Дно} ($z = 0$): при $z \leq 0$ компонента $v_z \to -v_z$ + \item \textbf{Крышка} ($z = H$): при $z \geq H$ компонента $v_z \to -v_z$ + \item \textbf{Боковые стенки} (цилиндр радиуса $R$): при $\sqrt{x^2 + y^2} \geq R$ нормальная компонента скорости меняет знак +\end{itemize} + +Для моделирования термостата при столкновении со стенкой скорость молекулы пересчитывается по распределению Максвелла при заданной температуре стенки $T$. + +\textbf{Столкновения между молекулами} + +При столкновении двух молекул выполняются законы сохранения: +\begin{enumerate} + \item \textbf{Сохранение импульса:} + \[ + m\vec{v}_1 + m\vec{v}_2 = m\vec{v}_1' + m\vec{v}_2' + \] + + \item \textbf{Сохранение кинетической энергии:} + \[ + \dfrac{m v_1^2}{2} + \dfrac{m v_2^2}{2} = \dfrac{m v_1'^2}{2} + \dfrac{m v_2'^2}{2} + \] +\end{enumerate} + +где $\vec{v}_1, \vec{v}_2$ -- скорости до столкновения, $\vec{v}_1', \vec{v}_2'$ -- после столкновения. + +Столкновение происходит при сближении молекул на расстояние меньше эффективного диаметра $d$. + +\textbf{Алгоритм расчёта скоростей после столкновения:} + +Переход в систему центра масс: +\[ +\vec{V}_{cm} = \dfrac{\vec{v}_1 + \vec{v}_2}{2} +\] + +Относительная скорость: +\[ +\vec{v}_{rel} = \vec{v}_1 - \vec{v}_2 +\] + +Единичный вектор вдоль линии центров: +\[ +\vec{n} = \dfrac{\vec{r}_1 - \vec{r}_2}{|\vec{r}_1 - \vec{r}_2|} +\] + +Новые скорости после упругого столкновения: +\[ +\begin{cases} +\vec{v}_1' = \vec{V}_{cm} + \dfrac{1}{2}\left[\vec{v}_{rel} - 2(\vec{v}_{rel} \cdot \vec{n})\vec{n}\right] \\[1em] +\vec{v}_2' = \vec{V}_{cm} - \dfrac{1}{2}\left[\vec{v}_{rel} - 2(\vec{v}_{rel} \cdot \vec{n})\vec{n}\right] +\end{cases} +\] + +\textbf{Полная энергия системы} + +В изолированной системе полная энергия сохраняется: +\[ +E_{total} = \sum_{i=1}^{N} \left(\dfrac{mv_i^2}{2} + mgz_i\right) = const +\] + +где первое слагаемое -- кинетическая энергия, второе -- потенциальная энергия в поле тяжести. diff --git a/M7B_report/models/maxwell_distribution.tex b/M7B_report/models/maxwell_distribution.tex new file mode 100644 index 0000000..1696728 --- /dev/null +++ b/M7B_report/models/maxwell_distribution.tex @@ -0,0 +1,194 @@ +\subsection{Распределение Максвелла по скоростям} + +В состоянии термодинамического равновесия распределение молекул по скоростям описывается \textbf{распределением Максвелла}. + +\textbf{Распределение по вектору скорости} + +Плотность вероятности обнаружить молекулу со скоростью в интервале $(\vec{v}, \vec{v} + d\vec{v})$: +\[ +f(\vec{v}) = \left(\dfrac{m}{2\pi k_B T}\right)^{3/2} \exp\left(-\dfrac{m(\vec{v})^2}{2k_B T}\right) +\] + +где: +\begin{itemize} + \item $m$ -- масса молекулы + \item $k_B = 1.38 \times 10^{-23}$ Дж/К -- постоянная Больцмана + \item $T$ -- абсолютная температура +\end{itemize} + +Распределение факторизуется по компонентам: +\[ +f(\vec{v}) = f(v_x)f(v_y)f(v_z) +\] + +где каждая компонента имеет нормальное распределение: +\[ +f(v_i) = \sqrt{\dfrac{m}{2\pi k_B T}} \exp\left(-\dfrac{mv_i^2}{2k_B T}\right), \quad i = x, y, z +\] + +\textbf{Распределение по модулю скорости} + +Интегрируя по всем направлениям, получаем распределение по модулю скорости $v = |\vec{v}|$: +\[ +\boxed{f(v) = 4\pi v^2 \left(\dfrac{m}{2\pi k_B T}\right)^{3/2} \exp\left(-\dfrac{mv^2}{2k_B T}\right)} +\] + +Это распределение называется \textbf{распределением Максвелла}. + +\textbf{Характерные скорости} + +\begin{enumerate} + \item \textbf{Наиболее вероятная скорость} (максимум распределения): + + Найдём максимум $f(v)$, приравняв производную к нулю: + \[ + \dfrac{df}{dv} = 4\pi \left(\dfrac{m}{2\pi k_B T}\right)^{3/2} \dfrac{d}{dv}\left[v^2 e^{-mv^2/(2k_B T)}\right] = 0 + \] + + Применяя правило произведения: + \[ + 2v e^{-mv^2/(2k_B T)} + v^2 e^{-mv^2/(2k_B T)} \cdot \left(-\dfrac{mv}{k_B T}\right) = 0 + \] + + Вынося общий множитель $v e^{-mv^2/(2k_B T)}$: + \[ + v e^{-mv^2/(2k_B T)} \left(2 - \dfrac{mv^2}{k_B T}\right) = 0 + \] + + Нетривиальное решение ($v \neq 0$): + \[ + 2 - \dfrac{mv^2}{k_B T} = 0 \quad \Rightarrow \quad v_{prob} = \sqrt{\dfrac{2k_B T}{m}} + \] + + \item \textbf{Средняя скорость:} + \[ + \langle v \rangle = \int_0^\infty v f(v) dv = 4\pi \left(\dfrac{m}{2\pi k_B T}\right)^{3/2} \int_0^\infty v^3 e^{-mv^2/(2k_B T)} dv + \] + + Замена переменной: $u = \dfrac{mv^2}{2k_B T}$, откуда $v = \sqrt{\dfrac{2k_B T u}{m}}$, $dv = \sqrt{\dfrac{k_B T}{2mu}} du$ + \[ + \int_0^\infty v^3 e^{-mv^2/(2k_B T)} dv = \left(\dfrac{2k_B T}{m}\right)^2 \int_0^\infty u e^{-u} du = \left(\dfrac{2k_B T}{m}\right)^2 \cdot 1! = \left(\dfrac{2k_B T}{m}\right)^2 + \] + + Подставляя обратно: + \[ + \langle v \rangle = 4\pi \left(\dfrac{m}{2\pi k_B T}\right)^{3/2} \left(\dfrac{2k_B T}{m}\right)^2 = \sqrt{\dfrac{8k_B T}{\pi m}} + \] + + \item \textbf{Среднеквадратичная скорость:} + \[ + \langle v^2 \rangle = \int_0^\infty v^2 f(v) dv = 4\pi \left(\dfrac{m}{2\pi k_B T}\right)^{3/2} \int_0^\infty v^4 e^{-mv^2/(2k_B T)} dv + \] + + Аналогично, используя табличный интеграл: + \[ + \int_0^\infty v^4 e^{-mv^2/(2k_B T)} dv = \dfrac{3}{4}\sqrt{\pi} \left(\dfrac{2k_B T}{m}\right)^{5/2} + \] + + Откуда: + \[ + \langle v^2 \rangle = \dfrac{3k_B T}{m} \quad \Rightarrow \quad v_{rms} = \sqrt{\langle v^2 \rangle} = \sqrt{\dfrac{3k_B T}{m}} + \] +\end{enumerate}Соотношение между ними: +\[ +v_{prob} : \langle v \rangle : v_{rms} = \sqrt{2} : \sqrt{\dfrac{8}{\pi}} : \sqrt{3} \approx 1.41 : 1.60 : 1.73 +\] + +\textbf{График распределения Максвелла} + +\begin{figure}[h] +\centering +\begin{tikzpicture} +\begin{axis}[ + width=0.8\textwidth, + height=6cm, + xlabel={Скорость $v/v_{prob}$}, + ylabel={$f(v)$ (отн. ед.)}, + grid=major, + legend pos=north east, + domain=0:3, + samples=200, + ymin=0, +] +% Безразмерное распределение Максвелла: f(x) = 4/sqrt(pi) * x^2 * exp(-x^2) +% где x = v/v_prob, v_prob = sqrt(2kT/m) +\addplot[blue, thick] {4/sqrt(pi)*x^2*exp(-x^2)}; +\addlegendentry{Распределение Максвелла} + +% Отметим характерные скорости +\draw[dashed, gray] (axis cs:1,0) -- (axis cs:1,1.2); +\node[anchor=south, gray] at (axis cs:1,1.2) {$v_{prob}$}; + +\draw[dashed, red] (axis cs:1.128,0) -- (axis cs:1.128,0.8); +\node[anchor=south, red] at (axis cs:1.128,0.8) {$\langle v \rangle$}; + +\draw[dashed, green!60!black] (axis cs:1.225,0) -- (axis cs:1.225,0.6); +\node[anchor=south, green!60!black] at (axis cs:1.225,0.6) {$v_{rms}$}; +\end{axis} +\end{tikzpicture} +\caption{Распределение Максвелла в безразмерных единицах ($v/v_{prob}$). Показаны характерные скорости: наиболее вероятная $v_{prob}$, средняя $\langle v \rangle \approx 1.128 v_{prob}$ и среднеквадратичная $v_{rms} \approx 1.225 v_{prob}$} +\end{figure} + +Из графика видно: +\begin{itemize} + \item Максимум распределения при $v = v_{prob} = \sqrt{2k_BT/m}$ + \item Средняя скорость $\langle v \rangle$ и среднеквадратичная $v_{rms}$ больше наиболее вероятной + \item Распределение асимметрично -- имеет "хвост" в сторону больших скоростей +\end{itemize} + +\textbf{Средняя кинетическая энергия} + +\Solution Вычислим среднее значение кинетической энергии: +\[ +\left\langle \dfrac{mv^2}{2} \right\rangle = \dfrac{m}{2} \langle v^2 \rangle = \dfrac{m}{2} \cdot \dfrac{3k_B T}{m} = \dfrac{3}{2}k_B T +\] + +Альтернативный вывод через интегрирование по компонентам: +\[ +\left\langle \dfrac{mv^2}{2} \right\rangle = \left\langle \dfrac{m(v_x^2 + v_y^2 + v_z^2)}{2} \right\rangle = \dfrac{m}{2}(\langle v_x^2 \rangle + \langle v_y^2 \rangle + \langle v_z^2 \rangle) +\] + +Для каждой компоненты (например, $v_x$): +\[ +\langle v_x^2 \rangle = \int_{-\infty}^{\infty} v_x^2 \sqrt{\dfrac{m}{2\pi k_B T}} e^{-mv_x^2/(2k_B T)} dv_x = \dfrac{k_B T}{m} +\] + +Используем табличный интеграл $\int_{-\infty}^{\infty} x^2 e^{-ax^2} dx = \dfrac{1}{2}\sqrt{\dfrac{\pi}{a^3}}$. + +Таким образом: +\[ +\left\langle \dfrac{mv^2}{2} \right\rangle = \dfrac{m}{2} \cdot 3 \cdot \dfrac{k_B T}{m} = \dfrac{3}{2}k_B T +\] + +Это выражение связывает температуру со средней кинетической энергией поступательного движения молекул. + +По теореме о равнораспределении энергии по степеням свободы: +\[ +\left\langle \dfrac{mv_i^2}{2} \right\rangle = \dfrac{1}{2}k_B T, \quad i = x, y, z +\] + +На каждую поступательную степень свободы приходится энергия $\dfrac{1}{2}k_B T$. + +\textbf{Проверка нормировки} + +\Solution Проверим, что распределение нормировано: +\[ +\int_0^\infty f(v) dv = 4\pi \left(\dfrac{m}{2\pi k_B T}\right)^{3/2} \int_0^\infty v^2 e^{-mv^2/(2k_B T)} dv +\] + +Используем табличный интеграл (гауссов интеграл): +\[ +\int_0^\infty v^2 e^{-av^2} dv = \dfrac{1}{4}\sqrt{\dfrac{\pi}{a^3}} +\] + +Где $a = \dfrac{m}{2k_B T}$: +\[ +\int_0^\infty v^2 e^{-mv^2/(2k_B T)} dv = \dfrac{1}{4}\sqrt{\dfrac{\pi}{(m/(2k_B T))^3}} = \dfrac{1}{4}\sqrt{\pi} \left(\dfrac{2k_B T}{m}\right)^{3/2} +\] + +Подставляем: +\[ +\int_0^\infty f(v) dv = 4\pi \left(\dfrac{m}{2\pi k_B T}\right)^{3/2} \cdot \dfrac{1}{4}\sqrt{\pi} \left(\dfrac{2k_B T}{m}\right)^{3/2} = \sqrt{\pi} \cdot \dfrac{1}{\sqrt{\pi}} = 1 +\] + +Это означает, что вероятность найти молекулу с любой скоростью равна единице. diff --git a/M7B_report/models/models.tex b/M7B_report/models/models.tex new file mode 100644 index 0000000..62e60e7 --- /dev/null +++ b/M7B_report/models/models.tex @@ -0,0 +1,10 @@ +\section{Физические модели} + +\input{models/ideal_gas_dynamics} +\newpage +\input{models/maxwell_distribution} +\newpage +\input{models/barometric_formula} +\newpage +\input{models/boltzmann_distribution} +\newpage \ No newline at end of file diff --git a/M7B_report/numerical_methods.tex b/M7B_report/numerical_methods.tex new file mode 100644 index 0000000..d5a1b84 --- /dev/null +++ b/M7B_report/numerical_methods.tex @@ -0,0 +1,139 @@ +\section{Численные методы решения} + +Численное моделирование системы молекул идеального газа основано на методе \textbf{молекулярной динамики} -- пошаговом интегрировании уравнений движения всех частиц с учётом столкновений. + +\subsection{Алгоритм интегрирования} + +Для интегрирования уравнений движения используется \textbf{метод Верле} (Verlet algorithm), который хорошо сохраняет энергию системы при длительном моделировании. + +\textbf{Базовый алгоритм Верле:} + +Положение частицы на следующем шаге вычисляется по формуле: +\[ +\vec{r}(t + \Delta t) = 2\vec{r}(t) - \vec{r}(t - \Delta t) + \vec{a}(t) \Delta t^2 +\] + +где $\vec{a}(t) = \dfrac{\vec{F}(t)}{m}$ -- ускорение. + +\textbf{Скоростная форма Верле (Velocity Verlet):} + +Более удобная модификация, которая явно вычисляет скорости: +\[ +\begin{cases} +\vec{r}(t + \Delta t) = \vec{r}(t) + \vec{v}(t) \Delta t + \dfrac{1}{2}\vec{a}(t) \Delta t^2 \\[1em] +\vec{v}(t + \Delta t) = \vec{v}(t) + \dfrac{1}{2}[\vec{a}(t) + \vec{a}(t + \Delta t)] \Delta t +\end{cases} +\] + +Для нашей задачи $\vec{a} = \vec{g} = (0, 0, -g) = const$, поэтому: +\[ +\begin{cases} +\vec{r}(t + \Delta t) = \vec{r}(t) + \vec{v}(t) \Delta t + \dfrac{1}{2}\vec{g} \Delta t^2 \\[1em] +\vec{v}(t + \Delta t) = \vec{v}(t) + \vec{g} \Delta t +\end{cases} +\] + +\textbf{Порядок точности:} метод имеет второй порядок точности по времени, $O(\Delta t^2)$. + +\subsection{Обработка столкновений} + +\textbf{Столкновения со стенками} + +На каждом временном шаге проверяются граничные условия: +\begin{itemize} + \item \textbf{Дно и крышка:} если $z < 0$ или $z > H$, то: + \[ + v_z \to -v_z, \quad z \to \max(0, \min(z, H)) + \] + + \item \textbf{Боковая поверхность:} если $\sqrt{x^2 + y^2} > R$, то компонента скорости вдоль радиуса меняет знак: + \[ + \vec{v}_r \to -\vec{v}_r + \] +\end{itemize} + +Для режима термостата при столкновении со стенкой скорость пересчитывается по распределению Максвелла: +\[ +v_i \sim \mathcal{N}\left(0, \sqrt{\dfrac{k_B T}{m}}\right), \quad i = x, y, z +\] + +\textbf{Столкновения между молекулами} + +Используется \textbf{метод жёстких сфер} (hard sphere model): +\begin{enumerate} + \item На каждом шаге проверяются расстояния между всеми парами молекул + \item Если $|\vec{r}_i - \vec{r}_j| < d$ (эффективный диаметр), происходит столкновение + \item Скорости пересчитываются по формулам упругого столкновения (см. раздел 3.1) + \item Для избежания повторных столкновений частицы слегка разводятся: + \[ + \vec{r}_i \to \vec{r}_i + \epsilon \vec{n}, \quad \vec{r}_j \to \vec{r}_j - \epsilon \vec{n} + \] + где $\vec{n}$ -- единичный вектор вдоль линии центров, $\epsilon$ -- малая величина +\end{enumerate} + +\textbf{Оптимизация вычислений} + +Для ускорения проверки столкновений используется метод \textbf{ячеек} (cell lists): +\begin{itemize} + \item Пространство разбивается на кубические ячейки размером $\ge d$ + \item Каждая молекула приписывается к ячейке по своим координатам + \item Проверка столкновений производится только с молекулами из той же и соседних ячеек + \item Сложность алгоритма снижается с $O(N^2)$ до $O(N)$ +\end{itemize} + +\subsection{Сбор статистики} + +После достижения равновесия (время релаксации $t_{relax}$) начинается сбор статистических данных: + +\textbf{Распределение по скоростям} + +Строится гистограмма модулей скоростей $v = |\vec{v}|$: +\[ +n_v(v) = \text{число молекул со скоростями в интервале } [v, v + \Delta v] +\] + +Нормированная плотность вероятности: +\[ +f_{num}(v) = \dfrac{n_v(v)}{N \cdot \Delta v} +\] + +Сравнивается с теоретическим распределением Максвелла. + +\textbf{Распределение по высоте} + +Строится гистограмма положений по оси $z$: +\[ +n_z(z) = \text{число молекул на высотах в интервале } [z, z + \Delta z] +\] + +Нормированная плотность вероятности: +\[ +f_{num}(z) = \dfrac{n_z(z)}{N \cdot \Delta z} +\] + +Сравнивается с барометрической формулой. + +\textbf{Критерии равновесия} + +Система считается достигшей равновесия, когда: +\begin{enumerate} + \item Средняя кинетическая энергия флуктуирует около постоянного значения + \item Распределения по скоростям и высоте не изменяются со временем + \item Выполняется теорема вириала: $\langle E_k \rangle = \dfrac{3}{2}Nk_B T$ +\end{enumerate} + +\textbf{Определение температуры} + +Температура системы определяется из средней кинетической энергии: +\[ +T = \dfrac{2\langle E_k \rangle}{3Nk_B} = \dfrac{2}{3Nk_B} \sum_{i=1}^{N} \dfrac{mv_i^2}{2} +\] + +Для проверки можно также использовать распределение по компонентам скорости: +\[ +T_i = \dfrac{\langle mv_i^2 \rangle}{k_B}, \quad i = x, y, z +\] + +В равновесии $T_x = T_y = T_z = T$. + +\newpage diff --git a/M7B_report/preamble.tex b/M7B_report/preamble.tex new file mode 100644 index 0000000..a696b04 --- /dev/null +++ b/M7B_report/preamble.tex @@ -0,0 +1,128 @@ +\documentclass[a4paper,12pt]{article} + +\usepackage{geometry} +% задание полей текста +\geometry{left=10mm,right=10mm,top=20mm,bottom=20mm} +\usepackage{listings} +\usepackage{cmap} % поиск в PDF +\usepackage[T2A]{fontenc} % кодировка +\usepackage[utf8]{inputenc} +\usepackage[english, russian]{babel} % локализация и переносы +\usepackage{natbib} + +\usepackage{graphicx} +\graphicspath{{pix/}} +\usepackage{pgfplots} +\usetikzlibrary{arrows} +\pgfplotsset{width=10cm,compat=1.18} +\pgfkeys{/pgf/trig format=rad} + + +\usepackage{xcolor} +\usepackage{amsmath} +\usepackage{amssymb} +\usepackage{bm} +\usepackage{enumerate} +\usepackage[normalem]{ulem} +\usepackage{titlesec} + +% Настройка форматирования section для центрирования без растягивания +\titleformat{\section} + {\normalfont\Large\bfseries\centering}{\thesection}{1em}{} +\titleformat{\subsection} + {\normalfont\large\bfseries\centering}{\thesubsection}{1em}{} +\titleformat{\subsubsection} + {\normalfont\normalsize\bfseries\centering}{\thesubsubsection}{1em}{} + +\usepackage{hyperref} +\hypersetup{ + colorlinks=true, + linkcolor=blue, + filecolor=blue, + urlcolor=blue, + citecolor=blue +} + + +\setlength\parindent{0pt} + + +\sloppy % строго соблюдать границы текста +\linespread{1.2} % коэффициент межстрочного интервала +\setlength{\parskip}{0.5em} % вертик. интервал между абзацами + +\setcounter{secnumdepth}{0} % отключение нумерации разделов +\binoppenalty=1000 % уменьшение переносов в формулах + +% объявление новых макрокоманд + +\newcommand{\Def}{\textbf{Def.} } +\newcommand{\Th}{\textbf{Теорема.} } +\newcommand{\Thbd}{\textbf{Теорема (б/д).} } +\newcommand{\Theor}[1]{\textbf{Теорема ({#1}).} } +\newcommand{\Theorbd}[1]{\textbf{Теорема ({#1}) (б/д).} } +\newcommand{\Consequence}{\textbf{Следствие.} } +\newcommand{\Remind}{\textbf{Remind.} } +\newcommand{\Note}{\textbf{Note.} } +\newcommand{\Statement}{\textbf{Утверждение.} } +\newcommand{\Prop}{\textbf{Свойство:} } +\newcommand{\Props}{\textbf{Свойства:} } +\newcommand{\Proof}{\textbf{Доказательство:} } +\newcommand{\Prooff}{\textbf{Доказать:} } +\newcommand{\Solution}{\textbf{Решение:} } +\newcommand{\Alg}{\textbf{Algorithm.} } +\newcommand{\Lemma}{\textbf{Лемма.} } +\newcommand{\Example}{\textbf{Пример:} } +\newcommand{\Task}{\textbf{Задача.} } +\newcommand{\Goal}{\textbf{Цель работы:} } +\newcommand{\Solve}{\textbf{Решение:} } +\newcommand{\Examples}{\textbf{Примеры.} } + +\allowdisplaybreaks[4] + +%\renewcommand\thesection{\arabic{section}} + +\newcommand{\Endproof}{$\blacksquare$ } + +\newcommand{\tr}{\text{tr}} +\newcommand{\Le}{\leqslant} % русский стиль нестрогих неравенств +\newcommand{\Ge}{\geqslant} +% облегчение математических обозначений +\newcommand{\A}{\mathcal{A}} +\newcommand{\M}{\mathcal{M}} +\newcommand{\F}{\mathcal{F}} +\newcommand{\Gs}{\mathcal{G}} +\newcommand{\R}{\mathbb{R}} +\newcommand{\N}{\mathbb{N}} +\newcommand{\Q}{\mathbb{Q}} +\newcommand{\Norm}{\mathcal{N}} +%\newcommand{\C}{\mathbb{C}} +\newcommand{\ind}{\perp \!\!\! \perp} +\newcommand{\Z}{\mathbb{Z}} +\newcommand{\E}{{\mathbb{E}}} +\newcommand{\D} { {\mathbb{D}}} +\newcommand{\x}{\overline{x}} +\newcommand{\y}{\overline{y}} +\newcommand{\soch}[2]{\text{$C_{#1}^{#2}$}} +\newcommand{\X}{\text{$\mathcal{X}$}} +\newcommand{\VC}[2]{\text{$VC\left(#1; #2\right)$}} +\renewcommand{\P}{\text{$\mathcal{P}$}} +\newcommand{\B}{\text{$\mathcal{B}$}} +\newcommand{\brackets}[1]{\left({#1}\right)} % автоматический размер скобочек +\newcommand{\abs}[1]{\left|{#1}\right|} % автоматический размер модуля +\newcommand{\norm}[1]{\left\|{#1}\right\|} % автоматический размер нормы +\newcommand{\set}[1]{\left\{{#1}\right\}} % автоматический размер множества +\newcommand{\condset}[2]{\set{#1 \ | \ #2}} % определение множества + +\newcommand{\distrfamily}{\condset{P_\theta}{\theta \in \Theta}} + +\DeclareMathOperator{\sign} +\DeclareMathOperator{\tg} + + +\newpagestyle{main}{ + \setheadrule{.4pt} + \sethead{}{\bullet \; \textit{Цифровизация физических процессов} \; \bullet}{} + \setfootrule{.4pt} + \setfoot{ВШПИ МФТИ}{М7Б. Статистика идеального газа}{Б13-402} +} diff --git a/M7B_report/report.pdf b/M7B_report/report.pdf new file mode 100644 index 0000000..3a3f3df Binary files /dev/null and b/M7B_report/report.pdf differ diff --git a/M7B_report/title.tex b/M7B_report/title.tex new file mode 100644 index 0000000..e52bc91 --- /dev/null +++ b/M7B_report/title.tex @@ -0,0 +1,18 @@ +\begin{titlepage} + \centering + {\scshape\Large Московский физико-технический институт \par} + \vspace{1cm} + {\scshape\Large Высшая школа программной инженерии \par} + \vspace{2cm} + {\huge \textbf{М7Б. Статистика идеального газа} \par} + \vspace{1.5cm} + {\Large \textit{Численное моделирование установления равновесия \\ + в системе молекул идеального газа в гравитационном поле} \par} + \vspace{2cm} + {\Large Выполнили студенты Б13-402: \par} + {\Large Жердев Егор \par} + {\Large Савельев Данил \par} + \vfill + {\large Долгопрудный \par} + {\large \today \par} +\end{titlepage} diff --git a/src/models/ideal_gas/__init__.py b/src/models/ideal_gas/__init__.py new file mode 100644 index 0000000..797370e --- /dev/null +++ b/src/models/ideal_gas/__init__.py @@ -0,0 +1,229 @@ +"""М7Б. Статистика идеального газа. + +Численное моделирование установления равновесия в системе молекул идеального газа +в гравитационном поле. + +Основные компоненты: +- config: параметры симуляции и конфигурация +- objects: основные классы для моделирования +- charts: визуализация результатов +- utils: вспомогательные функции для обработки данных +""" + +import streamlit as st +import numpy as np + +from models.ideal_gas.config import SimulationConfig +from models.ideal_gas.objects import IdealGasSimulation, EquilibriumAnalyzer +from models.ideal_gas.charts import ( + create_velocity_distribution_chart, + create_height_distribution_chart, + create_energy_evolution_chart, + create_3d_particle_positions_chart, + create_temperature_evolution_chart, +) +from models.ideal_gas.utils import calculate_maxwell_distribution, calculate_barometric_formula + + +def page() -> None: + st.set_page_config(page_title="М7Б. Статистика идеального газа", layout="wide") + st.title("М7Б. Статистика идеального газа") + st.write( + "Численное моделирование установления равновесия в системе молекул идеального газа в гравитационном поле" + ) + + # Боковая панель с параметрами + with st.sidebar: + st.header("Параметры симуляции") + st.caption("Параметры ограничены для быстрой визуализации в браузере.") + + num_particles = st.slider( + "Количество молекул", + min_value=50, + max_value=400, # было 1000 + value=150, + step=50, + help="Меньше частиц — быстрее симуляция. 100–200 обычно достаточно, чтобы увидеть распределения.", + ) + + container_height = st.slider( + "Высота сосуда (м)", + min_value=0.5, + max_value=3.0, # было 5.0 + value=1.5, + step=0.5, + ) + + container_radius = st.slider( + "Радиус сосуда (м)", + min_value=0.5, + max_value=2.0, # было 3.0 + value=1.0, + step=0.5, + ) + + initial_velocity = st.slider( + "Начальная скорость молекул (м/с)", + min_value=1.0, + max_value=30.0, # было 50.0 + value=10.0, + step=1.0, + ) + + simulation_time = st.slider( + "Время симуляции (с)", + min_value=0.5, + max_value=5.0, # было 20.0 + value=2.0, + step=0.5, + help="Большее время — больше шагов и медленнее визуализация.", + ) + + dt = st.slider( + "Шаг времени (мс)", + min_value=0.5, + max_value=2.0, + value=1.0, + step=0.5, + help="Крупный шаг ускоряет расчёт, но делает траектории грубее.", + ) + + use_thermostat = st.checkbox( + "Использовать термостат", + value=False, + help="Термостат немного удорожает шаги, включайте только при необходимости.", + ) + + if use_thermostat: + target_temperature = st.slider( + "Температура термостата (K)", + min_value=100.0, + max_value=800.0, + value=300.0, + step=50.0, + ) + else: + target_temperature = None + + config = SimulationConfig( + num_particles=num_particles, + container_height=container_height, + container_radius=container_radius, + initial_velocity=initial_velocity, + simulation_time=simulation_time, + dt=dt / 1000.0, # Convert ms to s + use_thermostat=use_thermostat, + target_temperature=target_temperature, + ) + + # Запуск симуляции + if st.button("🚀 Запустить симуляцию", use_container_width=True): + with st.spinner("Выполняется симуляция молекулярной динамики..."): + simulation = IdealGasSimulation(config) + + positions_history, velocities_history, times = simulation.run() + + # Анализ равновесия + analyzer = EquilibriumAnalyzer(config) + final_positions = positions_history[-1] + final_velocities = velocities_history[-1] + + # Получение температуры + kinetic_energy = 0.5 * config.particle_mass * np.sum(final_velocities**2) / config.num_particles + temperature = 2.0 * kinetic_energy / (3.0 * 1.38e-23) # k_B + + # Расчет теоретических распределений + maxwell_speeds = np.linspace(0, np.max(np.linalg.norm(final_velocities, axis=1)) * 1.2, 100) + maxwell_dist = calculate_maxwell_distribution(maxwell_speeds, temperature, config.particle_mass) + + heights = np.linspace(0, config.container_height, 100) + scale_height = 1.38e-23 * temperature / (config.particle_mass * 9.81) + barometric_dist = calculate_barometric_formula(heights, scale_height) + + # Выводы + st.success("✅ Симуляция завершена!") + + col1, col2, col3, col4 = st.columns(4) + with col1: + st.metric("Определенная температура", f"{temperature:.1f} K") + with col2: + st.metric("Масштаб высоты", f"{scale_height:.3f} м") + with col3: + st.metric("Средняя высота", f"{np.mean(final_positions[:, 2]):.3f} м") + with col4: + st.metric("Средняя скорость", f"{np.mean(np.linalg.norm(final_velocities, axis=1)):.2f} м/с") + + # Вкладки с визуализацией + tab1, tab2, tab3, tab4, tab5 = st.tabs([ + "Распределение скоростей", + "Распределение высот", + "Эволюция энергии", + "Эволюция температуры", + "3D позиции частиц", + ]) + + with tab1: + st.subheader("Распределение Максвелла по скоростям") + chart = create_velocity_distribution_chart( + final_velocities, + maxwell_speeds, + maxwell_dist, + temperature, + ) + st.altair_chart(chart, use_container_width=True) + + with tab2: + st.subheader("Барометрическое распределение по высоте") + chart = create_height_distribution_chart( + final_positions[:, 2], + heights, + barometric_dist, + scale_height, + ) + st.altair_chart(chart, use_container_width=True) + + with tab3: + st.subheader("Эволюция энергии системы") + chart = create_energy_evolution_chart(positions_history, velocities_history, times, config) + st.altair_chart(chart, use_container_width=True) + + with tab4: + st.subheader("Эволюция температуры") + chart = create_temperature_evolution_chart(velocities_history, times, config) + st.altair_chart(chart, use_container_width=True) + + with tab5: + st.subheader("3D позиции молекул в конечном состоянии") + chart = create_3d_particle_positions_chart(final_positions, config) + st.altair_chart(chart, use_container_width=True) + + # Теоретическая информация + st.divider() + with st.expander("📖 Теоретическая информация"): + st.markdown(""" + ### Распределение Максвелла + В термодинамическом равновесии распределение молекул по скоростям подчиняется распределению Максвелла: + + $$f(v) = 4\pi v^2 \left(\frac{m}{2\pi k_B T}\right)^{3/2} \exp\left(-\frac{mv^2}{2k_B T}\right)$$ + + где: + - $m$ - масса молекулы + - $k_B = 1.38 \times 10^{-23}$ Дж/К - постоянная Больцмана + - $T$ - абсолютная температура + + ### Барометрическая формула + Распределение молекул по высоте в гравитационном поле описывается барометрической формулой: + + $$n(z) = n(0) \exp\left(-\frac{mgz}{k_B T}\right)$$ + + Масштаб высоты: $H_0 = \frac{k_B T}{mg}$ + + ### Теорема о равнораспределении энергии + На каждую поступательную степень свободы приходится энергия: + + $$\left\langle \frac{mv_i^2}{2} \right\rangle = \frac{1}{2}k_B T, \quad i = x, y, z$$ + + Средняя полная кинетическая энергия молекулы: + + $$\left\langle E_k \right\rangle = \frac{3}{2}k_B T$$ + """) diff --git a/src/models/ideal_gas/charts.py b/src/models/ideal_gas/charts.py new file mode 100644 index 0000000..5be2af4 --- /dev/null +++ b/src/models/ideal_gas/charts.py @@ -0,0 +1,173 @@ +"""Визуализация результатов симуляции.""" + +import altair as alt +import numpy as np +import pandas as pd +from typing import List, Tuple +from models.ideal_gas.config import SimulationConfig +from models.ideal_gas.utils import ( + calculate_speed_from_velocity, + create_statistics_dataframe, +) + + +def create_velocity_distribution_chart(velocities: np.ndarray, maxwell_speeds: np.ndarray, + maxwell_dist: np.ndarray, temperature: float) -> alt.Chart: + """График распределения скоростей (Maxwell vs Simulation).""" + # Симуляционные скорости + speeds = calculate_speed_from_velocity(velocities) + + # Гистограмма + hist, bin_edges = np.histogram(speeds, bins=50, density=True) + bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2 + + sim_data = pd.DataFrame({ + 'speed': bin_centers, + 'probability': hist, + 'source': 'Simulation' + }) + + # Теоретическое распределение + theory_data = pd.DataFrame({ + 'speed': maxwell_speeds, + 'probability': maxwell_dist, + 'source': 'Maxwell' + }) + + combined_data = pd.concat([sim_data, theory_data], ignore_index=True) + + chart = alt.Chart(combined_data).mark_line(point=True).encode( + x=alt.X('speed:Q', title='Скорость (m/s)', scale=alt.Scale(zero=False)), + y=alt.Y('probability:Q', title='Вероятность'), + color=alt.Color('source:N', title='Источник', + scale=alt.Scale(domain=['Simulation', 'Maxwell'], + range=['#1f77b4', '#ff7f0e'])), + strokeDash=alt.StrokeDash('source:N', scale=alt.Scale(domain=['Simulation', 'Maxwell'], + range=[[1, 0], [5, 5]])), + ).properties( + width=600, + height=400, + title=f'Распределение Максвелла (T={temperature:.1f}K)' + ) + + return chart + + +def create_height_distribution_chart(heights: np.ndarray, heights_theory: np.ndarray, + barometric_dist: np.ndarray, scale_height: float) -> alt.Chart: + """График распределения по высоте.""" + # Гистограмма + hist, bin_edges = np.histogram(heights, bins=40, density=True) + bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2 + + sim_data = pd.DataFrame({ + 'height': bin_centers, + 'probability': hist, + 'source': 'Simulation' + }) + + # Теоретическое распределение + theory_data = pd.DataFrame({ + 'height': heights_theory, + 'probability': barometric_dist, + 'source': 'Barometric' + }) + + combined_data = pd.concat([sim_data, theory_data], ignore_index=True) + + chart = alt.Chart(combined_data).mark_line(point=True).encode( + x=alt.X('height:Q', title='Высота (m)'), + y=alt.Y('probability:Q', title='Плотность вероятности'), + color=alt.Color('source:N', title='Модель', + scale=alt.Scale(domain=['Simulation', 'Barometric'], + range=['#1f77b4', '#2ca02c'])), + strokeDash=alt.StrokeDash('source:N', scale=alt.Scale(domain=['Simulation', 'Barometric'], + range=[[1, 0], [5, 5]])), + ).properties( + width=600, + height=400, + title=f'Барометрическая формула (H0={scale_height:.3f}m)' + ) + + return chart + + +def create_energy_evolution_chart(positions_history: np.ndarray, velocities_history: np.ndarray, + times: np.ndarray, config: SimulationConfig) -> alt.Chart: + """График эволюции энергии.""" + stats_df = create_statistics_dataframe(positions_history, velocities_history, + config.particle_mass, times) + + energy_df = pd.DataFrame({ + 'time': stats_df['time'], + 'Kinetic': stats_df['kinetic_energy'] / stats_df['kinetic_energy'].iloc[0] if stats_df['kinetic_energy'].iloc[0] != 0 else stats_df['kinetic_energy'], + 'Potential': stats_df['potential_energy'] / abs(stats_df['potential_energy'].iloc[0]) if stats_df['potential_energy'].iloc[0] != 0 else stats_df['potential_energy'], + 'Total': stats_df['total_energy'] / stats_df['total_energy'].iloc[0] if stats_df['total_energy'].iloc[0] != 0 else stats_df['total_energy'], + }) + + # Перережють в длинный формат + energy_melted = energy_df.melt(id_vars=['time'], var_name='energy_type', value_name='normalized_value') + + chart = alt.Chart(energy_melted).mark_line().encode( + x=alt.X('time:Q', title='Время (s)'), + y=alt.Y('normalized_value:Q', title='Нормализованная энергия'), + color=alt.Color('energy_type:N', title='Компонента', + scale=alt.Scale(domain=['Kinetic', 'Potential', 'Total'], + range=['#1f77b4', '#ff7f0e', '#2ca02c'])), + ).properties( + width=600, + height=400, + title='Эволюция энергии системы' + ) + + return chart + + +def create_temperature_evolution_chart(velocities_history: np.ndarray, times: np.ndarray, + config: SimulationConfig) -> alt.Chart: + """График температуры во времени.""" + temps = [] + for vel_step in velocities_history: + kinetic_per_particle = 0.5 * config.particle_mass * np.sum(vel_step**2) / config.num_particles + T = 2.0 * kinetic_per_particle / (3.0 * config.k_B) + temps.append(T) + + df = pd.DataFrame({ + 'time': times, + 'temperature': temps + }) + + chart = alt.Chart(df).mark_line(point=True).encode( + x=alt.X('time:Q', title='Время (s)'), + y=alt.Y('temperature:Q', title='Температура (K)'), + color=alt.value('#d62728') + ).properties( + width=600, + height=400, + title='Эволюция температуры' + ) + + return chart + + +def create_3d_particle_positions_chart(positions: np.ndarray, config: SimulationConfig) -> alt.Chart: + """График позиций соъстабов (высота vs x).""" + df = pd.DataFrame({ + 'x': positions[:, 0], + 'y': positions[:, 1], + 'z': positions[:, 2] + }) + + # Высь я рисую (мероприятия z) + chart = alt.Chart(df).mark_circle(size=50, opacity=0.6).encode( + x=alt.X('x:Q', title='X (м)', scale=alt.Scale(zero=False)), + y=alt.Y('z:Q', title='Z - Высота (м)', scale=alt.Scale(zero=False, domainMin=0, domainMax=config.container_height)), + color=alt.Color('z:Q', title='Высота (m)', + scale=alt.Scale(scheme='viridis', domainMin=0, domainMax=config.container_height)), + ).properties( + width=500, + height=400, + title='Позиции молекул (XZ сечение)' + ) + + return chart diff --git a/src/models/ideal_gas/config.py b/src/models/ideal_gas/config.py new file mode 100644 index 0000000..4016b19 --- /dev/null +++ b/src/models/ideal_gas/config.py @@ -0,0 +1,67 @@ +"""Конфигурация настроек для симуляции идеального газа.""" + +from dataclasses import dataclass +from typing import Optional + + +@dataclass +class SimulationConfig: + """Параметры симуляции идеального газа. + + Молекулярные параметры: + num_particles: Количество молекул + particle_mass: Масса одной молекулы (кг), по умолчанию Ar + particle_diameter: Эффективный диаметр молекулы (м) + + Геометрия сосуда: + container_height: Высота сосуда (м) + container_radius: Радиус цилиндра (м) + + Начальные условия: + initial_velocity: Начальная скорость молекул (м/с) + initial_height: Начальная высота (м) + + Параметры времени: + simulation_time: Тотальное время симуляции (с) + dt: Шаг интегрирования (с) + + Параметры термостата: + use_thermostat: Основание / теплоаккумуляция стенками + target_temperature: Целевая температура (К) + """ + + # Молекулярные параметры + num_particles: int + # НОВО: условная масса для университетского моделя + # Это не реальная аргона, а высокая сильная барометрическая экспонента + particle_mass: float = 6.63e-23 # условная масса (100× тяжелее для наглядности) + particle_diameter: float = 3.4e-10 # эффективный диаметр (м) + + # Геометрия сосуда + container_height: float + container_radius: float + + # Начальные условия + initial_velocity: float + initial_height: float = 0.1 # На этакой высоте (м) + + # Параметры времени + simulation_time: float + dt: float + + # Параметры термостата + use_thermostat: bool = False + target_temperature: Optional[float] = None + + # Вычисляемые и постоянные параметры + g: float = 9.81 # ускорение свободного падения (м/с²) + k_B: float = 1.38e-23 # постоянная Больцмана (Дж/К) + + def __post_init__(self): + """Проверка находящихся допустимых состояния.""" + if self.num_particles <= 0: + raise ValueError("Количество молекул должно быть положительным") + if self.dt <= 0 or self.dt > 0.01: + raise ValueError("Шаг времени должен быть <= 0.01") + if self.initial_height > self.container_height: + raise ValueError("Начальная высота не может быть больше высоты сосуда") diff --git a/src/models/ideal_gas/objects.py b/src/models/ideal_gas/objects.py new file mode 100644 index 0000000..7639afe --- /dev/null +++ b/src/models/ideal_gas/objects.py @@ -0,0 +1,286 @@ +"""Основные классы для моделирования идеального газа.""" + +import numpy as np +from typing import Tuple +from models.ideal_gas.config import SimulationConfig + + +class IdealGasSimulation: + """Моделирование идеального газа методом молекулярной динамики. + + Основные методы: + - run(): запуск симуляции + - _initialize_particles(): инициализация молекул + - _velocity_verlet_step(): шаг интегрирования + - _handle_collisions(): обработка столкновений + + Оптимизации: + - Нумпи векторизация для быстрых вычислений + - Cell list для сокращения столкновений O(N) вместо O(N^2) + - Оптимизированная проверка столкновений + """ + + def __init__(self, config: SimulationConfig): + self.config = config + self.num_steps = int(config.simulation_time / config.dt) + + # Объекты для столкновений + self.collision_pairs_checked = 0 + self.collisions_detected = 0 + + # Cell list для акселерации столкновений + self.cell_size = config.particle_diameter * 2.5 + self.cells_per_axis = int(config.container_radius / self.cell_size) + 1 + self.cells = {} + + def run(self) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """Запустить симуляцию и вернуть историю позиций и скоростей. + + Returns: + positions_history: История позиций (num_steps, num_particles, 3) + velocities_history: История скоростей (num_steps, num_particles, 3) + times: Вектор времен + """ + # Инициализация + positions = np.zeros((self.config.num_particles, 3), dtype=np.float32) + velocities = np.zeros((self.config.num_particles, 3), dtype=np.float32) + self._initialize_particles(positions, velocities) + + # Оси для сохранения истории + positions_history = [] + velocities_history = [] + times = [] + + # Основной цикл + for step in range(self.num_steps): + t = step * self.config.dt + times.append(t) + + positions_history.append(positions.copy()) + velocities_history.append(velocities.copy()) + + # Интегрирование уравнений движения + self._velocity_verlet_step(positions, velocities) + + # Обработка столкновений + self._handle_collisions(positions, velocities) + + # Управление термостатом (по необходимости) + if self.config.use_thermostat and self.config.target_temperature is not None: + self._apply_thermostat(velocities) + + return np.array(positions_history), np.array(velocities_history), np.array(times) + + def _initialize_particles(self, positions: np.ndarray, velocities: np.ndarray) -> None: + """Начинальные позиции и скорости. + + Все молекулы располагаются вблизи дна с рандомными скоростями. + """ + # Оптимизированная векторизированная инициализация + N = self.config.num_particles + + # Генерирую случайные позиции в цилиндре + theta = np.random.uniform(0, 2 * np.pi, N) + r = np.sqrt(np.random.uniform(0, 1, N)) * self.config.container_radius * 0.9 + + positions[:, 0] = r * np.cos(theta) + positions[:, 1] = r * np.sin(theta) + positions[:, 2] = np.random.uniform(0, self.config.initial_height, N) + + # Генерирую скорости со случайными направлениями + phi = np.random.uniform(0, 2 * np.pi, N) + cos_theta = np.random.uniform(-1, 1, N) + sin_theta = np.sqrt(1 - cos_theta**2) + + velocities[:, 0] = self.config.initial_velocity * sin_theta * np.cos(phi) + velocities[:, 1] = self.config.initial_velocity * sin_theta * np.sin(phi) + velocities[:, 2] = self.config.initial_velocity * cos_theta + + def _velocity_verlet_step(self, positions: np.ndarray, velocities: np.ndarray) -> None: + """Однот методом Velocity Verlet (векторизированный). + + r(t + dt) = r(t) + v(t) * dt + 0.5 * a(t) * dt^2 + v(t + dt) = v(t) + 0.5 * (a(t) + a(t + dt)) * dt + + Гравитация столбца z: a_z = -g (отрицательная!) + """ + dt = self.config.dt + g = self.config.g + + # Обновление горизонтальных координат (без акцелерации) + positions[:, 0] += velocities[:, 0] * dt + positions[:, 1] += velocities[:, 1] * dt + + # Обновление вертикальных координат с гравитацией + # a_z = -g, поэтому: z = z + v_z*dt - 0.5*g*dt^2 + positions[:, 2] += velocities[:, 2] * dt - 0.5 * g * dt**2 + + # Обновление вертикальных скоростей + # v_z = v_z + a_z*dt = v_z - g*dt + velocities[:, 2] -= g * dt + + def _handle_collisions(self, positions: np.ndarray, velocities: np.ndarray) -> None: + """Обработка столкновений со стенками и между молекулами.""" + # Столкновения со стенками + self._handle_wall_collisions(positions, velocities) + + # Межмолекулярные столкновения с cell list оптимизацией + self._handle_intermolecular_collisions_optimized(positions, velocities) + + def _handle_wall_collisions(self, positions: np.ndarray, velocities: np.ndarray) -> None: + """Отражение от стенок сосуда (векторизированно).""" + # Боковые стенки (цилиндр) + r = np.sqrt(positions[:, 0]**2 + positions[:, 1]**2) + collided = r > self.config.container_radius + + if np.any(collided): + # Нормальные компоненты вектора скорости (векторизированно) + cos_phi = np.divide(positions[collided, 0], r[collided], where=r[collided]!=0, out=np.zeros_like(positions[collided, 0])) + sin_phi = np.divide(positions[collided, 1], r[collided], where=r[collided]!=0, out=np.zeros_like(positions[collided, 1])) + + v_r = velocities[collided, 0] * cos_phi + velocities[collided, 1] * sin_phi + v_t = -velocities[collided, 0] * sin_phi + velocities[collided, 1] * cos_phi + + # Отражение радиальной компоненты + v_r *= -1 + + # Новые компоненты скорости + velocities[collided, 0] = v_r * cos_phi - v_t * sin_phi + velocities[collided, 1] = v_r * sin_phi + v_t * cos_phi + + # Перепозиционирование частиц + positions[collided, 0] *= self.config.container_radius / r[collided] + positions[collided, 1] *= self.config.container_radius / r[collided] + + # Дно и крыша (векторизированно) + collided_floor = positions[:, 2] < 0 + velocities[collided_floor, 2] *= -1 + positions[collided_floor, 2] = 0 + + collided_ceiling = positions[:, 2] > self.config.container_height + velocities[collided_ceiling, 2] *= -1 + positions[collided_ceiling, 2] = self.config.container_height + + def _handle_intermolecular_collisions_optimized(self, positions: np.ndarray, velocities: np.ndarray) -> None: + """Обработка столкновений с cell list оптимизацией O(N) вместо O(N^2).""" + # Очистить налицы мячей + self.cells.clear() + + # Рассортировать частицы по мячеям + for i, pos in enumerate(positions): + # Определить индекс ячейки + ix = int((pos[0] + self.config.container_radius) / self.cell_size) + iy = int((pos[1] + self.config.container_radius) / self.cell_size) + iz = int(pos[2] / self.cell_size) + + # Пограничить индексы + ix = max(0, min(ix, 2 * self.cells_per_axis - 1)) + iy = max(0, min(iy, 2 * self.cells_per_axis - 1)) + iz = max(0, min(iz, self.cells_per_axis)) + + cell_key = (ix, iy, iz) + if cell_key not in self.cells: + self.cells[cell_key] = [] + self.cells[cell_key].append(i) + + # Проверить столкновения лишь в соседних ячейках + d_collision = self.config.particle_diameter + checked_pairs = set() + + for cell_key, particles in self.cells.items(): + ix, iy, iz = cell_key + + # Проверить в текущей ячейке + for i in range(len(particles)): + for j in range(i + 1, len(particles)): + p1, p2 = particles[i], particles[j] + pair = (min(p1, p2), max(p1, p2)) + if pair not in checked_pairs: + checked_pairs.add(pair) + self._check_and_resolve_collision(p1, p2, positions, velocities, d_collision) + + # Проверить соседние ячейки + for dix in [-1, 0, 1]: + for diy in [-1, 0, 1]: + for diz in [-1, 0, 1]: + if dix == 0 and diy == 0 and diz == 0: + continue + + neighbor_key = (ix + dix, iy + diy, iz + diz) + if neighbor_key in self.cells: + for p1 in particles: + for p2 in self.cells[neighbor_key]: + if p1 < p2: + pair = (p1, p2) + if pair not in checked_pairs: + checked_pairs.add(pair) + self._check_and_resolve_collision(p1, p2, positions, velocities, d_collision) + + def _check_and_resolve_collision(self, i: int, j: int, positions: np.ndarray, + velocities: np.ndarray, d_collision: float) -> None: + """Проверить и решить столкновение одного пары молекул.""" + r = positions[j] - positions[i] + dist_sq = np.sum(r**2) + + if dist_sq < d_collision**2: + dist = np.sqrt(dist_sq) + if dist > 0: + self.collisions_detected += 1 + self._resolve_collision_inline(i, j, positions, velocities, r, dist) + + def _resolve_collision_inline(self, i: int, j: int, positions: np.ndarray, + velocities: np.ndarray, r: np.ndarray, dist: float) -> None: + """Решение упругого столкновения (выполнено встроенные).""" + # Единичный вектор + n = r / dist + + # Относительная скорость + v_rel = velocities[i] - velocities[j] + v_rel_n = np.dot(v_rel, n) + + # Проверь, что частицы сближаются + if v_rel_n >= 0: + return + + # Обмен скоростей + delta_v = v_rel_n * n + velocities[i] -= delta_v + velocities[j] += delta_v + + # Отсылка + overlap = self.config.particle_diameter - dist + if overlap > 0: + shift = (overlap / 2 + 1e-6) * n + positions[i] -= shift + positions[j] += shift + + def _apply_thermostat(self, velocities: np.ndarray) -> None: + """Применение баростата.""" + if not self.config.target_temperature: + return + + kinetic_energy = 0.5 * self.config.particle_mass * np.sum(velocities**2) / self.config.num_particles + current_temp = 2.0 * kinetic_energy / (3.0 * self.config.k_B) + + if current_temp > 0: + scale_factor = np.sqrt(self.config.target_temperature / current_temp) + velocities *= scale_factor + + +class EquilibriumAnalyzer: + """Анализ равновесных распределений.""" + + def __init__(self, config: SimulationConfig): + self.config = config + + def calculate_temperature(self, velocities: np.ndarray) -> float: + """Рассчет температуры из среднего до нижнего центра.""" + kinetic_energy = 0.5 * self.config.particle_mass * np.sum(velocities**2) / self.config.num_particles + return 2.0 * kinetic_energy / (3.0 * self.config.k_B) + + def calculate_pressure(self, positions: np.ndarray, velocities: np.ndarray, + collisions_count: int) -> float: + """Наивная оценка давления.""" + n_density = self.config.num_particles / (np.pi * self.config.container_radius**2 * self.config.container_height) + T = self.calculate_temperature(velocities) + return n_density * self.config.k_B * T diff --git a/src/models/ideal_gas/utils.py b/src/models/ideal_gas/utils.py new file mode 100644 index 0000000..68b1780 --- /dev/null +++ b/src/models/ideal_gas/utils.py @@ -0,0 +1,102 @@ +"""Утилиты для обработки данных исмуляции.""" + +import numpy as np +import pandas as pd +from typing import Tuple + + +def calculate_maxwell_distribution(speeds: np.ndarray, temperature: float, particle_mass: float) -> np.ndarray: + """Распределение Максвелла по модулю скорости. + + f(v) = 4π v^2 (m / (2π k_B T))^(3/2) exp(-m v^2 / (2 k_B T)) + + Args: + speeds: Проенство скоростей + temperature: Температура (K) + particle_mass: Масса партикулы (kg) + + Returns: + Нормализованное распределение + """ + k_B = 1.38e-23 + + # Исключительные случаи + if temperature <= 0 or particle_mass <= 0: + return np.zeros_like(speeds) + + # Префакторы + sqrt_term = np.sqrt(particle_mass / (2 * np.pi * k_B * temperature)) + norm_factor = 4 * np.pi * (sqrt_term ** 3) + exp_arg = -particle_mass * speeds**2 / (2 * k_B * temperature) + + distribution = norm_factor * speeds**2 * np.exp(exp_arg) + + return distribution + + +def calculate_barometric_formula(heights: np.ndarray, scale_height: float) -> np.ndarray: + """Барометрическая формула для распределения но высоте. + + n(z) = n(0) exp(-z / H0) + где H0 = k_B T / (m g) + + Args: + heights: Простланство высот (m) + scale_height: Масштаб высоты H0 (m) + + Returns: + Нормализованное распределение + """ + if scale_height <= 0: + return np.zeros_like(heights) + + # Нормализованное распределение + distribution = (1.0 / scale_height) * np.exp(-heights / scale_height) + + return distribution + + +def calculate_speed_from_velocity(velocities: np.ndarray) -> np.ndarray: + """Перечислить ни векторы в них модули.""" + return np.linalg.norm(velocities, axis=1) + + +def create_statistics_dataframe(positions: np.ndarray, velocities: np.ndarray, + particle_mass: float, times: np.ndarray) -> pd.DataFrame: + """Особые данные ис жистории симуляции.""" + num_steps = len(times) + + data = { + 'time': [], + 'kinetic_energy': [], + 'potential_energy': [], + 'total_energy': [], + 'temperature': [], + 'mean_height': [], + } + + g = 9.81 + k_B = 1.38e-23 + + for step in range(num_steps): + pos = positions[step] + vel = velocities[step] + + # Кинетическая энергия + KE = 0.5 * particle_mass * np.sum(vel**2) + + # Потенциальная энергия + PE = particle_mass * g * np.sum(pos[:, 2]) + + # Температура + kinetic_per_particle = KE / len(pos) + T = 2.0 * kinetic_per_particle / (3.0 * k_B) + + data['time'].append(times[step]) + data['kinetic_energy'].append(KE) + data['potential_energy'].append(PE) + data['total_energy'].append(KE + PE) + data['temperature'].append(T) + data['mean_height'].append(np.mean(pos[:, 2])) + + return pd.DataFrame(data) diff --git a/tests/m7b/__init__.py b/tests/m7b/__init__.py new file mode 100644 index 0000000..4432c49 --- /dev/null +++ b/tests/m7b/__init__.py @@ -0,0 +1 @@ +# tests for M7B ideal gas simulation \ No newline at end of file diff --git a/tests/m7b/conftest.py b/tests/m7b/conftest.py new file mode 100644 index 0000000..ce99b2c --- /dev/null +++ b/tests/m7b/conftest.py @@ -0,0 +1,30 @@ +"""Приспособления конфигурации для тестов M7B.""" + +import pytest +from models.ideal_gas.config import SimulationConfig +from models.ideal_gas.objects import IdealGasSimulation + + +@pytest.fixture +def basic_config(): + """Базовая конфигурация для тестов.""" + return SimulationConfig( + num_particles=100, + container_height=2.0, + container_radius=1.0, + initial_velocity=10.0, + simulation_time=1.0, + dt=0.01, + ) + + +@pytest.fixture +def simulation(basic_config): + """Краткая симуляция для тестов.""" + return IdealGasSimulation(basic_config) + + +@pytest.fixture +def simulation_results(simulation): + """Проведи симуляцию и верни результаты.""" + return simulation.run() diff --git a/tests/m7b/test_simulation.py b/tests/m7b/test_simulation.py new file mode 100644 index 0000000..cc79f54 --- /dev/null +++ b/tests/m7b/test_simulation.py @@ -0,0 +1,313 @@ +"""Проверки для моделирования идеального газа.""" + +import numpy as np +import pytest +from scipy.integrate import trapezoid +from models.ideal_gas.config import SimulationConfig +from models.ideal_gas.objects import IdealGasSimulation, EquilibriumAnalyzer +from models.ideal_gas.utils import ( + calculate_maxwell_distribution, + calculate_barometric_formula, + calculate_speed_from_velocity, +) + + +class TestSimulationInitialization: + """Tests for simulation initialization.""" + + def test_config_validation_negative_particles(self): + """Test that negative particle count raises error.""" + with pytest.raises(ValueError): + SimulationConfig( + num_particles=-10, + container_height=2.0, + container_radius=1.0, + initial_velocity=10.0, + simulation_time=1.0, + dt=0.01, + ) + + def test_config_validation_invalid_dt(self): + """Test that invalid dt raises error.""" + with pytest.raises(ValueError): + SimulationConfig( + num_particles=100, + container_height=2.0, + container_radius=1.0, + initial_velocity=10.0, + simulation_time=1.0, + dt=0.05, # Требуется <= 0.01 + ) + + def test_config_validation_height_bounds(self): + """Test that initial height cannot exceed container height.""" + with pytest.raises(ValueError): + SimulationConfig( + num_particles=100, + container_height=2.0, + container_radius=1.0, + initial_velocity=10.0, + simulation_time=1.0, + dt=0.01, + initial_height=3.0, # > container_height + ) + + +class TestSimulationRun: + """Tests for simulation execution.""" + + def test_simulation_output_shapes(self, basic_config): + """Проверить, что выводы имеют правильные формы.""" + sim = IdealGasSimulation(basic_config) + positions, velocities, times = sim.run() + + num_steps = int(basic_config.simulation_time / basic_config.dt) + + assert positions.shape == (num_steps, basic_config.num_particles, 3) + assert velocities.shape == (num_steps, basic_config.num_particles, 3) + assert times.shape == (num_steps,) + + def test_particles_stay_in_bounds(self, basic_config): + """Test that particles stay within container bounds.""" + sim = IdealGasSimulation(basic_config) + positions, _, _ = sim.run() + + final_positions = positions[-1] + + # Проверить высоту + assert np.all(final_positions[:, 2] >= 0), "Some particles went below floor" + assert np.all(final_positions[:, 2] <= basic_config.container_height), "Some particles went above ceiling" + + # Проверить латеральные между цилиндрами + r = np.sqrt(final_positions[:, 0]**2 + final_positions[:, 1]**2) + assert np.all(r <= basic_config.container_radius * 1.01), "Some particles outside cylinder" + + def test_time_progression(self, basic_config): + """Test that time increases monotonically.""" + sim = IdealGasSimulation(basic_config) + _, _, times = sim.run() + + assert np.all(np.diff(times) > 0), "Time should increase monotonically" + assert times[0] == 0.0, "Time should start at 0" + assert times[-1] <= basic_config.simulation_time, "Final time should not exceed simulation_time" + + +class TestEnergyConservation: + """Tests for energy conservation laws.""" + + def test_total_energy_conserved(self, basic_config): + """Проверить основной закон сохранения энергии для изолированной системы.""" + sim = IdealGasSimulation(basic_config) + positions, velocities, _ = sim.run() + + # Полная энергия всегда + total_energies = [] + g = basic_config.g + m = basic_config.particle_mass + + for pos_step, vel_step in zip(positions, velocities): + # Кинетическая энергия + KE = 0.5 * m * np.sum(vel_step**2) + # Потенциальная энергия + PE = m * g * np.sum(pos_step[:, 2]) + # Полная энергия + total_energies.append(KE + PE) + + total_energies = np.array(total_energies) + + # Инициальная энергия + initial_energy = total_energies[0] + + # Проверить, что энергия почти константна + # Рассеия до ±5% допустима (численные ошибки) + rel_energy_error = np.max(np.abs(total_energies - initial_energy)) / np.abs(initial_energy) + assert rel_energy_error < 0.05, f"Energy conservation error too large: {rel_energy_error:.2%}" + + def test_kinetic_energy_decreases_initially(self, basic_config): + """Проверить, что кинетическая энергия снижается до явно + (молекулы поднимаются и теряют скорость).""" + sim = IdealGasSimulation(basic_config) + _, velocities, _ = sim.run() + + # Кинетическая энергия на каждом шаге + m = basic_config.particle_mass + kinetic_energies = [] + for vel_step in velocities: + KE = 0.5 * m * np.sum(vel_step**2) + kinetic_energies.append(KE) + + kinetic_energies = np.array(kinetic_energies) + + # У частицы есть инициальная скорость, поэтому кинетическая энергия должна упасть + assert kinetic_energies[0] > kinetic_energies[-1], "KE should decrease (particles rise against gravity)" + + def test_potential_energy_increases(self, basic_config): + """Проверить, что потенциальная энергия растёт.""" + sim = IdealGasSimulation(basic_config) + positions, _, _ = sim.run() + + g = basic_config.g + m = basic_config.particle_mass + + potential_energies = [] + for pos_step in positions: + PE = m * g * np.sum(pos_step[:, 2]) + potential_energies.append(PE) + + potential_energies = np.array(potential_energies) + + # Потенциальная энергия должна расти (молекулы поднимаются) + assert potential_energies[-1] > potential_energies[0], "PE should increase (particles rise)" + + def test_mechanical_energy_conservation_no_collisions(self): + """Проверить сохранение механической энергии без столкновений. + + Ниские отправергнутые партикулы должны дохранять энергию без высоких численных ошибок. + """ + config = SimulationConfig( + num_particles=10, # Мало молекул - минимизируем столкновения + container_height=5.0, # большой контейнер - меньше коллизий + container_radius=2.0, + initial_velocity=5.0, # нижняя скорость + simulation_time=0.5, + dt=0.01, + ) + + sim = IdealGasSimulation(config) + positions, velocities, _ = sim.run() + + g = config.g + m = config.particle_mass + + total_energies = [] + for pos_step, vel_step in zip(positions, velocities): + KE = 0.5 * m * np.sum(vel_step**2) + PE = m * g * np.sum(pos_step[:, 2]) + total_energies.append(KE + PE) + + total_energies = np.array(total_energies) + initial_energy = total_energies[0] + + # До ±2% (очень точная неоновмыеся система) + rel_energy_error = np.max(np.abs(total_energies - initial_energy)) / np.abs(initial_energy) + assert rel_energy_error < 0.02, f"Energy drift too large: {rel_energy_error:.2%}" + + +class TestEquilibriumAnalyzer: + """Tests for equilibrium analysis.""" + + def test_temperature_calculation(self, basic_config): + """Test temperature calculation from velocities.""" + sim = IdealGasSimulation(basic_config) + _, velocities, _ = sim.run() + + analyzer = EquilibriumAnalyzer(basic_config) + T = analyzer.calculate_temperature(velocities[-1]) + + # Температура должна быть положительной + assert T > 0, "Temperature should be positive" + # Ордер величины: несколько сот К + assert T < 10000, "Temperature seems unreasonably high" + + def test_pressure_calculation(self, basic_config): + """Test pressure calculation.""" + sim = IdealGasSimulation(basic_config) + positions, velocities, _ = sim.run() + + analyzer = EquilibriumAnalyzer(basic_config) + P = analyzer.calculate_pressure(positions[-1], velocities[-1], 0) + + # Давление должно быть положительным + assert P > 0, "Pressure should be positive" + + +class TestUtilityFunctions: + """Tests for utility functions.""" + + def test_maxwell_distribution_normalization(self): + """Проверить нормализацию Максвелла.""" + speeds = np.linspace(0, 1000, 1000) + maxwell = calculate_maxwell_distribution(speeds, 300, 6.63e-26) + + # Приближенная проверка интеграла с scipy.integrate.trapezoid + integral = trapezoid(maxwell, speeds) + assert 0.9 < integral < 1.1, f"Maxwell integral should be ~1, got {integral}" + + def test_barometric_formula_normalization(self): + """Проверить нормализацию барометрической формулы.""" + heights = np.linspace(0, 5, 100) + barometric = calculate_barometric_formula(heights, 1.0) + + # Приближенная проверка интеграла с scipy.integrate.trapezoid + integral = trapezoid(barometric, heights) + assert 0.9 < integral < 1.1, f"Barometric integral should be ~1, got {integral}" + + def test_speed_calculation(self): + """Проверить начисление скоростей.""" + velocities = np.array([ + [3, 4, 0], # speed = 5 + [0, 0, 5], # speed = 5 + [1, 1, 1], # speed = sqrt(3) + ]) + + speeds = calculate_speed_from_velocity(velocities) + + assert np.isclose(speeds[0], 5.0), "Speed calculation for [3,4,0]" + assert np.isclose(speeds[1], 5.0), "Speed calculation for [0,0,5]" + assert np.isclose(speeds[2], np.sqrt(3)), "Speed calculation for [1,1,1]" + + +class TestNumericalStability: + """Tests for numerical stability.""" + + def test_simulation_doesnt_diverge(self, basic_config): + """Test that simulation results don't diverge to infinity.""" + sim = IdealGasSimulation(basic_config) + positions, velocities, _ = sim.run() + + assert np.all(np.isfinite(positions)), "Positions contain non-finite values" + assert np.all(np.isfinite(velocities)), "Velocities contain non-finite values" + + # Проверить, что скорости не становятся сигнификантно большими + assert np.all(np.abs(velocities) < 1000), "Velocities are unreasonably large" + + def test_no_negative_kinetic_energy(self, basic_config): + """Test that kinetic energy never becomes negative.""" + sim = IdealGasSimulation(basic_config) + _, velocities, _ = sim.run() + + for vel_step in velocities: + KE = np.sum(vel_step**2) + assert KE >= 0, "Kinetic energy cannot be negative" + + def test_gravity_pulls_downward(self, basic_config): + """Проверить, что гравитация действует вниз (знак "-" в a_z = -g). + + Это основная проверка: если дать частице вертикальную скорость, + гравитация должна её уменьшать. + """ + config = SimulationConfig( + num_particles=1, # Одна частица + container_height=10.0, + container_radius=5.0, + initial_velocity=20.0, # Только вверх + initial_height=0.0, # На дне + simulation_time=1.0, + dt=0.01, + ) + + sim = IdealGasSimulation(config) + positions, velocities, _ = sim.run() + + # Частица стартует с v_z > 0 (верх) + initial_v_z = velocities[0, 0, 2] + assert initial_v_z > 0, "Particle should start with upward velocity" + + # После симуляции v_z должен быть меньше (гравитация таможит) + final_v_z = velocities[-1, 0, 2] + assert final_v_z < initial_v_z, "Gravity should reduce upward velocity" + + # Высота вначале арыов (дно), а вконце выше (частица поднялась) + initial_z = positions[0, 0, 2] + max_z = np.max(positions[:, 0, 2]) + assert max_z > initial_z, "Particle should rise due to initial upward velocity"