A&A Wykład 3

Struktura i własności sferycznych ciał niebieskich

Dynamic[{DateString[],Clock[{1,1},1]}//First]

Andrzej Odrzywołek

Zakład Teorii Względności i Astrofizyki, Instytut Fizyki Teoretycznej UJ
2025-03-25

Analiza struktury wewnętrznej dla logotropowego równania stanu

Zakładamy, że sferyczna ,,Ziemia’’ (bardziej asteroida) zbudowana jest z ,,metalu’’ o stałym wspólczynniku sprężystości objętościowej (ang. bulk modulus) K i gęstości nieskompresowanej (Subscript[ρ, 0]):

bulkModulus=K==ρ D[P[ρ],ρ]

rozw=DSolve[{bulkModulus,P[ρ0]==0},P,ρ]

P[ρ0]/.rozw[[1]]

eos=FullSimplify[P[ρ]/.rozw[[1]],ρ0>0]

P=eos/.ρ->ρ[r] (* Elastic EOS *)

eq1=1/ρ[r]D[P,r]==-G m[r]/r^2(* Hydrostatic equilibrium *)

eq1a=D[eq1[[1]]*r^2/G,r]+m'[r]/.m'[r]->4 Pi r^2 ρ[r] (* Use continuity equation *)

eq2=eq1a*G*ρ[r]^2/K/r^2//Expand

subs={ρ->Function[r,ρC u[r/λ]]}

eq2a=eq2/.subs/.r->λ x

eq3=eq2a*λ^2/ρC//Expand

eqλ= eq3/.{u->(1&)}

Λ=λ/.Solve[eqλ==1,λ][[-1]]

Λ^2

Λ/.{K->Quantity[170,"GPa"],ρC->Quantity[7874,"kg/m^3"],G->Quantity[1,"GravitationalConstant"]}

UnitConvert[%,"km"]

eq4=eq3/.λ->Λ

eq5=eq4/.u->Function[x,1/w[x]]

eq6=-eq5*w[x]^2//Expand

Funkcja logotropowa U[x]

(* Rozwiązanie numeryczne na funkcję U *)

ϵ=Sqrt[$MachineEpsilon]

rozw=NDSolve[{u''[x]+2/x u'[x]+u[x]^3-2 u'[x]^2/u[x]==0,u[ϵ]==1,u'[ϵ]==0},u,{x,ϵ/1024,1024*2},WorkingPrecision->32][[1]]

U=u/.rozw;

Export["U.wl",U]

SetDirectory[NotebookDirectory[]];

U = Import["http://th.if.uj.edu.pl/~odrzywolek/homepage/students/aa/U.wl"]

(*
ϵ=1/128;rozw=NDSolve[{u''[x]+2/x u'[x]+u[x]^3-2 u'[x]^2/u[x]==0,u[ϵ]==1-ϵ^2/6+13/360ϵ^4-25/3024 ϵ^6,u'[ϵ]==-(ϵ/3)+((13 (ϵ^3))/90)-((25 (ϵ^5))/504)},u,{x,ϵ/1024,1024},WorkingPrecision->64,MaxSteps->10^5][[1]];
Export["U.wl",u/.rozw];
U=Import["U.wl"]
*)

w=Plot[U[x],{x,0,12},PlotRange->{0,1.1},PlotStyle->{Black,AbsoluteThickness[1.5]},
AxesLabel->{"x","U(x)"},
Ticks->{{0,1,{InverseFunction[U][1/GoldenRatio],"(Subscript[x, t])"},3,4,5,10},{0,0.25,0.5,{1/GoldenRatio,"1/(Subscript[η, c])"},0.75,1}},
Prolog->
{Dashed,AbsoluteThickness[0.5],
Line[{{0,1/GoldenRatio},{InverseFunction[U][1/GoldenRatio],1/GoldenRatio}}],
Line[{{InverseFunction[U][1/GoldenRatio],1/GoldenRatio},{InverseFunction[U][1/GoldenRatio],0}}],
Point[{InverseFunction[U][1/GoldenRatio],1/GoldenRatio}]
}
]

inset=LogLogPlot[U[x]-Sqrt[2]/x//Abs,{x,1,2048},PlotRange->Automatic,PlotStyle->{Black,AbsoluteThickness[0.5]},
AxesLabel->{"x","|U(x)-(Sqrt[2])/x|"}
]

(* Rozwinięcie w szereg potęgowy *)

num=3 (* Zwiększyć, aby dostać rozwinięcie wyższego rzędu *)

szereg=Sum[a[i]x^i,{i,0,num}]

f=Function[x,Evaluate@szereg]

a[0]=1

a[1]=0

(u''[x]+2/x u'[x]+u[x]^3-2 u'[x]^2/u[x])/.u->f

rozwiniecie=Series[%,{x,0,num}]//Normal

c=CoefficientList[rozwiniecie,x]

sys=Table[c[[k]]==0,{k,1,num-1}]

niewiadome=CoefficientList[szereg,x][[3;;]]

rozw=Solve[sys,niewiadome]

approx=f[x]/.rozw[[1]]

Masa M, promień R i krzywa M(R)

wMR=ParametricPlot[{{x U[x],-3x^2 U'[x]},{x,x^3}},{x,0.001,1024},PlotStyle->{Black,AbsoluteThickness[1.5]},
PlotRange->{0,4.5},PlotPoints->164,GridLines->{{Sqrt[2]},{3Sqrt[2],9/2}},AxesLabel->{"R/(Subscript[R, 0])","M/(Subscript[M, 0])"},AspectRatio->1/GoldenRatio]

Integrate[D[K Log[ρ/ρ0],ρ]/ρ,ρ]

Kulistość Ziemi i innych ciał niebieskich

Wyobraźmy sobie gigantyczny wieżowiec o wysokości H postawiony na sferycznym ciele o promieniu R. Niech wszystko będzie zbudowane np. z żelaza o gęstości (Subscript[ρ, 0]). Jaka musi być wysokość H, aby wieżowiec rozpłynął się pod własnym ciężarem? Niech naprężenie niszczące (wytrzymałość na ściskanie) wynosi σ.

H ρ0 g == σ/.g->G M/R^2

%/.M->4/3Pi R^3 ρ0

%/.H->R

R^2/.Solve[%,R]

Analiza wymiarowa

KK=kg*m/s^2/m^2

ρ0=kg/m^3

G=m^3/s^2/kg

KK^α*ρ0^β*G^γ==m//PowerExpand (* Please check procedure for kg (mass scale) instead of m (length scale) *)

Log/@%

PowerExpand[%]

%/.{Log[kg]->lgkg,Log[m]->lgm,Log[s]->lgs}

rozw=SolveAlways[%,{lgm,lgs,lgkg}]

Clear[KK,ρ0,G]

KK^α*ρ0^β*G^γ/.rozw