Crank Nicolson Solution to 3d Heat Equation

Sharpybox · April 4, 2013, 08:08

Hello everyone. I am currently trying to create a Crank Nicolson solver to model the temperature distribution within a Solar Cell with heat sinking arrangement and have three question I would like to ask about my approach.
I am aiming to solve the 3d transient heat equation:
$\frac{\partial T}{\partial t}$ = $\alpha$ ( $\nabla^{2}$ T + $\frac{q}{\kappa}$ )
Where T = temperature (K), t = time (s), q is the rate of heat generation (W/m^3), $\kappa$ = thermal conductivity (W/mK) and $\alpha$ = thermal diffusivity (m^2/s).

Approximating the time derivative by a backward difference approximation and the spatial derivatives by central difference second order approximation yields my Crank Nicolson scheme for the bulk of the material:

(1+3 $\mu$ ) $T^{t+1}_{i, j, k}$ = $\frac{\mu}{2}$ ( $T^{t+1}_{i+1, j, k}$ + $T^{t+1}_{i-1, j, k}$ + $T^{t+1}_{i, j+1, k}$ + $T^{t+1}_{i, j-1, k}$ $T^{t+1}_{i, j, k+1}$ + $T^{t+1}_{i, j, k-1}$ ) + (1-3 $\mu$ ) $T^{t}_{i, j, k}$ + $\frac{\mu}{2}$ ( $T^{t}_{i+1, j, k}$ + $T^{t}_{i-1, j, k}$ + $T^{t}_{i, j+1, k}$ + $T^{t}_{i, j-1, k}$ $T^{t}_{i, j, k+1}$ + $T^{t}_{i, j, k-1}$ ) + $\frac{\alpha q k}{\kappa}$

Where i, j and k are unit vectors in the x, y, z directions respectively and t+1 and t refer to the current and previous time steps and $\mu$ = $\frac{k \alpha}{h^{2}}$ (with k = the step size between adjacent time steps and h = the step size for each spatial direction).

My three questions are as follows:
1.
I would like to model the boundary of my domain as losing heat to it's environment by convection and radiation. I am having a few problems implementing the boundary conditions so would like to check my approach.
As I understand it this is a Neumann condition and can be represented by (assuming a 1d case along the i=0 boundary):
- $\kappa$ $\frac{\partial^{} T}{\partial x^{}}$ = $h_{c}$ $(T^{t}_{i,j,k} - T_a)$ + $\epsilon$ $\sigma$ ( $T^{t4}_{i,j,k}$ - $T^{4}_{a}$ )

Where

is the convective heat transfer coefficient and

