Thermal Conduction Solver

!!Update 5.22.2015

In the newer version of the code, the unit of temperature in the code or Info%q(iTemp) was changed to Kelvin from C.U.. However the following equations and source code keep the same except the term "gamma7*rho" has to be changed to function getdEdT() which gives C.U energy/Kelvin.. and no gamma7 appears…

Before temperature unit was changed to Kelvin

The equation we are trying to solve in AstroBEAR is

begin{equation} \rho c_v \frac{\partial T}{\partial t} = \nabla \left ( \kappa T^{n} \nabla T \right ) \end{equation}

which can be "Laplacetasized" as

setcounter{equation} 1 \begin{equation} \rho c_v \frac{\partial T}{\partial t} =\nabla^2 \left ( \frac{1}{n+1} \kappa T^{n+1} \right ) \end{equation}

For backward euler we evaluate the RHS at time j+1 although we still have to linearize

T_i^{{j+1}^{n+1}}

So we can Taylor expand about T_i^j

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 )

now if we substitute that back into the discretized form of equation 2 we get

setcounter{equation} 2 \begin{eqnarray} \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 ) \\ & - & \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 ) \end{eqnarray}

Normal CK scheme

We can identify the left and right coefficients as

A_\pm = -\frac{\kappa \Delta t}{\Delta x^2}T_{i\pm1}^{j^n}

and the center coefficient as

A_0 = 2 \frac{\kappa \Delta t}{ \Delta x^2}T_{i}^{j^n} + \rho c_v

and the RHS left and right contributions

B_\pm = -\frac{n \kappa \Delta t}{ \left ( n+ 1 \right ) \Delta x^2} T_{i\pm 1}^{j^{n+1}}

and the RHS center term contribution as

B_0 = 2 \frac{n \kappa \Delta t}{ \left ( n+ 1 \right ) \Delta x^2} T_{i}^{j^{n+1}} + \rho c_v T_i^j

and finally in the code the values for A and B are negated:

    kx = kappa1*dt_diff/(dx**2)
    CALL getTemp(Info,i,j,k,T)
    stencil_fixed(0) = -2.0*REAL(nDim,8)*kx*T(0)**ndiff
    DO l = 1, 2*nDim
       stencil_fixed(l) = kx*T(l)**ndiff
    END DO
    source = ((ndiff)/(ndiff+1.0))*dot_product(T(0:2*nDim),stencil_fixed(0:2*nDim)) - T(0)*gamma7*Info%q(i,j,k,1)
    stencil_fixed(0) = stencil_fixed(0) - gamma7*Info%q(i,j,k,1)

also note that the equation is multipled by the mean molecular weight… which is why c_v \rightarrow c_v \times \chi = \gamma_7 and \kappa \rightarrow \kappa \times \chi = \kappa_1

a note about flux boundary conditions

If we have a specified flux at a boundary - then the equation is just

\rho c_v \frac{\partial T}{\partial t} = -\nabla Q

which can be discretized as

\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 )

which leads to a center coefficient of

A_0 = \rho c_v

and a left right source term contribution

B_\pm = \mp \frac{\Delta t}{\Delta x} Q_{i\pm 1/2}

Now if the left boundary is a flux boundary, but the right side is normal, then we have to adjust the coefficients

A_+ = -\frac{\kappa \Delta t}{\Delta x^2}T_{i+1}^{j^n}

A_- = 0

and the center coefficient as

A_0 = \frac{2\kappa \Delta t}{ \Delta x^2}T_{i}^{j^n} + \rho c_v

and the RHS left and right contributions

B_+ = -\frac{n \kappa \Delta t}{ 2 \left ( n+ 1 \right ) \Delta x^2} T_{i+1}^{j^{n+1}}

B_- = \frac{\Delta t}{\Delta x} Q_{i-1/2}

and the RHS center term contribution as

B_0 = \frac{2n \kappa \Delta t}{ \left ( n+ 1 \right ) \Delta x^2} T_{i}^{j^{n+1}} + \rho c_v T_i^j

Of course - remembering that in AstroBEAR everything is negated…

