Valversnelling in Chimay

[wnvl1]

Je rekent de integraal van \(Gdm/r^2\vec{e_r}\) uit over de volledige ellipsoide met r de afstand van een punt in de ellipsoide tot het beschouwde punt waar je g wilt berekenen. Om die integraal uit te rekenen werk je met de vergelijking van de ellipsoide. Je verrekent dus de volledige vorm van de ellips (alle tussenliggende punten worden in rekening gebracht).

Dat is een leuke oefening voor een cursus differentiaal en integraal rekening. Misschien een idee voor een volgend raadsel.

Hier een aanzet voor het soort van integraal dat je moet berekenen in geval van een sfeer.

https://www.mit.edu/~hlb/1802/notes-apm ... otes_g.pdf

De opmerking (remark) in de link is hier belangrijk bij ellipsoïdes.

[Xilvo]

hieruit blijkt dat je de verhouding hebt berekend tussen equatoriale en polaire diameter. maar er blijkt niet uit dat je niet ook tussenliggende punten berekend zou kunnen hebben. (als ik tenminste goed lees)

Nee. Die diameters en de verhouding ervan zijn bekend.

[Xilvo]

De opmerking (remark) in de link is hier belangrijk bij ellipsoïdes.

Klopt.
Ik heb er over gedacht de numerieke integratie wat slimmer aan te pakken, bijvoorbeeld door bij het bekende veld van een bol alleen de extra bijdrage van alles wat binnen de ellipsoïde maar buiten de ingeschreven bol valt numeriek te integreren.
Vanuit de pool levert dat mooi allemaal ringen op, vanuit de evenaar is dat weer lastiger.

Dus heb ik eerst maar eens gekeken of met de meest simpele manier van numeriek integreren een voldoende nauwkeurig resultaat binnen aanvaardbare tijd te bereiken was. Dat bleek mogelijk.

[HansH]

ik wil graag een poging doen om via mathcad met numerieke berekening het zwaartekrachtsveld van een bol en ellipsoide te vergelijken met de gegeven afplatting als functie van de hoek tov de polen. Ik begreep dat er in dit topic al wat werk was gedaan hiervoor dus zou mooi zijn als de code en globale beschrijving van die code dat gedeeld kan worden zodat ik daarop kan voortbouwen en om te voorkomen dat ik het wiel deels opnieuw moet uitvinden. (enige beschrijving van de gebruikte principes om reverse engineering te voorkomen)

[HansH]

Dat is een leuke oefening voor een cursus differentiaal en integraal rekening. Misschien een idee voor een volgend raadsel.

Hier een aanzet voor het soort van integraal dat je moet berekenen in geval van een sfeer.

https://www.mit.edu/~hlb/1802/notes-apm ... otes_g.pdf

in mathcad stelt deze berekening numeriek niet heel veel voor denk ik. verder optimaliseren door slim alleen het deel mee te nemen wat verschillend is tov een bol kan, maar kost waarschijnlijk meer tijd om het uit te werken dan om het met domme rekenkracht helemaal door te rekenen voor een complete ellipsoide. dus al ik tijd heb zal ik daar eens verder naar kijken.

[HansH]

https://www.mit.edu/~hlb/1802/notes-apm ... otes_g.pdf

We concentrate on finding just the k component of the gravitational attraction

wat 'the k component ' is wordt echter niet uitgelegd, en kan ik ook niet zo uit de context opmaken. dus weet iemand wat
'the k component ' is ?

[Xilvo]

ik wil graag een poging doen om via mathcad met numerieke berekening het zwaartekrachtsveld van een bol en ellipsoide te vergelijken met de gegeven afplatting als functie van de hoek tov de polen. Ik begreep dat er in dit topic al wat werk was gedaan hiervoor dus zou mooi zijn als de code en globale beschrijving van die code dat gedeeld kan worden zodat ik daarop kan voortbouwen en om te voorkomen dat ik het wiel deels opnieuw moet uitvinden. (enige beschrijving van de gebruikte principes om reverse engineering te voorkomen)

Neem een balk waar de ellipsoïde net in past, voorzie die van een rooster van massapunten, bijvoorbeeld 2000 per zijde.
Kies een punt net buiten de ellipsoïde waarvoor je de zwaartekrachtsversnelling wil weten.
Bepaal voor ieder massapunt dat binnen de ellipsoide ligt de afstand d tot het gekozen punt en daaruit de richting en grootte (∝ 1/d²) van de zwaartekrachtsversnellingsvector en sommeer die vectoren voor alle massapunten.
Als het alleen om de relatieve waardes gaat kun je de massa van de massapunten 1 kiezen, net als G.
Dat is, in principe, alles.

