Heat conduction with non-uniform conductivity
 

I am trying to solve this heat conduction problem :

(K T_x)_x + f =0 where "_" means derivative w.r.t x . K is the conductivity which is a function of x and f is a heat sink. I use the SOR method in solving Poisson equation and write :

T(1) = T(N) =0
for T= 2:N-1
T(j) = omega*(A1* T(j+1) + A2*T(j+2) + A3*f) + (1-omega)*T(j)

where A1 = 0.5 +0.25*K_x*dx/K , A2 = 0.5 -0.25*K_x*dx/K and A3 = -dx^2/2K

This works perfectly fine as long as K is a constant. But when K is not constant, (let's say K(x) =x , -> K_x =1) then it does not work.x

Does anyone have any idea?

Well, your indices are shifted. I doubt that you intended to update T(j) based on a stencil over j, j+1, and j+2. Diffusion should always be treated as symmetrically as possible.

Second, you need to read about harmonic averaging of diffusivities. Section 4 of this paper covers it:

http://web.cecs.pdx.edu/~gerry/class...iffusion2D.pdf

Thanks for your comments.

That was just a typo in my writing. What I actually have is

T(j) = omega*(A1* T(j+1) + A2*T(j-1) + A3*f) + (1-omega)*T(j)

The paper you send me is awesome. and I am looking into implementing it.
But by any chance do you know any similar paper but with FDM?

Thanks,

Most FD schemes are formulated as finite volume schemes (at least for real applications) because FV formulations can guarantee conservation that simple FD formulations do not. These are sometimes called "Conservative Finite Difference" or "Control Volume Finite Difference" methods. By just inserting FD approximations for the derivatives, you cannot know that, say, energy in this system is conserved. Mass, energy, and species conservations are VERY important, especially when dealing with reacting flows. But I digress...

I've verified your FD derivation (central differencing for the dT/dx term and standard 2nd-order differencing for the d^2T/dx*2 term), except I have an opposite sign for the source term coefficient A3...that won't effect stability of the solution, of course. With fixed boundary conditions at zero, it will just flip the sign of the overall solution.

My suspicion is that the matrix loses diagonal dominance as you change the conductivity in space. In the constant K case, the dk/dx terms all vanish and you get the standard {1 -2 1} second derivative FD stencil. That, of course, is diagonally dominate and SOR will work fine.

But if | K_x*dx/K | > 2.0, the off-diagonal coefficients will sum (in magnitude) to a value larger than the center coefficient and SOR is not guaranteed to converge. Even if the matrix is still diagonally dominate, the change in K in space will likely force you to choose a more conservative relaxation factor than you used for the constant K case.

There are implicit under-relaxation techniques that can allow GS or Jacobi iterative schemes to converge even when the underlying system is not diagonally dominate, but for a 1D system like this, you should just use the Thomas Algorithm and solve the system exactly.

Good luck.

I agree to check for the conditions for the convergence of the SOR matrix...furthermore, what about the value of omega? depending on both k(x) and omega, it could be that the spectral radius of the SOR matrix becomes >1