Changes between Version 173 and Version 174 of FluxLimitedDiffusion
- Timestamp:
- 04/03/13 13:49:21 (12 years ago)
Legend:
- Unmodified
- Added
- Removed
- Modified
-
FluxLimitedDiffusion
v173 v174 1 1 [[PageOutline()]] 2 2 3 test2 $$ 4 \pi $$4 5 test 3 \( 4\pi \)6 7 [[latex(4 \pi)]]8 3 = Physics of Radiation Transfer = 9 4 … … 15 10 * Instead of storing the phase space density of photons, the spectral intensity is the phase space density of energy flux... 16 11 17 Going between photon number and energy just involves a factor of [[latex(h \nu)]] and going from energy density to energy flux density just involves a factor of [[latex(c)]]so we have18 19 [[latex(I \left ( \nu, \mathbf{x}, \Omega, \right ) = h \nu c f \left ( \nu, \mathbf{x}, \Omega, \right ) )]]12 Going between photon number and energy just involves a factor of \( h \nu \) and going from energy density to energy flux density just involves a factor of \(c\) so we have 13 14 $$ I \left ( \nu, \mathbf{x}, \Omega, \right ) = h \nu c f \left ( \nu, \mathbf{x}, \Omega, \right ) $$ 20 15 21 16 This can also be seen by considering the differential energy … … 38 33 == Deriving the Transport Equation == 39 34 40 If we consider the Boltzmann transport equation for photons of a specific frequency [[latex(f_\nu)]]we have35 If we consider the Boltzmann transport equation for photons of a specific frequency \(f_\nu\) we have 41 36 42 37 [[latex(\frac{\partial}{\partial t} f_\nu + \mathbf{v} \cdot \nabla f + \mathbf{F} \cdot \nabla_p f = \left ( \frac{\partial f}{\partial t} \right )_{coll} )]] … … 46 41 [[latex(\frac{\partial}{\partial t} f_\nu + c \mathbf{n} \cdot \nabla f_\nu = A_{\nu} - \chi_\nu f_\nu c )]] 47 42 48 where [[latex(A_{\nu} )]] is the emission rate of photons of frequency [[latex(\nu)]] and the mean free path length is given by [[latex(\chi_\nu = \sigma_\nu n)]] where [[latex(\sigma_\nu)]] is the particle scattering cross section and [[latex(n)]]is the number density of particles.49 50 Now if we multiply through by [[latex( h\nu)]]43 where \(A_{\nu} \) is the emission rate of photons of frequency \(\nu\) and the mean free path length is given by \(\chi_\nu = \sigma_\nu n\) where \(\sigma_\nu\) is the particle scattering cross section and \(n\) is the number density of particles. 44 45 Now if we multiply through by \( h\nu\) 51 46 52 47 we have … … 55 50 56 51 where 57 [[latex(\eta_\nu = h\nu A_{\nu})]]is the radiative power52 \(\eta_\nu = h\nu A_{\nu}\) is the radiative power 58 53 59 54 If we solve the transport equation along a characteristic … … 65 60 [[latex(\frac{d I_\nu}{d s} = \frac{\partial I_\nu}{\partial x^i} \frac{\partial{x^i}}{\partial s} + \frac{\partial I_\nu}{\partial t}\frac{\partial t}{\partial s} = \mathbf{n} \cdot \nabla I_\nu + \frac{1}{c}\frac{\partial I_\nu}{\partial t} = \eta_\nu(s) - \chi_\nu(s) I_\nu (s) )]] 66 61 67 where [[latex(f(s) = f(\mathbf{x}(s), t(s)) = f \left ( \mathbf{x_0}+\mathbf{n} s, \frac{s}{c} \right ) )]]68 69 and then we can divide through by [[latex(\chi_\nu(s))]]we get62 where \(f(s) = f(\mathbf{x}(s), t(s)) = f \left ( \mathbf{x_0}+\mathbf{n} s, \frac{s}{c} \right ) \) 63 64 and then we can divide through by \(\chi_\nu(s)\) we get 70 65 71 66 [[latex(\frac{dI_\nu}{\chi_\nu(s) ds} = \frac{\eta_\nu(s)}{\chi_\nu(s)} - I_\nu(s) = S_\nu(s) - I_\nu(s))]] … … 73 68 Now if we define 74 69 75 [[latex(d\tau_\nu = \chi_\nu(s) ds )]]70 \(d\tau_\nu = \chi_\nu(s) ds \) 76 71 which gives 77 72 … … 88 83 [[latex( f \left ( \tau_\nu \right ) = f \left ( s \left ( \tau_\nu \right ) \right ) = f \left ( \mathbf{x} \left (s \left ( \tau_\nu \right ) \right ), t \left ( s \left ( \tau_\nu \right ) \right ) \right ) )]] 89 84 90 Also if we include scattering then the source function can depend on the mean radiative flux [[latex(\frac{cE}{4 \pi})]]and the transport equation becomes an integro-differential equation that must be solved iteratively...85 Also if we include scattering then the source function can depend on the mean radiative flux \(\frac{cE}{4 \pi}\) and the transport equation becomes an integro-differential equation that must be solved iteratively... 91 86 92 87 There are also a few important dimensionless numbers to consider: … … 121 116 || [[latex(c\mathbf{G} = \int_0^\infty d \nu \int d \Omega \left [ \kappa (\mathbf{n}, v)I(\mathbf{n}, \nu)-\eta(\mathbf{n},v) \right ] \mathbf{n})]] || 122 117 123 If we had a closure relation for the radiation pressure then we could solve this system. For gas particles, collisions tend to produce a Boltzmann Distribution which is isotropic and gives a pressure tensor that is a multiple of the identity tensor. Photons do not "collide" with each other and they all have the same velocity 'c' but in various directions. If the field were isotropic than [[latex(P^{ij}=\delta^{ij} 1/3 E)]] but in general [[latex(P^{ij}=f^{ij} E)]]where 'f' is the Eddington Tensor.118 If we had a closure relation for the radiation pressure then we could solve this system. For gas particles, collisions tend to produce a Boltzmann Distribution which is isotropic and gives a pressure tensor that is a multiple of the identity tensor. Photons do not "collide" with each other and they all have the same velocity 'c' but in various directions. If the field were isotropic than \(P^{ij}=\delta^{ij} 1/3 E\) but in general \(P^{ij}=f^{ij} E\) where 'f' is the Eddington Tensor. 124 119 125 120 [[CollapsibleEnd()]] … … 136 131 || [[latex(\kappa_{0F})]] || comoving-frame radiation flux weighted opacity || 137 132 138 * If the radiation has a blackbody spectrum then [[latex(\kappa_{0E}=\kappa_{0P})]]133 * If the radiation has a blackbody spectrum then \(\kappa_{0E}=\kappa_{0P}\) 139 134 * If the radiation is optically thick, then 140 135 [[latex(\mathbf{F_0}(\nu_0) \propto -\nabla E_0(\nu_0)/\kappa_0(\nu_0) \propto -[\partial B(\nu_0, T_0) / \partial T_0](\nabla T_0)/\kappa_0(\nu_0))]] 141 136 which implies that [[latex(\kappa_{0F}^{-1}=\kappa_{0R}^{-1}=\frac{\int_0^\infty d \nu_0 \kappa_0(\nu_0)^{-1}[\partial B(\nu_0,T_0)/\partial T_0]}{\int_0^\infty d \nu_0 [\partial B(\nu_0, T_0)/\partial T_0]})]] 142 * In the optically thin regime, [[latex(|\mathbf{F}_0(\nu_0)| \rightarrow cE_0(\nu_0))]] so we would have [[latex(\kappa_{0F}=\kappa_{0E})]] however assuming a blackbody temperature in the optically thin limit may be any more accurate than assuming that [[latex(\kappa_{0F}=\kappa_{0R})]]137 * In the optically thin regime, \(|\mathbf{F}_0(\nu_0)| \rightarrow cE_0(\nu_0)\) so we would have \(\kappa_{0F}=\kappa_{0E}\) however assuming a blackbody temperature in the optically thin limit may be any more accurate than assuming that \(\kappa_{0F}=\kappa_{0R}\) 143 138 [[CollapsibleEnd()]] 144 139 … … 147 142 148 143 The flux limited diffusion approximation drops the radiation momentum equation in favor of 149 [[latex(\mathbf{F}_0=-\frac{c\lambda}{\kappa_{0R}}\nabla E_0)]] 150 151 where [[latex(\lambda)]]is the flux-limiter144 \(\mathbf{F}_0=-\frac{c\lambda}{\kappa_{0R}}\nabla E_0\) 145 146 where \(\lambda\) is the flux-limiter 152 147 153 148 [[latex(\lambda = \frac{1}{R} \left ( \coth R - \frac{1}{R} \right ) )]] … … 252 247 == Implicit Update 1 == 253 248 254 For now we will assume that [[latex(\kappa_{0P})]] and [[latex(\kappa_{0R})]]are constant over the implicit update and we will treat the energy as the total internal energy ignoring kinetic and magnetic contributions. In this case we can solve the radiation energy equations:249 For now we will assume that \(\kappa_{0P}\) and \(\kappa_{0R}\) are constant over the implicit update and we will treat the energy as the total internal energy ignoring kinetic and magnetic contributions. In this case we can solve the radiation energy equations: 255 250 256 251 [[latex(\frac{\partial E}{\partial t} = \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E + \kappa_{0P} (4 \pi B(T)-cE) = \nabla \cdot \mathbf{F} + \kappa_{0P} (4 \pi B(T)-cE))]] … … 261 256 === Expanding about e,,0,, === 262 257 263 Of course even if the opacity is independent of energy and radiation energy, the above combined system of equations is still non-linear due to the dependence on Temperature of the Planck Function [[latex(B(T))]]258 Of course even if the opacity is independent of energy and radiation energy, the above combined system of equations is still non-linear due to the dependence on Temperature of the Planck Function \(B(T)\) 264 259 265 260 If we ignore the changes in the Temperature due to heating during the implicit step which would feed back into the source function. We can improve this by writing … … 276 271 [[latex(\frac{\partial e}{\partial t} = - \kappa_{0P} \left [ 4 \pi B_0 \left ( 1 + 4\Gamma \frac{e-e_0}{T_0} \right )-cE \right ] )]] 277 272 278 which will be accurate as long as [[latex(4\Gamma \frac{e-e_0}{T_0} < \xi << 1)]] or [[latex(\Delta e = e-e_0 < \xi \frac{T_0}{4 \Gamma})]]273 which will be accurate as long as \(4\Gamma \frac{e-e_0}{T_0} < \xi << 1\) or \(\Delta e = e-e_0 < \xi \frac{T_0}{4 \Gamma}\) 279 274 280 275 We can calculate the time scale for this to be true using the evolution equation for the energy density … … 324 319 or we can parameterize the solution 325 320 [[latex(E^*_i = \psi E^{n+1}_i + \bar{\psi}E^n_i)]] 326 where [[latex(\bar{\psi} = 1-\psi)]]327 328 Backward Euler has [[latex(\psi=1)]] and Crank Nicholson has [[latex(\psi=1/2)]]329 330 Forward Euler has [[latex(\psi=0)]]321 where \(\bar{\psi} = 1-\psi\) 322 323 Backward Euler has \(\psi=1\) and Crank Nicholson has \(\psi=1/2\) 324 325 Forward Euler has \(\psi=0\) 331 326 332 327 In any event in 1D we have the following matrix coefficients … … 345 340 346 341 347 If we ignore the change in the Planck function due to heating during the implicit solve, it is equivalent to setting [[latex(\psi \phi = 0)]]This gives the following equations:342 If we ignore the change in the Planck function due to heating during the implicit solve, it is equivalent to setting \(\psi \phi = 0\) This gives the following equations: 348 343 349 344 [[latex(\left [ 1 + \psi \left( \alpha^n_{i+1/2} + \alpha^n_{i-1/2} + \epsilon^n_i \right ) \right ] E^{n+1}_i - \left ( \psi \alpha^n_{i+1/2} \right ) E^{n+1}_{i+1} - \left ( \psi \alpha^n_{i-1/2} \right ) E^{n+1}_{i-1} =\left [ 1 - \bar{\psi} \left( \alpha^n_{i+1/2} + \alpha^n_{i-1/2} + \epsilon^n_i \right ) \right ] E^n_i+ \left ( \bar{\psi} \alpha^n_{i+1/2} \right ) E^{n}_{i+1} + \left ( \bar{\psi} \alpha^n_{i-1/2} \right ) E^{n}_{i-1} +\phi^n_i e^n_i + \theta^n_i)]] … … 376 371 [[latex(\frac{c \lambda}{\kappa_{0R}} \nabla E = \mathbf{F} = cE\mathbf{n} = \frac{c \lambda_g}{\kappa_g}\frac{ \left ( E-E_g \right ) }{\Delta x})]] 377 372 378 Clearly if we set [[latex(E_g = 0)]] and [[latex(\lambda_g =\kappa_g \Delta x)]]we should get the correct flux.373 Clearly if we set \(E_g = 0\) and \(\lambda_g =\kappa_g \Delta x\) we should get the correct flux. 379 374 380 375 This corresponds to an 381 [[latex(\alpha = c \frac{\Delta t}{\Delta x})]] 382 383 So we would just modify [[latex(\alpha)]]and zero out the matrix coefficient to the ghost zone376 \(\alpha = c \frac{\Delta t}{\Delta x}\) 377 378 So we would just modify \(\alpha\) and zero out the matrix coefficient to the ghost zone 384 379 385 380 === Constant Slope Boundary === 386 Here we want the flux to be constant so energy does not pile up near the boundary. If we cancel all derivative terms on both sides of the cell, this will effectively match the incoming flux with the outgoing flux. This can also be done by setting [[latex(\alpha_g = \alpha_i = 0)]]381 Here we want the flux to be constant so energy does not pile up near the boundary. If we cancel all derivative terms on both sides of the cell, this will effectively match the incoming flux with the outgoing flux. This can also be done by setting \(\alpha_g = \alpha_i = 0\) 387 382 388 383 === Periodic Boundary === … … 394 389 395 390 === Reflecting/ZeroSlope Boundary === 396 Reflecting boundary should be fairly straightforward. This an be achieved by setting [[latex(\alpha_g = 0)]] which zeros out any flux - and has the same effect as setting [[latex(E^{*}_g=E^{*}_i)]] or [[latex(E^{n+1}_g=E^{n+1}_i \mbox{ and } E^{n}_g=E^{n}_i)]]391 Reflecting boundary should be fairly straightforward. This an be achieved by setting \(\alpha_g = 0\) which zeros out any flux - and has the same effect as setting \(E^{*}_g=E^{*}_i\) or \(E^{n+1}_g=E^{n+1}_i \mbox{ and } E^{n}_g=E^{n}_i\) 397 392 398 393 === Constant radiative flux === 399 To have a constant radiative flux we must zero out terms involving the gradient and just add [[latex(F_0 \frac{\Delta t}{\Delta x})]]in the source vector394 To have a constant radiative flux we must zero out terms involving the gradient and just add \(F_0 \frac{\Delta t}{\Delta x}\) in the source vector 400 395 401 396 === Summary === … … 473 468 == Implicit Update 2 == 474 469 475 For now we will assume that [[latex(\kappa_{0P})]] and [[latex(\kappa_{0R})]]are constant over the implicit update and we will treat the energy as the total internal energy ignoring kinetic and magnetic contributions. In this case we can solve the radiation energy equations:470 For now we will assume that \(\kappa_{0P}\) and \(\kappa_{0R}\) are constant over the implicit update and we will treat the energy as the total internal energy ignoring kinetic and magnetic contributions. In this case we can solve the radiation energy equations: 476 471 477 472 [[latex(\frac{\partial e}{\partial t} = -\kappa_{0P}(4 \pi B-cE) + \lambda \left ( 2 \frac{\kappa_{0P}}{\kappa_{0R}}-1 \right ) \mathbf{v} \cdot \nabla E -\frac{3-R_2}{2}\kappa_{0P}\frac{v^2}{c}E )]] … … 536 531 [[latex(\frac{\partial e}{\partial t} = - \kappa_{0P} \left [ 4 \pi B_0 \left ( 1 + 4\Gamma \frac{e-e_0}{T_0} \right )-cE \right ] + \lambda \left ( 2 \frac{\kappa_{0P}}{\kappa_{0R}}-1 \right ) \mathbf{v} \cdot \nabla E )]] 537 532 538 which will be accurate as long as [[latex(4\Gamma \frac{e-e_0}{T_0} < \xi << 1)]] or [[latex(\Delta e = e-e_0 < \xi \frac{T_0}{4 \Gamma})]]533 which will be accurate as long as \(4\Gamma \frac{e-e_0}{T_0} < \xi << 1\) or \(\Delta e = e-e_0 < \xi \frac{T_0}{4 \Gamma}\) 539 534 540 535 We can calculate the time scale for this to be true using the evolution equation for the energy density … … 650 645 651 646 == Physical Boundary Conditions == 652 For each boundary type we need to specify [[latex(E^n_g, E^{n+1}_g, \kappa_g, \mbox { and } v_g)]]in terms of other quantities (ie including E^n+1^_i). The v,,g,, and kappa,,g,, terms should come from the hydro boundary conditions so we just need an equation for E^n^,,g,,647 For each boundary type we need to specify \(E^n_g, E^{n+1}_g, \kappa_g, \mbox { and } v_g\) in terms of other quantities (ie including E^n+1^_i). The v,,g,, and kappa,,g,, terms should come from the hydro boundary conditions so we just need an equation for E^n^,,g,, 653 648 654 649 … … 658 653 [[latex(\frac{c \lambda}{\kappa_{0R}} \nabla E = \mathbf{F} = cE\mathbf{n} = \frac{c \lambda_g}{\kappa_g}\frac{ \left ( E-E_g \right ) }{\Delta x})]] 659 654 660 Clearly if we set [[latex(E_g = 0)]] and [[latex(\kappa_g << \frac{1}{\Delta x} )]] we should get [[latex(\lambda = \kappa_g \Delta x)]] and a free streaming flux of [[latex(c E \mathbf{n})]]655 Clearly if we set \(E_g = 0\) and \(\kappa_g << \frac{1}{\Delta x} \) we should get \(\lambda = \kappa_g \Delta x\) and a free streaming flux of \(c E \mathbf{n}\) 661 656 662 657 which also corresponds to an 663 [[latex(\alpha = c \frac{\Delta t}{\Delta x})]] 658 \(\alpha = c \frac{\Delta t}{\Delta x}\) 664 659 665 660 === Constant Slope Boundary ===
