[UmbertoBitencourt]

Dear colleagues.

I've been struggling with this problem for quite some time. I appreciate any help!

I'm trying to find the temperature distribution between two parallel plates. The heat transfer in the plates is specified and known. I'm showing the version solved for 0<x<L, 0<y<H

I'm using the turbulent expression:
U(y)dT/dx=d/dy[(alfa+alfa_t(y))dT/dy]

Where U(y) and alfa_t(y) are known functions. I'm using the Van Driest Eq. for the alfa_t(y).

The BC's are for T(x,y):
T(0,y)=T_0;
dT/dy(x,0)=0;
dT/dy(x,H)=q_0/k;

I wrote a code on MATLAB as follows.

%SS Temperature profile with convection in the x-direction

clear
clc

%Geometric properties
dx=0.02;
dy=0.0005;
L=1;
H=0.15;
nx=round(L/dx);
ny=round(H/dy);
x=linspace(0,L,nx+1);
y=linspace(0,H,ny+1);
beta=dx/dy^2;

%Empirical constants
k=0.4;

%Data imported from ANSYS for comparison
num=xlsread('Temp_V_line2.xlsx');
m=length(num);
y1=abs(num(m:-1:1,1));
v=num(m:-1:1,3);
T_r=num(m:-1:1,2);

%Fluid properties
alfa=2.141*10^-5;
k_a=0.02551;
nu=1.562*10^-5;

%Velocity profile solved separately
U=ones(nx+1,ny+1);
v_max=11.7;
for i=1:1:nx+1
    for j=1:1:ny+1
            U(i,j)=v_max*(1-(y(j)/H)^20);
    end;
end;

%Initial condition
T_0=300;

%Temperature profile
T=T_0*ones(nx+1,ny+1);

%Maximum error
e=10^-6;

%Constants to initiate the iterations
q_0=500;
err=1;
w=0.01;
coun=0;
alfa_tn=alfa;
alfa_ts=alfa;
while err>e
    Ta=T;
for i=2:1:nx
    for j=2:1:ny
     
        %Calculating alfa_t in n
        dudy_n=0.5/dy*abs(U(i,j+1)-U(i,j)+U(i+1,j+1)-U(i+1,j));
        y_n=(y(j+1)+y(j))*0.5;
        alfa_tn=k^2*(H-y_n)^2*dudy_n*(1-exp(-1/26*((H-y_n)*sqrt(nu*dudy_n))));
        
        %Calculating alfa_t in s
        dudy_s=0.5/dy*abs(U(i,j)-U(i,j-1)+U(i+1,j)-U(i+1,j-1));
        y_s=(y(j)+y(j-1))*0.5;
        alfa_ts=k^2*(H-y_s)^2*dudy_s*(1-exp(-1/26*((H-y_s)*sqrt(nu*dudy_s))));
        alfa_n=alfa+alfa_tn;        
        alfa_s=alfa+alfa_ts;
        m_e=dy*(U(i+1,j+1)+U(i+1,j))*0.5;
        m_w=dy*(U(i,j)+U(i,j+1))*0.5;
        
        if j==ny
                T(i,j)=1/(m_w+dx/dy*(alfa_s))*(m_w*T(i-1,j)+dx*alfa_n*q_0/k_a+dx/dy*T(i,j-1)*alfa_s);
        else
        if j==2
                T(i,j)=1/(m_w+dx/dy*(alfa_n))*(m_w*T(i-1,j)+dx/dy*alfa_n*T(i,j+1));
        else
                const=1/(m_w+dx/dy*(alfa_n+alfa_s));
                T(i,j)=(m_w*T(i-1,j)+dx/dy*(alfa_n*T(i,j+1)+alfa_s*T(i,j-1)))*const;
        end;
        end;
    end;
end;
err=sqrt((sum(sum(abs(T-Ta)))))
coun=coun+1;
if coun>50000
    err=0;
end;
end;

figure
y2=0:dy:H-dy;
plot(T(round(nx/2),1:1:ny),y2,'b')
hold on
plot(T_r,y1,'r')

I'm updating my answer as well, where the redline is what the answer should be and the blue one is what I'm getting.

Does anyone spot a mistake? Or any suggestions of what I should do?

I really appreciate any help.