This is accomplished by

     source = source + (ndiff/(ndiff+1.0))*(T(0)*kx*T(0)**ndiff)
     source = source - (ndiff/(ndiff+1.0))*(T(p)*kx*T(p)**ndiff)
     stencil_fixed(0)=stencil_fixed(0)+kx*T(0)**ndiff
     stencil_fixed(p)=0.0
     source = source - flb*dt_diff/dx

If the right boundary is a zero-slope boundary, we have

T_{i+1}^{j+1}=T_{i}^{j+1} and T_{i+1}^{j}=T_{i}^{j}

So now on the right boundary, the equation becomes

setcounter{equation} 4 \begin{eqnarray} \rho_i c_v \left ( T_i^{j+1}-T_i^j \right ) & = & \frac{\kappa \Delta t}{\Delta x^2}\left (- T_{i}^{j^n} T_{i}^{j+1} + T_{i-1}^{j^n} T_{i-1}^{j+1} \right ) \\ & - & \frac{n\kappa \Delta t}{ \left ( n+ 1 \right ) \Delta x^2} \left (- T_{i}^{j^{n+1}} + T_{i-1}^{j^{n+1}} \right ) \end{eqnarray}

And the coeeficients becomes

A_+ = 0

A_-=-\frac{\kappa \Delta t}{\Delta x^2}T_{i-1}^{j^n}

A_0=\frac{\kappa \Delta t}{\Delta x^2}T_{i}^{j^n}+ \rho c_v

and the RHS left and right contributions

B_+=0

B_-=-\frac{n \kappa \Delta t}{ \left ( n+ 1 \right ) \Delta x^2} T_{i\pm 1}^{j^{n+1}}

B_0=\frac{n \kappa \Delta t}{ \left ( n+ 1 \right ) \Delta x^2} T_{i}^{j^n}+\rho C_vT_{i}^{j}

Again in AstroBEAR everything is negated

This is accomplished by

stencil_fixed(0) =stencil_fixed(0)+kx*info%q(iq(1),iq(2),iq(3),itemp)**ndiff ! store the stencil value to temp
ource = source -(ndiff)/(ndiff+1.0)*(stencil_fixed(p)*info%q(ip(1),ip(2),ip(3),itemp)+kx*info%q(iq(1),iq(2),iq(3),itemp)**(ndiff+1))
stencil_fixed(p) = 0.0 ! zero out the stencil

2nd order in Time CK scheme

Now to make the scheme 2nd order in time we can set T_i^{j+1}=\frac{T_i^{j}+T_i^{j+1}}{2} and we get

setcounter{equation} 4 \begin{eqnarray} \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 ) \\ & - & \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 ) \end{eqnarray}

Where we can identify the left and right coefficients as

A_\pm = -\frac{\kappa \Delta t}{2 \Delta x^2}T_{i\pm1}^{j^n}

and the center coefficient as

A_0 = 2 \frac{\kappa \Delta t}{ 2\Delta x^2}T_{i}^{j^n} + \rho c_v

and the RHS left and right contributions

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}}

and the RHS center term contribution as

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

and finally in the code the values for A and B are negated:

    kx = (0.5*kappa1*dt_diff/(dx**2))
    CALL getTemp(Info,i,j,k,T)
    stencil_fixed(0) = -2.0*REAL(nDim,8)*kx*T(0)**ndiff
    DO l = 1, 2*nDim
       stencil_fixed(l) = kx*T(l)**ndiff
    END DO  
    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)
    stencil_fixed(0) = stencil_fixed(0) - gamma7*Info%q(i,j,k,1)

also note that the equation is multipled by the mean molecular weight… which is why c_v \rightarrow c_v \times \chi = \gamma_7 and \kappa \rightarrow \kappa \times \chi = \kappa_1

a note about flux boundary conditions

If we have a specified flux at a boundary - then the equation is just

\rho c_v \frac{\partial T}{\partial t} = -\nabla Q

which can be discretized as

\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 )

which leads to a center coefficient of

A_0 = \rho c_v

and a left right source term contribution

B_\pm = \mp \frac{\Delta t}{\Delta x} Q_{i\pm 1/2}

