Valversnelling in Chimay

[Xilvo]

Kan je de richting waarin de gravitatie gaat staan op het aardoppervlak beredeneren vanuit symmetrie-overwegingen (uitgaande van een ellipsoïde)?

Loodrecht op het oppervlak, want dat moet een equipotentiaalvlak zijn.

[wnvl1]

Dat was ik vergeten. Maar dan is de centripetale component wel meegerekend, dus je maakt wel een foutje tov de integralen zoals ik ze nu opgeschreven heb.

[Xilvo]

Dat was ik vergeten. Maar dan is de centripetale component wel meegerekend, dus je maakt wel een foutje tov de integralen zoals ik ze nu opgeschreven heb.

Klopt, dat was ik dan weer vergeten.

Je kunt er wel omheen. Je weet hoe groot de totale versnelling bij benadering is, de centrifugale component kun je makkelijk berekenen dus daaruit is de richting van de zwaartekrachtcomponent behoorlijk nauwkeurig te bepalen.
Heel elegant is het natuurlijk niet.

[wnvl1]

Dat kan inderdaad.

Een andere bedenking:
De aarde is symmetrisch rondom de z-as. Als ik dan een punt op de aardbol neem kan ik via symmetriebeschouwingen voor een bepaald punt op de aarde uit de x-component van de gravitatie ook de y-component berekenen. De totale gravitatie moet loodrecht staan op het oppervlak, dus als ik de x-component en de y-component heb, dan heb ik direct ook de z-component. Dus 1 van mijn 3 integralen volstaat.

[Xilvo]

Als de evenaar in het x-y vlak ligt, dan kan je het punt op de aarde altijd in het x-z vlak kiezen. Wegens symmetrie middelen de y-componenten uit tot nul.
Dan kan je inderdaad met één component volstaan. Is dat de x-component, dan wordt de nauwkeurigheid bij de pool waarschijnlijk wel heel slecht; vice versa voor het kiezen voor de z-component en de nauwkeurigheid aan de evenaar.

[wnvl1]

Als ik mijn punt 2 aardstralen van het centrum van de aarde verwijderd leg, dan werkt mijn code heel snel en vermoedelijk is ze juist. De waarde ligt relatief dicht in de buurt van GMaarde/(2r)^2.

Numeriek zijn die integralen dan natuurlijk veeeeel stabieler. Het probleem dat de noemer naar nul gaat dicht bij het aardoppervlak is dan weg.

Code: Selecteer alles

import numpy as np
from scipy.integrate import tplquad

a=6356000 #straal bij pool
b=6378000 #straal bij evenaar
c=6378000 #straal bij evenaar

G=6.6743e-11 

r_loc=2 #positie Chimay
phi_loc=0 #positie Chimay
theta_loc=np.radians(90) #positie Chimay

rho = 5515 #massa dichtheid aarde

#gravitatie: x-component
gravx = tplquad(lambda phi, theta, r: rho*r**2 *a*b*c* np.sin(theta)*G*1/((r*a*np.sin(phi)*np.cos(theta) -  r_loc*a*np.sin(phi_loc)*np.cos(theta_loc))**2 +(r*b*np.sin(phi)*np.sin(theta) - r_loc*b*np.sin(phi_loc)*np.sin(theta_loc))**2+(r*c*np.cos(phi) - r_loc*c*np.cos(phi_loc))**2 )**1.5 *(r*a*np.sin(phi)*np.cos(theta) -  r_loc*a*np.sin(phi_loc)*np.cos(theta_loc))  ,
            0, 1,
            lambda theta: 0, np.pi,
            lambda theta, phi: 0, 2*np.pi,
            epsabs=0.000001, epsrel=0.000001)

print('Gravitatie x: ' +  str(gravx[0]) )

#gravitatie: y-component
gravy = tplquad(lambda phi, theta, r: rho*r**2 *a*b*c* np.sin(theta)*G*1/((r*a*np.sin(phi)*np.cos(theta) -  r_loc*a*np.sin(phi_loc)*np.cos(theta_loc))**2 +(r*b*np.sin(phi)*np.sin(theta) - r_loc*b*np.sin(phi_loc)*np.sin(theta_loc))**2+(r*c*np.cos(phi) - r_loc*c*np.cos(phi_loc))**2 )**1.5 *(r*b*np.sin(phi)*np.sin(theta) - r_loc*b*np.sin(phi_loc)*np.sin(theta_loc))  ,
            0, 1,
            lambda theta: 0, np.pi,
            lambda theta, phi: 0, 2*np.pi,
            epsabs=0.000001, epsrel=0.000001)

print('Gravitatie y: ' +  str(gravy[0]) )

#gravitatie: z-component
gravz = tplquad(lambda phi, theta, r: rho*r**2 *a*b*c* np.sin(theta)*G*1/((r*a*np.sin(phi)*np.cos(theta) -  r_loc*a*np.sin(phi_loc)*np.cos(theta_loc))**2 +(r*b*np.sin(phi)*np.sin(theta) - r_loc*b*np.sin(phi_loc)*np.sin(theta_loc))**2+(r*c*np.cos(phi) - r_loc*c*np.cos(phi_loc))**2 )**1.5 *(r*c*np.cos(phi) - r_loc*c*np.cos(phi_loc) )  ,
            0, 1,
            lambda theta: 0, np.pi,
            lambda theta, phi: 0, 2*np.pi,
            epsabs=0.1, epsrel=0.1)

print('Gravitatie z: ' +  str(gravz[0]) )

print((gravx[0]**2+gravy[0]**2+gravz[0]**2)**0.5)

Ruwe controle:

Code: Selecteer alles

print(G*4/3*np.pi*a**3*rho/(2*a)**2)

