{ "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 2

\n", "
\n", "

Wintersemester 23/24

\n", "\n", "\n", "**Website:** https://www.thp.uni-koeln.de/trebst/Lectures/2023-StatPhys.shtml\n", "\n", "**Abgabe**: keine Abgabe **Besprechung**: 23.10.2023 \n", "\n", "**Name**: Bitte geben Sie Ihren Namen an. " ] }, { "cell_type": "markdown", "id": "d8cd453a", "metadata": {}, "source": [ "\n", "\n", "
Aufgabe 7: Simulation eines zweidimensionalen Gases [0 Punkte]
" ] }, { "cell_type": "markdown", "id": "73c78fa0", "metadata": {}, "source": [ "In dieser Aufgabe sollen Sie sich mit der Simulation eines zweidimensionalen Gases vertraut machen. Für $N$ Teilchen sollen die klassischen Bewegungsgleichungen integriert, und elastische Kollisionen zwischen den Teilchen berücksichtigt werden. In diesem Notebook sind bereits alle wichtigen Funktionen implementiert. " ] }, { "cell_type": "markdown", "id": "778d53b3", "metadata": {}, "source": [ "**a)** Stellen Sie zunächst sicher, dass Sie die nachfolgenden Pakete installiert haben. Ergänzen Sie ggf. noch fehlende Pakete. Hilfe finden Sie z.B. hier." ] }, { "cell_type": "code", "execution_count": null, "id": "ba5a89c9", "metadata": {}, "outputs": [], "source": [ "using Random, LinearAlgebra, GLMakie" ] }, { "cell_type": "markdown", "id": "1b24c6ac", "metadata": {}, "source": [ "**b)** Führen Sie jetzt die nachfolgenden Zellen aus, in denen die für die Simulation wichtigen Funktionen definiert werden. Machen Sie sich die Wirkung der einzelnen Funktionen klar. Ergänzen Sie ggf. zusätzliche Kommentare." ] }, { "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", " return\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": "markdown", "id": "de83d77f", "metadata": {}, "source": [ "**c)** Führen Sie die nachfolgende Zelle aus, um die Simulation zu starten. Testen Sie anschließend auch andere Anfangsbedingungen, z.B.\n", "- Alle Teilchen starten in einer Ecke der Box.\n", "- Die Hälfte der Teilchen ruht zu Beginn.\n", "- Nur ein einziges Teilchen hat zu Beginn eine (sehr hohe) kinetische Energie.\n", "- Ein Teilchen ist 1000x schwerer oder 10x größer als alle anderen.\n", "\n", "Können Sie die unterschiedlichen Anfangsbedingungen bei großen Zeiten noch auseinanderhalten?" ] }, { "cell_type": "code", "execution_count": null, "id": "30cce5bf", "metadata": {}, "outputs": [], "source": [ "### Simulation Parameters ###\n", "N = 100 # Number of particles\n", "box = (0.0, 2, 0.0, 1.5) # Bounding box of simulation\n", "radius = 0.01 # Radius of particles\n", "v = 0.4 # Initial velocity\n", "dt = 0.01 # Time step in update.\n", "steps_per_frame = 10 # Number of updates before visualization is updated\n", "\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", "\n", "### Initialize particles ###\n", "initialize_particles_without_overlap!(x,y,vx,vy,r,v,box)\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", "\n", "### Set up animation ###\n", "fig = Figure()\n", "ax = Axis(fig[1,1], aspect = 1.333, spinewidth = 3)\n", "limits!(ax, box[1], box[2], box[3], box[4])\n", "hidedecorations!(ax)\n", "scatter!(ax, x, y, color = vabs, colormap = :viridis, markersize = 1000 * r, colorrange = (0, 2*v))\n", "display(fig)\n", "\n", "### Start/Stop Button ###\n", "run = Button(fig[2,1]; label = \"Start / Stop\", tellwidth = false)\n", "isrunning = Observable(false)\n", "on(run.clicks) do clicks; isrunning[] = !isrunning[]; 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.03)\n", " end \n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "551dd900", "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 }