[Xilvo]

Ik heb het integratieprogramma wat verfijnd.
Massapunten vlak onder het oppervlak (minder dan de afstand tussen de massapunten) krijgen een weegfactor die kleiner wordt naarmate ze dichterbij het oppervlak liggen. Dat verlaagt de "ruis" in de uitkomst.
Het dichtheidsverloop binnen de aarde verloopt ongeveer zo:

[bron Wikipedia]
Dat heb ik benaderd door een lineair afnemende dichtheid van 12,5 in de kern tot 10 op 55% van de diameter, vanaf daar een lineaire afname van 4,8 tot 4 aan het oppervlak. Het gaat om relatieve veranderingen dus eenheden van de dichtheden zijn niet belangrijk, alleen de onderlinge verhoudingen.
Dat geeft deze grafiek. x-as is breedtegraad, in rood de berekende zwaartekrachtsversnelling door integratie over de ellipsoïde, in blauw daar de centrifugale versnelling bij opgeteld.

Code: Selecteer alles

NB    0 g_g  9.81000 g_tot  9.77612
NB    5 g_g  9.81017 g_tot  9.77655
NB   10 g_g  9.81048 g_tot  9.77762
NB   15 g_g  9.81103 g_tot  9.77942
NB   20 g_g  9.81179 g_tot  9.78187
NB   25 g_g  9.81277 g_tot  9.78494
NB   30 g_g  9.81385 g_tot  9.78844
NB   35 g_g  9.81510 g_tot  9.79237
NB   40 g_g  9.81638 g_tot  9.79650
NB   45 g_g  9.81776 g_tot  9.80083
NB   50 g_g  9.81910 g_tot  9.80510
NB   55 g_g  9.82046 g_tot  9.80931
NB   60 g_g  9.82168 g_tot  9.81321
NB   65 g_g  9.82283 g_tot  9.81678
NB   70 g_g  9.82378 g_tot  9.81981
NB   75 g_g  9.82458 g_tot  9.82231
NB   80 g_g  9.82519 g_tot  9.82417
NB   85 g_g  9.82558 g_tot  9.82532
NB   90 g_g  9.82565 g_tot  9.82565

Als ik het verschil in versnelling tussen 50 en 55 graden door 5 deel en vermenigvuldig met 2,303 (het breedtegraadverschil tussen Amsterdam en Chimay), dan geeft dat 0,00194 m/s².
De massa van 500 g zou dan in Chimay 0,099 g minder moeten wegen dan in Amsterdam. Dat ligt binnen de foutenmarge van de meting.

Correctie: Hierbij is het hoogteverschil tussen Amsterdam en Chimay nog niet meegenomen.
Dat geeft een extra verschil van 0,000722 m/s².
Het verschil in versnelling wordt dan 0,00266 m/s² en het verschil in gewicht 0,136 g. Net buiten de foutenmarge.

Voor de verhouding tussen equatoriale en polaire diameter is 1,003364 gebruikt.
Voor de integratie is een 'box' van 1600 x 1600 x 1600 massapunten gekozen. Dat zijn er ca 4 miljard waarvan ongeveer de helft binnen de ellipsoïde valt.

[jkien]

Je kunt misschien nog nagaan of g(φ) van je model voor Amsterdam en Chimay voldoende nauwkeurig overeenkomt met de reeds genoemde IGF formule, die g(φ) op de ellipsoide geeft. De nauwkeurigheid van de IGF formule is ±10⁻⁶ m/s², volgens wikipedia. Wikipedia geeft bovendien een verbeterde formule die een nauwkeurigheid van ±10⁻⁹ m/s² heeft.

[Xilvo]

Dat levert deze grafiek op:

De geelgroene lijn is de versnelling volgens de IGF-formule.
De overeenstemming is goed maar zeker niet perfect.
Voor het gewichtsverschil tussen Amsterdam en Chimay zal het weinig uitmaken.

Waar het verschil door komt durf ik niet te zeggen. Misschien zou ik de dichtheidsverschillen binnen de aarde nauwkeuriger moeten meerekenen.
Tot nu heb ik een aardig resultaat verkregen door het model met bestaande kennis te bouwen zonder naar het resultaat te kijken.
Ik wil de valkuil vermijden te gaan sleutelen om maar een beter overeenstemming met een model als IGF te krijgen.

Mijn berekende waarde heb ik wat verschoven om beter bij die van IFG aan te sluiten. Mijn model berekent tenslotte alleen relatieve versnellingen, dus de waarde voor de versnelling van de zwaartekracht op de evenaar kies ik zelf.

[HansH]

