# Finite Difference Method (FDM) solution to heat equation in material having two different conductivity

I am not from a mechanical engineering background and I have not taken any courses in PDE so this may seem trivial for many. and I am writing a Matlab code with the objective to solve for the steady state temperature distribution in a 2D rectangular material that has 'two phases' of different conductivity.

I was able to do it considering the entire material as one single phase where at each iteration,the value of temperature is updated as the average of temperature of 4 nearest neighbors(central difference method) until the error is less than a specified value between consecutive iterations.

How can I calculate the temperature distribution when I have two phases, taking into consideration the conductivity of each material and how can I make sure that there is temperature continuity at the interface?

• You put the boundary half-way between two grid points. What do you think the difference equation would look like for a grid point immediately to the right of the boundary and for a grid point immediately to the left of the boundary? – Chet Miller May 13 '16 at 3:04
• If the two grid points are (0,j) and (i,j) the boundary lies at the center of the grid (i/2,j). Then dt/dx(@0,j)=Ti,j-Ti-1,j/dx and dt/dx(@i,j)=Ti+1,j-Ti,j/dx. I am sorry that I'm not able to show them as proper equations. – Dualphase May 13 '16 at 7:35
• There is a contact heat transfer coefficient that you need to figure before you march on programming. Also you need to answer if this is a steady state problem or a transient problem, if it has a fixed temperature boundary condition or it is a heat flux or it is insulated or it has fluid flow. Anyway, there is no short answer or solution. – user115350 May 13 '16 at 13:46
• @user115350 Thanks for your comment. I am first considering a steady state problem in 1D before moving on to 3D for my actual problem.Boundary conditions are fixed temperature. Two different materials of different conductivities, no insulation. I am a material science major and I did not think of some of the questions you have raised. Thanks for your insight. – Dualphase May 13 '16 at 18:25
• Doing a 1D model is a great idea to help you understand the heat transfer process. At the interface between the two blocks, the temperature gradients of two blocks are determined by their individual conductivity and the contact thermal conductivity, which is most often determined by test data. – user115350 May 13 '16 at 18:58

Let the discontinuity in thermal conductivity be located half-way between grid points i and i+1. Let the conductivity to the left of the discontinuity be $k_L$ and the conductivity to the right of the discontinuity be $k_R$. Then, for a grid point at i,j (immediately to the left of the discontinuity), the steady state heat balance equation (assuming a square grid) would be:
$$k_L(T_{i-1,j}-T_{i,j})+\left[\frac{2}{(\frac{1}{k_L}+\frac{1}{k_R})}\right](T_{i+1,j}-T_{i,j})+ k_L(T_{i,j+1}-2T_{i,j}+T_{i,j-1})=0$$ If you divide this equation by $k_L$ and solve for $T_{i,j}$, you will obtain the equation you desire.
SUPPLEMENT Let q be the heat flux from grid point i to grid point i + 1. This is given locally by the equation: $$-k\frac{\partial T}{\partial x}=q$$ Solving for $\frac{\partial T}{\partial x}$, we have:$$\frac{\partial T}{\partial x}=-\frac{q}{k}$$ Integrating this equation between i and i+1, we have: $$(T_{i+1,j}-T_{i.j})=-\int_{x_i}^{x_{i+1/2}}\frac{q}{k}dx-\int_{x_{i+1/2}}^{x_{i+1}}\frac{q}{k}dx=-\frac{q}{k_L}\frac{\Delta x}{2}-\frac{q}{k_R}\frac{\Delta x}{2}$$ So, $$q=-\left[\frac{2}{(\frac{1}{k_L}+\frac{1}{k_R})}\right]\frac{(T_{i+1,j}-T_{i,j})}{\Delta x}$$