2
$\begingroup$

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<<k2$. 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?enter image description here


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

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
$\endgroup$
2
  • 1
    $\begingroup$ Checking your code, you are not taking into account the value of the thermal diffusivity into account to compute the time step. $\endgroup$
    – nicoguaro
    Commented May 1, 2022 at 15:20
  • $\begingroup$ Thank you, nico. $\endgroup$
    – Kai Jiao
    Commented May 1, 2022 at 22:10

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.