Now if the left boundary is a flux boundary, but the right side is normal, then we have to adjust the coefficients

A_+ = -\frac{\kappa \Delta t}{2 \Delta x^2}T_{i+1}^{j^n}

A_- = 0

and the center coefficient as

A_0 = \frac{\kappa \Delta t}{ 2\Delta x^2}T_{i}^{j^n} + \rho c_v

and the RHS left and right contributions

B_+ = -\frac{\left (n - 1 \right ) \kappa \Delta t}{ 2 \left ( n+ 1 \right ) \Delta x^2} T_{i+1}^{j^{n+1}}

B_- = \frac{\Delta t}{\Delta x} Q_{i-1/2}

and the RHS center term contribution as

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

Of course - remembering that in AstroBEAR everything is negated…

This is accomplished by

     source = source + ((ndiff-1.0)/(ndiff+1.0))*(T(0)*kx*T(0)**ndiff)
     source = source - ((ndiff-1.0)/(ndiff+1.0))*(T(p)*kx*T(p)**ndiff)
     stencil_fixed(0)=stencil_fixed(0)+kx*T(0)**ndiff
     stencil_fixed(p)=0.0
     source = source - flb*dt_diff/dx

If the right boundary is a zero-slope boundary, we have

T_{i+1}^{j+1}=T_{i}^{j+1} and T_{i+1}^{j}=T_{i}^{j}

So now on the right boundary, the equation becomes

setcounter{equation} 4 \begin{eqnarray} \rho_i c_v \left ( T_i^{j+1}-T_i^j \right ) & = & \frac{\kappa \Delta t}{2 \Delta x^2}\left (- T_{i}^{j^n} T_{i}^{j+1} + T_{i-1}^{j^n} T_{i-1}^{j+1} \right ) \\ & - & \frac{\left (n - 1 \right ) \kappa \Delta t}{ 2 \left ( n+ 1 \right ) \Delta x^2} \left (- T_{i}^{j^{n+1}} + T_{i-1}^{j^{n+1}} \right ) \end{eqnarray}

And the coeeficients becomes

A_+ = 0

A_-=-\frac{\kappa \Delta t}{2 \Delta x^2}T_{i-1}^{j^n}

A_0=\frac{\kappa \Delta t}{2 \Delta x^2}T_{i}^{j^n}+ \rho c_v

and the RHS left and right contributions

B_+=0

B_-=-\frac{\left (n - 1 \right ) \kappa \Delta t}{ 2 \left ( n+ 1 \right ) \Delta x^2} T_{i\pm 1}^{j^{n+1}}

B_0=\frac{\left (n - 1 \right ) \kappa \Delta t}{ 2 \left ( n+ 1 \right ) \Delta x^2} T_{i}^{j^n}+\rho C_vT_{i}^{j}

Again in AstroBEAR everything is negated

This is accomplished by

stencil_fixed(0) =stencil_fixed(0)+kx*info%q(iq(1),iq(2),iq(3),itemp)**ndiff ! store the stencil value to temp
!source = source -((ndiff-1.0)/(ndiff+1.0))*(stencil_fixed(p)*info%q(ip(1),ip(2),ip(3),itemp)-kx*info%q(iq(1),iq(2),iq(3),itemp)**(ndiff+1))
source = source -((ndiff-1.0)/(ndiff+1.0))*(stencil_fixed(p)*info%q(ip(1),ip(2),ip(3),itemp)+kx*info%q(iq(1),iq(2),iq(3),itemp)**(ndiff+1))
stencil_fixed(p) = 0.0 ! zero out the stencil

where the commented out line is the old-version which got the sign before the "kx" wrong.

Limiting case

In the limit of T \rightarrow 0 with a flux from the left boundary:

our equation becomes \rho_i c_v T_i^{j+1} = \rho_i c_v T_i^j + \frac{\Delta t}{\Delta x} Q \delta_i^1

and the total energy E^{j+1}=\rho c_v \sum T_i^{j+1} \Delta x = E^{j} + \Delta t Q

or \dot{E} = Q