CFD Online Discussion Forums

CFD Online Discussion Forums (
-   Main CFD Forum (
-   -   Crank Nicolson Solution to 3d Heat Equation (

Sharpybox April 4, 2013 07:08

Crank Nicolson Solution to 3d Heat Equation
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:
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 h_c is the convective heat transfer coefficient and T_a
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?

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?

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,

FMDenaro April 4, 2013 07: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 08:05

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

All times are GMT -4. The time now is 02:19.