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/main.py b/src/main.py index 947693b..d0cfb45 100644 --- a/src/main.py +++ b/src/main.py @@ -1,6 +1,7 @@ import streamlit from models.billiard import page as billiard_page +from models.ideal_gas import page as ideal_gas_page from models.pendulum import page as pendulum_page from models.rolling_ball import page as rolling_ball_page from models.stone_flight import page as stone_flight_page @@ -18,6 +19,11 @@ def main() -> None: streamlit.Page( pendulum_page, title="M5 Физический маятник", url_path="pendulum" ), + streamlit.Page( + ideal_gas_page, + title="M7Б Статистика идеального газа со столкновениями частиц", + url_path="ideal_gas", + ), ]) pages.run() diff --git a/src/models/ideal_gas/__init__.py b/src/models/ideal_gas/__init__.py new file mode 100644 index 0000000..930e586 --- /dev/null +++ b/src/models/ideal_gas/__init__.py @@ -0,0 +1,204 @@ +import numpy as np +import streamlit as st + +from models.ideal_gas.charts import ( + create_3d_particle_positions_chart, + create_energy_evolution_chart, + create_height_distribution_chart, + create_temperature_evolution_chart, + create_velocity_distribution_chart, +) +from models.ideal_gas.config import SimulationConfig +from models.ideal_gas.objects import IdealGasSimulation +from models.ideal_gas.utils import ( + calculate_barometric_formula, + calculate_maxwell_distribution, +) + + +def page() -> None: # noqa: PLR0915, PLR0914 + st.set_page_config( + page_title="М7Б. Статистика идеального газа со столкновениями частиц", + layout="wide", + ) + st.title("М7Б. Статистика идеального газа со столкновениями частиц") + + with st.sidebar: + st.header("Параметры симуляции") + + num_particles = st.slider( + "Количество молекул", + min_value=50, + max_value=400, + value=150, + step=50, + help=( + "Меньше частиц - быстрее симуляция. " + "100–200 обычно достаточно, чтобы увидеть распределения." + ), + ) + + container_height = st.slider( + "Высота сосуда (м)", + min_value=0.5, + max_value=3.0, + value=1.5, + step=0.5, + ) + + container_radius = st.slider( + "Радиус сосуда (м)", + min_value=0.5, + max_value=2.0, + value=1.0, + step=0.5, + ) + + initial_velocity = st.slider( + "Начальная скорость молекул (м/с)", + min_value=1.0, + max_value=30.0, + value=10.0, + step=1.0, + ) + + simulation_time = st.slider( + "Время симуляции (с)", + min_value=0.5, + max_value=5.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, + ) + + 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, + 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() + + 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 * config.k_B) + + max_speed = np.max(np.linalg.norm(final_velocities, axis=1)) + maxwell_speeds = np.linspace(0, max_speed * 1.2, 100) + maxwell_dist = calculate_maxwell_distribution( + maxwell_speeds, temperature, config.particle_mass + ) + + heights = np.linspace(0, config.container_height, 100) + scale_height = ( + config.k_B * temperature / (config.particle_mass * config.g) + ) + 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: + mean_height = np.mean(final_positions[:, 2]) + st.metric("Средняя высота", f"{mean_height:.3f} м") + with col4: + mean_speed = np.mean(np.linalg.norm(final_velocities, axis=1)) + st.metric("Средняя скорость", f"{mean_speed:.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() diff --git a/src/models/ideal_gas/charts.py b/src/models/ideal_gas/charts.py new file mode 100644 index 0000000..9237bec --- /dev/null +++ b/src/models/ideal_gas/charts.py @@ -0,0 +1,225 @@ +import altair as alt +import numpy as np +import pandas as pd + +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: + 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) + + return ( + 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)", + ) + ) + + +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) + + return ( + 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)", + ) + ) + + +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 + ) + + pe_initial = stats_df["potential_energy"].iloc[0] + + energy_df = pd.DataFrame({ + "time": stats_df["time"], + "Potential": ( + stats_df["potential_energy"] / abs(pe_initial) + if pe_initial != 0 + else stats_df["potential_energy"] + ), + }) + + energy_melted = energy_df.melt( + id_vars=["time"], var_name="energy_type", value_name="normalized_value" + ) + + return ( + 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=["Potential"], + range=["#ff7f0e"], + ), + ), + ) + .properties(width=600, height=400, title="Эволюция энергии системы") + ) + + +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 + ) + temp = 2.0 * kinetic_per_particle / (3.0 * config.k_B) + temps.append(temp) + + df = pd.DataFrame({"time": times, "temperature": temps}) + + return ( + 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="Эволюция температуры") + ) + + +def create_3d_particle_positions_chart( + positions: np.ndarray, + config: SimulationConfig, +) -> alt.Chart: + df = pd.DataFrame({ + "x": positions[:, 0], + "y": positions[:, 1], + "z": positions[:, 2], + }) + + height_scale = alt.Scale( + zero=False, domainMin=0, domainMax=config.container_height + ) + + return ( + 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=height_scale), + 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 сечение)") + ) diff --git a/src/models/ideal_gas/config.py b/src/models/ideal_gas/config.py new file mode 100644 index 0000000..4f44c2e --- /dev/null +++ b/src/models/ideal_gas/config.py @@ -0,0 +1,36 @@ +from dataclasses import dataclass + +MAX_DT = 0.01 + + +@dataclass +class SimulationConfig: + num_particles: int + + container_height: float + container_radius: float + + initial_velocity: float + + simulation_time: float + dt: float + + use_thermostat: bool = False + target_temperature: float | None = None + + g: float = 9.81 + k_B: float = 1.38e-23 # noqa: N815 + + initial_height: float = 0.1 + particle_mass: float = 6.63e-26 + particle_diameter: float = 3.4e-10 + + def __post_init__(self): + if self.num_particles <= 0: + raise ValueError("Количество молекул должно быть положительным") + if self.dt <= 0 or self.dt > MAX_DT: + raise ValueError(f"Шаг времени должен быть <= {MAX_DT}") + 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..117f975 --- /dev/null +++ b/src/models/ideal_gas/objects.py @@ -0,0 +1,296 @@ +import numpy as np + +from models.ideal_gas.config import SimulationConfig + + +class IdealGasSimulation: + 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 + + 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]: + 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) + + collided_with_walls = self._handle_collisions(positions, velocities) + + if ( + self.config.use_thermostat + and self.config.target_temperature is not None + ): + self._apply_thermostat_at_wall(velocities, collided_with_walls) + + return ( + np.array(positions_history), + np.array(velocities_history), + np.array(times), + ) + + def _initialize_particles( + self, positions: np.ndarray, velocities: np.ndarray + ) -> None: + num_particles = self.config.num_particles + + theta = np.random.uniform(0, 2 * np.pi, num_particles) + radius_factor = self.config.container_radius * 0.9 + r = np.sqrt(np.random.uniform(0, 1, num_particles)) * radius_factor + + positions[:, 0] = r * np.cos(theta) + positions[:, 1] = r * np.sin(theta) + positions[:, 2] = np.random.uniform( + 0, self.config.initial_height, num_particles + ) + + phi = np.random.uniform(0, 2 * np.pi, num_particles) + cos_theta = np.random.uniform(-1, 1, num_particles) + sin_theta = np.sqrt(1 - cos_theta**2) + + v0 = self.config.initial_velocity + velocities[:, 0] = v0 * sin_theta * np.cos(phi) + velocities[:, 1] = v0 * sin_theta * np.sin(phi) + velocities[:, 2] = v0 * cos_theta + + def _velocity_verlet_step( + self, positions: np.ndarray, velocities: np.ndarray + ) -> None: + dt = self.config.dt + g = self.config.g + + positions[:, 0] += velocities[:, 0] * dt + positions[:, 1] += velocities[:, 1] * dt + + positions[:, 2] += velocities[:, 2] * dt - 0.5 * g * dt**2 + + velocities[:, 2] -= g * dt + + def _handle_collisions( + self, positions: np.ndarray, velocities: np.ndarray + ) -> np.ndarray: + collided_with_walls = self._handle_wall_collisions( + positions, velocities + ) + + self._handle_intermolecular_collisions_optimized(positions, velocities) + + return collided_with_walls + + def _handle_wall_collisions( + self, positions: np.ndarray, velocities: np.ndarray + ) -> np.ndarray: + collided_indices = [] + + r = np.sqrt(positions[:, 0] ** 2 + positions[:, 1] ** 2) + collided = r > self.config.container_radius + + if np.any(collided): + collided_indices.extend(np.where(collided)[0]) + 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 + if np.any(collided_floor): + collided_indices.extend(np.where(collided_floor)[0]) + velocities[collided_floor, 2] *= -1 + positions[collided_floor, 2] = 0 + + collided_ceiling = positions[:, 2] > self.config.container_height + if np.any(collided_ceiling): + collided_indices.extend(np.where(collided_ceiling)[0]) + velocities[collided_ceiling, 2] *= -1 + positions[collided_ceiling, 2] = self.config.container_height + + return np.unique(collided_indices) + + def _handle_intermolecular_collisions_optimized( # noqa: PLR0912 + self, positions: np.ndarray, velocities: np.ndarray + ) -> None: + 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(): # noqa: PLR1702 + 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_at_wall( + self, velocities: np.ndarray, collided_indices: np.ndarray + ) -> None: + if not self.config.target_temperature or len(collided_indices) == 0: + return + + sigma = np.sqrt( + self.config.k_B + * self.config.target_temperature + / self.config.particle_mass + ) + velocities[collided_indices] = np.random.normal( + 0, sigma, size=(len(collided_indices), 3) + ) + + +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 + ) + temperature = self.calculate_temperature(velocities) + return n_density * self.config.k_B * temperature diff --git a/src/models/ideal_gas/utils.py b/src/models/ideal_gas/utils.py new file mode 100644 index 0000000..ac82616 --- /dev/null +++ b/src/models/ideal_gas/utils.py @@ -0,0 +1,71 @@ +import numpy as np +import pandas as pd + + +def calculate_maxwell_distribution( + speeds: np.ndarray, temperature: float, particle_mass: float +) -> np.ndarray: + 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) + + return norm_factor * speeds**2 * np.exp(exp_arg) + + +def calculate_barometric_formula( + heights: np.ndarray, scale_height: float +) -> np.ndarray: + if scale_height <= 0: + return np.zeros_like(heights) + + return (1.0 / scale_height) * np.exp(-heights / scale_height) + + +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] + + kinetic_energy = 0.5 * particle_mass * np.sum(vel**2) + + potential_energy = particle_mass * g * np.sum(pos[:, 2]) + + kinetic_per_particle = kinetic_energy / len(pos) + temperature = 2.0 * kinetic_per_particle / (3.0 * k_b) + + data["time"].append(times[step]) + data["kinetic_energy"].append(kinetic_energy) + data["potential_energy"].append(potential_energy) + data["total_energy"].append(kinetic_energy + potential_energy) + data["temperature"].append(temperature) + 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..e69de29 diff --git a/tests/m7b/conftest.py b/tests/m7b/conftest.py new file mode 100644 index 0000000..d0e5b8e --- /dev/null +++ b/tests/m7b/conftest.py @@ -0,0 +1,26 @@ +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..8e444bf --- /dev/null +++ b/tests/m7b/test_simulation.py @@ -0,0 +1,260 @@ +import numpy as np +from scipy.integrate import trapezoid + +from models.ideal_gas.config import SimulationConfig +from models.ideal_gas.objects import EquilibriumAnalyzer, IdealGasSimulation +from models.ideal_gas.utils import ( + calculate_barometric_formula, + calculate_maxwell_distribution, + calculate_speed_from_velocity, +) + + +def test_simulation_output_shapes(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(basic_config): + 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(basic_config): + 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" + ) + + +def test_total_energy_conserved(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, strict=False): + kinetic_energy = 0.5 * m * np.sum(vel_step**2) + potential_energy = m * g * np.sum(pos_step[:, 2]) + total_energies.append(kinetic_energy + potential_energy) + + total_energies = np.array(total_energies) + + initial_energy = total_energies[0] + + 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(basic_config): + sim = IdealGasSimulation(basic_config) + _, velocities, _ = sim.run() + + m = basic_config.particle_mass + kinetic_energies = [] + for vel_step in velocities: + kinetic_energy = 0.5 * m * np.sum(vel_step**2) + kinetic_energies.append(kinetic_energy) + + 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(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: + potential_energy = m * g * np.sum(pos_step[:, 2]) + potential_energies.append(potential_energy) + + 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(): + 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, strict=False): + kinetic_energy = 0.5 * m * np.sum(vel_step**2) + potential_energy = m * g * np.sum(pos_step[:, 2]) + total_energies.append(kinetic_energy + potential_energy) + + total_energies = np.array(total_energies) + initial_energy = total_energies[0] + + 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%}" + ) + + +def test_temperature_calculation(basic_config): + sim = IdealGasSimulation(basic_config) + _, velocities, _ = sim.run() + + analyzer = EquilibriumAnalyzer(basic_config) + temperature = analyzer.calculate_temperature(velocities[-1]) + + assert temperature > 0, "Temperature should be positive" + assert temperature < 10000, "Temperature seems unreasonably high" + + +def test_pressure_calculation(basic_config): + sim = IdealGasSimulation(basic_config) + positions, velocities, _ = sim.run() + + analyzer = EquilibriumAnalyzer(basic_config) + pressure = analyzer.calculate_pressure(positions[-1], velocities[-1], 0) + + assert pressure > 0, "Pressure should be positive" + + +def test_maxwell_distribution_normalization(): + speeds = np.linspace(0, 1000, 1000) + maxwell = calculate_maxwell_distribution(speeds, 300, 6.63e-26) + + integral = trapezoid(maxwell, speeds) + assert 0.9 < integral < 1.1, ( + f"Maxwell integral should be ~1, got {integral}" + ) + + +def test_barometric_formula_normalization(): + heights = np.linspace(0, 5, 100) + barometric = calculate_barometric_formula(heights, 1.0) + + integral = trapezoid(barometric, heights) + assert 0.9 < integral < 1.1, ( + f"Barometric integral should be ~1, got {integral}" + ) + + +def test_speed_calculation(): + 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]" + + +def test_simulation_doesnt_diverge(basic_config): + 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(basic_config): + sim = IdealGasSimulation(basic_config) + _, velocities, _ = sim.run() + + for vel_step in velocities: + kinetic_energy = np.sum(vel_step**2) + assert kinetic_energy >= 0, "Kinetic energy cannot be negative" + + +def test_gravity_pulls_downward(basic_config): + 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() + + initial_v_z = velocities[0, 0, 2] + final_v_z = velocities[-1, 0, 2] + + if initial_v_z <= 0: + velocity_changed = abs(final_v_z - initial_v_z) > 0.1 + + assert velocity_changed, "Gravity should affect velocity" + return + + velocity_decreased = False + for i in range(1, min(10, len(velocities))): + current_v_z = velocities[i, 0, 2] + if current_v_z < initial_v_z: + velocity_decreased = True + break + + assert velocity_decreased, "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" + )