Het modelleren van de dichtheid van de Zon

Hoe modelleer je de dichtheid van de Zon?
De dichtheid van de Zon is natuurlijk allesbehalve homogeen, want het is evident dat de dichtheid in de kern van de Zon groter is dan in de buitenste lagen. Er zijn berekeningen waarbij die homogeniteit totaal niet relevant is, maar soms is het een doorslaggevende factor en dan moet die dichtheid op de een of andere manier wiskundig gemaakt worden.

Uit de combinatie van metingen en fysisch inzicht heeft men een zogenaamde best-fit-formula gemaakt, een vierdegraads polynoom, die te vinden is op de website van NASA:
Hierin is x de genormaliseerde straal, dus x = 0 is het middelpunt van de Zon en x = 1 is de rand van de Zon. Wanneer ik x = 0 invul dan krijg ik de dichtheid in het middelpunt van de Zon:
En wanneer ik x = 1 invul dan krijg ik de dichtheid aan de rand van de Zon:
Oeps, aan de rand van de Zon is de dichtheid niet nul, maar heeft zelfs een onfysische negatieve waarde. Ik ga de NASA-formule ook even integreren over het totale volume van de Zon, want dat leidt tot de massa van de Zon:
Omdat x genormaliseerd was moeten we dit nog vermenigvuldigen met R3 (R is de straal van de Zon) en dan vinden we m = 2.932 ∙ 1030 kg. Zo komen we bij de volgende oeps: deze waarde is bijna vijftig procent te hoog. Het NASA-polynoom is absoluut waardevol, maar er kleven een paar nadelen aan:
  1. de dichtheid aan de rand van de Zon is niet nul,
  2. de formule genereert ook onfysische negatieve waarden,
  3. het leidt tot een massa die bijna vijftig procent boven de werkelijke waarde ligt,
  4. de functie volgens vergelijking (1) is niet monotoon dalend hetgeen ook merkwaardig is (van het middelpunt van de Zon naar buiten toe neemt de dichtheid eerst af en later weer toe),
  5. de dichtheid in het middelpunt van de Zon heeft een verouderde waarde (155 in plaats van 162.2, zie de NASA Planetary Fact Sheet.
Iemand vertelde mij ooit: alle natuurlijke processen verlopen volgens e-machten. Dit bracht mij op het idee om de dichtheid van de Zon daarom volgens een e-macht te modelleren:
Hierin heb ik drie onbekende parameters, a, b en c, maar gelukkig heb ik ook drie randvoorwaarden waarmee ik die drie parameters kan bepalen. De eerste randvoorwaarde is dat voor r = 0 de dichtheid gelijk is aan de dichtheid in het middelpunt van de Zon:
De tweede randvoorwaarde is dat de dichtheid aan de rand van de Zon, waar r = R, gelijk is aan nul:
Uit de combinatie van deze twee randvoorwaarden kan ik c elimineren:
Maar ik had natuurlijk ook b kunnen elimineren:
En de derde randvoorwaarde is dat indien ik de dichtheid integreer over het totale volume van de Zon, dat dat de massa van de Zon als uitkomst moet opleveren. Ik maak hierbij gebruik van de tabel met integralen:
Hiervoor vond ik twee uitdrukkingen voor b en c als functie van a (de vergelijkingen (8) en (9)). Die ga ik nu gebruiken in vergelijking (10):
Voor de gemiddelde dichtheid geldt:
Daarmee wordt de vorige vergelijking:
Ik stel:

Waarmee vergelijking (13) uiteindelijk wordt:
Vervolgens roep ik de hulp in van Python om deze vergelijking op te lossen. Ik vind x = 8.7628642966 en daaruit volgen dan simpelweg de waarden van a (vergelijking (14)), b (vergelijking (8)) en c (vergelijking (6)). Tenslotte is het natuurlijk wel leuk om het NASA-polynoom en de e-macht even samen in een grafiek te zetten.

De grafiek van ρ (r) volgens het NASA-polynoom (de rode lijn) en volgens de e-macht (de groene lijn)
Dat ziet er toch best wel goed uit? En met de e-macht zijn de vijf nadelen van het NASA-polynoom, zoals ik die hiervoor opsomde, perfect opgeheven.

In één opzicht zijn beide functies echter heel onwerkelijk, want ze gaan niet naar een maximale waarde toe voor r = 0. We mogen ervan uitgaan dat de dichtheid in het middelpunt van de Zon het hoogst is en daarom verwachten we daar een horizontale raaklijn en zoals het hierbovenstaande plaatje laat zien is dat niet het geval. Ook daar is natuurlijk een oplossing voor te bedenken door in de exponent van de e-macht in plaats van r een r2 te plaatsen:
Vervolgens wordt het een herhaling van zetten om a, b en c te bepalen. De eerste randvoorwaarde is dat voor r = 0 de dichtheid gelijk is aan de dichtheid in het middelpunt van de Zon:
De tweede randvoorwaarde is dat de dichtheid aan de rand van de Zon, waar r = R, gelijk is aan nul:
Uit de combinatie van deze twee randvoorwaarden kan ik c elimineren:
Maar ik had natuurlijk ook b kunnen elimineren:
En de derde randvoorwaarde is dat indien ik de dichtheid integreer over het totale volume van de Zon, dat dat de massa van de Zon als uitkomst moet opleveren. Ik maak hierbij wederom gebruik van de tabel met integralen:
Hiervoor vond ik twee uitdrukkingen voor b en c als functie van a (de vergelijkingen (20) en (21)). Die ga ik nu gebruiken in vergelijking (22):
Ik stel:
Waarmee ik uiteindelijk kom tot:
Vervolgens is Python weer aan de beurt om deze vergelijking op te lossen. Ik vind x = 28.5928009589 en daaruit volgen dan weer de waarden van a (vergelijking (24)), b (vergelijking (20)) en c (vergelijking (18)). Ter vergelijking voeg ik deze e-macht toe aan de vorige grafiek.

De grafiek van ρ (r) volgens het NASA-polynoom (de rode lijn), volgens e−ar (de groene lijn)
en volgens e−ar2 (de blauwe lijn)
Om redenen die later duidelijk zullen worden doorloop ik dit hele spektakel nog een keer, maar nu met twee e-machten op elkaar gestapeld als volgt (en ik introduceer twee nieuwe parameters: p en q):
Het wordt nogmaals een herhaling van zetten om a, b en c te bepalen. De eerste randvoorwaarde is dat voor r = 0 de dichtheid gelijk is aan de dichtheid in het middelpunt van de Zon:
De tweede randvoorwaarde is dat de dichtheid aan de rand van de Zon, waar r = R, gelijk is aan nul:
Uit de combinatie van deze twee randvoorwaarden kan ik c elimineren:
Maar ik had natuurlijk ook b weer kunnen elimineren:
En de derde randvoorwaarde is dat indien ik de dichtheid integreer over het totale volume van de Zon, dat dat de massa van de Zon als uitkomst moet opleveren. Ik maak hierbij wederom gebruik van de tabel met integralen:
Hiervoor vond ik twee uitdrukkingen voor b en c als functie van a, p en q (de vergelijkingen (29) en (30)). Die ga ik nu gebruiken in vergelijking (31):
Ik stel:

Waarmee ik uiteindelijk kom tot:
Om goede redenen die later duidelijk zullen worden kies ik d = 0.228955 en p = 0.75, en ik roep nu voor de laatste keer Python te hulp om de bovenstaande vergelijking op te lossen. Ik vind x = 59.8961489749 en daaruit volgen dan weer de waarden van a (vergelijking (33)), b (vergelijking (29)) en c (vergelijking (27)). Ter vergelijking voeg ik deze e-macht ook toe aan de vorige grafiek.

De grafiek van ρ (r) volgens het NASA-polynoom (de rode lijn), volgens e−ar (de groene lijn)
volgens e−ar2 (de blauwe lijn) en volgens e−ar2 + e−qr2 (de oranje lijn)
In het volgende vraagstuk bereken ik het traagheidsmoment van de Zon met al deze verschillende modelleringen.