Changes between Version 14 and Version 15 of u/johannjc/scratchpad
- Timestamp:
- 01/10/14 14:46:10 (11 years ago)
Legend:
- Unmodified
- Added
- Removed
- Modified
-
u/johannjc/scratchpad
v14 v15 1 == Self-gravitating equilibrium profile using density profile and external pressure == 1 [[latex($\rho c_v \frac{\partial T}{\partial t} = \nabla \left ( \kappa T^{n} \nabla T \right )$)]] 2 2 3 Mathematically, we are solving the equation of hydrostatic equilibrium 3 [[latex($\rho c_v \frac{\partial T}{\partial t} =\nabla^2 \left ( \frac{1}{n+1} \kappa T^{n+1} \right )$)]] 4 4 5 [[latex($\nabla P(r) = -\rho(r) \nabla \phi(r)$)]] 6 7 subject to the constraint that [[latex($P(R)=P_0$)]]. 8 9 Since [[latex($\phi(r) = -\frac{G M(r)}{r}$)]] we can rewrite the equation for HSE as 10 {{{ 11 #!latex 12 \begin{equation} 13 \frac{dP}{dr}=-\frac{GM(r)\rho(r)}{r^2} 14 \end{equation} 15 }}} 5 [[latex($\rho c_v \frac{\partial T}{\partial t} =\nabla^2 \left ( \frac{1}{n+1} \kappa T^{n} T \right )$)]] 16 6 17 7 18 and we can express the enclosed mass as 8 [[latex($\rho_i c_v \left ( T_i^{j+1}-T_i^j \right ) = \frac{\kappa \Delta t}{\Delta x^2}\frac{1}{n+1} \left (T_{i+1}^{j^n} T_{i+1}^{j+1} - 2 T_{i}^{j^n} T_{i}^{j+1} + T_{i-1}^{j^n} T_{i-1}^{j+1} \right )$)]] 19 9 20 {{{ 21 #!latex 22 \setcounter{equation}{1} 23 \begin{equation} 24 M(r)=4\pi\displaystyle\int_0^r{\rho(r') r'^2dr'} 25 \end{equation} 26 }}} 10 [[latex($\alpha = \frac{\kappa \Delta t}{\Delta x^2 \left ( n+1 \right ) }$)]] 27 11 28 Now we can discretize equation 2 as 29 30 {{{ 31 #!latex 32 \setcounter{equation}{2} 33 \begin{equation} 34 M_{i+1}=M_i + 4\pi \rho_i r_i^2\Delta r_i 35 \end{equation} 36 }}} 37 38 where [[latex($\Delta r_j = r_{j+1}-r_j$)]] 39 40 and equation 1 as 41 42 {{{ 43 #!latex 44 \setcounter{equation}{3} 45 \begin{equation} 46 P_{i+1}=P_i-\Delta r_i \frac{G M_i \rho_i}{r_i^2} 47 \end{equation} 48 }}} 12 [[latex($ \left (\rho_i c_v + 2 \alpha T_i^{j^n} \right ) T_i^{j+1} - \alpha T_{i+1}^{j^n} T_{i+1}^{j+1} - \alpha T_{i-1}^{j^n} T_{i-1}^{j+1} = \rho_i c_v T_i^j$)]] 49 13 50 14 51 but a taylor series expansion shows that the error at each step in the integration goes as 52 53 [[latex($\frac{\Delta r_j}{r_j} + \frac{\Delta \rho_j}{\rho_j} + ...$)]] 54 55 so we need to limit out point spacing - or improve the accuracy of our discretization. 56 57 If we apply linear interpolation to our density profile, we can discretize equation 2 as 58 59 {{{ 60 #!latex 61 \setcounter{equation}{5} 62 \begin{equation} 63 M_i=M(r_i)=4\pi\displaystyle\sum_{j=0}^i{\rho_j r_j^2\Delta r_j \left ( 1 + \frac{\Delta r_j}{r_j} + \frac{\Delta \rho}{2\rho} + \frac{\Delta r_j^2}{3 r_j^2} + \frac{2\Delta \rho \Delta r_j}{3 \rho r_j} + \frac{\Delta \rho \Delta r_j^2}{4 \rho r_j^2}\right )} 64 \end{equation} 65 }}} 66 67 and equation 4 becomes a bit of a monster - and we have to worry about limits as r -> 0. 68 {{{ 69 if (r0 == 0d0) then 70 dp=drho*((pi*drho*dr**4)/4d0 + (4d0*pi*rho*dr**3)/9d0) + (rho*(pi*drho*dr**3 + 2d0*pi*rho*dr**2))/3d0 71 else 72 epsilon=log(r0/(r0+dr)) 73 dp=((-36d0*drho*r0)*epsilon*M0*dr + (-36d0*drho*r0**2)*epsilon*M0 + (- 12d0*pi*drho**2*r0**5 + 48d0*pi*rho*drho*r0**4)*epsilon*dr + (-12d0*pi*drho**2*r0**6 + 48d0*pi*rho*drho*r0**5)*epsilon + (36d0*rho - 36d0*drho*r0)*M0*dr + (9d0*pi*drho**2*r0)*dr**5 + (17d0*pi*drho**2*r0**2 + 28d0*pi*rho*drho*r0)*dr**4 + (2d0*pi*drho**2*r0**3 + 64d0*pi*drho*r0**2*rho + 24d0*pi*r0*rho**2)*dr**3 + (- 6d0*pi*drho**2*r0**4 + 24d0*pi*drho*r0**3*rho + 72d0*pi*r0**2*rho**2)*dr**2 +(- 12d0*pi*drho**2*r0**5 + 48d0*pi*rho*drho*r0**4)*dr)/(36d0*r0**2 + 36d0*dr*r0) 74 end if 75 data(i,i_Press)=data(i+1,i_Press)+ScaleGrav*dp 76 }}} 15 Now to make the scheme 2nd order in time we can replace [[latex($T_i^{j+1}$)]] with [[latex($\frac{T_i^{j}+T_i^{j+1}}{2}$)]] 77 16 78 17 18 [[latex($ \left (\rho_i c_v + 2 \beta T_i^{j^n} \right ) T_i^{j+1} - \beta T_{i+1}^{j^n} T_{i+1}^{j+1} - \beta T_{i-1}^{j^n} T_{i-1}^{j+1} = \rho_i c_v T_i^j - 2 \beta T_i^{j^n} T_i^{j} + \beta T_{i+1}^{j^n} T_{i+1}^{j} + \beta T_{i-1}^{j^n} T_{i-1}^{j} $)]] 19 20 21 Which simplifies to: 22 23 [[latex($\left (\rho_i c_v + 2 \beta T_i^{j^n} \right ) T_i^{j+1} - \beta T_{i+1}^{j^n} T_{i+1}^{j+1} - \beta T_{i-1}^{j^n} T_{i-1}^{j+1} = \rho_i c_v T_i^j - 2 \beta T_i^{j^{n+1}} + \beta T_{i+1}^{j^{n+1}} + \beta T_{i-1}^{j^{n+1}}$)]] 24 25