Changes between Version 3 and Version 4 of ThermalConduction

Timestamp:

01/31/14 10:42:37 (11 years ago)

Author:

Baowei Liu

Comment:

—

ThermalConduction

@@ -40 +40 @@
 }}}
 
-
-
+== Normal CK scheme ==
+We can identify the left and right coefficients as
+
+[[latex($A_\pm = -\frac{\kappa \Delta t}{\Delta x^2}T_{i\pm1}^{j^n} $)]]
+
+and the center coefficient as 
+
+[[latex($A_0 =  2 \frac{\kappa \Delta t}{ \Delta x^2}T_{i}^{j^n}  + \rho c_v$)]]
+
+and the RHS left and right contributions
+
+[[latex($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 
+
+[[latex($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 [[latex($A$)]] and [[latex($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 [[latex($c_v \rightarrow c_v \times  \chi = \gamma_7$)]] and [[latex($\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 
+
+[[latex($\rho c_v \frac{\partial T}{\partial t} = -\nabla Q $)]]
+
+which can be discretized as 
+
+[[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 ) $)]]
+
+which leads to a center coefficient of 
+
+[[latex($A_0 = \rho c_v$)]]
+
+and a left right source term contribution
+
+[[latex($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 ===
+
+[[latex($A_+ = -\frac{\kappa \Delta t}{\Delta x^2}T_{i+1}^{j^n} $)]]
+
+[[latex($A_- = 0 $)]]
+
+and the center coefficient as 
+
+[[latex($A_0 =  \frac{2\kappa \Delta t}{ \Delta x^2}T_{i}^{j^n}  + \rho c_v$)]]
+
+and the RHS left and right contributions
+
+[[latex($B_+ = -\frac{n \kappa \Delta t}{ 2  \left ( n+ 1 \right ) \Delta x^2} T_{i+1}^{j^{n+1}}$)]]
+
+[[latex($B_-  =  \frac{\Delta t}{\Delta x} Q_{i-1/2}$)]]
+
+and the RHS center term contribution as 
+
+[[latex($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 ==
+[[latex($T_{i+1}^{j+1}=T_{i}^{j+1}$)]]
+and 
+[[latex($T_{i+1}^{j}=T_{i}^{j}$)]]
+
+So now on the right boundary, the equation becomes
+{{{
+#!latex
+\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
+
+[[latex($A_+ = 0 $)]]
+
+[[latex($A_-=-\frac{\kappa \Delta t}{\Delta x^2}T_{i-1}^{j^n} $)]]
+
+[[latex($A_0=\frac{\kappa \Delta t}{\Delta x^2}T_{i}^{j^n}+ \rho c_v$)]]
+
+and the RHS left and right contributions
+
+[[latex($B_+=0 $)]]
+
+[[latex($B_-=-\frac{n \kappa \Delta t}{ \left ( n+ 1 \right ) \Delta x^2} T_{i\pm 1}^{j^{n+1}}$)]]
+
+[[latex($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 [[latex($T_i^{j+1}=\frac{T_i^{j}+T_i^{j+1}}{2}$)]] and we get
 
@@ -154 +275 @@
 [[latex($A_-=-\frac{\kappa \Delta t}{2 \Delta x^2}T_{i-1}^{j^n} $)]]
 
-[[latex($A_0=\frac{\kappa \Delta t}{2 \Delta x^2}T_{i}^{j^n}T_{i}^{j+1}+ \rho c_v$)]]
+[[latex($A_0=\frac{\kappa \Delta t}{2 \Delta x^2}T_{i}^{j^n}+ \rho c_v$)]]
 
 and the RHS left and right contributions