Applying a central difference to the derivative and rearranging for the 'ghost node' $T^{t}_{i-1, j, k}$ yields:
$T^{t}_{i-1, j, k}$ = $T^{t}_{i+1, j, k}$ - $\frac{2h (h_c(T^{t}_{i,j,k} - T_a) + \epsilon \sigma(T^{t4}_{i,j,k}-T^{4}_{a})}{\kappa}$

Am I correct in thinking that this can be substituted into the original Crank Nicolson Scheme in place of the 'ghost node' along the boundary? And should the same approach be implemented to approach the ghost nodes at both the t+1 and the t time steps?

2.
With the approach shown I am finding that my energy generation term produces the same temperature rise regardless of the time step I employ, which isn't consistent logically. With the approach above if I increase the time step by a factor x, then all terms in the equation are increased by the same factor (since $\mu$ is also increased by the same proportion) and this cancels out and produces the same temperature rise at a given node regardless of time step. Can anyone spot anywhere obvious that I have gone wrong?

3.
I am employing the Gauss-Seidel iterative technique to solve for each of the t+1th nodes (I appreciate it isn't the most efficient but is a useful starting point for me) and am finding that the number of iterations required increases consistently as the number of completed time steps increases. This is leading to a situation whereby each time step is taking an impractical number of iterations to converge (up to several thousand). I assume this is some issue with accumulation of errors but as a novice programmer am a bit perplexed as to how to stop this happening. If anyone can shed any light that would be fantastic, although I appreciate that the information given here is somewhat vague!

Apologies for the long winded post and thank you in advance for your help.
All the best,
Adam

FMDenaro · April 4, 2013, 08:35

concerning issue #1, I suggest to write the Laplacian as Div (Grad) so that you can substitute directly the Neumann condition at boundary....

I haven't read yet the other two questions

PS: in your CN formula you wrot "=" instead of "+" ?

Sharpybox · April 4, 2013, 09:05

Many thanks for the response. I have amended the original post

April 4, 2013, 08:08	Crank Nicolson Solution to 3d Heat Equation	#1
Sharpybox New Member Adam Sharpe Join Date: Apr 2013 Posts: 2 Rep Power: 0	Hello everyone. I am currently trying to create a Crank Nicolson solver to model the temperature distribution within a Solar Cell with heat sinking arrangement and have three question I would like to ask about my approach. I am aiming to solve the 3d transient heat equation: $\frac{\partial T}{\partial t}$ = $\alpha$ ( $\nabla^{2}$ T + $\frac{q}{\kappa}$ ) Where T = temperature (K), t = time (s), q is the rate of heat generation (W/m^3), $\kappa$ = thermal conductivity (W/mK) and $\alpha$ = thermal diffusivity (m^2/s). Approximating the time derivative by a backward difference approximation and the spatial derivatives by central difference second order approximation yields my Crank Nicolson scheme for the bulk of the material: (1+3 $\mu$ ) $T^{t+1}_{i, j, k}$ = $\frac{\mu}{2}$ ( $T^{t+1}_{i+1, j, k}$ + $T^{t+1}_{i-1, j, k}$ + $T^{t+1}_{i, j+1, k}$ + $T^{t+1}_{i, j-1, k}$ $T^{t+1}_{i, j, k+1}$ + $T^{t+1}_{i, j, k-1}$ ) + (1-3 $\mu$ ) $T^{t}_{i, j, k}$ + $\frac{\mu}{2}$ ( $T^{t}_{i+1, j, k}$ + $T^{t}_{i-1, j, k}$ + $T^{t}_{i, j+1, k}$ + $T^{t}_{i, j-1, k}$ $T^{t}_{i, j, k+1}$ + $T^{t}_{i, j, k-1}$ ) + $\frac{\alpha q k}{\kappa}$ Where i, j and k are unit vectors in the x, y, z directions respectively and t+1 and t refer to the current and previous time steps and $\mu$ = $\frac{k \alpha}{h^{2}}$ (with k = the step size between adjacent time steps and h = the step size for each spatial direction). My three questions are as follows: 1. I would like to model the boundary of my domain as losing heat to it's environment by convection and radiation. I am having a few problems implementing the boundary conditions so would like to check my approach. As I understand it this is a Neumann condition and can be represented by (assuming a 1d case along the i=0 boundary): - $\kappa$ $\frac{\partial^{} T}{\partial x^{}}$ = $h_{c}$ $(T^{t}_{i,j,k} - T_a)$ + $\epsilon$ $\sigma$ ( $T^{t4}_{i,j,k}$ - $T^{4}_{a}$ ) Where is the convective heat transfer coefficient and Applying a central difference to the derivative and rearranging for the 'ghost node' $T^{t}_{i-1, j, k}$ yields: $T^{t}_{i-1, j, k}$ = $T^{t}_{i+1, j, k}$ - $\frac{2h (h_c(T^{t}_{i,j,k} - T_a) + \epsilon \sigma(T^{t4}_{i,j,k}-T^{4}_{a})}{\kappa}$ Am I correct in thinking that this can be substituted into the original Crank Nicolson Scheme in place of the 'ghost node' along the boundary? And should the same approach be implemented to approach the ghost nodes at both the t+1 and the t time steps? 2. With the approach shown I am finding that my energy generation term produces the same temperature rise regardless of the time step I employ, which isn't consistent logically. With the approach above if I increase the time step by a factor x, then all terms in the equation are increased by the same factor (since $\mu$ is also increased by the same proportion) and this cancels out and produces the same temperature rise at a given node regardless of time step. Can anyone spot anywhere obvious that I have gone wrong? 3. I am employing the Gauss-Seidel iterative technique to solve for each of the t+1th nodes (I appreciate it isn't the most efficient but is a useful starting point for me) and am finding that the number of iterations required increases consistently as the number of completed time steps increases. This is leading to a situation whereby each time step is taking an impractical number of iterations to converge (up to several thousand). I assume this is some issue with accumulation of errors but as a novice programmer am a bit perplexed as to how to stop this happening. If anyone can shed any light that would be fantastic, although I appreciate that the information given here is somewhat vague! Apologies for the long winded post and thank you in advance for your help. All the best, Adam Last edited by Sharpybox; April 4, 2013 at 09:04. Reason: Mistake

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Exact solution for 2D inviscid burgers equation.	xiaowanzi01	Main CFD Forum	15	May 17, 2012 14:55
Constant velocity of the material	Sas	CFX	15	July 13, 2010 09:56
Analytical solution of poisson equation	kostaspan	Main CFD Forum	1	April 22, 2010 13:24
Concentric tube heat exchanger (Air-Water)	Young	CFX	5	October 7, 2008 00:17
Latent heat introduced by the TEM1 equation?	H.P.LIU	Phoenics	1	January 11, 2001 11:03

April 4, 2013, 08:35		#2
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,897 Rep Power: 73	concerning issue #1, I suggest to write the Laplacian as Div (Grad) so that you can substitute directly the Neumann condition at boundary.... I haven't read yet the other two questions PS: in your CN formula you wrot "=" instead of "+" ?

April 4, 2013, 09:05		#3
Sharpybox New Member Adam Sharpe Join Date: Apr 2013 Posts: 2 Rep Power: 0	Many thanks for the response. I have amended the original post