Changes between Version 21 and Version 22 of u/johannjc/scratchpad


Ignore:
Timestamp:
04/21/14 13:18:57 (11 years ago)
Author:
Jonathan
Comment:

—

Legend:

Unmodified
Added
Removed
Modified

  • u/johannjc/scratchpad

    v21 v22  
    1 [[latex($\rho c_v \frac{\partial T}{\partial t} = \nabla \left ( \kappa T^{n} \nabla T \right )$)]]
    21
    3 [[latex($\rho c_v \frac{\partial T}{\partial t} =\nabla^2 \left ( \frac{1}{n+1} \kappa T^{n+1} \right )$)]]
     2[[CollapsibleBegin(Single Temperature Solution)]]
     3If we plug the expressions for the radiation 4-momentum back into the gas equations and keep terms necessary to maintain accuracy we get:
    44
    5 For backward euler we evaluate the RHS at time [[latex($j+1$)]] although we still have to linearize
     5  [[latex($\frac{\partial }{\partial t} \left(\rho\mathbf{v}\right)+\nabla\cdot\left(\rho\mathbf{vv}\right)=\nabla  P\color{green}{-\lambda \nabla E}$)]]
     6   
     7  [[latex($\frac{\partial e}{\partial t}+\nabla\cdot\left[\left(e+P\right)\mathbf{v}\right]=\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{blue}{-\frac{3-R_2}{2}\kappa_{0P}\frac{v^2}{c}E}$)]]
    68
    7 [[latex($T_i^{{j+1}^{n+1}} = T_i^{j^{n+1}} + \left (n + 1 \right ) T_i^{j^n} \left ( T_i^{j+1} - T_i^j \right ) $)]]
    8 
    9 
    10 or equivalently
    11 
    12 [[latex($T_i^{{j+1}^{n+1}} = T_i^{j^n} \left (T_i^{j+1}  + n \left (T_i^{j+1} - T_i^{j} \right ) \right )$)]]
    13 
    14 
    15 
    16 {{{
    17 #!latex
    18 \begin{eqnarray}
    19 \rho_i c_v \left ( T_i^{j+1}-T_i^j \right ) & =  & \frac{\kappa \Delta t}{\Delta x^2}\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 ) \\
    20  & - & \frac{n \kappa \Delta t}{ \left ( n+ 1 \right ) \Delta x^2} \left (T_{i+1}^{j^{n+1}} - 2 T_{i}^{j^{n+1}} + T_{i-1}^{j^{n+1}} \right ) 
    21 \end{eqnarray}
    22 }}}
    23 
    24 
    25 
    26 Now to make the scheme 2nd order in time we can set [[latex($T_i^{j+1}=\frac{T_i^{j}+T_i^{j+1}}{2}$)]] and we get
    27 
    28 {{{
    29 #!latex
    30 \begin{eqnarray}
    31 \rho_i c_v \left ( T_i^{j+1}-T_i^j \right ) & =  & \frac{\kappa \Delta t}{2 \Delta x^2}\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 ) \\
    32  & - & \frac{\left (n - 1 \right ) \kappa \Delta t}{ 2 \left ( n+ 1 \right ) \Delta x^2} \left (T_{i+1}^{j^{n+1}} - 2 T_{i}^{j^{n+1}} + T_{i-1}^{j^{n+1}} \right ) 
    33 \end{eqnarray}
    34 }}}
    35 
    36 Where we can identify the left and right coefficients as
    37 
    38 [[latex($A_\pm = -\frac{\kappa \Delta t}{2 \Delta x^2}T_{i\pm1}^{j^n} $)]]
    39 
    40 and the center coefficient as
    41 
    42 [[latex($A_0 =  2 \frac{\kappa \Delta t}{ 2\Delta x^2}T_{i}^{j^n}  + \rho c_v$)]]
    43 
    44 and the RHS left and right contributions
    45 
    46 [[latex($B_\pm = -\frac{\left (n - 1 \right ) \kappa \Delta t}{ 2  \left ( n+ 1 \right ) \Delta x^2} T_{i\pm 1}^{j^{n+1}}$)]]
    47 
    48 and the RHS center term contribution as
    49 
    50 [[latex($B_0 = 2 \frac{\left (n - 1 \right ) \kappa \Delta t}{ 2  \left ( n+ 1 \right ) \Delta x^2} T_{i}^{j^{n+1}} + \rho c_v T_i^j$)]]
    51 
    52 and finally in the code the values for [[latex($A$)]] and [[latex($B$)]] are negated:
    53 
    54 {{{
    55     kx = (0.5*kappa1*dt_diff/(dx**2))
    56     CALL getTemp(Info,i,j,k,T)
    57     stencil_fixed(0) = -2.0*REAL(nDim,8)*kx*T(0)**ndiff
    58     DO l = 1, 2*nDim
    59        stencil_fixed(l) = kx*T(l)**ndiff
    60     END DO 
    61     source = ((ndiff-1.0)/(ndiff+1.0))*dot_product(T(0:2*nDim),stencil_fixed(0:2*nDim)) - T(0)*gamma7*Info%q(i,j,k,1)
    62     stencil_fixed(0) = stencil_fixed(0) - gamma7*Info%q(i,j,k,1)
    63 }}}
    64 
    65 also note that the equation is multipled by the mean molecular weight... which is why [[latex($c_v \rightarrow c_v \times  \chi = \gamma_7$)]] and [[latex($\kappa \rightarrow \kappa \times \chi = \kappa_1$)]]
    66 
    67 
    68 == a note about flux boundary conditions ==
    69 If we have a specified flux at a boundary - then the equation is just
    70 
    71 [[latex($\rho c_v \frac{\partial T}{\partial t} = -\nabla Q $)]]
    72 
    73 which can be discretized as
    74 
    75 [[latex($\rho_i c_v \left (T_i^{j+1}-T_i^j \right ) = \frac{-\Delta t}{\Delta x} \left ( Q_{i+1/2} - Q_{i-1/2} \right ) $)]]
    76 
    77 which leads to a center coefficient of
    78 
    79 [[latex($A_0 = \rho c_v$)]]
    80 
    81 and a left right source term contribution
    82 
    83 [[latex($B_\pm = \mp \frac{\Delta t}{\Delta x} Q_{i\pm 1/2}$)]]
    84 
    85 === Now if the left boundary is a flux boundary, but the right side is normal, then we have to adjust the coefficients ===
    86 
    87 [[latex($A_+ = -\frac{\kappa \Delta t}{2 \Delta x^2}T_{i+1}^{j^n} $)]]
    88 
    89 and the center coefficient as
    90 
    91 [[latex($A_0 =  \frac{\kappa \Delta t}{ 2\Delta x^2}T_{i}^{j^n}  + \rho c_v$)]]
    92 
    93 and the RHS left and right contributions
    94 
    95 [[latex($B_+ = -\frac{\left (n - 1 \right ) \kappa \Delta t}{ 2  \left ( n+ 1 \right ) \Delta x^2} T_{i+1}^{j^{n+1}}$)]]
    96 
    97 [[latex($B_-  =  \frac{\Delta t}{\Delta x} Q_{i-1/2}$)]]
    98 
    99 and the RHS center term contribution as
    100 
    101 [[latex($B_0 = \frac{\left (n - 1 \right ) \kappa \Delta t}{ 2  \left ( n+ 1 \right ) \Delta x^2} T_{i}^{j^{n+1}} + \rho c_v T_i^j$)]]
    102 
    103 
    104 Of course - remembering that in AstroBEAR everything is negated...
    105 
    106 This is accomplished by
    107 {{{
    108      source = source + ((ndiff-1.0)/(ndiff+1.0))*(T(0)*kx*T(0)**ndiff)
    109      source = source - ((ndiff-1.0)/(ndiff+1.0))*(T(p)*kx*T(p)**ndiff)
    110      stencil_fixed(0)=stencil_fixed(0)+kx*T(0)**ndiff
    111      stencil_fixed(p)=0.0
    112      source = source - flb*dt_diff/dx
    113 }}}
    114 
    115 == Limiting case ==
    116 
    117 In the limit of [[latex($T \rightarrow 0$)]] with a flux from the left boundary:
    118 
    119 our equation becomes [[latex($\rho_i c_v T_i^{j+1} = \rho_i c_v T_i^j + \frac{\Delta t}{\Delta x} Q \delta_i^1$)]]
    120 
    121 and the total energy [[latex($E^{j+1}=\rho c_v \sum T_i^{j+1} \Delta x = E^{j} + \Delta t Q$)]]
    122 
    123 or [[latex($\dot{E} = Q$)]]
     9  [[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}  $)]] 
     10[[CollasibleEnd()]]