Veel mogelijkheden om het op te lossen heb ik vermoedelijk niet met scipy. Je moet dan vermoedelijk uitwijken naar oplossing genre die van xilvo.

Misschien doet Mathcad beter?

[HansH]

Misschien doet Mathcad beter?

ik ben nu zo ver dat ik een functie F heb gedefinieerd die voor een punt p2 met massabijdrage dm de bijdrage aan g berekent in punt p1.
in die functie bereken ik dan eerst het kwadraat van de afstand tussen p1 en p2 en dan de bijdrage in de kracht F1 en dan de richtingsvector r waarin die kracht wijst. De uitkomst van die functie is dus een 3d vector met lengte gelijk aan de zwaartekrachtsbijdrage in p1 tgv massa element dm met coordinaten p2.
volgende stap is dan die bijdrage in de integraal zetten zodat de integraal dan voor elk stukje dm een bijdrage berekent aan de zwaartekrachtsvector en dan alles sommeert tot 1 totale vector die dan als het goed is de richting van de zwaartekracht van de hele aarde geeft en de grootte daarvan in punt p1.
ik ben alleen nog aan het puzzelen hoe ik dat in die integraal moet stoppen.

Mathcad - zwaartekracht3

(161.1 KiB) 23 keer gedownload

[wnvl1]

Je kan een aparte integraal opstellen zoals ik gedaan heb voor x-richting, y-richting en z-richting; en dan de 3 componenten met Pythagoras optellen. Dat gaat het simpelste zijn en is hoogstens maar een factor 3 trager dan alles in 1 keer optellen.

[HansH]

inmiddels lijkt het te werken.
in kan voor een willekurig punt buiten de aarde de zwaartekracht uitrekenen met grootte en richting.

Mathcad - zwaartekracht4

(164.17 KiB) 27 keer gedownload

berekening duurt maar een paar seconden per punt

[HansH]

wat me wel opvalt is dat bij punten op grotere afstand van de aarde de berekening maar heel kort duurt, maar als ik bv 2 km boven het oppervlak bereken of 100 meter of 1 meter dan wordt het steeds trager en op het oppervlak komt het algoritme er niet uit, blijft dan rekenen. Ik denk dat dat komt omdat er dan een stukje van de integraal voorbij komt waar de afstand naar 0 gaat dus je hele grote fouten kan maken als je dx,dy,dz te groot maakt. de tool bepaalt zelf hoeveel rekenpunten het moet nemen om de opgegeven nauwkeurigheid te halen, dus gaat dan heel lang duren denk ik en je kunt grote fouten maken.

maar dat probleem heb je ook bij de balk methode lijkt mij of je moet punten uit gaan sluiten die te dicht bij elkaar liggen.

[Xilvo]

Uiteraard krijg je problemen als je te dicht bij het oppervlak komt bij discrete massa-elementen.
Ik heb een hoogte van 0,1% van de aardstraal gekozen, ca 6 km.
Bij een globe met een diameter van 2 m zit je dan 1 mm boven het oppervlak.

Verder heb ik massa-elementen dicht bij het oppervlak gewogen meegerekend. Dat werkt prima. Rekentijd minder dan 30 s per punt.

[HansH]

De integraal is een black box berekening in mathcad. dus je kunt niet in dat algoritme ingrijpen vrees ik om bepaalde punten over te slaan. Dan zou ik de integtraal methode opnieuw moeten definieren via programming. dat kan wel maar is mogelijk trager?

Wat ik wel simpel kan doen is het zwaartekrachtsveld berekenen bv van noordpool (0,0,b) naar evenaar (a,0,0) dus bv in het x, z vlak op bv 1km, 2 km, 3km boven het aardoppervlak en op basis daarvan dan extrapoleren naar het aardoppervlak zelf. als ik even tijd heb zal ik daar verder naar kijken.

Voordeel van mathcad is dat het allemaal overzichtelijk modulair is dus je kunt ook heel makkelijk spelen met de massadichtheidsfunctie bijvoorbeeld door even een andere daarvoor in de plaats te zetten.

[Xilvo]

Wat ik wel simpel kan doen is het zwaartekrachtsveld berekenen bv van noordpool (0,0,b) naar evenaar (a,0,0) dus bv in het x, z vlak op bv 1km, 2 km, 3km boven het aardoppervlak en op basis daarvan dan extrapoleren naar het aardoppervlak zelf. als ik even tijd heb zal ik daar verder naar kijken.

Is niet nodig, daar kan je de hoogtecorrectie die hier gegeven wordt voor gebruiken.
Die volgt direct uit \(g(h)=\frac{G.M}{(R+h)^2}\)

[wnvl1]

In de documentatie van scipy staat:

https://docs.scipy.org/doc/scipy/refere ... lquad.html

For valid results, the integral must converge; behavior for divergent integrals is not guaranteed.

Ik vermoed dat ze bedoelen "improper integrals". Divergente integralen kan je toch nooit berekenen?
We hebben hier een oneigenlijke convergerende integraal.

[HansH]

ik heb de sheet nu uitgebreid met de berekening van de zwaartekracht (exclusief centripetale kracht) waarbij ik een pad volge van noordpool naar evenaar op een factor k1 keer de ellipsvorm.

als ik dan de zwaartekracht met k1^2 vermenigvuldig dan heb ik de vertaling naar zwaartekracht aan het oppervlak. het blijkt dat dat voor k1=1.005,1.01 en 1.015 vrijwel dezelfde grafiek oplevert.
op zich ook logisch want de zwaartekracht neemt af met het kwadraat dan de afstand tot het centrum voor een bol en dit is bijna een bol.

zwaartekracht

(531.93 KiB) 18 keer gedownload