CFD Online Logo CFD Online URL
Home > Forums > Main CFD Forum

Crank Nicolson Solution to 3d Heat Equation

Register Blogs Members List Search Today's Posts Mark Forums Read

LinkBack Thread Tools Display Modes
Old   April 4, 2013, 07:08
Unhappy Crank Nicolson Solution to 3d Heat Equation
New Member
Adam Sharpe
Join Date: Apr 2013
Posts: 2
Rep Power: 0
Sharpybox is on a distinguished road
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,

Last edited by Sharpybox; April 4, 2013 at 08:04. Reason: Mistake
Sharpybox is offline   Reply With Quote

Old   April 4, 2013, 07:35
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,586
Rep Power: 20
FMDenaro will become famous soon enough
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 "+" ?
FMDenaro is offline   Reply With Quote

Old   April 4, 2013, 08:05
New Member
Adam Sharpe
Join Date: Apr 2013
Posts: 2
Rep Power: 0
Sharpybox is on a distinguished road
Many thanks for the response. I have amended the original post
Sharpybox is offline   Reply With Quote


Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On

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

All times are GMT -4. The time now is 15:36.