| 7 | | [[latex($\frac{\partial E}{\partial t} \color{red}{ - \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E}=\color{red}{\kappa_{0P} (4 \pi B-cE)} \color{green}{-\lambda \left(2\frac{\kappa_{0P}}{\kappa_{0R}}-1\right)\mathbf{v}\cdot \nabla E} \color{green}{-\nabla \cdot \left ( \frac{3-R_2}{2}\mathbf{v}E\right )}\color{blue}{+\frac{3-R_2}{2}\kappa_{0P}\frac{v^2}{c}E} $)]] |
| | 5 | [[latex($V' = V \left ( \frac{1}{1+\epsilon} \right ) + $)]] |
| | 6 | |
| | 7 | [[latex($M'R' = M'R + \sum{\Delta m_i \left ( r_i - R \right )}=MR + \sum{\Delta m_i r_i }$)]] |
| | 8 | |
| | 9 | [[latex($S'=S+\sum{\Delta m_i \left ( r_i -R \right ) \times v_i}$)]] |
| | 10 | |
| | 11 | However, a better treatment is to conserve angular momentum. |
| | 12 | |
| | 13 | [[latex($L=MR \times V$)]] |
| | 14 | |
| | 15 | [[latex($J=L+S$)]] |
| | 16 | |
| | 17 | [[latex($J'=J +\sum{\Delta m_i r_i \times v_i}$)]] |
| | 18 | |
| | 19 | [[latex($S'=S + MR \times V - M'R' \times V' + \sum{\Delta m_i r_i \times v_i}$)]] |
| | 20 | |
| | 21 | [[latex($S'=S + MR \times V - M'R' \times V' + \sum{\Delta m_i r_i \times v_i}$)]] |
| 20 | | [[latex($\frac{\partial \left (e+E\right) }{\partial t}=\color{red}{ \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E}$)]] |
| 21 | | |
| 22 | | which simplifies to |
| 23 | | |
| 24 | | [[latex($\left(1+\frac{\rho c_v}{4 E}\left(\frac{E}{a_R}\right )^{1/4}\right) \frac{\partial E}{\partial t}= \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E$)]] |
| 25 | | |
| 26 | | or |
| 27 | | |
| 28 | | [[latex($\left(1+\frac{e}{4E}\right) \frac{\partial E}{\partial t}= \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E$)]] |
| 29 | | |
| 30 | | |
| 31 | | The second term in parenthesis represents the extra 'inertia' the radiation field has due to its coupling with the gas. It is non-linear and this limits the time step that can be taken. |
| 32 | | |
| 33 | | [[latex($\Delta t \approx \frac{E}{\frac{\partial E}{\partial t}} = \frac{E+\frac{e}{4}}{\nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E}$)]] |
| 34 | | |
| 35 | | == Changes to the discretization == |
| 36 | | |
| 37 | | For the coupled system of equations we had the following: |
| 38 | | |
| 39 | | [[latex($\color{purple}{\left [ 1 + \psi \left( \alpha_{i+1/2} + \alpha_{i-1/2} + \frac{\epsilon_i}{ 1 +\psi \phi_i}\right ) \right ] E^{n+1}_i - \left ( \psi \alpha_{i+1/2} \right ) E^{n+1}_{i+1} - \left ( \psi \alpha_{i-1/2} \right ) E^{n+1}_{i-1} =\left [ 1 - \bar{\psi} \left( \alpha_{i+1/2} + \alpha_{i-1/2} +\frac{\epsilon_i }{ 1 +\psi \phi_i} \right ) \right ] E^n_i + \left ( \bar{\psi} \alpha_{i+1/2} \right ) E^{n}_{i+1} + \left ( \bar{\psi} \alpha_{i-1/2} \right ) E^{n}_{i-1} + \frac{\theta_i}{ 1 +\psi \phi_i}}$)]] |
| 40 | | |
| 41 | | If the gas and radiation are in thermal equilibrium, then we have [[latex($\theta_i = \epsilon_i E_i$)]] and we also have that in the limit that [[latex($\kappa_P \rightarrow \infty$)]], we have [[latex($\epsilon \rightarrow \infty$)]] and [[latex($\phi \rightarrow \infty $)]] |
| 42 | | |
| 43 | | This simplifies the above equation to |
| 44 | | |
| 45 | | [[latex($\color{purple}{\left [ \left(1 + \frac{\epsilon_i}{\phi_i} \right) + \psi \left( \alpha_{i+1/2} + \alpha_{i-1/2} \right ) \right ] E^{n+1}_i - \left ( \psi \alpha_{i+1/2} \right ) E^{n+1}_{i+1} - \left ( \psi \alpha_{i-1/2} \right ) E^{n+1}_{i-1} =\left [ \left(1 + \frac{\epsilon_i}{\phi_i}\right) - \bar{\psi} \left( \alpha_{i+1/2} + \alpha_{i-1/2} \right ) \right ] E^n_i + \left ( \bar{\psi} \alpha_{i+1/2} \right ) E^{n}_{i+1} + \left ( \bar{\psi} \alpha_{i-1/2} \right ) E^{n}_{i-1} }$)]] |
| 46 | | |
| 47 | | And [[latex($\frac{\epsilon_i}{\phi_i} = \frac{\frac{T_i}{4\Gamma}}{E_i} = \frac{\frac{T_i}{4\frac{\partial T}{\partial e}}}{E_i} $)]] which if we use our equation of state where [[latex($e \propto T$)]] gives [[latex($\frac{\epsilon_i}{\phi_i} = \frac{e_i}{4E_i}$)]] |
| 48 | | |
| 49 | | Now if we go back and calculate [[latex($\frac{\partial e}{\partial E} = \frac{\partial e}{\partial T}\frac{\partial T}{\partial E} = \frac{\frac{1}{\Gamma}}{\frac{4 E}{T}}$)]] we arrive at [[latex($\frac{\epsilon_i}{\phi_i}$)]] instead of [[latex($\frac{e_i}{4E_i}$)]] which is consistent with our derivation above. |
| 50 | | |
| 51 | | Also our time equation should be |
| 52 | | |
| 53 | | [[latex($\Delta t \approx \frac{E}{\frac{\partial E}{\partial t}} = \frac{E+\frac{T_i}{4\Gamma}}{\nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E}$)]] |
| 54 | | |
| 55 | | |
| 56 | | So in principle combining the gas and energy equations in the limit that of high planck opacity, does not change the matrix or rhs vector, however it does limit the ability for there to be strong source terms on the right in regions where the gas and radiation have gotten out of equilibrium. It is not clear how this effects the ability of the elliptic solver to converge to given tolerances. |
| 57 | | |
| 58 | | == Modifications to time steps == |
| 59 | | |
| 60 | | More importantly is the recognition of the time scales over which the internal energy can change. |
| 61 | | |
| 62 | | Previously we looked at the decoupled equation for the gas energy density |
| 63 | | |
| 64 | | [[latex($\frac{\partial e}{\partial t}=-\kappa_{0P}(4 \pi B-cE)$)]] |
| 65 | | |
| 66 | | [[latex($\Delta e = \Delta t \kappa_{0P} \left | 4 \pi B_0 -cE \right | < \xi \frac{T_0}{4 \Gamma}$)]] |
| 67 | | |
| 68 | | which gives [[latex($\Delta t < \xi \frac{T_0}{4 \Gamma \kappa_{0P} \left | 4 \pi B_0 - cE \right | }$)]] |
| 69 | | |
| 70 | | however, if the gas is in equilibrium with the radiation, this does not limit the time step at all - even though diffusion may quickly move the gas out of equilibrium with the radiation. |
| 71 | | |
| 72 | | We can account for this by expanding our equation |
| 73 | | |
| 74 | | [[latex($\Delta t < \xi \frac{T_0}{4 \Gamma \kappa_{0P} \left | 4 \pi B_0 - c\left(E_0+\partial_tE\Delta t\right ) \right | }$)]] which gives us a quadratic for [[latex($\Delta t$)]] |
| 75 | | |
| 76 | | [[latex($4 \Gamma \kappa_{0P} c \left | \partial_tE \right | \Delta t^2 + \left(4 \Gamma \kappa_{0P} \left | 4 \pi B_0 - cE_0 \right | \right ) \Delta t < \xi T_0 $)]] |
| 77 | | |
| 78 | | where we have conservatively assumed that the diffusion always causes the gas to be more out of equilibrium. |
| 79 | | |
| 80 | | Now if we recognize that there are three time scales at play here: |
| 81 | | |
| 82 | | || Diffusion Time || [[latex($\tau_D=\frac{T}{4 \Gamma \left | \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E \right | }$)]] || |
| 83 | | || Coupling Time || [[latex($\tau_C=\frac{T}{ 4 \Gamma \kappa_{0P} \left |4 \pi B_0-cE \right |}$)]] || |
| 84 | | || Absorption Time || [[latex($\tau_A=\frac{1}{c \kappa_{0P}}$)]] || |
| 85 | | |
| 86 | | then this quadratic simplifies to |
| 87 | | |
| 88 | | [[latex($\frac{\Delta t^2}{\tau_A\tau_D} + \frac{\Delta t}{\tau_C} < \xi$)]] |
| 89 | | |
| 90 | | [[latex($\Delta t = \sqrt{\left ( \frac{\tau_A \tau_D}{2\tau_C}\right)^2+ \xi{\tau_D\tau_A}}-\frac{\tau_A \tau_D}{2\tau_C} $)]] |
| 91 | | |
| 92 | | [[latex($\Delta t = \frac{\tau_A \tau_D}{2\tau_C} \left ( \sqrt{1+ \frac{4\xi \tau_C^2}{\tau_D\tau_A}} - 1 \right) $)]] |
| 93 | | |
| 94 | | |
| 95 | | So we can choose |
| 96 | | * the diffusion time if we assume the gas and radiation are strongly coupled - optically thick, |
| 97 | | * the coupling time if we assume that the gas is optically thin (which should imply the radiation is fairly diffused) |
| 98 | | * Or we can solve the quadratic |
| | 34 | which implies our method is only accurate if $ V << v_i$ and $\sum{\Delta m_i} << M$ |
