{ "cells": [ { "cell_type": "markdown", "id": "c9115fd9", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "\n", "\n", "
Institut für Theoretische Physik
Universität zu Köln 		 Prof. Dr. Simon Trebst
Christoph Berke
\n", "
\n", "

Statistische Physik

\n", "

Übungsblatt 6

\n", "
\n", "

Wintersemester 23/24

\n", "\n", "\n", "**Website:** https://www.thp.uni-koeln.de/trebst/Lectures/2023-StatPhys.shtml\n", "\n", "**Abgabe**: 20.11.2023,10:00 Uhr **Besprechung**: 21.11.2023 \n", "\n", "**Name**: Bitte geben Sie Ihren Namen an. " ] }, { "cell_type": "markdown", "id": "8c379fc5", "metadata": {}, "source": [ "\n", "\n", "
Aufgabe 22: Maxwellsche Geschwindigkeitsverteilung [3 Bonuspunkte]
" ] }, { "cell_type": "markdown", "id": "73c78fa0", "metadata": {}, "source": [ "In dieser Aufgabe sollen Sie, aufbauend auf der Simulation des zweidimensionalen Gases, demonstrieren, dass die Verteilung der Geschwindigkeiten in einem Gas durch die Maxwellsche Geschwindigkeitsverteilung (auch Maxwell-Boltzmann-Verteilung) gegeben ist. Nachfolgend finden Sie die Ihnen bereits bekannten Funktionen, die zum Erzeugen einer Anfangskonfiguration und der Simulation elastischer Stöße benötigt werden." ] }, { "cell_type": "code", "execution_count": null, "id": "ba5a89c9", "metadata": {}, "outputs": [], "source": [ "using Random, LinearAlgebra, GLMakie, Statistics" ] }, { "cell_type": "code", "execution_count": null, "id": "a2c40f09", "metadata": {}, "outputs": [], "source": [ "# Initialize positions and velocities of particles in a 2D box.\n", "# x, y = Position arrays,\n", "# vx, vy = velocity arrays,\n", "# r = particle radius array\n", "# v = particle velocity (abs. value)\n", "# box = tuple with boundaries\n", "function initialize_particles!(x, y, vx, vy, r, v, box)\n", " \n", " # extract boundaries\n", " xmin, xmax, ymin, ymax = box\n", " \n", " # initialize each particle individually\n", " for i in 1:length(x)\n", " \n", " # random position with distance r to boundary\n", " x[i] = xmin+2*r[i] + rand() * (xmax - xmin - 4*r[i])\n", " y[i] = ymin+2*r[i] + rand() * (ymax - ymin - 4*r[i])\n", " \n", " # random velocity direction\n", " angle = rand()*2*pi\n", " vx[i], vy[i] = cos(angle)*v, sin(angle)*v\n", " \n", " end\n", " \n", "end\n", "\n", "# Same as above, but take account of the finite dimension of each particle. \n", "function initialize_particles_without_overlap!(x, y, vx, vy, r, v, box)\n", " \n", " # extract boundaries\n", " xmin, xmax, ymin, ymax = box\n", " \n", " # initialize each particle individually\n", " for i in 1:length(x)\n", " \n", " # reset counter of maximal tries\n", " tries, maxtries = 0, 1000\n", " \n", " # random position with distance r to boundary and no overlap to other particles\n", " overlapping = true\n", " while overlapping\n", " x[i] = xmin+2*r[i] + rand()*(xmax - xmin - 4*r[i])\n", " y[i] = ymin+2*r[i] + rand()*(ymax - ymin - 4*r[i])\n", " overlapping = false\n", " for j in 1:i-1\n", " if (x[i]-x[j])^2 + (y[i]-y[j])^2 < 1.01*(r[i]+r[j])^2\n", " overlapping = true\n", " end\n", " end\n", " tries += 1\n", " # maybe break loop because of maximum tries\n", " if tries > maxtries\n", " @error \"Could not place particle! Exceeded $(max_tries) tries\"\n", " return\n", " end\n", " end\n", " \n", " # random velocity direction\n", " angle = rand()*2*pi\n", " vx[i], vy[i] = cos(angle)*v, sin(angle)*v\n", " \n", " end\n", " \n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "b8a14600", "metadata": {}, "outputs": [], "source": [ "function move_particles!(x, y, vx, vy, dt, r, box)\n", " \n", " # move all particles by dt\n", " x .+= vx .* dt\n", " y .+= vy .* dt\n", " \n", " # check collision with bounding box\n", " for i in 1:length(x)\n", " # left / right\n", " if x[i]-r[i] < box[1]\n", " vx[i] *= -1\n", " elseif x[i]+r[i] > box[2]\n", " vx[i] *= -1\n", " end\n", " # up / down\n", " if y[i]-r[i] < box[3]\n", " vy[i] *= -1\n", " elseif y[i]+r[i] > box[4]\n", " vy[i] *= -1\n", " end\n", " end\n", " return 0\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "2a856eb1", "metadata": {}, "outputs": [], "source": [ "function collide_particles!(x, y, vx, vy, m, r)\n", " \n", " # check all particle pairs (i,j) for collision\n", " for i in 1:length(x)\n", " for j in 1:i-1\n", " # particles i and j are too close --> collision\n", " if (x[i]-x[j])^2 + (y[i]-y[j])^2 < (r[i]+r[j])^2\n", " er = [x[j]-x[i], y[j]-y[i]]\n", " er = er ./ norm(er)\n", " vi = er[1]*vx[i] + er[2]*vy[i]\n", " vj = er[1]*vx[j] + er[2]*vy[j]\n", " vip = 2*(m[i]*vi + m[j]*vj)/(m[i]+m[j]) - vi\n", " vjp = 2*(m[i]*vi + m[j]*vj)/(m[i]+m[j]) - vj\n", " vx[i] = vx[i] + (vip - vi)*er[1]\n", " vy[i] = vy[i] + (vip - vi)*er[2]\n", " vx[j] = vx[j] + (vjp - vj)*er[1]\n", " vy[j] = vy[j] + (vjp - vj)*er[2]\n", " end\n", " \n", " end\n", " end\n", " \n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "f80aac38", "metadata": {}, "outputs": [], "source": [ "# Maxwell-Boltzmann distribution in two dimensions.\n", "function maxwell(v, m, beta)\n", " return m * beta * v * exp(-m * v^2 * beta / 2)\n", "end" ] }, { "cell_type": "markdown", "id": "5fcc62ce", "metadata": {}, "source": [ "## Jetzt sind Sie dran:\n", "Nachfolgende Zelle erzeugt die Ihnen bereits bekannte Simulation. Zusätzlich ist jetzt auch noch die Verteilung der Geschwindigkeitsbeträge $|v|$ dargestellt. Die Anfangsbedingungen sind so gewählt, dass alle Teilchen zu Beginn dasselbe $|v|$ haben. Implementieren Sie andere Anfangsbedingungen etwa\n", "- jedes zweite Teilchen ruht,\n", "- zu Beginn sind die Geschwindigkeiten gleichförmig auf einem geeignet gewählten Intervall verteilt,\n", "- ...\n", "\n", "und überzeugen Sie sich davon, dass die Konvergenz gegen die Maxwellsche Geschwindigkeitsverteilung unabhängig von der Anfangskonfiguration erfolgt.\n", "Stellen Sie in dem Histogram zusätzlich die erwartete Verteilung ```maxwell(v, m, beta)``` dar (siehe obige Zelle). Überlegen Sie sich, wie Sie aus dem mittleren $|v|$ das entsprechende $\\beta$ berechnen können." ] }, { "cell_type": "code", "execution_count": null, "id": "30cce5bf", "metadata": {}, "outputs": [], "source": [ "### Simulation Parameters ###\n", "N = 1000 # Number of particles\n", "box = (0.0, 6, 0.0, 3) # Bounding box of simulation\n", "radius = 0.002 # Radius of particles\n", "v = 1.5 # Initial velocity\n", "dt = 0.01 # Time step in update.\n", "steps_per_frame = 10 # Number of updates before visualization is updated\n", "\n", "### Initial conditions ###\n", "x, y = zeros(N), zeros(N) # Arrays for positions.\n", "vx, vy = zeros(N), zeros(N) # Arrays for velocities\n", "m, r = ones(N), radius * ones(N) # Arrays for masses and radii\n", "\n", "### Initialize all particles with same |v| ###\n", "initialize_particles_without_overlap!(x,y,vx,vy,r,v,box)\n", "\n", "# Make the relevant arrays observables.\n", "x, y, vx, vy = Observable(x), Observable(y), Observable(vx), Observable(vy)\n", "vabs = @lift( sqrt.($vx.^2 + $vy.^2) ) # Observable that holds all abs. values of velocities.\n", "\n", "### Set up animation ###\n", "fig = Figure(resolution = (1300,600))\n", "ax = Axis(fig[1,1:2], aspect = 1.333, spinewidth = 3, backgroundcolor = :gray70)\n", "ax2 = Axis(fig[1,3:4], aspect = 1.333, xlabel = L\"v\", ylabel = L\"p(v)\")\n", "hist!(ax2, vabs, normalization = :pdf, bins = range(0,3*v,length = 30))\n", "\n", "limits!(ax, box[1], box[2], box[3], box[4])\n", "hidedecorations!(ax)\n", "scatter!(ax, x, y, color = vabs, colormap = :curl, markersize = 4500 * r, colorrange = (0, 3*v))\n", "display(fig)\n", "\n", "### Start/Stop Button ###\n", "run = Button(fig[2,1]; label = \"Start / Stop\", tellwidth = false, tellheight = true)\n", "isrunning = Observable(false)\n", "on(run.clicks) do clicks; isrunning[] = !isrunning[]; end\n", "\n", "# Restart Button\n", "restart = Button(fig[2,2]; label = \"Restart\", tellwidth = false, tellheight = true)\n", "on(restart.clicks) do _\n", " initialize_particles_without_overlap!(x[],y[],vx[],vy[],r,v,box)\n", " notify(x);notify(y);notify(vx);notify(vy)\n", "end\n", "\n", "### Link Button to Simulation ###\n", "on(run.clicks) do _\n", " @async while isrunning[]\n", " isopen(fig.scene) || break # Stoppen wenn Fenster geschlossen ist.\n", " for j in 1:steps_per_frame\n", " # let particles collide with each other\n", " collide_particles!(x[],y[],vx[],vy[],m,r)\n", " # let particles move\n", " move_particles!(x[],y[],vx[],vy[],dt/steps_per_frame,r,box)\n", " end\n", "\n", " notify(x);notify(y);notify(vx);notify(vy)\n", " sleep(0.01)\n", " end \n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "af8d2a76", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.9.3", "language": "julia", "name": "julia-1.9" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.9.3" } }, "nbformat": 4, "nbformat_minor": 5 }