Waar het verschil door komt durf ik niet te zeggen. Misschien zou ik de dichtheidsverschillen binnen de aarde nauwkeuriger moeten meerekenen.
Tot nu heb ik een aardig resultaat verkregen door het model met bestaande kennis te bouwen zonder naar het resultaat te kijken.

Ik vraag me wel af hoe ze aan die dichtheids gegevens komen. en of die IGF formule daar dan op is gebaseerd of gebaseerd op veel metingen en daar simpelweg de constanten uit bepaald voor de IGF formule. al met al denk ik dat Xilvo met de gegeven uitleg behoorlijk veel info heeft meegenomen en er eigenlijk niet zo veel extra meer te halen is. ik kan het in mathcad stoppen maar denk niet dat we daar nog veel verder mee komen.
voor mij zijn het inmiddels heel veel berichten , maar ik kan niet goed de vinger erachter krijgen of het verschil tussen amsterdam en Chimay nu volledig begrepen is.

belangrijkste inzicht voor mij uit dit topic is dat zwaartekracht van de aarde veel meer is dan alleen het puntmassa model en een afstand door een vervormde bol vanwege de aardrotatie. door de bijdrage evenredig met 1/r^2 heeft massa dichtbij relatief veel meer invloed dus zal een grotere straal op de evenaar veel minder afname van g geven dan op basis van een puntmassa model, omdat dat oppervlak met bijbehorende massa in de nabije diepte in de buurt van de waarnemer steeds even dichtbij de waarnemer zit.

[HansH]

wat ik denk te concluderen uit alle berichten is dat de oorspronkelijke berekening van Jkien op basis van elliptische straal en puntmassa een variatie geeft tussen pool en evenaar van ca 0.04m/s2 (exclusief centripetale kracht) en de nauwkeuriger berekeningen op basis van zwaartekracht van alle punten van een ellipsoide een variatie geeft tussen pool en evenaar van ca 0.01m/s2 (exclusief centripetale kracht) de variaties tgv lokale invloeden (dichtheid, hoogteprofiel) heeft teruggerekend naar Zeenivo een variatie < ca 0.001m/s2. centripetale kracht is helder dus hoeft niet verder onderzocht te worden.
het verschil tussen puntmassa methode en ellips is 0.03m/s2 op de 9.81m/s2=0.3% dar is op 500g dan 1.52gram. maar dat zou dan tussen pool en evenaar zijn. dus moet nog verrekend worden tussen Amsterdam en Chimay (verschil =2.3 graden)
dat levert minder dan 0.001m/s2 (grof afgeschat, niet berekend) dat is dan ca 50mg opde 500g gewicht op de weegschaal. er is 110mg verschil gemeten. dus zou mijn conclusie zijn dat het verschil tussen puntmassa methode en ellips methode dat redelijk verklaart en dat lokale verschillen in de marge weg vallen.

[wnvl1]

in mathcad stelt deze berekening numeriek niet heel veel voor denk ik. verder optimaliseren door slim alleen het deel mee te nemen wat verschillend is tov een bol kan, maar kost waarschijnlijk meer tijd om het uit te werken dan om het met domme rekenkracht helemaal door te rekenen voor een complete ellipsoide. dus al ik tijd heb zal ik daar eens verder naar kijken.

Ik zal het bij gelegenheid eens proberen met een drievoudige integraal in scipy (python package voor numerieke berekeningen).
Je kan de coördinaten kiezen zoals in het bovenste antwoord hier.

https://math.stackexchange.com/question ... -integrals

Geen idee wat de numerieke fout zal zijn.

[Xilvo]

Moet lukken met scipy.integrate.tplquad.
Ik ben benieuwd hoe snel dat is.

[HansH]

Ik ben bezig om de integraal op te stellen om uit te rekenen in mathcad.
ben even simpel begonnen met het oppervlak van een cirkel via eerst alle stukje dy optellen dus .

2dintegraal 247 keer bekeken

dat is de integraal van y=-wortel(1-x^2) tot y=wortel(1-x^2) en dan dat weer integreren voor x=-1 tot x=+1
daar komt mooi pi uit zoals verwacht.

maar hoe stel je de grenzen op voor een 3d integraal op dezelfde manier voor de elipsoide?

3dintegraal 247 keer bekeken

je integreert dan bv eerst in x richting voor een gegeven y en z coordinaat. de x range is dan een funcie van y en z
dan integreer je die som nog een keer in y richting en is de y range dus een functie van x en z. dat levert het opervlak van een curkel in ht x-y vlak voor de gegeven z coordinaat. dan moet je als die cirkel oppervlakts nog integreren in z richting dus van z=-b tot z=b. maar hoe zien de integratiegrenzen voor elk stukje van de integraal er dan uit?