# How to accelerate the computing of implicit finite difference method for heat conduction between two solids

Edit on May 3rd: I have found the problem. Because the difference of between $$k_1$$ and $$k_2$$ is huge, a very small time step need to be chosen so that the right green part can "feel" the heat from the left. Otherwise the temperature profile in the green region would barely change. However, this requres a lot of computation resources especially when I non-dimensionalize the problem. In the following figure, the blue line indicate the condiction that $$k_1=k_2$$, whereas the red line is $$k1<. In obtaining the blue line, I chose the non-dimensional time step $$\Delta\tau$$ to be 6e-11 to pass the grid sensitivity test, but as you can imagine this took a lot of time. So is there any technic to speed up the calculation? Can I compute the two sides separately?

I am solving a transient heat condcution problem involving 2 solids using implicit finite difference method.

At the interface, I soppose the finite difference equation at the interface should be $$\frac{k_1(T_{i-1}-T_i)}{\Delta x}=\frac{k_2(T_{i}-T_{i+1})}{\Delta x}$$ or $$\frac{1}{2}(\rho_1Cp_1+\rho_2Cp_2)\frac{T^{n+1}_i-T^n_i}{\Delta t}=k_1\frac{T^{n+1}_{i-1}-T^{n+1}_i}{\Delta x^2}+k_2\frac{T^{n+1}_{i+1}-T^{n+1}_i}{\Delta x^2}$$ The thernmal conducvivities of the 2 solids are 0.5 and 400, and $$\rho_1Cp_1$$ and $$\rho_2Cp_2$$ are assumed to be equal. But when I solved the problem using the following code (in Matlab), the result is unexpected. That is, the heat transfer performance of the 2 solids with $$k_1=0.5$$ and $$k_2=400$$ is worse than $$k_1=k_2=0.5$$.

So I am wondering if there is any special criteria need to be applied if the difference between $$k_1$$ and $$k_2$$ is huge?

clear all
L=1;
d=0.5;%model dimension
dx=d/40;%grid size

GN=d/dx+1;%interface position
N=L/dx+1; %total grid number

rho1=900;
cp1=2000;
k1=0.5;
rho2=900;
cp2=2000;
k2=400;
a1=k1/(rho1*cp1);
a2=k2/(rho2*cp2);%thermomechanical properties

Time=100000; %total time
dt=1;%time step
alpha=dx^2/dt;

T0=273.15+25;%initial temperature
Tb=273.15+100;%boundary temperature
Ntime=round(Time/dt);%time steps
theta=zeros(Ntime,N);%temperature field history
T=T0*ones(1,N);%initial temperature
T(:,1)=Tb;%LHS B.C.
T2=T;
R=1;

for p=1:1:Ntime
p
theta(p,:)=T;
while (R>10e-5)
T3=T2;
for i=2:1:GN-1
T2(i)=(a1*(T2(i-1)+T2(i+1))+alpha*T(i))/(alpha+2*a1);
end
T2(GN)=(k1*T2(GN-1)+k2*T2(GN+1))/(k1+k2);
for i=GN+1:1:N-1
T2(i)=(a2*(T2(i-1)+T2(i+1))+alpha*T(i))/(alpha+2*a2);
end
T2(N)=(a2*(T2(N-1)+T2(N-1))+alpha*T(N))/(alpha+2*a2);

R=max(abs(T3-T2)) ;

end
T=T2;
R=1;
end

• Checking your code, you are not taking into account the value of the thermal diffusivity into account to compute the time step. Commented May 1, 2022 at 15:20
• Thank you, nico. Commented May 1, 2022 at 22:10