Ik had beloofd dat ik wat zou vertellen over het scriptje dat ik geschreven heb om de zwaartekracht van twee objecten op elkaar te berekenen, dus bij dezen.
TL;DR: hak je object in stukjes, bereken de kracht van elk stukje op elk ander stukje, tel ze op, klaar. Het werkt, want de berekende kracht komt heel goed overeen met de theoretische kracht voor een bol en een holle bol.
Theorie:
Het basisidee is vrij eenvoudig. Je hebt een bepaalde massa, bijvoorbeeld een cilinder of een bol. Als je hiervan de zwaartekracht wil berekenen, deel je die massa op in (liefst zoveel mogelijk) kleine elementjes. Uiteraard moet je de positie en de massa van die elementjes zo kiezen dat ze een goede representatie zijn van het object dat je wil simuleren. Vervolgens neem je een andere puntmassa als 'testmassa'. De kracht tussen twee puntmassa's kun je simpelweg berekenen met de gravitatiewet van Newton.
Voor alle elementjes bereken je dan de kracht op je testmassa. Die krachten tel je bij elkaar op om de totale kracht van het voorwerp op je testmassa te krijgen.
Om de kracht van de ene massa op de andere te krijgen, moet je ze allebei in elementjes opdelen. Vervolgens bereken je voor elk elementje in massa 1 de totale kracht van alle elementen in massa 2. Vervolgens tel je de krachten van alle elementen in massa 1 bij elkaar op om de totale kracht van massa 2 op massa 1 te krijgen.
Ik gebruik cilindercoördinaten
\((r, \theta, z)\)
om de cilinder te beschrijven. De straal van de cilinder is
R, de hoogte is
h, de massa is
m. Het domein van de coördinaten is dan dus:
\(r \in [0, R] \quad \theta \in [0, 2\pi] \quad z \in [0,h]\)
Een volume-elementje van een cilinder is (zie
hier):
\(\mathrm{d}V = r\mathrm{d}r\,\mathrm{d}\theta \,\mathrm{d}z\)
De massa van dit elementje is de dichtheid
rho keer dat volume.
\(\mathrm{d}M= \rho\,\mathrm{d}V = \rho \, r \mathrm{d}r \,\mathrm{d}\theta\,\mathrm{d}z\)
met
\( \rho = \frac{m}{\pi r^2 h}.\)
De positie van een elementje is:
\(\vec{x} = \vec{x_0} + \left(\begin{matrix} r \cos(\theta) \\ r \sin(\theta) \\ z \end{matrix} \right)\)
De kracht op een puntmassa
m door dat elementje is:
\(\mathrm{d}\vec{F} = G \frac{m\,\mathrm{d}M}{s^2} \hat{s}\)
met s de afstand tussen de twee en s-dakje de eenheids-richtingsvector.
Stel ik neem N_r punten in de radiale richting, N_theta in de angulaire, en N_z in de longitudinale, dan wordt de totale kracht op een puntmassa:
\(\vec{F} = \sum\limits_{i=1}^{N_\theta}\sum\limits_{j=1}^{N_r}\sum\limits_{k=1}^{N_z}\mathrm {d}\vec{F}\)
Die kracht bereken ik dus voor alle punten in het voorwerp
waarop ik de kracht wil weten, en die sommeer ik dan.
Implementatie
De code is geschreven in Python 3.6.1 en maakt gebruik van o.a.de libraries numpy, scipy en matplotlib.
In
bericht 138 kun je een plotje zien van hoe deze verdeling in elementen eruit ziet voor verschillende N's.
Hieronder 2 voorbeelden van een resultaat van een berekening:
Code: Selecteer alles
============================= Simulation complete ==============================
N_r: 10 N_theta: 60 N_h: 10
Simulation force vector [ 1.20866740e-07 8.38270579e-24 -4.07534451e-26]
Simulation, total force: 1.2086674013576318e-07N
Newton, total force: 1.4288357777777776e-07N
Simulation completed in 174.2 seconds.
============================= Simulation complete ==============================
N_r: 10 N_theta: 100 N_h: 20
Simulation force vector [ 1.20684476e-07 1.04925070e-23 -2.22291519e-26]
Simulation, total force: 1.2068447606364614e-07N
Newton, total force: 1.4288357777777776e-07N
Simulation completed in 1912.9 seconds.
De 'Simulation force vector' is de vector van de totale kracht van cilinder 1 op cilinder 2. De cilinder staat op 8.4 cm in de x-richting van de ander. De kracht wijst dus alléén in de x-richting. In de y- en z-richting middelt de kracht uit naar (praktisch) 0. 'Simulation, total force' is de lengte van die vector.
'Newton, total force' is ter vergelijking de kracht van 2 puntmassa's met dezelfde totale massa als de cilinder op dezelfde (middelpunts)afstand. De ordegrootte van deze berekeningen lijkt dus in ieder geval in orde.
Merk op dat voor de laatste simulatie één cilinder bestaat uit 10*100*20 = 20.000 punten. Het totaal aantal krachten dat dus berekend werd voor deze simulatie is dan dus (20.000)
2 = 400 miljoen! Dat komt dus neer op ruim 20.000 berekeningen per seconde (en inmiddels heb ik nog wat optimalisatie doorgevoerd.)
Als iemand interesse heeft om de broncode in te zien, stuur me dan een pb.
Validatie
Allemaal leuk en aardig, maar klopt er ook iets van? Helaas heb ik geen referentie van de zwaartekracht van een cilinder (anders was deze berekening ook niet nodig geweest). Gelukkig weet ik wel de kracht van een bol: die word namelijk beschreven door het zgn. '
shell theorem'.
<b>Berekening van de zwaartekracht van een bol</b>
Met behulp van
bolcoördinaten kan ik een bol ook opdelen in elementjes. Vervolgens heb ik de kracht op een puntmassa voor verschillende posities op de z-as berekend. De bol had een massa van 1 kg, de puntmassa 10 g. De straal van de bol was 1 m en deze was opgedeeld in 40.000 elementjes (N_r = 20, N_phi = 20, N_theta = 100).
Hier is een voorbeeld van hoe de punten verdeeld waren. De groene lijn zijn de posities van de testmassa waarvoor de kracht berekend is.
- sphere_points 636 keer bekeken
Het resultaat zie je hier:
- sphere_force 636 keer bekeken
Op de x-as staat
r, in dit geval de afstand van de puntmassa tot het midden van de bol; op de y-as de kracht in Newton.
De resultaten van de theorie en de simulatie komen zéér goed overeen. Alleen in het midden van de bol klopt het niet helemaal, een gevolg van het eindig aantal elementen.
Berekening van de zwaartekracht van een holle bol
Om te controleren dat de berekening klopt, gebruiken we het shell theorem voor een holle bol. Binnen deze holle bol moet de totaalkracht op het deeltje 0 zijn. De bol is heeft weer een massa van 1 kg, en een straal van 1 m. De binnenstraal is 0.5 m. De bol is weer opgedeeld in 40.000 punten. Ook de puntmassa is hetzelfde.
De verdeling van de punten:
- hollow_sphere_points 636 keer bekeken
- hollow_sphere_force 636 keer bekeken
In de grafiek is duidelijk te zien dat de kracht inderdaad netjes naar 0 gaat binnen de holle bol.
- hollow_sphere_force_detail 636 keer bekeken
Een close up met wat meer berekende punten.
Conclusie
De resultaten van de simulatie voor een bol en een holle bol komen zeer goed overeen met de theorie. Het is dan ook aannemelijk dat het model ook goede resultaten geeft voor simulaties van een cilinder.
Outlook
Vanwege symmetrie hoef ik slechts een kwart van de cilinder waarop ik de kracht wil berekenen te simuleren: bijv. de bovenste helft van de rechterhelft. Door dit kwart te 2 keer te spiegelen kan ik dan de totale kracht berekenen. Dit moet ik nog implementeren. Voordeel hiervan is dat de simulatie met meer elementen gedraaid kan worden.
Daarnaast ben ik aan het nadenken over een andere manier om de cilinder in elementjes te verdelen, waarbij de elementen voor grotere
r niet zo 'grof' zijn.