Trying to solve Temperature profile in 2D Heat Transfer Through Composite Wall with Convective, Constant Temperature and Adiabatic Boundary Condition.

Hello,

I am attempting to code a solution for the temperature profile of a composite wall. The problem set up is the following, a composite wall made of two plastics and aluminum, the composite sections have the following thermal conductivities going from left to right, k1=0.2, k2=237, and k3=0.4 W/m*K. The boundaries are as follows, Left wall is at constant temperature 25C, the bottom wall considered adiabatic, the top and right walls are exposed to ambient atmoshpere at a temperature of T_inf=40C with natural convection. The convection coefficient of the air is assumed to be h_inf=10 W/m^2*K. Heat transfer into a composite interface is equal to the heat transferred away from the composite interface. The top and right side are considered to be facing "outside", and the left wall is considered to be faceing "inside". The composite wall interfaces are located at 1/3 the length of the total wall (m_12@L=1/3 and m_23@L=2/3). This is a practice problem that I created for myself, so I appologize if it dosen't seem that practical.

So there are a couple of issues I am having.

• One is, I am allowing the code to run for a maximum of 100000 iterations to solve. And the code will not solve completely within that allowance. If I do not put a limit it simply runs for a very long time. This is with an initialized guess temperature of 25C.
• Two is, it appears that the convective top wall of the alluminum composite is not calculating correctly, as if it is simply adiabatic, whenever the convective top wall of the other two plastics are functioning properly, and all the of equations used for the convective top wall BC are the same for all three composites. Given that the atmospheric temperature outside of 40 is greater than that of the inside temperature of 25, the middle wall (aluminum, since it is exposed to atmosphere) should be conducting heat from the outside, making it decently hot. This is not the case, the temperature of the middle wall (aluminum) is actually colder than the inside temperature, which should not be possible.
• Three, the interface BC that I am using do not seem to be working correctly, I am not sure if it is simply because the middle wall is the aluminum and seems to perform as some kind of heat sink, or if it is an error with my equation. I am using the following relation for the BC
• 
• And the following is the equation solved for T(m_12,2:n-1) or T(m_23,2:n-1).
• 
• 
• These equations are only dependent on heat transfer in the x-direction, which I am afraid may also be the issue, maybe I need to consider the y-direction heat transfer at the interfaces as well. That being the case, I am not sure how to implement this BC into the general heat equation d2T/dx2+d2T/dy2=0.

The following are my code, and the resulting figure I get when the code has exceeded its iteration count.

%% Initialize 
clear all 
close all
clc
%% Given Temperatures 
Twall=25; %degC
T_inf=40; 
%% Material Properties
k1=0.4; %W/m*K
k2=237;
k3=0.2;
%% Ambient Conditions 
h_inf=10; %W/m^2K
%% Defining Wall Dimensions
H=0.8; %m
L=0.3;
W=0.8;
Comp_length=L/3;
delx=0.001;
dely=0.001;
x=0:delx:L;
y=0:dely:H;
%% Defining Reference Values 
m=length(x);
n=length(y);
m_12=round((1/3)*m);
m_23=round((2/3)*m);
diff=1;
tolerance=0.000001;
k=0;
%% Initializing Temperature
T=ones(length(x),length(y))*25;

%% FDM
while diff > tolerance
    k=k+1;
    Told=T;
    % Constant Temperature @ Left Wall
    T(1,1:n)=Twall;
    % Convection @ Right Wall
    T(m,1:n)=(k3*T(m-1,1:n)/delx+h_inf*T_inf)/(k3/delx+h_inf);
    % Convection @ Top Wall of Composite 1
    T(2:m_12-1,n)=(k1*T(2:m_12-1,n-1)/dely+h_inf*T_inf)/(k1/dely+h_inf);
    % Convection @ Top Wall of Composite 2
    T(m_12:m_23-1,n)=(k2*T(m_12:m_23-1,n-1)/dely+h_inf*T_inf)/(k2/dely+h_inf);
    % Convection @ Top Wall of Composite 3
    T(m_23:m-1,n)=(k3*T(m_23:m-1,n-1)/dely+h_inf*T_inf)/(k3/dely+h_inf);
    % Adiabatic Condition @ Bottom Wall
    T(2:m-1,1)=T(2:m-1,2);
    % Composite 1-2 interface
    T(m_12,2:n-1)=(k1*T(m_12-1,2:n-1)+k2*T(m_12+1,2:n-1))/(k1+k2);
    % Composite 2-3 interface
    T(m_23,2:n-1)=(k2*T(m_23-1,2:n-1)+k3*T(m_23+1,2:n-1))/(k2+k3);
    % Internal Nodes for Composite 1
    T(2:m_12-1,2:n-1)=0.25*(T(3:m_12,2:n-1)+T(1:m_12-2,2:n-1)+T(2:m_12-1,3:n)+T(2:m_12-1,1:n-2));
    % Internal Nodes for Composite 2
    T(m_12+1:m_23-1,2:n-1)=0.25*(T(m_12+2:m_23,2:n-1)+T(m_12:m_23-2,2:n-1)+T(m_12+1:m_23-1,3:n)+T(m_12+1:m_23-1,1:n-2));
    % Internal Nodes for Composite 3
    T(m_23+1:m-1,2:n-1)=0.25*(T(m_23+2:m,2:n-1)+T(m_23:m-2,2:n-1)+T(m_23+1:m-1,3:n)+T(m_23+1:m-1,1:n-2));
    diff=max(max(T-Told));
    if k>100000
        "no good"
        break
    end
end

%% Plotting
T_act=T';
max(T(:))
min(T(:))
k
pcolor(x,y,T'),shading interp
title('Temperature Profile, S-S'), xlabel('x'), ylabel('y'), colorbar

[Torsten]

So we can look at the code, add a few things from which we think they might be important like heat transfer in the direction parallel to the interfaces and boundaries, get a different result and still don't know whether it's more correct than the previous or not.

What you need is a validation tool (like the PDE Toolbox, ANSYS, COMSOL or whatever) in which you set up exactly the same problem in order to compare the results. Then you can make well-directed changes to your code if necessary.